lme4:如何用随机截距指定 2 个相关性,而不添加随机斜率之间的相关性
lme4: how to specify 2 correlations with random intercept, without adding correlations between random slopes
我试图在 R 的 lme4 包中指定一个模型,其中我在随机截距和随机斜率之间有 2 个相关性,但不允许随机斜率相关。
lmer (Y ~ A + B + (1+A+B|Subject), data=mydata)
不好,因为它模拟了 A 和 B 的随机斜率之间的相关性。
鉴于:
lmer (Y ~ A + B + (1+A|Subject) + (1+B|Subject), data=mydata)
很糟糕,因为 Subject 的随机截距被引入模型两次。有没有第三种方法,也许更 hack-ish?
事实证明这比我想象的要难!
lme4
中的方差-协方差矩阵根据其Cholesky因子(本质上是矩阵平方根)进行参数化;因此,如果我们想建立一个特定相关性固定为零的模型,我们需要
t1 0 0 t1 t2 t3 t1^2 t1*t2 t1*t3
t2 t4 0 0 t4 t5 = t1*t2 t2^2 + t4^2 t2*t3 + t4*t5
t3 t5 t6 0 0 t6 = t1*t3 t2*t3 + t4*t5 t3^2 + t5^2 + t6^2
并求解 [3,2] 元素(A
和 B
之间的相关性)等于零;换句话说,我们将需要 t2 t3 + t4 t5 == 0
,或一个 6 元素向量,其中 t5 == -t2*t3/t4
;
tfun <- function(theta) {
theta5 <- -theta[2]*theta[3]/theta[4]
c(theta[1:4],theta5,theta[5])
}
模拟一些数据:
set.seed(101)
dd <- data.frame(A=rnorm(1000),B=rnorm(1000),
Subject=factor(rep(1:20,50)))
library("lme4")
dd$Y <- simulate(~A+B+(1+A+B|Subject),
newdata=dd,
family=gaussian(),
newparams=list(beta=c(1,2,3),
theta=tfun(c(1,0.2,0.3,2,3)),
sigma=1))[[1]]
现在按照 ?modular
中的步骤进行操作:
lmod <- lFormula(Y ~ A + B + (1+A+B|Subject), data=dd)
devfun <- do.call(mkLmerDevfun, lmod)
devfun()
的包装函数,它将采用 5 元素向量,计算相应的约束 theta 向量,并将其传递给 devfun()
:
devfun2 <- function(theta) {
devfun(tfun(theta))
}
从下界向量中删除一项:
lwr <- lmod$reTrms$lower
## [1] 0 -Inf -Inf 0 -Inf 0
lwr <- lwr[c(1:4,6)]
library("minqa")
## n.b. optwrap fails with minqa::bobyqa
opt <- lme4:::optwrap(optimizer=bobyqa,
par=ifelse(lwr==0,1,0),
fn=devfun2,
lower=lwr)
现在根据参数变换调整结果:
opt$par <- tfun(opt$par)
m1 <- mkMerMod(environment(devfun), opt, lmod$reTrms, fr = lmod$fr)
VarCorr(m1)
## Groups Name Std.Dev. Corr
## Subject (Intercept) 1.41450
## A 1.49374 0.019
## B 2.47895 0.316 0.000
## Residual 0.96617
所需的相关性现已固定为零。
我试图在 R 的 lme4 包中指定一个模型,其中我在随机截距和随机斜率之间有 2 个相关性,但不允许随机斜率相关。
lmer (Y ~ A + B + (1+A+B|Subject), data=mydata)
不好,因为它模拟了 A 和 B 的随机斜率之间的相关性。
鉴于:
lmer (Y ~ A + B + (1+A|Subject) + (1+B|Subject), data=mydata)
很糟糕,因为 Subject 的随机截距被引入模型两次。有没有第三种方法,也许更 hack-ish?
事实证明这比我想象的要难!
lme4
中的方差-协方差矩阵根据其Cholesky因子(本质上是矩阵平方根)进行参数化;因此,如果我们想建立一个特定相关性固定为零的模型,我们需要
t1 0 0 t1 t2 t3 t1^2 t1*t2 t1*t3
t2 t4 0 0 t4 t5 = t1*t2 t2^2 + t4^2 t2*t3 + t4*t5
t3 t5 t6 0 0 t6 = t1*t3 t2*t3 + t4*t5 t3^2 + t5^2 + t6^2
并求解 [3,2] 元素(A
和 B
之间的相关性)等于零;换句话说,我们将需要 t2 t3 + t4 t5 == 0
,或一个 6 元素向量,其中 t5 == -t2*t3/t4
;
tfun <- function(theta) {
theta5 <- -theta[2]*theta[3]/theta[4]
c(theta[1:4],theta5,theta[5])
}
模拟一些数据:
set.seed(101)
dd <- data.frame(A=rnorm(1000),B=rnorm(1000),
Subject=factor(rep(1:20,50)))
library("lme4")
dd$Y <- simulate(~A+B+(1+A+B|Subject),
newdata=dd,
family=gaussian(),
newparams=list(beta=c(1,2,3),
theta=tfun(c(1,0.2,0.3,2,3)),
sigma=1))[[1]]
现在按照 ?modular
中的步骤进行操作:
lmod <- lFormula(Y ~ A + B + (1+A+B|Subject), data=dd)
devfun <- do.call(mkLmerDevfun, lmod)
devfun()
的包装函数,它将采用 5 元素向量,计算相应的约束 theta 向量,并将其传递给 devfun()
:
devfun2 <- function(theta) {
devfun(tfun(theta))
}
从下界向量中删除一项:
lwr <- lmod$reTrms$lower
## [1] 0 -Inf -Inf 0 -Inf 0
lwr <- lwr[c(1:4,6)]
library("minqa")
## n.b. optwrap fails with minqa::bobyqa
opt <- lme4:::optwrap(optimizer=bobyqa,
par=ifelse(lwr==0,1,0),
fn=devfun2,
lower=lwr)
现在根据参数变换调整结果:
opt$par <- tfun(opt$par)
m1 <- mkMerMod(environment(devfun), opt, lmod$reTrms, fr = lmod$fr)
VarCorr(m1)
## Groups Name Std.Dev. Corr
## Subject (Intercept) 1.41450
## A 1.49374 0.019
## B 2.47895 0.316 0.000
## Residual 0.96617
所需的相关性现已固定为零。