Extract or Get Generalized Components from a Fitted Mixed Effects Model
getME.RdExtract (or “get”) “components” – in a
generalized sense – from a fitted mixed-effects model,
i.e., (in this version of the package) from an object of
class "merMod".
Usage
getME(object, name, ...)
# S3 method for class 'merMod'
getME(object,
name = c("X", "Z", "Zt", "Ztlist", "mmList", "y", "mu", "u", "b",
"Gp", "Tp", "L", "Lambda", "Lambdat", "Lind", "Tlist",
"A", "RX", "RZX", "sigma", "flist",
"fixef", "beta", "theta", "ST", "REML", "is_REML",
"n_rtrms", "n_rfacs", "N", "n", "p", "q",
"p_i", "l_i", "q_i", "k", "m_i", "m",
"cnms", "devcomp", "offset", "lower", "devfun", "glmer.nb.theta"),
...)Arguments
- object
a fitted mixed-effects model of class
"merMod", i.e., typically the result oflmer(),glmer()ornlmer().- name
a character vector specifying the name(s) of the “component”. If
length(name) > 1or ifname = "ALL", a namedlistof components will be returned. Possible values are:"X":fixed-effects model matrix
"Z":random-effects model matrix
"Zt":transpose of random-effects model matrix. Note that the structure of
Zthas changed sincelme4.0; to get a backward-compatible structure, usedo.call(Matrix::rBind,getME(.,"Ztlist"))"Ztlist":list of components of the transpose of the random-effects model matrix, separated by individual variance component
"mmList":list of raw model matrices associated with random effects terms
"y":response vector
"mu":conditional mean of the response
"u":conditional mode of the “spherical” random effects variable
"b":conditional mode of the random effects variable
"Gp":groups pointer vector. A pointer to the beginning of each group of random effects corresponding to the random-effects terms, beginning with 0 and including a final element giving the total number of random effects
"Tp":theta pointer vector. A pointer to the beginning of the theta sub-vectors corresponding to the random-effects terms, beginning with 0 and including a final element giving the number of thetas.
"L":sparse Cholesky factor of the penalized random-effects model.
"Lambda":relative covariance factor \(\Lambda\) of the random effects.
"Lambdat":transpose \(\Lambda'\) of \(\Lambda\) above.
"Lind":index vector for inserting elements of \(\theta\) into the nonzeros of \(\Lambda\).
"Tlist":vector of template matrices from which the blocks of \(\Lambda\) are generated.
"A":Scaled sparse model matrix (class
"dgCMatrix") for the unit, orthogonal random effects, \(U\), equal togetME(.,"Zt") %*% getME(.,"Lambdat")"RX":Cholesky factor for the fixed-effects parameters
"RZX":cross-term in the full Cholesky factor
"sigma":residual standard error; note that
sigma(object)is preferred."flist":a list of the grouping variables (factors) involved in the random effect terms
"fixef":fixed-effects parameter estimates
"beta":fixed-effects parameter estimates (identical to the result of
fixef, but without names)"theta":random-effects parameter estimates: these are parameterized as the relative Cholesky factors of each random effect term
"ST":A list of S and T factors in the TSST' Cholesky factorization of the relative variance matrices of the random effects associated with each random-effects term. The unit lower triangular matrix, \(T\), and the diagonal matrix, \(S\), for each term are stored as a single matrix with diagonal elements from \(S\) and off-diagonal elements from \(T\).
"n_rtrms":number of random-effects terms
"n_rfacs":number of distinct random-effects grouping factors
"N":number of rows of
X"n":length of the response vector,
y"p":number of columns of the fixed effects model matrix,
X"q":number of columns of the random effects model matrix,
Z"p_i":numbers of columns of the raw model matrices,
mmList"l_i":numbers of levels of the grouping factors
"q_i":numbers of columns of the term-wise model matrices,
ZtList"k":number of random effects terms
"m_i":numbers of covariance parameters in each term
"m":total number of covariance parameters, i.e., the same as
dims@nthbelow."cnms":the “component names”, a
list."REML":0indicates the model was fitted by maximum likelihood, any other positive integer indicates fitting by restricted maximum likelihood"is_REML":same as the result of
isREML(.)"devcomp":a list consisting of a named numeric vector,
cmp, and a named integer vector,dims, describing the fitted model. The elements ofcmpare:- ldL2
twice the log determinant of
L- ldRX2
twice the log determinant of
RX- wrss
weighted residual sum of squares
- ussq
squared length of
u- pwrss
penalized weighted residual sum of squares, “wrss + ussq”
- drsum
sum of residual deviance (GLMMs only)
- REML
REML criterion at optimum (LMMs fit by REML only)
- dev
deviance criterion at optimum (models fit by ML only)
- sigmaML
ML estimate of residual standard deviation
- sigmaREML
REML estimate of residual standard deviation
- tolPwrss
tolerance for declaring convergence in the penalized iteratively weighted residual sum-of-squares (GLMMs only)
The elements of
dimsare:- N
number of rows of
X- n
length of
y- p
number of columns of
X- nmp
n-p- nth
length of
theta- q
number of columns of
Z- nAGQ
see
glmer- compDev
see
glmerControl- useSc
TRUEif model has a scale parameter- reTrms
number of random effects terms
- REML
0indicates the model was fitted by maximum likelihood, any other positive integer indicates fitting by restricted maximum likelihood- GLMM
TRUEif a GLMM- NLMM
TRUEif an NLMM
"offset":model offset
"lower":lower bounds on random-effects model parameters (i.e, "theta" parameters). In order to constrain random effects covariance matrices to be semi-positive-definite, this vector is equal to 0 for elements of the
thetavector corresponding to diagonal elements of the Cholesky factor,-Infotherwise. (getME(.,"lower")==0can be used as a test to identify diagonal elements, as inisSingular.)"devfun":deviance function (so far only available for LMMs)
"glmer.nb.theta":negative binomial \(\theta\) parameter, only for
glmer.nb.
%% -- keep at the very end:
"ALL":get all of the above as a
list.
- ...
currently unused in lme4, potentially further arguments in methods.
Value
Unspecified, as very much depending on the name.
Details
The goal is to provide “everything a user may
want” from a fitted "merMod" object as far
as it is not available by methods, such as
fixef, ranef,
vcov, etc.
Examples
## shows many methods you should consider *before* using getME():
methods(class = "merMod")
#> [1] PBmodcomp PBrefdist VarCorr
#> [4] anova as.function coef
#> [7] confint cooks.distance deviance
#> [10] df.residual drop1 extractAIC
#> [13] family fitted fixef
#> [16] formula getData getL
#> [19] getME hatvalues influence
#> [22] isGLMM isLMM isNLMM
#> [25] isREML logLik model.frame
#> [28] model.matrix model2restriction_matrix ngrps
#> [31] nobs plot predict
#> [34] print profile qqmath
#> [37] ranef rePCA refit
#> [40] refitML residuals rstudent
#> [43] show sigma simulate
#> [46] summary terms update
#> [49] vcov weights
#> see '?methods' for accessing help and source code
(fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy))
#> Linear mixed model fit by REML ['lmerMod']
#> Formula: Reaction ~ Days + (Days | Subject)
#> Data: sleepstudy
#> REML criterion at convergence: 1743.628
#> Random effects:
#> Groups Name Std.Dev. Corr
#> Subject (Intercept) 24.741
#> Days 5.922 0.07
#> Residual 25.592
#> Number of obs: 180, groups: Subject, 18
#> Fixed Effects:
#> (Intercept) Days
#> 251.41 10.47
Z <- getME(fm1, "Z")
stopifnot(is(Z, "CsparseMatrix"),
c(180,36) == dim(Z),
all.equal(fixef(fm1), b1 <- getME(fm1, "beta"),
check.attributes=FALSE, tolerance = 0))
## A way to get *all* getME()s :
## internal consistency check ensuring that all work:
parts <- getME(fm1, "ALL")
str(parts, max=2)
#> List of 45
#> $ X : num [1:180, 1:2] 1 1 1 1 1 1 1 1 1 1 ...
#> ..- attr(*, "dimnames")=List of 2
#> ..- attr(*, "assign")= int [1:2] 0 1
#> ..- attr(*, "msgScaleX")= chr(0)
#> $ Z :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
#> $ Zt :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
#> $ Ztlist :List of 2
#> ..$ Subject.(Intercept):Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
#> ..$ Subject.Days :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
#> $ mmList :List of 1
#> ..$ Days | Subject: num [1:180, 1:2] 1 1 1 1 1 1 1 1 1 1 ...
#> .. ..- attr(*, "dimnames")=List of 2
#> .. ..- attr(*, "assign")= int [1:2] 0 1
#> $ y : num [1:180] 250 259 251 321 357 ...
#> $ mu : num [1:180] 254 273 293 313 332 ...
#> $ u : num [1:36] 2.34 39.68 -41.79 -34.58 -40.3 ...
#> $ b :Formal class 'dgeMatrix' [package "Matrix"] with 4 slots
#> $ Gp : int [1:2] 0 36
#> $ Tp : Named num [1:2] 0 3
#> ..- attr(*, "names")= chr [1:2] "beg__" "Subject"
#> $ L :Formal class 'dCHMsimpl' [package "Matrix"] with 11 slots
#> $ Lambda :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
#> $ Lambdat :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
#> $ Lind : int [1:54] 1 2 3 1 2 3 1 2 3 1 ...
#> $ Tlist :List of 1
#> ..$ Subject: num [1:2, 1:2] 0.9667 0.0152 0 0.2309
#> $ A :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
#> $ RX : num [1:2, 1:2] 3.79 0 2.3 16.56
#> ..- attr(*, "dimnames")=List of 2
#> $ RZX : num [1:36, 1:2] 3.022 0.269 3.022 0.269 3.022 ...
#> ..- attr(*, "dimnames")=List of 2
#> $ sigma : num 25.6
#> $ flist :List of 1
#> ..$ Subject: Factor w/ 18 levels "308","309","310",..: 1 1 1 1 1 1 1 1 1 1 ...
#> ..- attr(*, "assign")= int 1
#> $ fixef : Named num [1:2] 251.4 10.5
#> ..- attr(*, "names")= chr [1:2] "(Intercept)" "Days"
#> $ beta : num [1:2] 251.4 10.5
#> $ theta : Named num [1:3] 0.9667 0.0152 0.2309
#> ..- attr(*, "names")= chr [1:3] "Subject.(Intercept)" "Subject.Days.(Intercept)" "Subject.Days"
#> $ ST :List of 1
#> ..$ Subject: num [1:2, 1:2] 0.9667 0.0157 0 0.2309
#> $ REML : int 2
#> $ is_REML : logi TRUE
#> $ n_rtrms : int 1
#> $ n_rfacs : int 1
#> $ N : int 180
#> $ n : int 180
#> $ p : int 2
#> $ q : int 36
#> $ p_i : Named int 2
#> ..- attr(*, "names")= chr "Days | Subject"
#> $ l_i : Named int 18
#> ..- attr(*, "names")= chr "Subject"
#> $ q_i : Named int 36
#> ..- attr(*, "names")= chr "Days | Subject"
#> $ k : int 1
#> $ m_i : Named num 3
#> ..- attr(*, "names")= chr "Days | Subject"
#> $ m : int 3
#> $ cnms :List of 1
#> ..$ Subject: chr [1:2] "(Intercept)" "Days"
#> $ devcomp :List of 2
#> ..$ cmp : Named num [1:10] 7.60e+01 8.28 9.89e+04 1.77e+04 1.17e+05 ...
#> .. ..- attr(*, "names")= chr [1:10] "ldL2" "ldRX2" "wrss" "ussq" ...
#> ..$ dims: Named int [1:12] 180 180 2 178 36 3 1 1 0 2 ...
#> .. ..- attr(*, "names")= chr [1:12] "N" "n" "p" "nmp" ...
#> $ offset : num [1:180] 0 0 0 0 0 0 0 0 0 0 ...
#> $ lower : Named num [1:3] 0 -Inf 0
#> ..- attr(*, "names")= chr [1:3] "Subject.(Intercept)" "Subject.Days.(Intercept)" "Subject.Days"
#> $ devfun :function (theta)
#> $ glmer.nb.theta: logi NA
stopifnot(identical(Z, parts $ Z),
identical(b1, parts $ beta))