Referencias

Las notas se han elaborado a partir de vario material:

A.García Perez, Cap. 11 de Estadística aplicada avanzada con R (EAAR), Madrid, 2022.
Varias aportaciones personales.

Base teórica

Modelos lineales y no lineales

La relación lineal entre una variable dependiente de respuesta \(Y\) y unas variables idependientes o covariables \(\{X_1,X_2,...,X_k\}\), explicada por unos parámetro \(\{\beta_1,...,\beta_k\}\), para valores fijos de \(X_i\) hemos visto que es \[ Y=\beta_0+\beta_1X_1+...+\beta_kX_k \] En ocasiones, la relación entre la variable \(Y\) y las covariables como en \[ E[Y]=\text{exp}\{\beta_0+\beta_1x_1+...+\beta_kx_k\} \] pero puede linealizarse \[ \text{log}E[Y]=\beta_0+\beta_1x_1+...+\beta_kx_k \]

Estas relaciones no lineales que, después de una trasformación, se convierten en Modelos Lineales, son las que dan lugar a los Modelos Lineales Generalizados (Log-lineales, Regresión logística y Regresión de Posson) que ya vimos en otras ocasiones.

Los Modelos de Regresión no Lineal son no lineales en almenos un parámetro.

Carácterísticas de los modelos no lineales

Tenemos normalmente varios parámetros vinculados a cada covariable (en la regresión lineal detemos 1 por cada covariable + 1). Diremos por lo tanto que tenemos unos parámetros \(\theta\) (en otras ocasiones \(\beta\)) vinculados a cada covariable, siendo \(s\) el número de parámetros correspondientes.

En un modelo no lineal, la relación entre la variable dependiente \(Y\) y las independientes es de la forma \[ Y=\eta(x_1,..., x_k;\theta_0,...,\theta_k)+e \] siendo \(e\) una variable de error con distribución \(N(0,\sigma^2)\) y \(\eta\) una función de regresión no lineal en el sentido de que almeno una \(\frac{\partial^2 \eta}{\partial\theta^2}\neq0\).

Con lo que, \(Y\) será una variable aleatoria con media \(\eta(x_1,x_2,...,x_k;\theta_0,\theta_1,...\theta_k)\) que suele simplificarse con \(\eta (x,\theta)\).

Función sigmoide

También conocida como función de Richards o Logística generalizada, es del tipo

\[ \eta(x,\theta)=\theta_1+\frac{\theta_2-\theta_1}{1+\text{exp}\{\theta_3(x-\theta_4)\}}\tag{1.} \]

Función de crecimiento del tipo Weibull

\[ \eta(x,\theta)=\theta_1-\theta_2\;\text{exp}\{-\text{exp}(\theta_3+\theta_4\; \text{log}x)\}\tag{2.} \]

En la última parte del capítulo se analiza asimismo la Regresión Suavizada (smooth), donde no se presupone una forma concreta para la función a ajustar.

Regresión no lineal con R

La obtención de estimaciones de los parámetros del Modelo de Regresión no Lineal se realiza con R mediante la función \(\tt nls()\) (Nonlinear Least Squares) \[ \tt nls(model,data, start,trace) \]

\(\tt model\): es el modelo propuesto.
\(\tt data\): la matriz de datos en formato data frame.
\(\tt start\): un vector numérico que especifique los valores iniciales de las estimaciones de los parámetros.
\(\tt trace\): para indicar que queremos tener también los cálculos iterativos que han llevado a la solución.

Ejemplo 1

Test elisa para detectar anticuerpos del coronavirus de la neumonía asiática en suero de vacas y terneras, obteniéndose las siguientes densidade ópticas (\(\tt opti\)) en 8 disoluciones (\(\tt diso\)):

\[ \begin{array} {l|cccccccc} \text{diso} & 1/30& 1/90& 1/270& 1/810& 1/2430 &1/7290 &1/21869& 1/68609\\ \hline \text{opt}& 1.93 &1.86,1& 84& 1.58& 1.06& 0.56 &0.22 &0.09 \end{array} \]

df <- data.frame(
        diso=c(1/30, 1/90, 1/270, 1/810, 1/2430, 1/7290, 1/21869, 1/68609),
       opt= c(1.93,1.86,1.84,1.58,1.06,0.56,0.22,0.09))

df <- df |> 
  dplyr::mutate(x = -log10(diso)) |>
  dplyr::select(2,3)
df

##    opt        x
## 1 1.93 1.477121
## 2 1.86 1.954243
## 3 1.84 2.431364
## 4 1.58 2.908485
## 5 1.06 3.385606
## 6 0.56 3.862728
## 7 0.22 4.339829
## 8 0.09 4.836381

La variable \(\tt diso\) se convierte en su logaritmo con base 10. Considerando como variable independiente \(x= \log_{10} \tt diso\), la representación gráfica de \(x\) y \(\tt opt\) será

plot(df$x,df$opt, xlab = "x = -log10 diso", ylab = "y = opt",type = "p",pch=16, col= 2,
     main = "Nube de puntos de los datos - Ejemplo 1")

Establecer y ejecutar el modelo de regresión no lineal

Se pretende ajustar a estos datos una función de Regresión no Lineal de la forma \((1.)\), estimando los 4 parámetros \(\theta_1,...,\theta_4\). Mediante la función de R, \(nls()\).

El proceso para utilizar la función \(\tt nls()\) es el siguiente:

Entrar los valores inciales de los parámentros

start <- c(b1=0.04, b2=1.93, b3=2.59, b4=3.47)

start

##   b1   b2   b3   b4 
## 0.04 1.93 2.59 3.47

Establecer la fórmula de Regresión no Lineal que relaciona \(y\) y \(x\). En este caso optamos por una sigmoide:

formula <- opt ~ b1 + (b2 - b1) / (1 + exp(b3 * (x-b4)))

Ejecutar el modelo sobre los datos, con los parámetros elegidos

mod.nl1 <- nls(formula = formula, data = df, start = start)
mod.nl1

## Nonlinear regression model
##   model: opt ~ b1 + (b2 - b1)/(1 + exp(b3 * (x - b4)))
##    data: df
##      b1      b2      b3      b4 
## 0.04255 1.93254 2.58633 3.46541 
##  residual sum-of-squares: 0.003288
## 
## Number of iterations to convergence: 4 
## Achieved convergence tolerance: 6.194e-07

Ajustar el modelo

\[ y=0.043+\frac{1.933}{1+\exp\{2.586(x-3.465)\}} \]

Representar el modelo

Representemos juntos los datos y la función estimada

fnl <- function(x) {
  0.043 + 1.933 / (1 + exp(2.586 * (x - 3.465)))
}
z <- seq(1.5,5, len=100)
plot(df$x, df$opt, xlab = "x = -log10 diso", ylab = "y = opt",type = "p",pch=16, col= 2,
     main = "Nube de puntos de los datos y curva de regresión estimada")
points(z,fnl(z), type = "l", col=4, lwd=1.5)

En \(\tt ggplot\) disponemos de \(\tt stat\_function\) que puede contener la función, sin necesidad de simular valores con \(\tt z\). Aquí además accedemos directamente a los parámetros (\(\tt b1,b2,b3,b4\)) por medio de la función \(\tt coef(model)[index]\)

library(ggplot2)
ggplot(df, aes(x=x, y= opt)) +
  geom_point(colour= "red")+
  stat_function(fun = function(x) coef(mod.nl1)[1] + 
                  (coef(mod.nl1)[2] - coef(mod.nl1)[1]) / (1 + exp(coef(mod.nl1)[3] * (x - coef(mod.nl1)[4]))),
                colour="blue", linewidth = 1.2) +
  ylab("y = opt")+
  xlab("x = -log10 diso")+
  ggtitle("Nube de puntos de los datos y curva de regresión estimada") +
  theme_classic()

Cálculo del modelo utilizando la derivada parcial de la función de regresión no lieal

Derivada de la función con R:

mod1.d <- deriv(y ~ b1 + (b2 - b1) / (1 + exp(b3 * (x-b4))),
      namevec = c("b1","b2","b3","b4"),
      function.arg = function(b1,b2,b3,b4,x){}
)

Calculo del modelo:

nls(opt ~ mod1.d(b1,b2,b3,b4,x), data = df,start = start)

## Nonlinear regression model
##   model: opt ~ mod1.d(b1, b2, b3, b4, x)
##    data: df
##      b1      b2      b3      b4 
## 0.04255 1.93254 2.58633 3.46541 
##  residual sum-of-squares: 0.003288
## 
## Number of iterations to convergence: 4 
## Achieved convergence tolerance: 5.165e-07

Como se ve, se obtiene el mismo resultado

Cálculo de los valores iniciales

R dispone de funciones de autoarranque para poder obtener los valores de arranque a utilizar como \(\tt start\) en la función \(\tt nls()\).

Para la función sigmoide tenemos el modelo \(\tt SSfpl()\) (Self-Starting nls Four-Parameter)¹.

Los parámetros son:

\(\tt input\): un vector numérico de valores para evaluar el modelo.
\(\tt A\): un parámetro numérico que representa la asíntota horizontal del lado izquierdo (valores de entrada muy pequeños).
\(\tt B\): un parámetro numérico que representa la asíntota horizontal del lado derecho (valores de entrada muy grandes).
\(\tt xmid\): un parámetro numérico que representa el valor de entrada en el punto de inflexión de la curva. El valor de \(\tt SSfpl\) estará a medio camino entre \(\tt A\) y \(\tt B\) en \(\tt xmid\).
\(\tt scal\): un parámetro numérico de escala en el eje de entrada.

La relación que une esos parámetros es

\[ y=\tt A+\frac{B-A}{1+\text{exp}\{\tt (xmid-input)/scal\}} \] que, como vemos, es muy similar a la función de la sigmoide que hemos visto arriba. Los ajustes para obtener los valores que necesitamos serán: \[ \begin{array}{ll} \theta_1=\texttt{A}\\ \theta_2=\texttt{B}\\ \theta_3=1/\texttt{xmid}\\ \theta_4=-\texttt{xmid} \end{array} \]

a <- nls(opt ~ SSfpl(-x,A,B,xmid,scal), data=df)

print(c(b1=coef(a)[[1]], b2=coef(a)[[2]], b3=1/coef(a)[[4]], b4=-coef(a)[[3]]), digits = 1)

##   b1   b2   b3   b4 
## 0.04 1.93 2.59 3.47

Que son los valores iniciales que hemos entrado arriba.

Análisis del modelo ajustado

Gráficos de normalidad

par(mfrow=c(1,3))

hist(resid(mod.nl1), main = "Histograma de Residuals",xlab = "Residuals")

plot(resid(mod.nl1), type = "h", main = "Residuals según índice",ylab = "Residuals")
abline(h = 0, col = "red",lwd=2)

qqnorm(resid(mod.nl1),main = "Normalidad de residuals por cuantiles",pch=16, col= 2)

\(\text{Suma de } residuals^2 =\) 0.0033: es una suma muy pequeña. Lo que, junto con los gráficos anteriores, indica un buen ajuste del modelo a los datos.

Test de normalidad de Shapiro-Wilk

Como la muestra es pequeña, utilizamos el test de Shapiro-Wilk que es válido hasta a muestrs de 3, aunque se recomienda más de 20.

shapiro.test(resid(mod.nl1))

## 
##  Shapiro-Wilk normality test
## 
## data:  resid(mod.nl1)
## W = 0.96178, p-value = 0.8269

En cualquier caso, tanto el exámen visual como el test de normalidad muestran un nivel de normalidad más que aceptable.

Comparación con otros modelos posibles: el modelo de regresión lineal simple

También podríamos plantear un modelo de regresión lineal simple, que sería del tipo

\[y=\beta_0+\beta_1x\]

Ejecutando el modelo con la función de regresión lineal \(\tt lm()\) seria

mod.l1 <- lm(opt ~ x, data = df)
mod.l1

## 
## Call:
## lm(formula = opt ~ x, data = df)
## 
## Coefficients:
## (Intercept)            x  
##      3.1348      -0.6326

que generaría una recta que se podría representar en el plano cartesiano junto con la curva sigmoide que hemos dibujado arriba

library(ggplot2)
ggplot(df, aes(x=x, y= opt)) +
  geom_point(colour= "red")+
  stat_function(fun = function(x) coef(mod.nl1)[1] + 
                  (coef(mod.nl1)[2] - coef(mod.nl1)[1]) / (1 + exp(coef(mod.nl1)[3] * (x - coef(mod.nl1)[4]))),
                colour="blue", linewidth = 1.2) +
  stat_function(fun = function(x) coef(mod.l1)[1] + coef(mod.l1)[2]*x,
                colour = "green", linewidth = 1.2, linetype = 1) +
  #geom_smooth(method = "lm",se=F, linetype= 2)+
  ylab("y = opt")+
  xlab("x = -log10 diso")+
  ggtitle("Nube de puntos de los datos y curvas de regresión lineal y de regresión no lineal") +
  theme_classic()

Gráficos de normalidad de la regresión lineal

par(mfrow=c(1,3))

hist(resid(mod.l1), main = "Histograma de Residuals",xlab = "Residuals")

plot(resid(mod.l1), type = "h", main = "Residuals según índice",ylab = "Residuals")
abline(h = 0, col = "red",lwd=2)

qqnorm(resid(mod.l1),main = "Normalidad de residuals por cuantiles",pch=16, col= 2)

Test de Shapiro-Wilks de los residuos de la regresión lineal

shapiro.test(resid(mod.l1))

## 
##  Shapiro-Wilk normality test
## 
## data:  resid(mod.l1)
## W = 0.95739, p-value = 0.7849

La regresión lineal parece también aceptable a nivel de ajuste. Cabe sin embargo destacar que la \(\text{suma de } residuals^2 =\) 0.2657: 8 veces que en la regresión no lineal. Además, en la gráfica de resíduos por índice, los residuos siguen seguramente una pauta de tipo sinusoidal: lo que nos debería orientar hacia una solución diferente de la regresión lineal.

Ejemplo 2

Se observa el crecimiento de la masa verde en tierras de pastura desde el último día de pastoreo. Los datos fruto de la observación aparecen en las tabla

df2 <- data.frame(`dias t`= c(7,15,21,28,32,37,43,52,59),
           `crecimiento y`= c(9.85,11.1,18.84,27.85,39.55,56.32,61.01,64.55,65.10))

dias.t	7.00	15.0	21.00	28.00	32.00	37.00	43.00	52.00	59.0
crecimiento.y	9.85	11.1	18.84	27.85	39.55	56.32	61.01	64.55	65.1

df2 <- df2 |>
  dplyr::rename(t=dias.t, y=crecimiento.y)

plot(df2$t,df2$y, pch=16, col=2, xlab = "t", ylab = "y")

La pauta que se visualiza podría ser la de una función de crecimiento de tipo Weibull, que hemos visto en \((2.)\).

Para inicializar la función \(\tt nls()\), análogamente a como hemos hecho en el Ejemplo 1, hemos de encontrar una función para hacerlo. En R, para la curva de crecimiento de tipo Weibull existe una función específica \(\tt SSweibull()\) de \(\tt stats\).

Según la documentación de \(\tt stats\), los agumentos son:

\(\tt x\): Un vector numérico de valores para evaluar el modelo.
\(\texttt{Asym}=>\theta_1\): Un parámetro numérico que representa la asíntota horizontal del lado derecho (valores muy pequeños de x).
\(\texttt{Drop}=>\theta_2\): Un parámetro numérico que representa el cambio de \(\tt Asym\) a la intersección con el eje \(y\).
\(\texttt{lrc}=>\theta_3\): Un parámetro numérico que representa el logaritmo natural de la constante de velocidad.
\(\texttt{pwr}=>\theta_4\): Un parámetro numérico que representa la potencia a la que se eleva x.

La función sería:

\[ \tt Asym - Drop · exp\{-exp(lrc)·x^{pwr}\} \]

lo que, trasladado, a nuestro entorno, sería:²

\[ \eta(x,\theta)=\theta_1-\theta_2·\;\text{exp}\{-\text{exp}(\theta_3)·x^{\theta_4}\}= \theta_1-\theta_2\;\text{exp}\{-\text{exp}(\theta_3+\theta_4\; \text{log}x)\} \]

Cálculo de los parámetros iniciales

weib <- nls(y ~ SSweibull(t,Asym, Drop, lrc, pwr), data=df2)
start2 <- coef(weib)

start2 <- c(b1= start2[[1]], b2=start2[[2]], b3=start2[[3]], b4=start2[[4]])
start2

##         b1         b2         b3         b4 
##  64.511992  53.848662 -17.008080   4.847144

Cálculo del modelo

formula2 <- y ~ b1 - b2*exp(-exp(b3)*t^b4)

mod.nl2 <- nls(formula = formula2, data = df2, start = start2)

mod.nl2

## Nonlinear regression model
##   model: y ~ b1 - b2 * exp(-exp(b3) * t^b4)
##    data: df2
##      b1      b2      b3      b4 
##  64.512  53.849 -17.008   4.847 
##  residual sum-of-squares: 20.45
## 
## Number of iterations to convergence: 0 
## Achieved convergence tolerance: 5.186e-06

Ajustar el modelo

\[ y=64.51-53.85·\exp\{-\exp(-17)·x^{4.85}) \} \]

Representar el modelo

Representemos juntos los datos y la función estimada

ggplot(df2, aes(x=t, y= y)) +
  geom_point(colour= "red")+
  stat_function(fun = 
                  function(t) coef(mod.nl2)[[1]] - 
                  coef(mod.nl2)[[2]]*exp(-exp(coef(mod.nl2)[[3]])*t^coef(mod.nl2)[[4]]),
                colour="blue", linewidth = 1.2) +
  ylab("y = y")+
  xlab("x = x")+
  ggtitle("Nube de puntos de los datos y curva de regresión estimada") +
  theme_classic()

Análisis del modelo ajustado

Gráficos de normalidad

par(mfrow=c(1,3))
hist(resid(mod.nl2), main = "Histograma de Residuals", xlab = "Residuals")
plot(resid(mod.nl2), type = "h", main = "Residuals según índice", ylab = "Residuals")
abline(h = 0, col = "red",lwd=2)
qqnorm(resid(mod.nl2),main = "Normalidad de residuals por cuantiles",pch=16, col= 2)

Test de Shapiro-Wilks de los residuos de la regresión lineal

shapiro.test(resid(mod.nl2))

## 
##  Shapiro-Wilk normality test
## 
## data:  resid(mod.nl2)
## W = 0.88453, p-value = 0.1751

La regresión lineal parece también aceptable a nivel de ajuste. Cabe sin embargo destacar que la \(\text{suma de } residuals^2 =\) 20.4487: es importante, más que en el ejemplo anterior.

Regresión suavizada

En los modelos no lineales que hemos visto hasta ahora, hemos supuesto la forma de la solución a ajustar, salvo unos parámetros concretos a determinar en el proceso de ajuste. En esta sección vamos a ver como ajustar funciones suaves (smooth), cercanas a los datos, pero sin una forma funcional previamente especificada.

Ejemplo 3

Valores de \(HWC\) en la sangre de niños y niñas menores de 17 años

sangre <- read.table("datosr/sangre.txt", header = T)
head(sangre)

##   Edad HWC
## 1 0.00 230
## 2 0.01 238
## 3 0.01 169
## 4 0.01 186
## 5 0.01 179
## 6 0.01 259

plot(sangre$Edad,sangre$HWC,psc=16,ylab = "HWC",xlab = "Edad",pch=16, col=2,cex=0.75)

Cálculo del modelo: los polinomios

Los polinomios son unas funciones matemáticas que bien se adaptan a las nubes de puntos: a medida que se aumentará el grado, éste se adaptará mejor, y siempre serán modelos lineales.

Podemos por lo tanto imaginar de ajusar un modelo de la siguiente forma: \[ \small HWC=\beta_0+ \beta_1Edad+\beta_2Edad^2+\beta_3Edad^3+\beta_4Edad^4+\beta_5Edad^5+\beta_6Edad^6+\beta_7Edad^7+ \beta_8Edad^8 \]

Ejecución del modelo mediante el operador \(\tt lm()\) de regresión lineal. Notar que la notación \(\tt I()\) en la sintaxis de la fórmula en R significa “tal cual” (“as is”), es decir, \(\tt I(a+b)\) simplemente significa añadir la variable \(\tt a+b\) como predictor en el modelo \(\tt lm\). En su caso, \(\tt I(Viento * Temp^2)\) significa incluir como variable predictora el producto de \(\tt Viento\) y \(\tt Temp\) al cuadrado.

En nuestro caso \(\tt I(Edad^3)\) significa “añadir \(\tt I(Edad^3)\) como coeficiente”, para que no quede añadido al coeficiente \(\tt Edad\). Si conociéramos así, el modelo quedaría como una recta, ya que no entendería que los \(\tt Edad\) que hemos puesto en \(\tt I()\) son otros coeficientes diferentes de \(\tt Edad\).

Conclusión: La función \(\tt I()\) se utiliza para evitar confusiones con los operadores de la sintaxis de la fórmula.

poli <- lm(HWC ~ Edad+I(Edad^2)+I(Edad^3)+I(Edad^4)+I(Edad^5)+I(Edad^6)+I(Edad^7)+I(Edad^8),data = sangre)
poli

## 
## Call:
## lm(formula = HWC ~ Edad + I(Edad^2) + I(Edad^3) + I(Edad^4) + 
##     I(Edad^5) + I(Edad^6) + I(Edad^7) + I(Edad^8), data = sangre)
## 
## Coefficients:
## (Intercept)         Edad    I(Edad^2)    I(Edad^3)    I(Edad^4)    I(Edad^5)  
##   3.033e+02   -2.092e+02    1.134e+02   -3.443e+01    5.898e+00   -5.878e-01  
##   I(Edad^6)    I(Edad^7)    I(Edad^8)  
##   3.384e-02   -1.044e-03    1.336e-05

Bondad del ajuste y recálculo del modelo

Para ver si el ajuste planteado es el adecuado, se realiz un análisis de la varianza con \(\tt{ANOVA}\):

anova(poli)

## Analysis of Variance Table
## 
## Response: HWC
##            Df  Sum Sq Mean Sq  F value    Pr(>F)    
## Edad        1 1268941 1268941 593.4189 < 2.2e-16 ***
## I(Edad^2)   1  366218  366218 171.2613 < 2.2e-16 ***
## I(Edad^3)   1  136239  136239  63.7121 3.353e-14 ***
## I(Edad^4)   1   36707   36707  17.1660 4.494e-05 ***
## I(Edad^5)   1   25445   25445  11.8992 0.0006446 ***
## I(Edad^6)   1   16329   16329   7.6365 0.0060847 ** 
## I(Edad^7)   1    8838    8838   4.1330 0.0429610 *  
## I(Edad^8)   1    1715    1715   0.8019 0.3712694    
## Residuals 291  622262    2138                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

El test recomiendaría descartar el grado 8 y el grado 7 sería aceptable para niveles de significación inferiores a 0.430.

Volvemos por lo tanto a ejecutar la regresión al polinomio con los grado de máxima significación: 6.

\[ \small HWC=\beta_0+ \beta_1Edad+\beta_2Edad^2+\beta_3Edad^3+\beta_4Edad^4+\beta_5Edad^5+\beta_6Edad^6 \]

poli <- lm(HWC ~ Edad+I(Edad^2)+I(Edad^3)+I(Edad^4)+I(Edad^5)+I(Edad^6),data = sangre)

poli

## 
## Call:
## lm(formula = HWC ~ Edad + I(Edad^2) + I(Edad^3) + I(Edad^4) + 
##     I(Edad^5) + I(Edad^6), data = sangre)
## 
## Coefficients:
## (Intercept)         Edad    I(Edad^2)    I(Edad^3)    I(Edad^4)    I(Edad^5)  
##   2.949e+02   -1.563e+02    5.531e+01   -1.025e+01    9.781e-01   -4.551e-02  
##   I(Edad^6)  
##   8.166e-04

Representación gráfica

Con R base

par(mfrow=c(1,2))
z <- seq(0,18, len=500)
plot(sangre$Edad,sangre$HWC,psc=16,ylab = "HWC",xlab = "Edad",pch=16, col=2,cex=0.75,
     main = "Nube de puntos de los datos y curva \nde regresión suavizada")
points(z, predict(poli, data.frame(Edad=z)), type = "l", col="blue", lwd=2.5)

plot(resid(poli), main = "Residuals según índice", ylab = "Residuals",pch=16,cex=0.75)
abline(h = 0, col = "red",lwd=2)

Con \(\tt ggplot2\)

library(patchwork)
a <- ggplot(sangre, aes(x=Edad, y= HWC)) +
  geom_point(colour= "red")+
  stat_function(fun = function(x)
                 coef(poli)[1]+coef(poli)[2]*x+coef(poli)[3]*x^2+coef(poli)[4]*x^3+coef(poli)[5]*x^4+coef(poli)[6]*x^5+
                  coef(poli)[7]*x^6,
                colour="blue", linewidth = 0.8,linetype=1) +
  stat_smooth(method = "lm",
              formula = y ~ x + I(x^2) + I(x^3)+ I(x^4)+ I(x^5)+ I(x^6), # Equivalente a y ~ poly(x, 6), más sencillo, pero menos explícito
              colour= "green", se= F, linetype= 1,linewidth=1)+
  
  # geom_smooth(color= "violet", se=F)+
  #stat_smooth(method = "gam",
  #            formula=y ~ s(x, bs = "cs"), colour= "black")+ # Es el método y la formula que sigue geom_smooth por defecto
  ylab("y = y")+
  xlab("x = x")+
  ggtitle("Nube de puntos de los datos y \ncurva de regresión suavizada") +
  theme_classic()

b <- data.frame(residuals=resid(poli)) |>
  tibble::rownames_to_column(var = 'index') |>
  dplyr::mutate(index= as.numeric(index)) |>
  ggplot(aes(x=index,y=residuals)) +
    geom_point() +
    geom_hline(yintercept = 0, color= "red", lwd=1.2) +
    ggtitle("Residuals según índice") +
    theme_classic()

a + b

Notar la doble posibilidad que tenemos con \(\tt ggplot2\): podemos extraer con \(\tt geom\_stat()\) los coeficientes del modelo que hemos creado previamente o también, con \(\tt stat\_smooth\), podemos pintar la función sin generar previamente el modelo, especificando el tipo de función (polinómica y el número de coeficientes). También existe la función \(\tt geom\_smooth()\) que dibuja un smooth, con el método loess, -ligeramente diferente- por defecto. Pero no podemos controlar en este caso los coeficientes.

Ejemplo 4

Los datos

El set de datos \(\tt Wage\) del paquete \(\tt ISRL\) contiene información sobre 3000 trabajadores. Entre las 12 variables registradas se encuentra el salario (\(\tt wage\)) y la edad (\(\tt age\)). Dada la relación no lineal existente entre estas dos variables, se recurre , se pretende crear un modelo de regression splines que permita predecir el salario de un trabajador en función de la edad que tiene, utilizando una Regresión Syavizada.

library(ISLR)
data("Wage")

plot(Wage$age,Wage$wage, col="grey40", pch=16,cex=0.6, ylab ="Wage",xlab="Age",
     main = "Datos de ingresos salariales segúne edad - Dataset 'Wage'")

Cálculo inicial del modelo

Vamos a plantear inicialmente un modelo con un polínomio de grado 5.

modwage <- lm(wage ~ age+I(age^2)+I(age^3)+I(age^4)+I(age^5),data = Wage)

Bondad del ajuste y recálculo del modelo

Para ver si el ajuste planteado es el adecuado, se realiz un análisis de la varianza con \(\tt{ANOVA}\):

anova(modwage)

## Analysis of Variance Table
## 
## Response: wage
##             Df  Sum Sq Mean Sq  F value    Pr(>F)    
## age          1  199870  199870 125.4443 < 2.2e-16 ***
## I(age^2)     1  228786  228786 143.5931 < 2.2e-16 ***
## I(age^3)     1   15756   15756   9.8888  0.001679 ** 
## I(age^4)     1    6070    6070   3.8098  0.051046 .  
## I(age^5)     1    1283    1283   0.8050  0.369682    
## Residuals 2994 4770322    1593                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

El test recomiendaría descartar el grado 5 y el grado 4 por son superiores a niveles de significación de 0.05.

Volvemos por lo tanto a ejecutar la regresión al polinomio con los grado de máxima significación: 3.

modwage <- lm(wage ~ age+I(age^2)+I(age^3),data = Wage)

modwage

## 
## Call:
## lm(formula = wage ~ age + I(age^2) + I(age^3), data = Wage)
## 
## Coefficients:
## (Intercept)          age     I(age^2)     I(age^3)  
##  -7.524e+01    1.019e+01   -1.680e-01    8.495e-04

Representación gráfica

Con \(\tt ggplot2\)

library(patchwork)
a <- ggplot(Wage, aes(x=age, y= wage)) +
  geom_point(colour= "red", size= 1)+
  stat_function(fun = function(x)
                 coef(modwage)[1]+coef(modwage)[2]*x+coef(modwage)[3]*x^2+coef(modwage)[4]*x^3,
                colour="blue", linewidth = 1,linetype=1) +
  #stat_smooth(method = "lm",
  #            formula = y ~ x + I(x^2) + I(x^3), # Equivalente a y ~ poly(x, 6), más sencillo, pero menos explícito
  #            colour= "green", se= F, linetype= 1,linewidth=1)+
  
  # geom_smooth(color= "violet", se=F)+
  #stat_smooth(method = "gam",
  #            formula=y ~ s(x, bs = "cs"), colour= "black")+ # Es el método y la formula que sigue geom_smooth por defecto
  ylab("Wage")+
  xlab("Age")+
  ggtitle("Nube de puntos de los datos y \ncurva de regresión suavizada") +
  theme_classic()

b <- data.frame(residuals=resid(modwage)) |>
  tibble::rownames_to_column(var = 'index') |>
  #dplyr::mutate(index= as.numeric(index)) |> 
  ggplot(aes(x=index,y=residuals)) +
    geom_point(size=1) +
    geom_hline(yintercept = 0, color= "red", lwd=1.2) +
    ggtitle("Residuals según índice") +
    theme_classic() +
    theme(axis.text.x = element_blank(),
          axis.ticks = element_blank())

a + b

José Pinheiro y Douglas Bates.↩︎
Notar que \[\text{exp}(\theta_3)·x^{\theta_4} =\text{exp}(\theta_3+\theta_4\; \text{log}x)\] La igualdad \(\exp (b+c\cdot \ln (x))=\exp (b)\cdot x^{c}\) se basa en las propiedades de las potencias y logaritmos: la exponencial de una suma es el producto de las exponenciales \(\exp (A+B)=\exp (A)·\exp (B)\) y \(e^{c\cdot \ln (x)}\) se simplifica a \((e^{\ln (x)})^{c}\), lo cual es \(x^{c}\).↩︎

Regresión No Lineal y Regresión Suavizada

Notes d’estadística aplicada - Nota Núm. 11

Stefano Sanfilippo

Març 2026

Referencias

Base teórica

Modelos lineales y no lineales

Carácterísticas de los modelos no lineales

Función sigmoide

Función de crecimiento del tipo Weibull

Regresión no lineal con R

Ejemplo 1

Establecer y ejecutar el modelo de regresión no lineal

Entrar los valores inciales de los parámentros

Establecer la fórmula de Regresión no Lineal que relaciona \(y\) y \(x\). En este caso optamos por una sigmoide:

Ejecutar el modelo sobre los datos, con los parámetros elegidos

Ajustar el modelo

Representar el modelo

Cálculo del modelo utilizando la derivada parcial de la función de regresión no lieal

Cálculo de los valores iniciales

Análisis del modelo ajustado

Gráficos de normalidad

Test de normalidad de Shapiro-Wilk

Comparación con otros modelos posibles: el modelo de regresión lineal simple

Gráficos de normalidad de la regresión lineal

Test de Shapiro-Wilks de los residuos de la regresión lineal

Ejemplo 2

Cálculo de los parámetros iniciales

Cálculo del modelo

Ajustar el modelo

Representar el modelo

Análisis del modelo ajustado

Gráficos de normalidad

Test de Shapiro-Wilks de los residuos de la regresión lineal

Regresión suavizada

Ejemplo 3

Cálculo del modelo: los polinomios

Bondad del ajuste y recálculo del modelo

Representación gráfica

Con R base

Con \(\tt ggplot2\)

Ejemplo 4

Los datos

Cálculo inicial del modelo

Bondad del ajuste y recálculo del modelo

Representación gráfica

Con \(\tt ggplot2\)