r语言 - 从混合模型中获得随机效应矩阵

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，则有一行

参见lmeObject的组件

然后你做?lmeObject，你会发现两个有希望的条目:

apVar
为方差-协方差系数的近似协方差矩阵。如果apVar = FALSE在调用lme时使用的控制值，则该组件为NULL。

和

modelStruct从类lmeStruct继承的对象，表示混合效果模型组件的列表，例如reStruct、corStruct和varFunc对象。

嗯，我们实际上并不需要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的θ向量只是给定分组变量的下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限制集可以为不同的条目创建不同的下界(以及例如约束条目之间的关系)。