Still, "evidence" does not indicate "clinical importance". You should take these two into account (not a new insight, I know...).
Daniel 🕹️
@strengejacke.de
628 followers 113 following 522 posts
He/she/it - 's' muss mit. We're lower than the world! R easystats project: easystats.github.io/easystats
![library(modelbased)
data(mtcars)
m <- glm(am ~ vs, family = binomial, data = mtcars)
estimate_means(
m,
"vs",
transform = effectsize::probs_to_odds,
predict = "response",
post_process = difference ~ pairwise
)
#> Post-processing difference ~ pairwise...
#> Estimated Marginal Means
#>
#> Parameter | Probability | SE | 95% CI | z
#> ------------------------------------------------------
#> (b2) - (b1) | 0.50 | 0.17 | [0.16, 0.84] | 2.88
#>
#> Variable predicted: am
#> Predictors modulated: vs
#> Predictions are on the response-scale.](https://cdn.bsky.app/img/feed_thumbnail/plain/did:plc:d73oqmrkqkaoqdjbjsmkhxp2/bafkreib7kyvkl2kd3fns4in64i232h3literyaxjdltn5pzivoxddglpsq)
![m <- lm(mpg ~ wt + hp, data = mtcars)
ci <- 0.95
dof <- insight::get_df(m)
set.seed(123)
crit <- mvtnorm::qmvt(ci, df = dof, tail = "both.tails", corr = cov2cor(vcov(m)))$quantile
# 95% level%
confint(m)
#> 2.5 % 97.5 %
#> (Intercept) 33.95738245 40.49715778
#> wt -5.17191604 -2.58374544
#> hp -0.05024078 -0.01330512
parameters::model_parameters(m)
#> Parameter | Coefficient | SE | 95% CI | t(29) | p
#> ----------------------------------------------------------------------
#> (Intercept) | 37.23 | 1.60 | [33.96, 40.50] | 23.28 | < .001
#> wt | -3.88 | 0.63 | [-5.17, -2.58] | -6.13 | < .001
#> hp | -0.03 | 9.03e-03 | [-0.05, -0.01] | -3.52 | 0.001
#>
#> Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed
#> using a Wald t-distribution approximation.
# sup-t adjustment
confint(m, level = 1 - 2 * pt(-abs(crit), df = dof))
#> 1.04 % 98.96 %
#> (Intercept) 33.31930882 41.135231417
#> wt -5.42443900 -2.331222484
#> hp -0.05384452 -0.009701374
set.seed(123)
parameters::model_parameters(m, p_adjust = "sup-t")
#> Parameter | Coefficient | SE | 95% CI | t(29) | p
#> ----------------------------------------------------------------------
#> (Intercept) | 37.23 | 1.60 | [33.32, 41.14] | 23.28 | < .001
#> wt | -3.88 | 0.63 | [-5.42, -2.33] | -6.13 | < .001
#> hp | -0.03 | 9.03e-03 | [-0.05, -0.01] | -3.52 | 0.003
#>
#> p-value adjustment method: Simultaneous confidence bands
#>
#> Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed
#> using a Wald t-distribution approximation.](https://cdn.bsky.app/img/feed_thumbnail/plain/did:plc:d73oqmrkqkaoqdjbjsmkhxp2/bafkreihkibqoi7acbtdg36cxu4dvi6hernwdstzbm7u3t73szzxkebpv74)
![m <- lm(mpg ~ wt + hp, data = mtcars)
x <- 0.95
l <- mvtnorm::qmvnorm(x, tail = "both.tails", corr = cov2cor(vcov(m)))$quantile
1 - 2 * pnorm(-abs(l))
#> [1] 0.9797365
# 95% level%
confint(m)
#> 2.5 % 97.5 %
#> (Intercept) 33.95738245 40.49715778
#> wt -5.17191604 -2.58374544
#> hp -0.05024078 -0.01330512
# new "95%" level
confint(m, level = 1 - 2 * pnorm(-abs(l)))
#> 1.01 % 98.99 %
#> (Intercept) 33.30014759 41.154392645
#> wt -5.43202222 -2.323639269
#> hp -0.05395274 -0.009593154](https://cdn.bsky.app/img/feed_thumbnail/plain/did:plc:d73oqmrkqkaoqdjbjsmkhxp2/bafkreifg3dyowdtqb5h4sl3mwqcdaj7owcgpxpi2awun7dl52fyglly5zm)
![library(tinyplot)
tinytheme(
"clean2",
palette.qualitative = c("black", "sienna", "grey", "indianred", "goldenrod")
)
hec <- as.data.frame(proportions(HairEyeColor, 2:3))
# Add a new factor level "Dont know" to the variable "Hair"
hec$Hair <- factor(hec$Hair, levels = c(levels(hec$Hair)[1:2], "Dont know", levels(hec$Hair)[3:4]))
hec$Hair[5] <- "Dont know"
tinyplot(
Freq ~ Eye | Hair,
facet = ~Sex,
data = hec,
type = "barplot",
center = TRUE,
flip = TRUE,
facet.args = list(ncol = 1),
yaxl = "percent"
)](https://cdn.bsky.app/img/feed_thumbnail/plain/did:plc:d73oqmrkqkaoqdjbjsmkhxp2/bafkreic6vsklld7cye5pzi272j5g5mcxtlkxy7hn42turg3hrihv6atkrq)
![library(tinyplot)
tinytheme(
"clean2",
palette.qualitative = c("black", "sienna", "grey", "indianred", "goldenrod")
)
hec <- as.data.frame(proportions(HairEyeColor, 2:3))
# Add a new factor level "Dont know" to the variable "Hair"
hec$Hair <- factor(hec$Hair, levels = c(levels(hec$Hair)[1:2], "Dont know", levels(hec$Hair)[3:4]))
hec$Hair[5] <- "Dont know"
tinyplot(
Freq ~ Eye | Hair,
facet = ~Sex,
data = hec,
type = "barplot",
center = TRUE,
flip = TRUE,
facet.args = list(ncol = 1),
yaxl = "percent"
)](https://cdn.bsky.app/img/feed_thumbnail/plain/did:plc:d73oqmrkqkaoqdjbjsmkhxp2/bafkreig6fb5o664yvdpvuonrxtkdzjaxvejscirwnepbim7lnrqavepr6y)
![data(heart, package = "glm2")
start.p <- sum(heart$Deaths) / sum(heart$Patients)
fit.glm <- glm(
cbind(Deaths, Patients - Deaths) ~
factor(AgeGroup) + factor(Severity) + factor(Delay) + factor(Region),
family = binomial(),
data = heart
)
fit.logbin <- logbin::logbin(
formula(fit.glm),
data = heart,
start = c(log(start.p), rep(c(0.2, 0.4), 4)),
trace = 1
)
modelbased::estimate_contrasts(fit.logbin, "Delay")
#> Marginal Contrasts Analysis
#>
#> Level1 | Level2 | Difference | SE | 95% CI | t(65) | p
#> -------------------------------------------------------------------
#> 2 | 1 | 0.01 | 0.01 | [-0.01, 0.04] | 0.86 | 0.394
#> 3 | 1 | 0.03 | 0.02 | [ 0.00, 0.07] | 2.13 | 0.037
#> 3 | 2 | 0.02 | 0.02 | [-0.01, 0.05] | 1.56 | 0.124
#>
#> Contrasts are on the response-scale (in %-points).
parameters::model_parameters(fit.logbin, exponentiate = TRUE, keep = "Delay")
#> Parameter | Risk Ratio | SE | 95% CI | Statistic | df | p
#> ---------------------------------------------------------------------
#> Delay [2] | 1.06 | 0.07 | [0.93, 1.22] | 0.85 | 65 | 0.395
#> Delay [3] | 1.19 | 0.10 | [1.01, 1.39] | 2.13 | 65 | 0.034
#>
#> Uncertainty intervals (profile-likelihood) and p-values (two-tailed)
#> computed using a Wald distribution approximation.](https://cdn.bsky.app/img/feed_thumbnail/plain/did:plc:d73oqmrkqkaoqdjbjsmkhxp2/bafkreifo4hz5vwaraulk5vob36cydw4a6s44pg3we4v3srvmoxmdz3f4xy)
![library(easystats)
# sample data
data(contrast_example, package = "modelbased")
# center predictor of interaction term, for consistent results
contrast_example$score <- center(contrast_example$score)
# model with default factor-contrasts
model1 <- lm(outcome ~ score * tx, data = contrast_example)
# create custom contrasts for the "tx" factor
treat_vs_none <- c(-2/3, 1/3, 1/3)
short_vs_long <- c(0, -1/2, 1/2)
contrasts(contrast_example$tx) <- cbind(treat_vs_none, short_vs_long)
# model with custom factor-contrasts
model2 <- lm(outcome ~ score * tx, data = contrast_example)
# see parameter "tx [treat_vs_none]", value 1.94 - we got this results
# due to custom contrasts...
model_parameters(model2)
# instead of re-running the model each time with different contrasts, simply
# formulate what you want to do as "comparison" (or "hypothesis" in
# {marginaleffects}): we want to compare the average effect of short + long
# treatment against no treatment. We find these coefficients in rows 2+3
# (short+long) and row 1 (none).
estimate_means(model1, "tx")
# we now formulate this as comparison
estimate_contrasts(model1, "tx", comparison = "((b2+b3)/2) = b1")
# marginaleffects equivalent
marginaleffects::avg_predictions(
model1,
variables = "tx",
hypothesis = "((b2+b3)/2) = b1"
)](https://cdn.bsky.app/img/feed_thumbnail/plain/did:plc:d73oqmrkqkaoqdjbjsmkhxp2/bafkreigpi2rscuvpccfdzw6cz523htggdkiheahg2woyb4hejl74lsyx3q)
![library(modelbased)
# sample data and model
data(contrast_example, package = "modelbased")
model <- lm(outcome ~ score * tx, data = contrast_example)
# custom contrasts for factor "tx"
cond_tx <- cbind("no treatment" = c(1, 0, 0), "treatment" = c(0, 0.5, 0.5))
# comparison of slopes for factor "tx", taking factor-contrasts into account
estimate_slopes(model, "score", by = "tx", comparison = cond_tx)
#> Estimated Marginal Effects
#>
#> Parameter | Slope | SE | 95% CI | t(24) | p
#> -----------------------------------------------------------
#> no treatment | 0.76 | 0.27 | [ 0.21, 1.31] | 2.86 | 0.009
#> treatment | 0.30 | 0.22 | [-0.15, 0.75] | 1.37 | 0.184
#>
#> Marginal effects estimated for score
# marginaleffects-equivalent
marginaleffects::avg_slopes(
model,
variables = "score",
by = "tx",
hypothesis = cond_tx
)
#>
#> Term Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
#> no treatment 0.763 0.267 2.86 0.00429 7.9 0.239 1.287
#> treatment 0.301 0.220 1.37 0.17119 2.5 -0.130 0.732
#>
#> Type: response](https://cdn.bsky.app/img/feed_thumbnail/plain/did:plc:d73oqmrkqkaoqdjbjsmkhxp2/bafkreie2s3hilk6ewpc73s6y5nqlkg7qxi2hmkk62vuuakyhf3tyobe53u)
