Modelos No-Lineales

Tomado de:

Si luego de realizar el ajuste de muestros datos al modelo lineal simple, usar transformaciones o usar modelos polinomiales, no quedamos satisfechos (y las estadísticas son insaciables en cuanto a la búsqueda de un mejor modelo), podemos probar con modelos no-lineales.

Selección del modelo

Las opciones para formular un modelo no-lineal son infinitas, así que partir de cero no es una opción práctica; usualmente lo que se hace es consultar algún experto o la literatura, para ver si hay algún modelo ya probado en un sistema similar y a partir de este empezar la búsqueda.

Para el caso del desarrollo de la mandíbula de venado, estudios previos sugieren el siguiente modelo:
\[y = a - b*e^{(-c*x)}\] con los parámetros a, b y c.

Parámetros del modelo

La parte complicada al tratar de encontrar los parámetros finales de un modelo no-lineal, es encontrar unos valores iniciales que permitan una convergencia rápida y posible del algoritmo de computación. El procedimiento usual para establecer los valores iniciales es observar en los datos, como se comporta la ecuación para los valores \(x = 0\) y \(x = \infty\).

Estimación de valores iniciales de los parámetros

A continuación la gráfica de los datos de longitud de mandíbulas de venados y edad (años).

venado <- read.csv("data/jaws.csv")
library(ggplot2)
ggplot(venado, aes(age.y,bone.mm)) + 
  geom_point() + 
  ylab("Largo de quijada, mm") + xlab("Edad, años")

Si colocamos \(x = 0\) en la ecuación, tenemos:
\[y = a - b*e^0\] y por lo tanto \(y = a - b\). Para \(x = \infty\): \[y = a - b*e^{-\infty}\]
y por lo tanto \(y = a\).
Observando la gráfica podemos hacer un estimado de y cuando \(x = \infty\) de 125, y como \(y = a\), entonces \(a = 125\). Observando la gráfica para \(x = 0\), podemos estimar que \(y = 0\), y dado que \(y = a - b = 0\), por lo tanto \(a = b = 125\).
Para estimar c usamos la ecuación con algún punto de los datos y los estimados de a y b:
\[92.8 = 125 - 125*e^{(-c*10.8)}\] rearreglando: \[125*e^{(-c*10.8)} = 32.2\] tomando el logaritmo natural a ambos lados de la ecuación y despejando: \[c = \frac{ln125 - ln32.2}{10.8}\] por lo tanto \(c = 0.126\).

Modelo usando nls

Modelo completo (3 parámetros)

nlm3 <- nls(bone.mm ~ a-b*exp(-c*age.y),
           start=list(a=125,b=125,c=0.126), data = venado)
summary(nlm3)
## 
## Formula: bone.mm ~ a - b * exp(-c * age.y)
## 
## Parameters:
##    Estimate Std. Error t value Pr(>|t|)    
## a 115.05800    2.96285   38.83  < 2e-16 ***
## b 120.60280   10.07002   11.98 2.65e-16 ***
## c   0.12642    0.01978    6.39 5.41e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 13.33 on 50 degrees of freedom
## 
## Number of iterations to convergence: 2 
## Achieved convergence tolerance: 1.128e-06

Modelo con 2 parámetros

Como los parámetros a y b son parecidos (no difieren en más de dos errores estándares). Podemos simplificar a un modelo de dos parámetros:

\[y = a*(1 - e^{(-c*10.8)})\]

nlm2 <- nls(bone.mm ~ a*(1 - exp(-c*age.y)),
           start=list(a=125,c=0.126), data = venado)
summary(nlm2)
## 
## Formula: bone.mm ~ a * (1 - exp(-c * age.y))
## 
## Parameters:
##    Estimate Std. Error t value Pr(>|t|)    
## a 115.58052    2.87139  40.253  < 2e-16 ***
## c   0.11882    0.01245   9.542 6.22e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 13.23 on 51 degrees of freedom
## 
## Number of iterations to convergence: 3 
## Achieved convergence tolerance: 1.319e-06

Comparación de los modelos

Podemos comparar los dos modelos usando el comando anova (se usa para probar diferencias entre muestras)

anova(nlm3, nlm2)
## Analysis of Variance Table
## 
## Model 1: bone.mm ~ a - b * exp(-c * age.y)
## Model 2: bone.mm ~ a * (1 - exp(-c * age.y))
##   Res.Df Res.Sum Sq Df  Sum Sq F value Pr(>F)
## 1     50     8878.3                          
## 2     51     8929.1 -1 -50.839  0.2863  0.595

Al resultar no significativa la diferencia (p = 0.595), podemos quedarnos con el modelo más sencillo de dos parámetros.

Para los modelos no-lineales no tenemos un \(R^2\) para comparar modelos, pero podemos usar el estadístico conocido como S, Residual standard error en la tabla de resultados. Mientras más pequeño el valor de S, mejor el ajuste del modelo. En nuestro caso son muy semejantes, por lo tanto podemos seguir usando el modelo más sencillo.

Comparación con el modelo polinomial seleccionado

# datos del modelo en data.frame
venado$predictions <- predict(nlm2)
ggplot(venado, aes(age.y, bone.mm)) +
  geom_point() + 
  geom_smooth(method="lm", se = F, formula = y ~ poly(x,3)) +
  geom_line(aes(y=predictions), colour="red") +
  ylab("Largo de quijada, mm") + xlab("Edad, años")

¿Ventaja del modelo no-lineal?