Is this a valid method and does it have a common name? (covariance analysis) We measured the same persons four times (T1a,T1b,T2a,T2b). Our theory is that the differences between T2a and T2b are caused by the influence of T1a on T2a and of T1b on T2b. Our hypothesis is that if we remove the influence of T1a/b on T2a/b, the difference between T2a/b will be zero. This is more or less like an analysis of covariance for repeated measures but I couldn't find anywhere what the GLM for this case looks like (or how it is done in R). 
Hence, we calculated the regression of T2a~T1a and T2b~T1b. Then we calculated 
T2a' = T2a - beta1*T1a   and   T2b' = T2b - beta1*T1b and then a paired t-test for comparing T2a' and T2b'.
Is our method legit and if so does it have a name? 
Is it comparable to an analysis of covariance for repeated measures? 
And if not, what do instead (and how to do it)?
 A: I don't think your analysis is valid. The problem is that when you calculate T2a' and T2b' there is uncertainty in the beta1s and you pretend that those estimates  are fixed and correct. Hence your standard error of the difference will be two low, and your estimate too low. 
You can use a multilevel model or a structural equation model. 
For a structural equation model, set up your data as wide, so a person is a row. 
For the multilevel model, set it up so each person has a row for their t2a score and a row for their t2b score.
Here's some code to load the libraries, and make up some data:
library(lavaan)
library(lme4)
set.seed(12345)
d  <- data.frame(id = 1:100,
                 t1a = rnorm(100),
                 t1b = rnorm(100),
                 t2a = rnorm(100),
                 t2b = rnorm(100))

Here's the structural equation model:
model1 <- "
  # regressions
  t2a ~ t1a + t1b
  t2b ~ t1a + t1b
  # intercepts
  t2a ~ a * 1
  t2b ~ b * 1

  # compare intercepts
  d := a - b"

fit1 <- sem(model1, d)
summary(fit1)

The important bit of output is the parameter d:
Defined Parameters:
                   Estimate  Std.Err  z-value  P(>|z|)
    d                -0.217    0.153   -1.421    0.155

Now I can reshape the data:
a <- setNames(d[c("t1a", "t1b", "t2a")],
              c("t1a", "t1b", "t2"))
b <- setNames(d[c("t1a", "t1b", "t2b")],
              c("t1a", "t1b", "t2"))

a$group <- "a"
b$group <- "b"

a$id <- 1:100
b$id <- 1:100

dLong <- rbind(a, b)

And fit a mixed model:
summary(lmer(t2 ~ t1a * group + t1b * group + (1|id), data = dLong))

The important parameter is groupb
Fixed effects:
            Estimate Std. Error t value
groupb       0.21719    0.13731   1.582

The two parameters are the same, and the SEs differ by a little (as the sample size increases, the difference drops). 
