2018-09-23

Gráficos de An Introduction to Statistical Learning - Figura 2.6.

Gráfico a replicar

Continuamos con la serie iniciada sobre la creación de gráficos del libro An Introduction to Statistical Learning. En esta ocasión replicaremos los gráficos de la figura 2.6. Utiliza el conjunto de datos Income2 (renta). El gráfico representa los valores observados de renta (en miles de dólares) en función de los años de educación y la antigüedad de 30 individuos. Los puntos rojos representan los valores observados para esas tres variables. La superficie amarilla representa un modelo de suavizado con thin-plate spline. En este caso no hay líneas verticales negras que representaban el error asociado con cada observación, pues la superficie thin-plate spline se adapta perfectamente a los datos.

Solución

Como utilizamos el paquete mgcv que se utiliza para regresión no para interpolación por lo que tenemos que ajustarlo para que la superficie thin plate spline pase por todos los puntos.

  1. Desactivamos la aproximación de bajo rango que utiliza bs = 'tp' estableciendo el parámetro k con el número exacto de datos de la muestra.
  2. xt <- unique(income_2[c("Education", "Seniority")]) 
    nrow(xt)
    
    [1] 30
    
  3. Empleamos sp = 0 para desactivar la penalización de la spline.

  • Primera aproximación
  • model <- gam(Income ~ s(Education, Seniority, k = 30, sp = 0),
                               data = income_2)
    x <- range(income_2$Education)
    x <- seq(x[1], x[2], length.out=30)
    y <- range(income_2$Seniority)
    y <- seq(y[1], y[2], length.out=30)
    z <- outer(x, y, 
               function(Education,Seniority)
                         predict(model, data.frame(Education,Seniority)))
    p <- persp(x, y, z, theta = 30, phi = 30, 
               col = "yellow", expand = 0.5, shade = 0.2, 
               xlab = "Education", ylab = "Seniority", zlab = "Income")
    obs <-  trans3d(income_2$Education, income_2$Seniority, income_2$Income, p)
    pred <-  trans3d(income_2$Education, income_2$Seniority, fitted(model), p)
    points(obs, col = "red", pch = 16)
    segments(obs$x, obs$y, pred$x, pred$y)
    
    Como se puede apreciar en el gráfico, hemos eliminado los errores y se adapta perfectamente a los puntos. Sin embargo, se puede apreciar que la superficie es muy ondulada pues spline es isotrópica o radial. Como dos variables difieren considerablemente en escala, necesitaremos estandarizarlas.

    with(income_2, plot(Education, Seniority, asp = 1))
    
    Las estandarizamos y representamos nuevamente:

    xt_scaled <- scale(xt)
    dat <- data.frame(xt_scaled, Income = income_2$Income)
    with(dat, plot(Education, Seniority, asp = 1))
    
    # Ajustamos el modelo a los datos escalados
    interpolation_model <- gam(Income ~ s(Education, Seniority, k = 30, sp = 0),
                               data = dat)
    # Creamos las coordenadas para el gráfico
    x <- range(dat$Education)
    x <- seq(x[1], x[2], length.out=30)
    y <- range(dat$Seniority)
    y <- seq(y[1], y[2], length.out=30)
    z <- outer(x, y, 
               function(Education,Seniority)
                         predict(interpolation_model, data.frame(Education,Seniority)))
    
    # Volvemos a transformar las coordenadas x e y a su escala original.
    # No es necesario transformar los valores esperados (predicted values, pred)
    scaled_center <- attr(xt_scaled, "scaled:center")
    scaled_scale <- attr(xt_scaled, "scaled:scale")
    xx <- x * scaled_scale[1] + scaled_center[1]
    yy <- y * scaled_scale[2] + scaled_center[2]
    
    # Usamos `xx`, `yy` y `z`
    p <- persp(xx, yy, z, theta = 30, phi = 30,
               col = "yellow",expand = 0.5, shade = 0.2,
               xlab = "Education", ylab = "Seniority", zlab = "Income")
    obs <-  trans3d(income_2$Education, income_2$Seniority, income_2$Income, p)
    pred <-  trans3d(income_2$Education, income_2$Seniority, fitted(interpolation_model), p)
    points(obs, col = "red", pch = 16)
    segments(obs$x, obs$y, pred$x, pred$y)
    

    Entradas relacionadas

    Referencias

    No hay comentarios:

    Publicar un comentario

    Nube de datos