ggeffects :: ggpredict（）not nlme :: lme（）返回人口级别的预测间隔

Question

我正在尝试从ggeffects获取总体级别预测间隔（PI）：ggpredict（）使用type =“ re”从nlme：lme（）< /代码>模型。 GGPREDICT并未返回LME（）模型的预期数据，而等效LMER（）模型正常。我的数据是自相关的重复度量，因此我需要lme（）带有colorelation = corar1（）。

我不确定这是错误，还是我只是想做一些我使用的工具没有设计的工具？

library(lme4)
library(nlme)
library(ggeffects)

Data <- data.frame(
  Subject = factor(c(3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 
                     4, 4, 4, 4, 4, 4, 4, 4, 4, 
                     5, 5, 5, 5, 5, 5, 5, 5, 5, 
                     6, 6, 6, 6, 6, 6, 6, 6, 
                     7, 7, 7, 7, 7, 7, 7, 7, 
                     8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 
                     9, 9, 9, 9, 9, 9, 9, 9, 
                     13, 13, 13, 13, 13, 13, 13, 13, 
                     14, 14, 14, 14, 14, 14, 14, 14, 14, 
                     19, 19, 19, 19, 19, 19, 19)), 
  x = c(20.0, 28.5, 38.0, 47.5, 57.0, 66.5, 76.0, 85.5, 95.0, 100.0, 
           21.0, 31.5, 42.0, 53.0, 63.0, 73.5, 84.0, 95.0, 100.0, 
           20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 
           22.0, 33.0, 44.0, 56.0, 67.0, 78.0, 89.0, 100.0, 
           21.5, 32.0, 43.0, 54.0, 65.0, 76.0, 86.5, 100.0, 
           20.0, 29.0, 38.5, 48.5, 58.0, 67.5, 77.0, 87.0, 96.5, 100.0, 
           23.0, 33.0, 44.0, 56.0, 67.0, 78.0, 89.0, 100.0, 
           23.5, 34.5, 46.5, 57.5, 69.5, 80.5, 92.5, 100.0, 
           20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 
           25.0, 37.5, 50.0, 62.5, 75.0, 87.5, 100.0)/100,
  y = c(1.10, 1.00, 1.25, 1.60, 1.40, 1.20, 2.50, 4.60, 6.80, 10.40, 
          0.90, 1.00, 0.75, 0.90, 1.10, 1.70, 4.35, 9.95, 11.45, 
          1.20, 0.70, 1.30, 1.40, 0.70, 1.25, 2.30, 4.30, 8.20, 
          1.55, 1.15, 0.95, 1.10, 1.90, 3.25, 7.20, 14.30, 
          1.85, 2.00, 1.70, 2.00, 2.35, 3.30, 7.30, 12.10, 
          2.20, 1.95, 1.15, 1.55, 1.65, 3.00, 4.45, 9.05, 13.75, 15.85, 
          1.55, 1.20, 1.35, 1.60, 1.65, 4.70, 6.45, 10.80, 
          1.00, 0.90, 1.00, 1.10, 1.60, 3.60, 8.05, 12.30, 
          0.85, 1.00, 1.05, 1.00, 1.35, 2.00, 3.65, 6.75, 13.10, 
          2.25, 2.35, 2.40, 2.80, 4.90, 8.15, 13.50)
)

Model.lme4 <- lmer(
  y ~ x + (1 | Subject),
  data = Data
)

# first running lme() without autocorrelation
Model.nlme <- lme(
  fixed = y ~ x,
  random = ~ 1 | Subject, 
  data = Data,
)

# Expected data return fine from the lmer() model:
ggpredict(
  Model.lme4,
  terms = c("x [all]"),
  type = "re",
)
# Predicted values of y
# 
#    x | Predicted |         95% CI
# ---------------------------------
# 0.20 |     -1.09 | [-5.93,  3.74]
# 0.28 |     -0.08 | [-4.89,  4.73]
# 0.38 |      0.99 | [-3.80,  5.78]
# 0.47 |      2.06 | [-2.71,  6.84]
# 0.58 |      3.38 | [-1.39,  8.14]
# 0.67 |      4.51 | [-0.26,  9.28]
# 0.77 |      5.70 | [ 0.92, 10.48]
# 1.00 |      8.44 | [ 3.61, 13.27]
# 
# Adjusted for:
# * Subject = 0 (population-level)
# 
# Intervals are prediction intervals.

# When run on the lme() model, predicted values & PIs are missing:
ggpredict(
  Model.nlme,
  terms = c("x [all]"),
  type = "re",
)
# Predicted values of y
# 
#    x
# ----
# 0.20
# 0.28
# 0.38
# 0.47
# 0.58
# 0.67
# 0.77
# 1.00
# 
# Adjusted for:
# * Subject = 0 (population-level)
# 
# Intervals are prediction intervals.

如果我使用相关= corar1（）它会产生与上面相同的结果。 如果我明确调用enter = c（“ x [all]”，“ object [0]”）

当我添加预期的自相关结构时，我会得到预测＆amp; pi值，但仅适用于主题因素的第一级：

Model.nlme <- lme(
  fixed = y ~ x,
  random = ~ 1 | Subject, 
  correlation = corAR1(form = ~ x | Subject),
  data = Data,
)

ggpredict(
  Model.nlme,
  terms = c("x [all]"),
  type = "re",
)
# Predicted values of y
# 
#    x | Predicted |         95% CI
# ---------------------------------
# 0.20 |      1.64 | [-1.44,  4.71]
# 0.28 |      2.82 | [-0.23,  5.87]
# 0.38 |      4.07 | [ 1.03,  7.12]
# 0.47 |      5.33 | [ 2.29,  8.36]
# 0.58 |      6.86 | [ 3.82,  9.90]
# 0.67 |      8.18 | [ 5.13, 11.24]
# 0.77 |      9.58 | [ 6.50, 12.66]
# 1.00 |     12.78 | [ 9.61, 15.95]
# 
# Adjusted for:
# * Subject = 3
# 
# Intervals are prediction intervals.

我在某个地方犯错误吗？还是有更好的方法来获取我想要的PI？谢谢！

Answer 1

不幸的是，我认为您可能会被困。

潜在的问题是broom.mixed :: tidy不返回布鲁普斯的标准偏差：

library(broom.mixed)
tidy(Model.lme4, effects = "ran_vals") ## includes SDs
tidy(Model.nlme, effects = "ran_vals") ## doesn't

但是，这并不是真正的bramom.mixed的错误。问题在于，目前nlme软件包没有提供任何返回BLUPS标准偏差的选项（请参阅Eg [电子邮件＆nbsp; prected] 线程，它讨论了这个问题，但没有提供解决方案的情况下出现。

为了解决这个问题，我可以说，您（或某人）必须深入研究nlme的胆量，并弄清楚如何为蓝图构造SDS。如果我这样做，我将从等式56-59，看看是否可以从lme fit提取相应的组件。特别是，有关nlme：eq 2.22，p。 71，给出蓝图的表达式，希望一个人可以将其与lme code 和一起得出SDS ...

或者：您可以在glmmtmb中拟合AR1模型，该模型应为整理 -able/GGPREDICT -able。代码如下所示；我认为这些是等效的模型，但是您上面给出的示例数据不支持拟合AR1模型（我们获得单一拟合），因此我不确定。

## define a per-subject observation-number variable (factor)
## alt. tidyverse: Data |> group_by(Subject) |> mutate(obs = seq(n())) |>
##                  mutate(across(obs, factor))
Data2 <- (Data
    |> split(Data$Subject)
    |> lapply(FUN = \(x) transform(x, obs = seq(x)))
    |> do.call(what = rbind)
    |> transform(obs = factor(obs))
)
Model.nlme.acf <- update(Model.nlme,
                         data = Data2,
                         correlation = corAR1(form = ~obs))
library(glmmTMB)
Model.glmmTMB <-  glmmTMB(y ~ x + (1|Subject) + ar1(0 + obs|Subject),
                          data = Data2)