En este ejercicio, generaremos datos simulados y luego utilizaremos regresión no lineal para realizar la regresión

set.seed(1)

#200 observaciones de una normal con media 0 y desviación estandar 1
X<-rnorm(200,0,1)
#vector de errores (suponiendo normalidad)
e<-rnorm(200,0,1)

#definimos los predictores
x1<-X
x2<-X^2
x3<-X^3

#creamos el vector de respuestas
Y<-3+2*x1-7*x2+4*x3+e
#creamos nuestro data frame
df <- data.frame(Y, X1 = x1, X2 = x2, X3 = x3)

head(df)
##            Y         X1         X2           X3
## 1 -1.5740095 -0.6264538 0.39244438 -0.245848275
## 2  4.8448592  0.1836433 0.03372487  0.006193347
## 3 -4.3065899 -0.8356286 0.69827518 -0.583498718
## 4  4.2846614  1.5952808 2.54492084  4.059863355
## 5  0.7568581  0.3295078 0.10857537  0.035776429
## 6 -3.0647072 -0.8204684 0.67316837 -0.552313364
#definimos nuestra muestra de entrenamiento y nuestra muestra test
entrenamiento <- Y[1:150]
test <- Y[151:200]

df_entrenamiento <- df[1:150, ]
df_test <- df[151:200, ]
library(ggplot2)
#Diagrama de dispersión muestra de entrenamiento: X1 frente a Y
ggplot(df_entrenamiento, aes(x = Y, y = X1)) + geom_point(color = "blue")+
  labs(title = "Diagrama de Dispersión X1 frente a Y", x = "Y", y = "X1")

Como era de esperarse, existe una clara correlación entre las variables X1 e Y. Notemos que sabiamos que Y = 3+2X1-7X12+4X13+e. Los errores no son lo suficientemente grandes para que esta relación no sea visible.

Regresión polinómica

#polinomio de grado 2
fit_pol_grado2 <- lm(Y~X1+I(X1^2),data=df_entrenamiento)
summary(fit_pol_grado2)
## 
## Call:
## lm(formula = Y ~ X1 + I(X1^2), data = df_entrenamiento)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -28.219  -4.110   0.286   3.769  29.098 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   2.6197     0.5909   4.433  1.8e-05 ***
## X1           11.3325     0.5341  21.219  < 2e-16 ***
## I(X1^2)      -6.0936     0.4279 -14.242  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.861 on 147 degrees of freedom
## Multiple R-squared:  0.802,  Adjusted R-squared:  0.7993 
## F-statistic: 297.7 on 2 and 147 DF,  p-value: < 2.2e-16
#polinomio de grado 3
fit_pol_grado3 <- lm(Y~X1+I(X1^2)+I(X1^3),data=df_entrenamiento)
summary(fit_pol_grado3)
## 
## Call:
## lm(formula = Y ~ X1 + I(X1^2) + I(X1^3), data = df_entrenamiento)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.91451 -0.52378  0.07186  0.68544  2.40386 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.99295    0.10152   29.48   <2e-16 ***
## X1           1.68817    0.16609   10.16   <2e-16 ***
## I(X1^2)     -6.93233    0.07439  -93.19   <2e-16 ***
## I(X1^3)      4.12021    0.05918   69.62   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.006 on 146 degrees of freedom
## Multiple R-squared:  0.9942, Adjusted R-squared:  0.9941 
## F-statistic:  8358 on 3 and 146 DF,  p-value: < 2.2e-16
#polinomio de grado 4
fit_pol_grado4 <- lm(Y~X1+I(X1^2)+I(X1^3)+I(X1^4),data=df_entrenamiento)
summary(fit_pol_grado4)
## 
## Call:
## lm(formula = Y ~ X1 + I(X1^2) + I(X1^3) + I(X1^4), data = df_entrenamiento)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.8181 -0.5164  0.0891  0.6682  2.2545 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.88660    0.11749  24.569   <2e-16 ***
## X1           1.62509    0.16874   9.631   <2e-16 ***
## I(X1^2)     -6.59993    0.20261 -32.575   <2e-16 ***
## I(X1^3)      4.15385    0.06178  67.236   <2e-16 ***
## I(X1^4)     -0.08652    0.04911  -1.762   0.0802 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9984 on 145 degrees of freedom
## Multiple R-squared:  0.9943, Adjusted R-squared:  0.9942 
## F-statistic:  6359 on 4 and 145 DF,  p-value: < 2.2e-16

Notemos que el polinomio que mejor ajusta es el de grado 4, sin embargo la diferencia con el R2adj del polinomio de grado 3 es mínima, por lo cual nos quedaremos con el de grado 3 para el resto del trabajo.

library(splines)

#spline con 3 cortes
#obtenemos los cuantiles
qua_3 <- as.numeric(quantile(df_entrenamiento$X1,probs=seq(0,1,length.out=5)))
qua_3
## [1] -2.21469989 -0.58498722 -0.04936932  0.59393497  2.40161776
#obtenemos los intervalos 
cut_3 <- cut(df_entrenamiento$X1,breaks=qua_3,include.lowest=T)
table(cut_3)
## cut_3
##   [-2.21,-0.585] (-0.585,-0.0494]  (-0.0494,0.594]      (0.594,2.4] 
##               38               37               37               38
fit_x1_Y_spl_3 <- lm(Y~bs(X1,knots=qua_3[-c(1,5)]),data=df_entrenamiento)
summary(fit_x1_Y_spl_3)
## 
## Call:
## lm(formula = Y ~ bs(X1, knots = qua_3[-c(1, 5)]), data = df_entrenamiento)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.82227 -0.50459  0.07244  0.65908  2.24897 
## 
## Coefficients:
##                                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                      -80.1690     0.8283  -96.79   <2e-16 ***
## bs(X1, knots = qua_3[-c(1, 5)])1  51.6397     1.4840   34.80   <2e-16 ***
## bs(X1, knots = qua_3[-c(1, 5)])2  78.1429     0.8582   91.05   <2e-16 ***
## bs(X1, knots = qua_3[-c(1, 5)])3  83.8699     0.9232   90.85   <2e-16 ***
## bs(X1, knots = qua_3[-c(1, 5)])4  81.4860     0.9638   84.55   <2e-16 ***
## bs(X1, knots = qua_3[-c(1, 5)])5  81.1954     1.2271   66.17   <2e-16 ***
## bs(X1, knots = qua_3[-c(1, 5)])6 103.5211     1.1431   90.56   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.005 on 143 degrees of freedom
## Multiple R-squared:  0.9943, Adjusted R-squared:  0.9941 
## F-statistic:  4183 on 6 and 143 DF,  p-value: < 2.2e-16
#hacemos muchos más cortes k=12

qua_12 <- as.numeric(quantile(df_entrenamiento$X1,probs=seq(0,1,length.out=14)))
qua_12
##  [1] -2.21469989 -1.28181977 -0.74385365 -0.61715888 -0.46193483 -0.27118973
##  [7] -0.10079005  0.01984328  0.30190813  0.47734330  0.66284951  0.91085139
## [13]  1.39831177  2.40161776
cut_12 <- cut(df_entrenamiento$X1,breaks=qua_12,include.lowest=T)
table(cut_12)
## cut_12
##   [-2.21,-1.28]  (-1.28,-0.744] (-0.744,-0.617] (-0.617,-0.462] (-0.462,-0.271] 
##              12              11              12              11              12 
## (-0.271,-0.101] (-0.101,0.0198]  (0.0198,0.302]   (0.302,0.477]   (0.477,0.663] 
##              11              12              11              12              11 
##   (0.663,0.911]     (0.911,1.4]       (1.4,2.4] 
##              12              11              12
#instalamos la libreria
library(npreg)
## Package 'npreg' version 1.1.0
## Type 'citation("npreg")' to cite this package.
fit_x1_Y_spl_sua_12 <- ss(x=df_entrenamiento$X1,y=df_entrenamiento$Y,
                              method="OCV",knots=qua_12[-c(1,14)])
summary(fit_x1_Y_spl_sua_12)
## 
## Call:
## ss(x = df_entrenamiento$X1, y = df_entrenamiento$Y, method = "OCV", 
##     knots = qua_12[-c(1, 14)])
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -3.28368 -0.54725  0.04326  0.74258  2.60648 
## 
## Approx. Signif. of Parametric Effects:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -7.143     0.1209  -59.08        0 ***
## x             99.691     1.1581   86.08        0 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
## 
## Approx. Signif. of Nonparametric Effects:
##               Df  Sum Sq  Mean Sq F value Pr(>F)    
## s(x)        9.45 11836.1 1252.501    1023      0 ***
## Residuals 138.55   169.6    1.224                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
## 
## Residual standard error: 1.106 on 138.5 degrees of freedom
## Multiple R-squared:  0.9934,    Adjusted R-squared:  0.9928
## F-statistic: 1980 on 10.45 and 138.5 DF,  p-value: <2e-16

Veamos un gráfico que comparar los tres ajustes principales

ggplot(df_entrenamiento,aes(x=X1,y=Y)) + 
  theme_light(base_size=15) +
  geom_point(size=2,color='grey') +
  geom_line(aes(y=predict(fit_pol_grado3)),col='green',size=.8) +
  geom_line(aes(y=predict(fit_x1_Y_spl_3)),col='blue',size=.8) + 
  geom_line(aes(y=predict(fit_x1_Y_spl_sua_12,x=X1)$y),col='red',size=.8)
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Comparemos el ECM con la muestra test

#Regresion polinómica
pred_pol <- predict(fit_pol_grado3,df_test)
ECMT_pol <- mean((df_test$Y-pred_pol)^2)
ECMT_pol
## [1] 1.064443
#regresion spline
pred_spl <- predict(fit_x1_Y_spl_3,df_test)
ECMT_spl <- mean((df_test$Y-pred_spl)^2)
ECMT_spl
## [1] 1.058095
#regresion spline suavizado
pred_spl_sua <- predict(fit_x1_Y_spl_sua_12,x=df_test$X1)$y
ECMT_spl_sua <- mean((df_test$Y-pred_spl_sua)^2)
ECMT_spl_sua
## [1] 1.220433

La regresión con spline de grado 3 tiene el menor error cuadrático medio.

color_1 <- "turquoise2"
ECMT_tod <- c(ECMT_pol,ECMT_spl,ECMT_spl_sua)
names(ECMT_tod) <- names(c(ECMT_pol,ECMT_spl,ECMT_spl_sua))
names(ECMT_tod)[1] <- "Polinómica"
names(ECMT_tod)[2] <- "Spline"
names(ECMT_tod)[3] <- "Spline suave"
names(ECMT_tod)
## [1] "Polinómica"   "Spline"       "Spline suave"
ECMT_df <- data.frame(x=names(ECMT_tod),y1=ECMT_tod)
ggplot(ECMT_df,aes(x=x,y=y1)) + 
  theme_light(base_size=15) +
  geom_segment(aes(x=x,xend=x,y=0,yend=y1),color=color_1,lwd=1.5) +
  geom_point(size=2,color=color_1) + 
  xlab("Método") + 
  ylab("ECMT")

Las regresiones tienen R2 ajustados muy similares, a pesar de que la regresión con spline tuvo el error cuadratico medio más bajo, las tres se adaptan muy bien a los datos. Siempre hay que recordar que el modelo de spline es el mejor modelo para esta miuestra de entrenamoiento y muestra test. Si estas cambiaran, es probable que el resultado cambie.