TALLER REGRESION LINEAL Y NO LINEAL

Objetivos. Ajustar datos X,Y a modelos lineales y no lineales

Primeros pasos con regresión no lineal (nls) con R

PROCEDIMIENTO

Cargue los datos del dataset “HUMEDADES” del archivo de Excel llamado “Datos_2024_SEM_1.xlsx

library(readxl)
h <- read_excel("DATOS_2024_SEM_1.xlsx", sheet = "HUMEDADES")
h

#2. Seleccione los datos en Excel de tiempo y humedad para la temperatura de 40 grados y cree un dataset h40

h40 <- subset(h, temperatura == 40)
h40

3. Haga lo mismo para las demás temperaturas y cree los datasets h50, h60, h70 y h80

h50 <- subset(h, temperatura == 50)

Ajuste con Regresión Lineal Simple para mirar la trampa estadística

Intente ajustar los lineas suponiendo una relación lineal (Recuerde que el secado sigue una caída exponencial de acuerdo a la ley de enfriamiento de Newton)

modelo1 <- lm(humedad~tiempo, data = h40)
summary(modelo1)

## 
## Call:
## lm(formula = humedad ~ tiempo, data = h40)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.09104 -0.06348 -0.01044  0.04960  0.15036 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.849644   0.015197   55.91   <2e-16 ***
## tiempo      -0.054933   0.001895  -28.99   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.07058 on 47 degrees of freedom
## Multiple R-squared:  0.947,  Adjusted R-squared:  0.9459 
## F-statistic: 840.3 on 1 and 47 DF,  p-value: < 2.2e-16

plot(h40$tiempo, h40$humedad)
abline(modelo1)

Validando la linealidad de la regresión con el test de Durbin-Watson con la función dwtest() del paquete lmtest

Durbin-Watson Test Description Performs the Durbin-Watson test for autocorrelation of disturbances.

El test de Durbin-Watson es una prueba estadística utilizada para evaluar la presencia de autocorrelación en los residuos de un modelo de regresión.

library(lmtest)

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

dwtest(modelo1)

## 
##  Durbin-Watson test
## 
## data:  modelo1
## DW = 0.021528, p-value < 2.2e-16
## alternative hypothesis: true autocorrelation is greater than 0

Aquí están los resultados de la prueba aplicada a tu modelo de regresión lineal de humedad y tiempo:

DW (Durbin-Watson): El valor calculado del estadístico de Durbin-Watson es 0.021528. p-valor: El p-valor asociado con el estadístico es menor que 2.2e-16 (esencialmente cero). Interpretación:

La hipótesis nula del test de Durbin-Watson es que no hay autocorrelación en los residuos (es decir, los errores son independientes). La hipótesis alternativa es que existe autocorrelación positiva (los errores están correlacionados positivamente). Dado que el p-valor es extremadamente pequeño, tenemos evidencia significativa para rechazar la hipótesis nula. Esto sugiere que existe autocorrelación positiva en los residuos de tu modelo. En otras palabras, los errores no son completamente independientes y están correlacionados en el tiempo.

El valor de DW = 0.021528, indica que la regresión lineal simple no está capturando la totalidad de la variación de la humedad en el tiempo

Ajuste a cada uno de los siguientes modelos:

Regresiones no Lineales

Función nls del paquete stats para realizar Regresiones no Lineales

nls {stats} R Documentation Nonlinear Least Squares

Description Determine the nonlinear (weighted) least-squares estimates of the parameters of a nonlinear model.

MODELO DE HENDERSON & PABIS

\[ r1 <- nls(humedad \sim a * exp(-k * tiempo), trace=TRUE) \]

r1 <- nls(humedad ~ a*exp(k*tiempo), start = list(a=0.1, k=0.1), trace = TRUE, data = h40)

## 13.13646    (2.11e+00): par = (0.1 0.1)
## 4.731108    (4.25e+00): par = (0.6047614 -0.3104843)
## 0.5975009   (3.49e+00): par = (0.749213 -0.0912802)
## 0.04661037  (1.73e-01): par = (0.9380655 -0.1363996)
## 0.04526311  (2.29e-04): par = (0.9472385 -0.1370171)
## 0.04526311  (2.46e-07): par = (0.9472417 -0.1370139)

summary(r1)

## 
## Formula: humedad ~ a * exp(k * tiempo)
## 
## Parameters:
##    Estimate Std. Error t value Pr(>|t|)    
## a  0.947242   0.008662  109.36   <2e-16 ***
## k -0.137014   0.003069  -44.64   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.03103 on 47 degrees of freedom
## 
## Number of iterations to convergence: 5 
## Achieved convergence tolerance: 2.459e-07

MODELO DE LEWIS (a=1)

\[ r2 <- nls(humedad \sim exp(-k * tiempo), trace=TRUE) \]

r2 <- nls(humedad~exp(k*tiempo), start = list(k=0.1), trace = TRUE, data = h40)

## 236.4895    (5.14e+00): par = (0.1)
## 35.20167    (4.63e+00): par = (0.0295827)
## 4.797534    (4.28e+00): par = (-0.04115577)
## 0.4995013   (2.09e+00): par = (-0.1034074)
## 0.09285479  (3.87e-01): par = (-0.1392002)
## 0.08032712  (2.75e-02): par = (-0.1476878)
## 0.08026391  (1.20e-03): par = (-0.1483315)
## 0.08026379  (4.96e-05): par = (-0.1483596)
## 0.08026379  (2.05e-06): par = (-0.1483608)

summary(r2)

## 
## Formula: humedad ~ exp(k * tiempo)
## 
## Parameters:
##    Estimate Std. Error t value Pr(>|t|)    
## k -0.148361   0.003395   -43.7   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.04089 on 48 degrees of freedom
## 
## Number of iterations to convergence: 8 
## Achieved convergence tolerance: 2.052e-06

MODELO DE PAGE

\[ r3 <- nls(humedad \sim exp(-k * tiempo^n), trace=TRUE) \]

MODELO DE LOS DOS TÉRMINOS

\[ r4 <- nls(humedad \sim A1*exp(-B * tiempo) + A2*exp(B*tiempo), trace=TRUE) \]

MODELO DE PAGE MODIFICADO

\[ r5 <- nls(humedad \sim exp((k * tiempo)^n), trace=TRUE) \]

MODELO DE THOMPSON

\[ r6 <- nls(humedad \sim a*log(tiempo)+b*(log(tiempo))^2, trace=TRUE) \]

MODELO DE WANG Y SINGH

\[ r7 <- nls(humedad \sim 1+a*tiempo+b*tiempo^2, trace=TRUE) \]

ANOVA DE LOS MODELOS

hacemos un anova para ver cual modelo se ajusta mejor, significancia de los modelos

LA BONDAD DEL AJUSTE SERÁ PARA EL MAYOR R^2 Y LOS MENORES ssr y sse

a=anova(r1, r2, r3, r4,r7)

SIGNIFICANCIA DEL MODELO

LA BONDAD DEL AJUSTE SERÁ PARA EL MAYOR R² Y LOS MENORES ssr y sse

EL CALCULO DE R²: VER MONTGOMERY PAGINA 464 \[ R^2=SSr/SSt=1-SSe/SSt \]

\[ r_2r1=1-sum(resid(r1)^2)/sum((humedad-mean(humedad))^2) \]

R^2: VER MONTGOMERY PAGINA 464

\[ sser1=sum(resid(r1)^2) \]

\[ sstr1=sum((humedad-mean(humedad))^2) \]

\[ ssrr1=sstr1-sser1 \] VER MONTGOMERY PAGINA 451

UN EJERCICIO CON LOS PRECIOS DE LAS CASAS DE BOSTON Y EL ESTATUS

Visualice y analice la relación entre las variables:

lstat: lower status of the population (percent). medv: median value of owner-occupied homes in $1000s.

Visualice
Regresione con lm()
Diagnostique la linealidad con dwtest() de lmtest y bptest()

UN EJERCICIO CON LOS DATOS DE PRODUCCIÓN DE CLORO

Carque el archivo: produccion_de_cloro y realice una regresión lineal simple de la producción de cloro en función de las semanas de producción

Cargue el archivo en su carpeta de trabajo
Obtenga las variables x,y. x = semanas y y = cloro
Grafique x,y
Realice la regresión Lineal Simple
Diagnostique con dwtest() y bptest()

semanas=c(8,8,10,10,10,10,12,12,12,12,14,14,14,16,16,16,18,18,
        20,20,20,22,22,22,24,24,24,26,26,26,28,28,30,30,30,32,
        32,34,36,36,38,38,40,42)
cloro = c(0.49,0.49,0.48,0.47,0.48,0.47,0.46,0.46,0.45,0.43,0.45,
        0.43,0.43,0.44,0.43,0.43,0.46,0.45,0.42,0.42,0.43,0.41,
        0.41,0.40,0.42,0.40,0.40,0.41,0.40,0.41,0.41,0.40,0.40,
        0.40,0.38,0.41,0.40,0.40,0.41,0.38,0.40,0.40,0.39,0.39)
plot(semanas, cloro)

r <- lm(cloro~semanas)
summary(r)

## 
## Call:
## lm(formula = cloro ~ semanas)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.025741 -0.012042 -0.001608  0.012034  0.026224 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.4855103  0.0058907   82.42  < 2e-16 ***
## semanas     -0.0027168  0.0002431  -11.18 3.67e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01539 on 42 degrees of freedom
## Multiple R-squared:  0.7483, Adjusted R-squared:  0.7423 
## F-statistic: 124.9 on 1 and 42 DF,  p-value: 3.675e-14

library(lmtest)
dwtest(r)

## 
##  Durbin-Watson test
## 
## data:  r
## DW = 0.99208, p-value = 7.593e-05
## alternative hypothesis: true autocorrelation is greater than 0

Ajuste no lineal

Se desea ajustar el siguiente modelo a los datos:

\[ cloro = a+(0.49-a)*\exp^{(-b*(semanas-8))} \]

Este modelo, sugerido por un experto en el área, contiene dos incógnitas: a, el valor asintótico basal que se alcanza con valores grandes de semana, y b, la tasa exponencial de decaimiento.

modelocloro <- nls(cloro~a+(0.49-a)*exp(-b*(semanas-8)))

## Warning in nls(cloro ~ a + (0.49 - a) * exp(-b * (semanas - 8))): No starting values specified for some parameters.
## Initializing 'a', 'b' to '1.'.
## Consider specifying 'start' or using a selfStart model

summary(modelocloro)

## 
## Formula: cloro ~ a + (0.49 - a) * exp(-b * (semanas - 8))
## 
## Parameters:
##   Estimate Std. Error t value Pr(>|t|)    
## a 0.390140   0.005045  77.333  < 2e-16 ***
## b 0.101633   0.013360   7.607 1.99e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01091 on 42 degrees of freedom
## 
## Number of iterations to convergence: 8 
## Achieved convergence tolerance: 5.738e-07

Código de Python para resolver la regresión no lineal de la producción de cloro

NOTA: Para incrustar código de Python en un documento RMarkdown, se requiere instalar el paquete reticulate y además instalar todos los módulos de python que use el código de Python. Para instalar un paquete de R se debe lanzar en la consola

install.packages("reticulate")

La instalación de los módulos de Python se hacen en Posit Cloud en el panel inferior izquierdo, en la pestaña “Terminal”, utilizando el programa pip de Python

PIP (acrónimo de “Paquete de Instalación de Python” en inglés, o Python Package Installer) es una herramienta esencial para cualquier desarrollador que trabaje con Python.

Facilita la instalación, actualización y gestión de bibliotecas y módulos de Python desde el Python Package Index (PyPI) y otros repositorios. A través de PIP, los desarrolladores pueden acceder a una amplia variedad de paquetes de software que complementan y enriquecen la funcionalidad de Python.

Algunas acciones comunes que puedes realizar con PIP son:

Instalar paquetes: Utiliza pip install nombre_del_paquete para descargar e instalar un paquete específico.
Buscar paquetes disponibles: Utiliza pip search nombre_paquete para explorar los paquetes disponibles en PyPI.
Eliminar paquetes: Si ya no necesitas un paquete, puedes desinstalarlo con pip uninstall nombre_del_paquete.
Ver la lista de paquetes instalados: Ejecuta pip list para obtener una lista de los paquetes instalados en tu entorno Python. En resumen, PIP simplifica la gestión de paquetes y bibliotecas de terceros, lo que mejora la eficiencia del proceso de desarrollo en Python.

El código en Python es el siguiente:

import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import curve_fit

def f(x, a, b):
    return a + (0.49 -a)*np.exp(-b*(x-8))

semanas=[8,8,10,10,10,10,12,12,12,12,14,14,14,16,16,16,18,18,
        20,20,20,22,22,22,24,24,24,26,26,26,28,28,30,30,30,32,
        32,34,36,36,38,38,40,42]

cloro = [0.49,0.49,0.48,0.47,0.48,0.47,0.46,0.46,0.45,0.43,0.45,
        0.43,0.43,0.44,0.43,0.43,0.46,0.45,0.42,0.42,0.43,0.41,
        0.41,0.40,0.42,0.40,0.40,0.41,0.40,0.41,0.41,0.40,0.40,
        0.40,0.38,0.41,0.40,0.40,0.41,0.38,0.40,0.40,0.39,0.39]

x = semanas
y = cloro

res,cov = curve_fit(f,x,y)

print(res)

## [0.39014001 0.1016327 ]

print(cov)

## [[2.54520196e-05 5.98436879e-05]
##  [5.98436879e-05 1.78493841e-04]]


xx = np.linspace(8,40,50)
fig,axes=plt.subplots()
axes.scatter(x,y)
axes.plot(xx,f(xx,0.39014001,0.1016327))

NOTA: Utilice la función nls() de R y realice el ajuste no lineal del cloro en función de las semanas

Regresión Lineal y no Lineal

Justo

2024-04-23