机器算法验证 - 广义多级模型中令人费解的预测值 - 吾爱随笔录

广义多级模型中令人费解的预测值

机器算法验证物流混合模式重复测量多层次分析

2022-03-23 22:33:07

我使用逻辑链接函数运行了一个多级模型，其结果是二分变量。我使用lme4R 中的包进行了此操作。可以在此 GitHub Gist 中找到我的数据的近似值：https ://gist.github.com/markhwhiteii/3a954c4da82efdad5d8abee601b8a6a8

数据收集过程的设计非常简单：参与者被随机分配到治疗 ( t) 或控制 ( c) 条件。然后他们被问到三个问题；他们可能会回答其中的一个、两个或三个。对于每个问题（q1、q2、q3），他们可以正面回答（1）或负面回答（0）。我不一定关心变量的影响——我将每个问题视为或多或少相同结构的指标。

问题很简单：治疗是否增加了积极响应的概率（1）？重要的是，我想使用的效应量是两组之间模型隐含的百分点差异。这需要将预测值从 logits 转换为概率。

我指定的模型是一个仅拦截模型，响应为 1 级（表示为 $i$ ) 和第 2 级的人（由 $j$ ):

$\text{logit}(\pi_{ij}) = \beta_{0j}$ 和 $\beta_{0j} = \gamma_{00} + \gamma_{01}X_j + u_{0j}$

在的语法中lme4，这是：

mod_1 <- glmer(y ~ x + (1 | id), dat, binomial)

（如果您read.csv()在上面 Gist 链接中的 .csv 文件中，此代码将估计模型。）

这是混乱上升的地方。我可以将控制条件和治疗条件中响应的预测值转换为概率空间。他们给我：

inv_logit <- function(x) exp(x) / (1 + exp(x))
(control <- inv_logit(fixef(mod_1)[[1]]))
(treatment <- inv_logit(sum(fixef(mod_1))))
treatment - control

这将显示在控制条件和治疗条件下积极响应的概率约为 0.0003，或为 0.03% 的概率是积极的。条件之间的差异非常小。

每种情况下的这些条件概率似乎太小了，尤其是当我们查看两种情况下积极响应的天真百分比时：

with(dat, tapply(y, x, mean))

这将表明在控制条件下大约 22% 的反应是阳性的，而在治疗条件下 23% 是阳性的。我的基本问题是：从经验上看，22% 左右的积极结果与模型暗示的只有 0.03% 之间有什么如此大的差异？

如果我们首先按人平均（不允许那些反应更多的人有更多的权重）然后按条件进行平均，我们会得到同样的结果：

library(dplyr)
dat %>% 
  group_by(id) %>% 
  mutate(p = mean(y)) %>% 
  slice(1) %>% 
  group_by(x) %>% 
  summarise(p = mean(p))

两种情况都显示大约 24% 的反应是积极的。同样，这与 .03% 所暗示的完全不同mod_1。

其他一些可能有用的信息：

604 人有 1 条回复，301 人有 2 条回复，95 人有 3 条回复。我认为部分问题是大多数人只有一个观察结果？不过，据我所知，我们不需要为每个集群（在本例中为一个人）估计不同的响应正态分布，所以我不明白为什么只有 1 个响应会给我带来困扰—— UPS。
随机截距的分布显然不正常，如下所示hist(ranef(mod_1)$id[[1]], breaks = 50)：

该模型显然似乎是从随机截距的可怕分布中错误指定的，但我想知道：

为什么我天真地认为大约 22-24% 的正面响应与模型隐含的 0.03% 相差如此之大？我正在寻找“模型指定错误”之外的答案，因为我很好奇为什么会发生这种情况，而且这并没有让我更接近指定合适的模型。这导致我...
如何指定一个合适的模型来回答我关于实验条件的问题？我想我也可以尝试使用某种类型的三明治协方差估计来处理相关响应，例如sandwich::vcovCL()？在当前情况下，这种方法不会过多地改变 SE。这也不能满足我的好奇心，为什么混合模型会给我奇怪的预测。

1个回答

答案可能至少部分隐藏在您的答案中：

预测的随机效应都是非负的。将它们添加到固定效应的线性预测器中会将平均预测推高到所需的水平。

那么 eBLUP 的非中心分布的原因可能是什么？我认为它隐藏在数据极度不平衡的情况下：1000 个 id 中有 604 个只提供了一个单一的测量值，这很容易导致随机效应和固定效应之间的数值冲突。

在这种情况下，我经常运行一个正常的混合模型作为敏感性分析。如果我们对您的数据执行此操作，我们会得到：

mod_2 <- lmer(y ~ x + (1 | id), dat)
summary(mod_2)

# Output
Random effects:
 Groups   Name        Variance Std.Dev.
 id       (Intercept) 0.06961  0.2638  
 Residual             0.10825  0.3290  
Number of obs: 1491, groups:  id, 1000

Fixed effects:
            Estimate Std. Error t value
(Intercept) 0.228315   0.022716  10.051
xt          0.007224   0.026922   0.268

Correlation of Fixed Effects:
   (Intr)
xt -0.844

现在，固定效应看起来与描述性分析所预期的一样，与 GLMER 发现的非常不同。现在，eBLUP 甚至居中（当然离正常还很远）：

summary(unlist(ranef(mod_2, drop = TRUE)))

# Output
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
-0.15513 -0.12844 -0.09218  0.00000  0.14878  0.50347

根据您的确切研究问题和分析目标，您可以尝试 R 中的广义估计方程 GEE：

library(gee)
fit_gee <- gee(y ~ x, 
               id = id,         
               data = dat %>% arrange(id), 
               family = binomial, 
               corstr = "exchangeable")
summary(fit_gee)

# Output
Coefficients:
               Estimate Naive S.E.    Naive z Robust S.E.   Robust z
(Intercept) -1.22561668  0.1257257 -9.7483411   0.1235148 -9.9228358
xt           0.04325299  0.1484306  0.2914022   0.1479156  0.2924167

Estimated Scale Parameter:  0.9859774
Number of Iterations:  3

Working Correlation
          [,1]      [,2]      [,3]
[1,] 1.0000000 0.3133096 0.3133096
[2,] 0.3133096 1.0000000 0.3133096
[3,] 0.3133096 0.3133096 1.0000000

# Distribution of predictions
summary(fitted(fit_gee))

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.2269  0.2269  0.2346  0.2324  0.2346  0.2346

其它你可能感兴趣的问题

上一篇当 X 服从一般正态分布时，min{0,X} 的分布是什么？下一篇关于形状位置的时间序列聚类的合适距离度量