###############################################
#--- 2 fattori, ogni fattore 2 livelli
#--- in realtà sarebbe sufficiente calcolare F-value
#--- perchè i livelli sono meno di 3, ma la funzione anova (o aov)
#--- di R dà un risultato corretto anche in questo caso


# Dati:

F1<-c(1,1,1,1,-1,-1,-1,-1)
F2<-c(1,1,-1,-1,1,1,-1,-1)
R<-c(185,193,46,66,96,83,132,104)

Dati<-data.frame(F1=F1,F2=F2,R=R)


######
#-- Senza interazioni

Model_LM<-lm(R~F1+F2,data=Dati)

coefficients(Model_LM)

anova(Model_LM)

# Oppure, in alternativa, si esegue:

AOV<-aov(R~F1+F2,data=Dati)

summary(AOV)

######
#-- Con interazioni (sono significative! Osservare i p-value)

AOV_int<-aov(R~F1*F2,data=Dati)

summary(AOV_int)

# Per ottenere i coefficienti del modello di regressione:

Model_LM_int<-lm(R~F1*F2,data=Dati)

coefficients(Model_LM_int)


########################################################################################
### MULTI-WAY ANOVA


#http://openmv.net/info/blender-efficiency

#Blender efficiency
#Description:	The effect of 4 factors on blending efficiency.
#Data source:	

#    particle size
#    mixer diameter
#    mixer rotational speed
#    blending time

#Modified based on Box, Hunter, and Hunter, Statistics for Experimenters, 2nd edition, page 486.
#Data shape:	18 rows and 5 columns

originale<-read.table("Blending Efficiency.txt",sep=",",header=T)
colnames(originale)<-c("PSize","MixD","MixR","BlenT","Response")

str(originale)

#--Trasformo tutte le variabili indipendenti in vettori di tipo factor (in modo che R "veda" i livelli)

originale$PSize<-as.factor(originale$PSize)
originale$MixD<-as.factor(originale$MixD)
originale$MixR<-as.factor(originale$MixR)
originale$BlenT<-as.factor(originale$BlenT)

str(originale)


### TEST PRELIMINARI PER STABILIRE APPLICABILITA' DI ANOVA

#--Test normalità della risposta
#--L'ipotesi nulla è che la distribuzione dei valori sia normale (gaussiana)
#--quindi se p-value > 0.1 la distribuzione è normale

shapiro.test(originale$Response) #OK

#Approssimativamente si riesce a capire dalla visulaizzazione dell'istrogramma:

hist(originale$Response)


#--- Test omogeneità
#--L'ipotesi nulla è che la distribuzione dei valori degli errori sia omogenea per ogni livello
#--quindi se p-value > 0.1 c'è omogeneità (omoschedasticità)


bartlett.test(Response~PSize,data=originale) #OK
bartlett.test(Response~MixD,data=originale) #OK
bartlett.test(Response~MixR,data=originale) #OK (p-value=0.095, al limite)
bartlett.test(Response~BlenT,data=originale) #OK


### MULTI-WAY ANOVA

#-- ANOVA con interazioni

Blend_aov_int<-aov(Response~PSize*MixD*MixR*BlenT,data=originale)

summary(Blend_aov_int)


#-- ANOVA senza interazioni

Blend_aov<-aov(Response~PSize+MixD+MixR+BlenT,data=originale)

summary(Blend_aov)


### POST-HOC TEST 
# nel modello di "ANOVA senza interazioni" considero le due variabili che hanno p-value <0.05
# Per capire quali sono i livelli che danno valori più diversi rispetto alla Risposta
# si guardano i p-value (qui si chiamano p-adj)
# Nel grafico si guardano le coppie di livelli che si allontanano di più dallo zero

TukeyHSD(Blend_aov, "BlenT", conf.level = 0.95)
plot(TukeyHSD(Blend_aov, "BlenT", conf.level = 0.95))
# Per il fattore BlenT le Risposte più lontane tra di loro sono per i livelli 45 vs. 75

TukeyHSD(Blend_aov, "PSize", conf.level = 0.95)
plot(TukeyHSD(Blend_aov, "PSize", conf.level = 0.95))
# Per il fattore PSize, osservando sia i p-value (bassi)che il grafico (distanza dallo zero) si osserva che i valori della Risposta
# per ogni livello sono nettamante separati dagli altri livelli