tab <- read.table("ozone_fiche.txt",sep=" ")


tab$vent <- as.factor(tab$vent)
tab$pluie <- as.factor(tab$pluie)

nn <- table(tab$vent)
nleg <- paste(paste("n",1:4,sep="")," = ",nn,sep="")
png("../fig/ozone_boxplot.png")
  par(mar=c(5,5,5,2))
  boxplot(maxO3~vent,data=tab,cex.lab=1.6,
          xlab="Wind direction",ylab="Ozone level")
  mtext(side=3,nleg, at=1:4,line=1,cex=1.5)
dev.off()

mylm <- lm(maxO3~vent,data=tab)
anova(mylm)

par(mfrow=c(1,2))
plot(mylm,which=c(1,2))

png("../fig/ozone_residuals.png",width=920,height=460)
  par(mfrow=c(1,2))
  plot(mylm,which=c(1,2))
dev.off()


summary(mylm)

library(emmeans)
myemeans <- emmeans(mylm, ~ vent)
png("../fig/confidence_interval.png")
  plot(myemeans)
dev.off()  

cc <- contrast(Emeans, method="pairwise", adjust="none")
png("../fig/pairwise_tests.png")
  plot(cc)
dev.off()

contr <- list("NW vs SE" = c(-0.5, 0.5, 0.5, -0.5))
contrast(myemeans, contr, adjust="none")