# Dati tratti da:

# https://dasl.datadescription.com/datafile/carbon-footprint-2015/?_sf_s=car&_sfm_cases=4+59943

################################################################################################

originale<-read.table("carbon-footprint-2015.txt",header=T,sep="\t")
head(originale)

colnames(originale)<-c("car","X","Y","cyl","type")

# Disegnare un grafico con la funzione plot in cui le
# variabili siano le colonne X e Y del dataset originale

plot(originale$X,originale$Y)

# Creare un nuovo dataset chiamato Gas che contenga solo
# i dati di originale riferiti a type Gas (usare la funzione which)

Gas<-originale[which(originale$type=="Gas"),]

# Creare un nuovo dataset chiamato Gas che contenga solo
# i dati di originale riferiti a type Hybrid (usare la funzione which)

Hybrid<-originale[which(originale$type=="Hybrid"),]

# Aggiungere al grafico i punti di Gas e Hybrid con la funzione points.
# Hybrid in rosso, Gas in blu

points(Hybrid$X,Hybrid$Y,col="red")
points(Gas$X,Gas$Y,col="blue")

####################################################################
# Generare il modello bivariato Y contro X per il dataset originale:

Model<-lm(Y~X,data=originale)

#-- Eseguire quanto segue e osservare gli output:

# Risultati

Summary<-summary.lm(Model)

Coef<-coefficients(Model) # coefficenti beta-0 (Intercept) e beta-1 (slope) 
CIpar<-confint(Model, level=0.95) #Intervallo di confidenza per i coefficienti
Ymod<-fitted(Model)# Valori di Y modellati
Res<-residuals(Model) # Residui (Y-Ymod)


# Bande di intervallo di confidenza

ConfBand<-predict.lm(Model,interval="confidence",level=0.95)

ConfBand<-data.frame(ConfBand)


# Bande di intervallo di predizione

Xpred<-data.frame(X=seq(min(originale$X),max(originale$X),1))

PredBand<-predict.lm(Model,Xpred,interval="prediction",level=0.95)

PredBand<-data.frame(PredBand)


# Grafico

par(oma=c(0,0,0,0),mar=c(4,4,1,1))
plot(originale$X,originale$Y,type="n",xlab="X",ylab="Y",cex=0.7)
points(originale$X,originale$Y,pch=16,cex=1,col="gray")
abline(a=Coef[1],b=Coef[2],col="red",lwd=2)
points(originale$X,ConfBand$lwr,type="l",col="blue",lwd=2,lty=3)
points(originale$X,ConfBand$upr,type="l",col="blue",lwd=2,lty=3)
points(Xpred$X,PredBand$lwr,type="l",col="orange",lwd=2,lty=5)
points(Xpred$X,PredBand$upr,type="l",col="orange",lwd=2,lty=5)


# Predizione per un nuovo valore di X:

newdata<-data.frame(X=c(27,34,40))

PredNew<-predict.lm(Model,newdata,interval="prediction",level=0.95)
PredNew<-data.frame(PredNew)

#- Da aggiungere al grafico precedente, ci sono sia i valori predetti
#- che l'intervallo di predizione ("lwr","upr")

points(newdata$X,PredNew$fit,type="p",col="black",pch=3,cex=2,lwd=3)
points(newdata$X,PredNew$lwr,type="p",col="black",pch="-",cex=2,lwd=3)
points(newdata$X,PredNew$upr,type="p",col="black",pch="-",cex=2,lwd=3)

# Salvare il grafico in formato pdf

#################################################################
# Generare il modello bivariato Y contro X per il dataset Gas:

Model<-lm(Y~X,data=Gas)

# Risultati

Summary<-summary.lm(Model)

Coef<-coefficients(Model) # coefficenti beta-0 (Intercept) e beta-1 (slope) 
CIpar<-confint(Model, level=0.95) #Intervallo di confidenza per i coefficienti
Ymod<-fitted(Model)# Valori di Y modellati
Res<-residuals(Model) # Residui (Y-Ymod)

# Bande di intervallo di confidenza

ConfBand<-predict.lm(Model,interval="confidence",level=0.95)

ConfBand<-data.frame(ConfBand)


# Bande di intervallo di predizione

Xpred<-data.frame(X=seq(min(Gas$X),max(Gas$X),1))

PredBand<-predict.lm(Model,Xpred,interval="prediction",level=0.95)

PredBand<-data.frame(PredBand)


# Grafico

par(oma=c(0,0,0,0),mar=c(4,4,1,1))
plot(Gas$X,Gas$Y,type="n",xlab="X",ylab="Y",cex=0.7)
points(Gas$X,Gas$Y,pch=16,cex=1,col="gray")
abline(a=Coef[1],b=Coef[2],col="red",lwd=2)
points(Gas$X,ConfBand$lwr,type="l",col="blue",lwd=2,lty=3)
points(Gas$X,ConfBand$upr,type="l",col="blue",lwd=2,lty=3)
points(Xpred$X,PredBand$lwr,type="l",col="orange",lwd=2,lty=5)
points(Xpred$X,PredBand$upr,type="l",col="orange",lwd=2,lty=5)

# Predizione per un nuovo valore di X:

newdata<-data.frame(X=c(28,34,39))

PredNew<-predict.lm(Model,newdata,interval="prediction",level=0.95)
PredNew<-data.frame(PredNew)

#- Da aggiungere al grafico precedente, ci sono sia i valori predetti
#- che l'intervallo di predizione ("lwr","upr")

points(newdata$X,PredNew$fit,type="p",col="black",pch=3,cex=2,lwd=3)
points(newdata$X,PredNew$lwr,type="p",col="black",pch="-",cex=2,lwd=3)
points(newdata$X,PredNew$upr,type="p",col="black",pch="-",cex=2,lwd=3)

# Salvare il grafico in formato pdf