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-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On Behalf Of parkbomee Sent: Saturday, 7 February 2009 1:57 PM To: r-help@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@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. ______________________________________________ R-help@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.