从混合模型中获取随机效应矩阵

Question

在下面的代码中，我想知道如何从 library(nlme) 中的 lme() 对象中获取 out 和 Ts 的等价物？

dat <- read.csv("https://raw.githubusercontent.com/rnorouzian/v/main/mv.l.csv")

library(lme4)

x <- lmer(value ~0 + name+ (1| School/Student), data = dat,
          control = lmerControl(check.nobs.vs.nRE= "ignore"))

    lwr <- getME(x, "lower")
    theta <- getME(x, "theta")

    out = any(theta[lwr == 0] < 1e-4)  # find this from `x1` object below
     Ts = getME(x, "Tlist")            # find this from `x1` object below


# Fitting the above model using `library(nlme)`:

library(nlme)

x1 <- lme(value ~0 + name, random = ~1| School/Student, data = dat)

Answer 1

我强烈建议阅读文档！ lme4 和 nlme 使用本质上不同的方法来拟合混合模型——lme4 使用基于较低的 Cholesky 因子 (theta) 的惩罚最小二乘公式，nlme4 使用广义最小二乘公式，可以选择将其存储为 Cholesky 因子——但是他们的文档为您提供了从内部表示中获取所需信息的信息。之后，由您进行数学运算以在表示之间进行转换。

如果你这样做 ?lme 那么就会有行

See lmeObject for the components of the fit

然后你做 ?lmeObject，在那里你找到两个有前途的条目：

apVar
an approximate covariance matrix for the variance-covariance coefficients. If apVar = FALSE in the control values used in the call to lme, this component is NULL.

和

modelStruct an object inheriting from class lmeStruct, representing a list of mixed-effects model components, such as reStruct, corStruct, and varFunc objects.

好吧，我们实际上不需要 var-cov 系数，而是随机效应矩阵。所以我们可以看看reStruct。这在 nlme 中比 lme4 灵活得多，但它通常只是随机效应矩阵。要执行与 lme4 相当的任何操作，您需要将它们转换为较低的 Cholesky 因子。这是使用 sleepstudy 数据的示例：

> library("nlme")
> library("lme4")
> 
> data("sleepstudy")
> m_nlme <- lme(fixed=Reaction ~ 1 + Days,
+               random=~ 1 + Days | Subject,
+               data=sleepstudy,
+               method = "ML")
> m_lme4 <- lmer(Reaction ~ 1 + Days + (1 + Days|Subject),
+                data=sleepstudy,
+                REML=FALSE)
> 
> re_lme4 <- getME(m_lme4, "Tlist")$Subject
> print(re_lme4)
           [,1]      [,2]
[1,] 0.92919061 0.0000000
[2,] 0.01816575 0.2226432
> 
> re_nlme <- m_nlme$modelStruct$reStruct$Subject
> # entire covariance matrix
> print(re_nlme)
Positive definite matrix structure of class pdLogChol representing
            (Intercept)       Days
(Intercept)  0.86344433 0.01688228
Days         0.01688228 0.04990040
> # get the lower cholesky factor
> re_nlme <- t(chol(re_nlme)) # could also use pdMatrix(re_nlme, TRUE)
> print(re_nlme)
            (Intercept)      Days
(Intercept)  0.92921705 0.0000000
Days         0.01816829 0.2226439

lme4 的 theta 向量只是给定分组变量的下 Cholesky 因子的下三角形的行主表示。（对于具有多个分组变量的模型，您只需将它们连接在一起。）较低的 Cholesky 因子被限制为在对角线上没有小于零的条目（因为这将对应于负方差），否则不受限制。换句话说，对角线条目的下限为 0，所有其他条目的下限为 -Inf.

因此，在 lme4 中：

> re_lme4[lower.tri(re_lme4,diag = TRUE)]
[1] 0.92919061 0.01816575 0.22264321
> getME(m_lme4, "theta")
     Subject.(Intercept) Subject.Days.(Intercept)             Subject.Days 
              0.92919061               0.01816575               0.22264321 
> getME(m_lme4, "lower")
[1]    0 -Inf    0

我们可以为 nlme 实现这个（不是最有效的方法，但它显示了事情是如何构建的）：

> lowerbd <- function(x){
+   dd <- diag(0, nrow=nrow(x))
+   dd[lower.tri(dd)] <- -Inf
+   dd[lower.tri(dd, diag=TRUE)]
+ }
> lowerbd(re_nlme)
[1]    0 -Inf    0
> lowerbd(re_lme4)
[1]    0 -Inf    0

请注意，这是 nlme 实际上比 lme4 更强大的地方：整个 pdMatrix 限制集可以为不同的条目创建不同的下限（以及例如限制条目之间的关系）。

从混合模型中获取随机效应矩阵

Obtaining random-effects matrices from a mixed model

r

regression

lme4

mixed-models

nlme