data management
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
dta2 <- dta2 %>%
mutate(School = factor(School),
Gender = factor(Gender,
levels=0:1,
labels=c("Male","Female")),
US = factor(US),
TVtime = factor(TVtime,
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(dta2)
## School Gender US TVtime 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
descriptive
with(dta2, ftable(US, Gender, Enjoy))
## Enjoy 1 2 3 4
## US Gender
## 0 Male 41 19 58 84
## Female 9 13 31 103
## 1 Male 332 230 589 943
## Female 155 186 487 1334
with(dta2, ftable(US, Economic, Enjoy))
## Enjoy 1 2 3 4
## US Economic
## 0 >50% 30 13 50 112
## 26-50 10 8 22 30
## 11-25 5 4 3 15
## <10% 5 7 14 30
## 1 >50% 214 188 439 891
## 26-50 114 95 246 541
## 11-25 60 36 125 293
## <10% 99 97 266 552
with(dta2, ftable(US, Economic, TVtime, Enjoy))
## Enjoy 1 2 3 4
## US Economic TVtime
## 0 >50% no time 18 7 17 34
## up to 1 hour 1 0 10 11
## 1-3 hours 5 2 14 28
## 3-5 hours 5 4 8 32
## 5 hours or more 1 0 1 7
## 26-50 no time 8 3 6 8
## up to 1 hour 0 2 2 7
## 1-3 hours 2 1 10 10
## 3-5 hours 0 2 4 3
## 5 hours or more 0 0 0 2
## 11-25 no time 4 1 0 4
## up to 1 hour 1 0 1 2
## 1-3 hours 0 2 2 4
## 3-5 hours 0 1 0 4
## 5 hours or more 0 0 0 1
## <10% no time 2 2 3 5
## up to 1 hour 0 1 4 3
## 1-3 hours 2 3 3 12
## 3-5 hours 1 1 4 9
## 5 hours or more 0 0 0 1
## 1 >50% no time 119 73 140 268
## up to 1 hour 26 39 85 126
## 1-3 hours 31 45 135 239
## 3-5 hours 30 28 65 224
## 5 hours or more 8 3 14 34
## 26-50 no time 54 29 71 125
## up to 1 hour 19 24 39 83
## 1-3 hours 24 23 83 162
## 3-5 hours 15 18 51 145
## 5 hours or more 2 1 2 26
## 11-25 no time 24 10 27 55
## up to 1 hour 11 7 27 41
## 1-3 hours 13 10 47 130
## 3-5 hours 8 8 21 59
## 5 hours or more 4 1 3 8
## <10% no time 44 23 61 73
## up to 1 hour 15 21 59 81
## 1-3 hours 20 30 80 221
## 3-5 hours 14 20 61 147
## 5 hours or more 6 3 5 30
barplot(with(dta2, prop.table(table(Enjoy, US), m=2)), beside=T, legend=T, args.legend=list(x="topleft", cex=.5))
model
m0 <- ordinal::clmm(Enjoy ~ Gender + US + Reading + (1 | School) + (1 | TVtime) + (1 | Economic), data = dta2)
## Warning: Using formula(x) is deprecated when x is a character vector of length > 1.
## Consider formula(paste(x, collapse = " ")) instead.
## Warning: (2) Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?
## In addition: Absolute and relative convergence criteria were met
## Cumulative Link Mixed Model fitted with the Laplace approximation
##
## formula: Enjoy ~ Gender + US + Reading + (1 | School) + (1 | TVtime) +
## (1 | Economic)
## data: dta2
##
## link threshold nobs logLik AIC niter max.grad cond.H
## logit flexible 4614 -5095.06 10208.13 599(1850) 2.82e-01 4.9e+06
##
## Random effects:
## Groups Name Variance Std.Dev.
## School (Intercept) 0.166974 0.40863
## TVtime (Intercept) 0.093257 0.30538
## Economic (Intercept) 0.003375 0.05809
## Number of groups: School 180, TVtime 5, Economic 4
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## GenderFemale 0.6984740 0.0593244 11.774 <2e-16 ***
## US1 -0.1629514 0.1117947 -1.458 0.145
## Reading 0.0058490 0.0004666 12.535 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Threshold coefficients:
## Estimate Std. Error z value
## 1|2 1.0564 0.3011 3.509
## 2|3 1.8456 0.3015 6.122
## 3|4 3.1448 0.3041 10.343