# R-squared and F-stat in dummy variables regression vs panel FE model

When estimating a Fixed Effects model on panel data and an equivalent dummy variables regression, the coefficient estimates and associated SEs are identical. However, the R-squared and F-statistic are noticeably different (e.g. R-sq from dummy regression is usually much higher than R-sq from FE specification).

Once we obtain the R-squared & F-stat from the dummy variables regression, how can one adjust them to retrieve the same results as from the FE specification?

Consider this example:

library(foreign);library(plm);library(stargazer)

# Generate pdata.frame:
wagepan.p <- pdata.frame(wagepan, index=c("nr","year") )

# Estimate FE parameter in 3 different ways:
wagepan.p$yr<-factor(wagepan.p$year)

# Estimate dummy vars and FE models
reg.fe <-(plm(lwage~married+union+yr*educ,data=wagepan.p, model="within"))
reg.dum<-( lm(lwage~married+union+yr*educ+factor(nr), data=wagepan.p))

stargazer(reg.fe,reg.dum,type="text",model.names=FALSE,
keep=c("married","union"),omit.stat=c("ser"),
column.labels=c("Within","Dummies"))


Which will yield:

=================================================================
Dependent variable:
----------------------------------------------------
lwage
Within                    Dummies
(1)                       (2)
-----------------------------------------------------------------
married              0.055***                   0.055***
(0.018)                   (0.018)

union                0.083***                   0.083***
(0.019)                   (0.019)

-----------------------------------------------------------------
Observations           4,360                     4,360
R2                     0.171                     0.616
F Statistic  48.907*** (df = 16; 3799) 10.900*** (df = 560; 3799)
=================================================================
Note:                                 *p<0.1; **p<0.05; ***p<0.01


How can I adjust the Model 2 R2 and F-stat (0.616 and 10.9, respectively) to retrieve the same figures as in Model 1 (0.171 and 48.9)?

To get the same results for the F-Test and the R^2 with the lm() function you have to run the regression on the demeaned variables and adjust the available degrees of freedom. I have not looked up what plm() is doing in detail so results are not exact equal but pretty close (I do not know for example what is plm doing with 'educ' which is time-invariant and for which demeaning does not make sense). It does not matter really wether 'year' is demeaned or not since it is a balance panel here.

# Deman all variables (except 'educ' which is time-invariant)
wagepan$$lwage_de<-wagepan$$lwage-ave(wagepan$$lwage,factor(wagepan$$nr),FUN=mean)
wagepan$$married_de<-wagepan$$married-ave(wagepan$$married,factor(wagepan$$nr),FUN=mean)
wagepan$$union_de<-wagepan$$union-ave(wagepan$$union,factor(wagepan$$nr),FUN=mean)
wagepan$$year_de<-factor(wagepan$$year-ave(wagepan$$year,factor(wagepan$$nr),FUN=mean))
wagepan$$educ_de<-wagepan$$educ-ave(wagepan$$educ,factor(wagepan$$nr),FUN=mean)

# Define dummy variables for the years or 'yr80:educ' is also included
wagepan$$yr81<-as.numeric(wagepan$$year==1981)
wagepan$$yr82<-as.numeric(wagepan$$year==1982)
wagepan$$yr83<-as.numeric(wagepan$$year==1983)
wagepan$$yr84<-as.numeric(wagepan$$year==1984)
wagepan$$yr85<-as.numeric(wagepan$$year==1985)
wagepan$$yr86<-as.numeric(wagepan$$year==1986)
wagepan$$yr87<-as.numeric(wagepan$$year==1987)

# Run Regression
reg.fe2<-( lm(lwage_de~married_de+union_de+year_de+yr81:educ+yr82:educ+yr83:educ+yr84:educ+yr85:educ+yr86:educ+yr87:educ, data=wagepan))

# F-stat from the lm is wrong because it does not adjust degrees of freedom so calculate here
ss_res<-sum(reg.fe2$$residuals^2) ss_reg<-sum((reg.fe2$$fitted.values-mean(wagepan$$lwage_de))^2) ss_tot<-sum((wagepan$$lwage_de-mean(wagepan$$lwage_de))^2) p<-reg.fe2$$qr\$rank-1

# We have to adjust the degrees of freedom available
n<-length(wagepan$$lwage_de)-(reg.fe2$$qr$$rank+length(unique(wagepan$$nr))-1)

# R Squared
round(1 - (ss_res/ss_tot),2)

# F-Test
round((ss_reg/p)/(ss_res/n),2)


Here are the results:

> # R Squared
> round(1 - (ss_res/ss_tot),2)
[1] 0.17

> # F-Test
> round((ss_reg/p)/(ss_res/n),2)
[1] 48.59


For a nicer approach using matrix notation see for example

http://www.econ.uiuc.edu/~econ508/R/e-ta10_R.html

or

http://faculty.washington.edu/ezivot/econ582/fixedEffects.pdf

• I'm curious though if it is possible to adjust lm estimates without having to re-estimate the model on demeaned data. I am running a number of regressions in a loop, and re-estimating each model would prove time-consuming. Which is why I want to see if there is a way to somehow adjust 0.616 and 10.9 to retrieve the expected estimates... – landroni Feb 18 '16 at 22:32
• Of course, you can change the 'lm' function to account for this fact, but in this case it would probably make sense to also adjust standard error (to account for the clustering due to the panel structure of the data). Winthin estimation and Dummy variable approaches give the same results, so in the end why not use the plm function? You can even estimate coefficients for dummies using the within approach! – Arne Jonas Warnke Feb 19 '16 at 7:16
• I have issues with plm because it handles poorly cases of linear dependence in the model frame (see detect_lin_dep) which can generate errors in summary and/or vcovHC methods. Since the coefficient estimates are the same, I've switched to lm. For vcovHC.plm-like SEs the exact same Arellano estimator is implemented in multiwayvcov (e.g. cluster.vcov(mod.lm, ~ var1, df_correction=F)). Plus lm objs are more widely supported than plm in packages like car or lsmeans.. – landroni Feb 19 '16 at 12:30