library(tidyverse)
library(lme4)
library(marginaleffects)
# ?ChickWeight
# weight = body weight in grams
# Time = days since birth
# Chick = chick ID
# Diet = one of 4 possible diets
# This is a model similar to what you have
model <- lmer(weight ~ Time + I(Time^2) + Diet*Time + (1 | Chick), data = ChickWeight)
summary(model)
#> Linear mixed model fit by REML ['lmerMod']
#> Formula: weight ~ Time + I(Time^2) + Diet * Time + (1 | Chick)
#> Data: ChickWeight
#>
#> REML criterion at convergence: 5442.3
#>
#> Scaled residuals:
#> Min 1Q Median 3Q Max
#> -3.7484 -0.5503 -0.0059 0.5446 3.3252
#>
#> Random effects:
#> Groups Name Variance Std.Dev.
#> Chick (Intercept) 545.6 23.36
#> Residual 608.9 24.68
#> Number of obs: 578, groups: Chick, 50
#>
#> Fixed effects:
#> Estimate Std. Error t value
#> (Intercept) 41.04289 6.30890 6.506
#> Time 3.63070 0.60860 5.966
#> I(Time^2) 0.14567 0.02621 5.557
#> Diet2 -2.29496 10.47291 -0.219
#> Diet3 -12.67823 10.47291 -1.211
#> Diet4 -0.00622 10.48080 -0.001
#> Time:Diet2 1.84912 0.41691 4.435
#> Time:Diet3 4.66286 0.41691 11.184
#> Time:Diet4 2.93290 0.42224 6.946
#>
#> Correlation of Fixed Effects:
#> (Intr) Time I(T^2) Diet2 Diet3 Diet4 Tm:Dt2 Tm:Dt3
#> Time -0.413
#> I(Time^2) 0.272 -0.911
#> Diet2 -0.555 0.091 0.010
#> Diet3 -0.555 0.091 0.010 0.336
#> Diet4 -0.556 0.094 0.007 0.336 0.336
#> Time:Diet2 0.237 -0.230 -0.021 -0.423 -0.146 -0.146
#> Time:Diet3 0.237 -0.230 -0.021 -0.146 -0.423 -0.146 0.364
#> Time:Diet4 0.237 -0.239 -0.008 -0.144 -0.144 -0.423 0.359 0.359# The time coefficient is split across two different parts: Time and Time^2. You
# have to combine both parts somehow to get one overall slope. You can't
# interpret just one part on its own.
# You can see how the slope isn't consistent across time. It's shallower early
# on, then gets steeper as time progresses.
plot_predictions(model, condition = "Time")
# We can find the slope of that line at specific points, like 5, 10, 15, and 20
# (or any arbitrary points, really) If we feed hypothetical values into the
# model using datagrid(), marginaleffects will feed all the other necessary
# values (means or modes other coefficents, like Diet here)
slopes(model, newdata = datagrid(Time = c(5, 10, 15, 20)), variables = "Time")
#>
#> Term Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 % Diet Chick Time
#> Time 5.09 0.385 13.2 <0.001 129.8 4.33 5.84 1 1 5
#> Time 6.54 0.253 25.9 <0.001 487.1 6.05 7.04 1 1 10
#> Time 8.00 0.342 23.4 <0.001 399.1 7.33 8.67 1 1 15
#> Time 9.46 0.555 17.1 <0.001 214.2 8.37 10.54 1 1 20
#>
#> Columns: rowid, term, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, predicted, predicted_hi, predicted_lo, weight, Diet, Chick, Time
# ↑ Those are the coefficients that you want. At 5 days since birth, the effect
# of time is 5.09: a one-day increase in time is associated with a 5 gram
# increase in weight. At 15 days since birth, a one-day increase in time is
# associated with an 8 gram increase in weight; at 20 days, a one-day increase
# leads to 9 grams in weight, and so on. All those slopes are what you'd see in
# the plot of predictions:
plot_predictions(model, condition = "Time") +
geom_vline(xintercept = c(5, 10, 15, 20), linewidth = 0.5, color = "grey50")
# If you don't want to calculate the slope at specific points, you can find an
# overall average slope across the whole range of the line
avg_slopes(model, variables = "Time")
#>
#> Term Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
#> Time 8.7 0.153 56.8 <0.001 Inf 8.4 9
#>
#> Columns: term, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high
# ↑ Here the average marginal effect is 8.7. Across the average of days since
# birth, a one-day increase in time is associated with 8.7 more grams# The neat thing about this approach is that you can find marginal effects
# across all sorts of combinations of covariates. In your model you have
# MedianAge---you can use that to look at the slope of days_since for caes with
# a high median age and a low median age, something like this:
# slopes(
# Model_5,
# newdata = datagrid(
# days_since = c(10, 50, 100),
# Median_Age = c(30, 60)
# )
# )
# Or with the chick weights, we can look at the slopes for the different diets:
plot_predictions(model, condition = c("Time", "Diet")) +
geom_vline(xintercept = c(5, 10, 15, 20), linewidth = 0.5, color = "grey50")
slopes(model, newdata = datagrid(Time = c(5, 10, 15, 20), Diet = unique), variable = "Time")
#>
#> Term Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 % Chick Time Diet
#> Time 5.09 0.385 13.2 <0.001 129.8 4.33 5.84 1 5 1
#> Time 6.94 0.448 15.5 <0.001 176.9 6.06 7.82 1 5 2
#> Time 9.75 0.448 21.7 <0.001 345.8 8.87 10.63 1 5 3
#> Time 8.02 0.450 17.8 <0.001 233.8 7.14 8.90 1 5 4
#> Time 6.54 0.253 25.9 <0.001 487.1 6.05 7.04 1 10 1
#> Time 8.39 0.335 25.1 <0.001 458.6 7.74 9.05 1 10 2
#> Time 11.21 0.335 33.5 <0.001 814.2 10.55 11.86 1 10 3
#> Time 9.48 0.341 27.8 <0.001 562.9 8.81 10.15 1 10 4
#> Time 8.00 0.342 23.4 <0.001 399.1 7.33 8.67 1 15 1
#> Time 9.85 0.401 24.6 <0.001 441.3 9.07 10.64 1 15 2
#> Time 12.66 0.400 31.6 <0.001 726.6 11.88 13.45 1 15 3
#> Time 10.93 0.409 26.7 <0.001 520.0 10.13 11.74 1 15 4
#> Time 9.46 0.555 17.1 <0.001 214.2 8.37 10.54 1 20 1
#> Time 11.31 0.588 19.2 <0.001 271.0 10.15 12.46 1 20 2
#> Time 14.12 0.588 24.0 <0.001 420.3 12.97 15.27 1 20 3
#> Time 12.39 0.597 20.8 <0.001 315.6 11.22 13.56 1 20 4
#>
#> Columns: rowid, term, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, predicted, predicted_hi, predicted_lo, weight, Chick, Time, Diet
# ↑ This is some really rich detail! For chicks that are 5 days old on diet 1, a
# one-day increase in time is associated with a 5 gram increase in weight, but
# for chicks that are 5 days old on diet 3, a one-day increase is associated
# with a 9.75 gram increase in weight
# And finally, you can play with the re.form argument and set it to NULL or NA
# for conditional and marginal effects, where you control how to deal with the
# group-specific offsets in intercepts and slopes:
#
# Conditional effect = average chick; re.form = NA
# Marginal effect = chicks on average; re.form = NULL
#
# See this for more: https://www.andrewheiss.com/blog/2022/11/29/conditional-marginal-marginaleffects/
slopes(model, newdata = datagrid(Time = c(5, 10, 15, 20), Diet = unique), variable = "Time", re.form = NA)
#>
#> Term Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 % Chick Time Diet
#> Time 5.09 0.385 13.2 <0.001 129.8 4.33 5.84 1 5 1
#> Time 6.94 0.448 15.5 <0.001 176.9 6.06 7.82 1 5 2
#> Time 9.75 0.448 21.7 <0.001 345.8 8.87 10.63 1 5 3
#> Time 8.02 0.450 17.8 <0.001 233.8 7.14 8.90 1 5 4
#> Time 6.54 0.253 25.9 <0.001 487.1 6.05 7.04 1 10 1
#> Time 8.39 0.335 25.1 <0.001 458.6 7.74 9.05 1 10 2
#> Time 11.21 0.335 33.5 <0.001 814.1 10.55 11.86 1 10 3
#> Time 9.48 0.341 27.8 <0.001 562.9 8.81 10.15 1 10 4
#> Time 8.00 0.342 23.4 <0.001 399.1 7.33 8.67 1 15 1
#> Time 9.85 0.401 24.6 <0.001 441.3 9.07 10.64 1 15 2
#> Time 12.66 0.400 31.6 <0.001 726.6 11.88 13.45 1 15 3
#> Time 10.93 0.409 26.7 <0.001 520.0 10.13 11.74 1 15 4
#> Time 9.46 0.555 17.1 <0.001 214.2 8.37 10.54 1 20 1
#> Time 11.31 0.588 19.2 <0.001 271.0 10.15 12.46 1 20 2
#> Time 14.12 0.588 24.0 <0.001 420.4 12.97 15.27 1 20 3
#> Time 12.39 0.597 20.8 <0.001 315.6 11.22 13.56 1 20 4
#>
#> Columns: rowid, term, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, predicted, predicted_hi, predicted_lo, weight, Chick, Time, DietCreated on 2023-11-02 with reprex v2.0.2