What's the fundamental difference between these two regression models?

Question

Suppose I have a bivariate responses with significant correlation. I am trying to compare the two ways to model these outcomes. One way is to model the difference between the two outcomes: $$(y_{i2}-y_{i1}=\beta_0+X'\beta)$$ Another way is to use gls or gee to model them: $$(y_{ij}=\beta_0+\text{time}+X'\beta)$$

Here is a foo example:

#create foo data frame

require(mvtnorm)
require(reshape)
set.seed(123456)
sigma <- matrix(c(4,2,2,3), ncol=2)
y <- rmvnorm(n=500, mean=c(1,2), sigma=sigma)
cor(y)
x1<-rnorm(500)
x2<-rbinom(500,1,0.4)
df.wide<-data.frame(id=seq(1,500,1),y1=y[,1],y2=y[,2],x1,x2)
df.long<-reshape(df.wide,idvar="id",varying=list(2:3),v.names="y",direction="long")
df.long<-df.long[order(df.long$id),]
    df.wide$diff_y<-df.wide$y2-df.wide$y1


#regressions
fit1<-lm(diff_y~x1+x2,data=df.wide)
fit2<-lm(y~time+x1+x2,data=df.long)
fit3<-gls(y~time+x1+x2,data=df.long, correlation = corAR1(form = ~ 1 | time))

What's the fundamental difference between fit1 and fit2? And between fit2 and fit3, given they are so close on the $p$ values and estimates?

Answer 1

First, I will introduce yet a fourth model for the discussion in my answer:

fit1.5 <- lm(y_2 ~ x_1 + x_2 + y_1)

Part 0
The difference between fit1 and fit1.5 is best summarized as the difference between a constrained difference vs. an optimal difference.

I am going to use a simpler example to explain this than the one provided above. Let's start with fit1.5. A simpler version of the model would be $$y_2 = b_0 + b_1·x + b_2·y_1$$ Of course, when we obtain an OLS estimate, it will find the "optimal" choice for $b_2$. And, though it seems strange to write is as such, we could rewrite the formula as $$y_2 - b_2·y_1 = b_0 + b_1·x$$ We can think of this as the "optimal" difference between the two $y$ variables.

Now, if we decide to constraint $b_2=1$, then the formula/model becomes $$y_2 - y_1 = b_0 + b_1·x$$ which is just the (constrained) difference.

Note, in the above demonstration, if you let $x$ be a dichotomous variable, and $y_1$ be a pre-test and $y_2$ a post test score pairing, then the constrained difference model would just be the independent samples $t$-test for the gain in scores, whereas the optimal difference model would be the ANCOVA test with the pre-test scores being used as covariates.

Part 1
The model for fit2 can best be thought of in a similar fashion to the the difference approach used above. Though this is an oversimplification (as I am purposefully leaving out the error terms), the model could be presented as $$y = b_0 + b_1 · x + b_2 · t$$ where $t=0$ for the $y_1$ values and $t=1$ for the $y_2$ values. Here is the oversimplification...this let's us write $$\begin{align}y_1 & = b_0 + b_1 · x \\ y_2 & = b_0 + b_1 · x + b_2\end{align}$$ Written another way, $y_2 - y_1 = b_2$. Whereas model fit1.5 had $b_2$ as the value to make the optimal difference for the OLS analysis, here $b_2$ is essentially just the average difference between the $y$ values (after controlling for the other covariates).

Part 2
So what is the difference between models fit2 and fit3...actually, very little. The fit3 model does account for correlation in error terms, but this only changes the estimation process, and thus the differences between the two model outputs will be minimal (beyond the fact that the fit3 estimates the autoregressive factor).

Part 2.5
And I will include yet one more model in this discussion

fit4 <- lmer(y~time+x1+x2 + (1|id),data=df.long)

This mixed-effects model does a slightly different version of the autoregressive approach. If we were to include the time coefficient in the random effects, this would be comparable to calculating the difference between the $y$s for each subject. (But, this won't work...and the model won't run.)