2024-06-26

Cargar librerías

require(pacman)
pacman::p_load(easynls, easyreg, knitr, htmlTable, readr, dplyr,  openxlsx, flextable)

Fijar directorio de trabajo

setwd('~/ACTIVIDAD DE TITULACION/REGRESION NO LINEAL')

Limpiar memoria

g <- gc(reset = T)
rm(list = ls())
options(scipen = 999, warn = -1)

Librería “easynls”

## Libreria “easyreg”

Funcion “Self-starting”

Ejercicio 1

Los datos corresponden al crecimiento acumulado (altura) de una especie de planta desde el día 1 hasta el día 16. El objetivo es modelar la evolución del crecimiento a través de la altura acumulada en función del tiempo. La altura se expresa en milímetros.

df<-as.data.frame(read.xlsx("regnolineal.xlsx", sheet=1))

Estructura del set de datos

glimpse(df)

Tabla de datos

flextable(df)

Gráfico de dispersión

plot(df)

Nota: La forma o tendencia de la curva es del tipo sigmoidea

Modelos para curvas sigmoideas

Probar ajustar a los siguientes modelos, disponibles en los paquetes easynls y easyreg de R:

  • Gompertz
  • Logístico
  • Degradación ruminal
  • Curva de Lactación
  • Modelo de Michaelis-Menten

Ajuste a modelo Gompertz

  • Búsqueda de valores iniciales
mod_Gomp<-nls(y~SSgompertz(x,a,b,c), data = df)
summary(mod_Gomp)

a=723; b=12; c=0.6

Ajuste con easyreg

regplot(df, model=10, start=c(723, 12, 0.6),
        digits=2, position=8, xlab="x", ylab="y")

Ajuste con easynls

nlsplot(df, model=10, start=c(723, 12, 0.6))

Ajuste a modelo Logístico

  • Búsqueda de valores iniciales
mod_Log<-nls(y~SSlogis(x,a,b,c), data=df)
summary(mod_Log)

a= 703; b= 6.5; c=1.5

Ajuste a modelo logístico

  • Ajuste con easynls
regplot(df, model=7, start=c(703, 6.5, 1.5), digits=2,
        position=8, xlab="x", ylab="y")

Ajuste con easynls

nlsplot(df, model=10, start=c(723, 12, 0.6))

Ajuste a modelo Degradación Ruminal

  • Valores iniciales
mod_Deg <-nls(y~SSasymp(x,a,b,c), data=df)
summary(mod_Deg)

a= 1122; -186; -2

Ajuste con easyreg

regplot(df, model=12, start=c(1122,-186,-2),
        digits=2, position=8, xlab="x", ylab="y")

Ajuste con easynls

nlsplot(df, model=12, start=c(1122,-186,-2))

Ajuste a modelo de Lactación

  • Valores iniciales
mod_Lac <-nls(y~SSlogis(x,a,b,c), data=df)
summary(mod_Lac)

a= 703; b= 6.5; c= 1.5

Ajuste con easyreg

regplot(df, model=11, start=c(703, 6.5, 1.5),
        digits=2, position=8, xlab="x", ylab="y")

Ajuste con easynls

nlsplot(df, model=11, start=c(703, 6.5, 1.5))

Ajuste a modelo de Michaelis-Menten

  • Búsqueda de valores iniciales
mod_micmen <- nls(y~SSmicmen(x,a,b), data=df)
summary(mod_micmen)

a= 7382; b= 124

Ajuste con easyreg

MM<-nls(y~(a*x)/(b+x),start=list(a=7382, b=124), data=df)
summary(MM)

Se pide el valor del AIC

AIC(MM)
Modelos <-c("Gompertz", "Logístico", "Degradación ruminal",
            "Curva de lactación", "Michaelis-Menten")
AIC <-c ("152.72", "146.40", "176.63", "148.35", "181.75")
tabla <- data.frame(Modelos, AIC)
flextable(tabla)

Comentarios

De los cuatro modelos sigmoideos evaluados, el que presentó menor AIC (criterio de información de Akaike) fue el LOGISTICO.

Se ajusta, entonces a ese modelo, quedando:

\[y=702,87(1+84,99(e^{0,69x}))^{-1}\]

Limpiar Ambiente de Trabajo

rm(list = ls())

EJERCICIO 2

Se presentan datos de la relación entre dosis de antibióticos y % de daño de una enfermedad bacteriana.

df<-as.data.frame(read.xlsx("regnolineal.xlsx", sheet=2))

Tabla de datos

flextable(df)

Visualizar tendencia de los datos

plot(df)

Nota: Se ajustará a una curva de tipo exponencial

Método de linealización

Linealizar el modelo exponencial

\[y=a*e^{bx} / log\] \[log{(y)}=log(a)+bx\]

Linealización para estimar los parámetros iniciales usando un modelo lineal

model <- lm(log(Daño) ~ Dosis, data=df)  
alpha <- exp(coef(model))[1]
beta <- coef(model)[2]

Parámetros iniciales

start <- list(alpha = alpha, beta = beta)
start

a = 99.94; b = -0.04

Ajuste a modelo Exponencial

a = 99,94 y b = -0,04

regplot(df, model=6, start=c(99.94, -0.04), digits=2,
        position=8, xlab="Dosis", ylab="Daño")
nlsplot(df, model=6, start=c(99.94, -0.04))

Modelo ajustado

\[y=100,9778*e^{-0,0443x}\]

Limpiar Ambiente de Trabajo

rm(list = ls())

EJERCICIO 3

Pruebe ajustar los datos de la tabla a:

1.- Función lineal \((y=a+bx)\)

2.- Función parabólica \((y=a+bx+cx^2)\)

3.- Función exponencial \((y=ab^x)\)

Datos del ejercicio

df<-as.data.frame(read.xlsx("regnolineal.xlsx", sheet=3))

Visualizar tendencia de los puntos

plot(df)

Ajuste linear

regplot(df, model=1, digits=2, position=8,
        xlab="X", ylab="Y")
nlsplot(df, model=1)

Ajuste parabólico

regplot(df, model=2, digits=2, position=8,
        xlab="X", ylab="Y")
nlsplot(df, model=2)

Ajuste a modelo Exponencial

  • Linealización para estimar los parámetros iniciales
model <- lm(log(y) ~ x, data=df)  
alpha <- exp(coef(model))[1]
beta <- coef(model)[2]

Parámetros iniciales

start <- list(alpha = alpha, beta = beta)
start

Nota: a = 0.82; b = 0.78

Ajuste a modelo Exponencial utilizando los parámetros iniciales determinados

a = 0.82 y b = 0.78

regplot(df, model=6, start=c(0.82, 0.78),
        digits=2, position=8, xlab="X", ylab="Y")
nlsplot(df, model=6, start=c(0.82, 0.78))

Resumen de CME y AIC

resumen = data.frame(
"MODELO" = c("Función lineal", "Función parabólica",
             "Función exponencial"),
"CME" = c("5.76", "0.06", "3.38"), 
"AIC" = c("26.39", "3.30", "24.34"))
flextable(resumen)

¿Cuál de los tres modeles presenta el mejor ajuste?

  • Comentarios

De los tres modelos evaluados, el que presentó menor AIC (criterio de información de Akaike) fue la FUNCION PARABOLICA.

Se ajusta, entonces a ese modelo, quedando:

\[y=-0,45+0,49x+1,14x{^2}\]