Become an expert in R — Interactive courses, Cheat Sheets, certificates and more!
Get Started for Free

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, 
      ...)

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

A list of control parameters, with two possible elements:

  • nloptr, itself a list of control parameters to be copied in the opts argument of nloptr. Default controls are given by spaMM.getOption('nloptr')

  • 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). If $refit is a single boolean, it affects of parameters. By default only lambda is refitted, but this default may change in the future.

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.

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) 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. Non-boolean values of TRACE are meaningful, but the source code of spaMM:::.do_TRACE should be consulted for their meaning.

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.

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.

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.

spaMM

Mixed-Effect Models, with or without Spatial Random Effects

v3.10.0
CeCILL-2
Authors
François Rousset [aut, cre, cph] (<https://orcid.org/0000-0003-4670-0371>), Jean-Baptiste Ferdy [aut, cph], Alexandre Courtiol [aut] (<https://orcid.org/0000-0003-0637-2959>), GSL authors [ctb] (src/gsl_bessel.*)
Initial release
2022-02-06

We don't support your browser anymore

Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.