如何使用估算数据获取合并模型的 R^2、F 统计量和 p 值？答案

【问题标题】：How to get R^2, F-statistics, and p-value for pooled models with imputed data?如何使用估算数据获取合并模型的 R^2、F 统计量和 p 值？
【发布时间】：2022-11-19 19:08:18
【问题描述】：

我已经使用小鼠估算了具有估算数据的回归模型。

model1 <- with(imp, lm(outcome~ predictor1+ predictor2+ predictor3+ predictor4))).

在输出中我得到了一些信息

summary(pool(model1), conf.int = TRUE)

例如估计值、标准误差和 p 值。现在我想知道整个模型的 F 值和 R^2。

对于 R^2，我发现了以下代码：pool.r.squared(model1)。但我仍在寻找代码来显示 F 值。有没有人有这方面的经验？

【问题讨论】：

标签： r regression r-mice

【解决方案1】：

我们通过对 anova 的 F 值进行平均得到的常规 F 统计量，相比：

mean(anova(aov(bmi ~ hyp + chl, nhanes))[, 4], na.rm=TRUE)
summary(lm(bmi ~ hyp + chl, nhanes))$fstatistic[1]

对于合并分析，我们可以使用 miceadds::mi.anova 来获得 R^2 和 F 统计量。

library('miceadds')
nul <- capture.output(
  aov_fit <- miceadds::mi.anova(mi.res=imp, formula="bmi ~ hyp + chl" )
)

^{（不一定需要 capture.output，但可以防止控制台混乱。）}

所需信息现在存储在对象aov_fit 中。

aov_fit$r.squared  ## R-squared
# [1] 0.1158705

(fval <- mean(round(aov_fit$anova.table$`F value`, 2), na.rm=TRUE) ) ## F-statistic
# [1] 0.97

df_mod <- aov_fit$anova.table$df1[- nrow(aov_fit$anova.table)]  ## DF model
df_res <- el(fit$analyses)$df.residual  ## DF residual
c(df_mod, df_res)
# [1]  1  1 22

模型 p 值可以通过使用 F 分布pf() 的分布函数的右尾检验来计算。

pf(q=fval, df1=sum(df_mod), df_2=df_res, lower.tail=FALSE)  ## p-value
# [1] 0.3947152

我们现在可以使用 sprintf 来类似于 lm() 的 GOF 指标：

sprintf('Pooled R-squared: %s', round(aov_fit$r.squared, 4))
# [1] "Pooled R-squared: 0.1159"

tmp <- aov_fit$anova.table
sprintf('Pooled F-statistic: %s on %s and %s DF,  p-value: %s', 
        mean(round(tmp$`F value`, 2), na.rm=TRUE), 
        round(sum(tmp$df1[- nrow(aov_fit$anova.table)]), 2),
        round(el(fit$analyses)$df.residual, 2),
        format.pval(pf(fval, sum(df_mod), df_res, lower.tail=FALSE)))
# [1] "Pooled F-statistic: 0.97 on 2 and 22 DF,  p-value: 0.39472"

更新

得到r^2个_形容词，我们可以使用通常的公式，

adjR2 <- (r2, n, p) {
  1 - (n - 1)/(n - p - 1)*(1 - r2)
}

adjR2(aov_fit$r.squared, nrow(nhanes), sum(aov_fit$anova.table$df1, na.rm=TRUE))
# [1] 0.03549512

其中 n = 观察数，p = 参数数。

数据：

使用mice包的nhanes数据集。

library('mice')
set.seed(42)
imp <- mice(nhanes, m=100, printFlag=FALSE)
fit <- with(data=imp, exp=lm(bmi ~ hyp + chl))

【讨论】：

我测试了你的代码，它运行良好，谢谢。现在我想知道是否还有可能获得 F 统计量的 p 值，所以我知道整个模型是否显著地解释了我的数据中的一些差异。
@lanmi 谢谢你的好问题。请参阅 p 值的更新答案。
我有一个后续问题。是否也可以通过这种方式估计调整后的 R^2？
@lanmi 请参阅here如何调整 r^2。计算并更新我的答案。