library(Stat2Data)
library(ggplot2) 
library(dplyr)
library(nlme)
library(emmeans)
library(lme4)
library(agricolae)
#------------------------------------------------------------------------


#Example 1: Is there a within-factors difference for counting accuracy for the condition of music type?

#Get Data
data("MusicTime")
head(MusicTime)

#Test Assumptions

#Independence
#depends on sample design and taken to be valid here

#Normality
ggplot(data=MusicTime, aes(x=Accuracy, fill=Music))+
  geom_histogram(color=1, binwidth = 10)+
  facet_wrap(~Music)
ggplot(data=MusicTime, aes(sample=Accuracy, color=Music))+
  geom_qq(size=2) + geom_qq_line(size=2) +
  facet_grid(~Music)
#normality close enough

#Equal Variance
boxplot(Accuracy~Music, data=MusicTime, col="blue")
round(tapply(MusicTime$Accuracy, MusicTime$Music, var),2)
229.12/67.88 #3.37, okay

#Repeated Measures ANOVA
RM_AOV <-aov(Accuracy ~ Music + Error(Subject), data=MusicTime)
summary(RM_AOV)

#Post-Hoc Tests
posthoc.1<-emmeans(RM_AOV, ~ Music)
pairs(posthoc.1) #calm sign. different from control and upbeat

#Graph
MusicTime_sum<-MusicTime %>%
  group_by(Music) %>%
  summarise(mean=mean(Accuracy),
            sd = sd(Accuracy),
            error = qt(0.975,df=n()-1)*sd/sqrt(n()),
            ul = mean + error,
            ll = mean - error)
p1 <-ggplot(data=MusicTime_sum, aes(x=Music, y=mean)) +
  geom_bar(aes(fill=Music, color=1), stat="identity", position=position_dodge()) +
  geom_errorbar(aes(ymin=ll, ymax=ul, group=Music), width=0.1, position=position_dodge(0.9))+
  scale_fill_manual(values=c("lightblue","grey","orange"))+
  geom_text(aes(label="*"), x=c(1), y=c(29), size=7) 
p1
#------------------------------------------------------------------------


#Example 2: Is there a difference in pea yield across nitrogen (N), phosphate (P), or potassium (K)
#           application, while controlling for block?

#Get Data
head(npk)

#Test Assumptions

#Independence
#depends on sample design and taken to be valid here

#Normality
ggplot(data=npk, aes(x=yield, fill=N))+
  geom_histogram(color=1, binwidth = 10)+
  facet_wrap(~N)
ggplot(data=npk, aes(x=yield, fill=P))+
  geom_histogram(color=1, binwidth = 10)+
  facet_wrap(~P)
ggplot(data=npk, aes(x=yield, fill=K))+
  geom_histogram(color=1, binwidth = 10)+
  facet_wrap(~K)
#normality good

#Equal Variance
boxplot(yield~N + P + K, data=npk)
npk_var<-npk %>%
  group_by(N,P,K) %>%
  summarise(var=var(yield))
npk_var
88.6/5.59 #15.85, variance violated, compared between blocked or not

#Blocked ANOVA (no blocking)
Bl_AOV.1 <-lm(yield ~ factor(N) + factor(P) + factor(K), data=npk)
summary(Bl_AOV.1) #N significant

#Blocked ANOVA (with fixed blocking)
Bl_AOV.2 <-lm(yield ~ factor(N) + factor(P) + factor(K) + factor(block), data=npk)
summary(Bl_AOV.2) #N and K significant

#Blocked ANOVA (with random blocking)
Bl_AOV.3 <-lmer(yield ~ factor(N) + factor(P) + factor(K)+ (1|block), data=npk)
summary(Bl_AOV.3) #N and K significant


#Post-Hoc Tests
#not needed for 2 groups
boxplot(yield~N, data=npk) #higher yield with Nitrogen (N=1)
boxplot(yield~K, data=npk) #higher yield without Potassium (K=0)


#Graphs
npk_Nsum<- npk %>%
  group_by(N) %>%
  summarise(mean=mean(yield),
            sd = sd(yield),
            error = qt(0.975,df=n()-1)*sd/sqrt(n()),
            ul = mean + error,
            ll = mean - error)
p2a <-ggplot(data=npk_Nsum, aes(x=N, y=mean)) +
  geom_bar(aes(fill=N, color=1), stat="identity", position=position_dodge()) +
  geom_errorbar(aes(ymin=ll, ymax=ul, group=N), width=0.1, position=position_dodge(0.9))+
  scale_fill_manual(values=c("grey","blue"))
p2a

npk_Ksum<- npk %>%
  group_by(K) %>%
  summarise(mean=mean(yield),
            sd = sd(yield),
            error = qt(0.975,df=n()-1)*sd/sqrt(n()),
            ul = mean + error,
            ll = mean - error)
p2b <-ggplot(data=npk_Ksum, aes(x=K, y=mean)) +
  geom_bar(aes(fill=K, color=1), stat="identity", position=position_dodge()) +
  geom_errorbar(aes(ymin=ll, ymax=ul, group=K), width=0.1, position=position_dodge(0.9))+
  scale_fill_manual(values=c("grey","purple"))
p2b
#------------------------------------------------------------------------


#Example 3: Is there a difference in log-optical density across dilutions, 
#           where dilutions are nested by sample?

#Get Data
data("Assay")
head(Assay)

#Test Assumptions

#Independence
#dilutions are pseudo-replicates, so not independent

#Normality
ggplot(data=Assay, aes(x=logDens))+
  geom_histogram(color=1, binwidth = 1)+
  facet_wrap(~dilut) #too few numbers to separate out
ggplot(data=Assay, aes(x=logDens))+
  geom_histogram(color=1, bins=10)
#not normally distributed, but already log transformed, so go as is

#Equal Variance
boxplot(logDens~dilut, data=Assay)
Assay_var<-Assay %>%
  group_by(dilut) %>%
  summarise(var=var(logDens))
Assay_var
0.00992/0.00282 #3.52, variance good

#Nested ANOVA
Ne_aov <-lmer(logDens~factor(dilut) +0 +(1|sample), data=Assay)
summary(Ne_aov) #significant

#Post-Hoc Tests
posthoc.3<-emmeans(Ne_aov, ~ dilut)
pairs(posthoc.3) #all significantly different

#Graph
Assay_sum<- Assay %>%
  group_by(dilut) %>%
  summarise(mean=mean(logDens),
            sd = sd(logDens),
            error = qt(0.975,df=n()-1)*sd/sqrt(n()),
            ul = mean + error,
            ll = mean - error)
p3 <-ggplot(data=Assay_sum, aes(x=dilut, y=mean)) +
  geom_bar(aes(fill=dilut, col=1), stat="identity", position=position_dodge()) +
  geom_errorbar(aes(ymin=ll, ymax=ul, group=dilut), width=0.1, position=position_dodge(0.9))+
  scale_fill_manual(values=c("azure4","azure3","azure2","azure1","azure"))
p3