In Class Exercise 1
pacman::p_load(tidyverse, lme4)
#
dta <- read.table("family.txt", header = T)
# data structure
str(dta)
## 'data.frame': 12 obs. of 2 variables:
## $ family : int 1 1 1 2 2 2 2 3 3 4 ...
## $ liberalism: int 25 29 34 18 21 21 26 31 35 30 ...
# first 6 lines
head(dta)
## family liberalism
## 1 1 25
## 2 1 29
## 3 1 34
## 4 2 18
## 5 2 21
## 6 2 21
dta$family<- as.factor(dta$family)
# aov with fixed family effects
summary(m0 <- aov(liberalism ~ family, data = dta))
## Df Sum Sq Mean Sq F value Pr(>F)
## family 3 244.33 81.44 7.726 0.00951 **
## Residuals 8 84.33 10.54
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# aov with random family effects
summary(m1 <- aov(liberalism ~ Error(family), data = dta))
##
## Error: family
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 3 244.3 81.44
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 8 84.33 10.54
# regression **without intercept** term
summary(m2 <- lm(liberalism ~ family - 1, data = dta))
##
## Call:
## lm(formula = liberalism ~ family - 1, data = dta)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.333 -1.000 -0.500 1.500 4.667
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## family1 29.333 1.875 15.65 2.77e-07 ***
## family2 21.500 1.623 13.24 1.01e-06 ***
## family3 33.000 2.296 14.37 5.36e-07 ***
## family4 30.667 1.875 16.36 1.96e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.247 on 8 degrees of freedom
## Multiple R-squared: 0.9911, Adjusted R-squared: 0.9867
## F-statistic: 223.6 on 4 and 8 DF, p-value: 3.065e-08
anova(m2)
## Analysis of Variance Table
##
## Response: liberalism
## Df Sum Sq Mean Sq F value Pr(>F)
## family 4 9429.7 2357.42 223.63 3.065e-08 ***
## Residuals 8 84.3 10.54
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# lme4 package
summary(m3 <- lmer(liberalism ~ 1 + (1 | family), data = dta))
## Linear mixed model fit by REML ['lmerMod']
## Formula: liberalism ~ 1 + (1 | family)
## Data: dta
##
## REML criterion at convergence: 65.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.3116 -0.3850 -0.1135 0.5998 1.4772
##
## Random effects:
## Groups Name Variance Std.Dev.
## family (Intercept) 21.92 4.682
## Residual 10.48 3.238
## Number of obs: 12, groups: family, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 28.49 2.53 11.26
# show random effects
ranef(m3)
## $family
## (Intercept)
## 1 0.7271474
## 2 -6.2437983
## 3 3.6394896
## 4 1.8771613
##
## with conditional variances for "family"
# augument data with group mean and ID
dta_fy <- dta %>%
mutate( ID = paste0("S", 11:22), g_ave = mean(liberalism )) %>%
group_by(family) %>%
mutate( grp_ave = mean(liberalism) ) %>%
as.data.frame %>%
mutate( grp_ave_uw = mean(unique(grp_ave)) )
# output from model m3
dta_m3 <- data.frame(ID = paste0("S",11:22), fe = rep(summary(m3)$coef[1], 12),
re = rep(as.numeric(t(as.data.frame(ranef(m3)[1]))), c(3, 4, 2, 3)),
pred = fitted(m3))
# join two data frames
merge(dta_fy, dta_m3, by = "ID")
## ID family liberalism g_ave grp_ave grp_ave_uw fe re
## 1 S11 1 25 27.66667 29.33333 28.625 28.49027 0.7271474
## 2 S12 1 29 27.66667 29.33333 28.625 28.49027 0.7271474
## 3 S13 1 34 27.66667 29.33333 28.625 28.49027 0.7271474
## 4 S14 2 18 27.66667 21.50000 28.625 28.49027 -6.2437983
## 5 S15 2 21 27.66667 21.50000 28.625 28.49027 -6.2437983
## 6 S16 2 21 27.66667 21.50000 28.625 28.49027 -6.2437983
## 7 S17 2 26 27.66667 21.50000 28.625 28.49027 -6.2437983
## 8 S18 3 31 27.66667 33.00000 28.625 28.49027 3.6394896
## 9 S19 3 35 27.66667 33.00000 28.625 28.49027 3.6394896
## 10 S20 4 30 27.66667 30.66667 28.625 28.49027 1.8771613
## 11 S21 4 30 27.66667 30.66667 28.625 28.49027 1.8771613
## 12 S22 4 32 27.66667 30.66667 28.625 28.49027 1.8771613
## pred
## 1 29.21742
## 2 29.21742
## 3 29.21742
## 4 22.24648
## 5 22.24648
## 6 22.24648
## 7 22.24648
## 8 32.12976
## 9 32.12976
## 10 30.36744
## 11 30.36744
## 12 30.36744
# The end
In Class Exercise 2
#讀取資料
dta <- read.csv('Chicano.csv', stringsAsFactors = TRUE)
#載入所需套件
pacman::p_load(tidyverse, VCA, lme4, nlme)
#畫出nest design的變異
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)))
#take class as the random effect with aov function
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
#take class as the random effect with lmer function
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 ...
##
## 138 message(s): boundary (singular) fit: see ?isSingular
## 2.5 % 97.5 %
## .sig01 0.000000 3.648370
## .sigma 2.224671 4.176341
## (Intercept) 1.354633 6.809314
## TrtT 1.711791 10.026795
In Class Exercise 3
#Load package
pacman::p_load(tidyverse, lme4)
# input data
dta <- read.table("finger_ridges.txt", header = T)
# 將Pair指定為因子
dta$Pair <- factor(dta$Pair)
# show first 6 lines
head(dta)
## Pair Twin_1 Twin_2
## 1 1 71 71
## 2 2 79 82
## 3 3 105 99
## 4 4 115 114
## 5 5 76 70
## 6 6 83 82
# ICC
cor(c(dta$Twin_1, dta$Twin_2), c(dta$Twin_2, dta$Twin_1))
## [1] 0.9625626
# turn data from wide format to long format
dtaL <- dta %>%
gather(key = Twin, value = Ridges, Twin_1, Twin_2) %>%
mutate( Twin = as.factor(Twin))
# show mean and sd of ridges by twin pair (ordered)
dtaL %>% group_by(Pair) %>%
summarize( ridge_ave = mean(Ridges), ridge_sd = sd(Ridges)) %>%
arrange( desc(ridge_ave))
## # A tibble: 12 x 3
## Pair ridge_ave ridge_sd
## <fct> <dbl> <dbl>
## 1 4 114. 0.707
## 2 7 114. 0.707
## 3 9 114. 0.707
## 4 3 102 4.24
## 5 10 92.5 2.12
## 6 6 82.5 0.707
## 7 2 80.5 2.12
## 8 11 79 5.66
## 9 12 74 2.83
## 10 5 73 4.24
## 11 1 71 0
## 12 8 50.5 9.19
# dot plot
ggplot(dtaL, aes(reorder(Pair, Ridges, mean), Ridges, color = Pair)) +
geom_point() +
geom_hline(yintercept = mean(dtaL$Ridges), linetype = "dashed") +
labs(y = "Number of finger ridges", x = "Twin ID") +
theme_bw() +
theme(legend.position = "NONE")
# random effect anova for balanced designs
summary(m0 <- aov(Ridges ~ Error(Pair), data = dtaL))
##
## Error: Pair
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 11 8990 817.3
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 12 171.5 14.29
# random effects anova
summary(m1 <- lmer(Ridges ~ 1 + (1 | Pair), data = dtaL))
## Linear mixed model fit by REML ['lmerMod']
## Formula: Ridges ~ 1 + (1 | Pair)
## Data: dtaL
##
## REML criterion at convergence: 174.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.88917 -0.38618 -0.00834 0.37963 1.54959
##
## Random effects:
## Groups Name Variance Std.Dev.
## Pair (Intercept) 401.51 20.04
## Residual 14.29 3.78
## Number of obs: 24, groups: Pair, 12
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 87.208 5.836 14.94
#Extract the modes of the random effects
ranef(m1)
## $Pair
## (Intercept)
## 1 -15.924912
## 2 -6.591031
## 3 14.533018
## 4 26.814441
## 5 -13.959885
## 6 -4.626003
## 7 25.831927
## 8 -36.066447
## 9 25.831927
## 10 5.199136
## 11 -8.064801
## 12 -12.977371
##
## with conditional variances for "Pair"
Exercise 1
#install.packages("PairedData")
library(PairedData)
## Warning: package 'PairedData' was built under R version 4.0.2
## Loading required package: MASS
##
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
##
## select
## Loading required package: gld
## Warning: package 'gld' was built under R version 4.0.2
## Loading required package: mvtnorm
## Loading required package: lattice
##
## Attaching package: 'PairedData'
## The following object is masked from 'package:Matrix':
##
## summary
## The following object is masked from 'package:base':
##
## summary
#load data
data(Anorexia, package="PairedData")
#plot
with(Anorexia, plot(paired(Prior, Post), type="profile")) +
labs(x="Therapy", y="Weight (in lb)") +
theme_bw()
str(Anorexia)
## 'data.frame': 17 obs. of 2 variables:
## $ Prior: num 83.8 83.3 86 82.5 86.7 79.6 76.9 94.2 73.4 80.5 ...
## $ Post : num 95.2 94.3 91.5 91.9 100.3 ...
colMeans(Anorexia)
## Prior Post
## 83.22941 90.49412
# compute variable sd
print(apply(Anorexia, 2, sd), 3)
## Prior Post
## 5.02 8.48
# covariances
cov(Anorexia)
## Prior Post
## Prior 25.16721 22.88268
## Post 22.88268 71.82684
# correlations
print(cor(Anorexia),3)
## Prior Post
## Prior 1.000 0.538
## Post 0.538 1.000
# 2-sample t-test
t.test(Anorexia$Prior, Anorexia$Post)
##
## Welch Two Sample t-test
##
## data: Anorexia$Prior and Anorexia$Post
## t = -3.0414, df = 25.986, p-value = 0.005324
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -12.17472 -2.35469
## sample estimates:
## mean of x mean of y
## 83.22941 90.49412
# paired t-test
t.test(Anorexia$Prior, Anorexia$Post, pair=T)
##
## Paired t-test
##
## data: Anorexia$Prior and Anorexia$Post
## t = -4.1849, df = 16, p-value = 0.0007003
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -10.94471 -3.58470
## sample estimates:
## mean of the differences
## -7.264706
# turn wide format to long format
dtaL <- gather(Anorexia, key = "Therapy", value = "weight", Prior, Post)
# show first 6 lines
head(dtaL)
## Therapy weight
## 1 Prior 83.8
## 2 Prior 83.3
## 3 Prior 86.0
## 4 Prior 82.5
## 5 Prior 86.7
## 6 Prior 79.6
#
ggplot(dtaL, aes(Therapy, weight)) +
geom_point(shape = 1) +
stat_summary(fun.data = mean_cl_boot, geom = "pointrange") +
coord_flip() +
labs(x = "Therapy", y = "Weight") +
theme_bw()
# show results of anova
summary(lm(weight ~ Therapy, data=dtaL))
##
## Call:
## lm(formula = weight ~ Therapy, data = dtaL)
##
## Residuals:
## Min 1Q Median 3Q Max
## -15.294 -2.454 1.106 4.004 11.106
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 90.494 1.689 53.578 < 2e-16 ***
## TherapyPrior -7.265 2.389 -3.041 0.00467 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.964 on 32 degrees of freedom
## Multiple R-squared: 0.2242, Adjusted R-squared: 0.2
## F-statistic: 9.25 on 1 and 32 DF, p-value: 0.004673
# the same results
summary(aov(weight ~ Therapy, data=dtaL))
## Df Sum Sq Mean Sq F value Pr(>F)
## Therapy 1 448.6 448.6 9.25 0.00467 **
## Residuals 32 1551.9 48.5
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# generate sbj id for each group
dtaL <- mutate(dtaL, Subject=rep(paste0("S", 1:17), 2), Sbj=paste0("S", 1:34))
# within subject design
summary(aov(weight ~ Therapy + Error(Subject/Therapy), data = dtaL))
##
## Error: Subject
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 16 1142 71.38
##
## Error: Subject:Therapy
## Df Sum Sq Mean Sq F value Pr(>F)
## Therapy 1 448.6 448.6 17.51 7e-04 ***
## Residuals 16 409.8 25.6
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# within subject design
summary(aov(weight ~ Therapy + Subject, data = dtaL))
## Df Sum Sq Mean Sq F value Pr(>F)
## Therapy 1 448.6 448.6 17.513 0.0007 ***
## Subject 16 1142.1 71.4 2.787 0.0240 *
## Residuals 16 409.8 25.6
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# between design
summary(aov(weight ~ Therapy + Sbj, data = dtaL))
## Df Sum Sq Mean Sq
## Therapy 1 448.6 448.6
## Sbj 32 1551.9 48.5
#The end
Exercise 2
#Data: ergoStool{nlme}
data("ergoStool", package = "nlme")
head(ergoStool)
## Grouped Data: effort ~ Type | Subject
## effort Type Subject
## 1 12 T1 1
## 2 15 T2 1
## 3 12 T3 1
## 4 10 T4 1
## 5 10 T1 2
## 6 14 T2 2
faraway::sumary(m0 <- lmer(effort ~ Type + (1 | Subject), data = ergoStool, REML=TRUE))
## Fixed Effects:
## coef.est coef.se
## (Intercept) 8.56 0.58
## TypeT2 3.89 0.52
## TypeT3 2.22 0.52
## TypeT4 0.67 0.52
##
## Random Effects:
## Groups Name Std.Dev.
## Subject (Intercept) 1.33
## Residual 1.10
## ---
## number of obs: 36, groups: Subject, 9
## AIC = 133.1, DIC = 123.2
## deviance = 122.1
faraway::sumary(m1 <- lmer(effort ~ Type + (1 | Subject)-1, data = ergoStool, REML=TRUE))
## Fixed Effects:
## coef.est coef.se
## TypeT1 8.56 0.58
## TypeT2 12.44 0.58
## TypeT3 10.78 0.58
## TypeT4 9.22 0.58
##
## Random Effects:
## Groups Name Std.Dev.
## Subject (Intercept) 1.33
## Residual 1.10
## ---
## number of obs: 36, groups: Subject, 9
## AIC = 133.1, DIC = 123.2
## deviance = 122.1
Exercise 3
pacman::p_load(car, tidyverse, lme4, GGally)
dta_e3 <- read.table('thetaEEG.txt', header = T)
str(dta_e3)
## '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 ...
eeg <- as.factor(colnames(dta_e3[, -1]))
m0 <- lm(as.matrix(dta_e3[, -1]) ~ 1)
m0_aov <- Anova(m0, idata = data.frame(eeg), idesign = ~ eeg)
## Note: model has only an intercept; equivalent type-III tests substituted.
summary(m0_aov)
##
## 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: eeg
##
## Response transformation matrix:
## eeg1 eeg2 eeg3 eeg4 eeg5 eeg6 eeg7 eeg8
## 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:
## eeg1 eeg2 eeg3 eeg4 eeg5 eeg6
## eeg1 5.7475 1.650000e-02 1.485000000 0.506000000 7.70000000 1.930500000
## eeg2 0.0165 4.736842e-05 0.004263158 0.001452632 0.02210526 0.005542105
## eeg3 1.4850 4.263158e-03 0.383684211 0.130736842 1.98947368 0.498789474
## eeg4 0.5060 1.452632e-03 0.130736842 0.044547368 0.67789474 0.169957895
## eeg5 7.7000 2.210526e-02 1.989473684 0.677894737 10.31578947 2.586315789
## eeg6 1.9305 5.542105e-03 0.498789474 0.169957895 2.58631579 0.648426316
## eeg7 3.0635 8.794737e-03 0.791526316 0.269705263 4.10421053 1.028984211
## eeg8 -0.0165 -4.736842e-05 -0.004263158 -0.001452632 -0.02210526 -0.005542105
## eeg7 eeg8
## eeg1 3.063500000 -1.650000e-02
## eeg2 0.008794737 -4.736842e-05
## eeg3 0.791526316 -4.263158e-03
## eeg4 0.269705263 -1.452632e-03
## eeg5 4.104210526 -2.210526e-02
## eeg6 1.028984211 -5.542105e-03
## eeg7 1.632889474 -8.794737e-03
## eeg8 -0.008794737 4.736842e-05
##
## Multivariate Tests: eeg
## 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 *
## eeg 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
## eeg 0.00037434 9.7367e-11
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## eeg 0.31531 0.6559
##
## HF eps Pr(>F[HF])
## eeg 0.3709213 0.6854028
dta_e3L <-dta_e3 %>%
gather(key = Channel, value = Theta, 2:10) %>%
mutate( Channel = factor(Channel))
head(dta_e3L)
## ID Channel Theta
## 1 1 Ch3 -3.54
## 2 2 Ch3 5.72
## 3 3 Ch3 0.52
## 4 4 Ch3 0.00
## 5 5 Ch3 2.07
## 6 6 Ch3 1.67
# set theme to black and white and save parameters
ot <- theme_set(theme_bw())
ggplot(dta_e3L, aes(reorder(Channel, Theta, median), Theta)) +
geom_boxplot() +
geom_point(alpha = .3) +
geom_hline(yintercept = mean(dtaL$Theta)) +
labs(x ="Channel Number", y = "Theta (Power)")
## Warning in mean.default(dtaL$Theta): argument is not numeric or logical:
## returning NA
## Warning: Removed 1 rows containing missing values (geom_hline).
# Univariate approach
summary(m1 <- aov(Theta ~ Channel + Error(ID/Channel), data = dta_e3L))
##
## Error: ID
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 1 112.2 112.2
##
## Error: ID:Channel
## Df Sum Sq Mean Sq
## Channel 8 23.13 2.891
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## Channel 8 9.8 1.23 0.184 0.993
## Residuals 153 1020.5 6.67
# Theta Power over Channel with CIs
library(Hmisc)
## Warning: package 'Hmisc' was built under R version 4.0.2
## Loading required package: survival
## Loading required package: Formula
## Warning: package 'Formula' was built under R version 4.0.2
##
## Attaching package: 'Hmisc'
## The following objects are masked from 'package:dplyr':
##
## src, summarize
## The following objects are masked from 'package:base':
##
## format.pval, units
ggplot(dta_e3L, aes(reorder(Channel, Theta, median), Theta)) +
stat_summary(fun.data = "mean_cl_boot", geom = "pointrange", col = "firebrick") +
geom_hline(yintercept = mean(dtaL$Theta), linetype = "dashed") +
geom_point(alpha = .3, size = rel(.7)) +
labs(x = "Channel Number", y = "Theta (Power)")
## Warning in mean.default(dtaL$Theta): argument is not numeric or logical:
## returning NA
## Warning: Removed 1 rows containing missing values (geom_hline).
# Patients’ theta power with channels colors
ggplot(dta_e3L, aes(reorder(ID, Theta, mean), Theta, color = Channel)) +
geom_point() +
geom_hline(yintercept = mean(dtaL$Theta), linetype = "dashed") +
coord_flip() +
labs(y = "Theta (Power)", x = "Subject ID") +
theme(legend.position = c(.9, .2))
## Warning in mean.default(dtaL$Theta): argument is not numeric or logical:
## returning NA
## Warning: Removed 1 rows containing missing values (geom_hline).
Exercise 4
# read data
dta <- read.table("cognitive_task.txt", h=T)
# plot variance with Nested data
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
# 將每個班級的差異考量進去
m0 <- aov(Score ~ Method + Error(Class), data=dta)
#
m01 <- aov(Score ~ Method + Error(Method/Klass), 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
#
summary(m01)
##
## Error: Method
## Df Sum Sq Mean Sq
## Method 1 112.5 112.5
##
## Error: Method:Klass
## 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 <- lme4::lmer(Score ~ Method + (1 | Class), data=dta)
#
m11 <- lme4::lmer(Score ~ Method + (1 | Method:Klass), data=dta)
#
confint(m1, method="boot")
## Computing bootstrap confidence intervals ...
## 2.5 % 97.5 %
## .sig01 0.9245146 3.154740
## .sigma 0.6179640 1.131203
## (Intercept) 1.2612380 5.504713
## MethodI2 0.9716536 6.565271