profile-methods.RdMethods for profile() of [ng]lmer fitted
models.
The log() method and the more flexible logProf()
utility transform a lmer profile into one where logarithms of standard deviations
are used, while varianceProf converts from the
standard-deviation to the variance scale; see Details.
# S3 method for class 'merMod'
profile(fitted, which = NULL, alphamax = 0.01,
maxpts = 100, delta = NULL,
delta.cutoff = 1/8, verbose = 0, devtol = 1e-09,
devmatchtol = 1e-5,
maxmult = 10, startmethod = "prev", optimizer = NULL,
control=NULL, signames = TRUE,
parallel = c("no", "multicore", "snow"),
ncpus = getOption("profile.ncpus", 1L), cl = NULL,
prof.scale = c("sdcor","varcov"),
...)
# 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")
varianceProf(x, ranef = TRUE)a fitted model, e.g., the result of lmer(..).
NULL value, integer or character vector indicating which parameters to profile: default (NULL) is all parameters. For integer, i.e., indexing, the parameters are ordered as follows:
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.
residual standard deviation (or scale parameter for GLMMs where appropriate).
fixed effect (beta) parameters.
Alternatively, which may be a character, containing
"beta_" or "theta_" denoting the fixed or random
effects parameters, respectively, or also containing parameter
names, such as ".sigma" or "(Intercept)".
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.
maximum number of points (in each direction, for each parameter) to evaluate in attempting to construct the profile.
stepping scale for deciding on next point to profile.
The code uses the local derivative of the profile at the current
step to establish a change in the focal parameter that will lead
to a step of delta on the square-root-deviance scale.
If NULL, the delta.cutoff parameter will be used
to determine the stepping scale.
stepping scale (see delta)
expressed as a fraction of the
target maximum value of the profile on the square-root-deviance
scale. Thus a delta.cutoff setting of 1/n will
lead to a profile with approximately 2*n calculated points
for each parameter (i.e., n points in each direction,
below and above the estimate for each parameter).
level of output from internal calculations.
tolerance for fitted deviances less than baseline (supposedly minimum) deviance.
tolerance for match between original deviance computation and value returned from auxiliary deviance function
maximum multiplier of the original step size allowed, defaults to 10.
method for picking starting conditions for optimization (STUB).
(character or function) optimizer to use (see
lmer for details); default is to use the optimizer
from the original model fit.
a list of options controlling the
profiling (see lmerControl): default is to use the
control settings from the original model fit.
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., log() or
varianceProf) depends on signames==TRUE.
potential further arguments for various methods.
an object of class thpr (i.e., output of
profile)
the base of the logarithm. Defaults to natural logarithms.
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.
character string specifying the initial part of the sigma parameters to be log transformed.
The type of parallel operation to be used (if any).
If missing, the
default is taken from the option "profile.parallel" (and if that
is not set, "no").
integer: number of processes to be used in parallel operation: typically one would choose this to be the number of available CPUs.
An optional parallel or snow cluster for use if
parallel = "snow". If not supplied, a cluster on the
local machine is created for the duration of the profile
call.
whether to profile on the standard
deviation-correlation scale ("sdcor") or on
the variance-covariance scale ("varcov")
profile(<merMod>) returns an object of S3 class
"thpr",
which is data.frame-like.
Notable methods for such a profile object
confint(), which returns the
confidence intervals based on the profile,
and three plotting methods
(which require the lattice package),
xyplot, densityplot, and
splom.
In addition, the
log() (see above) and as.data.frame()
methods can transform "thpr" objects in useful ways.
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., the methods return a profile with all
of 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 random effect standard
deviations (.sigNN) can be reasonably be zero).
The forward and backward splines for the log-transformed parameters
are recalculated.
Note that correlation parameters are not handled sensibly at present
(i.e., they are logged rather than taking a more applicable
transformation such as an arc-hyperbolic tangent,
atanh(x)=\(\log((1+x)/(1-x))/2\)).
The varianceProf function works similarly, including
non-sensibility for correlation parameters, by squaring all
parameter values, changing the names by appending sq
appropriately (e.g. .sigNN to .sigsqNN).
Setting prof.scale="varcov" in the original
profile() call is a more computationally
intensive, but more correct, way to compute confidence
intervals for covariance parameters.
Methods for function profile (package
stats), here for profiling (fitted) mixed effect models.
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 -> 0.09s (2022)
#> user system elapsed
#> 0.071 0.011 0.077
## full profiling (default which means 'all') needs longer:
system.time( tpr <- profile(fm01ML, signames=FALSE))
#> user system elapsed
#> 0.385 0.284 0.261
## ~ 0.26s (2022) + possible warning about convergence
(confint(tpr) -> CIpr)
#> 2.5 % 97.5 %
#> sd_(Intercept)|Batch 12.19918 84.10119
#> sigma 38.25467 67.59099
#> (Intercept) 1486.44109 1568.53698
# too much precision (etc). but just FYI:
trgt <- array(c(12.19854, 38.22998, 1486.451,
84.06305, 67.6577, 1568.548), dim = 3:2)
stopifnot(all.equal(trgt, unname(CIpr), tol = .0001)) # had 3.1e-7
if (interactive()) {
library("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() > 2
if(doMore) { ## not typically, for time constraint reasons
## Batch and residual variance only
system.time(tpr2 <- profile(fm01ML, which=1:2)) # , optimizer="Nelder_Mead" gives warning
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 ~ 10 seconds on a modern machine {-> "verbose" while you wait}:
print( system.time(pr4 <- profile(gm1, verbose=TRUE)) )
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) )
## delta..: higher underlying resolution, only for 'sigma_1':
print( system.time(
pr4.hr <- profile(gm1, which="theta_", delta.cutoff=1/16)))
print( xyplot(pr4.hr) )
}
} # only if interactive()