在 lme4 中调整 object@pp$X

Adjusting object@pp$X in lme4

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

由于以下错误消息,我遇到的一个问题是更改 object@pp$X

错误:替换无效:引用 class 字段“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)!在它的函数中有两个调用:

  1. 是给object@beta(已经修改成功);
  2. 就是getME(object,"X"):(https://github.com/lme4/lme4/blob/master/R/lmer.R#L1932).

我愿意接受使用变量名称获取 fixef 系数的其他方法(能够直接调整 lmerMod)

提前致谢!

对象mod中槽pp的class为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"

class 个 lme4 个对象可以通过 lme4merPredD 命令生成。
这是使用 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