R help - Feb 2009 - fixed effects regression

Hi everyone,

I am running this fixed effects model in R.
The thing is, my own code and the canned routine keep returning different
estimates for the fixed effects factors and intercept.
(The other variables' estimates and R^2 are the same!!! weird!!)

I attach my own code, just in case.
I am pretty sure this would not have any fundamental error. Well, hopefully. I
have been using this code for a while..)

And does anyone know how I can include another fixed effect into this code? :(
Any help will be deeply appreciated....





feols = function(y,X,id) {
    n = length(y);
    uniq_id = unique(id);
    n_id = length(uniq_id);
    dummy = matrix(rep(0,n*(n_id-1)),nrow=n);

    for (i in 1:(n_id-1))
        dummy[id==uniq_id[i],i] = 1; 

    X = cbind(rep(1,n),X,dummy);
    k = ncol(X);

    invXX = solve(t(X) %*% X);
    b = as.numeric(invXX %*% (t(X) %*% y));
    res = y - X%*%b;
    s2 = as.numeric(t(res) %*% res /(n-k));
    omega = as.numeric(s2) * invXX;
    se = sqrt(diag(omega));
    dev = y - mean(y);
    r2 = 1 - as.numeric(t(res)%*%res)/as.numeric(t(dev) %*% dev);
    t = b/se;
    list(b=b,se=se,t=t,r2=r2,s2=s2,omega=omega,res=res);
}



B


_________________________________________________________________
[[elided Hotmail spam]]

	[[alternative HTML version deleted]]

I think your problem is that you are using SAS-style contrasts rather than
treatment contrasts.  That is, your 'dummy' matrix omits a final column,
whereas you should be omitting the first column.

Here is a version of your function which I find easier to follow.  You might
like to work through it.

feols <- function(y, X, id) {

  dummy <- outer(id, unique(id), "==") + 0
  X <- cbind(1, X, dummy[, -1])
  error_df <- nrow(X) - ncol(X)

  invXX <- solve(crossprod(X))
  b <- as.vector(invXX %*% crossprod(X, y))
  res <- as.vector(y - X %*% b)
  s2 <- sum(res^2)/error_df
  omega <- s2 * invXX
  se <- sqrt(diag(omega))
  r2 <- 1 - sum(res^2) / sum((y -mean(y))^2)

  list(b = b, se = se, t_stat = b/se, r2 = r2, s2 = s2,
       omega = omega, res = res, error_df = error_df,
       sig = 2 * pt(-abs(t_stat), error_df))
}

Warning: untested code.

Bill Venables.


Bill Venables
http://www.cmis.csiro.au/bill.venables/ 


-----Original Message-----
From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org]
On Behalf Of parkbomee
Sent: Saturday, 7 February 2009 1:57 PM
To: r-help at r-project.org
Subject: [R] fixed effects regression


Hi everyone,

I am running this fixed effects model in R.
The thing is, my own code and the canned routine keep returning different
estimates for the fixed effects factors and intercept.
(The other variables' estimates and R^2 are the same!!! weird!!)

I attach my own code, just in case.
I am pretty sure this would not have any fundamental error. Well, hopefully. I
have been using this code for a while..)

And does anyone know how I can include another fixed effect into this code? :(
Any help will be deeply appreciated....





feols = function(y,X,id) {
    n = length(y);
    uniq_id = unique(id);
    n_id = length(uniq_id);
    dummy = matrix(rep(0,n*(n_id-1)),nrow=n);

    for (i in 1:(n_id-1))
        dummy[id==uniq_id[i],i] = 1; 

    X = cbind(rep(1,n),X,dummy);
    k = ncol(X);

    invXX = solve(t(X) %*% X);
    b = as.numeric(invXX %*% (t(X) %*% y));
    res = y - X%*%b;
    s2 = as.numeric(t(res) %*% res /(n-k));
    omega = as.numeric(s2) * invXX;
    se = sqrt(diag(omega));
    dev = y - mean(y);
    r2 = 1 - as.numeric(t(res)%*%res)/as.numeric(t(dev) %*% dev);
    t = b/se;
    list(b=b,se=se,t=t,r2=r2,s2=s2,omega=omega,res=res);
}



B


_________________________________________________________________
[[elided Hotmail spam]]

	[[alternative HTML version deleted]]

______________________________________________
R-help at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.