1.ANÁLISIS DE TENDENCIAS

Lo primero a realizar es el análisis de tendencias. Este se divide el dos: tendencias de medición y tendencias de envasado. El primero busca evaluar si la medición estuvo bien y si no se requiere realizar algún ajuste o empleo de una técnica estadística diferente a la propuesta por ISO 35.

El segundo tiene como objetivo determinar si se realizó un proceso adecuado de envasado de cada una de las unidades.

1.1 Tendencia de medición

Para este caso, se construye una gráfica de la respuesta instrumental vs el número de medición. NOTA: en este caso omití los controles que ingresaste, pues se evidencia que tienen problemas, aparentemente de estabilidad… lo podemos discutir después.

1.1.1 THC

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## 
## Call:
## lm(formula = Respuesta ~ Medicion, data = datos_excel)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -41.754  -1.077   2.960   6.273  16.292 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 391.4754     4.1489  94.357   <2e-16 ***
## Medicion      0.1488     0.1449   1.027    0.311    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 12.67 on 37 degrees of freedom
## Multiple R-squared:  0.02771,    Adjusted R-squared:  0.001436 
## F-statistic: 1.055 on 1 and 37 DF,  p-value: 0.3111

El análisis de regresión, permite concluir que el modelo lineal no se ajusta a los datos (p value>0.05), por consiguiente se concluye que no hay tendencia de medición y se puede proceder a realizar la estimación de incertidumbre.

1.1.2 CBN

## 
## Call:
## lm(formula = RespuestaCBN ~ Medicion, data = datos_excel)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -12.2108  -1.5998   0.9495   1.7762   6.0302 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 133.47120    1.30374 102.375   <2e-16 ***
## Medicion      0.04109    0.04554   0.902    0.373    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.98 on 37 degrees of freedom
## Multiple R-squared:  0.02153,    Adjusted R-squared:  -0.004914 
## F-statistic: 0.8142 on 1 and 37 DF,  p-value: 0.3727

El análisis de regresión, permite concluir que el modelo lineal no se ajusta a los datos (p value>0.05), por consiguiente se concluye que no hay tendencia de medición y se puede proceder a realizar la estimación de incertidumbre.

** NOTA: En ambos casos parece existir problemas con las mediciones que se encuentran en la parte inferior de la gráfica, correspondientes a 18, 22 y 27. Sin embargo, en estos estudios esta casi prohibido realizar descarte de datos por pruebas estadísticas. Lo que se debe hacer, es revisar bien las integraciones,las señales cromatográficas y si hay algo mal se podría corregir o descartar los datos. Por ahora, debido al tiempo, voy a proceder y revisar si tal vez es necesario emplear métodos de estadística robusta.**

1.2 Tendencia de envasado

Para este caso, se construye una gráfica de la promedio para cada botella vs el número de botella

1.2.1 TCH

## # A tibble: 8 × 3
##   Botella Promedio Desviacion
##     <dbl>    <dbl>      <dbl>
## 1       7     388.      23.0 
## 2      20     396.       6.53
## 3      30     392.      16.9 
## 4      38     398.       6.19
## 5      45     400.       6.85
## 6      57     390.      18.1 
## 7      69     400.       7.08
## 8      76     399.       4.17

## 
## Call:
## lm(formula = Promedio ~ Botella, data = resumen_botellas)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.1670 -2.3637  0.9073  3.2219  4.0252 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 390.12909    3.20780 121.619 2.08e-11 ***
## Botella       0.12102    0.06652   1.819    0.119    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.198 on 6 degrees of freedom
## Multiple R-squared:  0.3555, Adjusted R-squared:  0.2481 
## F-statistic:  3.31 on 1 and 6 DF,  p-value: 0.1187

Nuevamente se encuentra que no hay tedencia de envasado significativa, por lo cual se concluye que el proceso de envasado esta ok.

1.2.2 CBN

## # A tibble: 8 × 3
##   Botella Promedio Desviacion
##     <dbl>    <dbl>      <dbl>
## 1       7     132.      23.0 
## 2      20     135.       6.53
## 3      30     134.      16.9 
## 4      38     134.       6.19
## 5      45     138.       6.85
## 6      57     133       18.1 
## 7      69     136.       7.08
## 8      76     136.       4.17

## 
## Call:
## lm(formula = Promedio ~ Botella, data = resumen_botellas2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.1462 -0.7178 -0.1407  0.3355  3.1796 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 132.88646    1.32238 100.491 6.54e-11 ***
## Botella       0.03964    0.02742   1.446    0.198    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.731 on 6 degrees of freedom
## Multiple R-squared:  0.2584, Adjusted R-squared:  0.1347 
## F-statistic:  2.09 on 1 and 6 DF,  p-value: 0.1984

Nuevamente se encuentra que no hay tedencia de envasado significativa, por lo cual se concluye que el proceso de envasado esta ok.

2. ESTIMACIÓN DE INCERTIDUMBRE

De acuerdo al ISO 35, se estima la incertidumbre a través de un ANOVA.

2.1 Estimación de incertidumbre para THC

##             Df Sum Sq Mean Sq F value Pr(>F)
## Botella      7    802   114.5   0.669  0.696
## Residuals   31   5304   171.1
##  Cuadrado Medio del Error (MSE): 171.0901 
##  Cuadrado Medio Entre Grupos (MSG): 114.5051 
##  Número de réplicas : 5

En este caso, da que el CME es mayor que CMD, es decir para estudios futuros es mejor realizar un mayor número de réplicas o en mejorar el método de medición. Aproximadamente, se requiere una precisión del método de unas 3 ó 4 veces mejor. Sin emabargo, la ISO 35 da una opción para estos casos que básicamente es atribuir la incertidumbre del método de medición a la homogeneidad del material. Vamos a tomar esta vía

A continuación coloco el código para los este escenarios:

promedioTHC<- mean(datos_excel$Respuesta)

   # Escenario 2: cuando la precisión del método no es la "adecuada", es decir CMentre < CMdentro
  
    uhom2 <- sqrt((CMdentro / Nreplicas)*(sqrt(2 / GLerror)))  # Realiza la estimación de la incertidumbre cuando la precisión es insuficiente
    cat(" Incertidumbre de Homogeneidad (%):",uhom2/promedioTHC*100, "\n")  # Muestra el resultado de incertiumbre 
##  Incertidumbre de Homogeneidad (%): 0.7459951

Esta sería la incertidumbre por homogeneidad.

Ahora procedemos a validar el modelo. Es decir, se evalua si los residuos siguen una distribución normal y son homocedasticos. PAra ello empleó la prueba de Bartlett y Shapiro-Wilk.

# Evaluación de supuestos:  Prueba de Bartlett para la homogeneidad de varianzas y  Shapiro-Wilk para normalidad

bartlett_result <- bartlett.test(Anova$residuals, g = datos_excel$Botella)
shapiro_test <- shapiro.test(Anova$residuals)
print(bartlett_result)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  Anova$residuals and datos_excel$Botella
## Bartlett's K-squared = 18.523, df = 7, p-value = 0.009821
print(shapiro_test)
## 
##  Shapiro-Wilk normality test
## 
## data:  Anova$residuals
## W = 0.90297, p-value = 0.002686

En ambos casos se encuentra que los datos no cumplen los supuestos, lo anterior puede ser ( ESTOY SEGURO) que es por los datos “anómalos” que se tienen. Nota: no necesariamente son datos anómalos, por ello las comillas. Es decir, durante el proceso de medición pasan estas cosas y es la variación natural de los datos.

Ahora, vamos a suponer que estos datos no son anomalos. Para este caso, vamos a usar un modelo de efectos aleatorios, el cual es robusto a este tipo de cosas.

## # A tibble: 8 × 3
##   Botella Promedio Varianza
##   <fct>      <dbl>    <dbl>
## 1 7           388.    530. 
## 2 20          396.     42.6
## 3 30          392.    285. 
## 4 38          398.     38.3
## 5 45          400.     46.9
## 6 57          390.    327. 
## 7 69          400.     50.1
## 8 76          399.     17.4
## 
## Random-Effects Model (k = 8; tau^2 estimator: DL)
## 
##   logLik  deviance       AIC       BIC      AICc   
## -25.6651    0.8963   55.3303   55.4892   57.7303   
## 
## tau^2 (estimated amount of total heterogeneity): 0 (SE = 30.5931)
## tau (square root of estimated tau^2 value):      0
## I^2 (total heterogeneity / total variability):   0.00%
## H^2 (total variability / sampling variability):  1.00
## 
## Test for Heterogeneity:
## Q(df = 7) = 0.8963, p-val = 0.9963
## 
## Model Results:
## 
## estimate      se      zval    pval     ci.lb     ci.ub      
## 398.2442  2.5246  157.7485  <.0001  393.2962  403.1923  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

En este caso, se observa que el valor de TAU es cero. Es decir se encuentra qu ela homogeneidad del MR es muy buena. De hecho, por ello es posible que el ANOVA no la pudiera detectar en la vía clásica.

Debido a todo el experimento de homogeneidad, el enfoque se tratar de “sobreestimar” la incertidumbre, tenemos que usar la desviación estándar del TAU (SE en la salida).

## uhom (%)= 1.399598

Este sería el resultado final.

2.2 Estimación de incertidumbre para CBN

##             Df Sum Sq Mean Sq F value Pr(>F)
## Botella      7  110.2   15.74   0.998  0.451
## Residuals   31  488.9   15.77
##  Cuadrado Medio del Error (MSE): 15.76974 
##  Cuadrado Medio Entre Grupos (MSG): 15.74396 
##  Número de réplicas : 5

Estamos en el mismo caso del THC, así que ya saben la historia

A continuación coloco el código para los este escenario:

promedioCBN<- mean(datos_excel$RespuestaCBN)

   # Escenario 2: cuando la precisión del método no es la "adecuada", es decir CMentre < CMdentro
  
    uhom2 <- sqrt((CMdentro / Nreplicas)*(sqrt(2 / GLerror)))  # Realiza la estimación de la incertidumbre cuando la precisión es insuficiente
    cat(" Incertidumbre de Homogeneidad (%):",uhom2/promedioTHC*100, "\n")  # Muestra el resultado de incertiumbre 
##  Incertidumbre de Homogeneidad (%): 0.2264831

Ahora procedemos a validar el modelo. Es decir, se evalua si los residuos siguen una distribución normal y son homocedasticos. PAra ello empleó la prueba de Bartlett y Shapiro-Wilk.

# Evaluación de supuestos:  Prueba de Bartlett para la homogeneidad de varianzas y  Shapiro-Wilk para normalidad

bartlett_result <- bartlett.test(Anova$residuals, g = datos_excel$Botella)
shapiro_test <- shapiro.test(Anova$residuals)
print(bartlett_result)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  Anova$residuals and datos_excel$Botella
## Bartlett's K-squared = 21.635, df = 7, p-value = 0.002935
print(shapiro_test)
## 
##  Shapiro-Wilk normality test
## 
## data:  Anova$residuals
## W = 0.89687, p-value = 0.001797

Misma situación, así que solucionemos el problema con un modelo de efectos aleatorios, se me olvidó contar antes que DL corresponde a DerSimmonian Laird, este es un método que permite estimar TAU.

## # A tibble: 8 × 3
##   Botella Promedio Varianza
##   <fct>      <dbl>    <dbl>
## 1 7           134.   35.6  
## 2 20          134.   10.6  
## 3 30          134.   32.2  
## 4 38          134.    0.315
## 5 45          134.    3.50 
## 6 57          134.   32.4  
## 7 69          134.    2.79 
## 8 76          134.    5.73
## 
## Random-Effects Model (k = 8; tau^2 estimator: DL)
## 
##   logLik  deviance       AIC       BIC      AICc   
## -15.2271   -0.0000   34.4541   34.6130   36.8541   
## 
## tau^2 (estimated amount of total heterogeneity): 0 (SE = 2.1946)
## tau (square root of estimated tau^2 value):      0
## I^2 (total heterogeneity / total variability):   0.00%
## H^2 (total variability / sampling variability):  1.00
## 
## Test for Heterogeneity:
## Q(df = 7) = 0.0000, p-val = 1.0000
## 
## Model Results:
## 
## estimate      se      zval    pval     ci.lb     ci.ub      
## 134.4974  0.4893  274.9032  <.0001  133.5385  135.4564  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Lo mismo que en el anterior.

## uhom CBN (%)= 1.101458

Este sería el resultado final para CBN

# 3. ESTABILIDAD

La estimación de estabilidad del MR si se realiza de una manera más sencilla. En principio, lo que vamos a evaluar si el modelo de regresión lineal se ajusta a los datos. Es decir una gráfica de respuesta vs tiempo.

Para este caso, normalicé todas las respuestas al tiempo cero, del manera que todo esta en % de concentración inicial.

library(readxl)
Datos <- read_excel("~/Desktop/Datos.xlsx", 
    sheet = "Estabilidad")

 # Previamente debemos  convertir la columna "Tiempo" en un factor para que R la interprete como una variable categórica.
    Datos$Tiempo <- as.factor(Datos$Tiempo)

En este caso, voy a trabajar directamente con minimos cuadrados ponderados, pues se evidencia que hay heterocedasticidad en los datos.

3.1 Estimación de incertidumbre por estabilidad para THC

# Para MCP requerimos estimar los pesos ( inverso de varianza)

      resumen_tiempo <- Datos %>% # %>% es el operador de paquete dplyr para encadenar varias 
      group_by(Tiempo) %>% # Esto es para agrupar por el código de la botella
      summarize(Promedio = mean(Respuesta1), Varianza = sd(Respuesta1)^2)
      Datos <- Datos %>%  #Usamos el operador %>% para encadenar operaciones en Datos. 
        left_join(resumen_tiempo, by = "Tiempo") %>%  
        mutate(pesos1 = Varianza)  # Creamos una nueva columna llamada "pesos" en Datos utilizando los valores de "Varianza" de resumen_tiempo.
      
   
      # Verificar que el Datframe cambio
      print(Datos, n=11)
## # A tibble: 30 × 8
##    Medición Etiquete Tiempo Respuesta1 Respuesta2 Promedio Varianza   pesos1
##       <dbl> <chr>    <fct>       <dbl>      <dbl>    <dbl>    <dbl>    <dbl>
##  1        1 14-2     48          1.00       1.01      1.00 0.000212 0.000212
##  2        2 54-2     0           0.969      0.980     1.00 0.000336 0.000336
##  3        3 1-1      72          0.961      0.968     1.00 0.000524 0.000524
##  4        4 47-2     48          1.01       1.01      1.00 0.000212 0.000212
##  5        5 73-3     24          0.993      0.990     1.00 0.000304 0.000304
##  6        6 70-2     120         1.02       1.02      1.00 0.000279 0.000279
##  7        7 55-3     72          1.01       1.03      1.00 0.000524 0.000524
##  8        8 47-1     48          0.972      0.999     1.00 0.000212 0.000212
##  9        9 51-3     120         1.02       1.02      1.00 0.000279 0.000279
## 10       10 14-3     48          1.01       1.01      1.00 0.000212 0.000212
## 11       11 35-3     0           1.02       1.01      1.00 0.000336 0.000336
## # ℹ 19 more rows

Con los pesos procedemos a estimar el modelo de regresión lineal ponderada:

 # Paso 5: Realizar el análisis de regresión ponderada (WLS)
      
      # Primero convertimos  la columna "Tiempo" en un factor númerico.
      
      Datos$Tiempo <- as.numeric(Datos$Tiempo)
      
      # Segundo estimamos los parámetros del modelo
      
      modelo_wls <- lm(Respuesta1 ~ Tiempo, data = Datos, weights = 1 / pesos1)

      # Resumen del modelo
      summary(modelo_wls)
## 
## Call:
## lm(formula = Respuesta1 ~ Tiempo, data = Datos, weights = 1/pesos1)
## 
## Weighted Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.0118 -0.4835  0.2605  0.8449  1.1946 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.9994708  0.0072042 138.735   <2e-16 ***
## Tiempo      0.0006492  0.0021872   0.297    0.769    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9499 on 28 degrees of freedom
## Multiple R-squared:  0.003137,   Adjusted R-squared:  -0.03247 
## F-statistic: 0.08811 on 1 and 28 DF,  p-value: 0.7688

En este caso, se encuentra que la regresión no es significativa, así que empleamos la desviación estándar de la pendiente para estimar la incertidumbre (ISO 35). Para esto, múltiplicamos la desviación estándar de la pendiente por el tiempo que se requería en transpote ( 72h que es lo que típicamente se demora una muestra en ser transportada), correspondería a el valor de % incertidumbre de estabilidad. A continuación muestro como hacerlo.

##  Extraer el valor p,  la pendiente y de la desviación estándar de la pendiente

      p_value <- summary(modelo_wls)$coefficients["Tiempo", "Pr(>|t|)"]
      pendiente <- coef(modelo_wls)["Tiempo"]
      desviacion_pendiente <- sqrt(vcov(modelo_wls)["Tiempo", "Tiempo"])
      
      # Mostrar los resultados
      cat(" Valor p:", p_value, "\n","Pendiente:", pendiente, "\n","Desviación estándar de la pendiente:", desviacion_pendiente, "\n")
##  Valor p: 0.7687847 
##  Pendiente: 0.0006492439 
##  Desviación estándar de la pendiente: 0.002187234
# Estimar el valor de la incertidumbre
      
      # Primero se verifica el valor de  p ( si es menor  o mayor que 0.05)
      
      if (p_value < 0.05) {
        # Calcular el factor uest con la pendiente y 10 días
        uest <- abs(pendiente) * 72 #Donde 72 corresponde al número de días
      } else {
        # Calcular el factor uest con la desviación estándar de la pendiente y 10 días
        uest <- desviacion_pendiente * 72  #Donde 72 corresponde al número de días
      }
            # Mostrar el valor de uest
      cat("Factor uest:", uest, "%\n")
## Factor uest: 0.1574809 %

Este es un valor de incertidumbre muy alto para esto, por lo cual es mejor transportarla a 4 ºC, para lo cual apoyados en la regla Q, establecemos que la inestabilidad. En esencia, acá se asume que por cada 10 ºC que disminuya la temperatura, la constante cinética disminuye cerca de 2 veces. En este caso sería de 7.2 veces la disminución (40 ºC a 4 ºC, es decir 36 ºC menos). 7.2 es genial porque es sólo dividir en 10 esta incertidumbre obtenida y la final sería 0.256%, para una temperatura de transporte de 4 ºC por máximo 3 días…

3.1 Estimación de incertidumbre por estabilidad para CBN

# Importamos nuevamente 
Datos <- read_excel("~/Desktop/Datos.xlsx", 
    sheet = "Estabilidad")

# Para MCP requerimos estimar los pesos ( inverso de varianza)

      resumen_tiempo <- Datos %>% # %>% es el operador de paquete dplyr para encadenar varias 
      group_by(Tiempo) %>% # Esto es para agrupar por el código de la botella
      summarize(Promedio = mean(Respuesta2), Varianza = sd(Respuesta2)^2)
      Datos <- Datos %>%  #Usamos el operador %>% para encadenar operaciones en Datos. 
        left_join(resumen_tiempo, by = "Tiempo") %>%  
        mutate(pesos2 = Varianza)  # Creamos una nueva columna llamada "pesos" en Datos utilizando los valores de "Varianza" de resumen_tiempo.
      
   
      # Verificar que el Datframe cambio
      print(Datos, n=11)
## # A tibble: 30 × 8
##    Medición Etiquete Tiempo Respuesta1 Respuesta2 Promedio  Varianza    pesos2
##       <dbl> <chr>     <dbl>      <dbl>      <dbl>    <dbl>     <dbl>     <dbl>
##  1        1 14-2         48      1.00       1.01      1.01 0.0000640 0.0000640
##  2        2 54-2          0      0.969      0.980     1.00 0.000261  0.000261 
##  3        3 1-1          72      0.961      0.968     1.02 0.000890  0.000890 
##  4        4 47-2         48      1.01       1.01      1.01 0.0000640 0.0000640
##  5        5 73-3         24      0.993      0.990     1.02 0.000327  0.000327 
##  6        6 70-2        120      1.02       1.02      1.01 0.0000711 0.0000711
##  7        7 55-3         72      1.01       1.03      1.02 0.000890  0.000890 
##  8        8 47-1         48      0.972      0.999     1.01 0.0000640 0.0000640
##  9        9 51-3        120      1.02       1.02      1.01 0.0000711 0.0000711
## 10       10 14-3         48      1.01       1.01      1.01 0.0000640 0.0000640
## 11       11 35-3          0      1.02       1.01      1.00 0.000261  0.000261 
## # ℹ 19 more rows

Con los pesos procedemos a estimar el modelo de regresión lineal ponderada:

 # Paso 5: Realizar el análisis de regresión ponderada (WLS)
      
 
      # Segundo estimamos los parámetros del modelo
      
      modelo_wls <- lm(Respuesta2 ~ Tiempo, data = Datos, weights = 1 / pesos2)

      # Resumen del modelo
      summary(modelo_wls)
## 
## Call:
## lm(formula = Respuesta2 ~ Tiempo, data = Datos, weights = 1/pesos2)
## 
## Weighted Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.63812 -0.84202 -0.00815  0.89695  1.79511 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 1.006e+00  3.957e-03  254.27   <2e-16 ***
## Tiempo      3.576e-05  4.898e-05    0.73    0.471    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9889 on 28 degrees of freedom
## Multiple R-squared:  0.01868,    Adjusted R-squared:  -0.01636 
## F-statistic: 0.5331 on 1 and 28 DF,  p-value: 0.4714

En este caso, se encuentra que la regresión no es significativa, así que empleamos la desviación estándar de la pendiente para estimar la incertidumbre (ISO 35). Para esto, múltiplicamos la desviación estándar de la pendiente por el tiempo que se requería en transpote ( 72h que es lo que típicamente se demora una muestra en ser transportada), correspondería a el valor de % incertidumbre de estabilidad. A continuación muestro como hacerlo.

##  Extraer el valor p,  la pendiente y de la desviación estándar de la pendiente

      p_value <- summary(modelo_wls)$coefficients["Tiempo", "Pr(>|t|)"]
      pendiente <- coef(modelo_wls)["Tiempo"]
      desviacion_pendiente <- sqrt(vcov(modelo_wls)["Tiempo", "Tiempo"])
      
      # Mostrar los resultados
      cat(" Valor p:", p_value, "\n","Pendiente:", pendiente, "\n","Desviación estándar de la pendiente:", desviacion_pendiente, "\n")
##  Valor p: 0.4713834 
##  Pendiente: 3.576005e-05 
##  Desviación estándar de la pendiente: 4.89785e-05
# Estimar el valor de la incertidumbre
      
      # Primero se verifica el valor de  p ( si es menor  o mayor que 0.05)
      
      if (p_value < 0.05) {
        # Calcular el factor uest con la pendiente y 10 días
        uest <- abs(pendiente) * 72 #Donde 72 corresponde al número de días
      } else {
        # Calcular el factor uest con la desviación estándar de la pendiente y 10 días
        uest <- desviacion_pendiente * 72  #Donde 72 corresponde al número de días
      }
            # Mostrar el valor de uest
      cat("Factor uest:", uest, "%\n")
## Factor uest: 0.003526452 %

Este es un valor de incertidumbre esta ok, por lo cual se podría transportar a 4 ºC, obviamente, si se debe transportar a 4 ºC ( por el THC) la incertidumbre es 10 veces menor.