Introducción a Regresión no lineal con R

---

Existen ciertos fenómenos que se encuentran mejor modelados mediante funciones no lineales que con funciones lineales. En este recurso, presentamos una introducción sobre la aplicación de regresión no lineal en R, una técnica para ajustar modelos no lineales a un conjunto de datos obtenidos en contextos experimentales u observacionales.

---

Librerías

Usaremos las siguientes librerías

library(tidyverse)

Literatura recomendada

• Ritz, C., & Streibig, J. C. (2008). Nonlinear regression with R. Springer Science & Business Media

• Ritz, C., Jensen, S. M., Gerhard, D., & Streibig, J. C. (2020). Dose-Response Analysis using R. CRC Press.

• Archontoulis, S. V., & Miguez, F. E. (2015). Nonlinear regression models and applications in agricultural research. Agronomy Journal, 107(2), 786–798

Introducción

Cuando la variable respuesta tiene una relación no lineal con ciertas predictoras, ajustar un modelo no lineal puede ser una alternativa más parsimoniosa que otras opciones (Archontoulis & Miguez, 2015) como regresión lineal con términos polinómicos o ajustar curvas loess. La Figura 1 presenta algunos ejemplos de datos reales que exhiben patrones no lineales.

Figure 1: Ejemplos de datos con patrones no lineales. En los tres casos, los puntos representan los datos observados y las líneas indican el compartamiento medio estimado por una función o modelo no lineal. (A) Experimento de dosis-respuesta en el cual la biomasa de plantas de lechuga cambia con un aumento en la concentración de alcohol isobutílico. (B) Tasa de crecimiento relativo en el tiempo para una especie vegetal. (C) Tasas de reacción enzimática en función de la cantidad de sustrato en cultivos celulares tratados y no tratados con Puromicina (En R, estos datos se activan como: data(Puromycin)).

En la Figura 1, las curvas en cada ejemplo representan el compartamiento medio estimado por un modelo o función no lineal. Por ejemplo, de acuerdo a Ritz & Streibig (2008), para los datos lettuce.csv de la Figura 1A, se ajustó el modelo no lineal de Brain Cousens dado por la ecuación:

\[\begin{equation} f(x) = \dfrac{\lambda + \theta\cdot x}{1+\left(\frac{x}{\tau}\right)^{\beta}} \tag{1} \end{equation}\]

El modelo (1) tiene cuatro parámetros (\(\lambda\), \(\theta\), \(\tau\) y \(\beta\)) y particularmente puede ser útil para describir relaciones de dosis - respuesta donde se espera que la respuesta exhiba un deacamiento cuando hay exposición a una sustancia potencialmente tóxica, pero cuando la caída se evidencia en altos valores de la sustancia (Ritz & Streibig, 2008). En otro ejemplo, el modelo ajustado para los datos de Puromicina de la Figura 1C es el de Michaelis-Menten, dado por:

\[\begin{equation} f(x) = \dfrac{(V_m + \beta \cdot z) \cdot x}{K + x} \tag{2} \end{equation}\]

En el modelo (2), \(x\) es la concentración del sustrato (una variable predictora), \(z\) es 1 si el cultivo es tratado con Puromicina o 0 si no lo es (una segunda variable predictora), mientras que \(V_m\), \(\beta\) y \(K\) son parámetros del modelo.

En este recurso revisaremos, de forma introductoria, como ejecutar la técnica de regresión no lineal en R, para estimar los parámetros de modelos no lineales como la función (1) o la función (2) desde datos o muestras obtenidas en escenarios como los ejemplos de la Figura 1. La teoría y ejemplos mostrados son basados principalmente en el libro de Ritz & Streibig (2008).

Conceptos sobre regresión no lineal

El modelo

Para un conjunto de pares de datos (x_i, y_i) con \(i = \{1,2,\ldots,n\}\) y donde \(x\) es la variable predictora y \(y\) es la variable respuesta, la técnica de regresión no lineal considera la siguiente relación entre la respuesta \(y\) y la predictora \(x\):

\[\begin{equation} y_i = f(x_i,\beta) + \varepsilon_i \tag{3} \end{equation}\]

De acuerdo a Ritz & Streibig (2008), una descripción del modelo (3) es la siguiente:

• \(f(x_i,\beta)\) es alguna función matemática que representa el comportamiento medio de \(y\) cuando la predictora \(x\) incrementa. Se aclara que la función \(f\) puede tener más de una predictora como el ejemplo de la Figura 1C.

• La función \(f\) también depende de \(\beta\) que representa los parámetros del modelo. Estos deben ser estimados desde el conjunto de datos. Si existen \(n\) datos y el modelo tiene \(p\) parámetros, debe cumplirse que \(n > p\) para poder realizar la estimación.

• Otro aspecto de la función \(f\) es que debe ser no lineal en al menos uno de los parámetros. El concepto de no lineal en los parámetros quiere decir, que la función \(f\) no se puede escribir como una combinación lineal de las variables predictoras y los parámetros. Si la función \(f\) fuera lineal en sus parámetros, podríamos/deberíamos usar la regresión lineal clásica.

• \(\varepsilon_i\) representa el error o diferencia que existe entre cada \(y_i\) y la media de \(y_i\). Dicho error es aleatorio y ocurre debido a diferentes fuentes de variación no contraladas en el proceso de medición. La técnica de regresión no lineal asume que los errores siguen una distribución normal con media 0 y cierta varianza \(\sigma^2\).

Mínimos cuadrados

La técnica de regresión no lineal estima los parámetros del modelo usando el criterio de mínimos cuadrados, en el cual se encuentran los valores para los parámetros (\(\beta\) en el modelo (3)) de tal forma que se minimice la siguiente suma:

\[\begin{equation} RSS(\beta) = \sum_{i = 1}^n (\ y_i - f(x_i,\beta)\ )^2 \tag{4} \end{equation}\]

La minimización de RSS es un problema no lineal debido a la no linealidad de la función \(f\) y por ende se requieren métodos de optimización numérica para su ejecución (Ritz & Streibig, 2008). Existen varios métodos para optimización numérica, siendo quizas el de Gauss-Newton el más común (Ritz & Streibig, 2008).

Valores iniciales

Un aspecto relevante para que los métodos de optimización se ejecuten con éxito es que se deben proveer valores iniciales para los parámetros (\(\beta\)) del modelo. Para esto tenemos dos opciones:

• Calcular (a mano) un estimado grueso de los parámetros que sirva como valores iniciales e indicarle directamente dichos valores al programa que realiza la estimación.

• Usar funciones de auto inicio (Self Start) que internamente calculan dichos valores iniciales y se los entregan al programa que realiza la estimación. Para cada modelo \(f\) que se quiera ajustar se requiere una función de auto inicio específica.

Comando nls de R

El comando nls realiza regresión no lineal. Sus argumentos principales son:

• formula: Aquí se específica el modelo. Existen al menos tres maneras de especificarla.
• data: El data.frame que contiene los datos para la respuesta y las predictoras. Los nombres de las columnas en este data.frame deben ser usados en la fórmula.
• start: Vector (o lista) nombrado con los valores iniciales de los parámetros.

Ejemplo 1: RGRcurve

Considere los datos sobre el crecimiento relativo de cierta especie vegetal en el tiempo de la Figura 1B:

# Codigo R para importar los datos:
RGRcurve <- read.csv("RGRcurve.csv")

'data.frame':   41 obs. of  2 variables:
 $ Day: int  0 0 0 0 0 0 1 1 1 1 ...
 $ RGR: num  0.169 0.164 0.21 0.215 0.183 0.181 0.188 0.229 0.268 0.202 ...

Use estos datos para ajustar el modelo exponencial de dos parámetros dado por:

\[\begin{equation} f(x) = y_0 e^{x/\theta} \tag{5} \end{equation}\]

donde \(x\) es la predictora (Day), y \(f\) es la respuesta (RGR). Los parámetros son \(y_0\) y \(\theta\). Observe que si \(x = 0\), \(f(0) = y_0\), de modo que \(y_0\) es el valor de la RGR en el Day = 0. De otro lado, \(\theta\) debe ser un tipo de pendiente. Observe que cuando \(x = \theta\), la RGR ha duplicado (aproximadamente) el valor \(y_0\).

Ejercicio 1

Escriba código R para realizar un diagrama de dispersión con los datos y, desde el gráfico, genere un estimado grueso (al ojo) del valor de \(y_0\) y \(\theta\) usando la descripción de estos parámetros dada anteriormente.

Código solución del gráfico

# Lectura de datos
RGRcurve <- read.csv("RGRcurve.csv")

# Grafico
ggplot(RGRcurve, aes(x = Day, y = RGR)) + geom_point() +
  scale_y_continuous(breaks = seq(0,2,0.1))

Razonamiento sobre los valores para y₀ y θ

El valor y0 debe estar alrededor de 0.2. Este valor se duplica (2*0.2 = 0.4) cuando
Day esta entre 3 y 4. Luego el valor de theta debe estar cerca a 3.5.

Ejercicio 2

Escriba código R que realice una diagrama de dispersión de los datos y agregue la curva del modelo utilizando los valores para los parámetros encontrados en el ejercicio anterior.

Código solución del gráfico

# Grafico base
gg <- ggplot(RGRcurve, aes(x = Day, y = RGR)) + geom_point() +
  scale_y_continuous(breaks = seq(0,2,0.1))

# Se agregan curvas para algunos valores de theta:
gg +
  geom_function(aes(color = "theta = 3.5"), 
                fun = function(x) 0.2*exp(x/3.5) ) +
  geom_function(aes(color = "theta = 4"), 
                fun = function(x) 0.2*exp(x/4) ) +
  geom_function(aes(color = "theta = 5"), 
                fun = function(x) 0.2*exp(x/5) )

El siguiente código R ejecuta el comando nls utilizando como valores iniciales para los parámetros, los estimados generados “al ojo” desde los ejercicios 1 y 2.

# Codigo R para ejecutar regresion no lineal con el comando nls:
m <- nls(RGR ~ y0*exp(Day / theta), data = RGRcurve, 
         start = c(y0 = 0.2, theta = 3.5))
summary(m)

Formula: RGR ~ y0 * exp(Day/theta)

Parameters:
      Estimate Std. Error t value Pr(>|t|)    
y0     0.16372    0.01421   11.52 4.03e-14 ***
theta  3.76454    0.17550   21.45  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.101 on 39 degrees of freedom

Number of iterations to convergence: 4 
Achieved convergence tolerance: 3.967e-09

Ejercicio 3

Escriba código R que realice una diagrama de dispersión de los datos y agregue la curva del modelo ajustado con el comando nls en el código anterior.

Código solución del gráfico

# Se extrae la estimacion de los parametros
mpar <- coef(m)
mpar  # se imprime para revisar

# Grafico base
gg <- ggplot(RGRcurve, aes(x = Day, y = RGR)) + geom_point() +
  scale_y_continuous(breaks = seq(0,2,0.1))

# Se agrega la curva del modelo ajustado:
gg + geom_function(fun = function(x) mpar[1]*exp(x/mpar[2]), 
                   color = "blue" )

Ejemplo 2: Puromicina

Considere los datos sobre la velocidad de reacción en función del sustrato para cultivos celulares tratados o no con Puromicina de la Figura 1C:

# Codigo R para activar los datos de Puromicina
puromy <- Puromycin
str(puromy)

'data.frame':   23 obs. of  3 variables:
 $ conc : num  0.02 0.02 0.06 0.06 0.11 0.11 0.22 0.22 0.56 0.56 ...
 $ rate : num  76 47 97 107 123 139 159 152 191 201 ...
 $ state: Factor w/ 2 levels "treated","untreated": 1 1 1 1 1 1 1 1 1 1 ...
 - attr(*, "reference")= chr "A1.3, p. 269"

Use estos datos para ajustar el modelo de Michaelis-Menten con cuatro parámetros dado por la función (2) ya presentada en la introducción:

\[\begin{equation*} f(x) = \dfrac{(V_m + \beta \cdot z) \cdot x}{K + x} \end{equation*}\]

donde \(x\) es la predictora que se refiere a la conc del sustrato, \(z\) es la predictora que se refiere a state, es decir, si el cultivo fue tratado o no con puromicina y \(f\) es la media de la respuesta (rate). Los parámetros son \(V_m\), \(K\) y \(\beta\). La predictora \(z\) se puede trabajar como una variable dummy (indicadora) tomando el valor de 1 si el cultivo es tratado y 0 si no lo es. Para crear esta variable indicadora en R use el siguiente código:

# Codigo R para crear la variable indicadora z:
library(tidyverse)
puromy <- mutate(puromy, z = ifelse(state == "treated", 1, 0))
str(puromy)

'data.frame':   23 obs. of  4 variables:
 $ conc : num  0.02 0.02 0.06 0.06 0.11 0.11 0.22 0.22 0.56 0.56 ...
 $ rate : num  76 47 97 107 123 139 159 152 191 201 ...
 $ state: Factor w/ 2 levels "treated","untreated": 1 1 1 1 1 1 1 1 1 1 ...
 $ z    : num  1 1 1 1 1 1 1 1 1 1 ...
 - attr(*, "reference")= chr "A1.3, p. 269"

Ejercicio 4

Al ojo anticipe valores aproximados para los parámetros del modelo (2); use estos como valores iniciales en el comando nls para ajustar el modelo (2).

Realice un gráfico que presente los puntos observados y las curvas del modelo para cada grupo de células (tratadas y no tratadas).

¿Podría realizar una prueba de hipótesis que compare la asintota de los dos grupos de células? Es decir:

\[H_0: V^{(z = 0)}_m = V^{(z = 1)}_m \quad \text{ contra } \quad H_1: V^{(z = 0)}_m \neq V^{(z = 1)}_m\]

Explique.

Código solución

# Se ejecuta el modelo para datos de puromy
m <- nls(rate ~ (Vm + beta*z)*conc / (K + conc), data = puromy, 
         start = c(Vm = 150, K = 0.1, beta = 50))
summary(m)

# Se realiza una tabla valores de las predictoras
nd <- expand.grid(
  conc = seq(0,1.3,0.05),
  z    = c(0,1)
)

# Se calculan las predicciones
nd$rate  <- predict(m, newdata = nd)

# Se agrega una columna de state
nd$state <- ifelse(nd$z == 1, "treated", "untreated")
head(nd, 3)

# Grafico
ggplot(puromy, aes(x = conc, y = rate, color = state)) + 
  geom_point() +
  geom_line(data = nd)

Referencias

Archontoulis, S. V., & Miguez, F. E. (2015). Nonlinear regression models and applications in agricultural research. Agronomy Journal, 107(2), 786–798. Retrieved from https://doi.org/10.2134/agronj2012.0506

Ritz, C., & Streibig, J. C. (2008). Nonlinear regression with R. Springer Science & Business Media.

Para saber más

Parámetros lineales y no lineales

La función \(f\) en (3) puede tener una mezcla de parámetros lineales y no lineales, y esta condición puede ser evaluada. Un parámetro es no lineal cuando la segunda derivada de la función \(f\) con respecto al parámetro no es cero (Archontoulis & Miguez, 2015).

Derivando en R

El comando D() del paquete mosaicCalc permite ejecutar derivadas con una sintaxis sencilla. Por ejemplo, suponga la siguiente función:

\(f(x) = a + bx^2\)

Recordando conceptos de cálculo diferencial, la derivada de esta función con respecto a \(x\) es:

\(f'(x) = \dfrac{df}{dx} = 2bx\)

Esta derivada se puede obtener con el comando D() (paquete mosaicCalc) como:

library(mosaicCalc)
dx <- D(a + b*x^2 ~ x)
dx

function (x, a = 0.8, b) 
b * (2 * x)

La 2da. derivada no es más que la derivada de la derivada. Así, la derivada de \(f'(x)\) está dada por:

\(f''(x) = \dfrac{d}{dx}\left( \dfrac{df}{dx}\right) = 2b\)

En R esto se logra aplicando una segunda vez el comando D() sobre el resultado previo de este mismo comando:

# Aplicando derivada a la derivada de la funcion anterior
D(dx(x) ~ x)

function (x, a = 0.8, b) 
b * 2

Para ahorrar escritura de código, el comando D() permite la siguiente sintaxis para ejecutar la 2da. derivada:

# Calculando directamente la 2da. derivada de la 
# funcion f original con respecto a x
D(a + b*x^2 ~ x & x)

function (x, a = 0.8, b) 
b * 2

Verificando linealidad sobre los parámetros

Podemos entonces usar el comando D() para verificar si un parámetro es lineal en cierta función \(f\). Por ejemplo en el modelo de Michaelis-Menten:

\[f(x) = \dfrac{V_m x}{K + x}\]

el parámetro \(V_m\) es lineal puesto que:

# 2da. derivada con respecto a Vm
D(Vm*x / (K + x) ~ Vm & Vm)

function (Vm, x, K) 
0

mientras que el parámetro \(K\) no es lineal puesto que:

# 2da. derivada con respecto a K
D(Vm*x / (K + x) ~ K & K)

function (K, Vm, x) 
Vm * x * (2 * (K + x))/((K + x)^2)^2

Aprovechando la condición de linealidad

Si al menos uno de los parámetros en un modelo es lineal de acuerdo a la regla de la 2da. derivada mostrada en la sección anterior, podemos aprovechar esta condición para facilitar el proceso de estimación por mínimos cuadrados.

El comando nls tiene el argumento algorithm en el cual se define el método de optimización numérica. Por defecto, dicho método es el de Gauss-Newton. Sin embargo, aquí también podemos indicar el método de Golub-Pereyra (plinear) de mínimos cuadrados parcialmente lineales. Al seleccionar este método, sólo debemos definir valores de inicio para aquellos parámetros no lineales, aunque al escribir la fórmula del modelo, también se debe omitir el parámetro. Revisar más sobre este aspecto en la sección 4.3 de Ritz & Streibig (2008).

Datos agrupados

Con datos agrupados nos refererimos a escenarios en donde se intenta comparar la función media entre dos o más grupos como en el ejemplo de la Puromycina de la Figura 1C. Observe que la comparación de los grupos se manejo usando el parámetro \(\beta\) y la variable dummy \(z\), la cual era 1 si state = treated y 0 si state = untreated. Otra manera de realizar la comparación o ajustar simultaneamente el mismo modelo para cada grupo (y luego tener la posibilidad de compararlos) es con la siguiente sintaxis:

# Codigo R para ejecutar un modelo separado para cada grupo:
# (Datos agrupados)
m <- nls(rate ~ Vm[state]*conc / (K[state] + conc), data = puromy,
         start = list(Vm = c(200, 160), K = c(0.01, 0.01)) )
summary(m)

Formula: rate ~ Vm[state] * conc/(K[state] + conc)

Parameters:
     Estimate Std. Error t value Pr(>|t|)    
Vm1 2.127e+02  6.608e+00  32.185  < 2e-16 ***
Vm2 1.603e+02  6.896e+00  23.242 2.04e-15 ***
K1  6.412e-02  7.877e-03   8.141 1.29e-07 ***
K2  4.771e-02  8.281e-03   5.761 1.50e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 10.4 on 19 degrees of freedom

Number of iterations to convergence: 8 
Achieved convergence tolerance: 3.623e-06

Observese que para cada parámetro, se reportan estimaciones para las dos categorías de la variable state. Revisar el capítulo 8 de Ritz & Streibig (2008) para más aspectos sobre este tema.