在 lme4 中调整 object@pp$X答案

【问题标题】：Adjusting object@pp$X in lme4在 lme4 中调整 object@pp$X
【发布时间】：2017-05-05 14:08:47
【问题描述】：

我正在测试 lmer 模型中的一些固定系数，但需要在进一步的过程中使用该模型（计算每个变量的贡献），因此需要更改 lmerMod 模型的某些部分。

由于以下错误消息，我正在努力更改 object@pp$X

错误：无效替换：引用类字段“X”是只读的

下面的可重现示例：

#Load package and data
library(lme4)
data(iris)

#build the model
mod<-lmer(Sepal.Length~Petal.Length + offset(Petal.Width*1) + (1|Species),data=iris)

fixef(mod) #not showing the offset coefficient

#apply changes to mod to get fixef(mod) to work with new coefficient
mod@beta <- c(mod@beta,1) #because model was offset by 1*Petal.Width
mod@pp$X <- matrix(data.frame(mod@pp$X, iris["Petal.Width"])) #causes the error

#check fixef:
fixef(mod) # should have Petal.Width at the end with a value of 1

注意：

[fixef]:(https://github.com/lme4/lme4/blob/master/R/lmer.R#L876)！在它的函数中有两个调用：

是object@beta（已经改成功了）；
是getME(object,"X"):(https://github.com/lme4/lme4/blob/master/R/lmer.R#L1932)。

我对使用变量名称获取 fixef 系数的替代方法持开放态度（能够直接调整 lmerMod）

提前致谢！

【问题讨论】：

标签： r lme4

【解决方案1】：

对象mod中槽pp的类为merPredD：

library(lme4)
data(iris)
mod <- lmer(Sepal.Length~Petal.Length + offset(Petal.Width*1) + (1|Species),data=iris)
class(slot(mod,"pp"))

[1] "merPredD"
attr(,"package")
[1] "lme4"

此类lme4对象可以通过lme4的merPredD命令生成。
下面是使用mod@pp的相同参数生成merPredD对象的示例：

obj1 <- merPredD(X=mod@pp$X, Zt=mod@pp$Zt, Lambdat=mod@pp$Lambdat, Lind=mod@pp$Lind, 
         theta=mod@pp$theta, n=nrow(mod@pp$X))
class(obj1)

使用mod@pp <- obj1 我们不会收到错误消息，我们可以用obj1 替换mod@pp 对象。
按照这种方式，我们可以例如更改mod@pp$X 的第二列：

mod <- lmer(Sepal.Length~Petal.Length + offset(Petal.Width*1) + (1|Species),data=iris)
Xold <-  Xnew <- mod@pp$X
set.seed(1)
Xnew[,2] <- rnorm(nrow(Xold))
mod@pp <- merPredD(X=Xnew, Zt=mod@pp$Zt, Lambdat=mod@pp$Lambdat, Lind=mod@pp$Lind, 
         theta=mod@pp$theta, n=nrow(mod@pp$X))
head(cbind(Xold[,2], mod@pp$X[,2]))

#######
  [,1]       [,2]
1  1.4 -0.6264538
2  1.4  0.1836433
3  1.3 -0.8356286
4  1.5  1.5952808
5  1.4  0.3295078
6  1.7 -0.8204684

我们还可以给mod@pp$X添加第三列：

mod <- lmer(Sepal.Length~Petal.Length + offset(Petal.Width*1) + (1|Species),data=iris)
Xnew <- as.matrix(data.frame(mod@pp$X, iris["Petal.Width"]))
colnames(Xnew) <- c(colnames(mod@pp$X),"Petal.Width")
mod@pp <- merPredD(X=Xnew, Zt=mod@pp$Zt, Lambdat=mod@pp$Lambdat, Lind=mod@pp$Lind, 
         theta=mod@pp$theta, n=nrow(mod@pp$X))
head(mod@pp$X)

#######
  (Intercept) Petal.Length Petal.Width
1           1          1.4         0.2
2           1          1.4         0.2
3           1          1.3         0.2
4           1          1.5         0.2
5           1          1.4         0.2
6           1          1.7         0.4

【讨论】：