Fitting function for fixed- and mixed-effect models with GLM response.
This is a common interface for fitting most models that spaMM can fit, from linear models to mixed models with non-gaussian random effects, therefore substituting to corrHLfit, HLCor and HLfit. By default, it uses ML rather than REML (differing in this respect from the other fitting functions). It may use “outer optimization”, i.e., generic optimization methods for estimating all dispersion parameters, rather than the iterative methods implemented in HLfit. The results of REML fits of non-gaussian mixed models by these different methods may (generally slightly) differ. Outer optimization should generally be faster than the alternative algorithms for large data sets when the residual variance model is a single constant term (no structured dispersion). For mixed models, fitme by default tries to select the fastest method when both can be applied, but precise decision criteria are subject to change in the future. corrHLfit (with non-default arguments to control the optimization method most suitable to a particular problem) may be used to ensure better consistency over successive versions of spaMM.
fitme(formula, data, family = gaussian(), init = list(), fixed = list(),
lower = list(), upper = list(), resid.model = ~1, init.HLfit = list(),
control = list(), control.dist = list(), method = "ML",
HLmethod = method, processed = NULL, nb_cores = NULL, objective = NULL,
...)formula |
Either a linear model |
data |
A data frame containing the variables in the response and the model formula. |
family |
|
init |
An optional list of initial values for correlation and/or dispersion parameters and/or response family parameters, e.g.
|
fixed |
A list similar to |
lower |
An optional (sub)list of values of the parameters specified through |
upper |
Same as |
resid.model |
See identically named |
init.HLfit |
See identically named |
control |
A list of control parameters, with two possible elements:
|
control.dist |
See |
method, HLmethod |
Character: the fitting method to be used, such as |
nb_cores |
Not yet operative, only for development purposes. Number of cores to use for parallel computations. |
processed |
For programming purposes, not documented. |
objective |
For development purpose, not documented. |
... |
Optional arguments passed to (or operating as if passed to) |
For approximations of likelihood, see method. For the possible structures of random effects, see random-effects,
For phi, lambda, and ranCoefs, fitme may or may not use the internal fitting methods of HLfit. The latter methods are well suited for structured dispersion models, but require computations which can be slow for large datasets. Therefore, fitme tends to outer-optimize by default for large datasets, unless there is a non-trivial resid.model. The precise criteria for selection of default method by fitme are liable to future changes.
Further, the internal fitting methods of HLfit also provide some more information such as the “cond. SE” (about which see warning in Details of HLfit). To force the evaluation of such information after an outer-optimization by a fitme call, use the control$refit argument (see Example). Alternatively (and possibly of limited use), one can force inner-optimization of lambda for a given random effect, or of phi, by setting it to NaN in init (see Example using ‘blackcap’ data). The same syntax may be tried for phi.
The return value of an HLCor or an HLfit call, with additional attributes. The HLCor call is evaluated at the estimated correlation parameter values. These values are included in the return object as its $corrPars member. The attributes added by fitme include the original call of the function (which can be retrived by getCall(<fitted object>), and information about the optimization call within fitme.
## Examples with Matern correlations
## A likelihood ratio test based on the ML fits of a full and of a null model.
data("blackcap")
(fullfit <- fitme(migStatus ~ means+ Matern(1|longitude+latitude),data=blackcap) )
(nullfit <- fitme(migStatus ~ 1 + Matern(1|longitude+latitude),data=blackcap))
## p-value:
1-pchisq(2*(logLik(fullfit)-logLik(nullfit)),df=1)
## See ?spaMM for examples of conditional autoregressive model and of non-spatial models.
## Contrasting different optimization methods:
# We simulate Gamma deviates with mean mu=3 and variance=2,
# ie. phi= var/mu^2= 2/9 in the (mu, phi) parametrization of a Gamma
# GLM; and shape=9/2, scale=2/3 in the parametrisation of rgamma().
# Note that phi is not equivalent to scale:
# shape = 1/phi and scale = mu*phi.
set.seed(123)
gr <- data.frame(y=rgamma(100,shape=9/2,scale=2/3))
# Here fitme uses HLfit methods which provide cond. SE for phi by default:
fitme(y~1,data=gr,family=Gamma(log))
# To force outer optimization of phi, use the init argument:
fitme(y~1,data=gr,family=Gamma(log),init=list(phi=1))
# To obtain cond. SE for phi after outer optimization, use the 'refit' control:
fitme(y~1,data=gr,family=Gamma(log),,init=list(phi=1),
control=list(refit=list(phi=TRUE))) ## or ...refit=TRUE...
## Outer-optimization is not necessarily the best way to find a global maximum,
# particularly when there is little statistical information in the data:
if (spaMM.getOption("example_maxtime")>1.6) {
data("blackcap")
fitme(migStatus ~ means+ Matern(1|longitude+latitude),data=blackcap) # poor
# Compare with the following two ways of avoiding outer-optimization of lambda:
corrHLfit(migStatus ~ means+ Matern(1|longitude+latitude),data=blackcap,
method="ML")
fitme(migStatus ~ means+ Matern(1|longitude+latitude),data=blackcap,
init=list(lambda=NaN))
}
## see help("COMPoisson"), help("negbin"), help("Loaloa"), etc., for further examples.Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.