## Trasformazioni dei dati per una migliore interpretazione dei coefficienti
## Esempio libro Gelman Hill
## Roma a.a. 2024--2025
setwd("/Users/mariagraziapittau/Library/CloudStorage/Dropbox-Universita'DegliStudiDiRomaLaSapienza/grazia.pittau@uniroma1.it/Lezioni/Statistica_Aziendale/R/Regressione")

library(foreign)
library(arm)
kidiq <- read.dta(file="kidiq.dta")

attach(kidiq)

##Standardizzazione
fit.std.1 <- lm (kid_score ~ mom_hs + mom_iq + mom_hs:mom_iq)
display(fit.std.1)


##Centro rispetto alla media
c.mom.hs <- mom_hs - mean(mom_hs)
c.mom.iq <- mom_iq - mean(mom_iq)

fit.std.2 <- lm (kid_score ~ c.mom.hs + c.mom.iq + c.mom.hs:c.mom.iq)
display(fit.std.2)


##standardizzo secondo la classica standardizzazione
fit.std.3 <- lm (kid_score ~ mom_hs + scale(mom_iq) + mom_hs:scale(mom_iq))
display(fit.std.3)

##il coefficiente dell'interazione e' pari a quello del modello non standardizzato moltiplicato per la deviazione standard-0.48*sd(mom_iq)

##Standardizzazione rispetto alla media e due volte la standard deviation
z.mom.hs <- (mom_hs - mean(mom_hs))/(2*sd(mom_hs))
z.mom.iq <- (mom_iq - mean(mom_iq))/(2*sd(mom_iq))

fit.std.4 <- lm (kid_score ~ mom_hs + z.mom.iq + mom_hs:z.mom.iq)
display(fit.std.4)

fit.4 <- lm (kid_score ~ mom_hs + mom_iq + mom_hs:mom_iq, data=kidiq)
display(fit.4)



##il coefficiente dell'interazione "-14.53" e' pari a quello del modello non standardizzato moltiplicato per la deviazione standard-0.48*2sd(mom_iq) 
 
## Trasformazioni logaritmiche 
## Esempio sui guadagni in funzione dell'altezza
library(foreign)
heights <- read.dta ("heights.dta")
attach(heights)

##Per semplicita' rimuovo i casi con missing data 
##e restringo il campione alle persone nate dopo il 1925
##SI OSSERVI BENE LA RICODIFICA DI SESSO IN MALE
table(sex)
male <- 2 - sex
table(male)
ok <- !is.na (earn+height+sex) & earn>0 & yearbn>2
heights.nuovo <- as.data.frame (cbind (earn, height, sex, male)[ok,])
n <- nrow (heights.nuovo)

# modello su scala originale
lm.earn0 <- lm (earn ~ height, data=heights.nuovo)
display (lm.earn0, digits=2)
hist(earn, nclass=30)
##1260$ a fronte di una variazione di un pollice nell'altezza

##Istogramma dei residui
residui <- resid(lm.earn0)
hist(residui, breaks=50, freq=F, ylab="densita'", main="Istogramma dei residui")  
lines(density(residui))
## forte asimmetria a sinistra


# modello su scala logaritmica
log.earn <- log (heights.nuovo$earn)
logmodel.1 <- lm (log.earn ~ scale(height), data=heights.nuovo)
display (logmodel.1) 
## A fronte della variazione di un pollice nell'altezza una variazione attesa dei salari del 6%

##exp(0.06)-1
residui <- resid(logmodel.1)
hist(residui, breaks=50, freq=F, ylab="densita'", main="Istogramma dei residui")  
lines(density(residui))


y.hat<-logmodel.1$fitted.values
plot (y.hat, residui, xlab="", ylab="residui")
abline(h=0)
abline(h=-sigma.hat(logmodel.1), lty=2)
abline(h=+sigma.hat(logmodel.1), lty=2)


#Aggiungo un predittore
logmodel.2 <- lm (log.earn ~ rescale(height) + male, data=heights.nuovo)
display (logmodel.2, digits=3)
summary(logmodel.2)
residui <- resid(logmodel.2)
hist(residui, breaks=50, freq=F, ylab="densita'", main="Istogramma dei residui")  
lines(density(residui))

## Si noti exp(0.42) gli uomini guadagnano una volta e mezzo in piu delle donne e si riduce la variazione attesa dei salari dal 6% al 2%
## ma anche un differenza del 2% nei salari non e' banale: per guadagni di $50000, diventa un rilevante $1000 per pollice di altezza in piu'! 
##Si noti la significativita' del coeff. dell'altezza
##Ma la deviazione standard dei resuidi e' 0.88 su scala logaritmica che vuol dire che i guadagni previsti espressi in log saranno piu' o meno precisi di 0.88. 
## Che succede ad una persona alta 70 pollici di sesso maschile?

x.new <- data.frame (height=70, male=1)
predict(logmodel.2, x.new)
pred <- coef(logmodel.2)[1] + coef(logmodel.2)[2] * 70 ##VALORE PREVISTO DEL SALARIO IN LOG
exp(pred)
int <- pred + c(-1, 1) *sigma.hat(logmodel.2)
## su scala originale
exp(int[1])
exp(int[2])

## Ci chiediamo: una previsione di questo tipo ha senso? Eppure il coeff. e' statisticamente significativo! Si noti anche il valore molto basso di R^2

#Aggiungo l'interazione
logmodel.3 <- lm (log.earn ~ rescale(height) + male + rescale(height):male, data=heights.nuovo)
display (logmodel.3, 3)

## male non possiamo interpretarlo--non ha senso! non e' stat.significativo e neanche economicamente significativo
## interazione invece non statisticamente significativa ma ha un significato importante
## un incremento di altezza pari a un pollice implica un aumento del 18% nei guadagni se l'individuo e' donna 
## se invece l'individuo e' un uomo l'incremento di altezza diventa 18%+0.7%=2.5%
##quindi anche se non e' stat. signifificativo decidiamo di tenerlo nel modello

#Standardizzo l'altezza al fine di poter interpretare l'intercetta e il coeff di male
z.height <- (height - mean(height))/(2*sd(height))
logmodel.3a <- lm (log.earn ~ z.height + male + z.height:male,data=heights.nuovo)
display (logmodel.3a, digits=3)

##13767 dollari il guadagno atteso per un uomo con altezza pari alla media pari a exp(9.527)
## gli uomini gaudagnano il 52% in piu' delle donne  exp(0.42)-1
## circa 14% in piu' sui salari attesi per due persone di esso femminile con differenza di altezza pari a 7.69 pollici pari a due volte la sd, ovvero per 20 cm di differenza nelle altezze 
## negli uomini questa differenza aumenta diventando circa 18% per gli uomini 



#Modello doppio logaritmo
##elasticita
loglogmodel.1 <- lm (log.earn ~ log(height) + male, data=heights.nuovo)
display (loglogmodel.1)


##Ricordare altre trasformazioni: per esempio quella quadatrica