Daniel 🕹️

@strengejacke.de

624 Followers 111 Following

He/she/it - 's' muss mit. We're lower than the world! R easystats project: https://easystats.github.io/easystats/

set.seed(123) n <- 200 d <- data.frame( outcome = rnorm(n), grp = as.factor(sample(c("treatment", "control"), n, TRUE)), episode = as.factor(sample(1:3, n, TRUE)), sex = as.factor(sample(c("female", "male"), n, TRUE, prob = c(0.4, 0.6))) ) model2 <- lm(outcome ~ grp * episode, data = d) avg_predictions( model2, by = c("episode", "grp"), hypothesis = "(b1 - b3) = (b2 - b4)" )
library(modelbased) library(marginaleffects) set.seed(123) n <- 200 d <- data.frame( outcome = rnorm(n), grp = as.factor(sample(c("treatment", "control"), n, TRUE)), episode = as.factor(sample(1:3, n, TRUE)), sex = as.factor(sample(c("female", "male"), n, TRUE, prob = c(0.4, 0.6))) ) model2 <- lm(outcome ~ grp * episode, data = d) estimate_contrasts(model2, c("episode", "grp"), comparison = "(b1 - b2) = (b4 - b5)", estimate = "average") avg_predictions( model2, by = c("episode", "grp"), hypothesis = "(b1 - b3) = (b2 - b4)" )
library(parameters) library(modelbased) library(marginaleffects) library(lme4) data("qol_cancer") qol_cancer <- datawizard::demean(qol_cancer, select = c("phq4", "QoL"), by = "ID") mixed_1 <- lmer( QoL ~ time + phq4_within + phq4_between + (1 | ID), data = qol_cancer ) # within- and between effects model_parameters(mixed_1) # context effect - difference between within- and between-effect estimate_contrasts(mixed_1, c("phq4_within", "phq4_between"), comparison = "slope") # context effect - using marginal effects avg_comparisons( mixed_1, variables = c("phq4_within", "phq4_between"), hypothesis = ~pairwise )
# variance components parameters::model_parameters(m, effects = "random") # BLUPs parameters::model_parameters(m, effects = "grouplevel") # fixef+ranef parameters::model_parameters(m, effects = "random_total") # same as modelbased::estimate_grouplevel(m) # dot plot of BLUPs parameters::model_parameters(m, effects = "grouplevel") |> plot() # dot plot of BLUPs modelbased::estimate_grouplevel(m) |> plot()
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.
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
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" )
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.
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 ) #> logbin parameterisation 1/81 #> Deviance = 149.321 Iterations - 9052 modelbased::estimate_means(fit.glm, "AgeGroup") |> print(select = "minimal") modelbased::estimate_means(fit.logbin, "AgeGroup") |> print(select = "minimal")
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" )
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
library(modelbased) data("nhanes2", package = "mice") imp <- mice::mice(nhanes2, printFlag = FALSE) # estimated marginal means predictions <- lapply(1:5, function(i) { m <- lm(bmi ~ age + hyp + chl, data = mice::complete(imp, action = i)) estimate_means(m, "age") }) out <- pool_predictions(predictions) out plot(out)
library(easystats) library(glmmTMB) data(qol_cancer) d <- cbind(qol_cancer, demean(qol_cancer, "phq4", "ID")) set.seed(1) d$x <- rnorm(nrow(d)) m1 <- glmmTMB( QoL ~ time + age + phq4_within + phq4_between + x + (1 + time | ID) + (0 + phq4_within | ID), data = d ) #> Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model convergence #> problem; non-positive-definite Hessian matrix. See vignette('troubleshooting') check_convergence(m1) #> [1] FALSE prior1 <- data.frame( prior = "gamma(1, 2.5)", class = "ranef" ) m2 <- glmmTMB( QoL ~ time + age + phq4_within + phq4_between + x + (1 + time | ID) + (0 + phq4_within | ID), data = d, priors = prior1 ) check_convergence(m2) #> [1] TRUE compare_parameters(m1, m2, effects = "all") |> print_html()
library(parameters) library(glmmTMB) data(sleepstudy, package = "lme4") m1 <- glmmTMB( Reaction ~ Days + (1 + Days | Subject), data = sleepstudy ) prior1 <- data.frame( prior = "t(0,3,3)", class = "fixef", coef = "Days" ) m2 <- glmmTMB( Reaction ~ Days + (1 + Days | Subject), data = sleepstudy, priors = prior1 ) prior2 <- data.frame( prior = "t(0,3,3)", class = "ranef", coef = "Subject" ) m3 <- glmmTMB( Reaction ~ Days + (1 + Days | Subject), data = sleepstudy, priors = prior2 ) compare_parameters(m1, m2, m3, effects = "all") |> print_html()