#Nested Designs - Chicano #input data
dta <- read.csv('C:/Users/Ching-Fang Wu/Documents/data/Chicano.csv', stringsAsFactors = TRUE)
#exam structure
str(dta)
## 'data.frame': 24 obs. of 4 variables:
## $ Pupil: Factor w/ 24 levels "P01","P02","P03",..: 1 2 3 4 5 6 7 8 9 10 ...
## $ Class: Factor w/ 6 levels "C1","C2","C3",..: 1 1 1 1 2 2 2 2 3 3 ...
## $ Trt : Factor w/ 2 levels "C","T": 2 2 2 2 2 2 2 2 2 2 ...
## $ Score: int 5 2 10 11 10 15 11 8 7 12 ...
pacman::p_load(tidyverse, VCA, lme4, nlme)
#plot data as variability-chart
VCA::varPlot(Score ~ Trt/Class/Pupil,
Data=dta,
YLabel=list(text="Score",
side=2,
cex=1),
MeanLine=list(var=c("Trt", "Class"),
col=c("darkred", "salmon"),
lwd=c(1, 2)))
Function varPlot determines the sequence of variables in the model formula and uses this information to construct the variability chart.(https://rdrr.io/cran/VCA/man/varPlot.html)
summary(m1 <- aov(Score ~ Trt + Error(Class), data=dta))
##
## Error: Class
## Df Sum Sq Mean Sq F value Pr(>F)
## Trt 1 216 216 9.818 0.0351 *
## Residuals 4 88 22
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 18 198 11
faraway::sumary(m2 <- lmer(Score ~ Trt + (1 | Class), data=dta))
## Fixed Effects:
## coef.est coef.se
## (Intercept) 4.00 1.35
## TrtT 6.00 1.91
##
## Random Effects:
## Groups Name Std.Dev.
## Class (Intercept) 1.66
## Residual 3.32
## ---
## number of obs: 24, groups: Class, 6
## AIC = 130.9, DIC = 131.8
## deviance = 127.4
confint(m2, method="boot")
## Computing bootstrap confidence intervals ...
##
## 130 message(s): boundary (singular) fit: see ?isSingular
## 2.5 % 97.5 %
## .sig01 0.000000 3.624950
## .sigma 2.223250 4.385830
## (Intercept) 1.213864 6.582716
## TrtT 2.178665 9.850228