#Title: 0928 Homework3
#Name: Szu-Yu Chen
#Date: 4 Oct,2020
pacman::p_load(car, tidyverse, lme4, GGally)
dta_EEG <- read.table("C:/Users/ASUS/Desktop/data/thetaEEG.txt", header=T)
str(dta_EEG)
## 'data.frame': 19 obs. of 10 variables:
## $ ID : int 1 2 3 4 5 6 7 8 9 10 ...
## $ Ch3 : num -3.54 5.72 0.52 0 2.07 1.67 9.13 -0.43 -0.56 1.28 ...
## $ Ch4 : num -3.11 5.07 -0.18 0.74 0.76 ...
## $ Ch5 : num -0.24 6.87 0.9 1.1 3.51 2.77 3.44 -0.31 -1.22 1.89 ...
## $ Ch6 : num 0.42 5.96 0.6 0.13 0.6 4.55 4.8 -0.61 0.67 1.77 ...
## $ Ch7 : num -0.49 8.2 1.27 0.19 3.71 1.8 0.48 -1.04 -0.97 1.83 ...
## $ Ch8 : num 2.13 4.87 1.28 0.07 1.86 4.79 1.63 -0.13 -0.98 0.91 ...
## $ Ch17: num -4.15 5.48 -0.95 0.8 1.49 4.51 9.94 -0.61 0 1.4 ...
## $ Ch18: num 2.87 5.57 1.74 0.25 3.11 3.24 1.34 -0.61 -1.22 1.1 ...
## $ Ch19: num 1.34 6.33 0.79 -0.66 1.8 3.99 1.53 -0.43 -0.91 -0.12 ...
head(dta_EEG)
## ID Ch3 Ch4 Ch5 Ch6 Ch7 Ch8 Ch17 Ch18 Ch19
## 1 1 -3.54 -3.11 -0.24 0.42 -0.49 2.13 -4.15 2.87 1.34
## 2 2 5.72 5.07 6.87 5.96 8.20 4.87 5.48 5.57 6.33
## 3 3 0.52 -0.18 0.90 0.60 1.27 1.28 -0.95 1.74 0.79
## 4 4 0.00 0.74 1.10 0.13 0.19 0.07 0.80 0.25 -0.66
## 5 5 2.07 0.76 3.51 0.60 3.71 1.86 1.49 3.11 1.80
## 6 6 1.67 4.18 2.77 4.55 1.80 4.79 4.51 3.24 3.99
ch <- as.factor(colnames(dta_EEG[,-1]))
m0 <- lm(as.matrix(dta_EEG[, -1]) ~ 1)
aov0 <- Anova(m0, idata = data.frame(ch), idesign = ~ ch)
## Note: model has only an intercept; equivalent type-III tests substituted.
summary(aov0)
##
## Type III Repeated Measures MANOVA Tests:
##
## ------------------------------------------
##
## Term: (Intercept)
##
## Response transformation matrix:
## (Intercept)
## Ch3 1
## Ch4 1
## Ch5 1
## Ch6 1
## Ch7 1
## Ch8 1
## Ch17 1
## Ch18 1
## Ch19 1
##
## Sum of squares and products for the hypothesis:
## (Intercept)
## (Intercept) 1761.616
##
## Multivariate Tests: (Intercept)
## Df test stat approx F num Df den Df Pr(>F)
## Pillai 1 0.2035907 4.601445 1 18 0.045839 *
## Wilks 1 0.7964093 4.601445 1 18 0.045839 *
## Hotelling-Lawley 1 0.2556358 4.601445 1 18 0.045839 *
## Roy 1 0.2556358 4.601445 1 18 0.045839 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------
##
## Term: ch
##
## Response transformation matrix:
## ch1 ch2 ch3 ch4 ch5 ch6 ch7 ch8
## Ch3 0 0 0 1 0 0 0 0
## Ch4 0 0 0 0 1 0 0 0
## Ch5 0 0 0 0 0 1 0 0
## Ch6 0 0 0 0 0 0 1 0
## Ch7 0 0 0 0 0 0 0 1
## Ch8 -1 -1 -1 -1 -1 -1 -1 -1
## Ch17 1 0 0 0 0 0 0 0
## Ch18 0 1 0 0 0 0 0 0
## Ch19 0 0 1 0 0 0 0 0
##
## Sum of squares and products for the hypothesis:
## ch1 ch2 ch3 ch4 ch5 ch6
## ch1 5.7475 1.650000e-02 1.485000000 0.506000000 7.70000000 1.930500000
## ch2 0.0165 4.736842e-05 0.004263158 0.001452632 0.02210526 0.005542105
## ch3 1.4850 4.263158e-03 0.383684211 0.130736842 1.98947368 0.498789474
## ch4 0.5060 1.452632e-03 0.130736842 0.044547368 0.67789474 0.169957895
## ch5 7.7000 2.210526e-02 1.989473684 0.677894737 10.31578947 2.586315789
## ch6 1.9305 5.542105e-03 0.498789474 0.169957895 2.58631579 0.648426316
## ch7 3.0635 8.794737e-03 0.791526316 0.269705263 4.10421053 1.028984211
## ch8 -0.0165 -4.736842e-05 -0.004263158 -0.001452632 -0.02210526 -0.005542105
## ch7 ch8
## ch1 3.063500000 -1.650000e-02
## ch2 0.008794737 -4.736842e-05
## ch3 0.791526316 -4.263158e-03
## ch4 0.269705263 -1.452632e-03
## ch5 4.104210526 -2.210526e-02
## ch6 1.028984211 -5.542105e-03
## ch7 1.632889474 -8.794737e-03
## ch8 -0.008794737 4.736842e-05
##
## Multivariate Tests: ch
## Df test stat approx F num Df den Df Pr(>F)
## Pillai 1 0.5243665 1.515882 8 11 0.25588
## Wilks 1 0.4756335 1.515882 8 11 0.25588
## Hotelling-Lawley 1 1.1024593 1.515882 8 11 0.25588
## Roy 1 1.1024593 1.515882 8 11 0.25588
##
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
##
## Sum Sq num Df Error SS den Df F value Pr(>F)
## (Intercept) 195.735 1 765.68 18 4.6014 0.04584 *
## ch 10.702 8 389.31 144 0.4948 0.85840
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## ch 0.00037434 9.7367e-11
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## ch 0.31531 0.6559
##
## HF eps Pr(>F[HF])
## ch 0.3709213 0.6854028
Long <- dta_EEG %>%
gather( key = channel, value = EEG, 2:10) %>%
mutate( channel = factor(channel)) %>%
mutate( ID = factor(ID))
str(Long)
## 'data.frame': 171 obs. of 3 variables:
## $ ID : Factor w/ 19 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
## $ channel: Factor w/ 9 levels "Ch17","Ch18",..: 4 4 4 4 4 4 4 4 4 4 ...
## $ EEG : num -3.54 5.72 0.52 0 2.07 1.67 9.13 -0.43 -0.56 1.28 ...
ot <- theme_set(theme_bw())
#Plots
ggplot(Long, aes(reorder(channel, EEG, median), EEG)) +
geom_boxplot() +
geom_point(alpha = .3) +
geom_hline(yintercept = mean(Long$EEG)) +
labs(x ="Channel effects", y = "EEG")
ggplot(Long, aes(reorder(ID, EEG, mean), EEG, color = channel)) +
geom_point() +
geom_hline(yintercept = mean(Long$EEG), linetype = "dashed") +
coord_flip() +
scale_color_manual(values = c("gray", "forestgreen", "peru", "navy", "blue", "red", "coral", "green", "pink"))+
labs(y = "EEG", x = "ID") +
theme(legend.position = c(.9, .2))
# Residual plot
summary(m1 <- lmer(EEG ~ channel + (1 | ID), data = Long))
## Linear mixed model fit by REML ['lmerMod']
## Formula: EEG ~ channel + (1 | ID)
## Data: Long
##
## REML criterion at convergence: 697
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.6502 -0.3882 -0.0002 0.4215 4.6394
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 4.426 2.104
## Residual 2.704 1.644
## Number of obs: 171, groups: ID, 19
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.4026 0.6126 2.290
## channelCh18 -0.5484 0.5335 -1.028
## channelCh19 -0.4079 0.5335 -0.765
## channelCh3 -0.5016 0.5335 -0.940
## channelCh4 0.1868 0.5335 0.350
## channelCh5 -0.3653 0.5335 -0.685
## channelCh6 -0.2568 0.5335 -0.481
## channelCh7 -0.5516 0.5335 -1.034
## channelCh8 -0.5500 0.5335 -1.031
##
## Correlation of Fixed Effects:
## (Intr) chnC18 chnC19 chnnC3 chnnC4 chnnC5 chnnC6 chnnC7
## channelCh18 -0.435
## channelCh19 -0.435 0.500
## channelCh3 -0.435 0.500 0.500
## channelCh4 -0.435 0.500 0.500 0.500
## channelCh5 -0.435 0.500 0.500 0.500 0.500
## channelCh6 -0.435 0.500 0.500 0.500 0.500 0.500
## channelCh7 -0.435 0.500 0.500 0.500 0.500 0.500 0.500
## channelCh8 -0.435 0.500 0.500 0.500 0.500 0.500 0.500 0.500
plot(m1, form = resid(., type = "pearson") ~ fitted(.) | channel,
abline = 0, lty = 2, layout = c(4,1), pch = 20,
xlab = "Fitted values", ylab = " Residuals")
