0.1 input Data

dta <- read.csv('Chicano.csv', stringsAsFactors = TRUE)

0.2 load package

pacman::p_load(tidyverse, VCA, lme4, nlme)

0.3 varPlot by Trt/Class/Pupil

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)))

Variability Chart for Hierarchical
Models.score依照Trt、Class、Pupil進行畫圖
左邊紅色線代表treament組的socre 平均
左邊橘色線代表c1、c2、c3的socre 平均
點代表每一個樣本的分數，右邊相同

結果圖顯示，treament組的socre高於control組 但左邊c1、c2、c3班級本身平均socre就高於右邊c4、c5、c6

r語法理解參考：https://rdrr.io/cran/VCA/man/varPlot.html

0.4 Error(Class) have p value

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

anova檢定實驗與控制組score平均有顯著差異(p=0.0351)

0.5 Caculating fix and random effect (1 | Class)

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

考慮class層級隨機效果，學生平均分數是4分，實驗組相比對照組分數多6分。
class的隨機效果變異量為1.66，看起來class間的隨機效果比個人間來的小(Std.Dev為3.32)

0.6 remove intercept effect

faraway::sumary(m2 <- lmer(Score ~ Trt -1 + (1 | Class), data=dta))
## Fixed Effects:
##      coef.est coef.se
## TrtC  4.00     1.35  
## TrtT 10.00     1.35  
## 
## 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

0.7 95% CI

confint(m2, method="boot")
## Computing bootstrap confidence intervals ...
## 
## 129 message(s): boundary (singular) fit: see ?isSingular
##           2.5 %    97.5 %
## .sig01 0.000000  3.589149
## .sigma 2.174806  4.288415
## TrtC   1.420648  6.557520
## TrtT   7.283613 12.528856

計算信賴區間
TrtT相對TrtC的score的effect信賴區間為:2.14~9.84

0.8 Question

1.在confint(m2, method=“boot”)的結果中.sig01和.sigma代表的是什麼?
2.m2 <- lmer(Score ~ Trt + (1 | Class)中沒有score fix effect的p值，fix effect的p值是否可以直接使用anova做出來的檢定結果?又或者應該用何種語法產出?