lme4 (version 1.1-7)

profile-methods: Profile method for merMod objects

Description

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

The log() method and the more flexible logProf() utility transform a lmer profile in one where logarithms of standard deviations are used, see Details.

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))
logProf(x, base = exp(1), ranef = TRUE,
           sigIni = if(ranef) "sig" else "sigma")

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: [object Object],[object Object],[object Object] In addition, which

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.

...

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.

ranef

logical indicating if the sigmas of the random effects should be log() transformed as well. If false, only $\sigma$ (standard deviation of errors) is transformed.

sigIni

character string specifying the initial part of the sigma parameters to be log transformed.

Value

• profile() returns an object of S3 class "thpr", data.frame-like. Methods for such a profile object are notably confint() and the three plotting methods (which require the lattice package), xyplot, densityplot, and splom.

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

Details

The log method and the more flexible logProf() function transform the profile into one where $\log(\sigma)$ is used instead of $\sigma$. By default all sigmas, including the standard deviations of the random effects are transformed, i.e., it returns a profile with all the .sigNN parameters replaced by .lsigNN. If ranef is false, only ".sigma", the standard deviation of the errors, is transformed as it should never be zero, whereas a random effect .sigNN is zero, when the corresponding random effect is entirely absent. The forward and backward splines for the log-transformed parameters are recalculated.

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

See Also

The plotting methods xyplot etc, for class "thpr"; varianceProf for transformation (from st.dev.) to variance scale.

For (more expensive) alternative confidence intervals: bootMer.

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)
xyplot(tpr, conf = c(0.95, 0.99), # (instead of all five 50, 80,...)
       main = "95% and 99% profile() intervals")
xyplot(logProf(tpr, ranef=FALSE),
       main = expression("lmer profile()s"~~ log(sigma)*"(only log)"))

densityplot(tpr, main="densityplot( profile(lmer(..)) )")
densityplot(varianceProf(tpr), main="varianceProf( profile(lmer(..)) )")
splom(tpr)
splom(logProf(tpr, ranef=FALSE))
doMore <- lme4:::testLevel() > 1if(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
 print( xyplot(logProf(tpr2, ranef=FALSE)) )

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