spaMM (version 1.4.1)

spaMM: Inference in spatial GLMMs (and other non-spatial mixed models)

Description

spaMM has been first designed to fit spatial mixed models. It can also fit a wider class of non-spatial mixed models. 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 (G)LMM response families gaussian, binomial, poisson, and Gamma are handled, as well as multinomial response (with some constraints) and negative binomial response (see Examples). 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 package fits spatial models following either a Matérn correlation model or an adjacency matrix model. It also handles a number of non-spatial models. It accepts several nested or crossed random effects (see Examples). Only one of these can be a spatial random effect. The spaMM code has not been optimized for speed in these models. The variance(s) of random effects ($u$) is (are) denoted $\lambda$ (lambda in input and output) and that of residual error is denoted $\phi$ (phi). A fixed-effects linear predictor for $\phi$, modeling heteroskedasticity, 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 Predictor 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). The package has been extensively tested mainly for analysis of spatial GLMMs (Rousset and Ferdy 2014), where the random effects are Gaussian. Other models have been checked against literature results and a few simulations.

References

Lee, Y., Nelder, J. A. and Pawitan, Y. (2006). Generalised 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

## Fit a Poisson GLMM with adjacency correlation model

data(scotlip) ## loads 'scotlip' data frame, but also 'Nmatrix'
corrHLfit(cases~I(prop.ag/10) +adjacency(1|gridcode)+offset(log(scotlip$expec)),
          data=scotlip,family=poisson(),
          adjMatrix=Nmatrix,lower=list(rho=0),upper=list(rho=0.1745))
          
## Adding a Gamma random effect to fit a negative-binomial response:
corrHLfit(cases~I(prop.ag/10) +(1|gridcode)+adjacency(1|gridcode)
                +offset(log(scotlip$expec)),
          data=scotlip,family=poisson(),rand.family=list(Gamma(log),gaussian()),
          adjMatrix=Nmatrix,lower=list(rho=0),upper=list(rho=0.1745))

## 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")

## Nested effects
# 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")

# A syntax described in ?formula
HLfit(cbind(Mate,1-Mate)~1+(1|Female)+(1|Male %in% Female),
       family=binomial(),rand.family=Beta(logit),data=salamander,HLmethod="ML")

## 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.formula = ~ X3+I(X3^2) ,data=wafers)

## Random-slope model (mind the output!)        
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.formula = ~ X3+I(X3^2) ,data=wafers)

## Fit of binary data using PQL/L. See ?arabidopsis
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)