从混合模型获取随机效应矩阵 [英] Obtaining random-effects matrices from a mixed model

问题描述

在下面的代码中，我想知道如何从 lme()对象中的 out 和 Ts library(nlme)?

In my below code, I was wondering how I can obtain the equivalent of out and Ts from an lme() object in the library(nlme)?

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)

推荐答案

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

I would strongly suggest reading the docs! lme4 and nlme use inherently different methods for fitting mixed models -- lme4 uses a penalized least-squares formulation based on the lower Cholesky factor (theta) and nlme4 uses a generalized least-squares formulation that can optionally be stored as a Cholesky factor -- but their docs give you the information to get what you need from the internal representation. After that, it's up to you to do the math to convert between the representations.

如果您执行 lme ，则有一行

有关适合的组件，请参见 lmeObject

然后您执行?lmeObject ，在其中找到两个有希望的条目:

Then you do ?lmeObject, where you find two promising entries:

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

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 从类 lmeStruct 继承的对象，表示混合效果模型组件的列表，例如 reStruct ， corStruct 和 varFunc 对象.

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 数据的示例:

Well, we don't actually want the var-cov coefficients, but rather the random effects matrices. So we can look at reStruct. That is a lot more flexible in nlme than lme4, but it's generally just the random effects matrices. To do anything comparable to lme4, you need to transform those to their lower Cholesky factor. Here's an example using the sleepstudy data:

> 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.

The theta vector for lme4 is just a row-major representation of the lower triangle of the lower Cholesky factor for a given grouping variable. (For models with multiple grouping variables, you just concatenate these together.) The lower Cholesky factor is constrained to not have entries smaller than zero on the diagonal (because that would correspond to a negative variance), and otherwise not constrained. In other words, diagonal entries get a a lower bound at 0, all other entries get a lower bound at -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实现此功能(不是最有效的方法，但它显示了事物是如何构建的):

We can implement this for nlme (not the most efficient way, but it shows how things are built up):

> 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 限制集可以为不同条目创建不同的下限(例如，限制条目之间的关系)

Note that this is one spot where nlme is actually more powerful than lme4: the whole pdMatrix set of restrictions can create different lower bounds for different entries (as well as e.g. constrain the relationship between entries).

这篇关于从混合模型获取随机效应矩阵的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！