#' title: "Homework 4 Cognitive task"
#' author: "Szu-Yu Chen"
#' date: "04 Oct 2020"


#載入資料
dta <- read.table("C:/Users/ASUS/Desktop/data/cognitive_task.txt", h=T)

#繪圖
VCA::varPlot(Score ~ Method/Class/ID, 
             Data=dta,
             YLabel=list(text="Score", side=2, cex=1),
             MeanLine=list(var=c("Method", "Class"),
                           col=c("darkred", "salmon"), lwd=c(1, 2)))

## Warning in min(x): min 中沒有無漏失的引數; 回傳 Inf

## Warning in max(x): max 中沒有無漏失的引數；回傳 -Inf

# define a ANOVA model containing  Method and error term of Class (Class as Random effect)
m0 <- aov(Score ~ Method + Error(Class), data=dta)
summary(m0)

## 
## Error: Class
##           Df Sum Sq Mean Sq F value Pr(>F)  
## Method     1  112.5  112.50   6.459  0.044 *
## Residuals  6  104.5   17.42                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Error: Within
##           Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 24   18.5  0.7708

#compute a model containing Method and error terms of Method, Method*Class 
m01 <- aov(Score ~ Method + Error(Method/Class), data=dta)

## Warning in aov(Score ~ Method + Error(Method/Class), data = dta): Error() model
## is singular

summary(m01)

## 
## Error: Method
##        Df Sum Sq Mean Sq
## Method  1  112.5   112.5
## 
## Error: Method:Class
##           Df Sum Sq Mean Sq F value Pr(>F)
## Residuals  6  104.5   17.42               
## 
## Error: Within
##           Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 24   18.5  0.7708

# m1 has a random intercept per Class 
m1 <- lme4::lmer(Score ~ Method + (1 | Class), data=dta)
summary(m1)

## Linear mixed model fit by REML ['lmerMod']
## Formula: Score ~ Method + (1 | Class)
##    Data: dta
## 
## REML criterion at convergence: 101.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.3028 -0.8416 -0.0126  0.3150  2.5753 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Class    (Intercept) 4.1615   2.040   
##  Residual             0.7708   0.878   
## Number of obs: 32, groups:  Class, 8
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)    3.500      1.043   3.355
## MethodI2       3.750      1.475   2.542
## 
## Correlation of Fixed Effects:
##          (Intr)
## MethodI2 -0.707

# m11 has a random intercept per Method*Klass 
m11 <- lme4::lmer(Score ~ Method + (1 | Method:Klass), data=dta)
summary(m11)

## Linear mixed model fit by REML ['lmerMod']
## Formula: Score ~ Method + (1 | Method:Klass)
##    Data: dta
## 
## REML criterion at convergence: 101.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.3028 -0.8416 -0.0126  0.3150  2.5753 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  Method:Klass (Intercept) 4.1615   2.040   
##  Residual                 0.7708   0.878   
## Number of obs: 32, groups:  Method:Klass, 8
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)    3.500      1.043   3.355
## MethodI2       3.750      1.475   2.542
## 
## Correlation of Fixed Effects:
##          (Intr)
## MethodI2 -0.707

# confidence intervals estimated by the bootstrapping method
confint(m1, method="boot")

## Computing bootstrap confidence intervals ...

##                 2.5 %   97.5 %
## .sig01      0.7983641 3.189876
## .sigma      0.6362361 1.133768
## (Intercept) 1.5120559 5.531759
## MethodI2    0.9332159 6.583586

# The end