#Ejercicio 1:
#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,  grados API (x1); 
#la presión de vapor del petróleo crudo, psi (x2); 
#el punto de 10% ASTM para el petróleo crudo, grados Fahrenheit (x3) y el punto final ASTM para la gasolina, 
#grados Farenheit (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?
  
  # Importamos y chequeamos los datos.
  
  prater<-read.table("./datos/prater.dat",head=T)
  colnames(prater)
  summary(prater)
  str(prater)
  
  # Ajuste del modelo
  
  ml1<-lm(y~x1+x2+x3+x4,data=prater)
  summary(ml1)
  attributes(ml1)
  ml1$coefficients
  ml1$fitted.values
  ml1$residuals
  attributes(summary(ml1))
  
  
  # Diagnóstico del modelo
  
  # gráficos de diagnóstico
  layout(matrix(1:4,2))
  plot(ml1)
  
  #linealidad
  #con gráficos
   
  # normalidad
  shapiro.test(rstandard(ml1))
  
  # homocedasticidad
  library(lmtest)
  bptest(ml1)
  
  # errores incorrelados
  library(lmtest)
  dwtest(ml1)
  bgtest(ml1,order=2)
  
 
  
  # multicolinealidad
  library(car)
  vif(ml1)
  
  #si todo ha funcionado, interpretamos los resultados del modelo
  summary(ml1)
  # Intervalos de confianza
  confint(ml1,level=0.99)
  
  # Selección de variables
  
  # basándonos en AIC
  library(MASS)
  ml0<-lm(y~1,data=prater)
  stepAIC(ml0,scope=list(upper=ml1,lower=ml0),direction="forward")
  stepAIC(ml1,scope=list(upper=ml1,lower=ml0),direction="backward")
  
  # basándonos en p-valor
  library(olsrr)
  ols_step_backward_p(ml1, p_val=0.1,print_plot=T)
  ols_step_forward_p(ml1, p_val = 0.1)
  ols_step_forward_p(ml1, p_val = 0.1,details=T)
  
  # predicciones
  ml2<-lm(y~x1+x2+x4,data=prater)
  predict(ml2,data.frame(x1=40, x2=5, x3=270,  x4=350), interval="confidence")
  predict(ml2,data.frame(x1=40, x2=5, x3=270,  x4=350), interval="prediction")
  
#  Ejercicio 2:
#    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.
  
  
  cemento<-read.table("./datos/cemento.dat",head=T)
  colnames(cemento)
  ml3<-lm(y~.,data=cemento)
  summary(ml3)
  
  # multicolinealidad
  library(car)
  vif(ml3)
  
  
  # selección de variables
  library(olsrr)
  ols_step_backward_p(ml3,prem=0.1,details=T)
  ols_step_backward_aic(ml3,details=T)
  ols_step_forward_p(ml3,prem=0.1,details=T)
  ols_step_forward_aic(ml3,details=T)
  
  
  # Regresión polinómica
  
#  Ejercicio 3:
#    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. Se observa la demanda del producto
#  (unidades) en una región en particular sobre un intervalo bastante amplio de precios
#  (dolares). Los datos que se encuentran en el archivo company.dat.
  
#  Se trata de realizar un análisis de dichos datos, indicando una ecuación que marque la
#  relación entre el precio del producto por unidad y la demanda del mismo.
  
  
  datos<-read.table("./datos/company.dat",head=T)
  
  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)
  
  
  attach(datos)
  
  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])))