library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
## ✔ ggplot2 3.3.6 ✔ purrr 0.3.4
## ✔ tibble 3.1.8 ✔ dplyr 1.0.9
## ✔ tidyr 1.2.0 ✔ stringr 1.4.0
## ✔ readr 2.1.2 ✔ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
library(ggpubr)
library(rstatix)
##
## Attaching package: 'rstatix'
## The following object is masked from 'package:stats':
##
## filter
library(readxl)
dat <- read_excel("/Users/riku/Documents/zxy.xlsx")
dat$id <- as.factor(dat$id)
dat$LH <- as.factor(dat$LH)
dat$PT <- as.factor(dat$PT)
dat$AB <- as.factor(dat$AB)
dat$LH <- ordered(dat$LH, levels=c("L","H"))
dat <- dat %>%
gather(key = "time", value = "score", "1st", "2rd") %>%
convert_as_factor(id, time)
dat$score <- as.numeric(dat$score)
# PS ----
## ver A + PS ----
dat_a <- dat %>% filter(dat$AB == "A")
dat_a_p <- dat_a %>% filter(dat_a$PT == "P")
bxp <- ggboxplot(
dat_a_p, x = "time", y = "score", add = "jitter",
color = "LH", palette = "lancet"
)
bxp + ggtitle("PS Ver A.")
### Check assumptions -------
dat_a_p %>%
group_by(time, LH) %>%
identify_outliers(score)
## # A tibble: 6 × 9
## LH time id PT AB sa score is.outlier is.extreme
## <ord> <fct> <fct> <fct> <fct> <dbl> <dbl> <lgl> <lgl>
## 1 H 1st 10 P A -1 18 TRUE FALSE
## 2 L 2rd 16 P A -2 9 TRUE FALSE
## 3 H 2rd 10 P A -1 17 TRUE FALSE
## 4 H 2rd 24 P A 0 11 TRUE FALSE
## 5 H 2rd 32 P A 2 11 TRUE FALSE
## 6 H 2rd 36 P A 5 16 TRUE FALSE
dat_a_p %>%
group_by(time, LH) %>%
shapiro_test(score)
## # A tibble: 4 × 5
## LH time variable statistic p
## <ord> <fct> <chr> <dbl> <dbl>
## 1 L 1st score 0.922 0.234
## 2 H 1st score 0.926 0.443
## 3 L 2rd score 0.887 0.0730
## 4 H 2rd score 0.927 0.454
ggqqplot(dat_a_p, "score", ggtheme = theme_bw()) +
facet_grid(time ~ LH)
dat_a_p %>%
group_by(time) %>%
levene_test(score ~ LH)
## # A tibble: 2 × 5
## time df1 df2 statistic p
## <fct> <int> <int> <dbl> <dbl>
## 1 1st 1 21 0.00260 0.960
## 2 2rd 1 21 0.000540 0.982
### Computation -------
res.aov <- anova_test(
data = dat_a_p, dv = score, wid = id,
between = LH, within = time
)
get_anova_table(res.aov)
## ANOVA Table (type III tests)
##
## Effect DFn DFd F p p<.05 ges
## 1 LH 1 21 0.723 0.405 0.025
## 2 time 1 21 2.900 0.103 0.032
## 3 LH:time 1 21 0.298 0.591 0.003
## ver B + PS -----
dat_b <- dat %>% filter(dat$AB == "B")
dat_b_p <- dat_b %>% filter(dat_b$PT == "P")
bxp <- ggboxplot(
dat_b_p, x = "time", y = "score", add = "jitter",
color = "LH", palette = "lancet"
)
bxp + ggtitle("PS Ver B.")
### Check assumptions -------
dat_b_p %>%
group_by(time, LH) %>%
identify_outliers(score)
## # A tibble: 1 × 9
## LH time id PT AB sa score is.outlier is.extreme
## <ord> <fct> <fct> <fct> <fct> <dbl> <dbl> <lgl> <lgl>
## 1 H 2rd 23 P B -2 8 TRUE TRUE
dat_b_p %>%
group_by(time, LH) %>%
shapiro_test(score)
## # A tibble: 4 × 5
## LH time variable statistic p
## <ord> <fct> <chr> <dbl> <dbl>
## 1 L 1st score 0.926 0.376
## 2 H 1st score 0.963 0.830
## 3 L 2rd score 0.923 0.347
## 4 H 2rd score 0.801 0.00971
ggqqplot(dat_b_p, "score", ggtheme = theme_bw()) +
facet_grid(time ~ LH)
dat_b_p %>%
group_by(time) %>%
levene_test(score ~ LH)
## # A tibble: 2 × 5
## time df1 df2 statistic p
## <fct> <int> <int> <dbl> <dbl>
## 1 1st 1 21 0.0505 0.824
## 2 2rd 1 21 1.24 0.278
### Computation -------
res.aov <- anova_test(
data = dat_b_p, dv = score, wid = id,
between = LH, within = time
)
get_anova_table(res.aov)
## ANOVA Table (type III tests)
##
## Effect DFn DFd F p p<.05 ges
## 1 LH 1 21 4.769 0.040 * 0.127
## 2 time 1 21 7.083 0.015 * 0.108
## 3 LH:time 1 21 0.079 0.781 0.001
# TS -----
## ver A + TS ----
dat_a_t <- dat_a %>% filter(dat_a$PT == "T")
bxp <- ggboxplot(
dat_a_t, x = "time", y = "score", add = "jitter",
color = "LH", palette = "lancet"
)
bxp + ggtitle("TS Ver A.")
### Check assumptions -------
dat_a_t %>%
group_by(time, LH) %>%
identify_outliers(score)
## # A tibble: 1 × 9
## LH time id PT AB sa score is.outlier is.extreme
## <ord> <fct> <fct> <fct> <fct> <dbl> <dbl> <lgl> <lgl>
## 1 L 2rd 9 T A 6 19 TRUE FALSE
dat_a_t %>%
group_by(time, LH) %>%
shapiro_test(score)
## # A tibble: 4 × 5
## LH time variable statistic p
## <ord> <fct> <chr> <dbl> <dbl>
## 1 L 1st score 0.909 0.238
## 2 H 1st score 0.958 0.758
## 3 L 2rd score 0.888 0.131
## 4 H 2rd score 0.932 0.402
ggqqplot(dat_a_t, "score", ggtheme = theme_bw()) +
facet_grid(time ~ LH)
dat_a_t %>%
group_by(time) %>%
levene_test(score ~ LH)
## # A tibble: 2 × 5
## time df1 df2 statistic p
## <fct> <int> <int> <dbl> <dbl>
## 1 1st 1 21 4.46 0.0468
## 2 2rd 1 21 0.110 0.744
### Computation -------
res.aov <- anova_test(
data = dat_a_t, dv = score, wid = id,
between = LH, within = time
)
get_anova_table(res.aov)
## ANOVA Table (type III tests)
##
## Effect DFn DFd F p p<.05 ges
## 1 LH 1 21 0.972 0.335000 0.036
## 2 time 1 21 17.544 0.000414 * 0.136
## 3 LH:time 1 21 0.379 0.545000 0.003
## ver B + PS -----
dat_b_t <- dat_b %>% filter(dat_b$PT == "T")
bxp <- ggboxplot(
dat_b_t, x = "time", y = "score", add = "jitter",
color = "LH", palette = "lancet"
)
bxp + ggtitle("TS Ver B.")
### Check assumptions -------
dat_b_t %>%
group_by(time, LH) %>%
identify_outliers(score)
## # A tibble: 1 × 9
## LH time id PT AB sa score is.outlier is.extreme
## <ord> <fct> <fct> <fct> <fct> <dbl> <dbl> <lgl> <lgl>
## 1 H 2rd 10 T B 0 12 TRUE FALSE
dat_b_t %>%
group_by(time, LH) %>%
shapiro_test(score)
## # A tibble: 4 × 5
## LH time variable statistic p
## <ord> <fct> <chr> <dbl> <dbl>
## 1 L 1st score 0.952 0.585
## 2 H 1st score 0.854 0.0819
## 3 L 2rd score 0.938 0.399
## 4 H 2rd score 0.837 0.0536
ggqqplot(dat_b_t, "score", ggtheme = theme_bw()) +
facet_grid(time ~ LH)
dat_b_t %>%
group_by(time) %>%
levene_test(score ~ LH)
## # A tibble: 2 × 5
## time df1 df2 statistic p
## <fct> <int> <int> <dbl> <dbl>
## 1 1st 1 21 1.89 0.184
## 2 2rd 1 21 1.10 0.306
### Computation -------
res.aov <- anova_test(
data = dat_b_t, dv = score, wid = id,
between = LH, within = time
)
get_anova_table(res.aov)
## ANOVA Table (type III tests)
##
## Effect DFn DFd F p p<.05 ges
## 1 LH 1 21 6.053 0.023000 * 0.174
## 2 time 1 21 15.381 0.000783 * 0.165
## 3 LH:time 1 21 0.296 0.592000 0.004
# PS/TS Between L and H ----
## ver A + PSTS ----
dat_2 <- read_excel("/Users/riku/Documents/zxy.xlsx")
dat_2$id <- as.factor(dat_2$id)
dat_2$LH <- as.factor(dat_2$LH)
dat_2$PT <- as.factor(dat_2$PT)
dat_2$AB <- as.factor(dat_2$AB)
dat_2$LH <- ordered(dat_2$LH, levels=c("L","H"))
dat_2$score <- as.numeric(dat_2$score)
dat_2_a <- dat_2 %>% filter(dat_2$AB == "A")
bxp <- ggboxplot(
dat_2_a, x = "LH", y = "sa", add = "jitter",
color = "PT", palette = "lancet"
)
bxp + ggtitle("PSTS Ver A.")
### Check assumptions -------
dat_2_a %>%
group_by(LH, PT) %>%
identify_outliers(sa)
## # A tibble: 1 × 11
## LH PT id `1st` `2rd` AB sa time score is.outlier is.extreme
## <ord> <fct> <fct> <chr> <dbl> <fct> <dbl> <lgl> <dbl> <lgl> <lgl>
## 1 L P 22 9 15 A 6 NA NA TRUE FALSE
dat_2_a %>%
group_by(LH, PT) %>%
shapiro_test(sa)
## # A tibble: 4 × 5
## LH PT variable statistic p
## <ord> <fct> <chr> <dbl> <dbl>
## 1 L P sa 0.906 0.138
## 2 L T sa 0.873 0.0844
## 3 H P sa 0.979 0.957
## 4 H T sa 0.916 0.251
ggqqplot(dat_2_a, "sa", ggtheme = theme_bw()) +
facet_grid(LH ~ PT)
dat_2_a %>%
group_by(PT) %>%
levene_test(sa ~ LH)
## # A tibble: 2 × 5
## PT df1 df2 statistic p
## <fct> <int> <int> <dbl> <dbl>
## 1 P 1 21 0.230 0.636
## 2 T 1 21 0.0624 0.805
### Computation -------
res.aov <- dat_2_a %>% anova_test(sa ~ LH * PT)
## Coefficient covariances computed by hccm()
res.aov
## ANOVA Table (type II tests)
##
## Effect DFn DFd F p p<.05 ges
## 1 LH 1 42 0.655 0.423 1.50e-02
## 2 PT 1 42 1.504 0.227 3.50e-02
## 3 LH:PT 1 42 0.003 0.955 7.52e-05
## ver B + PSTS ----
dat_2_b <- dat_2 %>% filter(dat_2$AB == "B")
bxp <- ggboxplot(
dat_2_b, x = "LH", y = "sa", add = "jitter",
color = "PT", palette = "lancet"
)
bxp + ggtitle("PSTS Ver B.")
### Check assumptions -------
dat_2_b %>%
group_by(LH, PT) %>%
identify_outliers(sa)
## # A tibble: 1 × 11
## LH PT id `1st` `2rd` AB sa time score is.outlier is.extreme
## <ord> <fct> <fct> <chr> <dbl> <fct> <dbl> <lgl> <dbl> <lgl> <lgl>
## 1 H T 24 13 17 B 4 NA NA TRUE FALSE
dat_2_b %>%
group_by(LH, PT) %>%
shapiro_test(sa)
## # A tibble: 4 × 5
## LH PT variable statistic p
## <ord> <fct> <chr> <dbl> <dbl>
## 1 L P sa 0.887 0.128
## 2 L T sa 0.966 0.818
## 3 H P sa 0.929 0.372
## 4 H T sa 0.941 0.595
ggqqplot(dat_2_b, "sa", ggtheme = theme_bw()) +
facet_grid(LH ~ PT)
dat_2_b %>%
group_by(PT) %>%
levene_test(sa ~ LH)
## # A tibble: 2 × 5
## PT df1 df2 statistic p
## <fct> <int> <int> <dbl> <dbl>
## 1 P 1 21 0.0152 0.903
## 2 T 1 21 4.66 0.0426
### Computation -------
res.aov <- dat_2_b %>% anova_test(sa ~ LH * PT)
## Coefficient covariances computed by hccm()
res.aov
## ANOVA Table (type II tests)
##
## Effect DFn DFd F p p<.05 ges
## 1 LH 1 42 0.012 0.914 0.000282
## 2 PT 1 42 0.114 0.738 0.003000
## 3 LH:PT 1 42 0.312 0.579 0.007000