Treatment effects with marginaleffects
Predictions, slopes and comparisons in a feeding experiment
Gavin Simpson
Aarhus University
A pig feeding experiment
In a simulated experiment, researchers randomly assign control , A or B to pens within indoor and outdoor housing. Each pen holds one pig. Researchers record starting weight and four-week gain, both in kg.
What gain does the model predict for an indoor pig starting at 30 kg?
How much would switching from control to diet A change predicted gain?
How quickly does predicted gain change near a starting weight of 30 kg?
library ("dplyr" )library ("marginaleffects" )
Opening, about 1 minute. Introduce the three questions in plain language. The first asks for an outcome, the second for a treatment contrast and the third for a local rate of change with a measured covariate. We can interpret the diet comparison in light of random assignment. Starting weight and housing are background characteristics in this example. Allow about 20 minutes in total.
Three questions for a model
What outcome does the model predict?
predictions()avg_predictions()plot_predictions()
How much do two settings differ?
comparisons()avg_comparisons()plot_comparisons()
How fast does the outcome change here?
slopes()avg_slopes()plot_slopes()
The rows describe the pigs or settings where we ask each question.
Pause on the difference between a finite change and a derivative. A finite change could switch a factor level or increase a numeric variable by a specified amount. A slope asks about a tiny change in a numeric variable at a particular point. For factors, ask for a comparison between named levels.
Simulated pen records
<- tidyr:: crossing (diet = factor (c ("control" , "A" , "B" ), levels = c ("control" , "A" , "B" )),housing = factor (c ("indoor" , "outdoor" ), levels = c ("indoor" , "outdoor" )),start_kg = seq (20 , 40 , by = 5 ),replicate = 1 : 8 |> filter (housing == "indoor" | replicate <= 2 )
The experiment includes eight indoor pens and two outdoor pens per combination of diet and starting weight.
About 2 minutes for setup and fitting. Diet and housing are factors with explicit levels. Crossing builds the full set of combinations. Replicate numbers index pens within each combination. Keep all eight pens per indoor combination and two per outdoor combination. The resulting 120 to 30 split is deliberate. Later it gives us two meaningful ways to average.
Simulating four-week gain
set.seed (148 )<- pens |> mutate (gain_kg = 8 + 0.20 * (start_kg - 30 ) - 0.018 * (start_kg - 30 )^ 2 + 1.2 * (diet == "A" ) + 1.8 * (diet == "B" ) + 2.0 * (diet == "A" & housing == "outdoor" ) - 0.8 * (diet == "B" & housing == "outdoor" ) + rnorm (n (), sd = 2.4 )
Expected gain follows a curve with starting weight. Diet and housing shift the curve. Each pen has its own random noise.
The simulated gain consists of a baseline curve, diet adjustments, extra adjustments outdoors for A and B, and random noise. The fitted model on the next slide estimates these patterns from the noisy observations.
Model for four-week gain
<- lm (~ diet * housing + start_kg + I ((start_kg - 30 )^ 2 ),data = pens
diet * housing lets diet differences vary with housing. The squared term lets the relationship with starting weight bend. The model treats the separate pens as independent.
This is an ordinary linear model for a continuous outcome with one record per independent pen. The response scale is kg of gain. The worked example assumes the model and its standard errors are suitable for these simulated data.
A grid of pig profiles
datagrid() crosses the supplied values. With its default settings, it fills in omitted numeric predictors at their means and omitted categorical predictors at their most common levels.
<- datagrid (model = fit,diet = distinct (pens, diet)$ diet,start_kg = 30
rowid housing diet start_kg
1 1 indoor control 30
2 2 indoor A 30
3 3 indoor B 30
The diet column keeps its fitted factor levels. Each row describes a 30 kg pig housed indoors.
About 2 minutes for the grid and the three row-level quantities. Point to the three rows and read one aloud. The housing column comes from the mode in the model data. We supplied the weight as 30 on purpose. When a grid is stored separately, pass the model explicitly. A datagrid call inside newdata can get the model from the enclosing function. Check that each row describes a realistic combination, especially with nested groups. See https://marginaleffects.com/man/r/datagrid.html
Predictions of weight gain
predictions (fit, newdata = grid, type = "response" )
diet start_kg Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
control 30 7.87 0.445 17.7 <0.001 230.5 7.00 8.74
A 30 8.94 0.445 20.1 <0.001 295.9 8.06 9.81
B 30 9.63 0.445 21.7 <0.001 343.2 8.76 10.51
Type: response
type = "response" gives kg of predicted gain. Each interval describes uncertainty about the mean gain for pigs in that setting. Individual pigs vary around that mean.
Read the row for diet A as expected four-week gain for an indoor pig starting at 30 kg. The three rows differ only in diet. In a binomial model the response scale could instead be a probability, so name the units when reporting a result. Prediction intervals for individual future observations answer a different question. See https://marginaleffects.com/man/r/predictions.html
Switching from control to A
<- c ("control" , "A" )comparisons (newdata = datagrid (start_kg = 30 ),variables = list (diet = diet_A),type = "response"
start_kg Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
30 1.06 0.535 1.99 0.047 4.4 0.0143 2.11
Term: diet
Type: response
Comparison: A - control
R predicts gain for A and for control at the same starting weight and housing, then subtracts A − control . The printed contrast gives the subtraction order.
The diet in the original row is replaced twice. Start by changing only diet to A, predict, then change only diet to control and predict again. Their difference should agree, up to rounding, with subtracting the corresponding two rows on the previous slide. Use the standard error for the difference to assess its precision, since both predictions use the same fitted model. See https://marginaleffects.com/man/r/comparisons.html
The slope at two weights
slopes (newdata = datagrid (start_kg = c (25 , 35 )),variables = "start_kg" ,type = "response"
start_kg Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
25 0.3097 0.0543 5.704 <0.001 26.4 0.203 0.4161
35 -0.0178 0.0543 -0.327 0.743 0.4 -0.124 0.0886
Term: start_kg
Type: response
Comparison: dY/dX
The slope is measured in kg of predicted gain per kg of starting weight . The two slopes can differ because the fitted curve bends.
Imagine laying a tangent along the gain curve at 25 kg, and then at 35 kg. Its steepness is the derivative. These numbers describe local rates. The difference in predicted gain between 25 and 35 kg would instead be a finite comparison over that whole interval. Starting weight was measured, so the slope describes a fitted association. See https://marginaleffects.com/man/r/slopes.html
Predicted gain by diet
Predicted gain. condition creates a grid over the named predictors. Other predictors use typical values unless you specify them.
plot_predictions (condition = c ("start_kg" , "diet" type = "response" + :: labs (x = "Starting weight (kg)" ,y = "Predicted gain (kg)"
About 3 minutes for the plots. Read the axes and trace one line. The lines describe the fitted means over starting weight in indoor housing, the most common housing level. The ribbons give pointwise uncertainty about each mean. Looking at a difference calls for a plot of the differences themselves. See https://marginaleffects.com/man/r/plot_predictions.html
Diet effects by housing
Here by averages the A − control differences separately over the observed pens in each housing group.
plot_comparisons (variables = list (diet = diet_Aby = "housing" ,type = "response" + :: labs (y = "A − control (kg)"
Read each point as an estimated average diet contrast for that housing group. The interval concerns that contrast. The question of whether the diet contrast differs between housing groups needs its own estimate and interval, shown later. See https://marginaleffects.com/man/r/plot_comparisons.html
The slope of predicted gain
Each point on the line gives the local slope at that starting weight. A zero slope means the fitted gain curve is flat there.
plot_slopes (variables = "start_kg" ,condition = "start_kg" ,type = "response" + :: labs (x = "Starting weight (kg)" ,y = "Slope (kg gain per kg)"
Connect this to the earlier gain curve. As the curve bends downward, its slope falls. The fitted slope is shared across diet and housing in this model. For each plot function, choose condition for a grid or by for averages over data rows. Use draw = FALSE if you want the values behind the plot. See https://marginaleffects.com/man/r/plot_slopes.html
Average effect across pens
For each observed pen , predict gain on A and on control at the same starting weight and housing. Subtract within each pen, then average the differences.
avg_comparisons (variables = list (diet = diet_A),type = "response"
Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
1.73 0.479 3.62 <0.001 11.7 0.796 2.67
Term: diet
Type: response
Comparison: A - control
With newdata omitted, each pen counts equally. Indoor pens account for 80% of the average.
About 4 minutes for averaging choices and factor contrasts. Emphasise that every pen contributes to both sides. The matched background gives the comparison its interpretation. An average of slopes or predictions also averages row-level results, but a slope or prediction is a different quantity from a diet comparison. Averaging individual standard errors would miss how the estimates share the same fitted coefficients. See https://marginaleffects.com/man/r/comparisons.html
A balanced housing average
avg_comparisons (newdata = "balanced" ,variables = list (diet = diet_A),type = "response"
Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
2.74 0.598 4.58 <0.001 17.7 1.57 3.91
Term: diet
Type: response
Comparison: A - control
The balanced grid crosses all diet and housing levels at mean starting weight. Each housing type contributes half of the average.
Compare this estimate with the previous one. Balancing also fixes starting weight at its mean. Here the starting-weight terms cancel in the A minus control contrast, so the change in estimate comes from the different mix of housing. In another model, such as a logistic model on the probability scale, setting a numeric predictor to its mean can change the contrast as well. The target population determines which average makes sense. See https://marginaleffects.com/man/r/datagrid.html
A counterfactual grid
A counterfactual grid copies the observed pens once per specified diet. Each copy keeps the original housing and starting weights.
<- datagrid (model = fit,diet = distinct (pens, diet)$ diet,grid_type = "counterfactual" count (all_diets, diet)
diet n
1 control 150
2 A 150
3 B 150
All three sets of 150 rows represent the same pens. Only the assigned diet changes between sets. A counterfactual grid lets us predict outcomes under each setting for the same population. See https://marginaleffects.com/man/r/datagrid.html
One predicted mean per diet
avg_predictions (newdata = all_diets,by = "diet" ,type = "response"
diet Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
control 6.94 0.339 20.5 <0.001 307.5 6.27 7.60
A 8.67 0.339 25.6 <0.001 478.2 8.01 9.33
B 8.48 0.339 25.0 <0.001 457.5 7.82 9.14
Type: response
by = "diet" averages each set of 150 predictions separately. Subtracting the control mean from the A mean gives the observed-pen average difference.
Show why the subtraction agrees: both means use the same pens with the same weights, and subtracting the means equals averaging the paired differences. In contrast, grouping the original rows by their observed diet asks for means in each actual treatment group, which can have different background characteristics. Here the simulated groups happen to have matching starting-weight and housing distributions. See https://marginaleffects.com/man/r/datagrid.html
Comparing three diets
avg_comparisons (variables = list (diet = "pairwise" ),type = "response"
Contrast Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
A - control 1.734 0.479 3.62 < 0.001 11.7 0.796 2.673
B - A -0.191 0.479 -0.40 0.68923 0.5 -1.130 0.747
B - control 1.543 0.479 3.22 0.00127 9.6 0.605 2.481
Term: diet
Type: response
"pairwise" requests all three pairs. "reference" compares each other diet with the first level, control . A vector such as c("control", "A") requests just one pair.
Point to each contrast label. The estimate is the first label in the printed subtraction minus the second label. Choose a planned pair when that answers the scientific question. When several tests form one family, decide how to address multiplicity. The pairwise call shown supplies separate p-values without a multiplicity adjustment. A Holm adjustment is one option for a family of planned tests. The displayed intervals are individual intervals. See https://marginaleffects.com/man/r/comparisons.html
A test against zero
avg_comparisons (variables = list (diet = diet_A),type = "response" ,hypothesis = 0
Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
1.73 0.479 3.62 <0.001 11.7 0.796 2.67
Term: diet
Type: response
Comparison: A - control
The null says A − control is 0 kg . The Wald statistic is the distance from zero in standard errors (1.73 / 0.479 ≈ 3.62). The two-sided p-value uses the normal distribution.
About 5 minutes for testing and interpretation. Zero is already the default here, so this repeats the earlier result with the null made explicit. The standard error measures the sampling uncertainty in the fitted contrast. The package uses the model covariance and the delta method, which approximates how coefficient uncertainty carries through to the quantity of interest. A two-sided p-value is the probability, if the null and test assumptions hold, of a statistic at least as far from zero as the one observed. Read the estimate and the 95% confidence interval before the p-value. The S column expresses the same p-value as minus log base two of p. The estimate, confidence limits and units matter for deciding whether the difference is useful. See https://marginaleffects.com/man/r/hypotheses.html
A test against 1 kg
Suppose the researchers use 1 kg as a practical benchmark for A − control.
avg_comparisons (variables = list (diet = diet_A),type = "response" ,hypothesis = 1
Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
1.73 0.479 1.53 0.125 3.0 0.796 2.67
Term: diet
Type: response
Comparison: A - control
The estimate and confidence interval still describe A − control. The test statistic now measures the distance from 1 kg in standard errors.
Compare the estimate, interval, statistic and p-value with the previous slide. The two-sided test asks whether the difference is above or below the benchmark. If the question specifically concerns exceeding a minimum, define that question and the direction in advance. For a ratio, equality between diets corresponds to a null of 1 instead of 0. Check the scale used for the test, especially if a transformation is applied for display. See https://marginaleffects.com/man/r/hypotheses.html
One effect per housing type
avg_comparisons (variables = list (diet = diet_A),by = "housing" ,type = "response"
housing Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
indoor 1.06 0.535 1.99 0.047 4.4 0.0143 2.11
outdoor 4.42 1.071 4.13 <0.001 14.7 2.3207 6.52
Term: diet
Type: response
Comparison: A - control
by gives one average difference per housing group. The first row is indoor and the second is outdoor . Each p-value tests that group’s difference against zero.
Connect these estimates to the dots on the comparison plot. The two separate p-values concern each contrast on its own. Use a direct test to estimate how much the diet effect changes across housing groups.
A test of how effects differ
avg_comparisons (variables = list (diet = diet_A),by = "housing" ,type = "response" ,hypothesis = "b2 - b1 = 0"
Hypothesis Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
b2-b1=0 3.36 1.2 2.8 0.00505 7.6 1.01 5.7
Type: response
The test subtracts the indoor diet effect from the outdoor diet effect. The standard error uses the covariance between the two fitted contrasts.
Read b1 and b2 from the output on the preceding slide before writing the equation. The estimate is 3.36 kg. Under the model the interval runs from about 1.01 to 5.70 kg. The standard error accounts for the fact that both estimates come from the same model. See https://marginaleffects.com/man/r/hypotheses.html
A joint test for the diet factor
<- avg_comparisons (variables = list (diet = "reference" ),type = "response" hypotheses (vs_control, joint = TRUE )
Joint hypothesis test:
A - control = 0
B - control = 0
F Pr(>|F|) Df 1 Df 2
7.89 <0.001 2 142
Type: response
One F test asks whether both average differences from control are zero.
This tests equality of all three average predictions in the target population of observed pens. It uses two independent restrictions because the third pairwise difference follows from the first two. A small joint p-value gives evidence that at least one average difference departs from zero. To identify individual differences, report the contrasts. For a family of separate pairwise decisions, Holm is one option for adjusting p-values. See https://marginaleffects.com/chapters/hypothesis.html
Reporting a diet comparison
Averaging paired predictions equally over the observed pens, we estimated that diet A increased four-week weight gain by 1.73 kg compared with control (95% CI 0.80 to 2.67 kg, two-sided p < 0.001).
For a mixed model, say how you handled the group effects. With an lme4 fit, re.form = NA sets group deviations to zero. In a nonlinear model this can differ from averaging predictions over the distribution of group effects.
Close by asking students to name the contrast, units, rows, weights and null before running a call in their own analysis. In repeated-measures data, equal row weight gives more influence to animals with more rows. Choose the grid and weights to reflect the intended population. Use the interpretation and reporting work in Unit 12 to communicate the resulting estimate.
Package references
The model and data in this deck are entirely simulated. The package references give details on other scales, model classes and ways to specify hypotheses.
