1 Data Management

  1. Load & Rename variables
dtaHW2 <- read.table("C:/Users/ASUS/Desktop/data/pirls_us.txt", h=T)
head(dtaHW2)
##   schoolid female bornUS watchTV econDisadv  reading enjoy
## 1        1      0      1       3       <10% 547.0412     1
## 2        1      0      1       4       <10% 545.9297     2
## 3        1      0      1       4       <10% 587.8460     2
## 4        1      1      1       3       <10% 450.8517     2
## 5        1      0      1       1       <10% 525.2201     3
## 6        1      0      1       1       <10% 555.4918     3
names(dtaHW2) <- c("ID","Gender","US", "TV", "Economic", "Reading", "Enjoy")
head(dtaHW2)
##   ID Gender US TV Economic  Reading Enjoy
## 1  1      0  1  3     <10% 547.0412     1
## 2  1      0  1  4     <10% 545.9297     2
## 3  1      0  1  4     <10% 587.8460     2
## 4  1      1  1  3     <10% 450.8517     2
## 5  1      0  1  1     <10% 525.2201     3
## 6  1      0  1  1     <10% 555.4918     3
str(dtaHW2)
## 'data.frame':    4614 obs. of  7 variables:
##  $ ID      : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ Gender  : int  0 0 0 1 0 0 1 0 0 0 ...
##  $ US      : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ TV      : int  3 4 4 3 1 1 3 3 3 3 ...
##  $ Economic: chr  "<10%" "<10%" "<10%" "<10%" ...
##  $ Reading : num  547 546 588 451 525 ...
##  $ Enjoy   : int  1 2 2 2 3 3 3 4 4 4 ...
pacman::p_load(magrittr, vcrpart, ggplot2, dplyr, tidyr, ltm, eRm, lme4,ordinal)
dtaHW2 <- dtaHW2 %>% 
  mutate(ID = factor(ID),
         Gender = factor(Gender, 
                         levels=0:1, 
                         labels=c("Male","Female")),
         US = factor(US),
         TV = factor(TV, levels=1:5,
                         labels=c("no time",
                                  "up to 1 hour", 
                                  "1-3 hours",
                                  "3-5 hours",
                                  "5 hours or more")), 
         Economic = factor(Economic,
                           levels = c(">50%", 
                                      "26-50", 
                                      "11-25", 
                                      "<10%")),
         Enjoy = factor(Enjoy, levels = 1:4))

head(dtaHW2)
##   ID Gender US        TV Economic  Reading Enjoy
## 1  1   Male  1 1-3 hours     <10% 547.0412     1
## 2  1   Male  1 3-5 hours     <10% 545.9297     2
## 3  1   Male  1 3-5 hours     <10% 587.8460     2
## 4  1 Female  1 1-3 hours     <10% 450.8517     2
## 5  1   Male  1   no time     <10% 525.2201     3
## 6  1   Male  1   no time     <10% 555.4918     3

2 GLMM

m0 <- ordinal::clmm(Enjoy ~ Gender + US + TV + (1 | ID) +(Reading| ID)+(Economic|ID), data = dtaHW2)
## Warning: Using formula(x) is deprecated when x is a character vector of length > 1.
##   Consider formula(paste(x, collapse = " ")) instead.
summary(m0)
## Cumulative Link Mixed Model fitted with the Laplace approximation
## 
## formula: Enjoy ~ Gender + US + TV + (1 | ID) + (Reading | ID) + (Economic |  
##     ID)
## data:    dtaHW2
## 
##  link  threshold nobs logLik   AIC      niter       max.grad cond.H 
##  logit flexible  4614 -5126.27 10298.53 3590(11052) 4.53e+02 6.8e+07
## 
## Random effects:
##  Groups Name          Variance  Std.Dev. Corr                 
##  ID     (Intercept)   0.000e+00 0.000000                      
##  ID     (Intercept)   1.459e+01 3.820145                      
##         Reading       4.357e-05 0.006601 -1.000               
##  ID     (Intercept)   1.917e-01 0.437892                      
##         Economic26-50 2.137e-01 0.462301 -0.999               
##         Economic11-25 8.030e-02 0.283368 -0.706  0.720        
##         Economic<10%  3.356e-01 0.579311 -0.846  0.859  0.480 
## Number of groups:  ID 180 
## 
## Coefficients:
##                   Estimate Std. Error z value Pr(>|z|)    
## GenderFemale       0.69340    0.06008  11.541  < 2e-16 ***
## US1               -0.07741    0.11383  -0.680  0.49649    
## TVup to 1 hour     0.27612    0.08917   3.097  0.00196 ** 
## TV1-3 hours        0.62729    0.07871   7.970 1.59e-15 ***
## TV3-5 hours        0.78165    0.08786   8.897  < 2e-16 ***
## TV5 hours or more  1.01146    0.18436   5.486 4.10e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Threshold coefficients:
##     Estimate Std. Error z value
## 1|2  -1.6219     0.1361 -11.919
## 2|3  -0.8247     0.1328  -6.210
## 3|4   0.4835     0.1316   3.674