wages <- read.table("https://stats.idre.ucla.edu/stat/r/examples/alda/data/wages_pp.txt", header=T, sep=",")
head(wages)
## id lnw exper ged postexp black hispanic hgc hgc.9 uerate ue.7
## 1 31 1.491 0.015 1 0.015 0 1 8 -1 3.215 -3.785
## 2 31 1.433 0.715 1 0.715 0 1 8 -1 3.215 -3.785
## 3 31 1.469 1.734 1 1.734 0 1 8 -1 3.215 -3.785
## 4 31 1.749 2.773 1 2.773 0 1 8 -1 3.295 -3.705
## 5 31 1.931 3.927 1 3.927 0 1 8 -1 2.895 -4.105
## 6 31 1.709 4.946 1 4.946 0 1 8 -1 2.495 -4.505
## ue.centert1 ue.mean ue.person.cen ue1
## 1 0.00 3.215 0.00 3.215
## 2 0.00 3.215 0.00 3.215
## 3 0.00 3.215 0.00 3.215
## 4 0.08 3.215 0.08 3.215
## 5 -0.32 3.215 -0.32 3.215
## 6 -0.72 3.215 -0.72 3.215
#install.packages("lme4")
library(lme4)
## Loading required package: Matrix
#install.packages("sjPlot")
library(sjPlot)
names(wages)
## [1] "id" "lnw" "exper" "ged"
## [5] "postexp" "black" "hispanic" "hgc"
## [9] "hgc.9" "uerate" "ue.7" "ue.centert1"
## [13] "ue.mean" "ue.person.cen" "ue1"
m0 <- lmer(data = wages, lnw ~ 1 + (1 | id))
summary(m0)
## Linear mixed model fit by REML ['lmerMod']
## Formula: lnw ~ 1 + (1 | id)
## Data: wages
##
## REML criterion at convergence: 5870.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.9641 -0.5806 -0.0419 0.5057 6.2096
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 0.06594 0.2568
## Residual 0.11804 0.3436
## Number of obs: 6402, groups: id, 888
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.882846 0.009838 191.4
exp(1.882846)# $6.57/hour in general
## [1] 6.572183
m1 <- lmer(data = wages, lnw ~ 1 + exper + (1 + exper) | id)
summary(m1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: lnw ~ 1 + exper + (1 + exper) | id
## Data: wages
##
## REML criterion at convergence: 5222.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.2634 -0.5066 -0.0327 0.4387 7.0537
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## id (Intercept) 0.066016 0.25694
## exper 0.004116 0.06415 -0.49
## Residual 0.094121 0.30679
## Number of obs: 6402, groups: id, 888
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.824025 0.009044 201.7
exp(1.715609)# $5.56/hour at the start of work
## [1] 5.560061
m2 <- lmer(data = wages, lnw ~ 1 + black + exper + (1 + exper | id))
summary(m2)
## Linear mixed model fit by REML ['lmerMod']
## Formula: lnw ~ 1 + black + exper + (1 + exper | id)
## Data: wages
##
## REML criterion at convergence: 4943.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.2067 -0.5170 -0.0322 0.4379 7.0332
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## id (Intercept) 0.054802 0.23410
## exper 0.001731 0.04161 -0.31
## Residual 0.095100 0.30838
## Number of obs: 6402, groups: id, 888
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.723253 0.012227 140.94
## black -0.027073 0.020027 -1.35
## exper 0.045588 0.002344 19.45
##
## Correlation of Fixed Effects:
## (Intr) black
## black -0.465
## exper -0.519 0.036
exp(1.723253)# $5.60/hour for non-blacks at point 1;
## [1] 5.602725
exp(1.723253-0.027073)# $5.45/hour for blacks at point 1.
## [1] 5.453077
sjp.lmer(m2, type = "fe.cor")
## Computing correlation using pearson-method with listwise-deletion...
m3 <- lmer(data = wages, lnw ~ 1 + black + exper + exper:black + (1 + exper | id))
summary(m3)
## Linear mixed model fit by REML ['lmerMod']
## Formula: lnw ~ 1 + black + exper + exper:black + (1 + exper | id)
## Data: wages
##
## REML criterion at convergence: 4941.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.1995 -0.5247 -0.0317 0.4441 7.0364
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## id (Intercept) 0.054513 0.23348
## exper 0.001678 0.04096 -0.30
## Residual 0.095115 0.30841
## Number of obs: 6402, groups: id, 888
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.712107 0.012684 134.98
## black 0.017350 0.024260 0.72
## exper 0.049705 0.002649 18.76
## black:exper -0.017993 0.005535 -3.25
##
## Correlation of Fixed Effects:
## (Intr) black exper
## black -0.523
## exper -0.567 0.296
## black:exper 0.271 -0.565 -0.479
exp(1.712107)# $5.54/hr for non-blacks at point 1
## [1] 5.540623
exp(1.71210+0.017350)# $5.64/hr for blacks at point 1
## [1] 5.637552
exp(1.712107+0.049705)# $5.82/hr for non-blacks at point 2
## [1] 5.822979
exp(1.712107+0.017350+0.049705
-0.017993)# $5.82/hr for blacks at point 2
## [1] 5.819236
sjp.int(m3, type = "eff", p.kr = F, show.ci = T) #cross-level interaction
sjp.lmer(m3, type = "fe")
## Computing p-values via Kenward-Roger approximation. Use `p.kr = FALSE` if computation takes too long.
m4 <- lmer(data = wages, lnw ~ 1 + black + exper:black + hgc.9 + exper + (1 + exper | id))
summary(m4)
## Linear mixed model fit by REML ['lmerMod']
## Formula: lnw ~ 1 + black + exper:black + hgc.9 + exper + (1 + exper |
## id)
## Data: wages
##
## REML criterion at convergence: 4915.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.2763 -0.5159 -0.0345 0.4430 7.0188
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## id (Intercept) 0.052112 0.22828
## exper 0.001659 0.04073 -0.31
## Residual 0.095166 0.30849
## Number of obs: 6402, groups: id, 888
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.717321 0.012565 136.68
## black 0.015200 0.023972 0.63
## hgc.9 0.038292 0.006444 5.94
## exper 0.049351 0.002641 18.69
## black:exper -0.018120 0.005515 -3.29
##
## Correlation of Fixed Effects:
## (Intr) black hgc.9 exper
## black -0.522
## hgc.9 0.074 -0.016
## exper -0.576 0.301 -0.026
## black:exper 0.274 -0.573 -0.003 -0.478
sjp.lmer(m4, type = "fe")
## Computing p-values via Kenward-Roger approximation. Use `p.kr = FALSE` if computation takes too long.
sjp.lmer(m4, type = "eff")
## Interaction terms in model have been ignored. Use `sjp.int()` to plot effects of interaction terms.
sjp.lmer(m4, type = "pred", vars = "hgc.9")
m5 <- lmer(data = wages, lnw ~ 1 + black + hgc.9 + exper:black + exper:hgc.9 + exper + (1 + exper | id))
summary(m5)
## Linear mixed model fit by REML ['lmerMod']
## Formula: lnw ~ 1 + black + hgc.9 + exper:black + exper:hgc.9 + exper +
## (1 + exper | id)
## Data: wages
##
## REML criterion at convergence: 4925.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.2696 -0.5172 -0.0332 0.4452 7.0194
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## id (Intercept) 0.052112 0.22828
## exper 0.001655 0.04068 -0.31
## Residual 0.095182 0.30852
## Number of obs: 6402, groups: id, 888
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.717146 0.012567 136.64
## black 0.015394 0.023973 0.64
## hgc.9 0.034932 0.007897 4.42
## exper 0.049345 0.002639 18.70
## black:exper -0.018214 0.005514 -3.30
## hgc.9:exper 0.001274 0.001728 0.74
##
## Correlation of Fixed Effects:
## (Intr) black hgc.9 exper blck:x
## black -0.523
## hgc.9 0.071 -0.020
## exper -0.575 0.301 -0.020
## black:exper 0.275 -0.573 0.011 -0.478
## hgc.9:exper -0.019 0.011 -0.578 -0.002 -0.023
sjp.lmer(m5, type = "pred", vars = "hgc.9")
sjt.lmer(m2,m3, m4, m5)
## Computing p-values via Kenward-Roger approximation. Use `p.kr = FALSE` if computation takes too long.
## Computing p-values via Kenward-Roger approximation. Use `p.kr = FALSE` if computation takes too long.
## Computing p-values via Kenward-Roger approximation. Use `p.kr = FALSE` if computation takes too long.
## Computing p-values via Kenward-Roger approximation. Use `p.kr = FALSE` if computation takes too long.
|
|
|
lnw
|
|
lnw
|
|
lnw
|
|
lnw
|
|
|
|
B
|
CI
|
p
|
|
B
|
CI
|
p
|
|
B
|
CI
|
p
|
|
B
|
CI
|
p
|
|
Fixed Parts
|
|
(Intercept)
|
|
1.72
|
1.70 – 1.75
|
<.001
|
|
1.71
|
1.69 – 1.74
|
<.001
|
|
1.72
|
1.69 – 1.74
|
<.001
|
|
1.72
|
1.69 – 1.74
|
<.001
|
|
black
|
|
-0.03
|
-0.07 – 0.01
|
.177
|
|
0.02
|
-0.03 – 0.06
|
.475
|
|
0.02
|
-0.03 – 0.06
|
.526
|
|
0.02
|
-0.03 – 0.06
|
.521
|
|
exper
|
|
0.05
|
0.04 – 0.05
|
<.001
|
|
0.05
|
0.04 – 0.05
|
<.001
|
|
0.05
|
0.04 – 0.05
|
<.001
|
|
0.05
|
0.04 – 0.05
|
<.001
|
|
black:exper
|
|
|
|
|
|
-0.02
|
-0.03 – -0.01
|
.001
|
|
-0.02
|
-0.03 – -0.01
|
.001
|
|
-0.02
|
-0.03 – -0.01
|
<.001
|
|
hgc.9
|
|
|
|
|
|
|
|
|
|
0.04
|
0.03 – 0.05
|
<.001
|
|
0.03
|
0.02 – 0.05
|
<.001
|
|
hgc.9:exper
|
|
|
|
|
|
|
|
|
|
|
|
|
|
0.00
|
-0.00 – 0.00
|
.461
|
|
Random Parts
|
|
σ2
|
|
0.095
|
|
0.095
|
|
0.095
|
|
0.095
|
|
τ00, id
|
|
0.055
|
|
0.055
|
|
0.052
|
|
0.052
|
|
ρ01
|
|
-0.310
|
|
-0.303
|
|
-0.313
|
|
-0.312
|
|
Nid
|
|
888
|
|
888
|
|
888
|
|
888
|
|
ICCid
|
|
0.366
|
|
0.364
|
|
0.354
|
|
0.354
|
|
Observations
|
|
6402
|
|
6402
|
|
6402
|
|
6402
|
|
R2 / Ω02
|
|
.573 / .558
|
|
.572 / .557
|
|
.571 / .556
|
|
.571 / .556
|
#install.packages("plm")
library(plm)
## Loading required package: Formula
# read data
load("us_prep2.Rdata")
# describe data
summary(usl)
## pidp sex age hiqual
## Min. : 68004087 female:29808 Min. :16.00 Degree :12762
## 1st Qu.:204919877 male :22704 1st Qu.:36.00 GCSE etc :10038
## Median :408090449 Median :49.00 A level etc : 9048
## Mean :385366595 Mean :48.93 No qual : 8118
## 3rd Qu.:545056897 3rd Qu.:62.00 Other higher: 6702
## Max. :818710487 Max. :95.00 Other qual : 5814
## (Other) : 30
## wave single urban
## Min. :1.0 No :35460 rural area:13409
## 1st Qu.:2.0 Yes:17052 urban area:39078
## Median :3.5 missing : 0
## Mean :3.5 NA's : 25
## 3rd Qu.:5.0
## Max. :6.0
##
## vote6 fimngrs
## not at all interested?:12306 Min. : 0.0
## not very :14683 1st Qu.: 745.1
## fairly :18430 Median : 1327.6
## very : 5581 Mean : 1688.8
## NA's : 1512 3rd Qu.: 2192.8
## Max. :10000.0
##
## sclfsato sf12pcs sf12mcs
## mostly satisfied :22283 Min. : 4.48 Min. : 1.61
## somewhat satisfied : 8003 1st Qu.:43.84 1st Qu.:45.44
## completely satisfied : 5406 Median :53.24 Median :52.32
## neither sat nor dissat: 4252 Mean :49.23 Mean :50.24
## somewhat dissatisfied : 3854 3rd Qu.:57.27 3rd Qu.:57.16
## (Other) : 3427 Max. :74.46 Max. :77.11
## NA's : 5287 NA's :5135 NA's :5135
## istrtdaty present high sati
## Min. :2009 Min. :1 Min. :0.0000 Min. :1.000
## 1st Qu.:2010 1st Qu.:1 1st Qu.:0.0000 1st Qu.:5.000
## Median :2012 Median :1 Median :0.0000 Median :6.000
## Mean :2012 Mean :1 Mean :0.3709 Mean :5.208
## 3rd Qu.:2013 3rd Qu.:1 3rd Qu.:1.0000 3rd Qu.:6.000
## Max. :2014 Max. :1 Max. :1.0000 Max. :7.000
## NA's :30 NA's :5287
## logincome
## Min. :4.605
## 1st Qu.:6.739
## Median :7.264
## Mean :7.163
## 3rd Qu.:7.738
## Max. :9.220
##
fe <- plm(data = usl, logincome ~ age + single + urban, index = "pidp", model = "within")
## This series is constant and has been removed: present
summary(fe)
## Oneway (individual) effect Within Model
##
## Call:
## plm(formula = logincome ~ age + single + urban, data = usl, model = "within",
## index = "pidp")
##
## Unbalanced Panel: n=8752, T=3-6, N=52487
##
## Residuals :
## Min. 1st Qu. Median 3rd Qu. Max.
## -3.8400 -0.1310 0.0121 0.1760 3.1500
##
## Coefficients :
## Estimate Std. Error t-value Pr(>|t|)
## singleYes 0.145067 0.015892 9.1284 <2e-16 ***
## urbanurban area 0.030316 0.025905 1.1703 0.2419
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 11607
## Residual Sum of Squares: 11584
## R-Squared: 0.0019292
## Adj. R-Squared: -0.19783
## F-statistic: 42.2664 on 2 and 43733 DF, p-value: < 2.22e-16
re <- plm(data = usl, logincome ~ age + single + urban, index = "pidp", model = "random")
## This series is constant and has been removed: present
summary(re)
## Oneway (individual) effect Random Effect Model
## (Swamy-Arora's transformation)
##
## Call:
## plm(formula = logincome ~ age + single + urban, data = usl, model = "random",
## index = "pidp")
##
## Unbalanced Panel: n=8752, T=3-6, N=52487
##
## Effects:
## var std.dev share
## idiosyncratic 0.2649 0.5147 0.34
## individual 0.5147 0.7174 0.66
## theta :
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.6173 0.7189 0.7189 0.7189 0.7189 0.7189
##
## Residuals :
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -3.4800 -0.1700 0.0683 0.0000 0.2680 2.7300
##
## Coefficients :
## Estimate Std. Error t-value Pr(>|t|)
## (Intercept) 7.24237333 0.02850305 254.0912 < 2.2e-16 ***
## age -0.00172702 0.00048076 -3.5923 0.0003281 ***
## singleYes 0.01236114 0.01190749 1.0381 0.2992294
## urbanurban area 0.00128380 0.01526287 0.0841 0.9329672
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 13954
## Residual Sum of Squares: 13944
## R-Squared: 0.00073846
## Adj. R-Squared: 0.00068134
## F-statistic: 12.8972 on 3 and 52483 DF, p-value: 2.0311e-08
phtest(fe, re)
##
## Hausman Test
##
## data: logincome ~ age + single + urban
## chisq = 159.81, df = 2, p-value < 2.2e-16
## alternative hypothesis: one model is inconsistent
#If Hausman test p-value is significant, retain the FE model.