# Method of alternating projections for linear fixed effects models [closed]

The standard fixed effects model (in econometrics, mostly) is $$y = \mathbf{X\beta} + \mathbf{D\alpha} + \epsilon$$ where $\mathbf{D}$ is a set of factors, potentially with thousands of levels. This is typical in longitudinal data -- $\mathbf{D}$ would be a matrix of dummies representing e.g. individuals.

When the dimensionality of $\mathbf{D}$ is high, and when there are more than one matrix of fixed effects, one way to project-out the fixed effects is the method of alternating projections, which follows the following algorithm:

This is implemented in the demeanlist function in the R package lfe. (All of this stuff is by Simen Gaure).

In the vignette for lfe, the author provides the following code snippet to illustrate how the method works:

demean <- function(x, flist) {
cx <- x
oldx <- x - 1
while(sqrt(sum((cx - oldx) ^ 2)) >= 1e-10) {
oldx <- cx
for(i in 1:length(flist)){
cx <- cx - ave(cx, flist[[i]])
}
}
return(cx)
}


To make sure that it works, I'm comparing the output against lfe's compiled function, and against projection via FWL, which is known to work but can be memory-prohibitive with large datasets:

demean.fwl <- function(x, flist){
xx <- model.matrix(~.-1, data= as.data.frame(flist))
Mx <- diag(rep(1, nrow(x))) - xx %*% solve(crossprod(xx)) %*% t(xx)
return(Mx %*% x)
}


The problem is that it doesn't work:

#fake data
x <- rnorm(1000)
x2 <- rnorm(length(x))
id <- factor(sample(20,length(x),replace=TRUE))
firm <- factor(sample(13,length(x),replace=TRUE))

#projecting...
Xdm <- demeanlist(X, list(id, firm))
Xdm.fwl <- demean.fwl(X, list(id, firm))
Xdm.r <- demean(X, list(id, firm))

#comparing...
1> all.equal(Xdm, Xdm.fwl, check.attributes = F)
[1] TRUE
1> all.equal(Xdm.r, Xdm.fwl, check.attributes = F)
[1] "Mean relative difference: 0.1290636"


The function demean provided in the vignette doesn't work, though it seems to hew to the algorithm.

Can anyone see why?

I'm interested in this because I want to add this functionality to an R packages that I'm writing, and I want to do it via a compiled function in Rcpp. demeanlist is already in C, but I want to embed it within a bunch of other compiled code (I am just learning C++).

## closed as off-topic by kjetil b halvorsen, Peter Flom♦Jul 13 at 13:50

This question appears to be off-topic. The users who voted to close gave this specific reason:

• "This question appears to be off-topic because EITHER it is not about statistics, machine learning, data analysis, data mining, or data visualization, OR it focuses on programming, debugging, or performing routine operations within a statistical computing platform. If the latter, you could try the support links we maintain." – kjetil b halvorsen, Peter Flom
If this question can be reworded to fit the rules in the help center, please edit the question.

I slightly modified your code (to remove run-time error) and ran it.

demean <- function(x, flist) {
cx <- x
oldx <- x - 1
while(sqrt(sum((cx - oldx) ^ 2)) >= 1e-10) {
oldx <- cx
for(i in 1:length(flist)){
cx <- cx - ave(cx, flist[[i]])
}
}
return(cx)
}

demean.fwl <- function(x, flist){
xx <- model.matrix(~.-1, data= as.data.frame(flist))
Mx <- diag(nrow(xx)) - xx %*% solve(crossprod(xx)) %*% t(xx)
return(Mx %*% x)
}

for(i in 1:10){
#fake data
set.seed(i)
x <- rnorm(1000)
id <- factor(sample(20, length(x), replace=TRUE))
firm <- factor(sample(13, length(x), replace=TRUE))

#projecting...
Xdm <- demeanlist(x, list(id, firm))
Xdm.fwl <- demean.fwl(x, list(id, firm))
Xdm.r <- demean(x, list(id, firm))
Xdm.fwl <- c(Xdm.fwl)

t1 <- all.equal(Xdm, Xdm.fwl, check.attributes = F)
t2 <- all.equal(Xdm.r, Xdm.fwl, check.attributes = F)
print(t1)
print(t2)
}


prints 20 TRUEs. It may be the problem has resolved recently.