I provide the original R code in case you want to check the logic of simulation.

Data Manipulation

Baseline Data

  • Assume we recruit 208 students.
n <- 60 # total number of participants
id <- c(1:n)

set.seed(123) # Do not change this number. This line of code ensures that the sampled data are the same when we use sample() function.

# Basic Information

group <- sample(c("None","Once", "Twice"), n, replace = TRUE) # randomly assign participants into three groups

## The following simulated data do not include "Prefer not to say"
gender <- sample(c("Male", "Female"), n, replace = TRUE)
age <- sample(c(18:22), n, replace = TRUE) # age can vary if we only include Oxford students
grade <- sample(c("Freshman", "Sophomore", "Junior","Senior"), n, replace = TRUE)
race <- sample(c("White", "Black","Indian_Alaska","Asian","Islander"), n, replace = TRUE, prob = c(.6,.2,.025,.15,.025))
ethnicity <- sample(c("Hisp","Non_Hisp"), n, replace = TRUE, prob = c(.1,.9))

# Academic Information
major <- sample(c("Psychology", "Pre-medical", "Biology", "QSS", "Sociology"), n, replace = TRUE, prob = c(.5,.2,.15,.1,.05))

credit_hr <- sample(c(15:21), n, replace = TRUE, prob = c(.1,.2,.25,.2,.1,.1,.05))
course_hd <- sample(c("Chemistry", "Biology", "Math", "Physics", "Philosophy"), n, replace = TRUE, prob = c(.2,.3,.3,.1,.1))

# Meditation Experiences
med_freq <- sample(c("Daily", "Weekly", "Monthly","Never", "Other"), n, replace = TRUE, prob = c(.05,.2,.1,.6,.05))
med_exp <- sample(c("Never", "1~10", "10~100","100~1000", "1000+"), n, replace = TRUE, prob = c(.05,.2,.1,.6,.05))
med_type <- sample(c("Mindfulness", "Loving", "FA", "Unknown", "Never"), n, replace = TRUE, prob = c(.1,.05,.05,.2,.6))

# Meditation Attitude
# I think the response to each question may be linked
relax_q <- sample(c(1:5), n, replace = TRUE)
focus_q <- sample(c(1:5), n, replace = TRUE)
score_q <- sample(c(1:5), n, replace = TRUE)

baseline <- data.frame(id, group, gender, age, grade, race, ethnicity, # basic
                       major, credit_hr, course_hd, # academic
                       med_freq, med_exp, med_type, # experience
                       relax_q, focus_q, score_q # attitude
                       )

knitr::kable(head(baseline,10)) # view data set
id group gender age grade race ethnicity major credit_hr course_hd med_freq med_exp med_type relax_q focus_q score_q
1 Twice Male 20 Freshman Black Non_Hisp Pre-medical 17 Philosophy Never 10~100 Never 1 2 4
2 Twice Male 19 Sophomore White Non_Hisp Pre-medical 20 Biology Weekly 10~100 Unknown 1 4 2
3 Twice Male 19 Freshman Black Non_Hisp Pre-medical 16 Philosophy Never 100~1000 Unknown 5 2 4
4 Once Male 19 Senior Black Non_Hisp Psychology 15 Biology Never 1~10 Mindfulness 5 1 3
5 Twice Female 21 Junior Black Non_Hisp QSS 19 Chemistry Never 100~1000 Unknown 4 3 1
6 Once Female 19 Freshman White Non_Hisp Pre-medical 17 Physics Monthly 100~1000 Unknown 4 5 3
7 Once Male 19 Freshman White Non_Hisp Psychology 17 Biology Never 1000+ Mindfulness 4 3 4
8 Once Female 21 Senior White Non_Hisp Psychology 15 Philosophy Monthly 100~1000 Unknown 5 1 2
9 Twice Female 21 Senior White Non_Hisp Psychology 16 Philosophy Never 1~10 Unknown 4 2 5
10 None Female 18 Sophomore White Non_Hisp QSS 20 Math Never 100~1000 Never 1 4 2

Daily Mood Survey (day 0 ~ 21)

daily.times <- 22
days <- daily.times - 1
id.mood <- rep(id,times=daily.times)
timepoint <- rep(c(0:days),each=n)
workload <- sample(c(T,F), n * daily.times, replace = TRUE,prob = c(.1,.9)) # question 1: use true or false to replace yes or no
# Whether the student has major assignment(s) or not


# completion rate of daily mood survey
probs.t <- seq(.8,.6,length.out = daily.times)
completion <- c()
for(i in probs.t){
    completion <- c(completion, sample(c(T,F),n,replace = TRUE, prob = c(i,1-i)))
}


# This function generates scores of likert scale questions by using normal distribution (can be used for daily mood survey and AEQ)
score_generator <- function(order, from, to,a=-Inf, b=Inf, times, negative = T){ # generate score from day 0 to day 21, on a scale from 1 to 5
    none = order %>%
              filter(group == "None")
    once = order %>%
              filter(group == "Once")
    twice = order %>%
              filter(group == "Twice")
    
    means.none = seq(from, from + 0.1, length.out=times) # nearly no change
    means.once = seq(from, to, length.out=times) # blocked meditation
    
    means.twice = c()
    if(negative){
        means.twice = seq(from, to * 0.7, length.out=times) # intermittent meditation: lower change than blocked meditation
    }
    else{
        means.twice = seq(from, to * 1.3, length.out=times) # intermittent meditation: higher change than blocked meditation
    }
    
    
    score = c()
    for(i in c(1:times)){ # 22 days
        score.none = rtruncnorm(n=length(none$id), a=a, b=b, mean=means.none[i], sd=1.5)
        score.once = rtruncnorm(n=length(once$id), a=a, b=b, mean=means.once[i], sd=1.5)
        score.twice = rtruncnorm(n=length(twice$id), a=a, b=b, mean=means.twice[i], sd=1.5)
        score = c(score, score.none, score.once, score.twice) # append a normally distributed score range from a to b
    }
    return(score)
}

# reorder the id number by group in ascending order
order.group = baseline %>%
                select(id, group) %>%
                arrange(group,id)

mood <- round(score_generator(order.group, 2.5, 4.5, 1, 5, times=daily.times,F)) 
study.length <- round(score_generator(order.group, 4, 7, 0, times=daily.times,F),1) # record in hours.
productivity <- round(score_generator(order.group, 1.5, 4.5, 1, 5, times=daily.times,F)) 


# final dataset

order.daily <- data.frame(id = rep(order.group$id, times=daily.times),group = rep(order.group$group, times=daily.times))

daily.mood <- data.frame(id = order.daily$id, group = order.daily$group, timepoint, workload, completion, mood, study.length, productivity) %>%
    arrange(id)
daily.mood$completion[daily.mood$group=="None"] = NA
knitr::kable(head(daily.mood,10)) # view data set
id group timepoint workload completion mood study.length productivity
1 Twice 0 FALSE TRUE 5 NA 3
1 Twice 1 FALSE TRUE 3 NA 1
1 Twice 2 FALSE FALSE 3 NA 2
1 Twice 3 FALSE TRUE 2 NA 2
1 Twice 4 TRUE TRUE 4 NA 1
1 Twice 5 FALSE FALSE 3 NA 4
1 Twice 6 FALSE TRUE 3 NA 4
1 Twice 7 FALSE TRUE 4 NA 4
1 Twice 8 FALSE TRUE 2 NA 3
1 Twice 9 TRUE TRUE 4 NA 4

AEQ (day 0, 6, 7, 14, 21)

AEQ.times = 5
id.AEQ <- rep(id,times = AEQ.times) # administer 5 times
timepoint <- factor(x = rep(c(0,6,7,14,21), each = n))

# The total score of each trait is 20.
enjoyment <- round(score_generator(order.group, 8, 13.5, 1, 20, times=AEQ.times,F))
hope <- round(score_generator(order.group, 12, 15.5, 1, 20, times=AEQ.times,F))
pride <- round(score_generator(order.group, 9, 15.5, 1, 20, times=AEQ.times,F))

anger <- round(score_generator(order.group, 17.5, 9, 1, 20, times=AEQ.times))
anxiety <- round(score_generator(order.group, 17.5, 10, 1, 20, times=AEQ.times))
shame <- round(score_generator(order.group, 8, 2, 1, 20, times=AEQ.times))
hopelessness <- round(score_generator(order.group, 10, 7, 1, 20, times=AEQ.times))
boredom <- round(score_generator(order.group, 9, 7, 1, 20, times=AEQ.times))


order.AEQ <- data.frame(id = rep(order.group$id, times=AEQ.times),group = rep(order.group$group, times=AEQ.times))

AEQ <- data.frame(id = order.AEQ$id, group = order.AEQ$group, timepoint, enjoyment, hope, pride, anger, anxiety, shame, hopelessness, boredom) %>% 
    pivot_longer(c(enjoyment, hope, pride, anger, anxiety, shame, hopelessness, boredom), names_to="trait", values_to="score") %>%
  arrange(id)

knitr::kable(head(AEQ,10)) # view data set
id group timepoint trait score
1 Twice 0 enjoyment 9
1 Twice 0 hope 12
1 Twice 0 pride 12
1 Twice 0 anger 20
1 Twice 0 anxiety 16
1 Twice 0 shame 11
1 Twice 0 hopelessness 9
1 Twice 0 boredom 9
1 Twice 6 enjoyment 12
1 Twice 6 hope 12
  • I haven’t include non-complying into this AEQ simulation.
  • The positive traits were simulated to have increasing scores and the negative traits were simulated to have decreasing scores.

Custom Measure (Expected Exam Score, Real Exam Score, total studying time, subjective study focus quality)

# expected
course.name <- sample(c("PSYC 110", "PSYC 111", "PSYC 205", "PSYC 210", "PSYC 222"), n, replace = TRUE) # several psychology courses provided at Oxford College
difficulty <- sample(1:5,n, replace=TRUE, prob = c(.05,.05,.1,.3,.5))

date <- sample(seq(as.Date('2021/10/01'), as.Date('2021/11/05'), by="day"), n, replace = TRUE)

score.exp <- rtruncnorm(n=n, a=60, b=100, mean=85, sd=8) # out of 100

# real score, generate by group

none = order.group %>%
    filter(group == "None") %>%
    summarise(length = n())

once = order.group %>%
    filter(group == "Once") %>%
    summarise(length = n())
twice = order.group %>%
    filter(group == "Twice") %>%
    summarise(length = n())

score.none <- rtruncnorm(n=none$length, a=60, b=100, mean=85, sd=8)
score.once <- rtruncnorm(n=once$length, a=60, b=100, mean=90, sd=8)
score.twice <- rtruncnorm(n=twice$length, a=60, b=100, mean=95, sd=8)
score <- c(score.none, score.once, score.twice) # out of 100


exam.score <- data.frame(id = order.group$id, group = order.group$group, course.name, difficulty, date, score.exp, score) %>%
    arrange(id)
fear <- round(rtruncnorm(n=n, a=1, b=5, mean=4.5, sd=1))
parent <- round(rtruncnorm(n=n, a=1, b=5, mean=3, sd=1))
peer <- round(rtruncnorm(n=n, a=1, b=5, mean=4, sd=1))
professor <- round(rtruncnorm(n=n, a=1, b=5, mean=4, sd=1))
profession <- round(rtruncnorm(n=n, a=1, b=5, mean=4.5, sd=1))
intelligence <- round(rtruncnorm(n=n, a=1, b=5, mean=4.5, sd=1))
effort <- round(rtruncnorm(n=n, a=1, b=5, mean=3, sd=1))

exam.fear <- data.frame(id, fear, parent, peer, professor, profession, intelligence, effort)
# same across 2 experimental groups
extra <- c(rep(NA,none$length), sample(seq(0.5,2, by = 1/6), size = n-none$length, replace = TRUE))
attention <- c(rep(NA,none$length),round(rtruncnorm(n=n-none$length, a=1, b=5, mean=4, sd=1)))

# differ by groups

perc.focus <- c(rep(NA,none$length),round(rtruncnorm(n=once$length, a=1, b=100, mean=80, sd=5)),
                round(rtruncnorm(n=twice$length, a=1, b=100, mean=87, sd=5))) # in percentage 
future <- c(rep(NA,none$length),sample(c(1:5), once$length, replace = TRUE, prob = c(.1,.2,.4,.2,.1)),
            sample(c(1:5), twice$length, replace = TRUE, prob = c(.1,.2,.2,.4,.1)))# 1 very Unlikely, 5 very likely
anxiety.sub <- c(rep(NA,none$length), sample(c("Yes", "No", "Not Sure"), once$length, replace = TRUE, prob = c(.5,.2,.3)),
                 sample(c("Yes", "No", "Not Sure"), twice$length, replace = TRUE, prob = c(.7,.1,.2)))

focus.sub <- c(rep(NA,none$length), sample(c("Yes", "No", "Not Sure"), once$length, replace = TRUE, prob = c(.4,.2,.4)),
                                           sample(c("Yes", "No", "Not Sure"), twice$length, replace = TRUE, prob = c(.7,.1,.2)))

final <- data.frame(id = order.group$id, group = order.group$group, extra, attention, perc.focus, future, anxiety.sub, focus.sub) %>%
    arrange(id)



# Clear Workspace
rm(list=setdiff(ls(), c("baseline","daily.mood","AEQ", "exam.score", "exam.fear", "final", "n", "id")))
# combine baseline, exam.score, exam fear survey, final survey data together for future use

combine <- baseline %>% full_join(exam.score) %>% full_join(exam.fear) %>% full_join(final)
## Joining, by = c("id", "group")
## Joining, by = "id"
## Joining, by = c("id", "group")
knitr::kable(head(combine, 10)) # preview data set
id group gender age grade race ethnicity major credit_hr course_hd med_freq med_exp med_type relax_q focus_q score_q course.name difficulty date score.exp score fear parent peer professor profession intelligence effort extra attention perc.focus future anxiety.sub focus.sub
1 Twice Male 20 Freshman Black Non_Hisp Pre-medical 17 Philosophy Never 10~100 Never 1 2 4 PSYC 110 3 2021-10-14 87.95563 92.22166 4 3 2 4 3 5 4 1.333333 4 82 3 Not Sure Not Sure
2 Twice Male 19 Sophomore White Non_Hisp Pre-medical 20 Biology Weekly 10~100 Unknown 1 4 2 PSYC 111 5 2021-10-24 90.29046 99.19380 4 4 2 5 4 4 2 1.166667 3 86 3 No Yes
3 Twice Male 19 Freshman Black Non_Hisp Pre-medical 16 Philosophy Never 100~1000 Unknown 5 2 4 PSYC 205 4 2021-10-09 84.47540 97.19650 4 2 3 4 5 4 3 1.333333 4 88 4 Yes Yes
4 Once Male 19 Senior Black Non_Hisp Psychology 15 Biology Never 1~10 Mindfulness 5 1 3 PSYC 111 4 2021-10-18 78.05886 82.76048 3 4 3 3 4 5 2 1.500000 4 77 1 Yes Yes
5 Twice Female 21 Junior Black Non_Hisp QSS 19 Chemistry Never 100~1000 Unknown 4 3 1 PSYC 222 5 2021-10-14 90.20686 78.30489 5 3 3 4 3 4 3 1.333333 4 91 2 Not Sure Not Sure
6 Once Female 19 Freshman White Non_Hisp Pre-medical 17 Physics Monthly 100~1000 Unknown 4 5 3 PSYC 205 4 2021-10-24 86.62935 88.18989 5 4 3 4 2 3 3 1.166667 4 83 2 Yes Yes
7 Once Male 19 Freshman White Non_Hisp Psychology 17 Biology Never 1000+ Mindfulness 4 3 4 PSYC 111 5 2021-10-13 75.10414 92.53632 4 4 4 2 3 4 2 1.500000 4 79 2 Yes Yes
8 Once Female 21 Senior White Non_Hisp Psychology 15 Philosophy Monthly 100~1000 Unknown 5 1 2 PSYC 210 2 2021-10-21 92.07442 87.39142 4 3 2 4 5 4 5 1.500000 4 82 3 Not Sure No
9 Twice Female 21 Senior White Non_Hisp Psychology 16 Philosophy Never 1~10 Unknown 4 2 5 PSYC 110 4 2021-10-27 94.02650 99.14277 4 1 4 4 5 4 3 0.500000 3 88 5 No Yes
10 None Female 18 Sophomore White Non_Hisp QSS 20 Math Never 100~1000 Never 1 4 2 PSYC 111 3 2021-10-25 90.57640 83.88125 4 1 4 3 4 3 2 NA NA NA NA NA NA

The data sets are combined according to the order: baseline, exam.score, exam.fear, final.

Hypothesis 1: FA meditation will reduce exam anxiety

Anxiety Score at 5 timepoints and Daily Mood Score from Day 0 to Day 21 (Animated Graph)

# Animated Graph
# error bar use standard error

anxiety.bar <- AEQ %>% 
    filter(trait == "anxiety") %>%
    left_join(select(baseline,id,group, gender)) %>%
    group_by(timepoint,group) %>%
    summarise(avg = mean(score,na.rm=T), n = n(), sd = sd(score,na.rm=T), se = sd/sqrt(n)) %>%
    plot_ly(x = ~group, y = ~avg, error_y = ~list(array = se,
                        color = '#000000'), 
            color = ~group, frame = ~timepoint, type = "bar", name = "FA Meditation X Anxiety") %>%
  layout(yaxis = list(range = c(0, Inf)), showlegend = FALSE) %>%
  animation_opts(500, easing = "linear", redraw = FALSE) %>%
  animation_slider(
    currentvalue = list(prefix = "Day ", font = list(color="darkblue"))
  )
## Joining, by = c("id", "group")
## `summarise()` has grouped output by 'timepoint'. You can override using the `.groups` argument.
data <- AEQ %>% 
    filter(trait == "anxiety") %>%
    left_join(select(baseline,id,group, gender)) %>%
    group_by(timepoint, group) %>%
    mutate(label.timepoint = as.factor(paste("Day", timepoint)))
## Joining, by = c("id", "group")
    # summarise(avg = mean(score,na.rm=T), n = n(), sd = sd(score,na.rm=T), se = sd/sqrt(n)) %>%

anxiety.dot <- ggplot(data, aes(x=group, y=score, fill = group, color = group
               )) + 
    geom_dotplot(alpha = 0.4, binaxis='y', stackdir='center', binwidth = 1, dotsize = 0.7) + 
    stat_summary(fun.data=mean_sdl, fun.args = list(mult=1), 
        geom="errorbar", color="black", width=0.2) +
  stat_summary(fun.y=mean, geom="point", color="black") + 
    theme_classic() + 
    facet_wrap(~label.timepoint) + 
    theme(legend.position = "none")
## Warning: `fun.y` is deprecated. Use `fun` instead.
anxiety.dot.animated <- ggplot(data, aes(x=group, y=score, fill = group, color = group
               )) + 
    geom_dotplot(alpha = 0.4, binaxis='y', stackdir='center', binwidth = 1, dotsize = 0.7) + 
    stat_summary(fun.data=mean_sdl, fun.args = list(mult=1), 
        geom="errorbar", color="black", width=0.2) +
    stat_summary(fun.y=mean, geom="point", color="black") + 
    # transition_time(as.numeric(timepoint)) + 
    transition_states(timepoint, transition_length = 1, state_length = 4) +
    ease_aes('cubic-in-out') +
    theme_classic() + 
    labs(title = 'Day: {closest_state}') + 
    # facet_wrap(~label.timepoint) + 
    theme(legend.position = "none")
## Warning: `fun.y` is deprecated. Use `fun` instead.
# Plot
anxiety.box <- data %>%
  mutate(text = paste("Score: ", score, "\nDay: ", timepoint)) %>%
  ggplot(aes(x=group, y=score, fill=group, frame = timepoint)) +
    # scale_fill_viridis(discrete = TRUE, alpha=0.6) +
    geom_jitter(color="black", size=0.4, alpha=0.5) +
    geom_boxplot(alpha=0.6) +
    theme_ipsum() +
    theme(
      legend.position="none",
      # plot.title = element_text(size=11)
    ) +
    ggtitle("FA Meditation X Anxiety boxplot with jitter") +
    xlab("")
anxiety.box.ani <- ggplotly(anxiety.box, tooltip = "text") %>%
    animation_opts(500, easing = "linear", redraw = FALSE) %>%
  animation_slider(
    currentvalue = list(prefix = "Day ", font = list(color="darkblue"))
  )
    
    
anxiety.bar
anxiety.dot

anxiety.dot.animated

anxiety.box.ani
rm(list=setdiff(ls(), c("baseline","daily.mood","AEQ", "exam.score", "exam.fear", "final", "n", "id", "combine")))

Four graphs are shown here. - Bar graph with error bars - Dot plot with error bars - Animated dot plot with error bars (one time point per frame) - Animated box plot with jitters (one time point per frame) - Double click the graph area when boxplot is not shown correctly

mood.bar <- daily.mood %>%
    group_by(timepoint,group) %>%
    rename(score = mood) %>%
    summarise(avg = mean(score,na.rm=T), n = n(), sd = sd(score,na.rm=T), se = sd/sqrt(n)) %>%
    plot_ly(x = ~group, y = ~avg, error_y = ~list(array = se,
                        color = '#000000'), color = ~group, frame = ~timepoint, type = "bar", name = "FA Meditation X Daily Mood (day 0 ~ 21)") %>%
  
  layout(yaxis = list(range = c(1, 5)), showlegend = FALSE) %>%
  animation_opts(500, easing = "linear", redraw = FALSE) %>%
  animation_slider(
    currentvalue = list(prefix = "Day ", font = list(color="darkblue"))
  )
## `summarise()` has grouped output by 'timepoint'. You can override using the `.groups` argument.
mood.bar
rm(list=setdiff(ls(), c("baseline","daily.mood","AEQ", "exam.score", "exam.fear", "final", "n", "id", "combine")))

ANOVA test

anova.mixed = function(data, dep.var = NULL){
# base dataset
data.name = deparse(substitute(data))
if(data.name == "daily.mood"){
    data = data %>%
        select(id, group, timepoint, dep.var) %>%
        rename(score = dep.var) %>%
        convert_as_factor(id, group, timepoint)
}

else if(data.name == "exam.score"){
    data = data %>%
        select(id, group, score.exp, score) %>%
        pivot_longer(c(score.exp, score), names_to = "timepoint", values_to = "score") %>%
        mutate(timepoint = ifelse(timepoint == "score.exp","expected","actual")) %>%
      mutate(timepoint = factor(timepoint, levels = c("expected", "actual")))
}

else{
    data = data %>% 
    filter(trait == dep.var) %>%
    select(-trait) %>%
    convert_as_factor(id, group, timepoint)
}

# View Summary Statistics
summary <- data %>%
    group_by(timepoint, group) %>%
    get_summary_stats(score, type = "mean_sd")
summary

bxp = NULL
if(n_distinct(as.numeric(data$timepoint))>5){
    bxp <- ggplot(data, mapping = aes(x = group, y = score)) + 
        geom_boxplot(aes(fill = group), alpha = 0.5) +
        # geom_jitter(width=0.1,alpha=0.2, aes(color = group)) + 
        # geom_dotplot(binaxis='y', stackdir='center', dotsize=1) + 
        facet_wrap(~timepoint)
}

else{
    bxp <- ggplot(data, mapping = aes(x = timepoint, y = score)) + 
        geom_boxplot(aes(fill = group), alpha = 0.5) 
        # geom_jitter(width=0.1,alpha=0.2, aes(color = group))
        # geom_dotplot(binaxis='y', stackdir='center', dotsize=1)
}
bxp <- bxp + geom_point(position=position_jitterdodge(0.2), mapping = aes(color = group)) + 
    # scale_color_grey() +
  scale_color_ft() + 
    theme_classic()


outliers <- data %>%
  group_by(timepoint, group) %>%
  identify_outliers(score)

## Normality assumption

# use this test if sample size <= 50
# anxiety %>%
#   group_by(timepoint, group) %>%
#   shapiro_test(score)

normality <- ggqqplot(data, "score", ggtheme = theme_bw()) +
  facet_grid(timepoint ~ group)
normality

## Homogneity of variance assumption p > 0.05
variance <- data %>%
  group_by(timepoint) %>%
  levene_test(score ~ group)
variance
## Homogeneity of covariances assumption
covariance <- box_m(data[, "score", drop = FALSE], data$group)$p.value > 0.001 # note the violation if p.value is false.
# In this simulation, p < 0.001 and the simulated data violate Homogeneity of covariances assumption.

covariance

# Anova test and Assumption of sphericity
# Two-way mixed ANOVA test

# automatically use type II test if each group has same number of participants, otherwise use type III test

res.aov <- anova_test(
  data = data, dv = score, wid = id,
  between = group, within = timepoint
  )
aov.tbl <- get_anova_table(res.aov)
aov.tbl

# Post-hoc Tests
## significant two-way interaction
# Simple main effect of group variable

# Effect of group at each time point
one.way <- data %>%
  group_by(timepoint) %>%
  anova_test(dv = score, wid = id, between = group) %>%
  get_anova_table() %>%
  adjust_pvalue(method = "bonferroni")
one.way
# Pairwise comparisons between group levels
pwc <- data %>%
  group_by(timepoint) %>%
  pairwise_t_test(score ~ group, p.adjust.method = "bonferroni")
pwc

# Simple main effects of time variable
# Effect of time at each level of exercises group
one.way2 <- data %>%
  group_by(group) %>%
  anova_test(dv = score, wid = id, within = timepoint) %>%
  get_anova_table() %>%
  adjust_pvalue(method = "bonferroni")
one.way2

# Pairwise comparisons between time points at each group levels
# Paired t-test is used because we have repeated measures by time
pwc2 <- data %>%
  group_by(group) %>%
  pairwise_t_test(
    score ~ timepoint, paired = TRUE, 
    p.adjust.method = "bonferroni"
    ) %>%
  select(-df, -statistic, -p) # Remove details
pwc2

# Procedure for non-significant two-way interaction
# time
non.pwc <- data %>%
    pairwise_t_test(
        score ~ timepoint, paired = TRUE, 
        p.adjust.method = "bonferroni"
    )

# group
non.pwc2 <- data %>%
    pairwise_t_test(
        score ~ group, 
        p.adjust.method = "bonferroni"
    )


# Visualization: boxplots with p-values
pwc.result <- pwc %>% add_xy_position(x = "timepoint")
pwc.filtered <- pwc.result %>% filter(timepoint != 0) # filter out timepoint 0

# if(length(distinct(pwc.result$timepoint))<=5){
    report <- bxp + 
    stat_pvalue_manual(pwc.filtered, label = "p.adj.signif", tip.length = 0, hide.ns = TRUE) + # show stars instead of actual p-values
    labs(
        subtitle = get_test_label(res.aov, detailed = TRUE),
        caption = get_pwc_label(pwc.result)
    )




print(report)

# assumptions <- list(summary, bxp, outliers, normality, variance, covariance)
stats.aov <- list(summary, bxp, outliers, normality, variance, covariance, 
                  aov.tbl, one.way, pwc, one.way2 ,pwc2, non.pwc, non.pwc2, report)

names(stats.aov) = c("summary", "bxp", "outliers", "normality", "variance", "covariance", 
                         "aov.tbl", "one.way", "pwc", "one.way2", "pwc2", "non.pwc", "non.pwc2", "report")
return(stats.aov)
}
anova.anxiety <- anova.mixed(AEQ, dep.var = "anxiety")
## Coefficient covariances computed by hccm()
## Coefficient covariances computed by hccm()
## Coefficient covariances computed by hccm()
## Coefficient covariances computed by hccm()
## Coefficient covariances computed by hccm()

# Automatically print the final boxplot for publication with p-value Asterisk (star)
# All data are saved for later use
anova.anxiety[[14]]

anova.mood <- anova.mixed(daily.mood, dep.var = "mood")
## Note: Using an external vector in selections is ambiguous.
## i Use `all_of(dep.var)` instead of `dep.var` to silence this message.
## i See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
## This message is displayed once per session.
## Coefficient covariances computed by hccm()
## Coefficient covariances computed by hccm()
## Coefficient covariances computed by hccm()
## Coefficient covariances computed by hccm()
## Coefficient covariances computed by hccm()
## Coefficient covariances computed by hccm()
## Coefficient covariances computed by hccm()
## Coefficient covariances computed by hccm()
## Coefficient covariances computed by hccm()
## Coefficient covariances computed by hccm()
## Coefficient covariances computed by hccm()
## Coefficient covariances computed by hccm()
## Coefficient covariances computed by hccm()
## Coefficient covariances computed by hccm()
## Coefficient covariances computed by hccm()
## Coefficient covariances computed by hccm()
## Coefficient covariances computed by hccm()
## Coefficient covariances computed by hccm()
## Coefficient covariances computed by hccm()
## Coefficient covariances computed by hccm()
## Coefficient covariances computed by hccm()
## Coefficient covariances computed by hccm()

rm(list=setdiff(ls(), c("baseline","daily.mood","AEQ", "exam.score", "exam.fear", "final", "n", "id", "combine", "anova.mixed")))

Hypothesis 2

anova.score <- anova.mixed(exam.score)
## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.

## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.
## Coefficient covariances computed by hccm()
## Coefficient covariances computed by hccm()

rm(list=setdiff(ls(), c("baseline","daily.mood","AEQ", "exam.score", "exam.fear", "final", "n", "id", "combine", "anova.mixed")))
# confidence

Hypothesis 3

Hypothesis 4