R: Profile method for merMod objects

profile-methods {lme4}

R Documentation

Profile method for merMod objects

Description

Methods for profile() of [ng]lmer fitted models.

The log() method transforms a lmer profile to the scale of the logarithm of the standard deviation of the random effects, i.e. it returns a profile with all the .sigNN parameters replaced by .lsigNN. The forward and backward splines for these parameters are recalculated.

Usage

## S3 method for class 'merMod'
profile(fitted, which = 1:nptot, alphamax = 0.01,
	maxpts = 100, delta = cutoff/8, verbose = 0, devtol = 1e-09,
        maxmult = 10, startmethod = "prev", optimizer = "bobyqa",
	signames = TRUE, ...)
## S3 method for class 'thpr'
 as.data.frame(x, ...)
## S3 method for class 'thpr'
 log(x, base = exp(1))

Arguments

`fitted`	a fitted model, e.g., the result of `lmer(..)`.
`which`	integer or character vector indicating which parameters to profile: default is all parameters. For integer, i.e., indexing, the parameters are ordered as follows: (1) random effects (theta) parameters; these are ordered as in `getME(.,"theta")`, i.e., as the lower triangle of a matrix with standard deviations on the diagonal and correlations off the diagonal. (2) residual standard deviation (or scale parameter for GLMMs where appropriate); (3) fixed effect (beta) parameters. In addition, `which` may be a character, containing `"beta_"` or `"theta_"` which means the fixed and random effects parameters, respectively, or also containing parameter names, such as `".sigma"` or `"(Intercept)"`.
`alphamax`	a number in (0,1), such that `1 - alphamax` is the maximum alpha value for likelihood ratio confidence regions; used to establish the range of values to be profiled.
`maxpts`	maximum number of points (in each direction, for each parameter) to evaluate in attempting to construct the profile.
`delta`	stepping scale for deciding on next point to profile.
`verbose`	level of output from internal calculations.
`devtol`	tolerance for fitted deviances less than baseline (supposedly minimum) deviance.
`maxmult`	maximum multiplier of the original step size allowed, defaults to 10.
`startmethod`	method for picking starting conditions for optimization (STUB).
`optimizer`	(character or function) optimizer to use (see `lmer` for details).
`signames`	logical indicating if abbreviated names of the form `.sigNN` should be used; otherwise, names are more meaningful (but longer) of the form `(sd\|cor)_(effects)\|(group)`. Note that some code for profile transformations (e.g., `varianceProf`) depends on `signames==TRUE`.
`...`	potential further arguments for various methods.
`x`	an object of class `thpr` (i.e., output of `profile`)
`base`	the base of the logarithm. Defaults to natural logarithms.

Details

Methods for function profile (package stats), here for profiling (fitted) mixed effect models.

Value

profile(<merMod>) returns an object of S3 class "thpr", data.frame-like. Methods for such a profile object are notably confint(), the three plots methods xyplot, densityplot, and splom.

Further, log() (see above) and as.data.frame().

Methods

signature(fitted = \"merMod\"): ...

Examples

fm01ML <- lmer(Yield ~ 1|Batch, Dyestuff, REML = FALSE)
system.time(
  tpr  <- profile(fm01ML, optimizer="Nelder_Mead", which="beta_")
)## fast; as only *one* beta parameter is profiled over
## full profiling (default which means 'all) needs
## ~2.6s (on a 2010 Macbook Pro)
system.time( tpr  <- profile(fm01ML))
## ~1s, + possible warning about bobyqa convergence
(confint(tpr) -> CIpr)

stopifnot(all.equal(CIpr,
  array(c(12.1985292, 38.2299848, 1486.4515,
          84.0630513, 67.6576964, 1568.54849), dim = 3:2,
        dimnames = list(c(".sig01", ".sigma", "(Intercept)"),
                        c("2.5 %", "97.5 %"))),
                    tol= 1e-07))# 1.37e-9 {64b}

require(lattice)
xyplot(tpr)
xyplot(tpr, absVal=TRUE) # easier to see conf.int.s (and check symmetry)
densityplot(tpr, main="densityplot( profile(lmer(..)) )")
splom(tpr)
doMore <- lme4:::testLevel() > 1 
if(doMore) { ## not typically, for time constraint reasons
 ## Batch and residual variance only
 system.time(tpr2 <- profile(fm01ML, which=1:2, optimizer="Nelder_Mead"))
 print( xyplot(tpr2) )
 print( xyplot(log(tpr2)) )# log(sigma) is better

 ## GLMM example
 gm1 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
	     data = cbpp, family = binomial)
 ## running time ~9 seconds on a modern machine:
 print( system.time(pr4 <- profile(gm1)) )
 print( xyplot(pr4,layout=c(5,1),as.table=TRUE) )
 print( xyplot(log(pr4), absVal=TRUE) ) # log(sigma_1)
 print( splom(pr4) )
 print( system.time( # quicker: only sig01 and one fixed effect
     pr2 <- profile(gm1, which=c("theta_", "period2"))))
 print( confint(pr2) )
}

[Package lme4 version 1.1-5 Index]