library(Stat2Data)
library(ggplot2) 
library(MASS)
library(agricolae)
library(dplyr)
#------------------------------------------------------------------------

#Example 1: Is alfalfa height different across acid treatment groups?

#Get data
data("Alfalfa")
head(Alfalfa)

#Test Assumptions

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

#Normality
hist(Alfalfa$Ht4, col="yellow") #good enough, too few observations by group to show

#Equal Variance
boxplot(Ht4~Acid, data=Alfalfa, col="yellow")
round(tapply(Alfalfa$Ht4, Alfalfa$Acid, var),2)
1.73/0.19 #9.11, violated, but continue for now

#ANOVA
AOV1.1 <-aov(Ht4~Acid, data=Alfalfa)
summary(AOV1.1) #significant, but assumption of equal variance violated

posthoc1 <-TukeyHSD(AOV1.1, "Acid")
posthoc1
plot(posthoc1) #water and 3.0 HCI just barely different

#Kruskal Wallic 
AOV1.2 <-kruskal.test(Ht4~Acid, data=Alfalfa)
AOV1.2 #not significant

#Graph
Alfalfa_sum<-Alfalfa %>% 
  group_by(Acid) %>%       
  summarise(mean=mean(Ht4),
            sd = sd(Ht4),
            error = qt(0.975,df=n()-1)*sd/sqrt(n()),
            ul = mean + error,
            ll = mean - error)
f1 <-ggplot(data=Alfalfa_sum, aes(x=Acid, y=mean, fill='red')) +
  geom_bar(stat="identity") +
  geom_errorbar(aes(ymin=ll, ymax=ul), width=0.1)+
  theme(legend.position="none")+
  labs(y="Mean height")
f1
#------------------------------------------------------------------------


#Example 2: Is ant number different across sandwich bread type?

#Get data
data("SandwichAnts")
head(SandwichAnts)

#Test Assumptions

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

#Normality
hist(SandwichAnts$Ants, col="red") #initially good
hist(SandwichAnts$Ants[SandwichAnts$Bread=="MultiGrain"], col="brown")
hist(SandwichAnts$Ants[SandwichAnts$Bread=="WholeWheat"], col="brown")
hist(SandwichAnts$Ants[SandwichAnts$Bread=="White"], col="brown")
hist(SandwichAnts$Ants[SandwichAnts$Bread=="Rye"], col="brown") #all good enough

#Equal Variance
boxplot(Ants~Bread, data=SandwichAnts, col="red")
round(tapply(SandwichAnts$Ants, SandwichAnts$Bread, var),2)
332.36/179.48 #1.85, good

#ANOVA
AOV2 <-aov(Ants~Bread, data=SandwichAnts)
summary(AOV2) #not significant

#Graph
SandwichAnts_sum<-SandwichAnts %>% 
  group_by(Bread) %>%       
  summarise(mean=mean(Ants),
            sd = sd(Ants),
            error = qt(0.975,df=n()-1)*sd/sqrt(n()),
            ul = mean + error,
            ll = mean - error)
f2 <-ggplot(data=SandwichAnts_sum, aes(x=Bread, y=mean, fill='red')) +
  geom_bar(stat="identity") +
  geom_errorbar(aes(ymin=ll, ymax=ul), width=0.1)+
  theme(legend.position="none")+
  labs(y="Mean Ant number")
f2
#------------------------------------------------------------------------


#Example 3: Is weight in chicks different across feed type?

#Get data
head(chickwts)

#Test Assumptions

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

#Normality
hist(chickwts$weight, col="orange") #good enough
ggplot(data=chickwts) +
  geom_histogram(aes(x=weight), bins=5, colour="black", fill="orange") +
  facet_wrap(~feed) #good enough

#Equal Variance
boxplot(weight~feed,data=chickwts, col='orange')
chickwts_var<-chickwts %>% 
  group_by(feed) %>%       
  summarise(var=var(weight))
chickwts_var
4212/1492 #2.82, good enough

#ANOVA
AOV3 <-aov(weight~feed,data=chickwts)
summary(AOV3) #feed is significant

posthoc3 <-TukeyHSD(AOV3, "feed")
posthoc3
plot(posthoc3)

HSD.test(AOV3, "feed", console=T) #need to switch over letters to organize
#b to c and c to b

#Graph
chickwts_sum<-chickwts %>% 
  group_by(feed) %>%       
  summarise(mean=mean(weight),
            sd = sd(weight),
            error = qt(0.975,df=n()-1)*sd/sqrt(n()),
            ul = mean + error,
            ll = mean - error)
f3 <-ggplot(data=chickwts_sum, aes(x=feed, y=mean, fill=feed)) +
  geom_bar(stat="identity") +
  geom_errorbar(aes(ymin=ll, ymax=ul), width=0.1)+
  scale_fill_manual(values=c("goldenrod", "goldenrod3", "khaki2","khaki3", "khaki4", "gold")) +
  theme(legend.position="none")+
  labs(y="mean chick weight", x="feed type")+
  scale_y_continuous(limits=c(0,400)) +
  geom_text(aes(label=c("a","b","bc","ac","c","a"), group=feed), vjust=-3, size=6, position=position_dodge(0.9))
f3