# Import Data
df <- read_excel(DATA)
# Define color
my_colours = brewer.pal(n = 9, PALETTE)[9:3]
# Determine data normality using the Shapiro test
shapiro_df <- df %>%group_by(Treatment, Line)%>%
summarise(statistic = shapiro.test(AB)$statistic,
p.value = shapiro.test(AB)$p.value)
flag_normal <- check_normality(shapiro_df)
# Data treatement according to normality status
if(flag_normal == TRUE) {
#print("Data follow a normal distribution")
# Summary
my_summary <- summarySE(df, measurevar="AB", groupvars=c("Line", "Treatment"))
# Plot
plot_normal(df, my_colours, my_summary)
# Stats
anova_results <- df %>% group_by(Treatment) %>% anova_test(AB ~ Line)
tukey_results <- df %>% group_by(Treatment) %>% tukey_hsd(AB ~ Line)
# Save plot
ggsave("my_ggplot1.pdf", width=4, height=5)
ggsave("my_ggplot1.svg", width=4, height=5)
} else {
#print("Data do not follow a normal distribution")
# Calculate Median and 95%CI by bootstrap (Largely based on code by Peter Kamerman:
#https://www.painblogr.org/2017-10-18-purrring-through-bootstraps.html)
#1. Set the seed
set.seed(12345)
#2. Create nested datafile to allow boostrapping on each subgroup.
my_data_nested <- df %>% group_by(Line,Treatment) %>% nest()
#3. Add boostrapped median values to my_nested_data
my_data_nested %<>% dplyr::mutate(booted = purrr::map(.x = data, ~ boot::boot(data = .x$AB,statistic = samplemedian, R = 10000, stype = "i")))
#4. Plot each <S3: boot> object to check boostrapping worked!
## Note: Saved to an object (plots) to stop the summary being printed
plots<-purrr::map(.x = my_data_nested$booted, ~ plot(.x))
print(plots)
#5. Add boostrapped 95%CI around median values to my_nested_data,
#95%CI determined based on "perc" & "bca"
my_data_nested %<>% mutate(booted_ci = purrr::map(.x = booted,
~ boot::boot.ci(.x,conf = 0.95,
type = c("perc", "bca"))))
#6. Table of each <S3: bootci> object to check
ci_table <-purrr::map(.x = my_data_nested$booted_ci, ~ print(.x))
#7. Get the stats (median and 95% CI) from the nested dataframe into a regular dataframe
#(my_data_booted_summary) with 95%CI determined by "bca" or "perc" method
booted_summary <- my_data_nested %>% mutate(AB= purrr::map(.x = booted_ci, ~ .x$t0),
lower_ci_bca = purrr::map(.x = booted_ci,~ .x$bca[[4]]),
upper_ci_bca = purrr::map(.x = booted_ci,~ .x$bca[[5]]),
lower_ci_perc = purrr::map(.x = booted_ci,~ .x$perc[[4]]),
upper_ci_perc = purrr::map(.x = booted_ci,~ .x$perc[[5]])) %>%
dplyr::select(-data, -booted, -booted_ci) %>% tidyr::unnest(cols = c(AB, lower_ci_bca, upper_ci_bca, lower_ci_perc, upper_ci_perc))
# Plot
plot_not_normal(df, my_colours, booted_summary)
# Stats
kruskal_pval <- (df %>% group_by(Treatment)%>%kruskal_test(AB ~ Line)) %>% select(Treatment, p)
pval_dunn <- test_dunn()
# Save plot
ggsave("my_ggplot1.pdf", width=4, height=5)
ggsave("my_ggplot1.svg", width=4, height=5)
}
