Add a `confint` method to use t-distribution based quantiles
Closed this issue · 2 comments
Currently we implicitly use Gaussian distributed based quantiles. Which does not make really sense when we are making so much effort to derive adjusted d.f. for the t-statistics.
In particular, this can lead to inconsistencies, where the p-value could be larger than 0.05 while the confidence interval would exclude 0.
Example:
library(dplyr)
library(mmrm)
df <- fev_data |>
filter(AVISIT %in% c("VIS1", "VIS2"), USUBJID %in% paste0("PT", 1:50)) |>
droplevels()
mm2 <- mmrm(data = df, formula = log(FEV1) ~ AVISIT + us(AVISIT | USUBJID))
mm2_sum <- summary(mm2)
mm2_sum
# This is Gaussian based:
confint(mm2)
# We want actually here:
mm2_coefs <- as.data.frame(mm2_sum$coefficients) |>
mutate(t_quantile = stats::qt(0.975, df = df)) |>
mutate(
ci_lower = Estimate - t_quantile * `Std. Error`,
ci_upper = Estimate + t_quantile * `Std. Error`
)
mm2_coefs
Hopefully we don't need to make this as complex as https://github.com/glmmTMB/glmmTMB/blob/master/glmmTMB/R/methods.R#L865 in our case.
Relatedly, we also need to test whether emmeans
confidence intervals are correctly calculated following this update.
when using df_1d
, is it also appropriate to report the confidence intervals?
ideas is that we use df_1d
and df_md
for testing, for df_1d
the d.f., estimate and distribution also provided, but lack the CI. it might be helpful to also report the CI
basically maybe we can create a function like this
confint.mmrm_tmb <- function(...) { #for confidence intervals of the coefficients
# call df_1d to obtain the d.f., est, etc and summarize them into a data.frame
}
for differences in marginal means, use
confint(pairs(emmeans(fit)))
@clarkliming hm good question, let me think about it