spaMM (version 3.8.23)

spaMM: Inference in mixed models, in particular spatial GLMMs

Description

Fits a range of mixed-effect models, including those with spatially correlated random effects. The random effects are either Gaussian (which defines GLMMs), or other distributions (which defines the wider class of hierarchical GLMs), or simply absent (which makes a LM or GLM). spaMM is designed to be used through the high-level fitting functions fitme (the most general function), HLfit (sometimes faster, for non-spatial models), HLCor (sometimes faster, for conditional-autoregressive models and fixed-correlation models), corrHLfit (now of lesser interest); and additional functions such as fixedLRT for likelihood-ratio testing, simulate and predict.

Both maximum likelihood (ML) and restricted likelihood (REML) can be used for linear mixed models, and extensions of these methods using Laplace approximations are used for non-Gaussian random response. Several variants of these methods discussed in the literature are included (see Details in HLfit), the most notable of which may be “PQL/L” for binary-response GLMMs (see Example for arabidopsis data). PQL methods implemented in spaMM are closer to (RE)ML methods than those implemented in MASS::glmmPQL.

Details

The standard response families gaussian, binomial, poisson, and Gamma are handled, as well as negative binomial (see negbin), zero-truncated poisson and negative binomial, and Conway-Maxwell-Poisson response (see Tpoisson, Tnegbin and COMPoisson). A multi family look-alike is also available for multinomial response, with some constraints.

The variance parameter of residual error is denoted \(\phi\) (phi): this is the residual variance for gaussian response, but for Gamma-distributed response, the residual variance is \(\phi\)\(\mu^2\) where \(\mu\) is expected response. A fixed-effects linear predictor for \(\phi\), modeling heteroscedasticity, can be considered (see Examples).

The package fits models including several nested or crossed random effects, including autocorrelated ones with the following correlation models: Matern, Cauchy, interpolated Markov Random Fields (IMRF), first-order autoregressive (AR1), conditional autoregressive as specified by an adjacency matrix, or any fixed correlation matrix (corrMatrix). GLMMs and HGLMs are fit via Laplace approximations for (1) the marginal likelihood with respect to random effects and (2) the restricted likelihood (as in REML), i.e. the likelihood of random effect parameters given the fixed effect estimates.

All handled models can be formulated in terms of a linear predictor of the traditional form offset+ X\(\beta\) + Z b, where X is the design matrix of fixed effects, \(\beta\) (beta) is a vector of fixed-effect coefficients, Z is a “design matrix” for the random effects (which is instead denoted M=ZAL elsewhere in the package documentation), and b a vector of random effect values. The general structure of Mb is described in random-effects.

The syntax for formulas extends that used in the lme4 package. In particular, non-autocorrelated random effects are specified using the (1|<block>) syntax, and Gaussian random-coefficient terms by (<lhs>|<block>). The double-vertical syntax, ( lhs || rhs ), is interpreted as in lme4 (see help('expandDoubleVerts') in lme4, including the comment that it will fail to break up factors into separate terms). Any such term is immediately converted to ( ( 1 | rhs ) + ( 0 + lhs | rhs ) ), and should be counted as two random effects for all purposes (e.g., for fixing the variances of the random effects). Autocorrelated random effects are specified by adding some prefix to this syntax, <prefix>(1|.), e.g. Matern(1|long+lat).

If one wishes to fit uncorrelated group-specific random-effects (say random effects with distinct variances for males and females or for different response variables), none of these syntaxes appears useful. In particular, using factors or boolean in the lhs of the term generally results in a correlated random-coefficient term being fitted, even when “0 + ”... is used. Instead, one should define numeric (0|1) variables such as numF <- as.numeric(<boolean for females>) and a term of the form (0 + numF| <block>). See “Elementary multivariate-response” example below.

Since version 2.6.0, it is possible to fit some “some autocorrelated random-coefficient” models by a syntax consistent with that of random-coefficient terms, <prefix>(<rhs>|.). For example, independent Mat<U+00E9>rn effects can be fitted for males and females by using the syntax Matern(male|.) + Matern(female|.), where male and female are TRUE/FALSE factors. A numeric variable z can also be considered, in which case the proper syntax is Matern(0+z|.), which represents an autocorrelated random-slope (only) term (or, equivalently, a direct specification of heteroscedasticy of the Mat<U+00E9>rn random effect). Thus, for boolean z, Matern(z|.) is de facto automatically converted to Matern(0+<as.numeric(z)>|.)). These effects are achieved by direct control of the elements of the incidence matrix Z through the <rhs> term: for numeric z, such elements are multiplied by z values, and thus provide a variance of order O(z squared).

By contrast, Matern(z|.) is not defined for numeric z. It could mean that a correlation structure between random intercepts and random slopes is to be combined with a Mat<U+00E9>rn correlation structure, but no way of combining them is yet defined and implemented (though this may change in some later version).

Since version 2.7.0, the syntax (z-1|.), for numeric z only, can also be used to fit some heteroscedastic non-Gaussian random effects. For example, a Gamma random-effect term (wei-1|block) specifies an heteroscedastic Gamma random effect \(u\) with constant mean 1 and variance wei^2 \(\lambda\), where \(\lambda\) is still the estimated variance parameter. See Details of negbin for a possible application. Here, this effect is not implemented through direct control of Z (multiplying the elements of an incidence matrix Z by wei), as this would have a different effect on the distribution of the random effect term. (z|.) is not defined for non-Gaussian random effects. It could mean that a correlation structure between random intercepts and random slopes for (say) Gamma-distributed random effects is considered, but such correlation structures are not well-specified by their correlation matrix.

References

Lee, Y., Nelder, J. A. and Pawitan, Y. (2006). Generalized linear models with random effects: unified analysis via h-likelihood. Chapman & Hall: London.

Rousset F., Ferdy, J.-B. (2014) Testing environmental and genetic effects in the presence of spatial autocorrelation. Ecography, 37: 781-790. 10.1111/ecog.00566

See Also

The test directory of the package provides many additional examples of spaMM usage beyond those from the formal documentation.

Examples

# NOT RUN {
data("wafers")
data("scotlip") ## loads 'scotlip' data frame, but also 'Nmatrix'

##     Linear model
fitme(y ~ X1, data=wafers)

##     GLM
fitme(y ~ X1, family=Gamma(log), data=wafers)
fitme(cases ~ I(log(population)), data=scotlip, family=poisson)

##     Non-spatial GLMMs
fitme(y ~ 1+(1|batch), family=Gamma(log), data=wafers)
fitme(cases ~ 1+(1|gridcode), data=scotlip, family=poisson)
#
# Random-slope model (mind the output!)        
fitme(y~X1+(X2|batch),data=wafers, method="REML")

## Spatial, conditional-autoregressive GLMM
if (spaMM.getOption("example_maxtime")>2) {   
  fitme(cases ~ I(log(population))+adjacency(1|gridcode), data=scotlip, family=poisson, 
        adjMatrix=Nmatrix) # with adjacency matrix provided by data("scotlip")
} 
# see ?adjacency for more details on these models 

## Spatial, geostatistical GLMM: 
# see e.g. examples in ?fitme, ?corrHLfit, ?Loaloa, or ?arabidopsis;
# see examples in ?Matern for group-specific spatial effects.

##     Hierachical GLMs with non-gaussian random effects
 data("salamander")
if (spaMM.getOption("example_maxtime")>1) {   
 # both gaussian and non-gaussian random effects
 fitme(cbind(Mate,1-Mate)~1+(1|Female)+(1|Male),family=binomial(),
        rand.family=list(gaussian(),Beta(logit)),data=salamander)
 
 # Random effect of Male nested in that of Female:
 fitme(cbind(Mate,1-Mate)~1+(1|Female/Male),
       family=binomial(),rand.family=Beta(logit),data=salamander)
 # [ also allowed is cbind(Mate,1-Mate)~1+(1|Female)+(1|Male %in% Female) ]
}

##    Modelling residual variance ( = structured-dispersion models)    
# GLM response, fixed effects for residual variance 
fitme( y ~ 1,family=Gamma(log),
      resid.model = ~ X3+I(X3^2) ,data=wafers)
#
# GLMM response, and mixed effects for residual variance
if (spaMM.getOption("example_maxtime")>1.5) {   
  fitme(y ~ 1+(1|batch),family=Gamma(log),
        resid.model = ~ 1+(1|batch) ,data=wafers)
}

## Elementary multivariate-response model (see fitmv() for less elementary ones) 
# Data preparation
fam <- rep(c(1,2),rep(6,2)) # define two biological 'families'
ID <- gl(6,2) # define 6 'individuals'
resp <- as.factor(rep(c("x","y"),6)) # distinguishes two responses per individual
set.seed(123)
toymv <- data.frame(
  fam = as.factor(fam), ID = ID, resp = resp, 
  y = 1 + (resp=="x") + rnorm(4)[2*(resp=="x")+fam] + rnorm(12)[6*(resp=="x")+as.integer(ID)],
  respX = as.numeric(resp=="x"),
  respY = as.numeric(resp=="y")
)
#
# fit response-specific variances of random effect and residuals: 
(fitme(y ~ resp+ (0+respX|fam)+ (0+respY|fam), resid.model = ~ 0+resp ,data=toymv))

# }