profilemethods {lme4}  R Documentation 
Methods 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
standarddeviation 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 = 1e09,
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)
fitted 
a fitted model, e.g., the result of 
which 
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:
Alternatively, 
alphamax 
a number in 
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.
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.cutoff 
stepping scale (see 
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

control 
a 
signames 
logical indicating if abbreviated names of the form

... 
potential further arguments for various methods. 
x 
an object of class 
base 
the base of the logarithm. Defaults to natural logarithms. 
ranef 
logical indicating if the sigmas of the random effects
should be 
sigIni 
character string specifying the initial part of the sigma parameters to be log transformed. 
parallel 
The type of parallel operation to be used (if any).
If missing, the
default is taken from the option 
ncpus 
integer: number of processes to be used in parallel operation: typically one would choose this to be the number of available CPUs. 
cl 
An optional parallel or snow cluster for use if

prof.scale 
whether to profile on the standard
deviationcorrelation scale ( 
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 logtransformed 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 archyperbolic tangent,
atanh(x)
=\log((1+x)/(1x))/2
).
The varianceProf
function works similarly, including
nonsensibility 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.
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 plotting methods xyplot
etc, for class
"thpr"
.
For (more expensive) alternative confidence intervals:
bootMer
.
fm01ML < lmer(Yield ~ 1Batch, 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)
## full profiling (default which means 'all') needs longer:
system.time( tpr < profile(fm01ML, signames=FALSE))
## ~ 0.26s (2022) + possible warning about convergence
(confint(tpr) > CIpr)
# 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.1e7
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()