rm(list=ls())
d.df <- read.table("exp1.csv", h=T, sep=";", dec=",")
d.df
## il comando seguente è un errore: si eliminano dei dati
## meglio non assegnare il risultato
na.omit(d.df)

d2.df <- data.frame(weight=as.vector(as.matrix(d.df[,1:3])))
d2.df
names(d.df)
d2.df$group <- rep(names(d.df), each=10)
d2.df
rep(c("M","F"),c(5,5) )
rep(rep(c("M","F"),c(5,5) ),3)
d2.df$sex<-rep(rep(c("M","F"),c(5,5) ),3)
d2.df
##a questo punto invece na.omit è corretto
d2.df <- na.omit(d2.df)
medie <- with(d2.df,(tapply(weight,list(group, sex), mean)))
medie

##funzione per il calcolo dell'errore standard
se <- function(x) sd(x)/sqrt(length(x))
ses <- with(d2.df,(tapply(weight,list(group, sex), se)))
ses

x <- 1:3
colori <- c("red", "blue")
names(colori) <-colnames(medie)
colori
max(medie + 2*ses)
min(medie - 2*ses)

plot(as.vector(medie) ~ rep(x,2), type="n", xaxt="n", ylab="Weight (g)", xlab="Treatment", ylim=c(2,8),xlim=c(0.5,3.5))
for(i in 1:ncol(medie)) {
    lines(medie[,i] ~c(x+(i-1)*0.04), col=colori[i])
    arrows(x+(i-1)*0.04, medie[,i]-2*ses[,i], x+(i-1)*0.04, medie[,i]+2*ses[,i], angle=90, code=3, length=0.05)
    points(medie[,i] ~c(x+(i-1)*0.04), pch=21, bg=colori[i], cex=2)
}
axis(1, at=x+0.02, lab=rownames(medie))
legend("topleft", pch=21, pt.bg=colori, legend=names(colori), pt.cex=2, title="Sex")


## MODO ALTERNATIVO USANDO dplyr
library(dplyr)
## questa può sembrare una sintassi strana
## ma è tipica dei pacchetti più nuovi di R ed è sostanzialmente
## semplice e concisa.
## La spiegheremo a lezione
## Il risultato è una oggetto leggermente diverso da un data.frame
## chiamato tibble.
d2.df %>% group_by(sex, group)%>% summarise(media=mean(weight),
                                            errstd=se(weight),
                                            enne=n()
                                            )
## qui converto la tibble in un data.frame normale. 
rd2.df <-as.data.frame(d2.df %>% group_by(sex, group)%>% summarise(
                                                             media=mean(weight),
                                                             errstd=se(weight),
                                                             enne=n()))
rd2.df
##il risultato ha una struttura diversa dalla matrice
##ma simile al dataframe d2.df
xc <-barplot(rd2.df$media, col=rep(colori,each=3), ylim=c(0,8),
             ylab="Weight (g)")
axis(1,at=xc,labels=paste(rd2.df$sex, rd2.df$group,sep="-"), las=3)
arrows(xc, rd2.df$media, xc, rd2.df$media+2*rd2.df$errstd,angle=90,
       code=3, length=0.04)

rd3.df <-as.data.frame(d2.df %>% group_by(group,sex)%>% summarise(
                                                            media=mean(weight),
                                                            errstd=se(weight),
                                                            enne=n()))
rd3.df
xc <-barplot(rd2.df$media, col=rep(colori,3), ylim=c(0,8), ylab="Weight (g)")
axis(1,at=xc,labels=paste(rd3.df$sex, rd3.df$group,sep="-"), las=3)
arrows(xc, rd2.df$media, xc, rd2.df$media+2*rd2.df$errstd,angle=90,
       code=3, length=0.04)
legend("topleft", fill=colori, legend=names(colori), title="Sex")