The function dglmstarma estimates a multivariate time series model based on double generalized linear models (DGLM) introduced by Smyth (1989). The primary application is for spatio-temporal data, but different applications, such as network data, are also feasible.
Conditionally on the past, each component of the multivariate time series is assumed to follow a distribution from the exponential dispersion family, see Jørgensen (1987).
In contrast to standard generalized linear models, the dispersion parameter of the distribution is allowed to vary.
The model framework links the mean of the time series conditional on the past, to a linear predictor. This linear predictor allows regression on past observations, past values of the linear predictor and covariates, as described in the details.
Additionally, the dispersion parameter of the distribution is modeled with an additional linear predictor, which can also include spatial and temporal dependencies as well as covariates.
Various distributions with several link-functions are available.
dglmstarma(
ts,
mean_model = list(),
dispersion_model = list(),
mean_family = NULL,
dispersion_link = c("log", "identity", "inverse"),
wlist = NULL,
mean_covariates = list(),
dispersion_covariates = list(),
pseudo_observations = c("deviance", "pearson"),
wlist_past_mean = NULL,
wlist_covariates = NULL,
wlist_pseudo_obs = NULL,
wlist_past_dispersion = NULL,
wlist_covariates_dispersion = NULL,
control = list()
)The function returns an object of class dglmstarma, which includes
mean A list containing information about the mean model, see glmstarma for details. Additionally, it contains:
param_history The sequence of parameters estimates of the mean model during the fitting process.
log_likelihood_history The sequence of log-likelihood evaluations of the mean model during the fitting process.
dispersion Information about the dispersion model, see glmstarma for details. In ts it stores the final pseudo-observations. Additionally, it contains:
pseudo_type Type of pseudo observations used ("deviance" or "pearson").
param_history The sequence of parameters estimates of the dispersion model during the fitting process.
log_likelihood_history The sequence of log-likelihood evaluations of the dispersion model during the fitting process.
target_dim Number of parameters in the model.
algorithm_info Information about the fitting algorithm for each iteration of the inner fitting loop.
convergence_info Information about the convergence of the inner fitting loop.
total_log_likelihood_history Evolution of the log-likelihood during the fitting process.
total_log_likelihood The final log-likelihood of the fitted model.
aic AIC of the (full) model based on the log-likelihood, see information_criteria.
bic BIC of the (full) model based on the log-likelihood, see information_criteria.
qic QIC of the (full) model based on the log-likelihood, see QIC.
call The function call.
control The control parameters used for fitting the model.
Multivariate time series. Rows indicate the locations and columns the time.
a named list specifying the model orders of the linear predictor, which can be of the following elements:
intercept : (Optional) character
'homogenous' (default) for a homogenous model, i.e. the same intercept for all components
'inhomogenous' for inhomogenous models, i.e. fitting an individual intercept for each component
past_obs : (Optional)
Integer vector with the maximal spatial orders for the time lags in past_obs_time_lags.
Alternatively: a binary matrix, with the entry in row \(i\) and column \(j\) indicating whether the \((i - 1)\)-spatial lag for the \(j\)-th time lag is included in the model.
past_obs_time_lags : (Optional) integer vector
indicates the time lags for past_obs. Defaults to seq(length(past_obs)) (for vectors) and seq(ncol(past_obs)) (for a matrix)
past_mean : (Optional)
Spatial orders for the regression on past values of (latent) linear process values.
Values can be entered in the same format as in past_obs. If not specified, no regression to the feedback process is performed.
past_mean_time_lags : (Optional) integer vector
Time lags for the regression on the (latent) linear process. Values can be entered in the same format as in past_obs_time_lags.
covariates : (Optional)
spatial orders for the covariate processes passed in the argument covariates. The values can be passed as in past_obs and past_means, where the \(j\)-th entry or column represents the \(j\)-th covariate.
Default is spatial order 0 for all covariates, which corresponds to the first matrix in argument wlist_covariates.
a named list specifying the model orders of the dispersion linear predictor, which can have the same elements as the mean_model argument.
A list generated by one of the family functions of this package, see stfamily. This argument specifies the marginal conditional distributions of the observations and the type of model fitted for the mean linear predictor.
Link function that is used for the dispersion model. Available options are "log" (default), "identity" and "inverse".
A list of quadratic matrices, with the same dimension as the time series has rows, which describe the spatial dependencies. Row-normalized matrices are recommended. See Details.
List of covariates for the mean linear predictor, containing matrices of same dimension as ts or returns of the covariate functions of this package (see also TimeConstant, SpatialConstant).
List of covariates for the dispersion linear predictor, containing matrices of same dimension as ts or returns of the covariate functions of this package (see also TimeConstant, SpatialConstant).
(character vector) Defines how pseudo observations for the past dispersion values are calculated. Options are "deviance" (default) and "pearson". See Details.
(Optional) List of matrices, which describes spatial dependencies for the values of the linear predictor. If this is NULL, the matrices from wlist are used.
(Optional) List of matrices, which describes spatial dependencies for the covariates. If this is NULL, the matrices from wlist are used.
(Optional) List of matrices, which describes spatial dependencies for past values of the pseudo observations. If this is NULL, the matrices from wlist are used.
(Optional) List of matrices, which describes spatial dependencies for the past dispersion values (latent process). If this is NULL, the matrices from wlist are used.
(Optional) List of matrices, which describes spatial dependencies for the covariates in the dispersion model. If this is NULL, the matrices from wlist are used.
A list of parameters for controlling the fitting process. This list is passed to dglmstarma.control.
For a multivariate time series \(\{Y_t = (Y_{1,t}, \ldots, Y_{p,t})'\}\), we assume that the (marginal) conditional components \(Y_{i,t} \mid \mathcal{F}_{t-1}\), on the past, follow a distribution that is a member of the exponential dispersion family.
The joint multivariate distribution of \(Y_t \mid \mathcal{F}_{t-1}\) is assumed to be generated by a process involving copulas. The distributional assumptions imply that the conditional mean \(\mathbf{\mu}_t := \mathbb{E}(Y_t \mid \mathcal{F}_{t-1})\) and the conditional variance \(\mathrm{Var}(Y_{i,t} \mid \mathcal{F}_{t-1}) = \phi_{i,t} V(\mu_{i,t})\), where \(V(\cdot)\) is the variance function of the chosen distribution and \(\phi_{i,t}\) is the dispersion parameter for location \(i\) at time \(t\).
The conditional mean is linked to a linear process by the link-function, i.e. \(g(\mathbf{\mu}_t) = \mathbf{\psi}_t\), which is applied elementwise.
A second linear process \(\mathbf{\zeta}_t\) is linked to the dispersion parameters of the distributions via a second link-function \(g_d\), i.e. \(g_d(\mathbf{\phi}_t) = \mathbf{\zeta}_t\).
The linear predictor for the mean process is defined by regression on past observations, past values of the linear predictor and covariates. It has the following structure:
$$\mathbf{\psi}_t = \mathbf{\delta} + \sum_{i = 1}^{q} \sum_{\ell = 0}^{a_i} \alpha_{i, \ell} W_{\alpha}^{(\ell)} h(\mathbf{\psi}_{t - i}) + \sum_{j = 1}^{r} \sum_{\ell = 0}^{b_j} \beta_{j, \ell} W_{\beta}^{(\ell)} \tilde{h}(\mathbf{Y}_{t - j}) + \sum_{k = 1}^{m} \sum_{\ell = 0}^{c_k} \gamma_{k, \ell} W_{\gamma}^{(\ell)} \mathbf{X}_{k, t},$$
where the matrices \(W_{\alpha}^{(\ell)}\), \(W_{\beta}^{(\ell)}\), and \(W_{\gamma}^{(\ell)}\) are taken from the lists wlist_past_mean, wlist, and wlist_covariates, respectively, and \(\ell\) denotes the spatial order.
If \(\delta = \delta_0 \mathbf{1}\) with a scalar \(\delta_0\), the model is called homogenous with respect to the intercept; otherwise, it is inhomogenous.
Spatial orders, intercept structure and time lags for the mean model are specified in the argument mean_model. If past_mean is specified, it is also required that past_mean is specified for identifiability.
The linear process of the dispersion model is defined similarly, but instead of direct observations it includes pseudo observations \(\mathbf{d}_t\), which are either defined based on deviance or Pearson residuals.
The linear process of the dispersion model has the following structure:
$$\mathbf{\zeta}_t = \mathbf{\tilde{\delta}} + \sum_{i = 1}^{\tilde{q}} \sum_{\ell = 0}^{\tilde{a}_i} \tilde{\alpha}_{i, \ell} W_{\alpha, \phi}^{(\ell)} \mathbf{\zeta}_{t - i} + \sum_{j = 1}^{\tilde{r}} \sum_{\ell = 0}^{\tilde{b}_j} \tilde{\beta}_{j, \ell} W_{\beta, \phi}^{(\ell)} \tilde{h}_{\phi}(\mathbf{d}_{t - j}) + \sum_{k = 1}^{\tilde{m}} \sum_{\ell = 0}^{\tilde{c}_k} \tilde{\gamma}_{k, \ell} W_{\gamma, \phi}^{(\ell)} \mathbf{\tilde{X}}_{k, t},$$
The model orders, neighborhood structures are specified analogously to the mean model via the argument dispersion_model and the wlist_ arguments.
The unknown parameters of the model are estimated with an iterative procedure, which alternates between estimating the mean and dispersion model until convergence.
Within each step, a quasi-maximum likelihood approach is used, where for the mean model the quasi-log-likelihood of the observations resulting from the mean_family argument is maximized. For the dispersion model, the quasi-likelihood resulting from a Gamma-Density with fixed dispersion parameter of 2 is maximized.
In case of a negative binomial family, the pseudo-observations are always calculated based on Pearson residuals using the Poisson variance function. The dispersion model is defined on these pseudo-observations, which have expectation \(1 + \phi_{i, t} \mu_{i,t}\).
Jørgensen, B. (1987), Exponential Dispersion Models. Journal of the Royal Statistical Society: Series B (Methodological), 49: 127-145. tools:::Rd_expr_doi("10.1111/j.2517-6161.1987.tb01685.x")
Smyth, G.K. (1989), Generalized Linear Models with Varying Dispersion. Journal of the Royal Statistical Society: Series B (Methodological), 51: 47-60. tools:::Rd_expr_doi("10.1111/j.2517-6161.1989.tb01747.x")
stfamily, glmstarma.control, dglmstarma, TimeConstant, SpatialConstant
# \donttest{
dat <- load_data("chickenpox", directory = tempdir())
chickenpox <- dat$chickenpox
population_hungary <- dat$population_hungary
W_hungary <- dat$W_hungary
model_autoregressive <- list(past_obs = rep(1, 7))
dglmstarma(chickenpox, model_autoregressive, dispersion_model = list(past_obs = 1),
mean_covariates = list(population = population_hungary),
wlist = W_hungary, mean_family = vquasipoisson("log"))
# }
Run the code above in your browser using DataLab