Transformaciones y Regresión Polinomial

Metadata y examen inicial de los datos

Vamos a examinar datos de mediciones de longitud de mandíbulas de venado cuya edad se conoce. Aunque estamos interesados en tener una posible curva de calibración, para muestras con edad desconocida, nos interesa tener un modelo de regresión que podamos comparar con modelos de desarrollo de otras especies, o animales creciendo en otros lugares.

Lectura de los datos

Los datos están en un archivo .csv, y contiene datos de edad en años (age.y) y una medida de longitud (mm), en la quijada de venados. [buscar referencias]

venado <- read.csv("jaws.csv")
head(venado)
##       age.y   bone.mm
## 1  5.112000  20.22000
## 2  1.320000  11.11130
## 3 35.240000 140.65000
## 4  1.632931  26.15218
## 5  2.297635  10.00100
## 6  3.322125  55.85634

Visualización de los datos

Mediante una gráfica de puntos (scatter plot) examinaremos si la posible relación entre edad (variable predictora) y longitud (variable respuesta) parece ser lineal.

library(ggplot2)
ggplot(venado, aes(age.y,bone.mm)) + 
  geom_point() + 
  ylab("Largo de quijada, mm") + xlab("Edad, años")

Figura 1. Relación entre edad (años) de venados y una medida de longitud de la quijada (mm).

Modelo de regresión lineal simple

El primer modelo a usar es una línea recta: \[Y = mX +b\] calculando los parámetros mediante el método de LSD.

reg1 <- lm(bone.mm ~ age.y, data = venado)
summary(reg1)
## 
## Call:
## lm(formula = bone.mm ~ age.y, data = venado)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -50.683  -8.625   2.026  16.378  32.238 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  57.1773     5.7363   9.968 1.44e-13 ***
## age.y         1.5264     0.1957   7.799 3.00e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 21.15 on 51 degrees of freedom
## Multiple R-squared:  0.5439, Adjusted R-squared:  0.535 
## F-statistic: 60.82 on 1 and 51 DF,  p-value: 3.003e-10
AIC(reg1)
## [1] 477.8322

Visualización de modelo lineal

ggplot(venado,aes(y=bone.mm,x=age.y)) + 
  geom_point() + 
  geom_smooth(method="lm") +
  ylab("Largo de quijada, mm") + xlab("Edad, años")
## `geom_smooth()` using formula 'y ~ x'

Evaluación del modelo

### evaluación del modelo
#### linearidad
plot(reg1, which = 1)

#### homocedasticidad
plot(reg1, which = 3)

#### normalidad
plot(reg1, which = 2)

Modelo de regresión lineal con transformación

Vamos a tratar de linearizar la relación, mediante una transformación de los datos, siguiendo la regla de Tukey y Mosteller:

Regla de Tukey y Mosteller

Transformación logarítmica de X

library(dplyr)
venado_log <- mutate(venado, age.log = log10(age.y))
head(venado_log)
##       age.y   bone.mm   age.log
## 1  5.112000  20.22000 0.7085908
## 2  1.320000  11.11130 0.1205739
## 3 35.240000 140.65000 1.5470359
## 4  1.632931  26.15218 0.2129679
## 5  2.297635  10.00100 0.3612811
## 6  3.322125  55.85634 0.5214160

Nuevo modelo

reg2 <- lm(bone.mm ~ age.log, data = venado_log)
summary(reg2)
## 
## Call:
## lm(formula = bone.mm ~ age.log, data = venado_log)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -36.858  -9.271   1.008  11.556  26.335 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    8.706      6.790   1.282    0.206    
## age.log       68.265      5.084  13.428   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 14.7 on 51 degrees of freedom
## Multiple R-squared:  0.7795, Adjusted R-squared:  0.7752 
## F-statistic: 180.3 on 1 and 51 DF,  p-value: < 2.2e-16
AIC(reg2)
## [1] 439.3099

Visualización de modelo con transformación

ggplot(venado_log,aes(y=bone.mm,x=age.log)) + 
  geom_point() + 
  geom_smooth(method="lm") +
  ylab("Largo de quijada, mm") + xlab("Edad, log(años)")
## `geom_smooth()` using formula 'y ~ x'

Evaluación del modelo 2

### evaluación del modelo
#### linearidad
plot(reg2, which = 1)

#### homocedasticidad
plot(reg2, which = 3)

#### normalidad
plot(reg2, which = 2)

Modelos polinomiales

Cuando las transformaciones no mejoran mucho los parámetros y evaluación del modelo, podemos probar modelos polinomiales. Estos son modelos lineales que contienen la variable X en un polinomio de potencias de X:

\[Y = \alpha + \beta_1X + \beta_2X^2 + \beta_3X^3 + ... + \beta_nX^n\]

Usualmente se llega hasta la potencia cúbica, para evitar caer en overfitting.

Cálculo de parámetros de los modelos cuadrático y cúbico

#cuadrático
rpol2 <- lm(bone.mm ~ poly(age.y, 2), data = venado)
summary(rpol2)
## 
## Call:
## lm(formula = bone.mm ~ poly(age.y, 2), data = venado)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -30.369  -8.310   0.205   8.819  36.986 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       95.752      2.053  46.635  < 2e-16 ***
## poly(age.y, 2)1  164.925     14.948  11.034 5.30e-15 ***
## poly(age.y, 2)2 -107.877     14.948  -7.217 2.75e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 14.95 on 50 degrees of freedom
## Multiple R-squared:  0.7766, Adjusted R-squared:  0.7677 
## F-statistic: 86.91 on 2 and 50 DF,  p-value: < 2.2e-16
AIC(rpol2)
## [1] 442.0019
# cúbico
rpol3 <- lm(bone.mm ~ poly(age.y, 3), data = venado)
summary(rpol3)
## 
## Call:
## lm(formula = bone.mm ~ poly(age.y, 3), data = venado)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -28.976  -7.941   2.256   9.214  27.087 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       95.752      1.888  50.705  < 2e-16 ***
## poly(age.y, 3)1  164.925     13.748  11.996 3.42e-16 ***
## poly(age.y, 3)2 -107.877     13.748  -7.847 3.28e-10 ***
## poly(age.y, 3)3   43.708     13.748   3.179  0.00256 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 13.75 on 49 degrees of freedom
## Multiple R-squared:  0.8148, Adjusted R-squared:  0.8035 
## F-statistic: 71.86 on 3 and 49 DF,  p-value: < 2.2e-16
AIC(rpol3)
## [1] 434.0621

Análisis de residuales de modelos polinomiales

Cuadrático

### evaluación del modelo
#### linearidad
plot(rpol2, which = 1)

#### homocedasticidad
plot(rpol2, which = 3)

#### normalidad
plot(rpol2, which = 2)

Cúbico

### evaluación del modelo
#### linearidad
plot(rpol3, which = 1)

#### homocedasticidad
plot(rpol3, which = 3)

#### normalidad
plot(rpol3, which = 2)

Gráfica comparativa de los modelos

library(sjPlot)
sjp.poly(venado$bone.mm, venado$age.y, poly.degree = c(1,2,3))
## Polynomial degrees: 1
## ---------------------
## p(x^1): 0.000
## 
## Polynomial degrees: 2
## ---------------------
## p(x^1): 0.000
## p(x^2): 0.000
## 
## Polynomial degrees: 3
## ---------------------
## p(x^1): 0.000
## p(x^2): 0.000
## p(x^3): 0.003

Comparación estadística de los modelos

Para determinar si hay diferencia en los modelos en su ajuste a los datos.

anova(rpol2, rpol3)
## Analysis of Variance Table
## 
## Model 1: bone.mm ~ poly(age.y, 2)
## Model 2: bone.mm ~ poly(age.y, 3)
##   Res.Df     RSS Df Sum of Sq      F   Pr(>F)   
## 1     50 11171.6                                
## 2     49  9261.2  1    1910.4 10.108 0.002559 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Gráfica del modelo escogido

library(ggpmisc)
ggplot(venado, aes(age.y, bone.mm)) +
  geom_point() + 
  geom_smooth(method="lm", se=TRUE, formula = y ~ poly(x,3)) +
  stat_poly_eq(formula = y ~ poly(x,3), 
      aes(label = paste(..eq.label..,..adj.rr.label..,sep="~~~")), 
      parse = TRUE) +
  ylab("Largo de quijada, mm") + xlab("Edad, años")