---

Modelos Estadísticos. Grado Biotecnología

---

Modelos no parámetricos

Abstract

En este tema se presentan los modelos no paramétricos mediante funciones de suavizado con respuesta Normal.

---

Modelos de suavizado

---

En este tema se presentan los modelos de suavizado lineales. Esto modelos surgen cuando la relación entre la predictora y la respuesta (en el caso de varaibles numéricas) no se puede escribir de forma lineal, sino más bien a través de una función desconocida. En temas anteriores utilizamos los modelos polinómicos para poder capturar comportamientos no lienales entre predictora y respuesta, pero en este caso utilizaremos funciones de suavizado que permiten capturar todo tipo de comportamiento entre ambas.

La mayor dificultad en este tipo de modelos es que no tenemos una forma explícita para la función de suavizado, y por tanto es necesario utilizar las funciones de predcción para obtnener el modelo resultante.

Para una variable predictora y una respuesta el modelo de suavizado viene dado por: \[Y = s(X) + \epsilon\] En situaciones con dos variables predictoras, \(X_1\) y \(X_2\), de tipo numérico se podrían plantear los modelos saturados siguientes:

\[\begin{array}{ll} M0: & Y \sim X_1 + X_2\\ M1: & Y \sim s(X_1) + X_2\\ M2: & Y \sim s(X_1) + s(X_2)\\ \end{array}\]

También resulta posible plantear este tipo de modelos donde se incluyen varaibles predictoras de tipo factor. En este caso debemos plantear un aecuación de suavizado para cada uno de los grupos determinados por el factor.

---

Ejemplos

---

Antes de comenzar necesitamos instalar la librería mgcv que nos permite utilizar las funciones de suavizado que utilizaremos a lo largo de este tema. Veamos los diferentes ejemplos con los que vamos a trabajar. Muchos de ellos ya los hemos utilizado en temas anteriores.

# Cargamos las librerías
library(tidyverse)
library(forcats)
library(broom)
library(reshape2)
library(lmtest)
library(mgcv)
library(MASS)
library(modelr)

airquality. Banco de datos airquality con varaibles Ozone, Solar.R, Wind, Temp, Month, Day. Representamos la relación entre la respuesta (Ozone) y las variables predictoras de tipo numérico (Solar.R, Wind, Temp) mediante modelos lienales y mediante suavizados.

datos <- airquality[,c("Ozone", "Solar.R", "Wind", "Temp")]
datacomp = melt(datos, id.vars='Ozone')
ggplot(datacomp) +
  geom_jitter(aes(value,Ozone, colour=variable),) + 
  geom_smooth(aes(value,Ozone, colour=variable), method=lm, se=FALSE) +
  facet_wrap(~variable, scales="free_x") +
  labs(x = "", y = "Ozono") +
  theme_bw()

# Seleccionamos varaibles de tipo numérico
datos <- airquality[,c("Ozone", "Solar.R", "Wind", "Temp")]
datacomp = melt(datos, id.vars='Ozone')
ggplot(datacomp) +
  geom_jitter(aes(value,Ozone, colour=variable),) + 
  geom_smooth(aes(value,Ozone, colour=variable), method=loess, se=FALSE) +
  facet_wrap(~variable, scales="free_x") +
  labs(x = "", y = "Ozono") +
  theme_bw()

¿Qué podemos decir de las tendencias observadas?

Ejemplo 14. Se desea estudiar la producción de cierto tipo de planta en dos localidades en función de la densidad de plantas. Representamos la información de los modelos saturados asumiendo un modelo lineal (temas anteriores) y mediante un modelo de suavizado

Densidad <- c(23.48, 26.22, 27.79, 32.88, 33.27, 36.79, 37.58, 37.58, 41.49, 
42.66, 44.23, 44.23, 51.67, 55.58, 55.58, 57.93, 58.71, 59.5, 
60.67, 62.63, 67.71, 70.06, 70.45, 73.98, 73.98, 78.67, 95.9, 
96.68, 96.68, 101.38, 103.72, 104.51, 105.68, 108.03, 117.82, 
127.21, 134.26, 137.39, 151.87, 163.61, 166.35, 184.75, 18.78, 
21.25, 23.23, 27.18, 30.15, 31.63, 32.12, 32.62, 32.62, 33.61, 
37.07, 38.55, 39.54, 39.54, 41.02, 42.5, 43.98, 45.47, 49.92, 
50.9, 53.87, 57.82, 61.78, 61.78, 63.75, 67.71, 71.66, 77.59, 
80.56, 86.49, 88.46, 89.45, 90.93, 92.91, 101.81, 103.78, 115.15, 
123.06, 144.31, 155.68, 158.15, 180.39)
Produccion <- c(5.41, 5.46, 5.4, 5.4, 5.29, 5.25, 5.35, 5.25, 5.05, 5.12, 5.29, 
5.04, 5.03, 4.96, 4.84, 5.12, 4.97, 5.02, 4.87, 4.83, 4.74, 4.76, 
4.79, 4.9, 4.74, 4.51, 4.62, 4.58, 4.62, 4.58, 4.47, 4.4, 4.34, 
4.47, 4.44, 4.24, 4.17, 4.2, 4.14, 4.02, 4.14, 4, 5.61, 5.46, 
5.2, 5.18, 4.95, 5.13, 4.93, 5.15, 4.72, 5.05, 4.92, 5.04, 4.82, 
4.99, 4.66, 4.94, 5, 4.7, 4.51, 4.63, 4.68, 4.53, 4.57, 4.55, 
4.6, 4.54, 4.5, 4.24, 4.3, 4.32, 4.29, 4.38, 4.37, 4.26, 4.11, 
4.31, 3.9, 4.04, 3.87, 3.69, 3.66, 3.37) 
Localidad <- c(rep("A",42), rep("B",42))
ejemplo14 <- data.frame(Densidad,Produccion,Localidad)
ggplot(ejemplo14,aes(x = Densidad, y = Produccion, color = Localidad)) +
  geom_point() + 
  geom_smooth(aes(x = Densidad, y = Produccion, color = Localidad), method=lm, se=FALSE) +
  labs(x = "Densidad", y = "Producción") +
  theme_bw()

ejemplo14 <- data.frame(Densidad,Produccion,Localidad)
ggplot(ejemplo14,aes(x = Densidad, y = Produccion, color = Localidad)) +
  geom_point() + 
  geom_smooth(aes(x = Densidad, y = Produccion, color = Localidad), method=loess, se=FALSE) +
  labs(x = "Densidad", y = "Producción") +
  theme_bw()

¿Qué podemos decir de las tendencias observadas?

Ejemplo 13. Se conoce como infiltración el proceso por el cual el agua (riego o lluvia) se va introduciendo bajo la superficie de un terreno cultivado. Este proceso es vital para determinar las cantidades de agua de riego necesarias, para mantener el terreno en condiciones óptimas. Un parámetro habitual que sirve para estudiar dicho proceso es la carga hidráulica. Este depende tanto de la profundidad de la infiltración como del procedimiento de riego usado. Se diseña un experimento para estudiar la carga hidráulica de un terreno bajo diferentes condiciones de riego (denominados tratamientos).

En este caso representamos unicamente el modelo suavizado, ya que el modelo lineal es de tipo polinómico.

tratamiento <- c(rep("A",15),rep("B",15),rep("C",15))
profundidad <-c(10,20,30,40,50,60,70,80,90,100,110,120,130,140,150,10,20,30,40,50,60,70,80,90,100,110,120,130,140,150,10,20,30,40,50,60,70,80,90,100,110,120,130,140,150)
cargahid <- c(-406.90,-345.70,-335.50,-315.10,-304.90,-315.10,-323.26,-335.50,-345.70,-362.02,-374.26,-386.50,-421.18,-435.46,-447.70,-896.50,-737.38,-653.74,-470.14,-406.90,-388.54,-396.70,-396.70,-396.70,-406.90,-419.14,-437.50,-468.10,-466.06,-492.58,-896.50,-855.70,-818.98,-788.38,-678.22,-590.50,-545.62,-515.02,-498.70,-496.66,-517.06,-555.82,-619.06,-623.14,-623.14)
ejemplo13 <- data.frame(tratamiento,profundidad,cargahid)
ggplot(ejemplo13,aes(x = profundidad, y = cargahid, color = tratamiento)) +
  geom_point() + 
  geom_smooth(aes(x = profundidad, y = cargahid, color = tratamiento), method=loess, se=FALSE) +
  labs(x = "Velocidad", y = "Vida") +
  theme_bw()

¿Qué podemos decir de las tendencias observadas?

---

Funciones de suavizado

---

Las funciones de suavizado son los denominados splines que consisten en funciones definidas sobre bases de polinomios. En nuestro caso utilizaremos los denominados splines penalizados o p-splines. Para el ajuste de este tipo de mosdelos utilizaremos la función gam de la libreria mgcv.

La función de suavizado tiene la estructura siguiente:

s(variable, k = , m = , bs ="ps", by = factor)

donde \(k\) es el tamaño de la base de polinomios, \(m\) es el orden de los polinomios, \(bs\) es el tipo de la base de splines utilizados y \(by\) identifica un factor para el ajuste de las curvas de suavizado (efecto de interacción).

Por defecto se toman valores \(k = 10\) cuando nuestra muestra es pequeña y \(k=20\) cuando nuestra muestra es grande. Para \(m\) se toma el valor de 2 y como base de splines se toman los splines penalizados (\(bs = ps\)) que son los que proporcionan mejor ajuste en la mayoría de situaciones prácticas.

---

Bondad del ajuste

---

La bondad del ajuste de este tipo de modelos se basa en los estadísticos AIC y GCV. El segundo de estos es específico de este tipo de modelos ya que se utiliza para estimar la curva de suavizado asociada con la varaible predictora. En ambos casos cuanto menor es el valor mejor será el modelo obtenido.

En esta situaciones se utilizan ambos criterios para comparar el modelo lineal con el suavizado y determinar cual es el mejor de los dos. una vez ajustado el estudio de los residuos propodionan el análisis diagnçostico del modelo en las mismas condiciones que los modelos lineales habituales, y bajo los mismos supuestos de independencia, normalidad y homogeneidad de varianzas.

Vemos al estudio de inferencia, bondad de ajuste y selección del mejor modelo para los diferentes ejemplos considerados.

Ejemplo airquality

Ya pudiomos ver en el análisis gráfico inicial que un modelos lienal sobre cada predictora no parecía muy adecuado y por ese motivo nos planteamos directamente el modelo aditivo. Utilizamos la función summary() para estudiar la capacidad explicativa del modelo y sus diferentes componentes o efectos.

# M1: modelo suavizado completo sin factores
M1 <- gam(Ozone ~ s(Solar.R, k = 20, m = 2, bs = "ps") + 
            s(Wind, k = 20, m = 2, bs = "ps") + 
            s(Temp, k = 20, m = 2, bs = "ps"), data = airquality)
# Resultados del modelo
summary(M1)

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Ozone ~ s(Solar.R, k = 20, m = 2, bs = "ps") + s(Wind, k = 20, 
##     m = 2, bs = "ps") + s(Temp, k = 20, m = 2, bs = "ps")
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   42.099      1.445   29.13   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##               edf Ref.df      F  p-value    
## s(Solar.R)  1.068  1.132 10.118  0.00173 ** 
## s(Wind)     3.588  4.416 20.900 5.49e-14 ***
## s(Temp)    16.392 17.430  7.082 3.05e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.791   Deviance explained = 83.1%
## GCV = 289.38  Scale est. = 231.9     n = 111

Todos los efectos del modelo resultan significativos. Además la capacidad explicativa del modelo (Deviance explained) es muy alta.

Ejemplo 14

El modelo saturado para este ejemplo viene dado por un efecto de la localidad y un efecto de suavizado de la densidad para cada localidad.

# M1: modelo suavizado completo sin factores
M1 <- gam(Produccion ~ Localidad + s(Densidad, k = 10, m = 2, bs = "ps", by = Localidad), data = ejemplo14)
summary(M1)

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Produccion ~ Localidad + s(Densidad, k = 10, m = 2, bs = "ps", 
##     by = Localidad)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.84394    0.01605  301.88   <2e-16 ***
## LocalidadB  -0.32929    0.02272  -14.49   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                          edf Ref.df     F p-value    
## s(Densidad):LocalidadA 2.541  3.101 216.4  <2e-16 ***
## s(Densidad):LocalidadB 4.550  5.333 165.3  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.951   Deviance explained = 95.6%
## GCV = 0.011773  Scale est. = 0.010499  n = 84

Tanto el efecto de localidad como el del suavizado por localidad resultan significativos. Además la capacidad explicativa del modelo (Deviance explained) es muy alta.

Ejemplo 13

El modelo saturado se obtiene asumiendo un efecto de tratamiento y un efecto de suavizado de la profundidad por cada tratamiento.

M1 <- gam(cargahid ~ tratamiento + s(profundidad, k = 10, m = 2, bs = "ps", by = tratamiento), data = ejemplo13)
summary(M1)

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## cargahid ~ tratamiento + s(profundidad, k = 10, m = 2, bs = "ps", 
##     by = tratamiento)
## 
## Parametric coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -363.652      3.238 -112.32   <2e-16 ***
## tratamientoB -131.920      4.579  -28.81   <2e-16 ***
## tratamientoC -277.848      4.579  -60.68   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                               edf Ref.df      F p-value    
## s(profundidad):tratamientoA 4.031  4.797  38.27  <2e-16 ***
## s(profundidad):tratamientoB 7.224  7.928 248.40  <2e-16 ***
## s(profundidad):tratamientoC 6.862  7.584 215.45  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.994   Deviance explained = 99.7%
## GCV = 296.27  Scale est. = 157.24    n = 45

Tanto el efecto de tratamiento como el del suavizado por tratamiento resultan significativos. Además la capacidad explicativa del modelo (Deviance explained) es muy alta.

---

Construcción del modelo

---

En este tipo de modelos tenemos varias situaciones posibles:

• Comparar el modelo lineal con el modelo suavizado.
• Comparar las componentes del modelo aditivo, es decir, comprobar si podemos elimianr algunos de los efectos presentes en el modelo. En este sentido, cuando tenemos más de dos varaibles predictoras se pueden trabajar con todos los modelos que presentamos al inicio de esta unidad.

En todos los casos utilizamos el valor del estadístico GCV para determianr el mejor modelo. Escogeremos aquel que tenga un valor de dicho estadístico más pequeño.

Veamos los diferentes ejemplos, empezando por el más simple.

Ejemplo 13

Comparamos el modelo polinómico con el modelo suavizado. Este ejemplo es el más sencillo ay que involucra unicamente una varaible predictora.

# M1: modelo lineal con interacción
M1 <- gam(cargahid ~ (profundidad+ I(profundidad^2)+I(profundidad^3))* tratamiento, data = ejemplo13)
# M2: modelo de suavizado con interaccion 
M2 <- gam(cargahid ~ tratamiento + s(profundidad, k = 10, m = 2, bs = "ps", by = tratamiento), data = ejemplo13)
data.frame(M1$gcv.ubre,M2$gcv.ubre)

El modelo suavizado es preferible al modelo polinómico.

# Resumen del modelo
summary(M2)

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## cargahid ~ tratamiento + s(profundidad, k = 10, m = 2, bs = "ps", 
##     by = tratamiento)
## 
## Parametric coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -363.652      3.238 -112.32   <2e-16 ***
## tratamientoB -131.920      4.579  -28.81   <2e-16 ***
## tratamientoC -277.848      4.579  -60.68   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                               edf Ref.df      F p-value    
## s(profundidad):tratamientoA 4.031  4.797  38.27  <2e-16 ***
## s(profundidad):tratamientoB 7.224  7.928 248.40  <2e-16 ***
## s(profundidad):tratamientoC 6.862  7.584 215.45  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.994   Deviance explained = 99.7%
## GCV = 296.27  Scale est. = 157.24    n = 45

Pasamos a obtener y representar la predicción del modelo suavizado. Dado que no tenemos una ecuación de predicción propiamente dicha, debemos aproximar los intervalos de predicción de acuerdo al procedimiento básico que establee que un intervalos de confianza al 95% de confianza se obtiene al calcular el intervalo: \[ [estimación - 1.96*sefit, estimación + 1.96*sefit]\]

A continuación vemos los resultados para este ejemplo:

# Creamos grid de predicción
# Utlizamos la opción each para indicar el número de repeticiones (tantas como niveles del factor)
newdata <- data.frame(profundidad = rep(seq(min(ejemplo13$profundidad), max(ejemplo13$profundidad), .01),each=3), tratamiento = factor(c("A", "B","C")))
# Obtenemos la predicción para el modelo ajusatdo
newdata <- data.frame(newdata, predict(M2, newdata, type = "response", se.fit=TRUE))
# Calculamos intervalos de confianza
newdata <- newdata %>% mutate(
  upr = fit + (1.96 * se.fit),
  lwr = fit - (1.96 * se.fit)
)
# Gráfico de la predicción
ggplot(newdata, aes(x = profundidad, y = fit, color = tratamiento)) +
  geom_line() +
  geom_ribbon(aes(ymax = upr, ymin = lwr, fill = tratamiento), alpha = 1/5) +
  labs(x = "Profundidad", y = "Carga Hidráulica", title = "Bandas de predicción") +
  theme_bw()

Pasmos a estudiar el modelo para los datos del ejemplo 14. Plantemaos modelos con ambas predictoras con y sin interacción.

Ejemplo 14

# M1: modelo lineal sin interacción
M1 <- gam(Produccion ~ Densidad + Localidad, data = ejemplo14)
# M2: modelo lineal con interacción
M2 <- gam(Produccion ~ Densidad + Localidad + Localidad*Densidad, data = ejemplo14)
# M3: modelo de suavizado sin interaccion 
M3 <- gam(Produccion ~ Localidad + s(Densidad, k = 10, m = 2, bs = "ps"), data = ejemplo14)
# M4: modelo de suavizado con interaccion 
M4 <- gam(Produccion ~ Localidad + s(Densidad, k = 10, m = 2, bs = "ps", by = Localidad), data = ejemplo14)
data.frame(M1$gcv.ubre,M2$gcv.ubre,M3$gcv.ubre,M4$gcv.ubre)

El mejor modelo es el M4 correspondiente al modelo suavizado con efecto de interacción entre localidad y densidad.

# Resumen del modelo
summary(M4)

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Produccion ~ Localidad + s(Densidad, k = 10, m = 2, bs = "ps", 
##     by = Localidad)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.84394    0.01605  301.88   <2e-16 ***
## LocalidadB  -0.32929    0.02272  -14.49   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                          edf Ref.df     F p-value    
## s(Densidad):LocalidadA 2.541  3.101 216.4  <2e-16 ***
## s(Densidad):LocalidadB 4.550  5.333 165.3  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.951   Deviance explained = 95.6%
## GCV = 0.011773  Scale est. = 0.010499  n = 84

Estudiamos ahora la predicción para dicho modelo.

# Creamos grid de predicción
# Utlizamos la opción each para indicar el número de repeticiones (tantas como niveles del factor)
newdata <- data.frame(Densidad = rep(seq(min(ejemplo14$Densidad), max(ejemplo14$Densidad), .01),each=2), Localidad = factor(c("A", "B")))
# Obtenemos la predicción para el modelo ajusatdo
newdata <- data.frame(newdata, predict(M4, newdata, type = "response", se.fit=TRUE))
# Calculamos intervalos de confianza
newdata <- newdata %>% mutate(
  upr = fit + (1.96 * se.fit),
  lwr = fit - (1.96 * se.fit)
)
# Gráfico de la predicción
ggplot(newdata, aes(x = Densidad, y = fit, color = Localidad)) +
  geom_line() +
  geom_ribbon(aes(ymax = upr, ymin = lwr, fill = Localidad), alpha = 1/5) +
  labs(x = "Densidad de plantas", y = "Producción", title = "Bandas de predicción") +
  theme_bw()

Veamos el último ejemplo.

Ejemplo airquality

# M1: modelo sin suavizado
M1 <- gam(Ozone ~ Solar.R + Wind + Temp, data = airquality)
# M2: modelo con suavizado para Solar.R
M2 <- gam(Ozone ~ s(Solar.R, k = 20, m = 2, bs = "ps") + Wind + Temp, data = airquality)
# M3: modelo con suavizado para Wind
M3 <- gam(Ozone ~ Solar.R + s(Wind, k = 20, m = 2, bs = "ps") + Temp, data = airquality)
# M4: modelo con suavizado para Temp
M4 <- gam(Ozone ~ Solar.R + Wind + s(Temp, k = 20, m = 2, bs = "ps"), data = airquality)
# M5: modelo suavizado en Solar.R y Wind
M5 <- gam(Ozone ~ s(Solar.R, k = 20, m = 2, bs = "ps") + 
            s(Wind, k = 20, m = 2, bs = "ps") + Temp, data = airquality)
# M6 suavizado en Solar.R y Temp
M6 <- gam(Ozone ~ s(Solar.R, k = 20, m = 2, bs = "ps") + 
            Wind + 
            s(Temp, k = 20, m = 2, bs = "ps"), data = airquality)
# M6 suavizado en Wind y Temp
M7 <- gam(Ozone ~ Solar.R + 
            s(Wind, k = 20, m = 2, bs = "ps") + 
            s(Temp, k = 20, m = 2, bs = "ps"), data = airquality)
# M8: modelo suavizado completo
M8 <- gam(Ozone ~ s(Solar.R, k = 20, m = 2, bs = "ps") + 
            s(Wind, k = 20, m = 2, bs = "ps") + 
            s(Temp, k = 20, m = 2, bs = "ps"), data = airquality)
data.frame(M1$gcv.ubre,M2$gcv.ubre,M3$gcv.ubre,M4$gcv.ubre,M5$gcv.ubre,M6$gcv.ubre,M7$gcv.ubre,M8$gcv.ubre)

El mejor modelos corresponde al M8. Sin embargo, en situaciones donde se involucran diferentes predictoras de tipo numérico resulta posible ajustar un suavizado conjunto para ambas variables. Para ello utilizamos la función te(). Comparamos los resultados el modeloa anterior con los suavizados dobles:

# M1: Doble en Wind y Temp
M1 <- gam(Ozone ~ s(Solar.R, k = 20, m = 2, bs = "ps") + te(Wind,Temp), data = airquality)
# M2: Doble en Solar.R y Wind 
M2 <- gam(Ozone ~ te(Solar.R, Wind) + s(Temp, k = 20, m = 2, bs = "ps"), data = airquality)
# M3: Doble en Solar.R y Temp
M3 <- gam(Ozone ~ te(Solar.R, Temp) + s(Wind, k = 20, m = 2, bs = "ps"), data = airquality)
# M4: modelo suavizadoindividual
M4 <- gam(Ozone ~ s(Solar.R, k = 20, m = 2, bs = "ps") + 
            s(Wind, k = 20, m = 2, bs = "ps") + 
            s(Temp, k = 20, m = 2, bs = "ps"), data = airquality)
# Resultados
data.frame(M1$gcv.ubre,M2$gcv.ubre,M3$gcv.ubre,M4$gcv.ubre)

El mejor modelo es el M1 correspondiente a un suavizado individual para Solar.R y un suavizado conjunto para Wind y Temp. Veamos el resumne del modelo

# Resumen del modelo
summary(M1)

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Ozone ~ s(Solar.R, k = 20, m = 2, bs = "ps") + te(Wind, Temp)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   42.099      1.411   29.85   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                  edf Ref.df      F p-value    
## s(Solar.R)     2.166  2.719  3.581  0.0235 *  
## te(Wind,Temp) 13.140 14.816 22.721  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.801   Deviance explained = 82.8%
## GCV = 258.86  Scale est. = 220.84    n = 111

Estudiamos ahora la predicción para dicho modelo.

# Creamos grid de predicción
g1 <- seq(min(airquality$Solar.R,na.rm = TRUE), max(airquality$Solar.R,na.rm = TRUE),length = 5)
g2 <- seq(min(airquality$Wind,na.rm = TRUE), max(airquality$Wind,na.rm = TRUE),length = 5)
g3 <- seq(min(airquality$Temp,na.rm = TRUE), max(airquality$Temp,na.rm = TRUE),length = 5)
newdata <- expand.grid(Solar.R = g1, Wind = g2, Temp = g3)
# Obtenemos la predicción para el modelo ajusatdo
newdata <- data.frame(newdata,predict(M1, newdata, type = "response", se.fit=TRUE))
# Calculamos intervalos de confianza
newdata <- newdata %>% mutate(
  upr = fit + (1.96 * se.fit),
  lwr = fit - (1.96 * se.fit)
)
# Gráfico para este tipo de modelos
par(mfrow = c(1,2)) # el 2 indica el número de efectos del modelo
plot(M1, shade = TRUE, seWithMean = TRUE, scale = 0)

# Para visualizar mejor la relación en el suavizado doble instalamos la libreria visreg
library(visreg)
visreg2d(M1, xvar='Wind', yvar='Temp', scale='response')

---

Diagnóstico del modelo

---

Las hipótesis de este modelo son las mismas que las del modelo de regresión, y por tanto, el diagnóstico se centra en los mismos procedimientos gráficos. Recordemos que debemos verificar la normalidad y varianza constante de los residuos del modelo. Sin embargo, en este tipo de modelos se añáde el diagnóstico sobre el grado de suavizado del modelo considerado para saber si es necesario modificarlo.

Para este tipo de modelos resulta posible obtener todos los gráficos de interés con la función gam.check(). Veamos el análisis de cada uno de los ejemplos.

Ejemplo 13

Utilizamos el mejor modleo obtenido en el paso anterior

# M2: modelo de suavizado con interaccion 
M2 <- gam(cargahid ~ tratamiento + s(profundidad, k = 10, m = 2, bs = "ps", by = tratamiento), data = ejemplo13)
# Gráficos de diagnóstico
gam.check(M2)

## 
## Method: GCV   Optimizer: magic
## Smoothing parameter selection converged after 7 iterations.
## The RMS GCV score gradient at convergence was 1.786611e-05 .
## The Hessian was positive definite.
## Model rank =  30 / 30 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##                               k'  edf k-index p-value
## s(profundidad):tratamientoA 9.00 4.03     1.1    0.73
## s(profundidad):tratamientoB 9.00 7.22     1.1    0.68
## s(profundidad):tratamientoC 9.00 6.86     1.1    0.68

Todos los gráficos de diagnóstico parece indicar que los residuos cumplen con als hipótesis del modelo (qq-plot, Resids vs linear pred., histogram of residuals, Response vs. ftted values). Por otro lado el análisis del parámetro de suavizado indica que el ajuste obtenido es adecuado, no encontrando problemas importantes.

Ejemplo 14

# Modelo
M4 <- gam(Produccion ~ Localidad + s(Densidad, k = 10, m = 2, bs = "ps", by = Localidad), data = ejemplo14)
# Diagnóstico
gam.check(M4)

## 
## Method: GCV   Optimizer: magic
## Smoothing parameter selection converged after 5 iterations.
## The RMS GCV score gradient at convergence was 3.261881e-07 .
## The Hessian was positive definite.
## Model rank =  20 / 20 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##                          k'  edf k-index p-value
## s(Densidad):LocalidadA 9.00 2.54    1.08    0.76
## s(Densidad):LocalidadB 9.00 4.55    1.08    0.67

No se detecta ningún problema en la fase de diagnóstico para este modelo.

airquality

# Modelo
M1 <- gam(Ozone ~ s(Solar.R, k = 20, m = 2, bs = "ps") + te(Wind,Temp), data = airquality)
# Diagnóstico
gam.check(M1)

## 
## Method: GCV   Optimizer: magic
## Smoothing parameter selection converged after 11 iterations.
## The RMS GCV score gradient at convergence was 3.152251e-05 .
## The Hessian was positive definite.
## Model rank =  44 / 44 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##                  k'   edf k-index p-value  
## s(Solar.R)    19.00  2.17    0.87    0.05 *
## te(Wind,Temp) 24.00 13.14    1.07    0.76  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

A pesar de que el modelo parece cumplir todas las hipótesis, el test sobre el suavizado indica que el modelo no es del todo correcto. Tomamos el segundo mejor modelo (el que tenía un suavizado individual para cada variable) y comprobamos de nuevo las hipótesis.

# Modelo
M1 <- gam(Ozone ~ s(Solar.R, k = 20, m = 2, bs = "ps") + 
            s(Wind, k = 20, m = 2, bs = "ps") +
            s(Temp, k = 20, m = 2, bs = "ps"), data = airquality)
# Diagnóstico
gam.check(M1)

## 
## Method: GCV   Optimizer: magic
## Smoothing parameter selection converged after 17 iterations.
## The RMS GCV score gradient at convergence was 5.831016e-05 .
## The Hessian was positive definite.
## Model rank =  58 / 58 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##               k'   edf k-index p-value
## s(Solar.R) 19.00  1.07    1.06    0.73
## s(Wind)    19.00  3.59    1.22    0.99
## s(Temp)    19.00 16.39    1.09    0.81

Podemos ver ahora (gráfico de Resids vs linear pred) que se incumple la hipótesis de varaina comstante ya que la variabilidad de los residuos aumenta cuando lo hace el predictor lineal. Este tipo de problemas se suele arraglar con la tranasformación raíz cuadrada. Aplicamos dicha transformación y ajustamos el nuevo modelo.

# Modelo
M1 <- gam(sqrt(Ozone) ~ s(Solar.R, k = 20, m = 2, bs = "ps") + 
            s(Wind, k = 20, m = 2, bs = "ps") +
            s(Temp, k = 20, m = 2, bs = "ps"), data = airquality)
# Bondad del ajuste
summary(M1)

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## sqrt(Ozone) ~ s(Solar.R, k = 20, m = 2, bs = "ps") + s(Wind, 
##     k = 20, m = 2, bs = "ps") + s(Temp, k = 20, m = 2, bs = "ps")
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   6.0165     0.1035   58.14   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##               edf Ref.df      F  p-value    
## s(Solar.R)  1.765  2.211  8.875 0.000214 ***
## s(Wind)     3.102  3.864 16.720 1.82e-10 ***
## s(Temp)    16.244 17.334  7.532 3.28e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =    0.8   Deviance explained = 83.9%
## GCV = 1.4844  Scale est. = 1.1887    n = 111

# Diagnóstico
gam.check(M1)

## 
## Method: GCV   Optimizer: magic
## Smoothing parameter selection converged after 12 iterations.
## The RMS GCV score gradient at convergence was 1.619835e-07 .
## The Hessian was positive definite.
## Model rank =  58 / 58 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##               k'   edf k-index p-value
## s(Solar.R) 19.00  1.76    1.11    0.89
## s(Wind)    19.00  3.10    1.11    0.89
## s(Temp)    19.00 16.24    1.12    0.89

Con este nuevo modelo se verifican todas las hipótesis y lo tomamos como modelo resultante. La predcción correspondiente viene dada por:

# Gráfico para este tipo de modelos
par(mfrow = c(1,3)) # el 3 indica el número de efectos del modelo
plot(M1, shade = TRUE, seWithMean = TRUE, scale = 0)

Aunque el modelo cumple todas las hipótesis el suavizado para la varaible temperatura es muy irregular y nos plantemos reducir el valor de \(k\) para conseguir un suavizado más adecuado. Probamos reduciendo en dos unidades. Tnede en cuenta que al tratarse un proceso iterativo de estimación los pvalores obtenidos son aproximados:

# Modelo
M1 <- gam(sqrt(Ozone) ~ s(Solar.R, k = 20, m = 2, bs = "ps") + 
            s(Wind, k = 20, m = 2, bs = "ps") +
            s(Temp, k = 18, m = 2, bs = "ps"), data = airquality)
# Bondad del ajuste
summary(M1)

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## sqrt(Ozone) ~ s(Solar.R, k = 20, m = 2, bs = "ps") + s(Wind, 
##     k = 20, m = 2, bs = "ps") + s(Temp, k = 18, m = 2, bs = "ps")
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   6.0165     0.1164   51.69   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##              edf Ref.df      F  p-value    
## s(Solar.R) 2.183  2.740  7.355 0.000309 ***
## s(Wind)    2.563  3.233 12.178 3.52e-07 ***
## s(Temp)    4.094  5.051 14.615 1.77e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.747   Deviance explained = 76.8%
## GCV = 1.6504  Scale est. = 1.504     n = 111

# Diagnóstico
gam.check(M1)

## 
## Method: GCV   Optimizer: magic
## Smoothing parameter selection converged after 11 iterations.
## The RMS GCV score gradient at convergence was 6.773148e-07 .
## The Hessian was positive definite.
## Model rank =  56 / 56 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##               k'   edf k-index p-value  
## s(Solar.R) 19.00  2.18    0.97   0.310  
## s(Wind)    19.00  2.56    0.98   0.405  
## s(Temp)    17.00  4.09    0.84   0.065 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Gráfico para este tipo de modelos
par(mfrow = c(1,3)) # el 3 indica el número de efectos del modelo
plot(M1, shade = TRUE, seWithMean = TRUE, scale = 0)

---

Copyright © 2018 Javier Morales. Universidad Miguel Hernández de Elche.