Regresion No Lineal:

library("MASS")
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ✖ dplyr::select() masks MASS::select()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(caret)
## Cargando paquete requerido: lattice
## 
## Adjuntando el paquete: 'caret'
## 
## The following object is masked from 'package:purrr':
## 
##     lift
library(ggplot2)
library(splines)
library(mgcv)
## Cargando paquete requerido: nlme
## 
## Adjuntando el paquete: 'nlme'
## 
## The following object is masked from 'package:dplyr':
## 
##     collapse
## 
## This is mgcv 1.9-1. For overview type 'help("mgcv-package")'.
data("Boston")
glimpse(Boston)
## Rows: 506
## Columns: 14
## $ crim    <dbl> 0.00632, 0.02731, 0.02729, 0.03237, 0.06905, 0.02985, 0.08829,…
## $ zn      <dbl> 18.0, 0.0, 0.0, 0.0, 0.0, 0.0, 12.5, 12.5, 12.5, 12.5, 12.5, 1…
## $ indus   <dbl> 2.31, 7.07, 7.07, 2.18, 2.18, 2.18, 7.87, 7.87, 7.87, 7.87, 7.…
## $ chas    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ nox     <dbl> 0.538, 0.469, 0.469, 0.458, 0.458, 0.458, 0.524, 0.524, 0.524,…
## $ rm      <dbl> 6.575, 6.421, 7.185, 6.998, 7.147, 6.430, 6.012, 6.172, 5.631,…
## $ age     <dbl> 65.2, 78.9, 61.1, 45.8, 54.2, 58.7, 66.6, 96.1, 100.0, 85.9, 9…
## $ dis     <dbl> 4.0900, 4.9671, 4.9671, 6.0622, 6.0622, 6.0622, 5.5605, 5.9505…
## $ rad     <int> 1, 2, 2, 3, 3, 3, 5, 5, 5, 5, 5, 5, 5, 4, 4, 4, 4, 4, 4, 4, 4,…
## $ tax     <dbl> 296, 242, 242, 222, 222, 222, 311, 311, 311, 311, 311, 311, 31…
## $ ptratio <dbl> 15.3, 17.8, 17.8, 18.7, 18.7, 18.7, 15.2, 15.2, 15.2, 15.2, 15…
## $ black   <dbl> 396.90, 396.90, 392.83, 394.63, 396.90, 394.12, 395.60, 396.90…
## $ lstat   <dbl> 4.98, 9.14, 4.03, 2.94, 5.33, 5.21, 12.43, 19.15, 29.93, 17.10…
## $ medv    <dbl> 24.0, 21.6, 34.7, 33.4, 36.2, 28.7, 22.9, 27.1, 16.5, 18.9, 15…

The Boston Housing Dataset

The Boston Housing Dataset is a derived from information collected by the U.S. Census Service concerning housing in the area of Boston MA. The following describes the dataset columns:

CRIM - per capita crime rate by town ZN - proportion of residential land zoned for lots over 25,000 sq.ft INDUS - proportion of non-retail business acres per town CHAS - Charles River dummy variable (1 if tract bounds river; 0 otherwise) NOX - nitric oxides concentration (parts per 10 million) RM - average number of rooms per dwelling AGE - proportion of owner-occupied units built prior to 1940 DIS - weighted distances to five Boston employment centres RAD - index of accessibility to radial highways TAX - full-value property-tax rate per 10,000 dollars PTRATIO - pupil-teacher ratio by town B - 1000(Bk - 0.63)^2 where Bk is the proportion of blacks by town LSTAT - % lower status of the population MEDV - Median value of owner-occupied homes in $1000’s

plot(Boston$lstat,Boston$medv)

set.seed(123)
training.samples<-Boston$medv %>%
   createDataPartition(p=0.8,list=FALSE)
train.data<-Boston[training.samples, ]
test.data<-Boston[-training.samples, ]

REGRESION LINEAL

model<-lm(medv~lstat,data=train.data)
predictions<-model %>% predict(test.data)
data.frame(
  RMSE=RMSE(predictions,test.data$medv),
  R2=R2(predictions,test.data$medv)
)
##       RMSE       R2
## 1 6.503817 0.513163
ggplot(train.data, aes(lstat,medv))+
  geom_point()+
  stat_smooth(method=lm, formula=y~x)

REGRESION POLINOMICA

El modelo sin variables vs el modelo con variables:

En regresion lineal hacemos 2 pruebas de hipotesis: 1. Para el modelo en general-> H0: modelo = modelo sin variables 2. Para cada coeficiente del modelo

A continuacion vemos que p-value: < 2.2e-16, por lo que rechazamos la H0: modelo = modelo sin variables

Grado 2

model<-lm(medv~poly(lstat,2),data=train.data)
summary(model)
## 
## Call:
## lm(formula = medv ~ poly(lstat, 2), data = train.data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -15.3654  -3.8250  -0.6439   2.2733  25.2922 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       22.5106     0.2727   82.56   <2e-16 ***
## poly(lstat, 2)1 -137.2577     5.5006  -24.95   <2e-16 ***
## poly(lstat, 2)2   55.3383     5.5006   10.06   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.501 on 404 degrees of freedom
## Multiple R-squared:  0.6418, Adjusted R-squared:   0.64 
## F-statistic: 361.9 on 2 and 404 DF,  p-value: < 2.2e-16

Grado 6

model<-lm(medv~poly(lstat,6),data=train.data)
summary(model)
## 
## Call:
## lm(formula = medv ~ poly(lstat, 6), data = train.data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -13.1962  -3.1527  -0.7655   2.0404  26.7661 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       22.5106     0.2572  87.533  < 2e-16 ***
## poly(lstat, 6)1 -137.2577     5.1881 -26.456  < 2e-16 ***
## poly(lstat, 6)2   55.3383     5.1881  10.666  < 2e-16 ***
## poly(lstat, 6)3  -24.4801     5.1881  -4.718 3.29e-06 ***
## poly(lstat, 6)4   22.9043     5.1881   4.415 1.30e-05 ***
## poly(lstat, 6)5  -15.3276     5.1881  -2.954  0.00332 ** 
## poly(lstat, 6)6    9.9068     5.1881   1.910  0.05691 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.188 on 400 degrees of freedom
## Multiple R-squared:  0.6845, Adjusted R-squared:  0.6798 
## F-statistic: 144.6 on 6 and 400 DF,  p-value: < 2.2e-16

Grado 3

model<-lm(medv~poly(lstat,3),data=train.data)
summary(model)
## 
## Call:
## lm(formula = medv ~ poly(lstat, 3), data = train.data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -14.6139  -3.7269  -0.5313   2.3783  26.3284 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       22.5106     0.2662  84.558  < 2e-16 ***
## poly(lstat, 3)1 -137.2577     5.3707 -25.557  < 2e-16 ***
## poly(lstat, 3)2   55.3383     5.3707  10.304  < 2e-16 ***
## poly(lstat, 3)3  -24.4801     5.3707  -4.558 6.86e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.371 on 403 degrees of freedom
## Multiple R-squared:  0.6594, Adjusted R-squared:  0.6568 
## F-statistic:   260 on 3 and 403 DF,  p-value: < 2.2e-16
model<-lm(medv~poly(lstat,5),data=train.data)
predictions<-model %>% predict(test.data)
data.frame(
  RMSE=RMSE(predictions,test.data$medv),
  R2=R2(predictions,test.data$medv)
)
##       RMSE        R2
## 1 5.270374 0.6829474
ggplot(train.data, aes(lstat,medv))+
  geom_point()+
  stat_smooth(method=lm,formula=y~poly(x,5))

TRANSFORMACION LOGARITMICA

model<-lm(medv~log(lstat),data=train.data)
predictions<-model %>% predict(test.data)
data.frame(
  RMSE=RMSE(predictions,test.data$medv),
  R2=R2(predictions,test.data$medv)
)
##       RMSE        R2
## 1 5.467124 0.6570091
ggplot(train.data, aes(lstat,medv))+
  geom_point()+
  stat_smooth(method=lm,formula=y~log(x))

REGRESION SPLINE

Otra forma de tratar las relaciones no lineales es la regresión spline. Aquí, como en la regresión polinómica, el modelo se crea tomando grados. Además, se obtienen los valores de los nudos (knots) y se puede crear la línea de regresión combinando adecuadamente los nudos con el grado especificado. En la regresión spline, se suele establecer un modelo con un 3er grado. Ahora, vamos a comprobar los resultados aplicando este método al conjunto de datos de entrenamiento:

NOSOTROS DEBEMOS ESCOGER DONDE VAMOS A PONER LOS ‘NUDOS’ DEL SPLINE

¿donde?

busco el percentil-10 de los datos en mis X. busco el percentil-50 busco el percentil-90

[0-10][10-50][50-90][90-100] <- 3 nudos

Acá ‘bs’ indica splines

knots<-quantile(train.data$lstat,p=c(0.10,0.50,0.90))
model<-lm(medv~bs(lstat, knots=knots),data=train.data)
predictions<-model %>% predict(test.data)
data.frame(
  RMSE=RMSE(predictions,test.data$medv),
  R2=R2(predictions,test.data$medv)
)
##       RMSE        R2
## 1 5.255264 0.6839899
ggplot(train.data, aes(lstat,medv))+
  geom_point()+
  stat_smooth(method=lm,formula=y~splines::bs(x,df=3))

MODELOS ADITIVOS GENERALIZADOS

En la regresión polinómica, se determinan los grados, y en la regresión spline, se determinan los nUdos y se crean los modelos. En los modelos aditivos generalizados, no ocurre lo mismo. La estimación de los datos puede obtenerse automáticamente. Para los modelos aditivos generalizados, debe utilizarse la función «gam» en R:

https://www.stat.cmu.edu/~ryantibs/advmethods/notes/smoothspline.pdf

model<-gam(medv~ s(lstat),data=train.data)
predictions<-model %>% predict(test.data)
data.frame(
  RMSE=RMSE(predictions,test.data$medv),
  R2=R2(predictions,test.data$medv)
)
##       RMSE        R2
## 1 5.318856 0.6760512
ggplot(train.data, aes(lstat,medv))+
  geom_point()+
  stat_smooth(method=gam,formula=y~s(x))

CONCLUSIÓN

Hemos analizado métodos alternativos que pueden utilizarse para modelizar relaciones no lineales y los hemos comparado en función de los criterios RMSE y coeficiente de determinación. Como resultado, podemos decir que la regresión polinómica, la regresión spline y los modelos aditivos generalizados pueden preferirse para este conjunto de datos.