问如何从R包小鼠运行的模型中获得拟合值
EN

Stack Overflow用户

提问于 2019-09-02 09:04:07

回答 2查看 887关注 0票数 1

我已经写了一个具有八个固定效果和两个随机效果的GLMM。我的两个修复效果包含缺失的数据，因此我使用R包MICE来估算缺失的值。

我想用我的模型中拟合的值和实际的观测值创建一个图表。如果我没有丢失数据，并且使用了lme4包来运行我的模型，我会简单地使用get ()函数来获得模型的拟合值。然而，由于我使用的是老鼠，我不确定如何为我的模型获得拟合值。当我使用函数the ()时，它返回“NULL”，而不是拟合值的向量。

我一直在互联网上搜寻，试图找到一个例子，在这个例子中，其他人在使用鼠标估计丢失的数据并运行GLMM后，获得了拟合值的向量，但我什么也找不到……

有没有人知道这个函数或从我的模型中计算拟合值的方法，该模型是使用MICE软件包运行的？或者推荐另一种可能对你有帮助的资源？

先谢谢你，奥利维亚

r-mice

电商安全场景方案

帮助电商行业客户一站式解决Web安全、黑产羊毛党对抗、高危漏洞、BOT工具、恶意入侵门户服务器等安全问题

回答 2

Stack Overflow用户

发布于 2019-09-02 10:10:20

如果没有一个有效的例子，就很难弄清楚你到底遇到了什么。

然而，这里是一个使用mice和lme4包的例子，用lmer拟合一个(无意义的)模型。

require(mice)
require(lme4)

dt <- mice(nhanes2, seed = 314)

mod <- with(dt, lme4::lmer(bmi ~ chl + (1 | hyp)))

summary(pool(mod))

提供：

Class: mipo    m = 5 
               estimate         ubar            b            t dfcom       df       riv    lambda       fmi
(Intercept) 21.74573850 14.857649384 2.799244e+00 1.821674e+01    21 13.85172 0.2260850 0.1843959 0.2811936
chl          0.02574043  0.000379629 7.879048e-05 4.741775e-04    21 13.36442 0.2490552 0.1993949 0.2972420

将模型拟合到列表列中，并提取每个补偿集合的拟合值。然后取拟合值的平均值，作为汇集这些拟合值的一种方法。不过，我不确定这是否是推荐的拟合值池化方式。另请参阅：https://github.com/stefvanbuuren/mice/issues/82获取一些专家建议。

dt %>%
  mice::complete(action = "long", include = FALSE) %>%
  group_by(.imp) %>%
  nest(.key = dt) %>%
  mutate(mod = map(dt, ~ lmer(formula =  bmi ~ chl + (1 | hyp), data = .x))) %>%
  mutate(fitted = map(mod, ~ data.frame(fitted = fitted(.x), id = seq_along(fitted(.x))))) %>%
  select(.imp, fitted) %>%
  unnest() %>%
  group_by(id) %>%
  summarise(fitted = mean(fitted))

提供：

# A tibble: 25 x 2
      id fitted
   <int>  <dbl>
 1     1   27.1
 2     2   26.6
 3     3   26.6
 4     4   27.1
 5     5   24.7
 6     6   26.5
 7     7   24.8
 8     8   26.6
 9     9   27.9
10    10   27.0
# ... with 15 more rows

票数 1

Stack Overflow用户

发布于 2019-09-02 12:24:46

MICE-imputation方法的不确定性由多个产生的推算数据集来表示。作为回归的计算需要与“鲁宾规则”合并。在mice::pool(.)方法中只实现了简单的lm或glm方法。您可能需要自己编写一些代码来汇集，例如随机效果的计算，就像您可能使用lme4所做的那样。你可以在唐纳德·鲁宾( Rubin )的唐纳德·B( Donald B. )中找到所需的公式。纽约:威利，1987 on page 76。

但是，如果您的模型不是那么复杂，您可以将不同估算的拟合值组合到一个图中，并用颜色将它们分开。

示例1

library(mice)
iris.mice <- complete(mice(iris.mis), "long")
with(iris.mice, plot(Petal.Length, lm(Petal.Length ~ Sepal.Width + Petal.Width)$fitted,
                     type="n", xlab="imp.actual", ylab="imp.yhat", main="Petal.Length"))
by(iris.mice, iris.mice$.imp, function(x) {
  with(x, points(Petal.Length, lm(Petal.Length ~ Sepal.Width + Petal.Width, x)$fitted,
                 col=x$.imp))
})
legend("bottomright", legend=unique(iris.mice$.imp), pch=1, col=unique(iris.mice$.imp),
       ncol=3, title="Imp.")

另一种可能性是使用不同的估算方法，例如MissForest，其仅产生一个具有误差容限的估算数据集。您可以在绘图中以文本形式报告误差值。

示例2

library(missForest)
iris.imp <- missForest(iris.mis, xtrue=iris)
with(iris.imp$ximp, plot(Petal.Length, 
                         lm(Petal.Length ~ Sepal.Width + Petal.Width)$fitted,
                         xlab="imp.actual", ylab="imp.yhat", main="Petal.Length"))
text(5.5, 1.7, paste("NRMSE=", round(iris.imp$error[1], 2)))

数据

iris.mis <- structure(list(Sepal.Length = c(5.1, NA, 4.7, 4.6, NA, 5.4, 4.6, 
5, 4.4, 4.9, 5.4, 4.8, 4.8, 4.3, 5.8, 5.7, 5.4, 5.1, 5.7, 5.1, 
5.4, 5.1, 4.6, NA, 4.8, 5, 5, 5.2, 5.2, 4.7, 4.8, 5.4, NA, NA, 
4.9, NA, NA, 4.9, 4.4, 5.1, 5, 4.5, 4.4, 5, 5.1, 4.8, 5.1, 4.6, 
5.3, 5, 7, 6.4, 6.9, 5.5, NA, 5.7, 6.3, 4.9, 6.6, 5.2, 5, 5.9, 
6, 6.1, 5.6, 6.7, 5.6, 5.8, 6.2, 5.6, 5.9, 6.1, NA, 6.1, 6.4, 
6.6, 6.8, 6.7, 6, 5.7, 5.5, 5.5, 5.8, NA, 5.4, 6, 6.7, NA, 5.6, 
5.5, 5.5, 6.1, 5.8, NA, 5.6, 5.7, NA, 6.2, 5.1, NA, 6.3, NA, 
7.1, 6.3, 6.5, 7.6, 4.9, 7.3, 6.7, NA, 6.5, NA, 6.8, 5.7, 5.8, 
6.4, 6.5, 7.7, 7.7, 6, 6.9, 5.6, 7.7, 6.3, 6.7, 7.2, 6.2, 6.1, 
NA, 7.2, NA, 7.9, 6.4, 6.3, 6.1, 7.7, 6.3, NA, NA, 6.9, 6.7, 
6.9, 5.8, 6.8, 6.7, NA, NA, 6.5, 6.2, NA), Sepal.Width = c(3.5, 
3, 3.2, NA, 3.6, 3.9, 3.4, 3.4, NA, 3.1, 3.7, NA, NA, 3, NA, 
4.4, 3.9, 3.5, 3.8, 3.8, 3.4, NA, 3.6, 3.3, 3.4, 3, 3.4, 3.5, 
3.4, 3.2, 3.1, NA, 4.1, 4.2, 3.1, NA, 3.5, 3.6, 3, NA, 3.5, 2.3, 
3.2, NA, 3.8, NA, 3.8, NA, 3.7, 3.3, 3.2, NA, NA, 2.3, NA, 2.8, 
3.3, NA, 2.9, 2.7, 2, 3, 2.2, 2.9, 2.9, 3.1, 3, NA, 2.2, 2.5, 
3.2, NA, NA, 2.8, 2.9, 3, NA, NA, 2.9, 2.6, 2.4, 2.4, NA, 2.7, 
3, 3.4, 3.1, 2.3, 3, 2.5, NA, NA, 2.6, 2.3, 2.7, NA, 2.9, 2.9, 
2.5, 2.8, 3.3, 2.7, 3, 2.9, 3, 3, 2.5, 2.9, 2.5, 3.6, 3.2, 2.7, 
3, 2.5, NA, NA, 3, 3.8, 2.6, NA, 3.2, 2.8, 2.8, 2.7, 3.3, 3.2, 
2.8, 3, 2.8, 3, 2.8, 3.8, NA, 2.8, 2.6, NA, 3.4, 3.1, 3, 3.1, 
3.1, 3.1, 2.7, 3.2, 3.3, 3, NA, 3, NA, 3), Petal.Length = c(NA, 
1.4, NA, NA, 1.4, 1.7, 1.4, 1.5, NA, 1.5, 1.5, 1.6, NA, 1.1, 
1.2, 1.5, 1.3, 1.4, 1.7, 1.5, 1.7, 1.5, 1, 1.7, NA, 1.6, NA, 
1.5, NA, 1.6, 1.6, NA, 1.5, 1.4, 1.5, NA, NA, NA, 1.3, 1.5, 1.3, 
1.3, NA, 1.6, 1.9, 1.4, 1.6, 1.4, 1.5, 1.4, NA, 4.5, 4.9, 4, 
4.6, 4.5, NA, 3.3, 4.6, 3.9, NA, NA, 4, 4.7, 3.6, NA, 4.5, 4.1, 
4.5, 3.9, 4.8, 4, NA, NA, 4.3, 4.4, 4.8, 5, 4.5, 3.5, NA, 3.7, 
3.9, 5.1, NA, 4.5, NA, 4.4, 4.1, 4, 4.4, 4.6, NA, 3.3, 4.2, 4.2, 
4.2, 4.3, NA, NA, 6, 5.1, 5.9, NA, 5.8, 6.6, 4.5, 6.3, NA, 6.1, 
5.1, NA, 5.5, 5, 5.1, 5.3, 5.5, 6.7, 6.9, 5, 5.7, 4.9, 6.7, 4.9, 
5.7, 6, 4.8, 4.9, 5.6, 5.8, NA, 6.4, 5.6, 5.1, 5.6, 6.1, 5.6, 
5.5, 4.8, NA, NA, 5.1, NA, 5.9, 5.7, 5.2, 5, 5.2, 5.4, 5.1), 
    Petal.Width = c(0.2, 0.2, 0.2, 0.2, 0.2, NA, 0.3, 0.2, NA, 
    0.1, 0.2, NA, 0.1, 0.1, 0.2, 0.4, 0.4, 0.3, 0.3, 0.3, 0.2, 
    0.4, 0.2, NA, 0.2, NA, 0.4, 0.2, 0.2, 0.2, 0.2, 0.4, NA, 
    0.2, 0.2, NA, 0.2, 0.1, 0.2, NA, 0.3, 0.3, 0.2, 0.6, 0.4, 
    0.3, 0.2, 0.2, 0.2, 0.2, NA, 1.5, 1.5, NA, 1.5, 1.3, NA, 
    1, 1.3, 1.4, 1, 1.5, 1, NA, 1.3, 1.4, 1.5, 1, 1.5, 1.1, NA, 
    NA, 1.5, 1.2, 1.3, 1.4, 1.4, NA, 1.5, NA, NA, 1, 1.2, NA, 
    1.5, 1.6, 1.5, 1.3, 1.3, NA, 1.2, NA, NA, 1, 1.3, 1.2, 1.3, 
    1.3, 1.1, 1.3, 2.5, NA, 2.1, 1.8, 2.2, 2.1, 1.7, NA, 1.8, 
    NA, 2, 1.9, 2.1, 2, 2.4, 2.3, NA, NA, 2.3, 1.5, 2.3, 2, 2, 
    1.8, NA, 1.8, 1.8, NA, 2.1, NA, 1.9, 2, 2.2, 1.5, 1.4, 2.3, 
    2.4, 1.8, 1.8, 2.1, 2.4, 2.3, NA, 2.3, 2.5, 2.3, 1.9, NA, 
    NA, 1.8), Species = structure(c(1L, 1L, NA, 1L, NA, 1L, NA, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, NA, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, NA, 1L, 1L, 1L, NA, 1L, 
    1L, 1L, 1L, NA, 1L, NA, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, NA, 
    2L, 2L, NA, 2L, NA, 2L, 2L, 2L, 2L, 2L, 2L, NA, NA, 2L, 2L, 
    2L, 2L, 2L, 2L, 2L, NA, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
    2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, NA, 2L, 2L, 2L, NA, NA, 
    NA, 2L, 2L, 3L, NA, NA, NA, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
    3L, 3L, 3L, NA, 3L, NA, 3L, 3L, 3L, 3L, 3L, 3L, NA, 3L, NA, 
    3L, 3L, NA, 3L, NA, 3L, 3L, 3L, 3L, NA, 3L, 3L, 3L, NA, 3L, 
    3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L), .Label = c("setosa", "versicolor", 
    "virginica"), class = "factor")), row.names = c(NA, -150L
), class = "data.frame")