rm(list=ls())
## Esercitazione
## Logistic Regression
## ANES 2024 Time Series Study
##STORIA ILLUSTRATA SULLA PENA DI MORTE: un bianco americano appartenente alla middle class la pensa come un afroamericano? E se sono ricchi? 
## E poi un ricco republicano e un ricco democratico la pensano diversamente?

## Load the data for the 2020 ANES. 
#load("/Users/mariagraziapittau/Library/CloudStorage/Dropbox-Universita'DegliStudiDiRomaLaSapienza/grazia.pittau@uniroma1.it/Lezioni/Statistica_Aziendale/R/Logistico/ANES_2020/nes20.Rdata")
load("nes20.RData") 
str(nes20)
dim(nes20)
head(nes20)
##FTF face to face

options(repos = c(CRAN = "http://cran.rstudio.com"))

ipak <- function(pkg){
  new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
  if (length(new.pkg)) 
    install.packages(new.pkg, dependencies = TRUE)
  sapply(pkg, require, character.only = TRUE)
}

packages <- c("arm", "car") 
ipak(packages)

#RICODIFICA ####

names(nes20)
# La variabile originale sulla pena di morte: death_penalty
table(nes20$death_penalty)
nes20$death_penalty<-recode(nes20$death_penalty, "-9 :-8= NA; else=nes20$death_penalty") 
table(nes20$death_penalty, exclude=NULL)
##La ricodifico considerando una variabile dicotomica 0, 1 
# 1 se l'individuo e' favorevole alla pena di morte, 0 altrimenti
nes20$death<-recode(nes20$death_penalty, "1=1; 2=0; else=NA")
table(nes20$death)

## Considero la variabile genere del rispondente GENDER_RESPONDENT_X
table(nes20$gender, exclude=NULL)
nes20$male<-recode(nes20$gender, "1=1; 2=0; else=NA")
table(nes20$male, exclude=NULL)

# A party identification variable 
## La variabile iniziale e' (V202439)
## LEFT-RIGHT SELF PLACEMENT -9   -8   -7   -6   -5 not applicable
## (0. Left;  10. Right)
## indica l'orientamento politico del rispondente (democratico o repub). 
table(nes20$party_id)
nes20$party_id<-1+nes20$party_id
table(nes20$party_id)
nes20$party_id<-recode(nes20$party_id, "-8:-1=NA; else=nes20$party_id")
table(nes20$party_id, exclude=NULL)

##INCGROUP_PREPOST_X income (22 classi) reddito famigliare
## Family income
##PRE: SUMMARY: Total (family) income |      Freq.
##-----------------------------------------+-----------
##                             -9. Refused |        584
##-5. Interview breakoff (sufficient parti |         32
##                         1. Under $9,999 |        719
##                       2. $10,000-14,999 |        282
##                       3. $15,000-19,999 |        210
##                       4. $20,000-24,999 |        325
##                       5. $25,000-29,999 |        263
##                       6. $30,000-34,999 |        327
##                       7. $35,000-39,999 |        242
##                       8. $40,000-44,999 |        321
##                       9. $45,000-49,999 |        230
##                      10. $50,000-59,999 |        546
##                      11. $60,000-64,999 |        325
##                      12. $65,000-69,999 |        182
##                      13. $70,000-74,999 |        264
##                      14. $75,000-79,999 |        233
##                      15. $80,000-89,999 |        426
##                      16. $90,000-99,999 |        304
##                    17. $100,000-109,999 |        506
##                    18. $110,000-124,999 |        343
##                    19. $125,000-149,999 |        347
##                    20. $150,000-174,999 |        404
##                    21. $175,000-249,999 |        416
##                    22. $250,000 or more

table(nes20$hh_income)

##La trasformo in numerica e la considero continua
nes20$hh_income=recode(nes20$hh_income,"-9=NA;-5=NA")
table(nes20$hh_income, exclude=NULL)

##Age 
table(nes20$age, exclude=NULL)

nes20$age <- recode(nes20$age,"-9=NA")
table(nes20$age,  exclude=NULL)

##Age group
nes20$generation<-ifelse(nes20$age<=35, 'giovani',
                  ifelse(nes20$age>35 & nes20$age<=65,  'adulti',
                  ifelse(nes20$age>65 & nes20$age<=75,  'anziani',
                          "vecchi")))
  
table(nes20$generation)
nes20$generation<-relevel(as.factor(nes20$generation), ref ="giovani")


##Race
table(nes20$race, exclude=NULL)
nes20$race<-recode(nes20$race, "-9:-8 = NA; else=nes20$race")

nes20$race<-ifelse(nes20$race==1, "white", 
             ifelse(nes20$race==2, "black", 
             ifelse(nes20$race==3, "hispanic", "other"))) 

             
table(nes20$race)
nes20$white<-ifelse(nes20$race=="white", 1, 0)

##Education
table(nes20$high_educ, exclude=NULL)
nes20$high_educ<-recode(nes20$high_educ, "-9:-2 = NA; else=nes20$high_educ")

nes20$educ<-ifelse(nes20$high_educ==1, 'Less than hs', 
            ifelse(nes20$high_educ>=2 & nes20$high_educ<=3, 'hs',
            ifelse(nes20$high_educ==4, 'Bachelor','Graduate degree')))
            
table(nes20$educ)
 
nes20$college<-ifelse(nes20$educ=="Bachelor" |  nes20$educ=="Graduate degree", 1, 0) 
  
 
## Origins of parents 
table(nes20$parents_bornUS)
nes20$parents_bornUS<-recode(nes20$parents_bornUS, "-9:-8 = NA; else=nes20$parents_bornUS") 
table(nes20$parents_bornUS, exclude=NULL)
nes20$parents_bornUS<-ifelse(nes20$parents_bornUS==1, 1, 0) 

##  Children
table(nes20$children)
nes20$children<-recode(nes20$children, "-9 = NA; else=nes20$children") 
table(nes20$children, exclude=NULL)

nes20$nochildren<-ifelse(nes20$children==0, 1, 0)
table(nes20$nochildren, exclude=NULL)


##Religion
table(nes20$relig_importance)
nes20$relig_importance<-recode(nes20$relig_importance, "-9 : -8= NA; else=6-nes20$relig_importance") 
table(nes20$relig_importance, exclude=NULL)

##save(nes20, file="~/Dropbox/mgrazia/univ/lezioni/2021-22/Statistica_Aziendale/R/Categorico/nes20_recode.Rdata")

round(prop.table(table(nes20$death)),2)
table(nes20$race,nes20$death)
round(prop.table(table(nes20$race,nes20$death),1),2)
##Etnia e death penalty ####

df <- data.frame(prop.table(table(nes20$race,nes20$death),1)) 
df <- df[rev(c(order(df$Freq[1:4]),order(df$Freq[1:4])+4)),]
df$Var1 <- factor(df$Var1, levels = df$Var1[1:4][order(df$Freq[1:4],decreasing = T)])
library(ggplot2)
library(tidyverse)
library(haven)
ggplot(df,aes(fill=Var2,y=Freq*100,x=Var1))+geom_bar(position = "stack",stat="identity")+ 
  theme_minimal()+theme(legend.position = "none")+ 
  labs(x="Etnia",y="Frequenze percentuali",title="% favorevoli per etnia")+ 
  scale_x_discrete(labels=c("Bianchi","Neri","Ispanici","Altri","Asiatici"))+ 
  scale_fill_manual(values=c("darkgrey","#3333CC"))
#favorevoli in blu
detach(package:ggplot2)

round(prop.table(table(nes20$hh_income,nes20$death),1),2)
round(prop.table(table(nes20$hh_income>11,nes20$death),1),2)

##OUTCOME DEATH PENALTY ####
round(prop.table(table(nes20$death, nes20$relig_importance),2)*100, 2)

### religiosità e death penalty ####
df <- data.frame(prop.table(table(nes20$relig_importance,nes20$death),1)) 
#df <- df[rev(c(order(df$Freq[1:5]),order(df$Freq[1:5])+5)),] 
#df$Var1 <- factor(df$Var1, levels = df$Var1[1:5][order(df$Freq[1:5],decreasing = T)]) 
library(ggplot2)
ggplot(df,aes(fill=Var2,y=Freq*100,x=Var1))+geom_bar(position = "stack",stat="identity")+ 
  theme_minimal()+theme(legend.position = "none")+ 
  labs(x="Religiosità",y="Frequenze percentuali",title="% favorevoli per religiosità")+ 
  scale_x_discrete(labels=c("Per nulla","Poco","Moderatamente","Abbastanza","Molto"))+ 
  scale_fill_manual(values=c("darkgrey","#3333CC"))
detach(package:ggplot2)
#favorevoli in blu

## 0 - modello nullo ####

fit.logit.0 <- glm(death ~ 1,family=binomial(link="logit"),data=nes20)
display(fit.logit.0)

 ## 1 - Stimo il modello in funzione dell'etnia ####
fit.logit.1 <-glm(death ~ white, family=binomial(link="logit"), data=nes20)
display(fit.logit.1)

## 2 - Stimo il modello in funzione del party.id ####
fit.logit.2 <- glm(death~rescale(hh_income), family=binomial(link="logit"),data=nes20)
display(fit.logit.2)

## 3 - Modello per domanda 1 senza interazione con controlli####
fit.logit.3 <-glm(death ~ white + male+ rescale(relig_importance) + rescale(party_id) + rescale(hh_income) + factor(generation) + factor(educ), family=binomial(link="logit"), data=nes20)
display(fit.logit.3)

## 4 - Modello per domanda 1 con interazione con controlli ####

fit.logit.4 <-glm(death ~ white + male+ rescale(relig_importance) + factor(generation) + factor(educ) + rescale(party_id) + rescale(hh_income) + rescale(hh_income) : white, family=binomial(link="logit"), data=nes20)
display(fit.logit.4)

## 4 - Modello per domanda 2 con interazione ####
fit.logit.5 <-glm(death ~ white + male + rescale(relig_importance) +rescale(party_id)  + rescale(hh_income) + factor(generation) + factor(educ) + rescale(party_id) : rescale(hh_income), family=binomial(link="logit"), data=nes20)
display(fit.logit.5, 3)
summary(fit.logit.5)



##ANALISI DEI RESIDUI ####
#Calcolo dei residui
##Valori stimati
pred <- fit.logit.4$fitted.values

##residui
res <- residuals(fit.logit.4)

#Grafico dei residui in funzione delle probabilita' stimate
binnedplot(pred, res, nclas=150, xlab="Prob. stimate di essere favorevoli alla pena di morte", ylab="residui medi", main="grafico dei residui binned")

## Error rate ####
## modello nullo
## Costruisco la matrice senza dati mancanti (cosi al denominatore ho i casi reali e non artificialmente grandi per NA)
data.complete<-na.omit(nes20[, c("death", "white", "party_id", "hh_income", "generation", "educ", "male","relig_importance")])                  

fit.null<-glm(death~1, data=data.complete, family=binomial(link="logit"))
display(fit.null)
invlogit(0.42)

## 0.60>0.5 assegna tutti 1
y<-data.complete[, "death"]
length(y)
table(y)

## l'errore nel modello nullo e' pari a
tax.error.null<-table(y)[1]/dim(data.complete)[1]
tax.error.null
                   
## The confusion matrix ####
table(y.oss=y$death, y.tilde=round(pred))
round(prop.table(table(y.oss=y$death, y.tilde=round(pred))),2)
tax.error<-mean(pred>0.5 & y==0)+mean(pred<0.5 & y==1)
tax.error

## Il modello non è male procedo con l'interpretazione ed eventuali interazioni
## ritorno sul dataset iniziale la cui interpretazione e' facilitata dalle labels

### Curve in funzione della identita' politica
jitter.binary <- function(a, jitt=.05){
  a + (1-2*a)*runif(length(a),0,jitt)
}

##Senza standardizzare per avere la scala originale (solo per le curve ma MAI per interpretazione) ####
  
fit.logit.5a <-glm(death ~ white + male + rescale(relig_importance) +party_id  + hh_income + factor(generation) + factor(educ) + party_id : hh_income, family=binomial(link="logit"), data=nes20)
display(fit.logit.5a)


##mi_income<-mean(nes20$hh_income, na.rm=TRUE) 
##sd_income<-sd(nes20$hh_income, na.rm=TRUE)  
par(mfrow=c(2,2),mar=c(5,3,3,3))
?par
##reddito medio basso (media==12, 7 per esempio)
inc_1=7 

plot(range(nes20$party_id, na.rm=TRUE)*1.02, c(0,1),  xlab="Identificazione politica (democratico vs repubblicano)", ylab="Pr (essere favorevole alla pena di morte)", type="n", xaxs="i", yaxs="i", mgp=c(2,.5,0), main="Reddito basso e titolo di studio basso")
points (nes20$party_id, jitter.binary(nes20$death), pch=20, cex=.1)


##Curva reddito basso uomo bianco non altamente istruito (hs) adulto ####
curve(invlogit(coef(fit.logit.5a)[1]+coef(fit.logit.5a)[2]*1+ #bianco
                 coef(fit.logit.5a)[3]*1+ #maschio
                 coef(fit.logit.5a)[5]*x + #party_id è la x della curva
                 coef(fit.logit.5a)[6] * inc_1 + #reddito basso(7)
                 coef(fit.logit.5a)[7]*1 + #adulto
               coef(fit.logit.5a)[11]+ #non altamente istruito
                 coef(fit.logit.5a)[13]*x*inc_1), lty=1, add=T, col="blue",lwd=2)
##Curva reddito basso di colore non altamente istruito adulto ####
curve(invlogit(coef(fit.logit.5a)[1]+coef(fit.logit.5a)[2]*0+ #non bianco
                 coef(fit.logit.5a)[3]*1+ #maschio
                 coef(fit.logit.5a)[5]*x + #party_id è la x della curva
                 coef(fit.logit.5a)[6] * inc_1 + #reddito basso(7)
                 coef(fit.logit.5a)[7]*1 + #adulto
                 coef(fit.logit.5a)[11]+ #non altamente istruito
                 coef(fit.logit.5a)[13]*x*inc_1), lty=1, add=T, col="red",lwd=2)
text(4, 0.8, "Uomo adulto reddito medio-basso caucasico", cex=0.7, col="blue")
text(4, 0.7, "Uomo adulto reddito medio-basso di colore", cex=0.7, col="red")


##reddito molto elevato (22 per esempio) ####
inc_2=22
plot(range(nes20$party_id, na.rm=TRUE)*1.02, c(0,1),  xlab="Identificazione politica (democratico vs repubblicano)", ylab="Pr (essere favorevole alla pena di morte)", type="n", xaxs="i", yaxs="i", mgp=c(2,.5,0), main="Reddito elevato e titolo di studio basso")
points (nes20$party_id, jitter.binary(nes20$death), pch=20, cex=.1)

##UOMO RICCO BIANCO ADULTO NON ISTRUITO ####
curve(invlogit(coef(fit.logit.5a)[1]+coef(fit.logit.5a)[2]*1+ #bianco
                 coef(fit.logit.5a)[3]*1+ #maschio
                 coef(fit.logit.5a)[5]*x + #party_id è la x della curva
                 coef(fit.logit.5a)[6] * inc_2 + #reddito alto(22)
                 coef(fit.logit.5a)[7]*1 + #adulto
                 coef(fit.logit.5a)[11]+ #non altamente istruito
                 coef(fit.logit.5a)[13]*x*inc_2), lty=1, add=T, col="blue",lwd=2)

##UOMO RICCO NERO ####
curve(invlogit(coef(fit.logit.5a)[1]+coef(fit.logit.5a)[2]*0+ #non-bianco
                 coef(fit.logit.5a)[3]*1+ #maschio
                 coef(fit.logit.5a)[5]*x + #party_id è la x della curva
                 coef(fit.logit.5a)[6] * inc_2 + #reddito alto(22)
                 coef(fit.logit.5a)[7]*1 + #adulto
                 coef(fit.logit.5a)[11]+ #non altamente istruito
                 coef(fit.logit.5a)[13]*x*inc_2), lty=1, add=T, col="red",lwd=2)

text(4, 0.7, "Uomo adulto reddito alto caucasico", cex=0.7, col="blue")
text(4, 0.8, "Uomo adulto reddito alto afroamericano", cex=0.7, col="red")


##Laureati readdito basso ####
inc_1=7
plot(range(nes20$party_id, na.rm=TRUE)*1.02, c(0,1),  xlab="Identificazione politica (democratico vs repubblicano)", ylab="Pr (essere favorevole alla pena di morte)", type="n", xaxs="i", yaxs="i", mgp=c(2,.5,0), main="Reddito basso e laureati")
points (nes20$party_id, jitter.binary(nes20$death), pch=20, cex=.1)

##College reddito basso uomo bianco istruito adulto anni (adulto) ####
curve(invlogit(coef(fit.logit.5a)[1]+coef(fit.logit.5a)[2]*1+ #bianco
                 coef(fit.logit.5a)[3]*1+ #maschio
                 coef(fit.logit.5a)[5]*x + #party_id è la x della curva
                 coef(fit.logit.5a)[6] * inc_1 + #reddito basso(7)
                 coef(fit.logit.5a)[7]*1 + #adulto
                 coef(fit.logit.5a)[10]+ #altamente istruito (college)
                 coef(fit.logit.5a)[13]*x*inc_1), lty=1, add=T, col="blue",lwd=2)

##College reddito basso di colore istruito adulto anni (adulto) ####
curve(invlogit(coef(fit.logit.5a)[1]+coef(fit.logit.5a)[2]*0+ #non-bianco
                 coef(fit.logit.5a)[3]*1+ #maschio
                 coef(fit.logit.5a)[5]*x + #party_id è la x della curva
                 coef(fit.logit.5a)[6] * inc_1 + #reddito basso(7)
                 coef(fit.logit.5a)[7]*1 + #adulto
                 coef(fit.logit.5a)[10]+ #altamente istruito (college)
                 coef(fit.logit.5a)[13]*x*inc_1), lty=1, add=T, col="red",lwd=2)

text(4, 0.8, "Uomo adulto reddito medio-basso caucasico", cex=0.7, col="blue")
text(4, 0.7, "Uomo adulto reddito medio-basso di colore", cex=0.7, col="red")

##Laureati reddito elevato ####
inc_2=22 
plot(range(nes20$party_id, na.rm=TRUE)*1.02, c(0,1),  xlab="Identificazione politica (democratico vs repubblicano)", ylab="Pr (essere favorevole alla pena di morte)", type="n", xaxs="i", yaxs="i", mgp=c(2,.5,0), main="Reddito elevato e laureati")
points (nes20$party_id, jitter.binary(nes20$death), pch=20, cex=.1)
##UOMO RICCO BIANCO ADULTO NON ISTRUITO
curve(invlogit(coef(fit.logit.5a)[1]+coef(fit.logit.5a)[2]*1+ #bianco
                 coef(fit.logit.5a)[3]*1+ #maschio
                 coef(fit.logit.5a)[5]*x + #party_id è la x della curva
                 coef(fit.logit.5a)[6] * inc_2 + #reddito alto(22)
                 coef(fit.logit.5a)[7]*1 + #adulto
                 coef(fit.logit.5a)[10]+ #altamente istruito (college)
                 coef(fit.logit.5a)[13]*x*inc_2), lty=1, add=T, col="blue",lwd=2)
##UOMO RICCO NERO ####
curve(invlogit(coef(fit.logit.5a)[1]+coef(fit.logit.5a)[2]*0+ #bianco
                 coef(fit.logit.5a)[3]*1+ #maschio
                 coef(fit.logit.5a)[5]*x + #party_id è la x della curva
                 coef(fit.logit.5a)[6] * inc_2 + #reddito alto(22)
                 coef(fit.logit.5a)[7]*1 + #adulto
                 coef(fit.logit.5a)[10]+ #altamente istruito (college)
                 coef(fit.logit.5a)[13]*x*inc_2), lty=1, add=T, col="red",lwd=2)

text(4, 0.7, "Uomo adulto reddito alto caucasico", cex=0.7, col="blue")
text(4, 0.8, "Uomo adulto reddito alto afroamericano", cex=0.7, col="red")





##https://www.pewresearch.org/politics/2021/06/02/most-americans-favor-the-death-penalty-despite-concerns-about-its-administration/