library(tidyverse)
## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --
## v ggplot2 3.3.5 v purrr 0.3.4
## v tibble 3.1.5 v dplyr 1.0.7
## v tidyr 1.1.4 v stringr 1.4.0
## v readr 2.0.2 v forcats 0.5.1
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
library(ggpubr)
library(rstatix)
##
## Attaching package: 'rstatix'
## The following object is masked from 'package:stats':
##
## filter
data("anxiety",package="datarium")
anxiety
## # A tibble: 45 x 5
## id group t1 t2 t3
## <fct> <fct> <dbl> <dbl> <dbl>
## 1 1 grp1 14.1 14.4 14.1
## 2 2 grp1 14.5 14.6 14.3
## 3 3 grp1 15.7 15.2 14.9
## 4 4 grp1 16 15.5 15.3
## 5 5 grp1 16.5 15.8 15.7
## 6 6 grp1 16.9 16.5 16.2
## 7 7 grp1 17 16.8 16.5
## 8 8 grp1 17 17.1 16.6
## 9 9 grp1 17.3 16.9 16.5
## 10 10 grp1 17.3 17.1 16.7
## # ... with 35 more rows
anxiety<- anxiety %>%
gather(key="time",value="score",t1,t2,t3) %>%
convert_as_factor(id,time)
anxiety
## # A tibble: 135 x 4
## id group time score
## <fct> <fct> <fct> <dbl>
## 1 1 grp1 t1 14.1
## 2 2 grp1 t1 14.5
## 3 3 grp1 t1 15.7
## 4 4 grp1 t1 16
## 5 5 grp1 t1 16.5
## 6 6 grp1 t1 16.9
## 7 7 grp1 t1 17
## 8 8 grp1 t1 17
## 9 9 grp1 t1 17.3
## 10 10 grp1 t1 17.3
## # ... with 125 more rows
anxiety %>% group_by(time) %>% get_summary_stats(score,type="mean_sd")
## # A tibble: 3 x 5
## time variable n mean sd
## <fct> <chr> <dbl> <dbl> <dbl>
## 1 t1 score 45 16.9 1.49
## 2 t2 score 45 16.1 1.77
## 3 t3 score 45 15.2 1.97
bxp <- ggboxplot(
anxiety, x="time", y="score",
color="group",pallete="jco"
)
bxp
anxiety %>% group_by(time,group) %>% identify_outliers(score)
## [1] group time id score is.outlier is.extreme
## <0 rows> (or 0-length row.names)
anxiety %>% group_by(time,group) %>% shapiro_test(score)
## # A tibble: 9 x 5
## group time variable statistic p
## <fct> <fct> <chr> <dbl> <dbl>
## 1 grp1 t1 score 0.964 0.769
## 2 grp2 t1 score 0.977 0.949
## 3 grp3 t1 score 0.954 0.588
## 4 grp1 t2 score 0.956 0.624
## 5 grp2 t2 score 0.935 0.328
## 6 grp3 t2 score 0.952 0.558
## 7 grp1 t3 score 0.949 0.506
## 8 grp2 t3 score 0.909 0.131
## 9 grp3 t3 score 0.925 0.232
ggqqplot(anxiety,"score",ggtheme=theme_bw()) + facet_grid(time~group)
anxiety %>% group_by(time) %>% levene_test(score~group)
## # A tibble: 3 x 5
## time df1 df2 statistic p
## <fct> <int> <int> <dbl> <dbl>
## 1 t1 2 42 0.176 0.839
## 2 t2 2 42 0.249 0.781
## 3 t3 2 42 0.335 0.717
res.aov<-anova_test(data=anxiety,dv=score,wid=id,between = group,within=time)
get_anova_table(res.aov)
## ANOVA Table (type II tests)
##
## Effect DFn DFd F p p<.05 ges
## 1 group 2 42 4.352 1.90e-02 * 0.168
## 2 time 2 84 394.909 1.91e-43 * 0.179
## 3 group:time 4 84 110.188 1.38e-32 * 0.108
one.way<- anxiety %>% group_by(time) %>% anova_test(dv=score,wid=id,between=group) %>%
get_anova_table() %>% adjust_pvalue(method="bonferroni")
## Coefficient covariances computed by hccm()
## Coefficient covariances computed by hccm()
## Coefficient covariances computed by hccm()
one.way
## # A tibble: 3 x 9
## time Effect DFn DFd F p `p<.05` ges p.adj
## * <fct> <chr> <dbl> <dbl> <dbl> <dbl> <chr> <dbl> <dbl>
## 1 t1 group 2 42 0.365 0.696 "" 0.017 1
## 2 t2 group 2 42 5.84 0.006 "*" 0.218 0.018
## 3 t3 group 2 42 13.8 0.0000248 "*" 0.396 0.0000744