Comenzamos por cargar los datos

library(latex2exp)
Datos = read.csv("ejemplo4.csv")

Luego, para usar funciones del paquete boot se ajusta usando la función glm en vez de utilizar la función para el ajuste del modelo lm, donde el paquete boot tiene funciones que nos permiten realizar la variación cruzada. Por default glm coincide con el caso normal de la lm pero tiene más argumentos para funcionar con modelos lineales generalizados.

library(boot)
fitglm=glm(y~x, data=Datos)
summary(fitglm)
## 
## Call:
## glm(formula = y ~ x, data = Datos)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -5.7847  -2.8954  -0.9503   2.3287   7.8916  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -9.9928     0.8039  -12.43   <2e-16 ***
## x             5.5151     0.1296   42.57   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 13.84727)
## 
##     Null deviance: 26451  on 99  degrees of freedom
## Residual deviance:  1357  on 98  degrees of freedom
## AIC: 550.58
## 
## Number of Fisher Scoring iterations: 2

en el cual podemos notar la similitud que hay entre la información que nos arrojaba el lm y lo que nos arroja el glm, así como las diferencias pues en este caso no se muestra la información de la prueba \(F\), el valor de \(R^{2}\), etcétera.

Posteriormente, recordemos que en estos datos no hay linealidad por lo que procediamos a transformar

# Transformamos los datos
Datos$Xprima=Datos$x^2

# Graficamos
par(mfrow=c(1,2)) 
par(mar=c(4, 5, 1, 1))

# Datos no transformados
plot(Datos$x, Datos$y, xlab = TeX("$x$"), ylab=TeX("$y$") )

# Datos transformados
plot(Datos$Xprima, Datos$y, xlab = TeX("$x^2$"), ylab=TeX("$y$") )

y volvemos a realizar el ajuste del modelo usando glm

fit2glm=glm(y~Xprima, data=Datos)
summary(fit2glm)
## 
## Call:
## glm(formula = y ~ Xprima, data = Datos)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -2.36940  -0.59665  -0.02778   0.58044   2.09823  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 1.056729   0.142366   7.423 4.23e-11 ***
## Xprima      0.500875   0.002829 177.078  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.8409146)
## 
##     Null deviance: 26450.68  on 99  degrees of freedom
## Residual deviance:    82.41  on 98  degrees of freedom
## AIC: 270.44
## 
## Number of Fisher Scoring iterations: 2

Así, con los modelos obtenidos en fitglm y fit2glm efectuaremos la comparación.

Método de validación cruzada

library(caret)
# colocamos una semilla puesto que realizaremos particiones aleatorias
set.seed(1234)

La función createFolds() del paquete caret ayuda a obtener la muestra en cada una de las \(K\) divisiones.

# Agregamos una nueva columna al dataframe con puros unos y la convertimos
# a una variable de tipo factor que nos servirá para la función createFolds
# para indicar cuántos elementos tenemos
Datos$unos=1
Datos$unos=factor(Datos$unos)

# Hacemos las diviones aleatorias en 10 particiones
Folds=createFolds(Datos$unos, k = 10, list = FALSE)
Folds
##   [1]  8 10  2  9  5  8  6  4  6 10 10  9  8  4  6  2  3  4  1 10  2  7  6  6  7
##  [26]  1  8  3  1  2  2  3  2  4  9  9  1  6  7  9  1  5  7  8  6  7  8  1 10  6
##  [51]  8  5  2  4  3 10 10  5  5  6  9  5  3  3 10  7  2  7  1  8  7  2  9  4  1
##  [76]  3  5  5  7  6  5  3  4  3  9  8  7  8  4  5  4  3  4  1  9  1  9  2 10 10

donde podemos ver 10 datos asignados a la partición 1, 10 datos a la partición 2, etcétera. Lo anterior puede verse en

table(Folds)
## Folds
##  1  2  3  4  5  6  7  8  9 10 
## 10 10 10 10 10 10 10 10 10 10

Una vez obtenido lo anterior, realizaremos un bucle for para ir probando los conjuntos de datos de las particiones como el conjunto para el test. Comenzamos definiendo las matrices

SE1=matrix(0, 10,1)
SE2=matrix(0, 10,1)

en las que almacenaremos la suma de los errores al cuadrado correspodientes al modelo 1 (fitglm) y 2 (fit2glm). Analicemos primero una iteración antes de correr el for completo. Consideraremos \(kj=1\).

  1. ``Datos.train=Datos[(Folds==kj)==FALSE,]```: realizamos el filtro de nuestros datos para aquellos que no tengan asignados el valor de 1. Éstos serán nuestros datos tomados para “entrenar” el modelo.

  2. Datos.test=Datos[Folds==kj,]: realizamos el filtro de nuestros datos para la partición 1, es decir, seleccionamos sólo aquellos datos que tengan asignados el valor de uno. Éstos serán nuestros datos para probar el modelo.

  3. fitkj1=lm(y~x, data=Datos.train): ajustamos el modelos sólo para los datos del entrenamiento.

  4. predkj1=predict(fitkj1, Datos.test): utilizando el objeto lm obtenido en el punto anterior realizamos la predicción correspondiente con los datos del test.

  5. SE1[kj,1]=sum((Datos.test$y-predkj1)^2): finalmente guardamos en la matriz SE1 la suma de los errores al cuadrado, donde Datos.test$y-predkj1 representa la diferencia entre el valor “real” de \(y\) en el dataframe Datos.test con el valor predecido por predkj1

Repetimos lo mismo para el segundo modelo

  # En este caso usamos los valores de la x transformada
  fitkj2=lm(y~Xprima, data=Datos.train)
  
  predkj2=predict(fitkj2, Datos.test)
  SE2[kj,1]=sum((Datos.test$y-predkj2)^2)

Por consiguiente

for(kj in 1:10){
  Datos.train=Datos[(Folds==kj)==FALSE,]
  Datos.test=Datos[Folds==kj,]
  fitkj1=lm(y~x, data=Datos.train)
  predkj1=predict(fitkj1, Datos.test)
  SE1[kj,1]=sum((Datos.test$y-predkj1)^2)
  
  fitkj2=lm(y~Xprima, data=Datos.train)
  predkj2=predict(fitkj2, Datos.test)
  SE2[kj,1]=sum((Datos.test$y-predkj2)^2)  
}

y vemos las sumas de los errores al cuadrao calculados para cada partición

SE1
##            [,1]
##  [1,] 147.23557
##  [2,] 134.96619
##  [3,]  78.89767
##  [4,] 102.69749
##  [5,] 133.67675
##  [6,] 184.89389
##  [7,] 116.90607
##  [8,] 172.39307
##  [9,] 134.02468
## [10,] 213.01157
SE2
##            [,1]
##  [1,]  8.939381
##  [2,] 16.977388
##  [3,]  3.288027
##  [4,] 11.122032
##  [5,]  3.698293
##  [6,]  4.627293
##  [7,]  5.718196
##  [8,] 13.257321
##  [9,]  9.247512
## [10,] 11.176527

Procedemos después a calcular el promedio, donde sumaremos los valores respectivos de la matriz y lo dividiremos entre 100. Recordemos que cada entrada de la matriz al ser promediada debía ser dividida entre 10, pues cada partición tenía 10 datos, y como queremos el promedio de la suma de los errores al cuadrado respectivas de cada partición, entonces dividiremos entre \(10\times 10\) porque hay \(k=10\) particiones. Luego

(MSEtest1=sum(SE1)/length(Folds))
## [1] 14.18703
(MSEtest2=sum(SE2)/length(Folds))
## [1] 0.8805197

donde en el modelo “sin hacer nada” el error estimado es de 14.26, que es mucho más grande que el obtenido para el error del modelo al hacer la transformación.

Notemos que lo hecho anteriormente se debe hacer \(B\) veces, pues antes sólo realizamos una forma de particionar los datos y en realidad debe probarse efectuando más de una forma de particionar, es decir, lo anterior se debe realizar \(B\) veces y se promedia. Existen dos opciones en R que ya hacen ese proceso.

Comenzamos utilizando la función cv.glm de la librería boot

# Valores necesarios
B=200
K=10

# 10-CV. Validación cruzada
library(boot)
set.seed(123)

# repeticiones
cvLRfit1=rep(NA,B)

for (i in 1:B){
  cvLRfit1[i]=cv.glm(Datos, fitglm, K=K)$delta[1]
}
(mean(cvLRfit1))
## [1] 14.31204

donde cvLRfit1[i]=cv.glm(Datos, fitglm, K=K)$delta[1] realiza una validación cruzada como lo hicimos líneas arriba, y agregamos el bucle for para efectuar las \(B\) repeticiones. Notemos que en el método de validación cruzada aplicado al modelo fitglm obtuvimos el valor de 14.2691, en contraste con el promedio obtenido de la repetición de \(B\) veces en el cual se obtuvo 14.31204. Repetimos para el otro modelo

set.seed(123)
cvLRfit2=rep(NA,B)
for (i in 1:B){
  cvLRfit2[i]=cv.glm(Datos, fit2glm, K=K)$delta[1]
}
(mean(cvLRfit2))
## [1] 0.8571143

La otra alternativa es utilizar la función repCV de la librería cvTools

library(cvTools)
set.seed(123)

# número de particiones y número de repeticiones
folds <- cvFolds(nrow(Datos), K = K, R = B)

# Modelo 1
fit=lm(y~x, data=Datos)
repCV(fit, cost = mspe, folds = folds)
## 10-fold CV results:
##       CV 
## 14.34317
# Modelo 2
fit2=lm(y~Xprima, data=Datos)
repCV(fit2, cost = mspe, folds = folds)
## 10-fold CV results:
##        CV 
## 0.8577174