3
$\begingroup$

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:

enter image description here

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++).

$\endgroup$

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.

0
$\begingroup$

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.

$\endgroup$

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