Introducción

En regresión no solo existen modelos de comportamiento lineal, si no que en ocasiones podemos encontrar no lineales en datos o en parámetros, en este caso, trataremos el de comportamiento no lineal de datos.

supongamos el siguiente comportamiento

# Datos
datos <- data.frame(x=c(1,2,3,4,5), y=c(1,1,2,2,4))

# Visualización de datos originales
plot(datos$x, datos$y, 
     main="Datos Simulados",
     xlab="x", ylab="y",
     pch=19, col="blue")

Librerias a usar para análisis

library(tidyverse)
library(lmtest)
library(ppcor)
library(lawstat)

si consideramos la posibilidad de que una curva, ajuste de mejor manera, podemos encontrar, lo siguiente:

\(Y = \beta_0 * e^{\beta_1X} +u_i\) modelo con error aditivo (funcion nls \(ln(Y) \ne ln(\beta_0 * e^{\beta_1X} +u_i)\) )

\(Y = \beta_0 * e^{\beta_1X} *u_i\) modelo con error multiplicativo ( funcion lm \(ln(Y) = ln(\beta_0)+\beta_1X+ ln(u_i)\))

# Modelo linealizado: ln(y) = ln(a) + b*x
mod_exp1 <- lm(log(y)~x, data = datos)

# Obtener parámetros del modelo exponencial
a <- exp(coef(mod_exp1)[1])  # intercepto
b <- coef(mod_exp1)[2]        # pendiente

# Predicciones
x_pred <- seq(min(datos$x), max(datos$x), length.out=100)
y_pred <- a * exp(b * x_pred)

# Gráfico con modelo ajustado
ggplot(datos, aes(x = x, y = y)) +
  geom_point(aes(color = "Datos"), size = 3) +
  geom_line(data = data.frame(x = x_pred, y = y_pred),
            aes(x = x, y = y, color = "Modelo ajustado"),
            linewidth = 1) +
  scale_color_manual(name = "",
                     values = c("Datos" = "blue", 
                               "Modelo ajustado" = "red")) +
  labs(title = "Modelo Exponencial",
       x = "x",
       y = "y") +
  theme_minimal() +
  theme(legend.position = c(0.02, 0.98),
        legend.justification = c(0, 1),
        legend.background = element_rect(fill = "white", color = "black"))

De acuerdo a esto, es posible que los datos presenten este comportamiento. Primero veamos el cumplimiento de supuestos

Supuestos

par(mfrow=c(2,2))
plot(mod_exp1)

par(mfrow=c(1,1))

El más interesante de observar es el gráfico que considera la distancia de COOK. En la figura resalta la posibilidad de una o 2 observaciones de alto impacto en el modelo. Se aconseja revisar y tomar decisión respecto de su consideración.

Por otra parte, el gráfico que presente el comportamiento Residuales vs Ajustados, si bien uno podría observar un patrón, este resalta demasiado al ser un numero reducido de datos.

Respecto del gráfico ajustados vs estandarizados, no se observa un comportamiento de embudo o símil que indique heterocedasticidad.

Supuesto de Normalidad

shapiro.test(residuals(mod_exp1))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(mod_exp1)
## W = 0.68403, p-value = 0.00647

dado que el p value es menor que \(\alpha\), quiere decir que no se puede aceptar la hipótesis de que los residuos de este modelo linealizado presente un comportamiento aproximadamente normal.

NOTA1: cuando ocurre esto debe pensar en:

  1. aumentar el tamaño de la muestra (siempre es lo más aconsejable, aunque muchas veces no practicable)

  2. transformar (no deja de tener sus dolores de cabeza), algunas opciones Box y Cox, WLS (mínimos cuadrados ponderados) entre otros.

  3. usar otros métodos de calculo de la regresión como NLS, o Métodos Robustos.

NOTA2: Impacto de la no normalidad

se estima que

  • Los estimadores de \(\beta_0\) y \(\beta_1\) siguen siendo insesgados (por mínimos cuadrados)
  • Las predicciones puntuales son válidas
  • \(R^2\) podría seguir ajustando bien, pero es difícil observar en pocos datos, revisar overfitting

Supuesto Homocedasticidad

bptest(mod_exp1)
## 
##  studentized Breusch-Pagan test
## 
## data:  mod_exp1
## BP = 2.4639e-29, df = 1, p-value = 1

dado que el p value es mayor que \(\alpha\), quiere decir que no se puede rechazar la hipótesis de que los residuales presenten homogeneidad en la varianza

Supuesto de independencia Temporal (autocorrelacion)

dwtest(mod_exp1)
## 
##  Durbin-Watson test
## 
## data:  mod_exp1
## DW = 3.3333, p-value = 0.9002
## alternative hypothesis: true autocorrelation is greater than 0

dado que el p value es mayor que \(\alpha\), quiere decir que no se puede rechazar la hipótesis de que los residuales presentan una Correlacion entre un valor presente y uno en algún instante anterior es 0, es decir \(R(e_t,e_{t-1})=0\)

Supuesto de independencia de medición (aleatoriedad en la secuencia de mediciones)

runs.test(mod_exp1$residuals)
## 
##  Runs Test - Two sided
## 
## data:  mod_exp1$residuals
## Standardized Runs Statistic = -0.43644, p-value = 0.6625

dado que el p value es mayor que \(\alpha\), no se puede rechazar la hipótesis de que los residuales presenten un comportamiento independiente en la secuencia de medición.

Conclusión análisis de supuestos del modelo 1.

existe una falta de normalidad que complica el continuar con el modelo. Como existe una presunción sobre la validez del modelo, podemos leer la implicancia de la existencia de los parámetros y analizar otro tipo de modelo si se estima conveniente. Entonces revisar los parámetros del modelo.

summary(mod_exp1)
## 
## Call:
## lm(formula = log(y) ~ x, data = datos)
## 
## Residuals:
##       1       2       3       4       5 
##  0.1386 -0.2079  0.1386 -0.2079  0.1386 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -0.48520    0.22989  -2.111   0.1253  
## x            0.34657    0.06931   5.000   0.0154 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2192 on 3 degrees of freedom
## Multiple R-squared:  0.8929, Adjusted R-squared:  0.8571 
## F-statistic:    25 on 1 and 3 DF,  p-value: 0.01539

de acuerdo con lo anterior, la recomendación pasa por plantear un nuevo modelo sin incluir \(\beta_0\)

Segundo modelo

\(Y = e^{\beta_1X} *ln(u_i)\) modelo con error aditivo (funcion nls \(ln(Y) \ne ln(\beta_0 * e^{\beta_1X} +u_i)\))

# Modelo linealizado: ln(y) = b*x+e
mod_exp2 <- lm(log(y)~x-1, data = datos)

# Obtener parámetros del modelo exponencial
b <- coef(mod_exp2)[1]        # pendiente

# Predicciones
x_pred <- seq(min(datos$x), max(datos$x), length.out=100)
y_pred <- exp(b * x_pred)

# Gráfico con modelo ajustado
ggplot(datos, aes(x = x, y = y)) +
  geom_point(aes(color = "Datos"), size = 3) +
  geom_line(data = data.frame(x = x_pred, y = y_pred),
            aes(x = x, y = y, color = "Modelo ajustado"),
            linewidth = 1) +
  scale_color_manual(name = "",
                     values = c("Datos" = "blue", 
                               "Modelo ajustado" = "red")) +
  labs(title = "Modelo Exponencial",
       x = "x",
       y = "y") +
  theme_minimal() +
  theme(legend.position = c(0.02, 0.98),
        legend.justification = c(0, 1),
        legend.background = element_rect(fill = "white", color = "black"))

De acuerdo a esto, es posible que los datos presenten este comportamiento. Primero veamos el cumplimiento de supuestos

Supuestos

par(mfrow=c(2,2))
plot(mod_exp2)

par(mfrow=c(1,1))

Sigue apareciendo en el gráfico que considera la distancia de COOK la posibilidad de una o más observaciones de alto impacto en el modelo. Vuelvo a insistir el impacto que tiene la poca cantidad de datos al generar el modelo. Se aconseja revisar y tomar decisión respecto de su consideración.

Supuesto de Normalidad

shapiro.test(residuals(mod_exp2))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(mod_exp2)
## W = 0.97558, p-value = 0.9097

dado que el p value es mayor que \(\alpha\), quiere decir que no se puede rechazar la hipótesis de que los residuos de este modelo (mod_exp2) linealizado presente un comportamiento aproximadamente normal .

Supuesto Homocedasticidad

bptest(mod_exp1)
## 
##  studentized Breusch-Pagan test
## 
## data:  mod_exp1
## BP = 2.4639e-29, df = 1, p-value = 1

dado que el p value es mayor que \(\alpha\), quiere decir que no se puede rechazar la hipótesis de que los residuales presenten homogeneidad en la varianza

Supuesto de independencia Temporal (autocorrelacion)

dwtest(mod_exp1)
## 
##  Durbin-Watson test
## 
## data:  mod_exp1
## DW = 3.3333, p-value = 0.9002
## alternative hypothesis: true autocorrelation is greater than 0

dado que el p value es mayor que \(\alpha\), quiere decir que no se puede rechazar la hipótesis de que los residuales presentan una Correlacion entre un valor presente y uno en algún instante anterior es 0, es decir \(R(e_t,e_{t-1})=0\)

Supuesto de independencia de medición (aleatoriedad en la secuencia de mediciones)

runs.test(mod_exp1$residuals)
## 
##  Runs Test - Two sided
## 
## data:  mod_exp1$residuals
## Standardized Runs Statistic = -0.43644, p-value = 0.6625

dado que el p value es mayor que \(\alpha\), no se puede rechazar la hipótesis de que los residuales presenten un comportamiento independiente en la secuencia de medición.

Conclusión análisis de supuestos.

Ya no existe una falta de normalidad, esto debido a la obtención de un modelo que responde a la información analizada respecto de los parámetros. Aunque su presunción puede estar sujeta a discución, vale la pena analizar otros modelos, como se planteo anteriormente.

Analisis Ajuste del Modelo

Existe una mayor tranquilidad a la hora de evaluar el modelo. Si bien el R de Pearson es una buena medida, es óptima al cumplir los supuestos sobre el modelo.

Uno de los elementos que se evalúa de forma conjunta es respecto de los parámetros del modelo así como de R

summary(mod_exp2)
## 
## Call:
## lm(formula = log(y) ~ x - 1, data = datos)
## 
## Residuals:
##        1        2        3        4        5 
## -0.21425 -0.42849  0.05041 -0.16383  0.31507 
## 
## Coefficients:
##   Estimate Std. Error t value Pr(>|t|)   
## x  0.21425    0.04035    5.31  0.00605 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2992 on 4 degrees of freedom
## Multiple R-squared:  0.8758, Adjusted R-squared:  0.8447 
## F-statistic:  28.2 on 1 and 4 DF,  p-value: 0.006047

Al respecto se observa la validación de la existencia de parámetros en el modelo, con un p-value de 0.00605 \(\beta_1\) es distinta de 0.

Luego al revisar \(R^2\) y \(R^2_{adj}\) se observa la diferencia entre ellos, esta se puede producir, por un efecto del tamaño de la muestra (pequeño) que adicionalmente a un aumento en la varianza de los residuales (\(\sigma^2_u\) = 0.08954, ver anova)

#analisis existencia de R H0: R=0

cor(predict(mod_exp2),log(datos$y), method = "pearson") # entrega la correlacion general
## [1] 0.9449112
cor.test(~log(y)+x-1,data=datos,method = "pearson")# evalua la existencia de relacion entre Y y X o lo que quiera
## 
##  Pearson's product-moment correlation
## 
## data:  log(y) and x
## t = 5, df = 3, p-value = 0.01539
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.3766144 0.9964629
## sample estimates:
##       cor 
## 0.9449112

Al evaluar la existen de correlación la \(H_0\) que se plantea es si \(\rho=0\), lo cual indica como negativo, al presentar un p value menor \(\alpha\)

#analisis existencia modelo
anova(mod_exp2)
## Analysis of Variance Table
## 
## Response: log(y)
##           Df  Sum Sq Mean Sq F value   Pr(>F)   
## x          1 2.52456 2.52456  28.195 0.006047 **
## Residuals  4 0.35816 0.08954                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Criterio Akaike (AIC)

establece que puedes evaluar como un modelo más adecuado el de menor valor AIC, no obstante es un criterio más

AIC(mod_exp1,mod_exp2)
##          df      AIC
## mod_exp1  3 2.457203
## mod_exp2  2 5.008261

El parámetro de intercepto en el modelo 1 no es significativo (p = 0.1253), por lo que no hay evidencia suficiente para mantenerlo en el modelo.

Entonces, el modelo reducido (sin intercepto) es más parsimonioso y cumple los supuestos de la regresión (normalidad, homocedasticidad de los residuales, independencia, etc.).

Aunque el modelo 1 presenta un AIC algo menor y una varianza del error más baja, la violación de los supuestos hace menos confiable la inferencia basada en él.

Nota: una falta de normalidad que se considera grave y que implique considerar un modelo con más varianza, presenta: - QQ-plot muestra curvatura fuerte (S o U marcada). - Asimetría notable (cola larga hacia un lado). - Kurtosis muy alta (colas pesadas). - Se observan outliers severos. - Cambian las conclusiones del modelo si usas métodos robustos. - La varianza no es constante (heterocedasticidad fuerte).

Esto sí afectaría la inferencia clásica. de lo contrario podria escojer el modelo 1 \(\hat{y} = 0.6156*e^{0.21425*x}\)

Modelo Final

\(\hat{y} = e^{\beta_1*x}\) entonces \(\hat{y} = e^{0.21425*x}\)