Robust Inference for Interaction Models

Why Robust Inference Needs the Whole Interaction

An interaction model allows the relationship between a quantitative predictor and an outcome to differ across groups. That flexibility is useful, but it means many practical questions are about combinations of coefficients, not single coefficients. For instance, the slope for one group may be the sum of a reference-group slope and an interaction coefficient. Its uncertainty therefore depends on both coefficient variances and their covariance.

When the error variance differs across observations, the usual least-squares covariance formula can give unreliable standard errors. A heteroskedasticity-consistent covariance estimate provides a robust alternative under appropriate conditions. This tutorial uses it to make inference about slopes and group differences in an interaction model. The regression coefficients themselves are still estimated by ordinary least squares; what changes is the estimated uncertainty around them.

The Interaction Model and Its Assumptions

Definition: Let $Y_i$ be a quantitative outcome, $x_i$ a quantitative predictor, and $G_i$ a group indicator coded 0 for a reference group and 1 for a comparison group. The two-group interaction model is $$ Y_i=\beta_0+\beta_1x_i+\beta_2G_i+\beta_3x_iG_i+\varepsilon_i. $$ Here $\beta_1$ is the reference-group slope, $\beta_1+\beta_3$ is the comparison-group slope, and $\beta_2+\beta_3x$ is the comparison-minus-reference difference in mean outcome at predictor value $x$.

This is a model for the conditional mean, with $E(\varepsilon_i\mid x_i,G_i)=0$. It assumes the mean is linear in the included terms, observations are independent under the sampling design, and the model matrix has full column rank so the coefficients can be estimated. For heteroskedasticity-robust inference, the error variance need not be constant: $\operatorname{Var}(\varepsilon_i\mid x_i,G_i)$ may vary across observations. The robust covariance estimate is justified by the zero conditional mean, independence, suitable finite error moments, and a sufficiently informative sample; it does not make small-sample inference exact.

Write the coefficient estimate vector as $\hat{\boldsymbol\beta}=(\hat\beta_0,\hat\beta_1,\hat\beta_2,\hat\beta_3)^\mathsf{T}$, and let $\widehat V_{\mathrm{HC}}$ be a heteroskedasticity-consistent estimate of its covariance matrix. For example, HC3 is one commonly used version. The diagonal entries estimate coefficient variances; off-diagonal entries estimate covariances. In particular, retaining the off-diagonal entries is essential when calculating uncertainty for a sum or difference of coefficients.

Linear Contrasts: The Calculation Behind the Inference

A linear contrast is a specified weighted combination of coefficients, written $L=\mathbf{c}^\mathsf{T}\boldsymbol\beta$. Its estimate is $\hat L=\mathbf{c}^\mathsf{T}\hat{\boldsymbol\beta}$. If $\widehat V_{\mathrm{HC}}$ estimates the covariance matrix of $\hat{\boldsymbol\beta}$, then the estimated variance and standard error of the contrast are

\widehat{\operatorname{Var}}_{\mathrm{HC}}(\hat L) =\mathbf{c}^\mathsf{T}\widehat V_{\mathrm{HC}}\mathbf{c}, \qquad SE_{\mathrm{HC}}(\hat L) =\sqrt{\mathbf{c}^\mathsf{T}\widehat V_{\mathrm{HC}}\mathbf{c}}.

This result follows from the covariance rule for a linear combination. If $A=\hat{\boldsymbol\beta}-\boldsymbol\beta$, then $\hat L-L=\mathbf{c}^\mathsf{T}A$. The variance of a weighted sum is

\operatorname{Var}(\mathbf{c}^\mathsf{T}A) =\sum_j\sum_k c_jc_k\operatorname{Cov}(A_j,A_k) =\mathbf{c}^\mathsf{T}\operatorname{Cov}(\hat{\boldsymbol\beta})\mathbf{c}.

Replacing the unknown covariance matrix by $\widehat V_{\mathrm{HC}}$ gives the estimated variance above. For a slope in the comparison group, $\mathbf{c}=(0,1,0,1)^\mathsf{T}$, so the variance is $\widehat{\operatorname{Var}}(\hat\beta_1)+\widehat{\operatorname{Var}}(\hat\beta_3)+2\widehat{\operatorname{Cov}}(\hat\beta_1,\hat\beta_3)$. Omitting the covariance term is generally incorrect.

For a large-sample two-sided test of $H_0:L=L_0$, the robust Wald statistic is $z=(\hat L-L_0)/SE_{\mathrm{HC}}(\hat L)$, compared with the standard normal distribution as an approximation. A corresponding approximate 95% confidence interval is $\hat L\pm1.96SE_{\mathrm{HC}}(\hat L)$. Small-sample performance can be limited, so the approximation and the analysis design should be considered rather than treating robust inference as an exact procedure.

Worked Example: The Comparison-Group Slope

In a hypothetical study, 72 independent clinic visits provide an outcome measured in symptom-score points, a dosage in milligrams, and a group indicator for daytime ($G=0$) versus evening ($G=1$) visits. The fitted interaction model has $\hat\beta_1=0.80$ points per milligram and $\hat\beta_3=-0.35$ points per milligram. An HC3 covariance estimate gives $\widehat{\operatorname{Var}}(\hat\beta_1)=0.0625$, $\widehat{\operatorname{Var}}(\hat\beta_3)=0.0400$, and $\widehat{\operatorname{Cov}}(\hat\beta_1,\hat\beta_3)=-0.0150$. Suppose the observations are independent clinic visits, the residual-versus-fitted plot shows no clear curvature in the mean, and the residual spread increases with fitted score. The spread pattern motivates robust uncertainty estimates; the assumed mean form and independence are plausible for this hypothetical setup.

The daytime slope is $\hat\beta_1=0.80$ points per milligram. The evening slope is the contrast $\hat\beta_1+\hat\beta_3$, so its estimate is $0.80+(-0.35)=0.45$ points per milligram. Its robust variance includes the covariance:

\widehat{\operatorname{Var}}_{\mathrm{HC}}(\hat\beta_1+\hat\beta_3) =0.0625+0.0400+2(-0.0150)=0.0725.

Thus $SE_{\mathrm{HC}}=\sqrt{0.0725}=0.2693$ points per milligram. The approximate 95% interval is $0.45\pm1.96(0.2693)=0.45\pm0.5278$, or $(-0.0778,0.978)$ points per milligram. It is reasonable, at this approximation level, to say that the data are compatible with a range of evening slopes that includes zero. The interval does not establish that the true slope is zero.

Contrasts for Group Differences at a Chosen Predictor Value

A different question is whether the groups have different expected outcomes at a particular predictor value $x_0$. From the model, the comparison-minus-reference difference at $x_0$ is $\beta_2+x_0\beta_3$. This contrast depends on the interaction coefficient, so its standard error changes with $x_0$. It also depends on the covariance between $\hat\beta_2$ and $\hat\beta_3$. The value $x_0$ and its units should be stated when reporting the result.

Worked Example: A Group Difference at Three Hours

Suppose a hypothetical sample contains 96 independent participants, with response measured in test-score points, predictor $x$ measured in hours of practice, and $G=1$ indicating a comparison program versus $G=0$ for the reference program. Assume the model's linear mean is plausible over the observed practice range and records are independent. The residual plot suggests changing spread but no obvious mean curvature. The fitted coefficients are $\hat\beta_2=1.60$ points and $\hat\beta_3=0.25$ points per hour. HC3 gives $\widehat{\operatorname{Var}}(\hat\beta_2)=0.25$, $\widehat{\operatorname{Var}}(\hat\beta_3)=0.04$, and $\widehat{\operatorname{Cov}}(\hat\beta_2,\hat\beta_3)=-0.02$.

At $x_0=3$ hours, the estimated mean difference, comparison minus reference, is $1.60+3(0.25)=2.35$ score points. The contrast uses $\mathbf{c}=(0,0,1,3)^\mathsf{T}$. Its robust variance and standard error are

\widehat{\operatorname{Var}}_{\mathrm{HC}}(\hat\beta_2+3\hat\beta_3) =0.25+3^2(0.04)+2(3)(-0.02) =0.49, \qquad SE_{\mathrm{HC}}=\sqrt{0.49}=0.700.

The approximate 95% interval is $2.35\pm1.96(0.700)$, or $(0.978,3.72)$ score points. Under the model, this estimates the difference in mean test score between programs among participants with three hours of practice. It is not a single group difference that applies at every practice duration; the estimated difference changes with $x$.

Testing Whether the Slopes Differ

The interaction coefficient directly represents the difference in slopes: comparison-group slope minus reference-group slope equals $\beta_3$. Therefore, testing whether the population slopes are equal is testing $H_0:\beta_3=0$ against an alternative such as $H_a:\beta_3\ne0$. This particular contrast has standard error $\sqrt{\widehat{\operatorname{Var}}_{\mathrm{HC}}(\hat\beta_3)}$. A test of the interaction term is not the same as a test that either group's slope is zero.

Worked Example: Evidence of a Slope Difference

In a hypothetical dataset of 120 independent participants, $x$ is weekly exercise time in hours and $Y$ is resting heart rate in beats per minute. Group $G=1$ is a comparison group and $G=0$ the reference group. Suppose the model has a plausible linear mean over the observed range, no strong residual evidence of mean curvature, and unequal residual spread; observations are independent by design. The fitted interaction coefficient is $\hat\beta_3=-0.30$ beats per minute per hour, with HC3 standard error $0.12$ beats per minute per hour.

For $H_0:\beta_3=0$, the robust statistic is $z=(-0.30-0)/0.12=-2.50$. The two-sided standard-normal p-value is approximately $0.0124$. Because this is below a significance level of $\alpha=0.05$, the data provide evidence that the groups' linear slopes differ. The estimated comparison-group slope is 0.30 beats per minute per hour lower than the reference-group slope. This is evidence of an association in the fitted model, not evidence that group membership caused the difference.

The approximate 95% interval for the slope difference is $-0.30\pm1.96(0.12)=(-0.535,-0.065)$ beats per minute per hour. It contains only negative values, consistent with the two-sided test result. These are large-sample robust calculations; their reliability still depends on the model and sampling conditions, not just on the choice of standard error.

Common Pitfall: Treating the Coefficients Separately

A common mistake is to report the interaction coefficient as if it were the slope for the comparison group. It is not: $\beta_3$ is the difference between the two slopes, while the comparison-group slope is $\beta_1+\beta_3$. Another mistake is to calculate the standard error of that sum by adding the two coefficient standard errors. The standard error of a contrast comes from its variance, which includes covariance terms; standard errors themselves do not add in that way.

Robust covariance estimation addresses uncertainty under nonconstant error variance, subject to its assumptions. It does not repair an incorrect conditional-mean form, remove confounding, make dependent observations independent, or turn an observational association into a causal effect. Before calculating a contrast, state exactly which slope or group difference is of interest, write it as a coefficient combination, and use the matching entries of the robust covariance matrix.

Key takeaway: In an interaction model, state the target contrast first. Calculate its estimate as a coefficient combination and its robust standard error as $\sqrt{\mathbf{c}^\mathsf{T}\widehat V_{\mathrm{HC}}\mathbf{c}}$, including all relevant covariances. Interpret the result at the predictor value and in the units that define the contrast.

Check Your Understanding

Use the interaction model and robust-contrast rules to answer these questions.

For the model $Y_i=\beta_0+\beta_1x_i+\beta_2G_i+\beta_3x_iG_i+\varepsilon_i$, what is the slope when $G=0$, and what is the slope when $G=1$?
At $x=5$, write the comparison-minus-reference difference in mean outcome as a linear contrast.
For a slope contrast $\hat\beta_1+\hat\beta_3$, which three entries of the covariance matrix determine its estimated variance?
If the estimated interaction coefficient is negative, what does its sign say about the comparison-group slope relative to the reference-group slope?
Name one important problem that a heteroskedasticity-consistent covariance estimate does not fix.