sites<-read.csv("bpa_sample_sites.csv")

#Sample size
str(sites)

#Normality assumptions
hist(sites$concentration, breaks=20) #normally distributed
qqnorm (sites$concentration)
qqline(sites$concentration, col = "red", lwd=3)
ks.test(sites$concentration, "pnorm", mean(sites$concentration), sd (sites$concentration))
shapiro.test(sites$concentration)

for (i in unique(sites$site)){
  par(mfcol = c(1,2))
  hist(sites$concentration[which(sites$site==i)], breaks=20) 
  qqnorm (sites$concentration[which(sites$site==i)])
  qqline(sites$concentration[which(sites$site==i)], col = "red", lwd=3)
  print(ks.test(sites$concentration[which(sites$site==i)], "pnorm", mean(sites$concentration[which(sites$site==i)]), sd (sites$concentration[which(sites$site==i)])))
  print(shapiro.test(sites$concentration[which(sites$site==i)]))
 
}

par(mfcol = c(1,1))

#Variance assumptions
boxplot(sites$concentration ~ sites$site)

#Anova itself
anova<-aov(sites$concentration ~ sites$site)

#Details of the model
summary(anova)

#Post-hoc
TukeyHSD(anova)


#II part (optional)


life_exp<-read.csv("country_lifeexp.csv")

#sample size
str(life_exp)

#normality assumptions 
hist(life_exp$life_expectancy, breaks=20) #is not normally distributed
qqnorm (life_exp$life_expectancy)
qqline(life_exp$life_expectancy, col = "red", lwd=3)
ks.test(life_exp$life_expectancy, "pnorm", mean(life_exp$life_expectancy), sd (life_exp$life_expectancy))
shapiro.test(life_exp$life_expectancy)

for (i in unique(life_exp$country)){
  par(mfcol = c(1,2))
  hist(life_exp$life_expectancy[which(life_exp$country==i)], breaks=20) 
  qqnorm (life_exp$life_expectancy[which(life_exp$country==i)])
  qqline(life_exp$life_expectancy[which(life_exp$country==i)], col = "red", lwd=3)
  print(ks.test(life_exp$life_expectancy[which(life_exp$country==i)], "pnorm", mean(life_exp$life_expectancy[which(life_exp$country==i)]), sd (life_exp$life_expectancy[which(life_exp$country==i)])))
  print(shapiro.test(life_exp$life_expectancy[which(life_exp$country==i)]))
  
}

par(mfcol = c(1,1))

#distribution shapes
par(mfcol = c(1,3))
for (i in unique(life_exp$country)){
  hist(life_exp$life_expectancy[which(life_exp$country==i)], breaks=20) 
}
par(mfcol = c(1,1))

#kruskal-wallis
kruskal.test(life_exp$life_expectancy ~ life_exp$country)

#Dunn's test

install.packages("dunn.test")
library(dunn.test)
dunn.test(life_exp$life_expectancy, life_exp$country, altp=TRUE)
boxplot(life_exp$life_expectancy ~ life_exp$country)