inner-vs-outer: Inner iterations versus outer optimization

In spaMM: Mixed-Effect Models, with or without Spatial Random Effects

Inner iterations versus outer optimization

Description

For estimation of phi and lambda, two approaches are implemented: generic optimization functions to maximize the log-likelihood (or restricted log-likelihood in REML case), and iterative algorithms described by Lee et al. (2006). This distinction is extended to random-coefficient parameters (ranCoefs) and to the correlation parameter of adjacency models in spaMM, but other correlation parameters of random effects are estimated only by generic optimization. Fixed-effect coefficients (\beta) are always estimated by iterative methods. The iterative method is well suited for estimating parameters of heteroscedastic (or “structured-dispersion”) models, but it then requires computations which can be slow for large datasets.

The HLfit and HLCor functions only use the iterative algorithms. Other fitting functions can call both approaches, internally using HLfit and HLCor to perform any iterative estimation. For this reason, all fitting functions, except HLfit and HLCor, are called “outer” optimization functions in this documentation, possibly calling “inner” iterative algorithms. Likewise, parameters may be “outer”- or “inner”-estimated.

fitme by default tries to select the fastest method when both can be applied. Therefore, fitme tends to use generic generic optimization, rather than iterative algorithms, for large datasets, unless there is a non-trivial resid.model. corrHLfit tends to use iterative methods as much as possible.

Compared to generic optimization, the iterative algorithms provide more information, such as the “cond. SE” of the dispersion parameter estimates (about which see warning in Details of HLfit). To force the evaluation of such information after an outer-optimization by an outer-optimizing function, use the control$refit argument (see Examples in fitme documentation). 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 further Example using ‘blackcap’ data). The same syntax may be used for phi.

Details

The results of REML fits of non-gaussian mixed models with both inner- and outer-estimated parameters may (generally slightly) differ according to which parameters are inner-estimated, when the fits include both some dispersion estimate(s) (say \hat{\lambda}), which maximize restricted likelihood for given fixed effects, and the fixed effects estimates (\hat{\beta}, which maximize marginal likelihood for given \hat{\lambda}). Fit results then depend whether \lambda is inner- or outer-estimated, if \hat{\lambda} depends on \hat{\beta} and conversely. The latter condition may have been ignored on the grounds of asymptotic arguments, but is nevertheless observed.

Explicit values of the init argument of the fitting functions may be used to control the use of inner iterations, and may be used to ensure better consistency of results over time if selection of methods changes changes over successive versions of spaMM: use NaN as noted above to enforce inner-optimization, and NA or a numeric value to enforce outer-optimization.

References

Lee, Y., Nelder, J. A. and Pawitan, Y. (2006). Generalized linear models with random effects: unified analysis via h-likelihood. Chapman & Hall: London.

Examples

###  Contrasting different optimization methods:

## Toy data:
set.seed(123)
gr <- data.frame(y=rgamma(100,shape=9/2,scale=2/3))
# => Gamma deviates with mean mu=3 and variance=2, 
#  ie. phi= var/mu^2= 2/9 

## Fits

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

spaMM documentation built on Sept. 10, 2026, 1:07 a.m.