## Regressione lineare semplice e multipla

# carichiamo i dati nel file "insulate" sul consumo di riscaldamento prima e dopo un intervento di coibentazione
# i dati contengono le osservazioni della temperatura esterna media (in gradi Celsius): temp
# e del consumo settimanale di gas (in migliaia di piedi cubi): cons
# per 26 settimane prima e 30 settimane dopo l'intervento: when

insulate <- read.table('insulate.dat', col.names=c("when","temp","cons"))

str(insulate)
insulate$when<-as.factor(insulate$when)
summary(insulate)
attach(insulate)  # per accedere direttamente alla variabili

plot(temp, cons, pch=19)

# modello lineare semplice 
ins.lm <- lm(cons ~ temp)   

ins.lm$coefficients
ins.lm$residuals
# # not run
# coef(ins.lm)
# fitted(ins.lm)
# resid(ins.lm)
confint(ins.lm)

# output completo
summary(ins.lm)

plot(temp, cons)
abline(ins.lm)

# analisi dei residui
par(mfrow=c(1,2), pty="s", mar=c(3,2,3,2))
plot(ins.lm, which = 1:2)

#predict: stima della media di Y o una previsione in corrispondenza di un nuovo valore di Y:
dev.off()
ins.mean <- predict(ins.lm, interval="none")
head(ins.mean)
plot(temp, cons, pch=16, ylab="predicted y")
points(temp, ins.mean, pch=8, col=2)

# nuove unità
new.data <- data.frame(temp=seq(-5, 15))
ins.meannew <- predict(ins.lm, newdata = new.data, interval="none")
head(ins.meannew)

plot(temp, cons, pch=16, xlim=range(new.data$temp), ylab="predicted y")
points(new.data$temp, ins.meannew, pch=8, col=2)
lines(new.data$temp, ins.meannew, col=2)

## intervali di confidenza: riflette l'incertezza attorno ai valori medi di previsione
ins.conf <- predict(ins.lm, interval="confidence")
head(ins.conf)
matplot(sort(temp), ins.conf[order(temp),], type="l", ylim=range(cons), 
        lty=c("solid","dashed","dashed"), col=c(2,4,4), 
        xlab="temp", ylab="predicted y")
points(temp, cons, pch=16)

# sui nuovi dati
#intervallo di confidenza e previsione

ins.confnew <- predict(ins.lm, newdata = new.data, interval="confidence")
ins.pred <- predict(ins.lm, newdata = new.data, interval="prediction")

matplot(new.data$temp, ins.pred, type="l", lty=c("solid","dashed","dashed"),
        col=c(2,6,6), xlab="temp", ylab="predicted y")
matlines(new.data$temp, ins.confnew, type="l", lty=c("solid","dashed","dashed"), 
         col=c(2,4,4))
points(temp, cons, pch=16)

## modello con più esplicative 

contrasts(when)

ins.lm.w <- lm(cons ~ temp*when)
summary(ins.lm.w)
# oppure
# summary(lm(cons ~ temp + when + temp:when))

plot(temp, cons)
abline(sum(ins.lm.w$coefficients[c(1,3)]), sum(ins.lm.w$coefficients[c(2,4)]), col=2)  # prima 
abline(ins.lm.w$coefficients[1],ins.lm.w$coefficients[2], col=3)  # dopo
legend("topright", c("prima", "dopo"), col=2:3, lty=1, bty="n", cex=1.5)

# confronto col modello lineare semplice 
anova(ins.lm,ins.lm.w)
# si rifiuta la nullità di beta3 e beta4

## intervalli di confidenza per il consumo medio con modello con interazione
new.data<-data.frame(temp=c(-2,-2,10,10),
           when=c("prima","dopo","prima","dopo"))
ins.confnew <- predict(ins.lm.w, newdata=new.data, interval="confidence")
ins.pred <- predict(ins.lm.w, newdata=new.data, interval="prediction")
ins.confnew 
ins.pred

detach(insulate)



###################################################
############# Orange Juice Regression #############
###################################################

oj <- read.csv("oj.csv") 
head(oj)
dim(oj)
# vendite (sales) settimanali per tre brand di OJ 
# costi (prices)
# brand: Tropicana, Minute Maid e Dominick’s 
# feat è una variabile dummy =1 se pubblicizzat, 0 altrimenti

oj$brand<-as.factor(oj$brand)
levels(oj$brand)

oj$feat<-as.factor(oj$feat)

bcol <- c("green","red","blue")
par(mfrow=c(1,4))
plot(sales ~ price, data=oj, col=bcol[oj$brand])
plot(log(sales) ~ price, data=oj, col=bcol[oj$brand])
plot(log(sales) ~ log(price), data=oj, col=bcol[oj$brand])
plot(log(price) ~ brand, data=oj, col=bcol)     

# modello lineare con log(Y) e log(x)
m1 <- lm(log(sales) ~ log(price), data=oj)
summary(m1) # coef, tests, fit

#un aumento dell'1% di prices corrisponde ad 
# una diminuzione delle vendite pari all'1.58%.
coef(m1)    

# due esplicative: aggiungiamo brand
m2 <- lm(log(sales) ~ log(price) + brand, data=oj)
summary(m2)
#coef(m2)

# le vendite calano di circa il $3\%$ 
# per ogni incremento dell'$1\%$ del prezzo
# a parità di prezzo, Tropicana realizza vendite 
# maggiori di Minute Maid
# che a sua volta ha in media vendite 
# superiori al brand Dominick's.

betav<-coef(m2)  # la pendenza è la stessa
layout(1:1)
plot(log(sales) ~ log(price), data=oj, col=bcol[oj$brand], 
     cex=.1, pch=20, bty="n")
abline(a=betav[1], b=betav[2], col=bcol[1], lwd=2)
abline(a=betav[1]+betav[3], b=betav[2], col=bcol[2], lwd=2)
abline(a=betav[1]+betav[4], b=betav[2], col=bcol[3], lwd=2)
legend("bottomleft", bty="n", lwd=2, col=bcol, legend=levels(oj$brand))

## Modello con Interazioni
m3 <- lm(log(sales) ~ log(price)*brand, data=oj)
summary(m3)

# Nel modello m3, per ogni brand si ottiene una retta 
# con intercetta e coefficiente angolare differente.
# il brand di riferimento è sempre Dominick's
# la differenza non è significativa per il brand minute.maid 
# per tropicana lo è

betav3<-coef(m3)
plot(log(sales) ~ log(price), data=oj, col=bcol[oj$brand], cex=.1, pch=20, bty="n")
# dominicks (green)
abline(a=betav3[1], b=betav3[2], col=bcol[1], lwd=2)
# minute maid (red)
abline(a=betav3[1]+betav3[3], b=betav3[2]+betav3[5], 
       col=bcol[2], lwd=2)
# tropicana (blue)
abline(a=betav3[1]+betav3[4], b=betav3[2]+betav3[6], col=bcol[3], lwd=2)
legend("bottomleft", bty="n", lwd=2, col=bcol, legend=levels(oj$brand))

# modello con la variabile "feat"
## OPTIONAL
m4 <- lm(log(sales) ~ log(price)*brand + feat, 
         data=oj)
summary(m4)