1

我正在尝试使用 lme4 拟合具有随机效应的泊松回归模型。我的响应变量 Y 表示双向表中的频率，但我只对协变量对交互的影响（表 Y 的对数线性表示的 alpha_aj 项）感兴趣，因此设计矩阵包括每一行的固定效应和在我的情况下不感兴趣的列。

我想将特定于单元格的随机效应与特定于列的差异包括在内。但是，当使用 glmer 时，我在模型中声明的第一列的估计方差分量始终为零，就好像存在某种多重共线性一样。我不明白为什么会发生这种情况，特别是因为当只假设一个共同方差时，预测的随机效应似乎不受任何限制。任何帮助将非常感激。

我的数据的一个例子是：

A <- 26 #number of rows 
J <- 3 #number of columns
y <- c(3117, 2426, 3048, 2830, 5151, 5399, 3331, 4631, 4419, 4434, 4356, 2161,
     6057, 3154, 5250, 4209, 2880, 4193, 5017, 7642, 3977, 5605, 4362, 5785,
     4044, 7259, 345,  371, 324,  393,  406,  303,  253,  284,  627,  520, 
     296,  428,  589,  381,  723,  506, 382,  825, 481,  572,  446,  486,
     603,  682,  526,  549, 1967, 1803, 1811, 1924, 1778, 1508, 1387, 1325,
     2259, 1737, 1209, 1699, 2247, 1463, 2548, 1469, 1504, 2610, 1703, 2444,
     1984, 1742, 2365, 2653, 1984, 2749)
x <- c(0.02, -0.26,  0.00, -0.20,  0.15,  0.21,  0.18,  0.22, -0.05,  0.06,
     0.28, -0.19, -0.02,  0.07, -0.15,  0.05, -0.04, -0.19,  0.13,  0.06,
     0.02,  0.08, -0.19, -0.19, -0.08,  0.05, -0.17,  0.13, -0.12,  0.10,
     -0.12, -0.05, -0.11, -0.06,  0.07, -0.01,  0.01,  0.06,  0.05, -0.13,
     0.15,  0.00, -0.09,  0.09, -0.04, -0.03, -0.05,  0.00,  0.13,  0.15,
     0.04,  0.01,  0.15,  0.14,  0.12, 0.10, -0.03, -0.17, -0.06, -0.16,
     -0.02, -0.05, -0.29,  0.13, -0.03,  0.06, 0.00, -0.05,  0.13,  0.11,
     -0.09, -0.03,  0.03, -0.08,  0.06,  0.05,  0.04, -0.06)

表 Y 按列堆叠。要定义我使用的设计矩阵，

zg <- kronecker(matrix(1,J,1),diag(A))  #area specific fixed effects
zl <- rbind(kronecker(diag(2),matrix(1,A,1)),matrix(-1,A,J-1)) #column specific fixed effects
z <- cbind(zg,zl,x)

具有共同方差的模型可以使用

library(lme4)
id <- seq(1,A*J)
model1 <- glmer(y ~ -1+z + (1|id),family = "poisson")
refit(model1)

为了使模型与特定于列的方差分量相匹配，我正在使用：

col <- kronecker(diag(3),matrix(1,A,1))
model2 <- glmer(y ~ -1+z +  
        (-1+col[,1]|id) + (-1+col[,2]|id)+(-1+col[,3]|id) ,family = "poisson")
refit(model2)

我得到结果：

Random effects:
 Groups Name     Std.Dev.
 id     col[, 1] 0.0000  
 id.1   col[, 2] 0.1973  
 id.2   col[, 3] 0.1231

如您所见，col1 的估计方差分量为零。如果我更改用于声明预测变量的顺序，

model3 <- glmer(y ~ -1+z +  
        (-1+col[,2]|id) + (-1+col[,3]|id)+(-1+col[,1]|id) ,family = "poisson")
refit(model3)

现在是第二列，它的方差分量估计为零：

Random effects:
 Groups Name     Std.Dev.
 id     col[, 2] 0.0000  
 id.1   col[, 3] 0.1422  
 id.2   col[, 1] 0.1861

从模型的定义来看，我不明白为什么我不能为三列提供特定于列的差异。我想我错过了一些东西，提前感谢任何提示！

已编辑（2）！！！

使用初始数据，我试图了解是否有任何理由应该限制方差。我使用 MCMCglmm 拟合模型：

model0 <- glm(y ~ -1+z ,family = "poisson")
prior <- list(B = list(mu=model0$coefficients,V=summary(model0)$cov.unscaled))
data <- data.frame(y,z,col,id)
data$column <- data$X1+2*data$X2+3*data$X3

library(MCMCglmm)
model3  <- MCMCglmm(y~-1+V1+V2+V3+V4+V5+V6+V7+V8+V9+V10+V11+V12+V13+V14+V15+V16+
        V17+V18+V19+V20+V21+V22+V23+V24+V25+V26+V27+V28+x, 
        random= ~idh(1):id, data = data, family="poisson", 
        rcov = ~idh(1):column,
        pr=TRUE, prior=prior,
verbose=TRUE, DIC=TRUE, singular.ok=FALSE, saveX=TRUE,
saveZ=TRUE, saveXL=TRUE, slice=FALSE, ginverse=NULL, 
nitt=50000, thin=10, burnin=10000)

std.dev <- c(sqrt(sum((model3$Sol[,-c(1:(A+J))]^2)[,1:A])/(nrow(model3$Sol)*A)),
sqrt(sum((model3$Sol[,-c(1:(A+J))]^2)[,(A+1):(2*A)])/(nrow(model3$Sol)*A)),
sqrt(sum((model3$Sol[,-c(1:(A+J))]^2)[,(2*A+1):(3*A)])/(nrow(model3$Sol)*A)))

估计的特定列方差分量为：

[1] 0.05413170 0.13366658 0.08073213

我仍然不明白为什么在与 lme4 拟合时不应该使用特定于列的差异。我认为在拟合模型之前会检查 lme4 是否正在执行，最终将方差分量设置为零。有人知道在 lme4 中指定这种模型的另一种方法吗？

非常感谢任何提示！

4

0 回答 0