如何在 lmer 中修改具有成对随机效应的插槽
How to modify slots with paired random effects in lmer
这是之前 post (How to modify slots lme4 >1.0) 的后续问题。我有一个类似的成对数据结构,并希望随机效应同时考虑成对的 "pops"。我有一个使用之前建议的代码的功能性随机拦截模型:
dat <- data.frame(pop1 = c(2,1,1,1,1,3,2,2,2,3,5,3,5,4,6),
pop2 = c(1,3,4,5,6,2,4,5,6,4,3,6,4,6,5),
X = c(20,25,18,40,36,70,68,72,78,76,97,100,115,110,108),
Y = c(18,16,15,40,22,18,18,18,18,45,10,47,67,5,6))
#build random effects matrix
Zl<-lapply(c("pop1","pop2"),function(nm)Matrix:::fac2sparse(dat[[nm]],"d",drop=FALSE))
ZZ<-Reduce("+",Zl[-1],Zl[[1]])
#specify model structure
mod<-lFormula(Y~X+(1|pop1),data=dat,REML=TRUE)
#replace slot
mod$reTrms$Zt <- ZZ
#fit model
dfun<-do.call(mkLmerDevfun,mod)
opt<-optimizeLmer(dfun)
mkMerMod(environment(dfun),opt,mod$reTrms,fr=mod$fr)
但是,当尝试添加随机斜率变量时:
mod2<-lFormula(Y~X+(1+X|pop1),data=dat,REML=TRUE)
mod2$reTrms$Zt <- ZZ
dfun<-do.call(mkLmerDevfun,mod2)
导致与之前 post 相同的错误(问题是调用了错误的数据帧):“Lambdat %*% Ut 错误:
文件 ../MatrixOps/cholmod_ssmult.c,第 82 行“Cholmod 错误 'A and B inner dimensions must match'”
查看每个 pop 的 lm
plot(1,type="n",xlim=c(0,150),ylim=c(0,75),ylab = "Y",xlab="X")
for(i in 1:length(unique(c(dat$pop1,dat$pop2)))){
subdat<-dat[which(dat$pop1==i | dat$pop2==i),]
out<-summary(lm(subdat$Y~subdat$X))
x=subdat$X
y=x*out$coefficients[2,1]+out$coefficients[1,1]
lines(x,y,col=i))
}
legend(125,60,1:6,col=1:6,lty=1,title="Pop")
dat <- data.frame(pop1 = c(2,1,1,1,1,3,2,2,2,3,5,3,5,4,6),
pop2 = c(1,3,4,5,6,2,4,5,6,4,3,6,4,6,5),
X = c(20,25,18,40,36,70,68,72,78,76,97,100,115,110,108),
Y = c(18,16,15,32,22,29,32,38,44,45,51,47,67,59,61))
有助于理解原始代码的实际作用:
## build random effects matrix
## 1. sparse dummy-variable matrices for each population ID
Zl <- lapply(dat[c("pop1","pop2")],
Matrix::fac2sparse,to="d",drop.unused.levels=FALSE)
## 2. take the sum of all components of the list of dummy-variable matrices ...
ZZ <- Reduce("+",Zl[-1],Zl[[1]])
如果列表很长,Reduce
形式通常很方便,但在这种情况下它只是 Zl[[1]]+Zl[[2]]
...
all.equal(Zl[[1]]+Zl[[2]],ZZ) ## TRUE
这个 RE 结构是什么样的?
library(gridExtra)
grid.arrange(
image(t(Zl[[1]]),main="pop 1",sub="",xlab="pop",ylab="obs"),
image(t(Zl[[2]]),main="pop 2",sub="",xlab="pop",ylab="obs"),
image(t(ZZ),main="combined",sub="",xlab="RE",ylab="obs"),
nrow=1)
对于随机斜率,我 认为 我们想要获取 ZZ
的每个填充元素并将其替换为相应 observation/row 观察到的 X
值dat
:这里的索引有点模糊 - 在这种情况下,它归结为 Z
的每一行/Zt
的列中有 2 个填充值(@p
稀疏矩阵的槽给出了一个指向每列中第一个非零元素的零索引指针...)
vals <- dat$X[rep(1:(length(ZZ@p)-1),diff(ZZ@p))]
ZZX <- ZZ
ZZX@x <- vals
image(t(ZZX))
library(lme4)
mod <- lFormula(Y~X+(X|pop1),data=dat,REML=TRUE)
## replace slot
mod$reTrms$Zt <- rbind(ZZ,ZZX)
## fit model
dfun <- do.call(mkLmerDevfun,mod)
opt <- optimizeLmer(dfun)
m1 <- mkMerMod(environment(dfun),opt,mod$reTrms,fr=mod$fr)
这 似乎 可以工作,但是您当然应该根据自己的知识来检查它,了解这里应该发生什么...
这是之前 post (How to modify slots lme4 >1.0) 的后续问题。我有一个类似的成对数据结构,并希望随机效应同时考虑成对的 "pops"。我有一个使用之前建议的代码的功能性随机拦截模型:
dat <- data.frame(pop1 = c(2,1,1,1,1,3,2,2,2,3,5,3,5,4,6),
pop2 = c(1,3,4,5,6,2,4,5,6,4,3,6,4,6,5),
X = c(20,25,18,40,36,70,68,72,78,76,97,100,115,110,108),
Y = c(18,16,15,40,22,18,18,18,18,45,10,47,67,5,6))
#build random effects matrix
Zl<-lapply(c("pop1","pop2"),function(nm)Matrix:::fac2sparse(dat[[nm]],"d",drop=FALSE))
ZZ<-Reduce("+",Zl[-1],Zl[[1]])
#specify model structure
mod<-lFormula(Y~X+(1|pop1),data=dat,REML=TRUE)
#replace slot
mod$reTrms$Zt <- ZZ
#fit model
dfun<-do.call(mkLmerDevfun,mod)
opt<-optimizeLmer(dfun)
mkMerMod(environment(dfun),opt,mod$reTrms,fr=mod$fr)
但是,当尝试添加随机斜率变量时:
mod2<-lFormula(Y~X+(1+X|pop1),data=dat,REML=TRUE)
mod2$reTrms$Zt <- ZZ
dfun<-do.call(mkLmerDevfun,mod2)
导致与之前 post 相同的错误(问题是调用了错误的数据帧):“Lambdat %*% Ut 错误: 文件 ../MatrixOps/cholmod_ssmult.c,第 82 行“Cholmod 错误 'A and B inner dimensions must match'”
查看每个 pop 的 lm
plot(1,type="n",xlim=c(0,150),ylim=c(0,75),ylab = "Y",xlab="X")
for(i in 1:length(unique(c(dat$pop1,dat$pop2)))){
subdat<-dat[which(dat$pop1==i | dat$pop2==i),]
out<-summary(lm(subdat$Y~subdat$X))
x=subdat$X
y=x*out$coefficients[2,1]+out$coefficients[1,1]
lines(x,y,col=i))
}
legend(125,60,1:6,col=1:6,lty=1,title="Pop")
dat <- data.frame(pop1 = c(2,1,1,1,1,3,2,2,2,3,5,3,5,4,6),
pop2 = c(1,3,4,5,6,2,4,5,6,4,3,6,4,6,5),
X = c(20,25,18,40,36,70,68,72,78,76,97,100,115,110,108),
Y = c(18,16,15,32,22,29,32,38,44,45,51,47,67,59,61))
有助于理解原始代码的实际作用:
## build random effects matrix
## 1. sparse dummy-variable matrices for each population ID
Zl <- lapply(dat[c("pop1","pop2")],
Matrix::fac2sparse,to="d",drop.unused.levels=FALSE)
## 2. take the sum of all components of the list of dummy-variable matrices ...
ZZ <- Reduce("+",Zl[-1],Zl[[1]])
如果列表很长,Reduce
形式通常很方便,但在这种情况下它只是 Zl[[1]]+Zl[[2]]
...
all.equal(Zl[[1]]+Zl[[2]],ZZ) ## TRUE
这个 RE 结构是什么样的?
library(gridExtra)
grid.arrange(
image(t(Zl[[1]]),main="pop 1",sub="",xlab="pop",ylab="obs"),
image(t(Zl[[2]]),main="pop 2",sub="",xlab="pop",ylab="obs"),
image(t(ZZ),main="combined",sub="",xlab="RE",ylab="obs"),
nrow=1)
ZZ
的每个填充元素并将其替换为相应 observation/row 观察到的 X
值dat
:这里的索引有点模糊 - 在这种情况下,它归结为 Z
的每一行/Zt
的列中有 2 个填充值(@p
稀疏矩阵的槽给出了一个指向每列中第一个非零元素的零索引指针...)
vals <- dat$X[rep(1:(length(ZZ@p)-1),diff(ZZ@p))]
ZZX <- ZZ
ZZX@x <- vals
image(t(ZZX))
library(lme4)
mod <- lFormula(Y~X+(X|pop1),data=dat,REML=TRUE)
## replace slot
mod$reTrms$Zt <- rbind(ZZ,ZZX)
## fit model
dfun <- do.call(mkLmerDevfun,mod)
opt <- optimizeLmer(dfun)
m1 <- mkMerMod(environment(dfun),opt,mod$reTrms,fr=mod$fr)
这 似乎 可以工作,但是您当然应该根据自己的知识来检查它,了解这里应该发生什么...