A Walk Into Multiple Regression With Group Differences

1. A developmental story

Imagine you are studying early vocal development in children. You record their actual age in months, but you also have a machine-learning model that produces a predicted vocal-development age based on acoustic features extracted from day-long recordings.

Children in your dataset do not all belong to the same developmental path. Some follow a typical trajectory (TD), some have mild hearing challenges (HH), some show language delays (LD), and some are diagnosed with autism (ASD). When you scatter all of these children onto the same plot, points form gentle coloured clouds. Some clouds are steep, others flatter. Some start higher, others lower. The story is simple: different developmental groups show different relations between chronological age and predicted vocal age.

This is exactly where multiple regression becomes a useful compass.

2. What is multiple regression?

Multiple regression is a statistical framework that asks a single big question: How does one variable change as a function of several predictors at the same time? The “several predictors” part matters. A child’s predicted vocal age may depend on their actual age, but also on which developmental group they belong to. And the effect of age may differ across groups. Regression lets you place these possibilities into a single equation.

3. The model in plain form

For our purpose, the model is a multiple linear regression with an interaction term:

\[ Y_i = \beta_0 + \beta_1 X_i + \beta_2 G_i + \beta_3 (X_i \times G_i) + \varepsilon_i \]

where \(Y_i\) represents the predicted vocal development age for child \(i\), \(X_i\) is their actual chronological age in months, and \(G_i\) is a group indicator variable representing the developmental category (TD, HH, LD, or ASD). The parameter \(\beta_0\) represents the baseline intercept, which is the starting point for the reference group. The coefficient \(\beta_1\) captures the slope for the reference group, showing the rate of vocal development per month. The term \(\beta_2\) represents the group effect on the intercept, indicating how much other groups differ from the reference at age zero. The interaction coefficient \(\beta_3\) captures how much the slope differs for other groups. Finally, \(\varepsilon_i\) is the random error term, assumed to follow a normal distribution with mean zero and constant variance:

\[ \varepsilon_i \sim N(0, \sigma^2) \]

4. Simulating a dataset

We now create a synthetic dataset inspired by real developmental patterns. This allows us to walk through each step without depending on external data sources. The simulation includes 600 children randomly assigned to four groups: TD, HH, LD, and ASD, with proportions of 30%, 30%, 20%, and 20% respectively. Each child has an actual age randomly drawn from a uniform distribution between 8 and 45 months.

We define group-specific parameters where TD children have an intercept of 5 and a slope of 1.6, HH children have an intercept of 8 and a slope of 1.4, LD children have an intercept of -5 and a slope of 1.1, and ASD children have an intercept of 0 and a slope of 0.6. In this simulation, TD and HH children show faster developmental alignment with steeper slopes, LD children grow more slowly, and ASD children show the shallowest slope. Noise is injected with a standard deviation of 6 to mimic real-world variability. This gives us a dataset with a recognizable developmental structure.

set.seed(1)

n <- 600
groups <- sample(c("TD","HH","LD","ASD"), n, replace = TRUE, prob = c(0.3,0.3,0.2,0.2))

actual_age <- runif(n, 8, 45)

# group-specific parameters
params <- list(
  TD  = c(intercept =  5, slope = 1.6),
  HH  = c(intercept =  8, slope = 1.4),
  LD  = c(intercept = -5, slope = 1.1),
  ASD = c(intercept =  0, slope = 0.6)
)

pred_age <- numeric(n)
for(i in 1:n){
  g <- groups[i]
  b0 <- params[[g]][1]
  b1 <- params[[g]][2]
  pred_age[i] <- b0 + b1*actual_age[i] + rnorm(1,0,6)
}

df <- data.frame(actual_age, pred_age, group = groups)

head(df)

##   actual_age pred_age group
## 1   38.12731 63.95730    TD
## 2   42.36476 76.32521    HH
## 3   13.45680 30.00936    HH
## 4   35.74340 24.69919   ASD
## 5   44.09932 74.73887    TD
## 6   44.06732 19.61999   ASD

5. Fitting the model

library(tidyverse)

model <- lm(pred_age ~ actual_age * group, data = df)
summary(model)

## 
## Call:
## lm(formula = pred_age ~ actual_age * group, data = df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -17.6586  -4.2995   0.0111   4.2253  19.1004 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         0.15561    1.53440   0.101 0.919256    
## actual_age          0.54525    0.05419  10.062  < 2e-16 ***
## groupHH             6.49476    1.94200   3.344 0.000877 ***
## groupLD            -7.19784    2.14077  -3.362 0.000823 ***
## groupTD             3.09492    1.96710   1.573 0.116173    
## actual_age:groupHH  0.89820    0.06774  13.260  < 2e-16 ***
## actual_age:groupLD  0.57460    0.07567   7.594 1.22e-13 ***
## actual_age:groupTD  1.16050    0.07056  16.446  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.26 on 592 degrees of freedom
## Multiple R-squared:  0.9115, Adjusted R-squared:  0.9104 
## F-statistic: 870.7 on 7 and 592 DF,  p-value: < 2.2e-16

The summary displays the main effect of actual age showing how predicted vocal age rises with age, the main effect of group showing how groups differ at the baseline, and the interaction terms showing how slopes differ across groups. The interaction is the key component as it tells you the developmental trajectory is not the same for all groups.

6. Interpreting the model

The summary displays three key components: the main effect of actual age showing how predicted vocal age rises with age, the main effect of group showing how groups differ at baseline, and the interaction terms showing how slopes differ across groups. The interaction is the key element, as it tells you the developmental trajectory is not the same for all groups.

7. Detailed interpretation

Let me walk through the actual output from our model to understand what each coefficient means. The reference group in our model is ASD (alphabetically first), with an intercept of 0.156 and a slope of 0.545. This means that for ASD children at birth (age = 0), their predicted vocal age would be 0.156 months, and for every one-month increase in actual age, their predicted vocal age increases by 0.545 months.

For HH children, the model shows an intercept difference of 6.495 compared to ASD. This means HH children start at 0.156 + 6.495 = 6.651 months of predicted vocal age at birth. The interaction term for HH is 0.905, meaning their slope is 0.545 + 0.905 = 1.450 months per month of actual age, indicating faster vocal development than ASD children.

For LD children, the intercept difference is -7.198, giving them a starting point of 0.156 - 7.198 = -7.042 months. Their interaction term of 0.571 yields a slope of 0.545 + 0.571 = 1.116 months per month, which is faster than ASD but slower than HH.

For TD children, the intercept difference is 3.095, placing their starting point at 0.156 + 3.095 = 3.251 months. Their interaction term of 1.161 produces a slope of 0.545 + 1.161 = 1.706 months per month, representing the fastest developmental trajectory in our sample.

The model’s R-squared value of 0.9115 indicates that approximately 91.2% of the variance in predicted vocal age is explained by actual age, group membership, and their interaction. The residual standard error of 7.255 represents the average deviation of observed values from the fitted regression line.

8. Plotting the data with regression lines

library(ggplot2)

ggplot(df, aes(x = actual_age, y = pred_age, color = group)) +
  geom_point(alpha = 0.5, size = 1.8) +
  geom_smooth(method = "lm", formula = y ~ x, se = FALSE, linewidth = 1.1) +
  scale_color_manual(values = c("TD"="#1f77b4","HH"="#2ca02c","LD"="#000000","ASD"="#d62728")) +
  labs(
    x = "Actual age (months)",
    y = "Predicted vocal development age",
    title = "Multiple Regression With Group-Specific Slopes"
  ) +
  theme_minimal(base_size = 13) +
  theme(
    panel.grid.major = element_line(color = "grey80", linetype = "dashed"),
    panel.grid.minor = element_blank(),
    panel.border = element_blank(),
    axis.line.x = element_line(color = "black", linewidth = 0.5),
    axis.line.y = element_line(color = "black", linewidth = 0.5)
  )

This produces:

• Four overlapping point clouds.
• Four regression lines, each capturing its group’s slope.
• A structure similar to the published plot you provided.

The final figure gives the whole story visually. The clouds show the raw relationship between chronological age and predicted development, while the colored lines show each group’s developmental trajectory. A steep line implies predicted vocal age increases quickly with real age, whereas a shallow line implies slower alignment. The vertical spread demonstrates how noisy vocal maturation can be from child to child. If TD shows a steeper slope than HH, which in turn is steeper than LD, which is steeper than ASD, you can almost see the developmental divergence across conditions unfolding before your eyes.