Modularity

Following along with the paper where \(b = (X^TX)^{-1}X^TY\) (beta values) and \(H = X(X^TX)^{-1}X^T\) (projection matrix), so \(\hat{y} = Hy\)

Lets take a small dataset where it’s easy to see what could happen. In this scenario we will look at how “hp” and “disp” are related to “mpg”. Think of “hp” and “disp” as something like movement and spiking where we want to remove those artifacts from our data. These variables are correlated as I’d imagine many artifacts are that you are trying to remove.

data("mtcars")
cor(mtcars$hp,mtcars$disp)

[1] 0.7909486

y<-mtcars$mpg

let’s get rid of the effect of “hp” on “mpg” first

x1<-cbind(rep(1,nrow(mtcars)),mtcars$hp)
b<-solve(t(x1)%*%x1)%*%t(x1)%*%y
h<-x1%*%solve(t(x1)%*%x1)%*%t(x1)
y.hat<-h%*%y
resids<-mtcars$mpg - y.hat

The residuals (in the context of fmri, our signal of interest with motion removed) should be uncorrelated with our nuisance variable “hp” and that is indeed the case.

round(cor(resids, mtcars$hp))

     [,1]
[1,]    0

OK, let’s move on to the next step of preprocessing. Let’s take our signal we care about (resids) and remove the effect of “disp”, this could be something like spiking.

x2<-cbind(rep(1,nrow(mtcars)),mtcars$disp)
b.2<-solve(t(x2)%*%x2)%*%t(x2)%*%resids
h.2<-x2%*%solve(t(x2)%*%x2)%*%t(x2)
y.hat.2<-h.2%*%resids
resids.2<-resids - y.hat.2

In theory, “resids.2” would be our data with the effects of motion (hp) and spiking (disp) removed, so let’s check.

cor(resids.2,mtcars$disp)

             [,1]
[1,] 3.300913e-16

Disp and resids.2 are uncorrelated which is good. What about resids.2 and hp?. Nope, because hp and disp are correlated (and as a result their projection matrices) we have actually reintroduced some of the covariance by projecting this data back into a space that shares some overlap with the hp (i.e, the motion regressor)

cor(resids.2,mtcars$hp)

          [,1]
[1,] 0.3155504

However, if you get rid of modularization and do everything in one step, you’re good to go.

x3<-cbind(rep(1,nrow(mtcars)),mtcars$hp,mtcars$disp)
b.3<-solve(t(x3)%*%x3)%*%t(x3)%*%y
h.3<-x3%*%solve(t(x3)%*%x3)%*%t(x3)
y.hat.3<-h.3%*%y
resids.3<-y - y.hat.3

Correlation between resids and nuisance regressors is 0.

round(cor(data.frame(resids.3,mtcars$hp,mtcars$disp)),2)

            resids.3 mtcars.hp mtcars.disp
resids.3           1      0.00        0.00
mtcars.hp          0      1.00        0.79
mtcars.disp        0      0.79        1.00