修正lme4/lmer中特定值的差异
我正在做一个混合效应模型的模拟研究(三个层次；观察嵌套在学校的科目中)：
f <- lmer(measurement ~ time + race + gender + s_ses +
fidelity + (1 + time|school/subject), mydata_long, REML=0)该模型允许截距和时间斜率随学科和学校的不同而变化。我想知道如何才能将方差修正为特定值。我知道如何在只有随机截获的情况下做到这一点：
VarCorr(f)['subject:school']<-0.13
VarCorr(f)['school']<-0.20然而，当存在随机斜率时，这些代码不起作用，因为在方差方面有不同的分量(见附图)。
在这种情况下，如何将subject:学校(Intercept)、subject:学校时间、学校(Intercept)和学校时间的差异固定为特定值。有什么建议吗？
回答 1
Stack Overflow用户
发布于 2020-07-30 02:45:11
一个模拟的例子。最困难的部分是正确指定随机效果参数:您需要知道的关键事情是(1)内部随机效果方差矩阵由残差方差缩放；(2)对于矢量值随机效果(如此随机斜率模型)，方差协方差矩阵是根据其乔列斯基因子指定的:如果我们想要协方差矩阵V，则存在一个下三角矩阵，使得C %*% t(C) == V。我们使用chol()计算C，然后以列为主的顺序读出下三角形(包括对角线)的元素(参见下面的辅助函数)。
设置实验设计(简化自您的设计，但使用相同的随机效果组件)：
mydata_long <- expand.grid(time=1:40,
school=factor(letters[1:25]),
subject=factor(LETTERS[1:25]))从转换的帮助器函数
- 标准差向量、一个或多个相关参数(按下三角/列主顺序)和残差标准差
到
- 一个由
lme4内部使用的"theta“参数的向量(参见上面的描述)
..。从另一个方向返回(conv_chol)
conv_sc <- function(sdvec,cor,sigma) {
## construct symmetric matrix with cor in lower/upper triangles
cormat <- matrix(1,nrow=length(sdvec),ncol=length(sdvec))
cormat[lower.tri(cormat)] <- cor
cormat[upper.tri(cormat)] <- t(cormat)[upper.tri(cormat)]
## convert to covariance matrix and scale by 1/sigma^2
V <- outer(sdvec, sdvec)*cormat/sigma^2
## extract lower triangle in column-major order
return(t(chol(V))[lower.tri(V,diag=TRUE)])
}
conv_chol <- function(ch, s) {
m <- matrix(NA,2,2)
m[lower.tri(m,diag=TRUE)] <- ch
m[upper.tri(m)] <- 0
V <- m %*% t(m) * s^2
list(sd=sqrt(diag(V)), cor=cov2cor(V)[1,2])
}如果您想从协方差矩阵开始，而不是从标准差和相关性开始，您可以修改代码以跳过一些步骤(以V开始和结束)。
选取一些值并进行转换(和反向转换，以进行检查)
tt1 <- conv_sc(c(0.7, 1.2), 0.3, 0.5)
tt2 <- conv_sc(c(1.4, 0.2), -0.2, 0.5)
tt <- c(tt1, tt2)
conv_chol(tt1, s=0.5)
conv_chol(tt2, s=0.5)设置公式并进行模拟：
form <- m ~ time + (1 + time|school/subject)
set.seed(101)
mydata_long$m <- simulate(form[-2], ## [-2] drops the response
family=gaussian,
newdata=mydata_long,
newparams=list(theta=tt,
beta=c(1,1),
sigma=0.5))[[1]]
f <- lmer(form, data=mydata_long, REML=FALSE)
VarCorr(f)拟合的结果接近我们上面要求的结果。
Groups Name Std.Dev. Corr
subject:school (Intercept) 0.66427
time 1.16488 0.231
school (Intercept) 1.78312
time 0.22459 -0.156
Residual 0.49772 现在做同样的事情200次，来探索估计的分布：
simfun <- function() {
mydata_long$m <- simulate(form[-2],
family=gaussian,
newdata=mydata_long,
newparams=list(theta=tt,
beta=c(1,1),
sigma=0.5))[[1]]
f <- lmer(form, data=mydata_long, REML=FALSE)
return(as.data.frame(VarCorr(f))[,"sdcor"])
}
set.seed(101)
res <- plyr::raply(200,suppressMessages(simfun()),.progress="text")为了方便起见，这里使用了循环，您可以随意执行此操作(for plyr::raply()、lapply()、replicate()、purrr::map() ...)
par(las=1)
boxplot(res)
## add true values to the plot
points(1:7,c(0.7,1.2,0.3,1.4,0.2,-0.3,0.5),col=2,cex=3,lwd=3)
https://stackoverflow.com/questions/63145801
复制相似问题
腾讯云开发者
