添加随机效应会影响系数估计

机器算法验证 r 混合模式 随机效应模型
2022-03-28 10:20:59

我一直被教导随机效应只影响方差(误差),而固定效应只影响平均值。但我发现了一个例子,其中随机效应也会影响平均值 - 系数估计:

require(nlme)
set.seed(128)
n <- 100
k <- 5
cat <- as.factor(rep(1:k, each = n))
cat_i <- 1:k # intercept per kategorie
x <- rep(1:n, k)
sigma <- 0.2
alpha <- 0.001
y <- cat_i[cat] + alpha * x + rnorm(n*k, 0, sigma)
plot(x, y)

# simulate missing data
y[c(1:(n/2), (n*k-n/2):(n*k))] <- NA

m1 <- lm(y ~ x)
summary(m1)

m2 <- lm(y ~ cat + x)
summary(m2)

m3 <- lme(y ~ x, random = ~ 1|cat, na.action = na.omit)
summary(m3)

您可以看到模型的估计系数为x-0.013780 m1,而模型的估计系数为m30.0011713 - 两者都与零显着不同。

请注意,当我删除模拟缺失数据的行时,结果是相同的(它是完整的矩阵)。

这是为什么?

PS:请注意我不是专业的统计学家,所以如果你要回答很多数学问题,那么请给傻瓜做一些简单的总结:-)

2个回答

“我一直被教导随机效应只会影响方差(误差),而固定效应只会影响平均值。”

正如您所发现的,这仅适用于没有连续预测变量的平衡、完整(即没有缺失数据)数据集。换句话说,对于经典 ANOVA 文本中讨论的数据/模型类型。在这些理想情况下,可以独立估计固定效应和随机效应。

当这些条件不成立时(因为它们在“现实世界”中经常不成立),固定效应和随机效应就不是独立的。有趣的是,这就是为什么“现代”混合模型是使用迭代优化方法来估计的,而不是像经典的混合方差分析那样用一点矩阵代数来精确求解:为了估计固定效应,我们必须知道随机效应,但是为了估计随机效应,我们必须知道固定效应!与当前问题更相关的是,这也意味着当数据不平衡/不完整和/或模型中有连续的预测变量时,调整混合模型的随机效应结构可以改变模型固定部分的估计,反之亦然。

编辑 2016-07-05。来自评论:“您能否详细说明或引用为什么连续预测变量会影响模型固定部分的估计?

模型固定部分的估计值将取决于模型随机部分的估计值——即估计的方差分量——如果(但不仅是)预测变量的方差在集群中不同。如果任何预测变量是连续的(至少在“现实世界”数据中 - 理论上这可能不正确,例如在构建的数据集中),这几乎肯定是正确的。

在第一层,我认为你们都忽略了人口价值的收缩;混合效应模型中的每个受试者的斜率和截距比受试者内最小二乘估计值更接近总体估计值。 ” [参考。1]。以下链接可能也会有所帮助(对于我的混合模型,要查看哪些正确的描述?),请参阅 Mike Lawrence 的回答)。

此外,我认为您在玩具示例中有点不走运,因为您有一个完美平衡的设计,导致您在没有缺失值的情况下获得完全相同的估计。

现在尝试以下具有相同过程且没有缺失值的代码:

 cat <- as.factor(sample(1:5, n*k, replace=T) ) #This should be a bit unbalanced.
 cat_i <- 1:k # intercept per kategorie
 x <- rep(1:n, k)
 sigma <- 0.2
 alpha <- 0.001
 y <- cat_i[cat] + alpha * x + rnorm(n*k, 0, sigma) 

 m1 <- lm(y ~ x)  
 m3 <- lme(y ~ x, random = ~ 1|cat, na.action = na.omit) 

 round(digits= 7,fixef(m3)) ==  round(digits=7, coef(m1)) #Not this time lad.
 #(Intercept)           x 
 #      FALSE       FALSE 

现在在哪里,因为您的设计不是完全平衡的,所以您没有相同的系数估计值。

实际上,如果您以一种愚蠢的方式使用您的缺失值模式(例如:),y[ c(1:10, 100 + 1:10, 200 + 1:10, 300 + 1:10, 400 +1:10)] <- NA那么您的设计仍然是完美平衡的,您将再次获得相同的系数。

 require(nlme)
 set.seed(128)
 n <- 100
 k <- 5
 cat <- as.factor(rep(1:k, each = n))
 cat_i <- 1:k # intercept per kategorie
 x <- rep(1:n, k)
 sigma <- 0.2
 alpha <- 0.001
 y <- cat_i[cat] + alpha * x + rnorm(n*k, 0, sigma)
 plot(x, y)

 # simulate missing data in a perfectly balanced way
 y[ c(1:10, 100 + 1:10, 200 + 1:10, 300 + 1:10, 400 +1:10)] <- NA

 m1 <- lm(y ~ x)  
 m3 <- lme(y ~ x, random = ~ 1|cat, na.action = na.omit) 

 round(digits=7,fixef(m3)) ==  round(digits=7, coef(m1)) #Look what happend now...
 #(Intercept)           x 
 #       TRUE        TRUE 

您对原始实验的完美设计略有误导。当您将 NA 插入不平衡的距离时,您改变了个体受试者可以相互借用多少“强度”的模式。

简而言之,您看到的差异是由于收缩效应造成的,更具体地说,是因为您用非完美平衡的缺失值扭曲了原始的完美平衡设计。

参考 1:Douglas Bates lme4:使用 R 进行混合效果建模,第 71-72 页