# Modelos de Regresión Múltiple.


# N.H. Prater  desarrolló una ecuación de regresión para
# estimar la producción de gasolina como una función de las
# propiedades de destilación de cierto tipo de petróleo crudo. Se
# identificaron cuatro variables de predicción: la graduación del
# petróleo crudo, $^o$API (x1); la presión de vapor del petróleo
# crudo, psi (x2); el punto de 10\% ASTM para el petróleo crudo,
# $^o$F (x3) y el punto final ASTM para la gasolina, $^o$F (x4). Los
# dos primeros miden la graduación y la presión de vapor del
# petróleo crudo. El punto de 10\% ASTM es la temperatura para la
# cual se ha evaporado cierta cantidad de líquido, y el punto final
# para la gasolina es la temperatura para la cual se ha evaporado
# todo el líquido. La variable respuesta (y) fue la cantidad de
# gasolina producida expresada como un porcentaje respecto al total
# de petróleo crudo. Los datos de laboratorio obtenidos por Prater
# se muestran en el archivo PRATER.DAT. Determinar una ecuación de
# regresión para la producción de gasolina como una función lineal
# de las propiedades de destilación de cierto tipo de petróleo crudo
# x1, x2, x3 y x4. ¿Podemos eliminar alguna de las variables
# predictoras del modelo?

dat<-read.table("prater.dat",head=T)

res<-lm(y~.,data=dat)
summary(res)

# 1. Métodos de selección de variables.

# 1.1 Métodos paso a paso.

# A. Método Backward.
# 1. Partimos del modelo con todas las variables predictoras.
# 2. Quitamos la variable predictora con el mayor p-valor que supere un alfa-crítico (p.ej. 0.1,0.15)
# 3. Reajustamos el modelo y volvemos al punto 2.
# 4. Paramos cuando todos los p-valores sean menores que el alfa-crítico.

res1<-lm(y~.-x2, data=dat)


# B. Método Forward.
# 1. Partimos del modelo básico sin variables predictoras.
# 2. Para todas las variables predictoras que no están en el modelo calculamos su p-valor si 
#    son añadidas al modelo. Elegimos aquella variable con el menor p-valor menor que el 
#    alfa-crítico (p. ej. 0.1,0.15)
# 3. Continuamos hasta que ninguna variable pueda ser añadida al modelo.

coefficients(summary(lm(y~x1,data=dat)))[,4][2]
coefficients(summary(lm(y~x2,data=dat)))[,4][2]
coefficients(summary(lm(y~x3,data=dat)))[,4][2]
coefficients(summary(lm(y~x4,data=dat)))[,4][2]

coefficients(summary(lm(y~x4+x1,data=dat)))[,4][3]
coefficients(summary(lm(y~x4+x2,data=dat)))[,4][3]
coefficients(summary(lm(y~x4+x3,data=dat)))[,4][3]

coefficients(summary(lm(y~x4+x3+x1,data=dat)))[,4][4]
coefficients(summary(lm(y~x4+x3+x2,data=dat)))[,4][4]

coefficients(summary(lm(y~x4+x3+x1+x2,data=dat)))[,4][5]


# 1.2 Métodos globales.

# Si tenemos p variables predictoras, entonces hay 2^p modelos posibles.
# Elegiremos entre el mejor según un criterio.

# A. Criterio de Información de Akaike (AIC).
# AIC=-2(log-verosimilitud)+2(p+1)
# Deviance=-2(log-verosimilitud)=n(log(SCR/n))     SCR=Suma de Cuadrados Residual
# Elegimos el modelo de menor AIC

step(res)


# B. Criterio R^2-ajustado.
# Elegimos el modelo de mayor R^2-ajustado.

summary(res)$adj.r.squared

summary(lm(y~.-x1,data=dat))$adj.r.squared
summary(lm(y~.-x2,data=dat))$adj.r.squared
summary(lm(y~.-x3,data=dat))$adj.r.squared
summary(lm(y~.-x4,data=dat))$adj.r.squared


# 2. Métodos de diagnóstico.

layout(matrix(1:4,2))
plot(res)

# Normalidad.
library(ctest)
shapiro.test(rstandard(res))

# Incorrelación en las observaciones.
# Test de Durbin-Watson

library(car)
durbin.watson(res)


# 3. Predicciones.

predict(res,data.frame(x1=32,x2=6,x3=240,x4=350),se.fit=T,interval="confidence")
predict(res,data.frame(x1=32,x2=6,x3=240,x4=350),se.fit=T,interval="prediction")


# 4. Multicolinealidad.

# Se realiza un experimento para determinar el calor
# desarrollado en la fabricación de cemento en función de los
# porcentajes de 4 compuestos activos que se utilizan en la
# fabricación. Para ello se elige una muestra de 13 cementos,
# midiéndose el calor desarrollado de cals/g (y) y los porcentajes
# de los compuestos referidos ( x1, x2, x3, x4). Los datos se
# encuentran en el archivo CEMENTO.DAT. Calcula la ecuación pedida,
# investigando la bondad del ajuste y eliminando de modo razonado
# alguna de las variables predictoras del modelo, si ello fuera
# necesario.

datos<-read.table("cemento.dat",head=T)
res<-lm(y~.,data=datos)
summary(res)

cor(datos)

# Factor de inflacción de la varianza: 1/(1-Ri^2)
# Ri^2 coeficiente de determinación entre xi y el resto de variables explicativas.

library(car)
vif(res)

step(res)

res1<-lm(y~.-x3,data=datos)
summary(res1)

vif(res1)

res2<-lm(y~.-x3-x4,data=datos)
summary(res2)

vif(res2)

plot(res2)


# 5. Regresión polinómica.

# Una compañía observa que la demanda de uno de sus productos cambió
# debido a una variación rápida de su precio por unidad. Supóngase que la
# demanda Y del producto se observa en una región en particular sobre un
# intervalo bastante amplio de precios x. Dados los datos que se encuentran en
# COMPANY.DAT, determinar el grado de un polinomio que mejor ajuste estos
# datos.

datos<-read.table("company.dat",head=T)
attach(datos)

res<-lm(unidades~dolares,data=datos)
summary(res)

layout(matrix(1:4,2))
plot(res)


res1<-lm(unidades~dolares+I(dolares^2),data=datos)
summary(res1)

plot(res1)


res2<-lm(unidades~dolares+I(dolares^2)+I(dolares^3),data=datos)
summary(res2)

plot(res2)


res3<-lm(unidades~dolares+I(dolares^2)+I(dolares^3)+I(dolares^4),data=datos)
summary(res3)


par(mfrow=c(1,1))

plot(dolares,unidades)
abline(res)

z<-seq(min(dolares),max(dolares),len=1000)
lines(z,1330.407-155.467*z+4.866*z^2,col=2)
lines(z,as.vector(res2$coefficients%*%rbind(1,z,z^2,z^3)),col=4)


a<-predict(res2,data.frame(dolares=z),interval="confidence")
lines(z,a[,3],col=4,lty=2)
lines(z,a[,2],col=4,lty=2)

b<-predict(res2,data.frame(dolares=z),interval="prediction")
lines(z,b[,3],col=4,lty=3)
lines(z,b[,2],col=4,lty=3)

plot(dolares,unidades,ylim=c(min(b[,2]),max(b[,3])))