Help laying out and interpreting multiple regression

Question

Let's say I want to model children's reading scores as a function of family income, type of school (K-6 or K-8), and an interaction. Therefore my model is

\begin{align} Y = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \beta_3 X_1 X_2 + \epsilon \end{align}

where $X_1$ is for income and $X_2$ is a dummy for school type (baseline K-6).

Question: I want to understand what kind of p-values in my R output would tell me that I have the following scenario:

• For K-6 schools, the change in reading score per unit change of income is not statistically significant.
• For K-8 schools, the change in reading score per unit change of income is statistically significant.

My attempt:

I have in expectation

\begin{align} \mathbb E \left[ Y \mid X_2 = 0 \right] = f(X_1) = \beta_0 + \beta_1 X_1\\ \mathbb E \left[ Y \mid X_2 = 1 \right] = g(X_1) = (\beta_0 + \beta_2) + (\beta_1 + \beta_3) X_1 \end{align}

So, the expected/predicted reading score, as a function of income given a particular school type, is two different lines with different intercepts and slopes. Because the slope coefficients for each line are intermingled, I am not sure any p-values in the R output would help me answer my question. For example,

\begin{align} p_{\beta_0} < 0.05\\ p_{\beta_1} > 0.05\\ p_{\beta_2} < 0.05\\ p_{\beta_3} < 0.05 \end{align}

would not tell me that I have the scenario in my question, would it?

Answer 1

Statistical significance will not tell you the magnitude of the effect in any way, and will consequently be uninformative if you use it as a metric of your coefficient effects. Even the word "effect" has to be chosen carefully here, as this is behavioral data and so you need to be a bit careful about the wording with relation to causality.

As to actual responses to your question:

So, the expected/predicted reading score, as a function of income given a particular school type, is two different lines with different intercepts and slopes.

Your intercept, $\beta_0$, will be the conditional mean of the reference group (if you had a treatment group and a control group, one of these is selected as a reference group to compare against). Your following coefficient $\beta_1$ will be adjustments to this conditional mean. So if the K-6 school's conditional mean is $90$ and K-8 school's slope is $2$, then the conditional mean for K-8 school is equal to $92$ after controlling for all other predictors (I'm assuming you mean two distinct schools and not literally one school that has kindergarten kids, another with primary, and so on). A visualized example straight from Regression and Other Stories is shown below (though here I use a frequentist version rather than the Bayesian example they give for your purposes):

#### Fit ####
fit <- lm(kid_score ~ mom_hs + mom_iq, data=kidiq)
summary(fit)

#### Plot ####  
colors <- ifelse(kidiq$mom_hs==1, "black", "gray")
plot(kidiq$mom_iq,
     kidiq$kid_score,
  xlab="Mother IQ score",
  ylab="Child test score",
  col=colors,
  pch=20)
b_hat <- coef(fit)
abline(b_hat[1] + b_hat[2], b_hat[3] + b_hat[4], col="black")
abline(b_hat[1], b_hat[3], col="gray")

You can see that the association between mother IQ and student test scores are positive in general. However, the conditional mean (intercept) shifts by whether or not the mom graduated, and the magnitude of the relationship between each mother changes as well (the slope).

Consequently, if you only have one continuous predictor and one categorical factor as main effects, along with their interaction, you should have 4 coefficients and their respective $p$ values, which match the number of $p$ values you have listed already:

• $\beta_0$: Intercept (School 1)
• $\beta_1$: School 2
• $\beta_2$: $X$
• $\beta_3$: School 2 $\times$ $X$

For comparison to our earlier regression, this can be easily shown by running the summary(fit) code above:

Call:
lm(formula = kid_score ~ mom_hs + mom_iq + mom_hs:mom_iq, data = kidiq)

Residuals:
    Min      1Q  Median      3Q     Max 
-52.092 -11.332   2.066  11.663  43.880 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)   -11.4820    13.7580  -0.835 0.404422    
mom_hs         51.2682    15.3376   3.343 0.000902 ***
mom_iq          0.9689     0.1483   6.531 1.84e-10 ***
mom_hs:mom_iq  -0.4843     0.1622  -2.985 0.002994 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 17.97 on 430 degrees of freedom
Multiple R-squared:  0.2301,    Adjusted R-squared:  0.2247 
F-statistic: 42.84 on 3 and 430 DF,  p-value: < 2.2e-16

We can see everything noted in the plot reflects the point estimates shown here, however the regression isn't centered, so it is likely not super interpretable. Better to center the predictors in order to make it easier to read:

#### Center Variable ####
kidiq$iq_center <- kidiq$mom_iq - mean(kidiq$mom_iq)

#### Fit Again ####
fit2 <- lm(formula = kid_score ~ mom_hs + iq_center + mom_hs:iq_center,
           data = kidiq)

#### Summarize ####
summary(fit2)

We now have something easier to interpret:

Call:
lm(formula = kid_score ~ mom_hs + iq_center + mom_hs:iq_center, 
    data = kidiq)

Residuals:
    Min      1Q  Median      3Q     Max 
-52.092 -11.332   2.066  11.663  43.880 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)       85.4069     2.2182  38.502  < 2e-16 ***
mom_hs             2.8408     2.4267   1.171  0.24239    
iq_center          0.9689     0.1483   6.531 1.84e-10 ***
mom_hs:iq_center  -0.4843     0.1622  -2.985  0.00299 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 17.97 on 430 degrees of freedom
Multiple R-squared:  0.2301,    Adjusted R-squared:  0.2247 
F-statistic: 42.84 on 3 and 430 DF,  p-value: < 2.2e-16

The conditional mean test score is around $85$ for those who have mothers that did not graduate. We increase our mean by about $3$ if they have graduated. For the slope, we interpret as the average change in test score with each one unit increase in IQ. Our interaction estimate here adjusts the slope by the number given, shown by the different slopes for each line in the plot.

Notice that the key here in the interaction is the slope change. If we simply fit the data like so with just the main effects:

#### Fit Again ####
fit3 <- lm(formula = kid_score ~ mom_hs + mom_iq,
           data = kidiq)

#### Plot ####  
colors <- ifelse(kidiq$mom_hs==1, "black", "gray")
plot(kidiq$mom_iq,
     kidiq$kid_score,
     xlab="Mother IQ score",
     ylab="Child test score",
     col=colors,
     pch=20)
b_hat <- coef(fit3)
abline(b_hat[1] + b_hat[2], b_hat[3], col="black")
abline(b_hat[1], b_hat[3], col="gray")

We get our positive trends, but the magnitude does not change by group: