Sun. Dec. 18, 2022

STEP 1: Test ANOVA assumptions

ANOVA Assumptions (in order of importance):

Independence - Data are independent (not sampled from the same individual, etc.)
Normality - The response variable (and any numeric factors) has a normal distribution
Homoscedasticity - The sample variance across factors is similar

# Test normality; transform if necessary

# Define dataset
dat_stat <- dat_stat

# Assign factors
dat_stat$treatment  <- factor(dat_stat$treatment, levels = c("control","SS","MS"))
dat_stat$ploidy     <- factor(dat_stat$ploidy)
dat_stat$timepoint  <- factor(dat_stat$timepoint, levels = c("0","11","12","15","20"), ordered = "TRUE")

# Assign response variable
test_me <- dat_stat$ATPase

# Test for normality
qqnorm(test_me, main = "Q-Q Plot: untransformed") # check linearity
qqline(test_me)
norm_test <- shapiro.test(test_me) # p-value > 0.05 = good, don't need transformation
print(paste("shapiro test p-value, untransformed:", norm_test$p.value))

# Normalize response variable if normality test failed
if(norm_test$p.value<0.05)     {
        normalized <- bestNormalize(test_me, main = "Q-Q Plot: transformed")
        test_me <- normalized$x.t # overwrite
        qqnorm(test_me) # check linearity of transformed response
        qqline(test_me)
        norm_test <- shapiro.test(test_me) # p-value > 0.05 = good
        print(paste("shapiro test p-value, transformed:", norm_test$p.value))}
dat_stat$response <- test_me # overwrite

untransformed	transformed

STEP 2: Run ANOVA

# Run ANOVA
my_test <- aov(response ~ ploidy * timepoint * treatment, data = dat_stat)
my_test_summary <- summary(my_test)
summary(my_test)

STEP 3: compare model AIC scores

# Compare model AIC scores (lowest score wins)
other <- aov(response ~ ploidy * timepoint, data = dat_stat)

model.set <- list(my_test, other)
model.names <- c("ploidy:timepoint:treatment", "ploidy:timepoint")

aictab(model.set, modnames = model.names)

STEP 4: Check for Homoscedasticity across factors in model

leveneTest(response ~ interaction(ploidy,timepoint,treatment), dat_stat)

STEP 5: Run post-hoc test if interaction is significant

tx <- with(dat_stat, interaction(timepoint,treatment,ploidy)) # build interaction
amod <- aov(response ~ tx, data = dat_stat) # run model
mult_comp <- HSD.test(amod, "tx", group=TRUE, console=TRUE) # run HSD test