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)
wagepan<-read.dta("http://fmwww.bc.edu/ec-p/data/wooldridge/wagepan.dta")

# 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           
Adjusted R2            0.049                     0.560           
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)? 
 A: 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
