## Dati European Social Survey 2022
## code per Regressione esercitazione--a.a. 2024--2025
require(arm)
require(car)
load("/Users/mariagraziapittau/Library/CloudStorage/Dropbox-Universita'DegliStudiDiRomaLaSapienza/grazia.pittau@uniroma1.it/Lezioni/Statistica_Aziendale/Assignments_2024-2025/EuropeanSocialSurvey_2020.rda")

table(ess_data$country)

##Considero solo alcune nazioni Italia, Germania, Spagna e Rep. Ceca
##Si restringa l’analisi alle sole nazioni Italia (IT), Germania (DE), Francia (FR), Regno Unito (GB), Bulgaria (BG) e Olanda (NL) 

ess = subset(ess_data, ess_data$country == "IT" | ess_data$country == "FR" |ess_data$country == "CZ" | ess_data$country == "GB"|ess_data$country == "BG" |ess_data$country == "NL" )
table(ess$country)
 
hist(ess$age)
Rescale=function(x){
	
	(x-mean(x))/(2*sd(x))
	
}


# fit
ess$country <-  relevel(as.factor(ess$country), ref = "IT")
ess$health <- relevel(as.factor(ess$health), ref = "Pessima")


## Education
table(ess$edu)

istr<- as.numeric(ess$edu)
istr <- recode(istr, "1=1; 2=2; 3=4; 4=4; 5=5")
table(istr)
istr<-factor(istr, levels=c(1,2,4,5), labels=c("Licenza Elementare", "Licenza Media", "Licenza Superiore/Laurea Triennale", "Laurea/Master/Dottorato"))
ess$istr<-istr


ess$happy<-as.numeric(ess$fact.happy)
boxplot(ess$happy~ess$country)

z.age<-Rescale(ess$age)

fit.1 = lm(happy ~ country + female + z.age + I(z.age^2) + health + istr + employed + safe + 
            Rescale(partyid) + Rescale(income), data = ess)
summary(fit.1)
display(fit.1)


fit.2 = lm(happy ~ country + female + z.age + I(z.age^2) + health + istr + employed + safe + 
            Rescale(partyid) + Rescale(income) + Rescale(income):Rescale(partyid), data = ess)

display(fit.2)
summary(fit.2)




##Un po' di diagnostica
par(mfrow=c(2,2))
plot(fit.2$fitted.values, fit.2$residuals, xlab = "Valori Stimati", ylab = "Residui", cex = 0.6)
abline(h = 0) 
abline(h=-sigma.hat(fit.2), lty=2)
abline(h=+sigma.hat(fit.2), lty=2)

qqnorm(fit.2$residuals, cex = 0.6, main = "", xlab = "Quantili", ylab = "Residui")
qqline(fit.2$residuals)

hist(fit.2$residuals, freq = FALSE, col = "gray", main = "", xlab = "Residui",ylab = "Densita'", ylim=c(0, 1.2)) 
lines(density(fit.2$residuals))
curve(dnorm(x, 0, sigma.hat(fit.2)), -3 * sigma.hat(fit.2), 4 * sigma.hat(fit.2), add = TRUE, col="red")



x <- predict(fit.1)
y <- resid(fit.1)
binnedplot(x, y, nclass=80, xlab="Valori medi stimati", ylab="Residui medi", main="Grafico dei residui binned")

##
## Ricco =9 
z.ricco= (9 - mean(ess$income))/(2*sd(ess$income))
## Povero =2 
z.povero= (2 - mean(ess$income))/(2*sd(ess$income))
##ITALIA

dev.off()
pdf("/Users/mariagraziapittau/Library/CloudStorage/Dropbox-Universita'DegliStudiDiRomaLaSapienza/grazia.pittau@uniroma1.it/Lezioni/Statistica_Aziendale/aa.pdf")
par(mfrow=c(2,2))
plot (range(rescale(ess$partyid)), range(ess$happy), xlab="Orientamento Politico (riscalato)", ylab="Felicita'", pch=20, type="n", main="", ylim=c(4.5,6.2))
curve (coef(fit.2)[1] + coef(fit.2)[20] * rescale(x) + coef(fit.2)[21] * z.ricco + coef(fit.2)[22]  * rescale(x) * z.ricco, add=TRUE, col="blue", lty=2) ## Ricco 
curve (coef(fit.2)[1] + coef(fit.2)[20] * rescale(x) + coef(fit.2)[21] * z.povero + coef(fit.2)[22]  * rescale(x) * z.povero, add=TRUE, col="red", lty=2) ## Povero

text(.2, 6.0, "Ricchi", cex=1, col="blue")
text(.2, 5.3, "Povero", cex=1, col="red")



##BULGARIA
plot (range(rescale(ess$partyid)), range(ess$happy), xlab="Orientamento Politico (riscalato)", ylab="Felicita'", pch=20, type="n", main="", ylim=c(3.5,5.5))
curve (coef(fit.2)[1] +coef(fit.2)[2] + coef(fit.2)[20] * rescale(x) + coef(fit.2)[21] * z.ricco + coef(fit.2)[22]  * rescale(x) * z.ricco, add=TRUE, col="blue", lty=2) ## Ricco 
curve (coef(fit.2)[1]+ coef(fit.2)[2] + coef(fit.2)[20] * rescale(x) + coef(fit.2)[21] * z.povero + coef(fit.2)[22]  * rescale(x) * z.povero, add=TRUE, col="red", lty=2) ## Povero

text(.2, 5.2, "Ricchi", cex=1, col="blue")
text(.2, 4.5, "Povero", cex=1, col="red")


##Curva della felicita'
plot (range(rescale(ess$age)), range(ess$happy), xlab="Eta' (riscalata)", ylab="Felicita'", pch=20, type="n", main="Italia", ylim=c(5.4,5.7))
curve (coef(fit.2)[1] + coef(fit.2)[8] * rescale(x) + coef(fit.2)[9] * I(rescale(x)^2), add=TRUE, col="blue", lty=2)