Add a `confint` method to use t-distribution based quantiles

[danielinteractive]

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.

[clarkliming]

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)))

[danielinteractive]

@clarkliming hm good question, let me think about it