#-----------------------------------------------------#
#                                                     #
#        LABORATORIO DI METODI STATISTICI             #
#                PER LA FINANZA                      #
#               PROF F. BARTOLUCCI                    #
#                                                     #
#                    LAB-06                          #
#-----------------------------------------------------#

##########################################################################
# Lezione 6
##########################################################################1
#                                                                           
#   Indice:                                                                 
#           
#   * Regressione lineare semplice: forma matriciale
#   * Regressione lineare semplice: funzione lm()
#   * Regressione lineare semplice: tabella ANOVA
#   * Regressione lineare semplice: analisi dei residui
#   * Regressione lineare multipla: introduzione
# 
##########################################################################


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

names(dati)

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

# Diamo alle variabili dei nomi piu semplici
names(dati) = c("cranio","altezza","gest","peso","eta.mamma","tossemia")


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




# Regressione lineare semplice: forma matriciale 
# -----------------------------------------------
# matrice del disegno associata al modello
n = length(y)
X = cbind(rep(1,n),x)

# Stima dei parametri
# ------------------------------------------------
betahat = solve(t(X) %*% X) %*% t(X) %*% y
betahat 
rownames(betahat) = c("Intercept", "Peso")
betahat

# valori predetti di y
# --------------------------
yhat = X %*% betahat
yhat

# residui
# --------
e = y - yhat

# standard error del vettore beta
# --------------------------------------------------------

# stima della varianza dell'errore sigma2
RSS = sum(e^2)
s2 = RSS / (n-2)
s2 

# stima della matrice di var e covar del vettore beta
V = s2 * solve(t(X)%*%X)
V

# stima dello standard error del vettore beta
se = sqrt(diag(V))
se

# statistica-test t per verificare la significatività dei singoli parametri beta: H: beta = 0, H1: beta != 0
t = betahat / se

# calcoliamo il p-value
pvalue = 2*(pt(abs(t), n-2, lower.tail = FALSE))
pvalue
# R: Rifiutiamo entrambe le ipotesi nulle


# creiamo una tabella riepilogativa delle stime
coeff = cbind(betahat, se, t, pvalue)
coeff
# diamo un nome alle colonne
colnames(coeff) = c("Estimates", "Se", "t", "p-value")
coeff

# intervalli di confidenza per i beta
t = qt(0.05/2, df = n-2, lower.tail = FALSE)

l1.be0 = betahat[1,] - t * se[1]
l2.be0 = betahat[1,] + t * se[1]
l1.be1 = betahat[2,] - t * se[2]
l2.be1 = betahat[2,] + t * se[2]
int.conf = matrix(c(betahat[1], betahat[2], l1.be0, l1.be1, l2.be0, l2.be1), ncol = 3)
int.conf
colnames(int.conf) = c("Estimate", "Lower", "Upper")
int.conf



# Previsione di un valore di Y
x.new = c(1, 790)
yhat.new = x.new %*% betahat
yhat.new



# Regressione lineare semplice: funzione lm()
# ------------------------------------------------
mod.peso = lm(y ~ x)
mod.peso

# qualcosa in più sulle stime e sul fit del modello
summary(mod.peso)


# rappresentazione grafica
# --------------------------------
plot(x, y)
abline(mod.peso, col = "red")


## Estrazione delle informazioni dall'oggetto lm
# ------------------------------------------------
# coefficienti
coefficients(mod.peso) # oppure:
mod.peso$coefficients

# fitted values
yhat = fitted(mod.peso) 
yhat
# oppure:
yhat = mod.peso$fitted.values
yhat

# residui
res = residuals(mod.peso) 
res 
# oppure:
res = mod.peso$residuals
res

# varianza dell'errore
s2 = summary(mod.peso)$sigma^2
s2

# r^2
summary(mod.peso)$r.squared

# r^2 aggiustato che tiene conto del numero dei parametri del modello
summary(mod.peso)$adj.r.squared

# intervalli di confidenza per i parametri
confint(mod.peso, level = 0.95)



# Previsione di un nuovo valore di y
x.new = data.frame(x = 790)
predict(mod.peso, newdata = x.new, se.fit = TRUE)




# Scomposizione della devianza
# ----------------------------------
# Explained sum of squares
ESS = sum((yhat - mean(y))^2)
ESS
dfESS = ncol(as.matrix(x))

# residual sum of squares
RSS = sum((y - yhat)^2)
RSS
dfRSS = length(y) - ncol(as.matrix(x))-1

# total sum of squares
TSS = RSS + ESS
TSS
dfTSS = length(y) -1

# costruiamo la tabella di analisi della varianza (ANOVA)
df = c(dfESS, dfRSS)
sumSQ = c(ESS, RSS)
meanSQ = c(sumSQ / df)
f = meanSQ[1]/meanSQ[2]
pvalue = 2*pf(f, df[1], df[2], lower.tail = FALSE)

anova = cbind(sumSQ, df, meanSQ)
anova

# manca la statistica test F ed il pvalue.. proviamo ad avere una tabella ANOVA "completa"
cbind(anova,c(f, NA), c(pvalue, NA))


# più velocemente
anova(mod.peso)



# intervallo  di confidenza per la media
# ---------------------------------------------------------------------
pred.med = predict(mod.peso, data.frame(x), int = "confidence")
head(pred.med)

plot(x,y)
matlines(x[order(x)], pred.med[order(x),], col = c(1,2,2), lty = c(1,2,2))


# intervallo predittivo
# ------------------------------------------------------
pred.y = predict(mod.peso, newdata = data.frame(x), int = "p")

# sovrapponiamo quest'altro intervallo al grafico (che chiaramente sarà più ampio)
matlines(x[order(x)],pred.y[order(x),], col=c(1,3,3), lty=c(1,3,3))





# residui del modello di regressione
# ----------------------------------

# valori stimati vs valori osservati
plot(mod.peso$fitted.values, y)

# valori stimati vs residui
plot(mod.peso$fitted.values, mod.peso$residuals, xlab = 'Fitted values', ylab = 'Residuals')

# Normal Probability Plot
qqnorm(mod.peso$residuals, ylab = 'Residuals')

# più velocemente...
par(mfrow = c(1,2))
plot(mod.peso, which = 1:2)


# Modello di regressione lineare multipla: introduzione
# -------------------------------------------------------

# consideriamo la variabile "cranio" come variabile dipendente (y)
# e "peso" e "gest" come variabili esplicative (X)

y = dati$cranio
X = cbind(dati$peso, dati$gest)

# matrice del disegno associata al modello
n = length(y)
X = cbind(rep(1,n),X)

# Stima dei parametri
# ------------------------------------------------
betahat = solve(t(X) %*% X) %*% t(X) %*% y
betahat 
rownames(betahat) = c("Intercept", "Peso", "Gest")
betahat

# valori predetti di y
# --------------------------
yhat = X %*% betahat
yhat

# residui
# --------
e = y - yhat


# standard error del vettore beta
# --------------------------------------------------------

# stima della varianza dell'errore sigma2
RSS = sum(e^2)
s2 = RSS / (n-ncol(X))
s2 

# stima della matrice di var e covar del vettore beta
V = s2 * solve(t(X)%*%X)
V

# stima dello standard error del vettore beta
se = sqrt(diag(V))
se

# coefficiente di determinazione R^2
R2 = 1-RSS / TSS

# statistica-test t per verificare la significatività dei singoli parametri beta: H: beta = 0, H1: beta != 0
t = betahat / se

# calcoliamo il p-value
pvalue = 2*(pt(abs(t), n-ncol(X), lower.tail = FALSE))
pvalue
# R: Rifiutiamo l'ipotesi di coefficienti nulli

# più velocemente
mod.peso.gest = lm(cranio ~ peso + gest, data = dati)
summary(mod.peso.gest)


savehistory()