How to get R^2, F-statistics, and p-value for pooled models with imputed data?

user_19836486 0
Viewed 54

I have estimated regression models with imputed data using mice.

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

In the output I get some information with

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

like estimate, standard error and p-value. Now I'd like to know the F-value and R^2 of the whole model.

For R^2 I found the following code: pool.r.squared(model1). But I'm still looking for a code to show the F-value. Does anyone have experience with that?

r-mice regression r
1 Answers

Conventional F-statistics we get by averaging the F-values from an anova, compare:

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

For pooled analyses, we may use miceadds::mi.anova to get both R^2 and F-statistic.

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

(The capture.output isn't necessarily needed but prevents the console from cluttering.)

The desired information is now stored in the object 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

The model p-value can be calculated by a right-tailed test using the distribution function for the F distribution pf().

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

We now could use sprintf to resemble somewhat the GOF metrics of lm():

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"

Data:

Using the nhanes data set of the mice package.

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