spaMM (version 4.5.0)

fitme: Fitting function for fixed- and mixed-effect models with GLM response.

Description

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.

Usage

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, 
      weights.form = NULL, ...)

Value

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.

Arguments

formula

Either a linear model formula (as handled by various fitting functions) or a predictor, i.e. a formula with attributes (see Predictor and examples below). See Details in spaMM for allowed terms in the formula.

data

A data frame containing the variables in the response and the model formula.

family

Either a response family or a multi value.

init

An optional list of initial values for correlation and/or dispersion parameters and/or response family parameters, e.g. list(rho=1,nu=1,lambda=1,phi=1) where rho and nu are parameters of the Matérn family (see Matern), and lambda and phi are dispersion parameters (see Details in spaMM for the meaning of these parameters). All are optional, but giving values for a dispersion parameter changes the ways it is estimated (see Details and Examples). rho may be a vector (see make_scaled_dist) and, in that case, it is possible that some or all of its elements are NA, for which fitme substitutes automatically determined values.

fixed

A list similar to init, but specifying fixed values of the parameters not estimated. See fixed for further information; and keep in mind that fixed fixed-effect coefficients can be passed as the etaFix argument as part of the ‘...’.

lower

An optional (sub)list of values of the parameters specified through init, in the same format as init, used as lower values in calls to optim. See Details for default values.

upper

Same as lower, but for upper values.

resid.model

See identically named HLfit argument.

init.HLfit

See identically named HLfit argument.

control.dist

See control.dist in HLCor

method, HLmethod

Character: the fitting method to be used, such as "ML", "REML" or "PQL/L". "ML" is the default, in contrast to "REML" for HLfit, HLCor and corrHLfit. Other possible values of HLfit's method argument are handled. method=c(<"ML" or "REML">,"exp") can be distinctly useful for slow fits of models with Gamma(log) response family (see method).

weights.form

Specification of prior weights by a one-sided formula: use weights.form = ~ pw instead of prior.weights = pw. The effect will be the same except that such an argument, known to evaluate to an object of class "formula", is suitable to enforce safe programming practices (see good-practice).

control

A list of (rarely needed) control parameters, with possible elements:

• refit, a boolean, or a list of booleans with possible elements phi, lambda and ranCoefs. If either element is set to TRUE, then the corresponding parameters are refitted by the internal HLfit methods (see Details), unless these methods were already selected for such parameters in the main fitting step. If refit is a single boolean, it affects all parameters. By default no parameter is refitted.

• optimizer, the numerical optimizer, specified as a string and whose default is controlled by the global spaMM option "optimizer". Possible values are "nloptr", "bobyqa", "L-BFGS-B" and ".safe_opt", whose meanings are detailed in the documentation for the optimizer argument of spaMM.options. Better left unchanged unless suspect fits are obtained.

• nloptr, itself a list of control parameters to be copied in the opts argument of nloptr. Default value is given by spaMM.getOption('nloptr') and possibly other global spaMM options. Better left unchanged unless you are ready to inspect source code.

• bobyqa, optim, lists of controls similar to nloptr but for methods "bobyqa" and "L-BFGS-B", respectively.

nb_cores

For development purpose, not documented.

processed

For programming purpose, not documented.

objective

For development purpose, not documented.

...

Optional arguments passed to (or operating as if passed to) HLCor, HLfit or mat_sqrt, for example rand.family, control.HLfit , verbose or the distMatrix argument of HLCor (so that estimation of Matern or Cauchy parameters can be combined with use of an ad hoc distance matrix). In a fitme call, the verbose vector of booleans may include a TRACE=TRUE element, in which case information is displayed for each set of correlation and dispersion parameter values considered by the optimiser (see verbose for further information, mostly useless except for development purposes).

Details

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.

Examples

## 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.