13.6 Modelos no lineales
En ocasiones un modelo lineal es insuficiente para explicar un fenómeno: la variable objetivo puede no depender linealmente de una variable explicativa. Aunque este tipo de modelos no suelen abordarse en cursos introductorios (ni siquiera de estadística), es ilustrativo mostrar cómo pueden aplicarse de manera sencilla en R.
En el siguiente trozo de código, vamos a cargar un conjunto de datos y ajustador un modelo lineal idéntico al que hubiésemos obtenido con la función lm, aunque usando la función gam del paquete mgcv.
library(mgcv)
pisa <- read.table("data/pisasci2006.csv", header = T, sep = ",")
mod.gam.0 <- gam(Overall ~ Income, data = pisa)
summary(mod.gam.0)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Overall ~ Income
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 204.32 35.37 5.777 4.32e-07 ***
## Income 355.85 46.79 7.606 5.36e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## R-sq.(adj) = 0.518 Deviance explained = 52.7%
## GCV = 1504.5 Scale est. = 1448.8 n = 54
Compara los resultados del modelo anterior con los que se obtendrían con lm.
Lo interesante es que gam permite introducir términos especiales para modelar efectos no lineales. Por ejemplo, en estos datos, los ingresos. El gráfico generado ilustra el impacto (no lineal) del ingreso sobre la variable objetivo.
mod.gam.1 <- gam(Overall ~ s(Income, bs = "cr"), data = pisa)
summary(mod.gam.1)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Overall ~ s(Income, bs = "cr")
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 470.444 4.082 115.3 <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(Income) 6.895 7.741 16.67 1.59e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.7 Deviance explained = 73.9%
## GCV = 1053.7 Scale est. = 899.67 n = 54plot(mod.gam.1)
Usa predict con ambos modelos para estimar su correspondiente error cuadrático medio (la media de la diferencia entre los valores reales y estimados al cuadrado).