# install.packages(c("wooldridge", "plm", "texreg"), dependencies = TRUE)
library(wooldridge)
data(wagepan)
library(plm)

POLS <- plm(lwage ~ educ + black + hisp + exper+I(exper^2)+ married + union +
            factor(year), data = wagepan, index=c("nr","year") , model="pooling")

RE <- plm(lwage ~ educ + black + hisp + exper + I(exper^2) + married + union +
          factor(year), data = wagepan, index = c("nr","year") , model = "random")

FE <- plm(lwage ~ I(exper^2) + married + union + factor(year),
          data = wagepan, index = c("nr","year"), model="within")

# library(texreg)
texreg::screenreg(list(POLS, RE, FE), custom.coef.map = list('married' = 'Marrtied', 'union' = 'Union'))
#> ================================================
#>            Model 1      Model 2      Model 3
#> ------------------------------------------------
#> Marrtied      0.11 ***     0.06 ***     0.05 *
#>              (0.02)       (0.02)       (0.02)
#> Union         0.18 ***     0.11 ***     0.08 ***
#>              (0.02)       (0.02)       (0.02)
#> ------------------------------------------------
#> R^2           0.19         0.18         0.18
#> Adj. R^2      0.19         0.18         0.06
#> Num. obs.  4360         4360         4360
#> ================================================
#> *** p < 0.001, ** p < 0.01, * p < 0.05

#> ------------------------------------------------
#> Obs.  (N)   4360         4360         4360
#> Indiv.(n)    545          545          545
#> Time  (T)      8            8            8
#> R^2           0.19         0.18         0.18
#> Adj. R^2      0.19         0.18         0.06
#> F-stat       72.458       68.4124      83.8515
#> P-value        (2.22e-16)   (2.22e-16)   (2.22e-16)
#> ================================================
#> *** p < 0.001, ** p < 0.01, * p < 0.05

library(texreg)

fv.1 <- summary(POLS)$fstatistic$statistic
pv <- summary(POLS)$fstatistic$p.value
rbind(f, p)

ex.1 <- extract(POLS)
ex.1@gof.names <- c(ex.1@gof.names, "F-stat", "P-value")
ex.1@gof <- c(ex.1@gof, fv.1, pv)
ex.1@gof.decimal <- c(ex.1@gof.decimal, TRUE, TRUE)


fv.2 <- summary(RE)$fstatistic$statistic
pv <- summary(RE)$fstatistic$p.value
rbind(f, p)

ex.2 <- extract(RE)
ex.2@gof.names <- c(ex.2@gof.names, "F-stat", "P-value")
ex.2@gof <- c(ex.2@gof, fv.2, pv)
ex.2@gof.decimal <- c(ex.2@gof.decimal, TRUE, TRUE)


fv.3 <- summary(FE)$fstatistic$statistic
pv <- summary(FE)$fstatistic$p.value
rbind(f, p)

ex.3 <- extract(FE)
ex.3@gof.names <- c(ex.3@gof.names, "F-stat", "P-value")
ex.3@gof <- c(ex.3@gof, fv.3, pv)
ex.3@gof.decimal <- c(ex.3@gof.decimal, TRUE, TRUE)

> screenreg(list(ex.1, ex.2, ex.2))

=======================================================
                  Model 1      Model 2      Model 3
-------------------------------------------------------
(Intercept)          0.09         0.02         0.02
                    (0.08)       (0.15)       (0.15)
educ                 0.09 ***     0.09 ***     0.09 ***
                    (0.01)       (0.01)       (0.01)
black               -0.14 ***    -0.14 **     -0.14 **
                    (0.02)       (0.05)       (0.05)
hisp                 0.02         0.02         0.02
                    (0.02)       (0.04)       (0.04)
exper                0.07 ***     0.11 ***     0.11 ***
                    (0.01)       (0.02)       (0.02)
I(exper^2)          -0.00 **     -0.00 ***    -0.00 ***
                    (0.00)       (0.00)       (0.00)
married              0.11 ***     0.06 ***     0.06 ***
                    (0.02)       (0.02)       (0.02)
union                0.18 ***     0.11 ***     0.11 ***
                    (0.02)       (0.02)       (0.02)
[...]
-------------------------------------------------------
R^2                  0.19         0.18         0.18
Adj. R^2             0.19         0.18         0.18
Num. obs.         4360         4360         4360
F-stat              72.46        68.41        68.41
P-value              0.00         0.00         0.00
=======================================================
*** p < 0.001, ** p < 0.01, * p < 0.05