Learn R Programming

spaMM (version 4.5.0)

fixed: Fixing some parameters

Description

The fitting functions allow all parameters to be fixed rather than estimated:
* Fixed-effect coefficients can be set by way of the etaFix argument (linear predictor coefficients) for all fitting functions.
* Random-effect parameters and the phi parameter of the gaussian and Gamma response families can be set for all fitting function by the fixed argument, or for some fitting functions by an alternative argument with the same effect (see Details for this confusing feature, but using fixed uniformly is simpler).
* The ad-hoc dispersion parameter of some response families (COMPoisson, negbin1, negbin2, beta_resp, betabin and possibly future ones) can be fixed using the ad-hoc argument of such families rather than by fixed.

Arguments

Details

etaFix is a list with single documented element beta, which should be a vector of (a subset of) the coefficients (\(\beta\)) of the fixed effects, with names as shown in a fit without such given values. If REML is used to fit random effect parameters, then etaFix affects by default the REML correction for estimation of dispersion parameters, which depends only on which \(\beta\) coefficients are estimated rather than given. This default behaviour will be overridden whenever a non-null REMLformula is provided to the fitting functions (see Example). Alternatively, with a non-NULL etaFix$beta, REML can also be performed as if all \(\beta\) coefficients were estimated, by adding attribute keepInREML=TRUE to etaFix$beta; in that case the REML computation is by default that implied by the fixed effects in the full model formula, unless a non-default REMLformula is also used.

The older equivalent for the fixed argument is ranFix for HLfit and corrHLfit, and ranPars for HLCor. Do not use both one such argument and fixed in a call. This older diversity of names was confusing, but its logic was that ranFix allows one to fix parameters that HLfit and corrHLfit would otherwise estimate, while ranPars can be used to set correlation parameters that HLCor does not estimate but nevertheless requires (e.g., Matérn parameters).

Theses arguments for fixing random-effect parameters all have a common syntax. They is a list, with the following possible elements, whose nature is further detailed below:
* phi (variance of residual error, for gaussian and Gamma HGLMs),
* lambda (random-effect variances, except for random-coefficient terms),
* ranCoefs (random-coefficient parameters),
* corrPars (correlation parameters, when handled by the fitting function).
* Individual correlation parameters such as rho, nu, Nugget, ARphi... are also possible top-level elements of the list when there is no ambiguity as to which random effect these correlation parameters apply. This syntax was conceived when spaMM handled a single spatial random effect, and it is still convenient when applicable, but it should not be mixed with corrPars element usage.

phi may be a single value or a vector of the same length as the response vector (the number of rows in the data, once non-informative rows are removed).

lambda may be a single value (if there is a single random effect, or a vector allowing to specify unambiguously variance values for some random effect(s). It can thus take the form lambda=c(NA,1) or lambda=c("2"=1) (note the name) to assign a value only to the variance of the second of two random effects.

ranCoefs is a list of numeric vectors, each numeric vector specifying the variance and correlation parameters for a random-coefficient term. As for lambda, it may be incomplete, using names to specify the random effect to which the parameters apply. For example, to assign variances values 3 and 7, and correlation value -0.05, to the second random effect in a model formula, one can use ranCoefs=list("2"=c(3,-0.05,7)) (note the name). The elements of each vector are variances and correlations, matching those of the printed summary of a fit. The order of these elements must be the order of the lower.tri of a covariance matrix, as shown e.g. by
m2 <- matrix(NA, ncol=2,nrow=2); m2[lower.tri(m2,diag=TRUE)] <- seq(3); m2.
fitme accepts partially fixed parameters for a random coefficient term, e.g.,
ranCoefs=list("2"=c(NA,-0.05,NA)), although this may not mix well with some obscure options, such as
control=list(refit=list(ranCoefs=TRUE)) which will ignore the fixed values. GxE shows how to use partially-fixed ranCoefs to fit different variances for different levels of a factor.

corrPars is a list, and it may also be incomplete, using names to specify the affected random effect as shown for lambda and ranCoefs. For example, ranFix=list(corrPars=list("1"=list(nu=0.5))) makes explicit that nu=0.5 applies to the first ("1") random effect in the model formula. Its elements may be the correlation parameters of the given random effect. For the Matérn model, these are the correlation parameters rho (scale parameter(s)), nu (smoothness parameter), and (optionally) Nugget (see Matern). The rho parameter can itself be a vector with different values for different geographic coordinates. For the adjacency model, the only correlation parameter is a scalar rho (see adjacency). For the AR1 model, the only correlation parameter is a scalar ARphi (see AR1). Consult the documentation for other types of random effects, such as Cauchy or IMRF, for any information missing here.

Examples

Run this code
if (FALSE) {
data("wafers")
# Fixing random-coefficient parameters:
fitme(y~X1+(X2|batch), data=wafers, fixed=list(ranCoefs=list("1"=c(2760, -0.1, 1844))))
##  HLfit syntax for the same effect (except that REML is used here)
# HLfit(y~X1+(X2|batch), data=wafers, ranFix=list(ranCoefs=list("1"=c(2760, -0.1, 1844))))


### Fixing coefficients of the linear predictor:
#
## ML fit
#
fitme(y ~X1+X2+X1*X3+X2*X3+I(X2^2)+(1|batch), data=wafers, family=Gamma(log), 
      etaFix=list(beta=c("(Intercept)"=5.61208)))
#      
## REML fit
# Evaluation of restricted likelihood depends on which fixed effects are estimated,
# so simply fixing the coefficients to their REML estimates will not yield 
# the same REML fits, as see by comparing the next two fits:
#
unconstr <- fitme(y ~X1+X2+X1*X3+X2*X3+I(X2^2)+(1|batch), data=wafers, 
                  family=Gamma(log), method="REML")
#  
# Second fit is different from 'unconstr' despite the same fixed-effects:     
naive <- fitme(y ~X1+X2+X1*X3+X2*X3+I(X2^2)+(1|batch), data=wafers, family=Gamma(log), 
               method="REML", etaFix=list(beta=fixef(unconstr)))    
#
# Using REMLformula to obtain the same REML fit as the unconstrained one:
fitme(y ~X1+X2+X1*X3+X2*X3+I(X2^2)+(1|batch), data=wafers, family=Gamma(log), 
      method="REML", etaFix=list(beta=fixef(unconstr)),
      REMLformula=y ~X1+X2+X1*X3+X2*X3+I(X2^2)+(1|batch))

data("Loaloa")
# Fixing some Matern correlation parameters, in fitme():
fitme(cbind(npos,ntot-npos) ~ elev1 +Matern(1|longitude+latitude),
             data=Loaloa,family=binomial(),fixed=list(nu=0.5,Nugget=2/7))
# Fixing all mandatory Matern correlation parameters, in HLCor():
HLCor(cbind(npos,ntot-npos) ~ elev1 + Matern(1|longitude+latitude),
             data=Loaloa,family=binomial(),ranPars=list(nu=0.5,rho=0.7))
}

Run the code above in your browser using DataLab