####    Regressione Lineare Semplice    ###
########    prof. Mini Valentina    #######
#    appunti per STATISTICA MULTIVARIATA  #
##############    2019-2020   #############

#########creare un database semplice in R: 

X=c(1213,1518,3050,852,1550,1215,2120,2207,2175)
Y=c(13474,16497,29349,11314,17224,14459,22186,23483,24095)

#dove X= numero di pezzi di laminati prodotti (dimensione del processo di produzione)

#Y= costo di produzione in serie di laminati

#########caricare dall'esterno un database semplice in R
#cambio directory e indico quella di riferimento

dati=read.csv2("produzione.csv",header=T)

#facciamo riferimento al database semplice "produzione" da Pearson

dati #per visualizzare il dataset "dati"

dim(dati) #per visualizzare il numero di unità e di var

n=dim(dati)[1] #n dimensione del campione

names(dati) #nome delle variabili

attach(dati) #siamo sicuri che in R ci sia il dataset 


############VISUALIZZAZIONE GRAFICA:SCATTERPLOT

#può essere utile visualizzare graficamente i dati 
#mediante scatterplot
#tra costo di produzione e numero di pezzi prodotti 

#scatterplot di Y vs X

plot(X,Y,main="scatterplot di Y vs X", lwd=2, xlab="Numero pezzi prodotti", ylab="Costo di produzione")


########## MODELLO DI REGRESSIONE LINEARE SEMPLICE 

#ora richiamiamo il comando per creare un "linear Model" lm (ovvero una regressione lineare semplice)
#l'output è vasto e complesso, quindi suggerisco di salvare l'analisi 
#in una variabile per poi estrarre le informazioni necessarie per analizzare i dati

#output del modello di regressione
#si usa il simbolo tilde. Generalmente: 
#tilde ~ in windows è alt+126
#tilde ~ in MC è alt+5


result=lm(Y~X)
result

#la variabile result è una variabile creata da lm ed è un oggetto di classe lm
#in tale classe sono contenuti: 
# a) la formula in mase alla quale l'oggetto è stimato
# b) le stime puntuali dei coefficienti di regressione
# c) una lista conetenete coefficients, residuals, fitted.values, rank, weights, df.residuals, call ...

#tutte queste informazioni le useremo nell'analisi di lm, per avere info aggiuntive help(lm)
#per "estrarre" le informazioni si utilizza la specifica $

########SCATTERPLOT CON MODELLO STIMATO#############

#per sovrapporre allo scatterplot il nostro modello stimato è necessario
#estrarre i valori stimati di Y in corrispondenza delle osservazioni della var. esplicativa X

#valori stimati per la variabile Y
Ycapp=result$fitted

#ovvero estraiamo dall'analisi contenuta nella variabile "result" 
# i valori del modello stimato (o fitted model)


plot(X,Y,main="scatterplot di Y vs X", lwd=2, xlab="numero pezzi prodotti", ylab="costo di produzione")
lines(X,Ycapp,lwd=2,col="red")
lines(X,Ycapp,lwd=2,col="green")
lines(X,Ycapp,lwd=2,col="blue")

############## OSSERVAZIONE E COMMENTO DEL GRAFICO OTTENUTO

# dal grafico si può osservare che la relazione tra numero di pezzi prodotti e 
# costo di produzione è approssimativamente lineare; 
# sembra esserci una correlazione positiva (al crescere diel numero di pezzi prodotti
# cresce il costo di produzione)
# inoltre non sono evidenti outlier (o valori anomali)

############# OSSERVARE L'OUTPUT DELL'ANALISI

#il comando summay() ci permette di ottenere maggiori risultati sull'analisi condotta

summary(result)

#il risultato ottenuto da summary contiene: 

######## 1) call : la formula in base alla quale è stato stimato il modello
#  vediamo 
#  Call:
#  lm(formula = Y ~ X)


######### 2) residuals : minimo, primo quartile, mediana, terzo quartile e massimo dei residui del modello
#  vediamo: 
#   Residuals:
#    Min      1Q      Median      3Q      Max 
#   -926.29  -461.61   -5.41    144.49  1425.51 



#########3) coefficients: per ciascun coefficiente sono riportati:
#		- la stima puntuale
#		- il corrispondente errore standard
#		- il valore della statistica test t associata al test bilatero H0:Bi=0
#		  contro l'alternativa H1:B1 diverso da 0; 
#		- p-value relativo  questo test di ipotesi
#vediamo: 

# Coefficients:
#              Estimate Std. Error  t value Pr(>|t|)    
#(Intercept)  3763.6291   750.9216   5.012  0.00154 ** 
# X              8.6923     0.3996  21.752  1.1e-07 ***


########## 4) signif. codes : la legenda da associare ad ogni p-value 
#vediamo: 
#Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1



# 5) ANOVA : il valore dell'errore standard dei residui 
# vediamo: 
# Residual standard error: 767.5 on 7 degrees of freedom

# Multiple R-squared:  0.9854,    Adjusted R-squared:  0.9833 

# F-statistic: 473.1 on 1 and 7 DF,  p-value: 1.095e-07

#####SPIEGAZIONE DI ANOVA: 
# il valore dell'errore standard dei residui con i relativi gradi di libertà
# il valore di R2 e del valore aggiustato (utile per comparare mdelli diversi)
# il valore della statistica test F dell ANOVA con i relativi gradi di libertà e il p-value associato
#relativi al seguente test di ipotesi: 
#   H0: modello senza predittore Y=B0+Err
#   H1: modello con predittore Y=B0+B1X+Err

#R2 è la proporzione di variabilità nella Y osservata, spiegata 
#   dal modello di regressione lineare basato su X

#si può notare che ESSENDO UNA REG. LIN. SEMPL. il test F coincide con il test t eseguito su B1
# ovvero F=473.1 = quadrato di t su B1=(21.752)^2 e
# il p-value dei due test coincidono 

############# IL MODELLO STIMATO

# nel nostro caso il modello stimato è: Y= 4587.3854+8.3503X

#entrambi i test t sui coefficienti b0 e b1 presentano un p-value <0
# quindi si ha alta evidenza statistica per rifiutare l'ipotesi nulla che il singolo coefficiente sia uguale a zero
# ovvero: 
#	 - l'intercetta è statisticamente diversa da zero
#	 - il predittore risulta statisticamente rilevante nel modello 

#R^2 = 0.97, indica che il modello spiega il 97% della variabilità di Y
# si può quindi ritenere un buon modello

########### TABELLA ANOVA PIU' DETTAGLIATA

# per ottenere una tabella ANOVA più dettagliata possiamo usare due funzioni
# anova () oppure
# aov ()   questa con sintassi analoga ad lm
#entrambe producono medesimo risultato

result.aov=anova(result)
result.aov

anova.table=aov(Y~X,data=dati)
summary(anova.table)

#otteniamo: 
# - valori delle somme degli scarti al quadrato (sum sq)
# - errori quadratici medi (mean sq)
# - il valore della statistica F (F value)
# - il p-value associato al relativo test (Pr>F)

################ estrazione del ERRORE STANDARD DEI RESIDUI
# Si osservi che la stima dello scarto quadratico medio, ovvero l'errore standard dei residui, 
# può essere ricavata dalla radice del Mean Square Error(relativo alla regressione)
# che si trova nella tabella di Anova e che possiamo estrarre nel seguente modo: 

#stima dell'errore standard dei residui

result.aov
sqrt(result.aov[2,3])

#ovvero sqrt(589072)