x1:Conversion del n-heptano en acetileno(%)

x2:Temperatura del reactor

x3:Reaccion(molar)de H2 a heptano

x4:Tiempo de contacto(seg)

x1<-c(49,50.2,50.5,48.5,47.5,44.5,28,31.5,34.5,35,38,38.5,15,17,20.5,29.5) #Variable respuesta
x2<-c(1300,1300,1300,1300,1300,1300,1200,1200,1200,1200,1200,1200,1100,1100,1100,1100)
x3<-c(7.5,9,11,13.5,17,23,5.3,7.5,11,13.5,17,23,5.3,7.5,11,17)
x4<-c(0.012,0.012,0.0115,0.013,0.0135,0.012,0.04,0.038,0.032,0.026,0.034,0.041,0.084,0.098,0.092,0.086)
datos=data.frame(x1,x2,x3,x4) #x2,x3,x4, son variables explicativas
datos#El data frame corresponde a los datos sobre el acetileno(ejemplo 10.1,introduccion al analisis de regresion lineal)
##      x1   x2   x3     x4
## 1  49.0 1300  7.5 0.0120
## 2  50.2 1300  9.0 0.0120
## 3  50.5 1300 11.0 0.0115
## 4  48.5 1300 13.5 0.0130
## 5  47.5 1300 17.0 0.0135
## 6  44.5 1300 23.0 0.0120
## 7  28.0 1200  5.3 0.0400
## 8  31.5 1200  7.5 0.0380
## 9  34.5 1200 11.0 0.0320
## 10 35.0 1200 13.5 0.0260
## 11 38.0 1200 17.0 0.0340
## 12 38.5 1200 23.0 0.0410
## 13 15.0 1100  5.3 0.0840
## 14 17.0 1100  7.5 0.0980
## 15 20.5 1100 11.0 0.0920
## 16 29.5 1100 17.0 0.0860
z<-aov(x1~x2+x3+x4)#Modelo ajustado(regresion lineal clasica)
print(coef(z))     # Estimacion de los parametros
##  (Intercept)           x2           x3           x4 
## -121.2696214    0.1268539    0.3481576  -19.0216969
print(summary(z))  #ANOVA
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## x2           1 1896.7  1896.7 133.655 7.32e-08 ***
## x3           1   56.3    56.3   3.967   0.0696 .  
## x4           1    0.4     0.4   0.031   0.8630    
## Residuals   12  170.3    14.2                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
print(vif(z))      #Detectar multicolinealidad.Se puede observar que efectivamente hay multicolinealidad,valores de VIF mayores a 10
##        x2        x3        x4 
## 12.225045  1.061838 12.324964
m<-model.matrix(x1~.,-1,data=datos)#Eliminar el intercepto que contiene unos,y se especifica que x1,es la variable de respuesta y datos es el data frame creado anteriormente
m
##    (Intercept)   x2   x3     x4
## 1            1 1300  7.5 0.0120
## 2            1 1300  9.0 0.0120
## 3            1 1300 11.0 0.0115
## 4            1 1300 13.5 0.0130
## 5            1 1300 17.0 0.0135
## 6            1 1300 23.0 0.0120
## 7            1 1200  5.3 0.0400
## 8            1 1200  7.5 0.0380
## 9            1 1200 11.0 0.0320
## 10           1 1200 13.5 0.0260
## 11           1 1200 17.0 0.0340
## 12           1 1200 23.0 0.0410
## 13           1 1100  5.3 0.0840
## 14           1 1100  7.5 0.0980
## 15           1 1100 11.0 0.0920
## 16           1 1100 17.0 0.0860
## attr(,"assign")
## [1] 0 1 2 3
n<-x1#Variable de respuesta

#Regresion ridge

fit.ridge=glmnet(m,n,alpha=0) #Como se realiza regresion ridge,alpha debe ser igual a cero.Donde m es una matriz y n es la variable respuesta
fit.ridge#Por defecto encuentra 100 modelos y el lambda asociado a cada uno
## 
## Call:  glmnet(x = m, y = n, alpha = 0) 
## 
##     Df  %Dev  Lambda
## 1    3  0.00 10890.0
## 2    3  0.43  9920.0
## 3    3  0.47  9039.0
## 4    3  0.52  8236.0
## 5    3  0.57  7504.0
## 6    3  0.63  6838.0
## 7    3  0.69  6230.0
## 8    3  0.75  5677.0
## 9    3  0.83  5173.0
## 10   3  0.91  4713.0
## 11   3  0.99  4294.0
## 12   3  1.09  3913.0
## 13   3  1.19  3565.0
## 14   3  1.31  3249.0
## 15   3  1.43  2960.0
## 16   3  1.57  2697.0
## 17   3  1.72  2457.0
## 18   3  1.89  2239.0
## 19   3  2.07  2040.0
## 20   3  2.27  1859.0
## 21   3  2.49  1694.0
## 22   3  2.72  1543.0
## 23   3  2.98  1406.0
## 24   3  3.26  1281.0
## 25   3  3.57  1167.0
## 26   3  3.91  1064.0
## 27   3  4.28   969.2
## 28   3  4.68   883.1
## 29   3  5.11   804.7
## 30   3  5.59   733.2
## 31   3  6.11   668.1
## 32   3  6.67   608.7
## 33   3  7.28   554.6
## 34   3  7.94   505.4
## 35   3  8.66   460.5
## 36   3  9.44   419.6
## 37   3 10.28   382.3
## 38   3 11.18   348.3
## 39   3 12.16   317.4
## 40   3 13.22   289.2
## 41   3 14.35   263.5
## 42   3 15.56   240.1
## 43   3 16.86   218.8
## 44   3 18.25   199.3
## 45   3 19.72   181.6
## 46   3 21.29   165.5
## 47   3 22.96   150.8
## 48   3 24.72   137.4
## 49   3 26.57   125.2
## 50   3 28.51   114.1
## 51   3 30.55   103.9
## 52   3 32.67    94.7
## 53   3 34.87    86.3
## 54   3 37.14    78.6
## 55   3 39.48    71.6
## 56   3 41.87    65.3
## 57   3 44.31    59.5
## 58   3 46.78    54.2
## 59   3 49.27    49.4
## 60   3 51.77    45.0
## 61   3 54.26    41.0
## 62   3 56.73    37.4
## 63   3 59.16    34.0
## 64   3 61.54    31.0
## 65   3 63.86    28.2
## 66   3 66.10    25.7
## 67   3 68.24    23.5
## 68   3 70.29    21.4
## 69   3 72.24    19.5
## 70   3 74.07    17.7
## 71   3 75.78    16.2
## 72   3 77.38    14.7
## 73   3 78.85    13.4
## 74   3 80.21    12.2
## 75   3 81.45    11.1
## 76   3 82.58    10.2
## 77   3 83.60     9.3
## 78   3 84.51     8.4
## 79   3 85.34     7.7
## 80   3 86.07     7.0
## 81   3 86.72     6.4
## 82   3 87.30     5.8
## 83   3 87.81     5.3
## 84   3 88.26     4.8
## 85   3 88.65     4.4
## 86   3 89.00     4.0
## 87   3 89.30     3.6
## 88   3 89.57     3.3
## 89   3 89.80     3.0
## 90   3 90.00     2.8
## 91   3 90.18     2.5
## 92   3 90.34     2.3
## 93   3 90.48     2.1
## 94   3 90.61     1.9
## 95   3 90.72     1.7
## 96   3 90.82     1.6
## 97   3 90.91     1.4
## 98   3 91.00     1.3
## 99   3 91.07     1.2
## 100  3 91.14     1.1
plot(fit.ridge,"lambda",label=TRUE)

#Validacion cruzada:Para encontrar el mejor lambda,es decir el mejor K.
cv.ridge=cv.glmnet(m,n,alpha=0)
cv.ridge
## 
## Call:  cv.glmnet(x = m, y = n, alpha = 0) 
## 
## Measure: Mean-Squared Error 
## 
##     Lambda Index Measure    SE Nonzero
## min  1.089   100   20.75 6.675       3
## 1se  7.681    79   27.40 6.685       3
plot(cv.ridge)
abline(v=log(fit.ridge$lambda[100]),col="blue",lwd=4,lty=3)#Permite observar cual deberia ser el valor de lambda,o en que rango se encuentra.

mejorlambda<-cv.ridge$lambda.min # El mejor lambda, es decir donde los coeficientes se estabilizan, el RMSE disminuye, y aumenta el coeficiente de determinacion.
mejorlambda
## [1] 1.088771
log(mejorlambda)
## [1] 0.08504987
coef(fit.ridge)[,which(fit.ridge$lambda==mejorlambda)]#Modelo ajustado asociado con el mejor lambda
##   (Intercept)   (Intercept)            x2            x3            x4 
##  -63.07470183    0.00000000    0.08235242    0.32670125 -117.50061031
#Datos de entrenamiento

split = sample.split(datos, SplitRatio = 0.7)#Del dataframe escojo el 70 por ciento de las observaciones,sin reemplazo.
test1 = subset(datos, split == FALSE)
test1
##      x1   x2   x3    x4
## 2  50.2 1300  9.0 0.012
## 4  48.5 1300 13.5 0.013
## 6  44.5 1300 23.0 0.012
## 8  31.5 1200  7.5 0.038
## 10 35.0 1200 13.5 0.026
## 12 38.5 1200 23.0 0.041
## 14 17.0 1100  7.5 0.098
## 16 29.5 1100 17.0 0.086
#Predecir con el modelo ridge

x.test1=model.matrix(x1~.,-1,data=test1)
pred=predict(fit.ridge, s=mejorlambda,newx = x.test1)
pred #Predicion de la variable de respuesta x1
##          s1
## 2  45.51375
## 4  46.86640
## 6  50.08756
## 8  33.73344
## 10 37.10365
## 12 38.44480
## 14 18.44816
## 16 22.96183
data.frame(RMSE=RMSE(pred,test1$x1),Rsquare = R2(pred,test1$x1))#RMSE y R^2,del mejor modelo,es decir con el lambda mas pequeño.
##       RMSE        s1
## 1 3.709986 0.8788269