Why Plot a Factor-by-Numeric Interaction?
When a model includes an interaction between a categorical predictor and a numeric predictor, the relationship between the numeric predictor and the fitted response can differ across categories. A fitted-line plot makes those group-specific relationships visible. It also provides a direct check on the algebra: each plotted line should have the slope and category differences implied by the model coefficients.
The plot shows model-implied fitted values, not necessarily observed values. Points can be added to show the data, but the lines themselves should be calculated from the fitted model. This distinction matters: observed points vary around the fitted values, while the line represents the conditional mean specified by the model.
From Coefficients to Fitted Lines
Consider an ordinary least squares linear regression with a numeric response \(Y\), numeric predictor \(x\), and factor \(G\) with levels A and B. Suppose A is the treatment-coded reference level. The model for the conditional mean is
where \(I(G=B)\) equals 1 for level B and 0 otherwise. The model-implied fitted value is obtained by replacing the coefficients with their estimates. For level A, both indicator terms are zero, so its line has intercept \(\beta_0\) and slope \(\beta_1\). For level B, the indicator terms are active, so its intercept is \(\beta_0+\gamma_B\) and its slope is \(\beta_1+\delta_B\).
The fitted difference between B and A at a particular value \(x\) is not generally just the factor coefficient. It is
Thus, \(\widehat{\gamma}_B\) is the fitted difference at \(x=0\), and \(\widehat{\delta}_B\) is the change in the B-versus-A difference for each one-unit increase in \(x\). If \(x=0\) is not a useful comparison point, the difference at a relevant value of \(x\) is usually more meaningful.
A Short Derivation: Why the Lines Have These Slopes
The fitted line for each category follows by substituting its indicator value into the model. For A, \(I(G=B)=0\), so the model reduces to \(\beta_0+\beta_1x\). For B, \(I(G=B)=1\), so it reduces to \(\beta_0+\beta_1x+\gamma_B+\delta_Bx\). Collecting the constant terms and the terms multiplied by \(x\) gives \((\beta_0+\gamma_B)+(\beta_1+\delta_B)x\). Subtracting the A expression from the B expression cancels the shared reference terms, leaving \(\gamma_B+\delta_Bx\). This is why the interaction coefficient changes the slope, while the category coefficient contributes to the difference at \(x=0\).
Worked Example: Check Three Lines Against Their Coefficients
Worked Example: Check Three Lines Against Their Coefficients
A hypothetical pricing analysis uses numeric predictor \(x\), measured in units, and a three-level factor with reference level A. The fitted linear model is \(\widehat{Y}=50+0.4x-5I(B)+0.3xI(B)+2I(C)-0.2xI(C)\). The observed \(x\)-values range from 0 to 20. Find each fitted line, compare the fitted values at \(x=10\), and identify the steepest line.
Method. Substitute the indicator values for each factor level, then calculate the intercept, slope, and fitted value. This avoids treating an interaction coefficient as the group’s entire slope.
For A, both indicators are zero:
For B, the B indicator is 1 and the C indicator is 0. Its intercept is \(50-5=45\), and its slope is \(0.4+0.3=0.7\), giving \(\widehat{Y}_B(x)=45+0.7x\). For C, its intercept is \(50+2=52\), and its slope is \(0.4-0.2=0.2\), giving \(\widehat{Y}_C(x)=52+0.2x\).
At \(x=10\), the fitted values are:
The B-versus-A difference is \(52-54=-2\); the C-versus-A difference is \(54-54=0\). These agree with the coefficient-based checks: for B, \(-5+0.3(10)=-2\); for C, \(2-0.2(10)=0\). B is steepest, with slope \(0.7\). A and C have the same fitted value at \(x=10\), but their slopes differ.
As a further plot check, at \(x=0\) the fitted values are 50, 45, and 52 for A, B, and C; at \(x=20\) they are 58, 59, and 56. Each pair of endpoint values is consistent with its line and slope. The values are checked by direct substitution into the three fitted-line equations.
Plotting Fitted Lines in R
A reliable plotting workflow is to fit the interaction model, create a prediction grid spanning the observed range of \(x\), calculate fitted values for each factor level, and draw one line per level. For example, with a data frame named dat and factor variable group, the R formula `lm(y ~ group * x, data = dat)` includes the main effects and their interaction. Ensure that group is a factor and that its reference level is the one intended for interpretation.
Create a grid with values of \(x\) from the observed minimum to maximum and one row for every combination of those \(x\)-values and the levels of group. Use `predict(fit, newdata = grid)` to add fitted values to the grid. Plot the fitted values against \(x\), using group for both line color and grouping. In ggplot2, the key layer is `geom_line(aes(color = group, group = group))`. Because the values being connected come from the fitted model, the displayed lines correspond to the model’s group-specific fitted relationships.
A grid covering the observed range is generally more useful than a range far beyond the data. A line can always be extended algebraically, but its fitted values outside the range of observed predictors may be extrapolations and should not be presented as well-supported comparisons. If observed data points are shown as well, distinguish them visually from the fitted lines.
Worked Example: Interpret a Centered Predictor Plot
Worked Example: Interpret a Centered Predictor Plot
A hypothetical staffing model uses workload \(x\) and factor levels A and B. Workload is centered at 40, so \(z=x-40\). The fitted model is \(\widehat{Y}=120+1.5z-8I(B)+0.5zI(B)\). Workload values from 36 to 44 are in the observed range. Find the group-specific lines in terms of \(z\), compare the fitted values at workloads 36, 40, and 44, and state what the coefficient \(-8\) represents.
Method. First use \(z=x-40\) to identify the centered predictor value at each workload. Then substitute the factor indicator and calculate each group’s intercept and slope in \(z\).
For A, \(I(B)=0\), so \(\widehat{Y}_A=120+1.5z\). For B, \(I(B)=1\), so the intercept is \(120-8=112\), and the slope is \(1.5+0.5=2.0\). Therefore, \(\widehat{Y}_B=112+2.0z\).
At workload 36, \(z=36-40=-4\). The fitted values are \(\widehat{Y}_A=120+1.5(-4)=120-6=114\) and \(\widehat{Y}_B=112+2(-4)=112-8=104\). At workload 40, \(z=0\), giving A = 120 and B = 112. At workload 44, \(z=44-40=4\), giving A \(=120+6=126\) and B \(=112+8=120\).
The B-versus-A differences are \(-10\), \(-8\), and \(-6\) at workloads 36, 40, and 44. The formula \(-8+0.5z\) gives the same results: \(-8+0.5(-4)=-10\), \(-8+0=-8\), and \(-8+0.5(4)=-6\). The coefficient \(-8\) is the fitted difference at workload 40, where \(z=0\), not at workload zero. B’s slope is 0.5 units per workload unit greater than A’s.
A plot with workload on the horizontal axis should show two lines. At workload 40, the vertical separation is 8 fitted-response units, with B below A. Moving to the right, the separation narrows because B’s line has the greater slope. The calculated values confirm the direction and size of that change.
Worked Example: Use Plot Coordinates to Catch a Coefficient Error
Worked Example: Use Plot Coordinates to Catch a Coefficient Error
A hypothetical processing-time model has factor levels A, B, and C, with A as reference. The fitted equation is \(\widehat{Y}=30+2x+6I(B)-0.5xI(B)-4I(C)+xI(C)\), and the plot is intended to cover \(x=0\) to \(x=10\). Calculate each line and its fitted values at \(x=0,5,10\). A draft plot shows the B line with slope 0.5 and the C line with slope 1. Check whether it matches the equation.
Method. Add each level’s interaction coefficient to the reference slope. Then calculate line coordinates by substitution; comparing these coordinates with the draft makes a mistaken slope or intercept visible.
For A, the line is \(\widehat{Y}_A=30+2x\). For B, the intercept is \(30+6=36\), and the slope is \(2-0.5=1.5\), so \(\widehat{Y}_B=36+1.5x\). For C, the intercept is \(30-4=26\), and the slope is \(2+1=3\), so \(\widehat{Y}_C=26+3x\).
At \(x=0\), the fitted values are A \(=30\), B \(=36\), and C \(=26\). At \(x=5\), they are A \(=30+2(5)=40\), B \(=36+1.5(5)=36+7.5=43.5\), and C \(=26+3(5)=41\). At \(x=10\), they are A \(=50\), B \(=36+15=51\), and C \(=26+30=56\).
The draft does not match: the B slope must be 1.5, not 0.5, and the C slope must be 3, not 1. For another check, at \(x=5\) the B-versus-A difference is \(43.5-40=3.5\), which equals \(6-0.5(5)=3.5\). The C-versus-A difference is \(41-40=1\), which equals \(-4+1(5)=1\). Both differences agree with the interaction equation.
Common Plotting and Interpretation Traps
- Using an interaction coefficient as the whole group slope. The slope for a non-reference level is the numeric-predictor coefficient plus that level’s interaction coefficient.
- Comparing categories using only their factor coefficient. That coefficient is the difference at \(x=0\); at another value, include the interaction term.
- Reading an intercept at the wrong origin. If the predictor is centered, \(x=0\) in the model may correspond to a nonzero value of the original predictor.
- Plotting a line inconsistent with the fitted model. Check at least two coordinates per line and confirm its slope. A line with the right intercept but the wrong slope still misrepresents the interaction.
- Confusing fitted values with observed outcomes. The line represents model-implied conditional means; individual observations need not lie on it.
- Extrapolating without noting the data range. A plotted line beyond the observed predictor range is algebraically defined, but its behavior there may be poorly supported by data.
Check Your Understanding
For each problem, derive the relevant lines before checking the answer.
- A treatment-coded model has reference-level line \(\widehat{Y}_A=18+0.6x\), with B factor coefficient 4 and B interaction coefficient \(-0.2\). Find both fitted values at \(x=5\), and state the B slope. (Answer: A \(=21\), B \(=24\), B slope \(=0.4\).)
- A model has \(\widehat{Y}=70+1.2x-3I(B)+0.4xI(B)\). Find A and B fitted values at \(x=10\) and check the B-versus-A difference by both methods. (Answer: A \(=82\), B \(=83\); difference \(=1=-3+0.4(10)\).)
- A centered-predictor model is \(\widehat{Y}=9+0.5z+2I(B)-0.1zI(B)\), where \(z=x-12\). Find A and B fitted values at \(x=10\) and the category difference there. (Answer: A \(=8\), B \(=10.2\); B-minus-A \(=2.2\).)
- For a two-level model, the B-minus-A fitted difference is \(-6+0.3x\). At what \(x\)-value do the lines meet, and what does that equality mean? (Answer: \(x=20\); the fitted values for A and B are equal there.)
- A three-level model has lines \(\widehat{Y}_A=40+0.8x\), \(\widehat{Y}_B=45+0.4x\), and \(\widehat{Y}_C=38+1.4x\). Calculate all fitted values at \(x=5\), identify the steepest line, and compare B and C with A at that value. (Answer: A \(=44\), B \(=47\), C \(=45\); C is steepest; B-minus-A \(=3\), C-minus-A \(=1\).)