# Esempio dalle slides regressione (spesa alimentare di un camiopne di 20 famiglie)

spesa <- c(5.2, 5.1, 5.6, 4.6, 11.3,8.1, 7.8, 5.8,5.1, 18,4.9, 11.8, 5.2,4.8,7.9,6.4,20,13.7,5.1,2.9)

reddito <- c(28,26,32,24,54,59,44,30,40,82,42,58,28,20,42,47,112,85,31,26)

nc <- c(3,3,2,1,4,2,3,2,1,6,3,4,1,5,3,1,6,5,2,2)

d <- data.frame(spesa, reddito, nc)# creo la matrice dati            


install.packages("questionr") # pacchetto per suddividere una variabile quantitativa in classi
library(questionr)
d$reddito_rec <- cut(d$reddito, include.lowest=FALSE,  right=TRUE, breaks=3) ## Cutting d$reddito into d$reddito_rec
class(d$reddito_rec) # dopo la suddivisione in classi la var. quantitativa diventa una variabile qualitativa con tante modalita' quante sono le classi scelte (in R e' un factor)
d$reddito_rec

levels(d$reddito_rec)

levels(d$reddito_rec) <- c("basso",  "medio", "alto")

nrow(d)

#########################################
# analisi descrittiva univariata (riportato solo per spesa ma va fatto per tutte le variabili

attach(d) # carico in memoria il dataframe d (non uso piu il simbolo $ per richiamare le variabili)

summary(spesa) # indici descrittivi
sd(spesa) # deviazione standard
hist(spesa) # istogramma
boxplot(spesa) # boxplot

hist(reddito)

#########################################
# analisi bivariata

#grafici
plot(reddito, spesa) # Diagramma di dispersione
plot(nc, spesa) # Diagramma di dispersione
plot(nc, reddito) # Diagramma di dispersione

plot(d[,c(1,2,3)])

cor(reddito, spesa) # coefficiente di regressione
cor(nc, spesa)
cor(nc, reddito)

cor(d[,1:3])

#install.packages("corrplot")
library(corrplot)

corrplot(cor(d[,1:3]))

#########################################
# Modello di regressione semplice

# model.redd: spesa = B0 + B1*reddito + e

model.redd <- lm(spesa ~ reddito, data=d)
summary(model.redd)

#########################################
# Analisi multivariata


#########################################
# Modello di regressione multipla
# model.redd.nc: spesa = B0 + B1*reddito + B2*nc + e

# STEP1: matrice di correlazione (se non c'e' multicollinearita' si puo' stimare il model.redd.nc)
cor(d[ ,1:3]) # moderata correlazione tra le due variabili esplicative

# Stima del modello (ipotizzando che tutte le assunzioni del modello siano rispettate, vedi slides)
model.redd.nc <- lm(spesa ~ reddito + nc, data = d)
summary(model.redd.nc)


# STEP 2: commentare p-value del test congiunto sui coefficienti (p-value: 8.558e-11 --> OK)
# STEP 3: Commentare R^2 (Adjusted R-squared:  0.9269 Ok anche se non migliora tantissimo rispetto al valore di R^2 del model.redd)
# STEP 4: Commentare p-value dei test sui singoli coefficienti B (eccetto intercetta) --> sembrerebbe ragionevole ritenere che tutti i coefficienti siano significativamente diversi da 0 nella popolazione (forma breve: i coeff. sono statisticamente significativi)
# STEP 5: Commentare i valori dei coefficienti (vedi slides del corso)

# Eventuale STEP 6: usare il modello per stimare la spesa alimentare familiare annua sapendo reddito annuo e numero componenti famiglia


######################################################################################
# Ulteriori diagnostiche (analisi residui)

plot(model.redd.nc, which=1:2)

# Esempio utilizzo del modello

nc_i <- 4
reddito_i <- 28

spesa_hat_i <- -1.11829+0.14821*reddito_i+0.79311*nc_i

coef(model.redd.nc)[1]

spesa_hat_i <- coef(model.redd.nc)[1] +coef(model.redd.nc)[2]*reddito_i+ coef(model.redd.nc)[3]*nc_i


#########################################
# Usiamo la variabile qualitativa al posto della var quantitativa
reddito_rec
model.reddrec <- lm(spesa ~ reddito_rec, data = d)
summary(model.reddrec)


detach(d)