library(car)
library(merTools)
## Warning: package 'merTools' was built under R version 3.2.5
## Loading required package: arm
## Warning: package 'arm' was built under R version 3.2.5
## Loading required package: MASS
## Loading required package: Matrix
## Loading required package: lme4
##
## arm (Version 1.9-3, built: 2016-11-21)
## Working directory is C:/Users/s-das/Syncplicity Folders/MY_Projects/My_Papers/TRB 2018/17-01 Teen Drivers
##
## Attaching package: 'arm'
## The following object is masked from 'package:car':
##
## logit
## Loading required package: dplyr
## Warning: package 'dplyr' was built under R version 3.2.5
##
## Attaching package: 'dplyr'
## The following object is masked from 'package:MASS':
##
## select
## The following object is masked from 'package:car':
##
## recode
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(lme4)
setwd("C:/Users/s-das/Syncplicity Folders/MY_Projects/My_Papers/TRB 2018/17-01 Teen Drivers/Data")
dat1 <- read.csv("Tn_Final_Dataset.csv")
str(dat1)
## 'data.frame': 109267 obs. of 24 variables:
## $ X : int 1 2 3 4 5 6 7 8 9 10 ...
## $ School : int 1 1 1 2 2 2 2 2 2 2 ...
## $ Year : int 2008 2008 2009 2009 2009 2009 2009 2009 2009 2009 ...
## $ Q03_ : Factor w/ 3 levels "","Female","Male": 3 3 2 2 2 2 1 3 3 2 ...
## $ Q05_ : int NA 12 12 12 12 12 12 11 12 10 ...
## $ Q07_ : Factor w/ 5 levels "Blank","learner",..: 1 2 2 2 3 3 4 3 4 2 ...
## $ Q07a_ : int 1 3 3 3 2 2 4 2 4 3 ...
## $ Q09_ : Factor w/ 3 levels "Blank","No","Yes": 1 1 2 1 2 2 3 2 2 2 ...
## $ Q10_ : Factor w/ 3 levels "Blank","No","Yes": 1 2 2 2 2 2 2 3 2 3 ...
## $ Q11_ : Factor w/ 3 levels "Blank","No","Yes": 1 2 2 2 2 2 3 2 3 2 ...
## $ Q12a_ : Factor w/ 4 levels "A lot","Blank",..: 2 3 1 3 3 3 1 3 4 1 ...
## $ Q12b_ : Factor w/ 4 levels "A lot","Blank",..: 2 3 1 4 3 3 1 3 4 1 ...
## $ Q12c_ : Factor w/ 4 levels "A lot","Blank",..: 2 3 4 3 3 3 3 3 3 3 ...
## $ Q12d_ : Factor w/ 4 levels "A lot","Blank",..: 2 3 1 3 3 3 4 3 3 3 ...
## $ Q12e_ : Factor w/ 4 levels "A lot","Blank",..: 2 4 1 4 3 3 4 4 3 3 ...
## $ Q12f_ : Factor w/ 4 levels "A lot","Blank",..: 2 4 1 4 3 3 4 4 3 3 ...
## $ Q12g_ : Factor w/ 4 levels "A lot","Blank",..: 2 1 1 3 3 3 1 4 4 4 ...
## $ Q12h_ : Factor w/ 4 levels "A lot","Blank",..: 2 3 4 3 4 3 1 3 3 3 ...
## $ Q12i_ : Factor w/ 4 levels "A lot","Blank",..: 2 3 4 3 3 3 4 3 3 3 ...
## $ Q12j_ : Factor w/ 4 levels "A lot","Blank",..: 2 3 4 3 3 3 3 3 3 3 ...
## $ Q12k_ : Factor w/ 4 levels "A lot","Blank",..: 2 4 1 3 3 3 1 4 3 4 ...
## $ Q12l_ : Factor w/ 4 levels "A lot","Blank",..: 2 4 1 4 1 3 1 1 3 4 ...
## $ Q12m_ : Factor w/ 4 levels "A lot","Blank",..: 2 4 1 3 3 3 1 3 3 1 ...
## $ schoolState: Factor w/ 11 levels "CA","CO","CT",..: 11 11 11 11 11 11 11 11 11 11 ...
head(dat1,2)
## X School Year Q03_ Q05_ Q07_ Q07a_ Q09_ Q10_ Q11_ Q12a_ Q12b_
## 1 1 1 2008 Male NA Blank 1 Blank Blank Blank Blank Blank
## 2 2 1 2008 Male 12 learner 3 Blank No No Never Never
## Q12c_ Q12d_ Q12e_ Q12f_ Q12g_ Q12h_ Q12i_ Q12j_ Q12k_ Q12l_ Q12m_
## 1 Blank Blank Blank Blank Blank Blank Blank Blank Blank Blank Blank
## 2 Never Never Some Some A lot Never Never Never Some Some Some
## schoolState
## 1 TX
## 2 TX
dim(dat1)
## [1] 109267 24
dat1$Q05_ <- as.factor(dat1$Q05_)
model1 <- lmer(Q07a_ ~ Q03_ + Q05_ + Q09_ + Q10_ + Q11_ + Q12a_+ Q12b_+ Q12c_+ Q12d_+
Q12e_+ Q12f_+ Q12g_+ Q12f_+ Q12g_+Q12h_+Q12i_+ Q12j_+ Q12k_+ Q12l_+(1|School) + (1|Year)+
(1|schoolState), data=dat1)
model1
## Linear mixed model fit by REML ['lmerMod']
## Formula:
## Q07a_ ~ Q03_ + Q05_ + Q09_ + Q10_ + Q11_ + Q12a_ + Q12b_ + Q12c_ +
## Q12d_ + Q12e_ + Q12f_ + Q12g_ + Q12f_ + Q12g_ + Q12h_ + Q12i_ +
## Q12j_ + Q12k_ + Q12l_ + (1 | School) + (1 | Year) + (1 |
## schoolState)
## Data: dat1
## REML criterion at convergence: 265416.2
## Random effects:
## Groups Name Std.Dev.
## School (Intercept) 0.2714
## Year (Intercept) 0.1632
## schoolState (Intercept) 0.1175
## Residual 0.8205
## Number of obs: 108201, groups: School, 281; Year, 11; schoolState, 11
## Fixed Effects:
## (Intercept) Q03_Female Q03_Male Q05_10 Q05_11
## 2.318370 0.131448 0.164231 0.250889 0.527421
## Q05_12 Q09_No Q09_Yes Q10_No Q10_Yes
## 0.742437 0.039233 0.815000 0.121851 0.083850
## Q11_No Q11_Yes Q12a_Blank Q12a_Never Q12a_Some
## 0.069398 0.219649 -0.604259 -0.505271 -0.155797
## Q12b_Blank Q12b_Never Q12b_Some Q12c_Blank Q12c_Never
## -0.219988 -0.199648 -0.074172 0.074474 0.007754
## Q12c_Some Q12d_Blank Q12d_Never Q12d_Some Q12e_Blank
## 0.167594 -0.002932 0.139468 0.031851 0.041352
## Q12e_Never Q12e_Some Q12f_Blank Q12f_Never Q12f_Some
## 0.124003 0.078373 0.025031 0.142397 0.080009
## Q12g_Blank Q12g_Never Q12g_Some Q12h_Blank Q12h_Never
## -0.316447 -0.406454 -0.147701 0.085336 0.189231
## Q12h_Some Q12i_Blank Q12i_Never Q12i_Some Q12j_Blank
## 0.090965 0.201099 0.069290 0.106950 0.074082
## Q12j_Never Q12j_Some Q12k_Blank Q12k_Never Q12k_Some
## 0.207269 0.077243 -0.249488 -0.315113 -0.187781
## Q12l_Blank Q12l_Never Q12l_Some
## -0.071798 0.046030 0.045288
summary(model1)
## Linear mixed model fit by REML ['lmerMod']
## Formula:
## Q07a_ ~ Q03_ + Q05_ + Q09_ + Q10_ + Q11_ + Q12a_ + Q12b_ + Q12c_ +
## Q12d_ + Q12e_ + Q12f_ + Q12g_ + Q12f_ + Q12g_ + Q12h_ + Q12i_ +
## Q12j_ + Q12k_ + Q12l_ + (1 | School) + (1 | Year) + (1 |
## schoolState)
## Data: dat1
##
## REML criterion at convergence: 265416.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.3263 -0.5795 -0.0471 0.5674 4.8437
##
## Random effects:
## Groups Name Variance Std.Dev.
## School (Intercept) 0.07368 0.2714
## Year (Intercept) 0.02665 0.1632
## schoolState (Intercept) 0.01381 0.1175
## Residual 0.67316 0.8205
## Number of obs: 108201, groups: School, 281; Year, 11; schoolState, 11
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.318370 0.087255 26.57
## Q03_Female 0.131448 0.040229 3.27
## Q03_Male 0.164231 0.040232 4.08
## Q05_10 0.250889 0.007506 33.42
## Q05_11 0.527421 0.007801 67.61
## Q05_12 0.742437 0.008507 87.27
## Q09_No 0.039233 0.011568 3.39
## Q09_Yes 0.815000 0.012999 62.70
## Q10_No 0.121851 0.018414 6.62
## Q10_Yes 0.083850 0.018349 4.57
## Q11_No 0.069398 0.017328 4.00
## Q11_Yes 0.219649 0.019803 11.09
## Q12a_Blank -0.604259 0.039174 -15.43
## Q12a_Never -0.505271 0.013087 -38.61
## Q12a_Some -0.155797 0.011132 -14.00
## Q12b_Blank -0.219988 0.040976 -5.37
## Q12b_Never -0.199648 0.012481 -16.00
## Q12b_Some -0.074172 0.011029 -6.73
## Q12c_Blank 0.074474 0.036212 2.06
## Q12c_Never 0.007754 0.019621 0.40
## Q12c_Some 0.167594 0.020601 8.14
## Q12d_Blank -0.002932 0.036553 -0.08
## Q12d_Never 0.139468 0.012813 10.89
## Q12d_Some 0.031851 0.013073 2.44
## Q12e_Blank 0.041352 0.024599 1.68
## Q12e_Never 0.124003 0.009252 13.40
## Q12e_Some 0.078373 0.008857 8.85
## Q12f_Blank 0.025031 0.029598 0.85
## Q12f_Never 0.142397 0.010216 13.94
## Q12f_Some 0.080009 0.009901 8.08
## Q12g_Blank -0.316447 0.030170 -10.49
## Q12g_Never -0.406454 0.010629 -38.24
## Q12g_Some -0.147701 0.009630 -15.34
## Q12h_Blank 0.085336 0.033263 2.57
## Q12h_Never 0.189231 0.014454 13.09
## Q12h_Some 0.090965 0.015177 5.99
## Q12i_Blank 0.201099 0.032735 6.14
## Q12i_Never 0.069290 0.018708 3.70
## Q12i_Some 0.106950 0.018846 5.67
## Q12j_Blank 0.074082 0.031149 2.38
## Q12j_Never 0.207269 0.017145 12.09
## Q12j_Some 0.077243 0.018184 4.25
## Q12k_Blank -0.249488 0.030614 -8.15
## Q12k_Never -0.315113 0.009892 -31.85
## Q12k_Some -0.187781 0.009251 -20.30
## Q12l_Blank -0.071798 0.026256 -2.73
## Q12l_Never 0.046030 0.008782 5.24
## Q12l_Some 0.045288 0.008117 5.58
##
## Correlation matrix not shown by default, as p = 48 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need it
##draw(m1, type = "random")
#3(p1 <- plotFEsim(FEsim(m1)))
##(p1 <- plotREsim(REsim(m1)))
RMSE.merMod(model1)
## [1] 0.8193262
fastdisp(model1)
## lmer(formula = Q07a_ ~ Q03_ + Q05_ + Q09_ + Q10_ + Q11_ + Q12a_ +
## Q12b_ + Q12c_ + Q12d_ + Q12e_ + Q12f_ + Q12g_ + Q12f_ + Q12g_ +
## Q12h_ + Q12i_ + Q12j_ + Q12k_ + Q12l_ + (1 | School) + (1 |
## Year) + (1 | schoolState), data = dat1)
## coef.est coef.se
## (Intercept) 2.32 0.09
## Q03_Female 0.13 0.04
## Q03_Male 0.16 0.04
## Q05_10 0.25 0.01
## Q05_11 0.53 0.01
## Q05_12 0.74 0.01
## Q09_No 0.04 0.01
## Q09_Yes 0.82 0.01
## Q10_No 0.12 0.02
## Q10_Yes 0.08 0.02
## Q11_No 0.07 0.02
## Q11_Yes 0.22 0.02
## Q12a_Blank -0.60 0.04
## Q12a_Never -0.51 0.01
## Q12a_Some -0.16 0.01
## Q12b_Blank -0.22 0.04
## Q12b_Never -0.20 0.01
## Q12b_Some -0.07 0.01
## Q12c_Blank 0.07 0.04
## Q12c_Never 0.01 0.02
## Q12c_Some 0.17 0.02
## Q12d_Blank 0.00 0.04
## Q12d_Never 0.14 0.01
## Q12d_Some 0.03 0.01
## Q12e_Blank 0.04 0.02
## Q12e_Never 0.12 0.01
## Q12e_Some 0.08 0.01
## Q12f_Blank 0.03 0.03
## Q12f_Never 0.14 0.01
## Q12f_Some 0.08 0.01
## Q12g_Blank -0.32 0.03
## Q12g_Never -0.41 0.01
## Q12g_Some -0.15 0.01
## Q12h_Blank 0.09 0.03
## Q12h_Never 0.19 0.01
## Q12h_Some 0.09 0.02
## Q12i_Blank 0.20 0.03
## Q12i_Never 0.07 0.02
## Q12i_Some 0.11 0.02
## Q12j_Blank 0.07 0.03
## Q12j_Never 0.21 0.02
## Q12j_Some 0.08 0.02
## Q12k_Blank -0.25 0.03
## Q12k_Never -0.32 0.01
## Q12k_Some -0.19 0.01
## Q12l_Blank -0.07 0.03
## Q12l_Never 0.05 0.01
## Q12l_Some 0.05 0.01
##
## Error terms:
## Groups Name Std.Dev.
## School (Intercept) 0.27
## Year (Intercept) 0.16
## schoolState (Intercept) 0.12
## Residual 0.82
## ---
## number of obs: 108201, groups: School, 281; Year, 11; schoolState, 11
## AIC = 265520
feexd <- FEsim(model1, 1000)
cbind(feexd[,1] , round(feexd[, 2:4], 3))
## feexd[, 1] mean median sd
## 1 (Intercept) 2.324 2.322 0.088
## 2 Q03_Female 0.131 0.130 0.039
## 3 Q03_Male 0.164 0.163 0.039
## 4 Q05_10 0.251 0.251 0.007
## 5 Q05_11 0.527 0.527 0.008
## 6 Q05_12 0.742 0.742 0.009
## 7 Q09_No 0.039 0.040 0.011
## 8 Q09_Yes 0.815 0.815 0.013
## 9 Q10_No 0.122 0.122 0.019
## 10 Q10_Yes 0.084 0.084 0.019
## 11 Q11_No 0.069 0.069 0.018
## 12 Q11_Yes 0.219 0.219 0.020
## 13 Q12a_Blank -0.603 -0.602 0.041
## 14 Q12a_Never -0.504 -0.504 0.013
## 15 Q12a_Some -0.155 -0.154 0.011
## 16 Q12b_Blank -0.222 -0.221 0.043
## 17 Q12b_Never -0.200 -0.200 0.013
## 18 Q12b_Some -0.075 -0.075 0.011
## 19 Q12c_Blank 0.074 0.074 0.036
## 20 Q12c_Never 0.008 0.007 0.019
## 21 Q12c_Some 0.167 0.167 0.021
## 22 Q12d_Blank -0.002 -0.002 0.037
## 23 Q12d_Never 0.139 0.139 0.013
## 24 Q12d_Some 0.031 0.031 0.013
## 25 Q12e_Blank 0.041 0.041 0.024
## 26 Q12e_Never 0.125 0.125 0.009
## 27 Q12e_Some 0.079 0.078 0.008
## 28 Q12f_Blank 0.026 0.027 0.031
## 29 Q12f_Never 0.143 0.143 0.010
## 30 Q12f_Some 0.081 0.081 0.010
## 31 Q12g_Blank -0.317 -0.318 0.030
## 32 Q12g_Never -0.407 -0.407 0.010
## 33 Q12g_Some -0.148 -0.149 0.010
## 34 Q12h_Blank 0.085 0.086 0.032
## 35 Q12h_Never 0.190 0.189 0.014
## 36 Q12h_Some 0.092 0.092 0.015
## 37 Q12i_Blank 0.199 0.199 0.032
## 38 Q12i_Never 0.068 0.068 0.019
## 39 Q12i_Some 0.106 0.106 0.019
## 40 Q12j_Blank 0.075 0.075 0.030
## 41 Q12j_Never 0.207 0.207 0.017
## 42 Q12j_Some 0.077 0.077 0.018
## 43 Q12k_Blank -0.250 -0.249 0.031
## 44 Q12k_Never -0.315 -0.314 0.010
## 45 Q12k_Some -0.188 -0.188 0.009
## 46 Q12l_Blank -0.072 -0.071 0.027
## 47 Q12l_Never 0.046 0.045 0.009
## 48 Q12l_Some 0.045 0.045 0.008
library(ggplot2)
ggplot(feexd[feexd$term!= "(Intercept)", ]) +
aes(x = term, ymin = median - 1.96 * sd,
ymax = median + 1.96 * sd, y = median) +
geom_pointrange() +
geom_hline(yintercept = 0, size = I(1.1), color = I("red")) +
coord_flip() +
theme_bw() + labs(title = "Coefficient Plot",
x = "Median Effect Estimate", y = "License Rating")
plotFEsim(feexd) +
theme_bw() + labs(title = "Coefficient Plot",
x = "Median Effect Estimate", y = "License Rating")
reexd <- REsim(model1)
head(reexd)
## groupFctr groupID term mean median sd
## 1 School 1 (Intercept) -0.07822608 -0.07415088 0.07244622
## 2 School 2 (Intercept) -0.06108175 -0.06350232 0.05533199
## 3 School 3 (Intercept) -0.11390221 -0.09601669 0.16442154
## 4 School 4 (Intercept) -0.21920848 -0.21423055 0.20909131
## 5 School 5 (Intercept) -0.14074314 -0.13964380 0.05556122
## 6 School 6 (Intercept) -0.12946944 -0.12636285 0.14645082
table(reexd$term)
##
## (Intercept)
## 303
table(reexd$groupFctr)
##
## School schoolState Year
## 281 11 11
plot1 <- plotREsim(reexd)
plot1
lattice::dotplot(ranef(model1, condVar=TRUE))
## $School
##
## $Year
##
## $schoolState
