spaMM (version 2.6.1)

spaMM: Inference in mixed models, in particular spatial GLMMs

Description

Fits a range of mixed 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 GLM).

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 package fits models including several nested or crossed random effects, some of which can be following a Mat<U+00E9>rn correlation model (see Matern), or a Cauchy correlation model (see CauchyCorr), or an adjacency matrix model (see adjacency), or an AR1 model, or a gven 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.

The variance(s) of random effects (\(u\)) is (are) denoted \(\lambda\) (lambda in input and output). 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). Fixed effects are described in the standard form X\(\beta\) where X is the design matrix of fixed effects and \(\beta\) (beta) is a vector of fixed effect parameters.

The structure of the random effects can generally be described by the following steps. First, independent and identically distributed (iid) random effects \(u\) are drawn from one of the following distributions: gaussian, Beta-distributed, Gamma and inverse-Gamma distributed random effects, implemented as detailed in the HLfit documentation. Second, a transformation v\(=f\)(u) is applied (v elements are still iid). Third, correlated random effects are obtained as Mv, where the matrix M can describe spatial correlation between observed locations, block effects (or repeated observations in given locations), and correlations involving unobserved locations. See Details in covStruct for the general form of M. In most cases M is determined from the model formula, but it can also be input directly (e.g., to describe genetic correlations).

Since version 2.6.0, it is possible to fit “autocorrelated random-coefficient” models by a syntax analogous to that of other random-coefficient terms, providing direct control of the elements of the Z matrix. 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). By contrast, Matern(z|.) is not defined. 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.

The double-vertical syntax for random effects, ( rhs || lhs ), is handled as in lme4. Any such term is immediately converted to ( ( 1 | lhs ) + ( 0 + lhs | rhs ) ), and should be counted as two random effects for all purposes (e.g., for fixing the variances of the random effects).

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. http://dx.doi.org/10.1111/ecog.00566

See Also

spaMM is designed to be used through the high-level functions corrHLfit, HLCor, HLfit, and fixedLRT

Examples

# NOT RUN {
## Fit a Poisson GLMM with adjacency (CAR) correlation model
# see ?adjacency for how to fit efficiently such model models 
data("scotlip") ## loads 'scotlip' data frame, but also 'Nmatrix'
HLCor(cases~I(prop.ag/10) +adjacency(1|gridcode)+offset(log(expec)),
          adjMatrix=Nmatrix,family=poisson(),data=scotlip) 


if (spaMM.getOption("example_maxtime")>2.1) {          
 ## Adding a Gamma random effect to fit a negative-binomial response:
 HLCor(cases~I(prop.ag/10) +(1|gridcode)+adjacency(1|gridcode)
                +offset(log(expec)),
          data=scotlip,family=poisson(),rand.family=list(Gamma(log),gaussian()),
          adjMatrix=Nmatrix)
}

 
# }
# NOT RUN {
<!-- % tested in test-spaMM.R -->
# }
# NOT RUN {
 ## fit non-spatial crossed random effects with distinct families
 data("salamander")
 HLfit(cbind(Mate,1-Mate)~1+(1|Female)+(1|Male),family=binomial(),
        rand.family=list(gaussian(),Beta(logit)),data=salamander,HLmethod="ML")
# }
# NOT RUN {
 ## Nested effects
 
# }
# NOT RUN {
<!-- % tested in test-spaMM.R -->
# }
# NOT RUN {
 # lmer syntax allowing several degrees of nesting
 HLfit(cbind(Mate,1-Mate)~1+(1|Female/Male),
       family=binomial(),rand.family=Beta(logit),data=salamander,HLmethod="ML")
 # [ also allowed is cbind(Mate,1-Mate)~1+(1|Female)+(1|Male %in% Female) ]
# }
# NOT RUN {
## fit a non-spatial, Gamma GLMM:
data("wafers")
HLfit(y ~X1*X3+X2*X3+I(X2^2)+(1|batch),family=Gamma(log),
          data=wafers)

## Same with fixed-effects predictor for residual variance 
## ( = structured-dispersion model):
HLfit(y ~X1*X3+X2*X3+I(X2^2)+(1|batch),family=Gamma(log),
          resid.model = ~ X3+I(X3^2) ,data=wafers)

## Random-slope model (mind the output!)        
if (spaMM.getOption("example_maxtime")>1) {          
  HLfit(y~X1+(X2|batch),data=wafers)
}

## fit a GLM (not mixed) with structured dispersion:
HLfit( y ~X1+X2+X1*X3+X2*X3+I(X2^2),family=Gamma(log),
      resid.model = ~ X3+I(X3^2) ,data=wafers)

## Fit of binary data using PQL/L. See ?arabidopsis
# }
# NOT RUN {
data("arabidopsis")
HLCor(cbind(pos1046738,1-pos1046738)~seasonal+Matern(1|LAT+LONG),
                   ranPars=list(rho=0.129,lambda=4.28,nu=0.291),
                   family=binomial(),HLmethod="PQL/L",data=arabidopsis)
# }