Fixed and Mixed Effect Models
This code was adapted from those provided by Gelman and Hill and Heiss.
Reading in Gelman and Hill’s radon data
urlfile="http://www.stat.columbia.edu/~gelman/arm/examples/radon/srrs2.dat"
srrs2=read.csv(url(urlfile),header=T)
mn <- srrs2$state=="MN"
radon <- srrs2$activity[mn]
log.radon <- log (ifelse (radon==0, .1, radon))
floor <- srrs2$floor[mn] # 0 for basement, 1 for first floor
n <- length(radon)
y <- log.radon
x <- floor
hist(radon)
Creating county index
# get county index variable
county.name <- as.vector(srrs2$county[mn])
uniq <- unique(county.name)
uniq## [1] "AITKIN " "ANOKA " "BECKER "
## [4] "BELTRAMI " "BENTON " "BIG STONE "
## [7] "BLUE EARTH " "BROWN " "CARLTON "
## [10] "CARVER " "CASS " "CHIPPEWA "
## [13] "CHISAGO " "CLAY " "CLEARWATER "
## [16] "COOK " "COTTONWOOD " "CROW WING "
## [19] "DAKOTA " "DODGE " "DOUGLAS "
## [22] "FARIBAULT " "FILLMORE " "FREEBORN "
## [25] "GOODHUE " "HENNEPIN " "HOUSTON "
## [28] "HUBBARD " "ISANTI " "ITASCA "
## [31] "JACKSON " "KANABEC " "KANDIYOHI "
## [34] "KITTSON " "KOOCHICHING " "LAC QUI PARLE "
## [37] "LAKE " "LAKE OF THE WOODS " "LE SUEUR "
## [40] "LINCOLN " "LYON " "MAHNOMEN "
## [43] "MARSHALL " "MARTIN " "MCLEOD "
## [46] "MEEKER " "MILLE LACS " "MORRISON "
## [49] "MOWER " "MURRAY " "NICOLLET "
## [52] "NOBLES " "NORMAN " "OLMSTED "
## [55] "OTTER TAIL " "PENNINGTON " "PINE "
## [58] "PIPESTONE " "POLK " "POPE "
## [61] "RAMSEY " "REDWOOD " "RENVILLE "
## [64] "RICE " "ROCK " "ROSEAU "
## [67] "SCOTT " "SHERBURNE " "SIBLEY "
## [70] "ST LOUIS " "STEARNS " "STEELE "
## [73] "STEVENS " "SWIFT " "TODD "
## [76] "TRAVERSE " "WABASHA " "WADENA "
## [79] "WASECA " "WASHINGTON " "WATONWAN "
## [82] "WILKIN " "WINONA " "WRIGHT "
## [85] "YELLOW MEDICINE "J <- length(uniq)
county <- rep (NA, J)
for (i in 1:J){
county[county.name==uniq[i]] <- i
}Completely pooled versus partial pooling
library(texreg)## Version: 1.38.6
## Date: 2022-04-06
## Author: Philip Leifeld (University of Essex)
##
## Consider submitting praise using the praise or praise_interactive functions.
## Please cite the JSS article in your publications -- see citation("texreg").# Linear model of log(radon) on dummy of first floor (ref: basement)
## Complete pooling regression
lm.pooled <- lm (y ~ x)
summary(lm.pooled)##
## Call:
## lm(formula = y ~ x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.6293 -0.5383 0.0342 0.5603 2.5486
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.32674 0.02972 44.640 <2e-16 ***
## x -0.61339 0.07284 -8.421 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8226 on 917 degrees of freedom
## Multiple R-squared: 0.07178, Adjusted R-squared: 0.07077
## F-statistic: 70.91 on 1 and 917 DF, p-value: < 2.2e-16## Partial pooling regression . )-1 takes out the intercept)
lm.unpooled <- lm (y ~ x + factor(county) -1)
summary(lm.unpooled)##
## Call:
## lm(formula = y ~ x + factor(county) - 1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.14595 -0.45405 0.00065 0.45376 2.65987
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## x -0.72054 0.07352 -9.800 < 2e-16 ***
## factor(county)1 0.84054 0.37866 2.220 0.026701 *
## factor(county)2 0.87482 0.10498 8.333 3.23e-16 ***
## factor(county)3 1.52870 0.43946 3.479 0.000530 ***
## factor(county)4 1.55272 0.28897 5.373 1.00e-07 ***
## factor(county)5 1.43257 0.37866 3.783 0.000166 ***
## factor(county)6 1.51301 0.43672 3.464 0.000558 ***
## factor(county)7 2.01216 0.20243 9.940 < 2e-16 ***
## factor(county)8 1.98958 0.37999 5.236 2.08e-07 ***
## factor(county)9 1.00304 0.23931 4.191 3.07e-05 ***
## factor(county)10 1.56391 0.31099 5.029 6.04e-07 ***
## factor(county)11 1.40113 0.33828 4.142 3.80e-05 ***
## factor(county)12 1.73025 0.37821 4.575 5.49e-06 ***
## factor(county)13 1.03872 0.30881 3.364 0.000804 ***
## factor(county)14 1.98838 0.20325 9.783 < 2e-16 ***
## factor(county)15 1.33797 0.37999 3.521 0.000453 ***
## factor(county)16 0.66486 0.53487 1.243 0.214204
## factor(county)17 1.27480 0.38221 3.335 0.000890 ***
## factor(county)18 1.12155 0.21913 5.118 3.83e-07 ***
## factor(county)19 1.33831 0.09541 14.026 < 2e-16 ***
## factor(county)20 1.80032 0.43672 4.122 4.13e-05 ***
## factor(county)21 1.73399 0.25227 6.873 1.23e-11 ***
## factor(county)22 0.63679 0.30905 2.060 0.039663 *
## factor(county)23 1.39999 0.53613 2.611 0.009183 **
## factor(county)24 2.10162 0.25267 8.318 3.64e-16 ***
## factor(county)25 1.95072 0.20243 9.636 < 2e-16 ***
## factor(county)26 1.36058 0.07422 18.332 < 2e-16 ***
## factor(county)27 1.77336 0.30978 5.725 1.45e-08 ***
## factor(county)28 1.24159 0.34115 3.639 0.000290 ***
## factor(county)29 1.05600 0.43672 2.418 0.015818 *
## factor(county)30 0.92576 0.22807 4.059 5.39e-05 ***
## factor(county)31 2.02057 0.33828 5.973 3.45e-09 ***
## factor(county)32 1.23629 0.37821 3.269 0.001124 **
## factor(county)33 2.06187 0.37821 5.452 6.58e-08 ***
## factor(county)34 1.59044 0.43946 3.619 0.000314 ***
## factor(county)35 0.81920 0.28897 2.835 0.004695 **
## factor(county)36 2.95897 0.53613 5.519 4.55e-08 ***
## factor(county)37 0.40209 0.25227 1.594 0.111345
## factor(county)38 1.86772 0.37999 4.915 1.07e-06 ***
## factor(county)39 1.74807 0.33860 5.163 3.05e-07 ***
## factor(county)40 2.31580 0.37866 6.116 1.48e-09 ***
## factor(county)41 1.96715 0.26759 7.351 4.69e-13 ***
## factor(county)42 1.36098 0.75642 1.799 0.072343 .
## factor(county)43 1.60224 0.25543 6.273 5.69e-10 ***
## factor(county)44 1.04099 0.28609 3.639 0.000291 ***
## factor(county)45 1.29541 0.21101 6.139 1.28e-09 ***
## factor(county)46 1.21461 0.33828 3.591 0.000349 ***
## factor(county)47 0.88393 0.53613 1.649 0.099583 .
## factor(county)48 1.14812 0.25227 4.551 6.13e-06 ***
## factor(county)49 1.70211 0.21010 8.102 1.93e-15 ***
## factor(county)50 2.49321 0.75642 3.296 0.001022 **
## factor(county)51 2.16504 0.37821 5.724 1.45e-08 ***
## factor(county)52 1.92769 0.43672 4.414 1.15e-05 ***
## factor(county)53 1.25080 0.43741 2.860 0.004348 **
## factor(county)54 1.30676 0.15802 8.270 5.28e-16 ***
## factor(county)55 1.61799 0.26885 6.018 2.64e-09 ***
## factor(county)56 1.10110 0.43946 2.506 0.012415 *
## factor(county)57 0.76218 0.30905 2.466 0.013855 *
## factor(county)58 1.86092 0.37866 4.915 1.07e-06 ***
## factor(county)59 1.72178 0.37999 4.531 6.73e-06 ***
## factor(county)60 1.27939 0.53487 2.392 0.016979 *
## factor(county)61 1.15873 0.13389 8.654 < 2e-16 ***
## factor(county)62 1.98301 0.33860 5.856 6.80e-09 ***
## factor(county)63 1.67070 0.43741 3.820 0.000144 ***
## factor(county)64 1.84784 0.22817 8.099 1.97e-15 ***
## factor(county)65 1.29912 0.53487 2.429 0.015357 *
## factor(county)66 1.66574 0.20648 8.067 2.50e-15 ***
## factor(county)67 1.80312 0.21101 8.545 < 2e-16 ***
## factor(county)68 1.09002 0.26743 4.076 5.02e-05 ***
## factor(county)69 1.24245 0.37821 3.285 0.001062 **
## factor(county)70 0.86763 0.07096 12.227 < 2e-16 ***
## factor(county)71 1.49184 0.15174 9.832 < 2e-16 ***
## factor(county)72 1.57990 0.23920 6.605 7.08e-11 ***
## factor(county)73 1.79176 0.53487 3.350 0.000845 ***
## factor(county)74 0.98704 0.37821 2.610 0.009223 **
## factor(county)75 1.72372 0.43741 3.941 8.80e-05 ***
## factor(county)76 2.00844 0.37866 5.304 1.45e-07 ***
## factor(county)77 1.82168 0.28609 6.367 3.17e-10 ***
## factor(county)78 1.28569 0.33956 3.786 0.000164 ***
## factor(county)79 0.61488 0.37866 1.624 0.104785
## factor(county)80 1.32952 0.11181 11.890 < 2e-16 ***
## factor(county)81 2.70953 0.43946 6.166 1.09e-09 ***
## factor(county)82 2.23001 0.75642 2.948 0.003286 **
## factor(county)83 1.62292 0.21048 7.711 3.57e-14 ***
## factor(county)84 1.64535 0.20987 7.840 1.38e-14 ***
## factor(county)85 1.18652 0.53487 2.218 0.026801 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7564 on 833 degrees of freedom
## Multiple R-squared: 0.7671, Adjusted R-squared: 0.7431
## F-statistic: 31.91 on 86 and 833 DF, p-value: < 2.2e-16lm.unpooled2 <- lm (y ~ x + factor(county))
screenreg(list(lm.pooled, lm.unpooled, lm.unpooled2))##
## ====================================================
## Model 1 Model 2 Model 3
## ----------------------------------------------------
## (Intercept) 1.33 *** 0.84 *
## (0.03) (0.38)
## x -0.61 *** -0.72 *** -0.72 ***
## (0.07) (0.07) (0.07)
## factor(county)1 0.84 *
## (0.38)
## factor(county)2 0.87 *** 0.03
## (0.10) (0.39)
## factor(county)3 1.53 *** 0.69
## (0.44) (0.58)
## factor(county)4 1.55 *** 0.71
## (0.29) (0.47)
## factor(county)5 1.43 *** 0.59
## (0.38) (0.53)
## factor(county)6 1.51 *** 0.67
## (0.44) (0.58)
## factor(county)7 2.01 *** 1.17 **
## (0.20) (0.43)
## factor(county)8 1.99 *** 1.15 *
## (0.38) (0.54)
## factor(county)9 1.00 *** 0.16
## (0.24) (0.45)
## factor(county)10 1.56 *** 0.72
## (0.31) (0.49)
## factor(county)11 1.40 *** 0.56
## (0.34) (0.51)
## factor(county)12 1.73 *** 0.89
## (0.38) (0.54)
## factor(county)13 1.04 *** 0.20
## (0.31) (0.49)
## factor(county)14 1.99 *** 1.15 **
## (0.20) (0.43)
## factor(county)15 1.34 *** 0.50
## (0.38) (0.54)
## factor(county)16 0.66 -0.18
## (0.53) (0.66)
## factor(county)17 1.27 *** 0.43
## (0.38) (0.54)
## factor(county)18 1.12 *** 0.28
## (0.22) (0.44)
## factor(county)19 1.34 *** 0.50
## (0.10) (0.39)
## factor(county)20 1.80 *** 0.96
## (0.44) (0.58)
## factor(county)21 1.73 *** 0.89 *
## (0.25) (0.45)
## factor(county)22 0.64 * -0.20
## (0.31) (0.49)
## factor(county)23 1.40 ** 0.56
## (0.54) (0.66)
## factor(county)24 2.10 *** 1.26 **
## (0.25) (0.45)
## factor(county)25 1.95 *** 1.11 **
## (0.20) (0.43)
## factor(county)26 1.36 *** 0.52
## (0.07) (0.39)
## factor(county)27 1.77 *** 0.93
## (0.31) (0.49)
## factor(county)28 1.24 *** 0.40
## (0.34) (0.51)
## factor(county)29 1.06 * 0.22
## (0.44) (0.58)
## factor(county)30 0.93 *** 0.09
## (0.23) (0.44)
## factor(county)31 2.02 *** 1.18 *
## (0.34) (0.51)
## factor(county)32 1.24 ** 0.40
## (0.38) (0.54)
## factor(county)33 2.06 *** 1.22 *
## (0.38) (0.54)
## factor(county)34 1.59 *** 0.75
## (0.44) (0.58)
## factor(county)35 0.82 ** -0.02
## (0.29) (0.47)
## factor(county)36 2.96 *** 2.12 **
## (0.54) (0.66)
## factor(county)37 0.40 -0.44
## (0.25) (0.45)
## factor(county)38 1.87 *** 1.03
## (0.38) (0.54)
## factor(county)39 1.75 *** 0.91
## (0.34) (0.51)
## factor(county)40 2.32 *** 1.48 **
## (0.38) (0.53)
## factor(county)41 1.97 *** 1.13 *
## (0.27) (0.46)
## factor(county)42 1.36 0.52
## (0.76) (0.85)
## factor(county)43 1.60 *** 0.76
## (0.26) (0.46)
## factor(county)44 1.04 *** 0.20
## (0.29) (0.47)
## factor(county)45 1.30 *** 0.45
## (0.21) (0.43)
## factor(county)46 1.21 *** 0.37
## (0.34) (0.51)
## factor(county)47 0.88 0.04
## (0.54) (0.66)
## factor(county)48 1.15 *** 0.31
## (0.25) (0.45)
## factor(county)49 1.70 *** 0.86 *
## (0.21) (0.43)
## factor(county)50 2.49 ** 1.65
## (0.76) (0.85)
## factor(county)51 2.17 *** 1.32 *
## (0.38) (0.54)
## factor(county)52 1.93 *** 1.09
## (0.44) (0.58)
## factor(county)53 1.25 ** 0.41
## (0.44) (0.58)
## factor(county)54 1.31 *** 0.47
## (0.16) (0.41)
## factor(county)55 1.62 *** 0.78
## (0.27) (0.46)
## factor(county)56 1.10 * 0.26
## (0.44) (0.58)
## factor(county)57 0.76 * -0.08
## (0.31) (0.49)
## factor(county)58 1.86 *** 1.02
## (0.38) (0.53)
## factor(county)59 1.72 *** 0.88
## (0.38) (0.54)
## factor(county)60 1.28 * 0.44
## (0.53) (0.66)
## factor(county)61 1.16 *** 0.32
## (0.13) (0.40)
## factor(county)62 1.98 *** 1.14 *
## (0.34) (0.51)
## factor(county)63 1.67 *** 0.83
## (0.44) (0.58)
## factor(county)64 1.85 *** 1.01 *
## (0.23) (0.44)
## factor(county)65 1.30 * 0.46
## (0.53) (0.66)
## factor(county)66 1.67 *** 0.83
## (0.21) (0.43)
## factor(county)67 1.80 *** 0.96 *
## (0.21) (0.43)
## factor(county)68 1.09 *** 0.25
## (0.27) (0.46)
## factor(county)69 1.24 ** 0.40
## (0.38) (0.54)
## factor(county)70 0.87 *** 0.03
## (0.07) (0.38)
## factor(county)71 1.49 *** 0.65
## (0.15) (0.41)
## factor(county)72 1.58 *** 0.74
## (0.24) (0.45)
## factor(county)73 1.79 *** 0.95
## (0.53) (0.66)
## factor(county)74 0.99 ** 0.15
## (0.38) (0.54)
## factor(county)75 1.72 *** 0.88
## (0.44) (0.58)
## factor(county)76 2.01 *** 1.17 *
## (0.38) (0.53)
## factor(county)77 1.82 *** 0.98 *
## (0.29) (0.47)
## factor(county)78 1.29 *** 0.45
## (0.34) (0.51)
## factor(county)79 0.61 -0.23
## (0.38) (0.53)
## factor(county)80 1.33 *** 0.49
## (0.11) (0.39)
## factor(county)81 2.71 *** 1.87 **
## (0.44) (0.58)
## factor(county)82 2.23 ** 1.39
## (0.76) (0.85)
## factor(county)83 1.62 *** 0.78
## (0.21) (0.43)
## factor(county)84 1.65 *** 0.80
## (0.21) (0.43)
## factor(county)85 1.19 * 0.35
## (0.53) (0.66)
## ----------------------------------------------------
## R^2 0.07 0.77 0.29
## Adj. R^2 0.07 0.74 0.21
## Num. obs. 919 919 919
## ====================================================
## *** p < 0.001; ** p < 0.01; * p < 0.05– Being on the first floor versus basement is estimated to reduce radon exposure by 61%.
Estimating random-intercept models
library(lme4) # linear fixed effects model## Loading required package: Matrixlibrary(sjstats)## Registered S3 method overwritten by 'parameters':
## method from
## format.parameters_distribution datawizard## Varying-intercept model w/ no predictors
?lmer## starting httpd help server ...## doneM0 <- lmer (y ~ 1 + (1 | county))
summary(M0)## Linear mixed model fit by REML ['lmerMod']
## Formula: y ~ 1 + (1 | county)
##
## REML criterion at convergence: 2259.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.4661 -0.5734 0.0441 0.6432 3.3516
##
## Random effects:
## Groups Name Variance Std.Dev.
## county (Intercept) 0.09581 0.3095
## Residual 0.63662 0.7979
## Number of obs: 919, groups: county, 85
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.31258 0.04891 26.84# from sjstats package
performance::icc(M0)## # Intraclass Correlation Coefficient
##
## Adjusted ICC: 0.131
## Conditional ICC: 0.131screenreg(M0)##
## =====================================
## Model 1
## -------------------------------------
## (Intercept) 1.31 ***
## (0.05)
## -------------------------------------
## AIC 2265.44
## BIC 2279.91
## Log Likelihood -1129.72
## Num. obs. 919
## Num. groups: county 85
## Var: county (Intercept) 0.10
## Var: Residual 0.64
## =====================================
## *** p < 0.001; ** p < 0.01; * p < 0.05Question: How much of the variance is accounted for by the county level?
Answer: 13% of the variance in model 2 is attributed to level 2 of the county. “(0.09581/(0.09581 + 0.63662))”
Random-intercept model with one predictor
## Including x as a predictor
M1 <- lmer (y ~ x + (1 | county))
summary(M1)## Linear mixed model fit by REML ['lmerMod']
## Formula: y ~ x + (1 | county)
##
## REML criterion at convergence: 2171.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.3989 -0.6155 0.0029 0.6405 3.4281
##
## Random effects:
## Groups Name Variance Std.Dev.
## county (Intercept) 0.1077 0.3282
## Residual 0.5709 0.7556
## Number of obs: 919, groups: county, 85
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.46160 0.05158 28.339
## x -0.69299 0.07043 -9.839
##
## Correlation of Fixed Effects:
## (Intr)
## x -0.288# estimated regression coefficicents
coef (M1)## $county
## (Intercept) x
## 1 1.1915004 -0.6929937
## 2 0.9276468 -0.6929937
## 3 1.4792143 -0.6929937
## 4 1.5045012 -0.6929937
## 5 1.4461503 -0.6929937
## 6 1.4801817 -0.6929937
## 7 1.8581255 -0.6929937
## 8 1.6827736 -0.6929937
## 9 1.1600747 -0.6929937
## 10 1.5086099 -0.6929937
## 11 1.4322449 -0.6929937
## 12 1.5771520 -0.6929937
## 13 1.2370518 -0.6929937
## 14 1.8380232 -0.6929937
## 15 1.4024982 -0.6929937
## 16 1.2432992 -0.6929937
## 17 1.3723633 -0.6929937
## 18 1.2209415 -0.6929937
## 19 1.3462611 -0.6929937
## 20 1.5840332 -0.6929937
## 21 1.6311136 -0.6929937
## 22 1.0211903 -0.6929937
## 23 1.4409443 -0.6929937
## 24 1.8605721 -0.6929937
## 25 1.8135585 -0.6929937
## 26 1.3626875 -0.6929937
## 27 1.6222663 -0.6929937
## 28 1.3467692 -0.6929937
## 29 1.3149879 -0.6929937
## 30 1.0999775 -0.6929937
## 31 1.7329562 -0.6929937
## 32 1.3646863 -0.6929937
## 33 1.7197950 -0.6929937
## 34 1.5015319 -0.6929937
## 35 1.0870316 -0.6929937
## 36 1.8680899 -0.6929937
## 37 0.7928242 -0.6929937
## 38 1.6303574 -0.6929937
## 39 1.5979922 -0.6929937
## 40 1.8260564 -0.6929937
## 41 1.7636308 -0.6929937
## 42 1.4456250 -0.6929937
## 43 1.5404841 -0.6929937
## 44 1.2199767 -0.6929937
## 45 1.3375197 -0.6929937
## 46 1.3416955 -0.6929937
## 47 1.2995481 -0.6929937
## 48 1.2623708 -0.6929937
## 49 1.6294468 -0.6929937
## 50 1.6253580 -0.6929937
## 51 1.7641693 -0.6929937
## 52 1.6300755 -0.6929937
## 53 1.3820836 -0.6929937
## 54 1.3328317 -0.6929937
## 55 1.5494611 -0.6929937
## 56 1.3246513 -0.6929937
## 57 1.0877739 -0.6929937
## 58 1.6303984 -0.6929937
## 59 1.5675867 -0.6929937
## 60 1.4116740 -0.6929937
## 61 1.1995431 -0.6929937
## 62 1.7120492 -0.6929937
## 63 1.5338618 -0.6929937
## 64 1.7205672 -0.6929937
## 65 1.4170797 -0.6929937
## 66 1.5982671 -0.6929937
## 67 1.6981946 -0.6929937
## 68 1.2380854 -0.6929937
## 69 1.3673371 -0.6929937
## 70 0.8899487 -0.6929937
## 71 1.4829168 -0.6929937
## 72 1.5389227 -0.6929937
## 73 1.5520593 -0.6929937
## 74 1.2574763 -0.6929937
## 75 1.5530274 -0.6929937
## 76 1.6938490 -0.6929937
## 77 1.6642922 -0.6929937
## 78 1.3708520 -0.6929937
## 79 1.0944378 -0.6929937
## 80 1.3404792 -0.6929937
## 81 1.9060482 -0.6929937
## 82 1.5835784 -0.6929937
## 83 1.5716875 -0.6929937
## 84 1.5906331 -0.6929937
## 85 1.3862294 -0.6929937
##
## attr(,"class")
## [1] "coef.mer"# fixed and random effects
fixef (M1)## (Intercept) x
## 1.4615979 -0.6929937ranef (M1)## $county
## (Intercept)
## 1 -0.27009750
## 2 -0.53395107
## 3 0.01761646
## 4 0.04290331
## 5 -0.01544759
## 6 0.01858386
## 7 0.39652761
## 8 0.22117571
## 9 -0.30152320
## 10 0.04701206
## 11 -0.02935302
## 12 0.11555412
## 13 -0.22454606
## 14 0.37642528
## 15 -0.05909970
## 16 -0.21829868
## 17 -0.08923455
## 18 -0.24065633
## 19 -0.11533676
## 20 0.12243536
## 21 0.16951571
## 22 -0.44040759
## 23 -0.02065353
## 24 0.39897419
## 25 0.35196058
## 26 -0.09891042
## 27 0.16066845
## 28 -0.11482866
## 29 -0.14661001
## 30 -0.36162037
## 31 0.27135836
## 32 -0.09691159
## 33 0.25819715
## 34 0.03993398
## 35 -0.37456628
## 36 0.40649202
## 37 -0.66877366
## 38 0.16875953
## 39 0.13639435
## 40 0.36445852
## 41 0.30203292
## 42 -0.01597290
## 43 0.07888623
## 44 -0.24162113
## 45 -0.12407819
## 46 -0.11990235
## 47 -0.16204982
## 48 -0.19922711
## 49 0.16784892
## 50 0.16376017
## 51 0.30257145
## 52 0.16847763
## 53 -0.07951423
## 54 -0.12876622
## 55 0.08786324
## 56 -0.13694655
## 57 -0.37382399
## 58 0.16880050
## 59 0.10598884
## 60 -0.04992385
## 61 -0.26205479
## 62 0.25045135
## 63 0.07226395
## 64 0.25896934
## 65 -0.04451816
## 66 0.13666918
## 67 0.23659674
## 68 -0.22351243
## 69 -0.09426076
## 70 -0.57164915
## 71 0.02131895
## 72 0.07732481
## 73 0.09046140
## 74 -0.20412162
## 75 0.09142952
## 76 0.23225113
## 77 0.20269435
## 78 -0.09074583
## 79 -0.36716012
## 80 -0.12111867
## 81 0.44445030
## 82 0.12198051
## 83 0.11008961
## 84 0.12903524
## 85 -0.07536845
##
## with conditional variances for "county"Question: What are the estimated intercept for counties 1 and 2?
Answer: The estimated intercepts for county one is 1.2 (1.462 + (-0.27009750)) and county two is 0.93 (1.462 + (-0.53395107))
Adding a group-level variable
## Get the county-level predictor
srrs2.fips <- srrs2$stfips*1000 + srrs2$cntyfips
urlfile="http://www.stat.columbia.edu/~gelman/arm/examples/radon/cty.dat"
cty=read.csv(url(urlfile),header=T)
usa.fips <- 1000*cty[,"stfips"] + cty[,"ctfips"]
usa.rows <- match (unique(srrs2.fips[mn]), usa.fips)
uranium <- cty[usa.rows,"Uppm"]
u <- log (uranium)Adding a group-level predictor
library(lmtest)## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric## Varying-intercept model w/ group-level predictors
u.full <- u[county]
M2 <- lmer (y ~ x + u.full + (1 | county))
summary(M2)## Linear mixed model fit by REML ['lmerMod']
## Formula: y ~ x + u.full + (1 | county)
##
## REML criterion at convergence: 2134.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.9673 -0.6117 0.0274 0.6555 3.3848
##
## Random effects:
## Groups Name Variance Std.Dev.
## county (Intercept) 0.02446 0.1564
## Residual 0.57523 0.7584
## Number of obs: 919, groups: county, 85
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.46576 0.03794 38.633
## x -0.66824 0.06880 -9.713
## u.full 0.72027 0.09176 7.849
##
## Correlation of Fixed Effects:
## (Intr) x
## x -0.357
## u.full 0.145 -0.009# Goodness of fit
screenreg(list(M0,M1,M2))##
## =================================================================
## Model 1 Model 2 Model 3
## -----------------------------------------------------------------
## (Intercept) 1.31 *** 1.46 *** 1.47 ***
## (0.05) (0.05) (0.04)
## x -0.69 *** -0.67 ***
## (0.07) (0.07)
## u.full 0.72 ***
## (0.09)
## -----------------------------------------------------------------
## AIC 2265.44 2179.31 2144.19
## BIC 2279.91 2198.60 2168.30
## Log Likelihood -1129.72 -1085.65 -1067.09
## Num. obs. 919 919 919
## Num. groups: county 85 85 85
## Var: county (Intercept) 0.10 0.11 0.02
## Var: Residual 0.64 0.57 0.58
## =================================================================
## *** p < 0.001; ** p < 0.01; * p < 0.05lrtest(M0,M1)lrtest(M1,M2)# Looking at model results
coef (M2)## $county
## (Intercept) x u.full
## 1 1.445120 -0.6682448 0.7202676
## 2 1.477009 -0.6682448 0.7202676
## 3 1.478185 -0.6682448 0.7202676
## 4 1.576891 -0.6682448 0.7202676
## 5 1.473999 -0.6682448 0.7202676
## 6 1.439566 -0.6682448 0.7202676
## 7 1.593872 -0.6682448 0.7202676
## 8 1.509045 -0.6682448 0.7202676
## 9 1.397556 -0.6682448 0.7202676
## 10 1.466362 -0.6682448 0.7202676
## 11 1.531245 -0.6682448 0.7202676
## 12 1.475555 -0.6682448 0.7202676
## 13 1.486617 -0.6682448 0.7202676
## 14 1.562814 -0.6682448 0.7202676
## 15 1.449652 -0.6682448 0.7202676
## 16 1.431496 -0.6682448 0.7202676
## 17 1.396746 -0.6682448 0.7202676
## 18 1.499187 -0.6682448 0.7202676
## 19 1.383197 -0.6682448 0.7202676
## 20 1.482112 -0.6682448 0.7202676
## 21 1.507354 -0.6682448 0.7202676
## 22 1.252279 -0.6682448 0.7202676
## 23 1.435133 -0.6682448 0.7202676
## 24 1.593848 -0.6682448 0.7202676
## 25 1.591105 -0.6682448 0.7202676
## 26 1.432151 -0.6682448 0.7202676
## 27 1.451023 -0.6682448 0.7202676
## 28 1.471547 -0.6682448 0.7202676
## 29 1.480673 -0.6682448 0.7202676
## 30 1.445935 -0.6682448 0.7202676
## 31 1.524017 -0.6682448 0.7202676
## 32 1.437994 -0.6682448 0.7202676
## 33 1.540931 -0.6682448 0.7202676
## 34 1.476560 -0.6682448 0.7202676
## 35 1.456290 -0.6682448 0.7202676
## 36 1.563196 -0.6682448 0.7202676
## 37 1.307618 -0.6682448 0.7202676
## 38 1.591774 -0.6682448 0.7202676
## 39 1.488871 -0.6682448 0.7202676
## 40 1.540850 -0.6682448 0.7202676
## 41 1.519206 -0.6682448 0.7202676
## 42 1.457093 -0.6682448 0.7202676
## 43 1.492752 -0.6682448 0.7202676
## 44 1.339205 -0.6682448 0.7202676
## 45 1.363380 -0.6682448 0.7202676
## 46 1.418700 -0.6682448 0.7202676
## 47 1.429966 -0.6682448 0.7202676
## 48 1.394820 -0.6682448 0.7202676
## 49 1.480131 -0.6682448 0.7202676
## 50 1.495954 -0.6682448 0.7202676
## 51 1.541446 -0.6682448 0.7202676
## 52 1.484985 -0.6682448 0.7202676
## 53 1.417856 -0.6682448 0.7202676
## 54 1.297171 -0.6682448 0.7202676
## 55 1.536860 -0.6682448 0.7202676
## 56 1.426593 -0.6682448 0.7202676
## 57 1.344875 -0.6682448 0.7202676
## 58 1.471195 -0.6682448 0.7202676
## 59 1.471316 -0.6682448 0.7202676
## 60 1.435281 -0.6682448 0.7202676
## 61 1.459521 -0.6682448 0.7202676
## 62 1.508361 -0.6682448 0.7202676
## 63 1.455963 -0.6682448 0.7202676
## 64 1.541677 -0.6682448 0.7202676
## 65 1.422890 -0.6682448 0.7202676
## 66 1.586234 -0.6682448 0.7202676
## 67 1.563953 -0.6682448 0.7202676
## 68 1.495340 -0.6682448 0.7202676
## 69 1.408462 -0.6682448 0.7202676
## 70 1.246716 -0.6682448 0.7202676
## 71 1.431690 -0.6682448 0.7202676
## 72 1.441835 -0.6682448 0.7202676
## 73 1.464737 -0.6682448 0.7202676
## 74 1.363079 -0.6682448 0.7202676
## 75 1.496793 -0.6682448 0.7202676
## 76 1.490651 -0.6682448 0.7202676
## 77 1.520897 -0.6682448 0.7202676
## 78 1.515393 -0.6682448 0.7202676
## 79 1.317927 -0.6682448 0.7202676
## 80 1.442140 -0.6682448 0.7202676
## 81 1.587608 -0.6682448 0.7202676
## 82 1.490002 -0.6682448 0.7202676
## 83 1.398639 -0.6682448 0.7202676
## 84 1.551352 -0.6682448 0.7202676
## 85 1.423816 -0.6682448 0.7202676
##
## attr(,"class")
## [1] "coef.mer"fixef(M2)## (Intercept) x u.full
## 1.4657628 -0.6682448 0.7202676ranef(M2)## $county
## (Intercept)
## 1 -0.0206423630
## 2 0.0112457688
## 3 0.0124220276
## 4 0.1111284341
## 5 0.0082358459
## 6 -0.0261963570
## 7 0.1281093537
## 8 0.0432824990
## 9 -0.0682068835
## 10 0.0005992601
## 11 0.0654820919
## 12 0.0097919700
## 13 0.0208539173
## 14 0.0970507335
## 15 -0.0161106390
## 16 -0.0342668097
## 17 -0.0690169594
## 18 0.0334239015
## 19 -0.0825654644
## 20 0.0163495692
## 21 0.0415913432
## 22 -0.2134841935
## 23 -0.0306298131
## 24 0.1280847611
## 25 0.1253422380
## 26 -0.0336118129
## 27 -0.0147402932
## 28 0.0057840506
## 29 0.0149105894
## 30 -0.0198281643
## 31 0.0582539647
## 32 -0.0277690233
## 33 0.0751682451
## 34 0.0107974655
## 35 -0.0094725278
## 36 0.0974331099
## 37 -0.1581443594
## 38 0.1260113268
## 39 0.0231083697
## 40 0.0750874972
## 41 0.0534435511
## 42 -0.0086695336
## 43 0.0269893894
## 44 -0.1265576167
## 45 -0.1023824445
## 46 -0.0470625166
## 47 -0.0357963317
## 48 -0.0709428996
## 49 0.0143685026
## 50 0.0301915894
## 51 0.0756831087
## 52 0.0192220801
## 53 -0.0479063547
## 54 -0.1685913456
## 55 0.0710968071
## 56 -0.0391702639
## 57 -0.1208880435
## 58 0.0054318336
## 59 0.0055532559
## 60 -0.0304815132
## 61 -0.0062415829
## 62 0.0425977970
## 63 -0.0097996359
## 64 0.0759142226
## 65 -0.0428729911
## 66 0.1204715381
## 67 0.0981904668
## 68 0.0295769306
## 69 -0.0573007541
## 70 -0.2190464434
## 71 -0.0340728075
## 72 -0.0239274353
## 73 -0.0010260185
## 74 -0.1026833598
## 75 0.0310303716
## 76 0.0248879657
## 77 0.0551338148
## 78 0.0496306720
## 79 -0.1478361415
## 80 -0.0236229693
## 81 0.1218447545
## 82 0.0242393918
## 83 -0.0671242614
## 84 0.0855889222
## 85 -0.0419463713
##
## with conditional variances for "county"Question: Interpret the effects of the best-fitting model.
Answer: Model 2 is the best fitting model. The model included “county” as random effect. The model’s total explanatory power is moderate (conditional R2 = 0.23) and the part related to the fixed effects alone (marginal R2) is of 0.09. The model’s intercept, corresponding to x = 0, is at 1.46 (with 95% confidence, p < 0.05)
The effect of x is statistically significant and negative (beta = -0.69, p < 0.05)
Panel data
library(plm);library(stargazer)##
## Please cite as:## Hlavac, Marek (2022). stargazer: Well-Formatted Regression and Summary Statistics Tables.## R package version 5.2.3. https://CRAN.R-project.org/package=stargazerdata(wagepan, package='wooldridge')
# Generate pdata.frame:
wagepan.p <- pdata.frame(wagepan, index=c("nr","year") )
pdim(wagepan.p)## Balanced Panel: n = 545, T = 8, N = 4360# Check variation of variables within individuals
pvar(wagepan.p)## no time variation: nr black hisp educ
## no individual variation: year d81 d82 d83 d84 d85 d86 d87Fixed and random effects models
library(texreg)
# Estimate different models
wagepan.p$yr<-factor(wagepan.p$year)
reg.ols<- (plm(lwage~educ+black+hisp+exper+I(exper^2)+married+union+yr,
data=wagepan.p, model="pooling") )
reg.fd <- (plm(lwage~exper + I(exper^2)+married+union+yr,
data=wagepan.p, model="fd") )
reg.fe <- (plm(lwage~ exper+I(exper^2)+married+union+yr,
data=wagepan.p, model="within") )
reg.re <- (plm(lwage~ educ+black+hisp+exper+I(exper^2)+married+union+yr,
data=wagepan.p, model="random") )
screenreg(list(reg.ols,reg.fd,reg.fe,reg.re))##
## ===============================================================
## Model 1 Model 2 Model 3 Model 4
## ---------------------------------------------------------------
## (Intercept) 0.09 0.14 *** 0.02
## (0.08) (0.03) (0.15)
## educ 0.09 *** 0.09 ***
## (0.01) (0.01)
## black -0.14 *** -0.14 **
## (0.02) (0.05)
## hisp 0.02 0.02
## (0.02) (0.04)
## exper 0.07 *** 0.13 *** 0.11 ***
## (0.01) (0.01) (0.02)
## exper^2 -0.00 ** -0.01 ** -0.01 *** -0.00 ***
## (0.00) (0.00) (0.00) (0.00)
## married 0.11 *** 0.04 0.05 * 0.06 ***
## (0.02) (0.02) (0.02) (0.02)
## union 0.18 *** 0.04 * 0.08 *** 0.11 ***
## (0.02) (0.02) (0.02) (0.02)
## yr1981 0.06 0.02 0.02 0.04
## (0.03) (0.02) (0.02) (0.02)
## yr1982 0.06 -0.02 -0.01 0.03
## (0.03) (0.03) (0.02) (0.03)
## yr1983 0.06 -0.05 -0.04 * 0.02
## (0.04) (0.04) (0.02) (0.04)
## yr1984 0.09 * -0.04 -0.04 0.04
## (0.04) (0.04) (0.02) (0.05)
## yr1985 0.11 * -0.05 -0.04 * 0.06
## (0.04) (0.03) (0.02) (0.06)
## yr1986 0.14 ** -0.03 -0.03 0.09
## (0.05) (0.02) (0.02) (0.07)
## yr1987 0.17 *** 0.13
## (0.05) (0.08)
## ---------------------------------------------------------------
## R^2 0.19 0.01 0.18 0.18
## Adj. R^2 0.19 0.00 0.06 0.18
## Num. obs. 4360 3815 4360 4360
## s_idios 0.35
## s_id 0.32
## ===============================================================
## *** p < 0.001; ** p < 0.01; * p < 0.05# Hausman test of RE versus FE; H0: RE is consistent
phtest(reg.fe,reg.re)##
## Hausman Test
##
## data: lwage ~ exper + I(exper^2) + married + union + yr
## chisq = 31.707, df = 10, p-value = 0.000448
## alternative hypothesis: one model is inconsistentHere, the random effects model is not consistent.