I am trying to do a repeated-measures multivariate ANOVA on some data, however it is not normally distributed therefore, I am trying to use bootstrapping to be more confident about my results. The error I get is when I try to apply my model to the bootstrapped data.
Thank you very much
library(tidyverse)
library(ggpubr)
library(rstatix)
#>
#> Attaching package: 'rstatix'
#> The following object is masked from 'package:stats':
#>
#> filter
library(rsample)
library(reprex)
#Load data
allA<-read.csv('alldata_A_R.csv')
allA<-as.data.frame(allA)
head(allA) #to show you what my raw data looks like
#> subject exp Oddball PvN fitA
#> 1 1 1 E N 7.998360e-02
#> 2 2 1 E N 5.507097e-01
#> 3 3 1 E N 7.541230e-11
#> 4 4 1 E N 1.512817e-01
#> 5 5 1 E N 1.162404e-01
#> 6 6 1 E N 9.120629e-02
allA %>%
group_by(exp,Oddball,PvN) %>%
get_summary_stats(fitA,type='mean_sd')
#> # A tibble: 8 x 7
#> exp Oddball PvN variable n mean sd
#> <int> <chr> <chr> <chr> <dbl> <dbl> <dbl>
#> 1 1 E N fitA 70 0.117 0.069
#> 2 1 E P fitA 70 0.158 0.115
#> 3 1 P N fitA 70 0.111 0.045
#> 4 1 P P fitA 70 0.163 0.069
#> 5 2 E N fitA 30 0.099 0.047
#> 6 2 E P fitA 30 0.153 0.116
#> 7 2 P N fitA 30 0.113 0.042
#> 8 2 P P fitA 30 0.144 0.066
bxp<-ggboxplot(allA, x='Oddball',y='fitA',
color='PvN', palette='jco',
facet.by='exp', short.panel.labs=FALSE)
bxp
[image]
#check assumptions
allA %>%
group_by(exp,Oddball,PvN) %>%
identify_outliers(fitA) #test for outliers
#> # A tibble: 12 x 7
#> exp Oddball PvN subject fitA is.outlier is.extreme
#> <int> <chr> <chr> <int> <dbl> <lgl> <lgl>
#> 1 1 E N 2 5.51e- 1 TRUE TRUE
#> 2 1 E N 3 7.54e-11 TRUE FALSE
#> 3 1 E N 62 2.30e- 1 TRUE FALSE
#> 4 1 E P 7 4.15e- 1 TRUE TRUE
#> 5 1 E P 39 4.38e- 1 TRUE TRUE
#> 6 1 E P 57 8.86e- 1 TRUE TRUE
#> 7 1 P P 11 3.00e- 1 TRUE FALSE
#> 8 1 P P 17 2.70e- 1 TRUE FALSE
#> 9 1 P P 39 5.63e- 1 TRUE TRUE
#> 10 1 P P 68 3.18e- 1 TRUE FALSE
#> 11 2 E P 12 6.77e- 1 TRUE TRUE
#> 12 2 P N 28 2.19e- 1 TRUE FALSE
#test for normality
normality<-allA %>%
group_by(exp,Oddball,PvN) %>%
shapiro_test(fitA)
normality
#> # A tibble: 8 x 6
#> exp Oddball PvN variable statistic p
#> <int> <chr> <chr> <chr> <dbl> <dbl>
#> 1 1 E N fitA 0.713 2.16e-10
#> 2 1 E P fitA 0.634 6.23e-12
#> 3 1 P N fitA 0.971 1.10e- 1
#> 4 1 P P fitA 0.753 1.60e- 9
#> 5 2 E N fitA 0.986 9.60e- 1
#> 6 2 E P fitA 0.674 6.72e- 7
#> 7 2 P N fitA 0.976 7.16e- 1
#> 8 2 P P fitA 0.945 1.26e- 1
#Create QQ plot for each cell of design
ggqqplot(allA,'fitA',ggtheme=theme_bw())+
facet_grid(exp+Oddball~PvN, labeller='label_both')
[image]
res.aov<-anova_test(
data=allA,dv=fitA,wid=subject,between=c(exp),
within=c(Oddball,PvN)
)
#> Warning: The 'wid' column contains duplicate ids across between-subjects
#> variables. Automatic unique id will be created
get_anova_table(res.aov) #this give the EXACT same results as JASP.
#> ANOVA Table (type III tests)
#>
#> Effect DFn DFd F p p<.05 ges
#> 1 exp 1 98 1.105 2.96e-01 4.00e-03
#> 2 Oddball 1 98 0.021 8.84e-01 3.80e-05
#> 3 PvN 1 98 22.963 5.88e-06 * 6.60e-02
#> 4 exp:Oddball 1 98 0.049 8.25e-01 8.74e-05
#> 5 exp:PvN 1 98 0.041 8.40e-01 1.26e-04
#> 6 Oddball:PvN 1 98 0.167 6.84e-01 3.42e-04
#> 7 exp:Oddball:PvN 1 98 1.385 2.42e-01 3.00e-03
##test bootstrapping here
set.seed(123)
dataA.boot<-bootstraps(allA,times=1000,apparent=TRUE)
#dataA.boot<-as.data.frame(dataA.boot)
#fit a model to all of these bootstraps
dataA.boot %>%
mutate(model=map(splits,~anova_test(
data=.,dv=fitA,wid=subject,between=c(exp),
within=c(Oddball,PvN))),
coef_info=map(model,tidy))
#> Error: Problem with `mutate()` input `model`.
#> x Can't subset columns that don't exist.
#> x Column `fitA` doesn't exist.
#> ℹ Input `model` is `map(...)`.
#post-hoc compairson on the main effect
PvN.pwc<-allA %>%
group_by(exp,Oddball)%>%
pairwise_t_test(fitA~PvN, paired=TRUE, p.adjust.method='bonferroni') %>%
select(-df,-statistic) #remove details
PvN.pwc
#> # A tibble: 4 x 10
#> exp Oddball .y. group1 group2 n1 n2 p p.adj p.adj.signif
#> <int> <chr> <chr> <chr> <chr> <int> <int> <dbl> <dbl> <chr>
#> 1 1 E fitA N P 70 70 1.40e-2 1.40e-2 *
#> 2 1 P fitA N P 70 70 1.59e-6 1.59e-6 ****
#> 3 2 E fitA N P 30 30 2.70e-2 2.70e-2 *
#> 4 2 P fitA N P 30 30 3.80e-2 3.80e-2 *