如何从R包小鼠运行的模型中获得拟合值
如何从R包小鼠运行的模型中获得拟合值
提问于 2019-09-02 09:04:07
我已经写了一个具有八个固定效果和两个随机效果的GLMM。我的两个修复效果包含缺失的数据，因此我使用R包MICE来估算缺失的值。
我想用我的模型中拟合的值和实际的观测值创建一个图表。如果我没有丢失数据，并且使用了lme4包来运行我的模型，我会简单地使用get ()函数来获得模型的拟合值。然而，由于我使用的是老鼠，我不确定如何为我的模型获得拟合值。当我使用函数the ()时，它返回“NULL”，而不是拟合值的向量。
我一直在互联网上搜寻，试图找到一个例子，在这个例子中，其他人在使用鼠标估计丢失的数据并运行GLMM后，获得了拟合值的向量，但我什么也找不到……
有没有人知道这个函数或从我的模型中计算拟合值的方法，该模型是使用MICE软件包运行的？或者推荐另一种可能对你有帮助的资源？
先谢谢你，奥利维亚
回答 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 rowsStack 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")页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接：
https://stackoverflow.com/questions/57754162
复制相关文章
相似问题
腾讯云开发者
