缺少混合模型和GGEFFECT的标准错误和信心

Question

我正在尝试使用ggeffects :: ggpredict为我的模型制作一些效果图。我发现许多结果都缺少标准错误和置信度限制。我可以通过一些模拟数据来重现问题。似乎是专门用于观察到标准错误将预测概率接近0或1的观察结果。

我试图在链接刻度上获得预测，以诊断是否是从链接到响应的翻译问题的问题，但我不相信这是由包裹支持。

有什么想法解决这个问题吗？非常感谢。

library(tidyverse)
library(lme4)
library(ggeffects)

# number of simulated observations
n <- 1000

# simulated data with a numerical predictor x, factor predictor f, response y
# the simulated effects of x and f are somewhat weak compared to the noise, so expect high standard errors
df <- tibble(
  x = seq(-0.1, 0.1, length.out = n),
  g = floor(runif(n) * 3),
  f = letters[1 + g] %>% as.factor(),
  y = pracma::sigmoid(x + (runif(n) - 0.5) + 0.1 * (g - mean(g))),
  z = if_else(y > 0.5, "high", "low") %>% as.factor()
)

# glmer model
model <- glmer(z ~ x + (1 | f), data = df, family = binomial)
print(summary(model))
#> Generalized linear mixed model fit by maximum likelihood (Laplace
#>   Approximation) [glmerMod]
#>  Family: binomial  ( logit )
#> Formula: z ~ x + (1 | f)
#>    Data: df
#> 
#>      AIC      BIC   logLik deviance df.resid 
#>   1373.0   1387.8   -683.5   1367.0      997 
#> 
#> Scaled residuals: 
#>     Min      1Q  Median      3Q     Max 
#> -1.3858 -0.9928  0.7317  0.9534  1.3600 
#> 
#> Random effects:
#>  Groups Name        Variance Std.Dev.
#>  f      (Intercept) 0.0337   0.1836  
#> Number of obs: 1000, groups:  f, 3
#> 
#> Fixed effects:
#>             Estimate Std. Error z value Pr(>|z|)    
#> (Intercept)  0.02737    0.12380   0.221    0.825    
#> x           -4.48012    1.12066  -3.998 6.39e-05 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Correlation of Fixed Effects:
#>   (Intr)
#> x -0.001

# missing standard errors
ggpredict(model, c("x", "f")) %>% print()
#> Data were 'prettified'. Consider using `terms="x [all]"` to get smooth plots.
#> # Predicted probabilities of z
#> 
#> # f = a
#> 
#>     x | Predicted |       95% CI
#> --------------------------------
#> -0.10 |      0.62 | [0.54, 0.69]
#>  0.00 |      0.51 |             
#>  0.10 |      0.40 |             
#> 
#> # f = b
#> 
#>     x | Predicted |       95% CI
#> --------------------------------
#> -0.10 |      0.62 | [0.56, 0.67]
#>  0.00 |      0.51 |             
#>  0.10 |      0.40 |             
#> 
#> # f = c
#> 
#>     x | Predicted |       95% CI
#> --------------------------------
#> -0.10 |      0.62 | [0.54, 0.69]
#>  0.00 |      0.51 |             
#>  0.10 |      0.40 |
ggpredict(model, c("x", "f")) %>% as_tibble() %>% print(n = 20)
#> Data were 'prettified'. Consider using `terms="x [all]"` to get smooth plots.
#> # A tibble: 9 x 6
#>       x predicted std.error conf.low conf.high group
#>   <dbl>     <dbl>     <dbl>    <dbl>     <dbl> <fct>
#> 1  -0.1     0.617     0.167    0.537     0.691 a    
#> 2  -0.1     0.617     0.124    0.558     0.672 b    
#> 3  -0.1     0.617     0.167    0.537     0.691 c    
#> 4   0       0.507    NA       NA        NA     a    
#> 5   0       0.507    NA       NA        NA     b    
#> 6   0       0.507    NA       NA        NA     c    
#> 7   0.1     0.396    NA       NA        NA     a    
#> 8   0.1     0.396    NA       NA        NA     b    
#> 9   0.1     0.396    NA       NA        NA     c

^由

Answer 1

我认为这可能是由于单一模型拟合所致。

我深入研究了代码的核心 here ，似乎存在不匹配的地方预测协方差矩阵的维度 (3x3) 和预测值的数量 (15) 之间。

我进一步怀疑问题可能发生在此处 ：

rows_to_keep <- as.numeric(rownames(unique(model_matrix_data[
             intersect(colnames(model_matrix_data), terms)])))

也许该函数变得混乱，因为每个组的条件模式/BLUP 都是相同（一般来说，当随机效应方差为零时，这才是正确的）...？

这似乎值得在 ggeffects 问题列表 上提出一个问题？