I am running multiple regression models testing whether a change from T1 to T2 in my IV predicts a change in my DV (also from T1 to T2). I am running these analyses in R using the glm function.

I ran the same model again in a slightly different manner, predicting change in my DV directly rather than by controlling for my DV at T1. However, this yields slightly different results. Why is this? I assumed the two models would be doing exactly the same thing, but apparently not. Both models also contain both my IV and my DV at two measurement points, so this can't be due to more missings in of the two models. Here is my R-code:

### Version 1: 
M0 <- glm(IV.t2 ~ IV.t1,
          data = data,
          na.action = na.omit)

Res_IV <- resid(M0)

M1 <- glm(DV.t2 ~ DV.t1 + Res_IV,
          data = data,
          na.action = na.omit)

### Version 2:
M0.1 <- glm(IV.t2 ~ IV.t1,
          data = data,
          na.action = na.omit)

Res_IV <- resid(M0.1)

M0.2 <- glm(DV.t2 ~ DV.t1,
          data = data,
          na.action = na.omit)

Res_DV <- resid(M0.2)

M2 <- glm(Res_DV ~ Res_IV,
          data = coronaUK_T1_T2_cost,
          na.action = na.omit)

What might be the reason for these differences? Do the two analyses mathematically do something different?


1 Answer 1


I'm not sure I follow your code, and how it relates to what you have in mind.

For starters, it's slightly odd that you are using glm() but aren't setting the family argument. By default, family=gaussian, in which case, glm() fits a standard linear model, like what you would get by using lm(). The differences are that glm() estimates the parameters by searching using a variant of the Newton-Raphson algorithm, instead of using ordinary least squares (OLS). In addition, it tests the parameters using Wald $z$-tests, instead of $t$-tests. In both cases, there are (slight) reasons to prefer the methods used by lm(). Moreover, there are methods that operate over the output (e.g., summary()) that will provide more useful information if lm() was used. Was this supposed to be a model for non-normal data (e.g., a binary response or a count)?

Second, neither of your models uses the difference between measures at time 1 and time 2, as far as I can tell. Instead, you seem to be using the part of the time 2 measure that cannot be explained by knowledge of the time 1 measure. That's not the same thing. Consider this simple simulation (coded in R):

set.seed(5758)  # this makes the example exactly reproducible

t1  =             round(rnorm(6, mean=50, sd=10), digits=1)
t2  = 10 + 2*t1 + round(rnorm(6, mean=0,  sd=1),  digits=1)
dif = t2-t1
m   = lm(t2~t1)
round(coef(summary(m)), digits=2)
#             Estimate Std. Error t value Pr(>|t|)
# (Intercept)    10.00       1.66    6.03        0
# t1              2.01       0.03   64.71        0
res = residuals(m)
cor(dif,res)  # [1] 0.06151258
data.frame(t1, t2, dif, res)
#     t1    t2  dif         res
# 1 62.3 134.9 72.6 -0.07585836
# 2 64.6 140.4 75.8  0.81043618
# 3 58.2 125.4 67.2 -1.35142689
# 4 52.6 115.9 63.3  0.38194292
# 5 42.6  95.8 53.2  0.34153188
# 6 33.9  77.9 44.0 -0.10662574
windows(width=7, height=4)
  layout(matrix(1:2, nrow=1))
  dotchart(t1, xlim=c(min(t1,t2),max(t1,t2)), cex=.8,
           main="Values\n(w/ differences implied)")
  axis(side=2, at=1:6, cex=.8)
  points(t2, 1:6, pch="+")
  plot(t1, t2, main="Scatterplot\n(w/ residuals implied)", cex.axis=.8)
  abline(coef(m), col="gray")
  plot(dif, res)

enter image description here

enter image description here

Setting that aside, we can ask if the answers have to / should be the same. The answer is no. Consider these simulations:


y1c = rnorm(10, mean=100, sd=15)
y2c = y1c + rnorm(10, mean=0, sd=1)

y1d = rnorm(10, mean=0, sd=1)
y2d = rnorm(10, mean=5, sd=1)

round(coef(summary(lm(y2c~y1c))), digits=2)
#             Estimate Std. Error t value Pr(>|t|)
# (Intercept)     5.00       4.81    1.04     0.33
# y1c             0.95       0.05   20.76     0.00
round(coef(summary(lm(y2c-y1c~1))), digits=2)
#             Estimate Std. Error t value Pr(>|t|)
# (Intercept)      0.1       0.35    0.27     0.79

round(coef(summary(lm(y2d~y1d))), digits=2)
#             Estimate Std. Error t value Pr(>|t|)
# (Intercept)     4.97       0.48   10.29     0.00
# y1d            -0.08       0.55   -0.15     0.89
round(coef(summary(lm(y2d-y1d~1))), digits=2)
#             Estimate Std. Error t value Pr(>|t|)
# (Intercept)     5.43       0.49   11.18        0

The short answer is that the two methods can give different answers because they correspond to different questions. In general, regression models treat $X$ and $Y$ variables asymmetrically. Specifically, $x$-values are taken to be fixed and known, whereas $y$-values are considered to be a mix of true signal and random noise. Using change scores as your response, means that you are assuming the differences are a combination of input from $X$ and some pure randomness. On the other hand, when you model time 2 values controlling for time 1 values, the time 1 values are treated as constants measured with perfect accuracy, and the time 2 values are partly stochastic error.

This leads to a third question, which should you use? You should be aware that, as @rolando2 points out, this is a very contentious topic in statistics. There have been various camps in favor of each that have been fussing at each other for decades. The threads he links to are worth reading to get some perspective. If pressed, I suppose I would say you should use the method that best corresponds to your research question.

  • $\begingroup$ Recent article (semanticscholar.org/paper/…) by Tennant et al. (2021) on this topic arguing in favor of Y1 as a regressor based on causal graphs. $\endgroup$
    – RobertF
    Commented Aug 17, 2022 at 17:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.