# Help laying out and interpreting multiple regression

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?

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:

• Hello Shawn, grazie for your thorough answer. Forgive me for possibly misunderstanding, but I do not see an answer to the core of my question. I will try to rephrase my question and remove unnecessary details. Suppose I set up a multiple regression with one continuous variable, one (2-level) categorical variable, and their interaction. Commented Sep 27, 2023 at 3:15
• This will result in four parameter estimates and two different lines with different slopes and intercepts, one line for each value of the categorical variable. Now, could I deduce from the p-values in the model output that I. the slope for the line corresponding to the baseline value of the categorical variable is not statistically significant, but II. the slope for the line corresponding to the other value of the categorical is statistically significant? Commented Sep 27, 2023 at 3:15
• I think not, because the slope estimates for each line are intertwined. But I may be wrong. Thank you for your help. Commented Sep 27, 2023 at 3:15
• (I may well ask my revised version as a separate question, just because I worry that my attempt above is too confusing.) Commented Sep 27, 2023 at 3:21
• You can examine the statistical significance of the slope beta_1 to see if income is significantly related to reading in the group (school type) that is coded as zero. To find out whether income is significantly related to reading in the other school type (i.e., the group that is coded as 1), you could simply flip the dummy coding and rerun the analysis. Then, beta_1 and its p value will refer to the relationship in the second group/second school type. Commented Sep 27, 2023 at 10:18