####################################################################
##### Laboratorio 5- Analisi con  modelli per dati gerarchici- Aprile 2025######


## Dati 
df<- read.csv("sales1.csv")
str(df)

##rendiamo fattore
df$store <- factor(df$store)

#
# creazione di variabili

#logaritmo del valore del cliente/spesa
y <- df$value 

#cassa automatica
x <- df$auto

#negozio (variabile di gruppo)
store<- df$store

######

## Modello 'Complete pooling' - un'unica intercetta per tutti i negozi, ignorando la variabile gerarchica

lm_pooled <- lm(y ~ x) 
summary(lm_pooled)
AIC(lm_pooled)
#
## Modello 'No pooling' -1 intercetta per ciascun negozio (togliamo l'intercetta) 
lm_unpooled <- lm(y ~ x + store -1) 
## lm(y ~ x + store) 

summary(lm_unpooled)
coef(lm_unpooled)
AIC(lm_unpooled)
#
## Confronto tra complete pooling e no-pooling- solo per il grafico
plot(x, y, xlab = "cassa automatica", ylab = "valore")


plot(x[store==70], y[store==70], xlab = "cassa automatica", ylab = "valore")


## aggiungo una minima variabilità alla x per facilitare la visualizzazione
x_jitter <- x + runif(nrow(df),-.05,.05) 
plot(x_jitter[store==70], y[store==70], xlab = "cassa automatica", ylab = "valore")

# negozi da mostrare
alcuni.negozi <- c (70, 36, 1, 35, 21, 14, 71, 61)  


par(mfrow = c(2,4))
for (j in alcuni.negozi){
  plot(x_jitter[store==j], y[store==j], xlim = c(-.05,1.05),
       xlab = "cassa automatica", ylab = "valore",
       pch = 20,
       main = j)
  abline(coef(lm_pooled), lty = 2)
  abline(coef(lm_unpooled)[j+1], coef(lm_unpooled)[1], col = "red")
}
##coef(lm_unpooled)[62]


#la retta tratteggiata (pooled) è la stessa su tutti i negozi
#cambiamo la scala della y in modo che questo si veda (specifico ylim)
par(mfrow = c(2,4))
for (j in alcuni.negozi){
  plot(x_jitter[store==j], y[store==j], xlim = c(-.05,1.05),
       xlab = "cassa automatica", ylab = "valore", ylim=c(min(y), max(y)),
       pch = 20, 
       main = j)
  abline(coef(lm_pooled), lty = 2)
  abline(coef(lm_unpooled)[j+1], coef(lm_unpooled)[1], col = "red")
}
###
###
#grafico di alcune rette unpooled (sono parallele tra loro, ma _non_ con la pooled)
par(mfrow=c(1,1))
plot(x_jitter, y, xlab = "cassa automatica", ylab = "valore")
abline(coef(lm_pooled), lty = 2)
for(j in alcuni.negozi) {
	abline(coef(lm_unpooled)[j+1], coef(lm_unpooled)[1], col = "red")
    text(0.4,0.1+coef(lm_unpooled)[j+1] + coef(lm_unpooled)[1]*0.4, j, cex=0.5)
}
legend(0.6,-1, legend=c("unpooled", "pooled"), lty=c(1,2), col=2:1)


#####Modelli gerarchici######
### Partial pooling senza e con predittore (modello gerarchico)
## Libreria 
library(lme4)
#
## Modello a intercetta variabile senza predittori

M0 <- lmer(y ~ 1 + (1 | store)) 

###effetto fisso + effetto casuale

#M0b <- lmer(y ~ (1 | store)) ####alternativa di scrittura possibile 

summary(M0)

# coefficienti stimati
coef(M0)

# effetti fissi e casuali
fixef(M0)
ranef(M0)

## Modello a intercetta variabile e con predittore x (partial pooling con predittore)
M1 <- lmer(y ~ x + (1 | store))
summary(M1)

#confrontiamo i due modelli rispetto al valore di AIC
AIC(M1)
AIC(M0)

#guardiamo i coefficienti
# coefficienti stimati
coef(M1)

# effetti fissi e casuali
fixef(M1)
ranef(M1)

## Grafici
a_hat_M1 <- coef(M1)$store[,1]                
# primo elemento, intercetta

b_hat_M1 <- coef(M1)$store[,2]                
# secondo elemento, pendenza

par(mfrow = c(2,4))
for(j in alcuni.negozi){
  plot(x_jitter[store==j], y[store==j], xlim = c(-.05,1.05), ylim=c(min(y), max(y)),
       xlab = "cassa automatica", ylab = "valore", main = j, pch = 20 )
  abline(coef(lm_pooled), lty = 2, col = "gray10")
  abline(coef(lm_unpooled)[j+1], coef(lm_unpooled)[1], col = "red")
  abline(a_hat_M1[j], b_hat_M1[j],  col = "green")
}  
#la retta verde e' in mezzo tra la rossa e la nera: compromesso tra due estremi, complete pooling e no pooling


#grafico di alcune rette partial pooled (sono parallele tra loro, ma _non_ con la pooled)
par(mfrow=c(1,1))
plot(x_jitter, y, xlab = "automatic", ylab = "sales level")
abline(coef(lm_pooled), lty = 2)
for(j in alcuni.negozi) {
  abline(a_hat_M1[j], b_hat_M1[j],  col = "green")
    text(0.4,0.1+a_hat_M1[j] + b_hat_M1[j]*0.4, j, cex=0.5)
}
legend(0.6,-1, legend=c("pooled", "partial pooling"), lty=c(2,1), col=c(1,"green"))


#Inserimento di effetto casuale anche sulla pendenza
#modello con intercetta e coefficiente angolare come effetti casuali 

M2 <- lmer(y ~ x  + (x | store))
summary(M2)
AIC(M2)
#
coef(M2)
fixef(M2)
ranef(M2)


## Grafici 
a_hat_M2 <- fixef(M2)[1] + ranef(M2)$store[,1] 
b_hat_M2 <- fixef(M2)[2]+ ranef(M2)$store[,2]

par(mfrow = c(2,4))
for(j in alcuni.negozi){
  plot(x_jitter[store==j], y[store==j], xlim = c(-.05,1.05), ylim=c(min(y), max(y)),
       xlab = "cassa automatica", ylab = "valore", main = j, pch = 20 )
  abline(coef(lm_pooled), lty = 2, col = "gray10")
  abline(coef(lm_unpooled)[j+1], coef(lm_unpooled)[1], col = "red")
  abline(a_hat_M2[j], b_hat_M2[j],  col = 6)
}  


#grafico di alcune rette unpooled (non sono parallele tra loro, e _non_ con la pooled)
par(mfrow=c(1,1))
plot(x_jitter, y, xlab = "cassa automatica", ylab = "valore")
abline(coef(lm_pooled), lty = 2)
for(j in alcuni.negozi) {
  abline(a_hat_M2[j], b_hat_M2[j],  col = 6)
#	abline(coef(lm_unpooled)[j+1], coef(lm_unpooled)[1], col = "red")
    text(0.4,0.1+a_hat_M2[j] + b_hat_M2[j]*0.4, j, cex=0.5)
}
legend(0.6,-1, legend=c("pooled", "random a and b"), lty=c(2,1), col=c(1,6))

##Nota: 
#AIC e' superiore a quello con solo intercetta casuale 
#sembra ragionevole stimare un unico coefficiente angolare per tutti i negozi e non uno per negozio######################

#########################################################################################
### Inserimento di predittore a livello di gruppo
#########################################################################################
reddito <- df$income

####Modello a intercetta casuale con x come predittore e reddito come predittore a livello di gruppo

M3 <- lmer(y ~ x + reddito + (1 | store))
summary(M3)
AIC(M3)
##
coef(M3)
fixef(M3)
ranef(M3)

## Grafici 

a_hat_M1 <- fixef(M1)[1] + ranef(M1)$store                
b_hat_M1 <- fixef(M1)[2]


#creo reddito di lunghezza 85, uno per ogni store (ricorda che è reddito medio della zona in cui si trova lo store)

red.store <- unique(reddito)
str(red.store)


a_hat_M3 <- fixef(M3)[1] + fixef(M3)[3]*red.store +ranef(M3)$store
#a_hat_M3
#str(fixef(M3)[3]*reddito) i primi 4 sono uguali

b_hat_M3 <- fixef(M3)[2] 
#b_hat_M3


par(mfrow = c(2,4), mar = c(4,4,3,1), oma = c(1,1,2,1))
for(j in alcuni.negozi){
  plot(x_jitter[store==j], y[store==j], xlim = c(-.05,1.05), ylim=c(min(y), max(y)),
       xlab = "cassa automatica", ylab = "valore", cex.lab = 1.2, cex.axis = 1.1,
       pch = 20, mgp = c(2,.7,0), xaxt = "n", yaxt = "n", cex.main = 1.1, main = j)
  axis(1, c(0,1), mgp = c(2,.7,0), cex.axis = 1.1)
  axis(2, seq(-1,3,2), mgp = c(2,.7,0), cex.axis = 1.1)
  abline(a_hat_M1[j,1], b_hat_M1, lwd=.5, col="gray10")
  abline(a_hat_M3[j,1], b_hat_M3, lwd=1, col="blue")
}



####Modello con intercetta e coefficiente angolare come effetti casuali + reddito

M4 <- lmer(y ~  x + reddito + (x | store))
summary(M4)
AIC(M4)
#
coef(M4)
fixef(M4)
ranef(M4)


## Grafici di confronto, intercetta e pendenza casuale, senza e con la variabile reddito 


a_hat_M4 <- fixef(M4)[1] + fixef(M4)[3]*red.store +ranef(M4)$store[,1]
b_hat_M4 <- fixef(M4)[2] +ranef(M4)$store[,2]

par(mfrow = c(2,4), mar = c(4,4,3,1), oma = c(1,1,2,1))
for(j in alcuni.negozi){
  plot(x_jitter[store==j], y[store==j], xlim = c(-.05,1.05), ylim=c(min(y), max(y)),
       xlab = "cassa automatica", ylab = "valore", cex.lab = 1.2, cex.axis = 1.1,
       pch = 20, mgp = c(2,.7,0), xaxt = "n", yaxt = "n", cex.main = 1.1, main = j)
  axis(1, c(0,1), mgp = c(2,.7,0), cex.axis = 1.1)
  axis(2, seq(-1,3,2), mgp = c(2,.7,0), cex.axis = 1.1)
abline(a_hat_M4[j], b_hat_M4[j], lwd=1, col="blue")
  abline(a_hat_M2[j], b_hat_M2[j],  col = 6)
}