library(help="datasets")	# ottenere informazioni su dataset di esempio
?airquality			# informazioni su dataset air quality
air.cleaned <- na.omit(airquality)	# filtro righe contenenti valori NA


par(mfrow=c(3,2))		# impostazione finestra output grafico di R (3 righe x 2 colonne : 6 grafici)

# Multiple linear regression : fitting del modello
# Effettua una regressione lineare (multipla) della concentrazione di ozono in funzione di radiazione solare, vento e temperatura. Viene utilizzata la funzione lm(). Questa funzione utilizza il metodo 'Ordinary least squares' (OLS). OLS stima i parametri sconosciuti di una funzione lineare di un set di variabili esplicative (regressori) minimizzando la somma dei quadrati delle differenze tra i valori predetti dal modello e i valori della variabile dipendente (quella che vogliamo predire utilizzando il modello) nel dataset in esame.    

model1 <- lm(air.cleaned$Ozone ~ air.cleaned$Solar.R + air.cleaned$Wind + air.cleaned$Temp, data=air.cleaned)
summary(model1)			# stampa informazioni sul modello ottenuto

plot(model1, which=c(1:4))	# stampa plot diagnostici a partire al modello (di default ne crea 4)
				# specifichiamo quali posizioni vogliamo che occupino i grafici nella
				# griglia 3 x 2 definita in precedenza (posizioni da 1 a 4)
				# abbiamo ancora 2 slot liberi ...
hist(residuals(model1), breaks=10)	# grafico istogramma dei residui
plot(predict(model1) ~ air.cleaned$Ozone)	# plot valori predetti dal modello vs valori osservati
abline(0,1)


# Ipotizziamo che ci sia una interazione tra i valori di temperatura e radiazione solare nel determinare la concentrazione di ozono. Per testare velocemente se i nostri dati supportan questa ipotesi creiamo uno scatterplot di radiazione solare vs temperatura in cui la dimensione dei simboli e' proporzionale ai valori osservati di concentrazione di ozono. Dopo aver esaminato il grafico che conclusioni potete trarre? (Motivate la risposta)

symbols(air.cleaned$Temp, air.cleaned$Solar.R, circles=air.cleaned$Ozone/100, ylab="Solar Radiation", xlab="Temperature", main="Interaction Plot", inches=FALSE)

formula2 <- "Ozone ~ Wind + Solar.R * Temp"	# Creiamo una nuova formula che tenga conto dell'ipotetica interazione 
						# tra temperatura e radiazione solare nell'influenzare concentrazione ozono
model2 <- lm(formula2 , data=air.cleaned)	# Usiamo la formula (che e' salvata in una stringa) per fittare un modello

# Confrontiamo le informazioni di base dei modelli ottenuti ignorando e tenendo conto dell'ipotetica interazione Temp. Rad. solare 
summary(model1)
summary(model2)

# Come posso confrontare i modelli ottenuti? 
# Potrei pensare che e' sufficiante esaminare il valore di R^2 ... purtroppo esso aumenta sempre in seguito all'aggiunta di parametri e, di conseguenza, questa non e' una strategia adeguata. In questi casi e' possibile effettuare un F test (o, se preferite, un Likelihood Ratio Test, LRT) per testare in modo formale se il modello piu' complesso (quello che include il termine di interazione) descrive i dati osservati significativamente meglio del modello che ne e' sprovvisto.  

anova(model1, model2, test="F")