##########################################################################
                                                                          
# Laboratorio di "Metodi Statistici avanzati" con R                         
# Prof.ssa Lea Petrella                                                                                                                                
# Facolta'  di Economia 
# Dipartimento Metodi e Modelli per l’'Economia il Territorio e la Finanza 
   
# Sapienza Universita'  di Roma 
# a.a. 2015-2016   


# alberto.diiorio@uniroma1.it
                                                         
##########################################################################

##########################################################################
# Lezione 4: MODELLO di REGRESSIONE LINEARE SEMPLICE                                   
##########################################################################
#                                                                           
#   Indice:                                                                 
#                                                                           
#   * richiami teorici
#   * dati simulati dal "vero modello", stima della regressione lineare semplice
#     e inferenza sui parametri
#   * analisi dei residui
#	* previsione
#   * applicazione a dati reali (circonferenza cranica)
#
##########################################################################


#----------------------------------------------------------
#   * richiami teorici
#----------------------------------------------------------

# Y[i]=a+bX[i]+epsilon[i]   

# Y variabile risposta-dipendente
# X variabile predittiva-indipendente-regressore 
# (deterministica...non aleatoria!)
# epsilon e' il termine di errore...(aleatorio) ---> Y aleatoria

# Parametri del modello:

# a: intercetta 
# b: coefficiente di regressione
# sigma^2 : varianza dell'errore

# assunzioni di base (deboli):
# E(epsilon[i])=0  ;  V(epsilon[i])=sigma^2  ;  COV(epsilon[i],epsilon[j])=0

# ---> 

# E(Y[i])=a+bX[i]  ;  V(Y[i])=sigma^2  ;  COV(Y[i],Y[j])=0

# stimatore MQO: min (sum(y[i]-(a+bx[i]))^2)

# ---> hat(a) e hat(b)  (Gauss-Markov: Stimatori BLUE)

# NOTA: finora nessuna ipotesi sulla distribuzione di epsilon!

# per fare inferenza (intervalli di confidenza-test di ipotesi)
# necessarie assunzioni sulla distribuzione della v.a. epsilon

# epsilon[i]~N(0,sigma^2) ---> Y[i]~N(a+bX[i], sigma^2)


#----------------------------------------------------------
#   * dati simulati dal "vero modello", stima della 
#     regressione lineare semplice e inferenza sui parametri  
#----------------------------------------------------------

#x=c(0,1,2,5,7,10,14,20)
#y=c(375,450,500,725,800,950,1025,1200)

n=100
set.seed(123)
x=round(runif(n,5,25),1)

true.alpha=5
true.beta=3
true.sigma2=100

y=rnorm(n,true.alpha+true.beta*x,sqrt(true.sigma2))

# Ispezione grafica preliminare attraverso lo
# Scatter plot (o grafico di dispersione)
# Ogni punto e' relativo ad un'unita'  statistica
# individuato dall'ascissa, pari al valore della variabile indipendente,
# e dall'ordinata, pari al valore della variabile dipendente.
# E' importante osservare l'orientamento della nuvola dei punti!

plot(x,y)

# dal grafico si evince come ci sia una certa "dipendenza lineare" di segno pisitivo...

# Verifichiamo attraverso il coefficiente di correlazione 
# l'eventuale presenza di dipendenza LINEARE

cor(x,y)

# Per implementare il modello di regressione lineare
# si impiega il comanda lm...

mod1=lm(y~x)

# NOTA: per digitare ~ 
# con mac: alt+5
# con Windows: alt+126
# altrimenti...copia ed incolla...

summary(mod1)

# Visto che il comando summary ci fornisce molte informazioni utili
# salviamo il suo output in un oggetto a parte, vedremo poi
# come sfruttarlo 

sintesi1=summary(mod1) 


# Vediamo la struttura dell'output 
# del comando lm...

str(mod1)

# ...e' piuttosto articolata ma possiamo estrarre elementi
# importanti per l'analisi di regressione attraverso il $

# ad esempio, per ottenere le stime dei paramteri della retta di regressione,
# i valori di y stimati (o predetti) dalla retta di regressione
# e i residui...

mod1$coefficients
mod1$fitted.values
mod1$residuals

# oppure usiamo le funzioni specifiche 

coef(mod1)
fitted(mod1) 
residuals(mod1)

# Salviamo le stime dei parametri della retta in oggetti separati

stima.alpha=mod1$coefficients[1]
stima.beta=mod1$coefficients[2]

# Possiamo aggiungere alla nuvola dei punti la retta di regressione stimata

plot(x,y,xlim=c(0,30),ylim=c(0,100),main="Scatter plot e la retta di regressione stimata")
abline(a=stima.alpha,b=stima.beta,lty=2)

# oppure

plot(x,y,xlim=c(0,30),ylim=c(0,100),main="Scatter plot e la retta di regressione stimata")
abline(mod1,lty=2)


# Scomposizione della devianza

# SST= SSR + SSE

SST=var(y)*(n-1)		                # devianza totale
SSR=var(mod1$fitted.values)*(n-1)		# devianza della regressione
SSE=var(mod1$residuals)*(n-1)			# devianza residua

# Verifichiamo che i calcoli siano giusti

SST
SSE+SSR

# ma allora l'indice di determinazione lo potevo ottenere così...

coefdet=SSR/SST
coefdet

# oppure, nel caso della regressione semplice, anche col
# quadrato del coeffciente di correlazione

cor(x,y)^2

# stima MQO (corretta) della varianza dell'errore

stima.cor.var.sigma2=SSE/(n-2)

# stima MV (non corretta) della varianza dell'errore

stima.ncor.var.sigma2=SSE/n


# Intervalli di confidenza

# Intervallo di confidenza per alpha al livello 1-q
# [stima.alpha-t(n-2;1-q/2)*SE(stima.alpha),stima.alpha-t(n-2;1-q/2)*SE(stima.alpha)]

# Intervallo di confidenza per beta al livello 1-q
# [stima.beta-t(n-2;1-q/2)*SE(stima.beta),stima.beta-t(n-2;1-q/2)*SE(stima.beta)]

# Intervallo di confidenza per sigma2 al livello 1-q
# [(n-2)*stima.var.sigma2/chisq(n-2;1-q/2),(n-2)*stima.var.sigma2/chisq(n-2;q/2)]

# Dunque, se ci viene richiesto l'intervallo al livello del 95%
# abbiamo che q=0.05 ma ci servono gli standard error delle stime di alpha e beta 
# che possiamo rubare al comando summary(mod1) precedentemente salvato

str(sintesi1)
sintesi1$coefficients
class(sintesi1$coefficients)
se.stima.alpha=sintesi1$coefficients[1,2]
se.stima.alpha
se.stima.beta=sintesi1$coefficients[2,2]
se.stima.beta


alpha.int95=c(stima.alpha-qt(p=1-0.025,df=n-2)*se.stima.alpha,stima.alpha+qt(p=1-0.025,df=n-2)*se.stima.alpha)
alpha.int95

beta.int95=c(stima.beta-qt(p=1-0.025,df=n-2)*se.stima.beta,stima.beta+qt(p=1-0.025,df=n-2)*se.stima.beta)
beta.int95

# per gli intervalli di alpha e beta esiste anche il comando diretto  

confint(mod1)

# Di default calcola gli intervalli al 95% ma possiamo cambiare
# livello con l'opzione level

# ora vediamo l'intervallo per la varianza dell'errore

sigma2.int95=c((n-2)*stima.cor.var.sigma2/qchisq(p=0.975,df=n-2),(n-2)*stima.cor.var.sigma2/qchisq(p=0.025,df=n-2))
sigma2.int95


# Verificate sempre che la stima puntuale appartenga
# all'intervallo di confidenza, altrimenti c'e' qualcosa che non va!

#----------------------------------------------------------
#   * analisi dei residui
#----------------------------------------------------------

predetti=mod1$fitted.values
residui=mod1$residuals


# Verifichiamo che la media sia nulla con un opportuno test

t.test(residui)		  

# Per valutare l'adeguatezza dell'ipotesi di omoschedasticita'
#plottiamo residui (in ordinata) verso valori predetti (in ascissa)
# e vediamo come si dispongono attorno alla retta di riferimento
# in corrispondenza dello zero

plot(predetti, residui, xlab="Valori predetti", ylab="Residui")
abline(h=0, col="red",lwd=2)

# Oscillazioni non costanti attorno allo zero evidenziano
# la violazione dell'omoschedasticita' .
# Inoltre l'evidenza di particolari pattern in questo grafico 
# segnalano la violazione all'ipotesi di linearita'  nella relazione
# funzionale tra X e Y oppure la presenza di dipendenza dei residui
# Punti isolati con residui molto elevati indicano le osservazioni
# descritte male dal modello e potenzialmente anomale

# Per l'omoschedasticita'  possiamo effettuare anche il
# test di Breusch-Pagan

library(lmtest)

bptest(mod1)

# Per l'incorrelazione si puo' costruire il
# correlogramma

acf(residui)

# e anche effettuare il test di 
# Durbin-Watson che sottopone a verifica l'ipotesi 
# di assenza di correlazione seriale (autocorrelazione)
# del primo ordine
# si trova sempre nello stesso pacchetto lmtest!

dwtest(mod1)
 
# Verifichiamo la conformita' dei residui alla distribuzione normale
# dando un'occhiata alla distribuzione dei residui 

hist(residui,prob=T,main="Distribuzione dei residui",xlab="Residui")
abline(v=0,lwd=3)

plot(residui)
abline(h=0,lty=2)

# Costruiamo il qqplot per i residui standardizzati

qqnorm(scale(residui))
abline(0,1)

# Test di normalita'  di Jarque-Bera

library(tseries)
jarque.bera.test(residui)


# Test di normalita'  di Shapiro-Wilks

shapiro.test(residui)  

#----------------------------------------------------------
#   * previsione 
#---------------------------------------------------------- 

# riprendiamo i dati di iniziali

set.seed(123)
x=round(runif(n,5,25),1)

true.alpha=5
true.beta=3
true.sigma2=100

y=rnorm(n,true.alpha+true.beta*x,sqrt(true.sigma2))

d <- data.frame(x,y)
mod1 <- lm(y~x,data=d)

# Fare previsione significa utilizzare il modello stimato
# per predire il valore di y in corrispondenza
# di nuove osservazioni della x, ovvero di dati che non hanno contribuito
# alla stima del modello 

 
# Primo passo: creiamo dei "nuovi" dati (nuove x) su cui prevedere i valori di y

x.new=seq(0,40)
nd <- data.frame(x=x.new)

# Previsione di tipo 1: siamo interessati a stimare il valore MEDIO di Y 
# per un dato nuovo valore di X
# e a calcolare il corrispondente intervallo di confidenza

p_conf <- predict(mod1,interval="confidence",newdata=nd)

# Previsione di tipo 2: siamo interessati a stimare il valore singolo di Y
# per un dato nuovo valore di X
# e a calcolare il corrispondente intervallo di confidenza

p_pred <- predict(mod1,interval="prediction",newdata=nd)

# Il valore individuale predetto
# e' uguale a quello ottenuto in precedenza ma l'errore standard e'
# diverso in quanto tiene conto anche della dispersione di y intorno
# alla media. 

p_conf[,1]
p_pred[,1]

# il comando da utilizzare è predict

 
# Plottiamo tutto insieme...otteniamo lo stesso grafico di prima
 
plot(y~x,xlim=c(0,30),ylim=c(0,100)) ## data

matlines(nd, p_conf, col=c(1,3,3), lty=c(1,3,3))	# col=3 corrisponde al green
matlines(nd, p_pred, col=c(1,2,2), lty=c(1,2,2))	# col=2 corrisponde al red

#----------------------------------------------------------
#   * applicazione a dati reali (circonferenza cranica)
#----------------------------------------------------------

dati=read.csv("Circonferenza-Cranica.csv",header=T,sep=";")
str(dati)
names(dati)
n=nrow(dati)

# headcirc: circonferenza cranica
# lenght: lunghezza / altezza del bambino
# gestage: eta'  gestazionale 
# birthwt: peso alla nascita
# momage: eta' della madre al momento del parto
# toxemia: presenza/assenza di tossemia della madre durante la gravidanza

# anzitutto codifichiamo la tossemia come fattore

dati$toxemia=factor(dati$toxemia)

# rinominare le variabili con nomi piu semplici

names(dati)=c("cranio","altezza","gest","peso","eta.mamma","tossemia")

# vediamo la struttura di correlazione nei dati
cor(dati[,-6])

# cerchiamo di facilitare la lettura della matrice di correlazione

library(corrplot)

cor.mat=cor(dati[,-6])
corrplot(cor.mat,method="circle")

# consideriamo la variabile "cranio" come variabile dipendente (y)
# e "peso" come variabile esplicativa (x)

plot(dati)

# e attacchiamo il file in modo dal lavorare in modo più agevole

attach(dati)

mod.peso=lm(cranio~peso)
summary(mod.peso)
sintesi.mod.peso=summary(mod.peso)

# Nuvola dei punti piu' retta di regressione stimata

plot(peso,cranio,main="Scatterplot con retta di regressione stimata")
abline(mod.peso)


# Scomposizione della devianza

# SST= SSR + SSE

SST=var(cranio)*(n-1)		                # devianza totale
SSR=var(mod.peso$fitted.values)*(n-1)		# devianza della regressione
SSE=var(mod.peso$residuals)*(n-1)			# devianza residua

# Verifichiamo che i calcoli siano giusti

SST
SSE+SSR
SSR/SST
cor(peso,cranio)^2

# stima MQO (corretta) della varianza dell'errore

stima.cor.var.sigma2=SSE/(n-2)

# stima MV (non corretta) della varianza dell'errore

stima.ncor.var.sigma2=SSE/n

# Intervalli di confidenza

stima.alpha=mod.peso$coefficients[1]
stima.beta=mod.peso$coefficients[2]

se.stima.alpha=sintesi.mod.peso$coefficients[1,2]
se.stima.alpha
se.stima.beta=sintesi.mod.peso$coefficients[2,2]
se.stima.beta


alpha.int95=c(stima.alpha-qt(p=1-0.025,df=n-2)*se.stima.alpha,stima.alpha+qt(p=1-0.025,df=n-2)*se.stima.alpha)
alpha.int95

beta.int95=c(stima.beta-qt(p=1-0.025,df=n-2)*se.stima.beta,stima.beta+qt(p=1-0.025,df=n-2)*se.stima.beta)
beta.int95

# per gli intervalli di alpha e beta esiste anche il comando diretto  

confint(mod.peso)

# Di default calcola gli intervalli al 95% ma possiamo cambiare
# livello con l'opzione level

# ora vediamo l'intervallo per la varianza dell'errore

sigma2.int95=c((n-2)*stima.cor.var.sigma2/qchisq(p=0.975,df=n-2),(n-2)*stima.cor.var.sigma2/qchisq(p=0.025,df=n-2))
sigma2.int95

# analisi dei residui

predetti=mod.peso$fitted.values
residui=mod.peso$residuals


# Verifichiamo che la media sia nulla con un opportuno test

t.test(residui)		   


# Omoschedasticita' 

plot(predetti, residui, xlab="Valori predetti", ylab="Residui")
abline(h=0, col="red",lwd=2)


library(lmtest)

bptest(mod.peso)

# Incorrelazione

acf(residui)

dwtest(mod.peso)
 
# Verifichiamo la conformita' dei residui alla distribuzione normale
# dando un'occhiata alla distribuzione dei residui 

hist(residui,prob=T,main="Distribuzione dei residui",xlab="Residui")
abline(v=0,lwd=3)

plot(residui)
abline(h=0,lty=2)


# Costruiamo il qqplot per i residui standardizzati

qqnorm(scale(residui))
abline(0,1)

# Test di normalita'  di Jarque-Bera

library(tseries)
jarque.bera.test(residui)


# Test di normalita'  di Shapiro-Wilks

shapiro.test(residui)