Fit spatial generalized autoregressive models

Fit spatial generalized linear models for areal data (i.e., spatial generalized autoregressive models) using a variety of estimation methods, allowing for random effects, partition factors, and row standardization.

Usage

spgautor(
  formula,
  family,
  data,
  spcov_type,
  spcov_initial,
  dispersion_initial,
  estmethod = "reml",
  random,
  randcov_initial,
  partition_factor,
  W,
  row_st = TRUE,
  M,
  range_positive = TRUE,
  ...
)

Arguments

formula: A two-sided linear formula describing the fixed effect structure of the model, with the response to the left of the ~ operator and the terms, separated by + operators, on the right.
family: The generalized linear model family describing the distribution of the response variable to be used. Available options "poisson", "nbinomial", "binomial", "beta", "Gamma", and "inverse.gaussian". Can be quoted or unquoted. Note that the family argument only takes a single value, rather than the list structure used by stats::glm. See Details for more.
data: A data frame or sf object that contains the variables in fixed, random, and partition_factor, as well as potentially geographical information.
spcov_type: The spatial covariance type. Available options include "car" and "sar". Parameterizations of each spatial covariance type are available in Details. When spcov_type is specified, relevant spatial covariance parameters are assumed unknown, requiring estimation. spcov_type is not required (and is ignored) if spcov_initial is provided. Multiple values can be provided in a character vector. Then spgautor() is called iteratively for each element and a list is returned for each model fit. The default for spcov_type is "car".
spcov_initial: An object from spcov_initial() specifying initial and/or known values for the spatial covariance parameters. Not required if spcov_type is provided. Multiple spcov_initial() objects can be provided in a list. Then spgautor() is called iteratively for each element and a list is returned for each model fit.
dispersion_initial: An object from dispersion_initial() specifying initial and/or known values for the dispersion parameter for the "nbinomial", "beta", "Gamma", and "inverse.gaussian" families. family is ignored if dispersion_initial is provided.
estmethod: The estimation method. Available options include "reml" for restricted maximum likelihood and "ml" for maximum likelihood The default is "reml".
random: A one-sided linear formula describing the random effect structure of the model. Terms are specified to the right of the ~ operator. Each term has the structure x1 + ... + xn | g1/.../gm, where x1 + ... + xn specifies the model for the random effects and g1/.../gm is the grouping structure. Separate terms are separated by + and must generally be wrapped in parentheses. Random intercepts are added to each model implicitly when at least one other variable is defined. If a random intercept is not desired, this must be explicitly defined (e.g., x1 + ... + xn - 1 | g1/.../gm). If only a random intercept is desired for a grouping structure, the random intercept must be specified as 1 | g1/.../gm. Note that g1/.../gm is shorthand for (1 | g1/.../gm). If only random intercepts are desired and the shorthand notation is used, parentheses can be omitted.
randcov_initial: An optional object specifying initial and/or known values for the random effect variances.
partition_factor: A one-sided linear formula with a single term specifying the partition factor. The partition factor assumes observations from different levels of the partition factor are uncorrelated.
W: Weight matrix specifying the neighboring structure used. Not required if data is an sf polygon object, as W is calculated internally using queen contiguity. If calculated internally, W is computed using sf::st_intersects().
row_st: A logical indicating whether row standardization be performed on W. The default is TRUE.
M: M matrix satisfying the car symmetry condition. The car symmetry condition states that \((I - range * W)^{-1}M\) is symmetric, where \(I\) is an identity matrix, \(range\) is a constant that controls the spatial dependence, W is the weights matrix, and \(^{-1}\) represents the inverse operator. M is required for car models when W is provided and row_st is FALSE. When M, is required, the default is the identity matrix. M must be diagonal or given as a vector or one-column matrix assumed to be the diagonal.
range_positive: Whether the range should be constrained to be positive. The default is TRUE.
...: Other arguments to stats::optim().

Value

A list with many elements that store information about the fitted model object. If spcov_type or spcov_initial are length one, the list has class spgautor. Many generic functions that summarize model fit are available for spgautor objects, including AIC, AICc, anova, augment, coef, cooks.distance, covmatrix, deviance, fitted, formula, glance, glances, hatvalues, influence, labels, logLik, loocv, model.frame, model.matrix, plot, predict, print, pseudoR2, summary, terms, tidy, update, varcomp, and vcov. If spcov_type or spcov_initial are length greater than one, the list has class spgautor_list and each element in the list has class spgautor. glances can be used to summarize spgautor_list

objects, and the aforementioned spgautor generics can be used on each individual list element (model fit).

Details

The spatial generalized linear model for areal data (i.e., spatial generalized autoregressive model) can be written as \(g(\mu) = \eta = X \beta + \tau + \epsilon\), where \(\mu\) is the expectation of the response (\(y\)) given the random errors, \(g(.)\) is called a link function which links together the \(\mu\) and \(\eta\), \(X\) is the fixed effects design matrix, \(\beta\) are the fixed effects, \(\tau\) is random error that is spatially dependent, and \(\epsilon\) is random error that is spatially independent.

There are six generalized linear model families available: poisson assumes \(y\) is a Poisson random variable nbinomial assumes \(y\) is a negative binomial random variable, binomial assumes \(y\) is a binomial random variable, beta assumes \(y\) is a beta random variable, Gamma assumes \(y\) is a gamma random variable, and inverse.gaussian assumes \(y\) is an inverse Gaussian random variable.

The supports for \(y\) for each family are given below:

family: support of \(y\)
poisson: \(0 \le y\); \(y\) an integer
nbinomial: \(0 \le y\); \(y\) an integer
binomial: \(0 \le y\); \(y\) an integer
beta: \(0 < y < 1\)
Gamma: \(0 < y\)
inverse.gaussian: \(0 < y\)

The generalized linear model families and the parameterizations of their link functions are given below:

family: link function
poisson: \(g(\mu) = log(\eta)\) (log link)
nbinomial: \(g(\mu) = log(\eta)\) (log link)
binomial: \(g(\mu) = log(\eta / (1 - \eta))\) (logit link)
beta: \(g(\mu) = log(\eta / (1 - \eta))\) (logit link)
Gamma: \(g(\mu) = log(\eta)\) (log link)
inverse.gaussian: \(g(\mu) = log(\eta)\) (log link)

The variance function of an individual \(y\) (given \(\mu\)) for each generalized linear model family is given below:

family: \(Var(y)\)
poisson: \(\mu \phi\)
nbinomial: \(\mu + \mu^2 / \phi\)
binomial: \(n \mu (1 - \mu) \phi\)
beta: \(\mu (1 - \mu) / (1 + \phi)\)
Gamma: \(\mu^2 / \phi\)
inverse.gaussian: \(\mu^2 / \phi\)

The parameter \(\phi\) is a dispersion parameter that influences \(Var(y)\). For the poisson and binomial families, \(\phi\) is always one. Note that this inverse Gaussian parameterization is different than a standard inverse Gaussian parameterization, which has variance \(\mu^3 / \lambda\). Setting \(\phi = \lambda / \mu\) yields our parameterization, which is preferred for computational stability. Also note that the dispersion parameter is often defined in the literature as \(V(\mu) \phi\), where \(V(\mu)\) is the variance function of the mean. We do not use this parameterization, which is important to recognize while interpreting dispersion parameter estimates. For more on generalized linear model constructions, see McCullagh and Nelder (1989).

Together, \(\tau\) and \(\epsilon\) are modeled using a spatial covariance function, expressed as \(de * R + ie * I\), where \(de\) is the dependent error variance, \(R\) is a matrix that controls the spatial dependence structure among observations, \(ie\) is the independent error variance, and \(I\) is an identity matrix. Note that \(de\) and \(ie\) must be non-negative while \(range\) must be between the reciprocal of the maximum eigenvalue of W and the reciprocal of the minimum eigenvalue of W. Recall that \(\tau\) and \(\epsilon\) are modeled on the link scale, not the inverse link (response) scale. Random effects are also modeled on the link scale.

spcov_type Details: Parametric forms for \(R\) are given below:

car: \((I - range * W)^{-1}M\), weights matrix \(W\), symmetry condition matrix \(M\)
sar: \([(I - range * W)(I - range * W)^T]^{-1}\), weights matrix \(W\), \(^T\) indicates matrix transpose

If there are observations with no neighbors, they are given a unique variance parameter called extra, which must be non-negative.

estmethod Details: The various estimation methods are

reml: Maximize the restricted log-likelihood.
ml: Maximize the log-likelihood.

Note that the likelihood being optimized is obtained using the Laplace approximation.

By default, all spatial covariance parameters except ie as well as all random effect variance parameters are assumed unknown, requiring estimation. ie is assumed zero and known by default (in contrast to models fit using spglm(), where ie is assumed unknown by default). To change this default behavior, specify spcov_initial (an NA value for ie in spcov_initial to assume ie is unknown, requiring estimation).

random Details: If random effects are used, the model can be written as \(y = X \beta + Z1u1 + ... Zjuj + \tau + \epsilon\), where each Z is a random effects design matrix and each u is a random effect.

partition_factor Details: The partition factor can be represented in matrix form as \(P\), where elements of \(P\) equal one for observations in the same level of the partition factor and zero otherwise. The covariance matrix involving only the spatial and random effects components is then multiplied element-wise (Hadmard product) by \(P\), yielding the final covariance matrix.

Observations with NA response values are removed for model fitting, but their values can be predicted afterwards by running predict(object). This is the only way to perform prediction for spgautor() models (i.e., the prediction locations must be known prior to estimation).

Note

This function does not perform any internal scaling. If optimization is not stable due to large extremely large variances, scale relevant variables so they have variance 1 before optimization.

References

McCullagh P. and Nelder, J. A. (1989) Generalized Linear Models. London: Chapman and Hall.

Examples

spgmod <- spgautor(I(log_trend^2) ~ 1, family = "Gamma", data = seal, spcov_type = "car")
summary(spgmod)
#> 
#> Call:
#> spgautor(formula = I(log_trend^2) ~ 1, family = "Gamma", data = seal, 
#>     spcov_type = "car")
#> 
#> Deviance Residuals:
#>    Min     1Q Median     3Q    Max 
#> -4.484 -2.461 -1.030  0.386  2.823 
#> 
#> Coefficients (fixed):
#>             Estimate Std. Error z value Pr(>|z|)    
#> (Intercept)  -3.7975     0.3396  -11.18   <2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Coefficients (car spatial covariance):
#>       de    range    extra 
#> 0.001738 0.995833 0.002374 
#> 
#> Coefficients (Dispersion for Gamma family):
#> dispersion 
#>     0.3051