Aplicación de Programas Informáticos 2025

Análisis de la relación entre dos variables cuantitativas

Author
Affiliation

Eva Maure

Cátedra de Cálculo Estadístico y Biometría, FCA-UNCUYO 

library(tidyverse)
library(rstatix)
library(flextable)
library(agricolae)

Caso de estudio

Una estrategia para el control de plagas en frutales es el uso de la confusión sexual, que libera feromonas en el cultivo. Se trata de un método específico capaz de reducir efectivamente las poblaciones de insectos sin riesgos de toxicidad para el hombre y los ecosistemas. Con el objetivo de evaluar la efectividad, para condiciones ambientales y sistemas de producción locales, se desea determinar la influencia de la temperatura en la cantidad de feromonas liberadas por difusores en un monte duraznero.

Datos

library(readxl)
feromona <- read_excel("feromona.xlsx")

Análisis exploratorio

  1. Identificar las dos variables que intervienen en el estudio bivariado. 

Variable 1: Temperatura media en °C

Variable 2: Dosis media de feromona liberada

  1. Realizar una representación gráfica de los datos (diagrama de dispersión).
ggplot(data = feromona) +
  geom_point(mapping = aes(x = T_media, y = Dosis))+
  labs(x ="Temperatura media °C", y="Dosis media de feromona liberada por difusor ")+
theme_light()

  1. Describir numéricamente cada una de las variables
res<- feromona %>% 
  get_summary_stats(
    T_media, Dosis,
    show = c("mean", "sd", "min", "max")
  )
flextable(res)

variable

n

mean

sd

min

max

T_media

22

21.554

3.354

16.020

26.970

Dosis

22

0.013

0.005

0.007

0.025

¿Existe asociación entre las variables?

Coeficiente de Correlación Simple Lineal

Para el cálculo del coeficiente de correlación simple lineal muestral, se utilizan los comandos base de R:

r=cor(feromona$T_media, feromona$Dosis)
r
[1] 0.6421292

La función cor() tiene varios argumentos, entre ellos el método para el cálculo del coeficiente de CSL muestral que por defecto utiliza Pearson.

TipNota

  • recordar que r puede tomar valores entre -1 y 1

  • un valor cercano a -1 indica asociación fuerte y negativa

  • un valor cercano a 1 indica asociación furete y positiva

Note💻 Actividad

interpretar el valor de \(r\) para las variables en estudio

Prueba para el coeficiente de CSL

Hipótesis científica: “Existe correlación simple lineal entre Temperatura media (°C) y la dosis de feromona aplicada”

\(H_0 : \rho =0\)

\(H_1 : \rho \neq 0\)

\(\alpha = 0,05\)

La prueba de hipótesis para el Coeficiente de Correlación Simple Lineal \(\rho\) también sale el cálculo de la estimación:

r1_test = cor.test(feromona$T_media, feromona$Dosis)
r1_test

    Pearson's product-moment correlation

data:  feromona$T_media and feromona$Dosis
t = 3.746, df = 20, p-value = 0.001273
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.3023844 0.8371094
sample estimates:
      cor 
0.6421292

Por lo obtenido, se rechaza la hipótesis nula bajo un \(\alpha = 0,05\).

Note💻 Actividad

Completa la interpretación estádistica de la prueba de hipotesis

“Existe evidencia muestral suficiente para decir que existe asociación entre …..✏️”

¿Existe una relación funcional entre las variables?

Análisis de Regresión Simple Lineal

Sospechamos que puede existir una relación funcional entre las variables en estudio. Por ello se plantea el modelo:

\[Y_i = 𝛽_0 + 𝛽_1 x_i + 𝜀_i\]
\[i=1,...,n\] Para la estimación de los parámetros y las pruebas de hipótesis correspondientes, se utiliza:

# modelo = Y ~ X
modelo1 = lm(Dosis ~ T_media, 
             data=feromona)
modelo1

Call:
lm(formula = Dosis ~ T_media, data = feromona)

Coefficients:
(Intercept)      T_media  
 -0.0052514    0.0008679

En la salida principal del modelo tenemos las estimaciones de los parámetros. Con ellos podemos graficar la recta que ajusta a los datos:

ggplot(feromona) +
  aes(x = T_media, y = Dosis) +
  geom_smooth(method="lm", se=FALSE, color="#4A4C4F")+ 
  geom_point(colour = "#0C4C8A") +
  labs(
    x = "Temperatura media °C",
    y = "Dosis media de feromona liberada por difusor "
  ) +
  theme_light()

Para las pruebas debemos llamar a los diferentes elementos que se generaron con la función lm() en el objeto lista “modelo1”.

names(modelo1)
 [1] "coefficients"  "residuals"     "effects"       "rank"         
 [5] "fitted.values" "assign"        "qr"            "df.residual"  
 [9] "xlevels"       "call"          "terms"         "model"

Para realizar la inferencia de los coeficientes:

summary(modelo1)

Call:
lm(formula = Dosis ~ T_media, data = feromona)

Residuals:
       Min         1Q     Median         3Q        Max 
-0.0059303 -0.0022763  0.0002058  0.0020585  0.0081383 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)   
(Intercept) -0.0052514  0.0050509  -1.040  0.31089   
T_media      0.0008679  0.0002317   3.746  0.00127 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.003561 on 20 degrees of freedom
Multiple R-squared:  0.4123,    Adjusted R-squared:  0.3829 
F-statistic: 14.03 on 1 and 20 DF,  p-value: 0.001273

Para la del modelo, es decir teniendo en cuenta la partición de la variabilidad del modelo:

anova(modelo1)
Analysis of Variance Table

Response: Dosis
          Df     Sum Sq    Mean Sq F value   Pr(>F)   
T_media    1 0.00017790 1.7790e-04  14.033 0.001273 **
Residuals 20 0.00025355 1.2678e-05                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Supuestos

  • Independencia entre observaciones.

  • Homocedasticidad.

  • Normalidad de residuos.

Q-Q Plot

Este gráfico muestra cómo se acumulan los residuos respecto de los cuantiles teóricos de una distribución normal. Si la distribución de residuos es normal, los veremos cercanos a la recta. Desviaciones de la recta indican que la distribución de los residuos no es normal.

# 1. Extraer los residuos y meterlos en un dataframe simple
datos_residuos <- data.frame(mis_residuos = residuals(modelo1))

# 2. Graficar leyendo directamente ese dataframe
ggplot(datos_residuos, aes(sample = mis_residuos)) +
  stat_qq() +
  stat_qq_line() +
  xlab("Cuartiles teóricos (normales)") +
  ylab("Cuartiles de los residuos")

Residuos estandarizados vs Predichos

Este gráfico es similar al primero, pero los residuos están estandarizados. Esperamos no ver ningún patrón, con los residuos distribuidos normalmente alrededor de cero.

ggplot(data=modelo1)+
  geom_point(mapping=aes(x=modelo1$fitted.values, y=rstandard(modelo1)))+
  geom_hline(yintercept = c(-3,0,3),col="red")+
  xlab("Predichos")+
  ylab("Residuos estandarizados")
Warning: `fortify(<lm>)` was deprecated in ggplot2 4.0.0.
ℹ Please use `broom::augment(<lm>)` instead.
ℹ The deprecated feature was likely used in the ggplot2 package.
  Please report the issue at <https://github.com/tidyverse/ggplot2/issues>.

Ajustamiento de modelos intrínsecamente lineales

Seguims usando la funicón lm() pero moficiacmos los argumentos según el modelo a evaluar.

Modelo polinómico de 2° grado

Teniendo en cuenta que el modelo es. \[Y_i = \beta _0 + \beta _1 x_i + \beta _2 x_i^2 + \varepsilon _i\]

\[i=1,...,n\] La función se escribe:

modelo2 = lm(feromona$Dosis~feromona$T_media + I(feromona$T_media^2))
summary(modelo2)

Call:
lm(formula = feromona$Dosis ~ feromona$T_media + I(feromona$T_media^2))

Residuals:
       Min         1Q     Median         3Q        Max 
-0.0054085 -0.0019800 -0.0003171  0.0017204  0.0077466 

Coefficients:
                        Estimate Std. Error t value Pr(>|t|)
(Intercept)            3.209e-02  3.380e-02   0.950    0.354
feromona$T_media      -2.681e-03  3.184e-03  -0.842    0.410
I(feromona$T_media^2)  8.236e-05  7.371e-05   1.117    0.278

Residual standard error: 0.003539 on 19 degrees of freedom
Multiple R-squared:  0.4486,    Adjusted R-squared:  0.3905 
F-statistic: 7.728 on 2 and 19 DF,  p-value: 0.003501

Modelo potencial

Teniendo en cuenta que el modelo es. \[Y_i = \beta _0 * x_i^{\beta_1}*\varepsilon _i\]

\[i=1,...,n\] Aplicamos logarirmo en ambos miembros del modelo y obtenemos la linealización del modelo exponencial: \[ln(Y_i) = ln(\beta _0) +ln(x_i)\beta_1+ln(\varepsilon _i)\] \[Y_i^* = \beta _0^* +\beta_1 x_i^*+\varepsilon _i^*\]

Se escribe:

modelo3 = lm(log(feromona$Dosis)~log(feromona$T_media))
summary(modelo3)

Call:
lm(formula = log(feromona$Dosis) ~ log(feromona$T_media))

Residuals:
     Min       1Q   Median       3Q      Max 
-0.57977 -0.13483  0.05252  0.17481  0.44511 

Coefficients:
                      Estimate Std. Error t value Pr(>|t|)    
(Intercept)            -8.2375     1.1415  -7.216 5.52e-07 ***
log(feromona$T_media)   1.2673     0.3727   3.400  0.00284 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.2707 on 20 degrees of freedom
Multiple R-squared:  0.3663,    Adjusted R-squared:  0.3346 
F-statistic: 11.56 on 1 and 20 DF,  p-value: 0.002841

Modelo exponencial

Teniendo en cuenta que el modelo es. \[Y_i = \beta _0 * e^{\beta_1x_i}*\varepsilon _i\]

\[i=1,...,n\] Aplicamos logarirmo en ambos miembros del modelo y obtenemos la linealización del modelo exponencial: \[ln(Y_i) = ln(\beta _0) +\beta_1x_i+ln(\varepsilon _i)\] \[Y_i^* = \beta _0^* +\beta_1x_i+\varepsilon _i^*\]Se escribe:

modelo4 = lm(log(feromona$Dosis)~feromona$T_media)
summary(modelo4)

Call:
lm(formula = log(feromona$Dosis) ~ feromona$T_media)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.56353 -0.16086  0.05162  0.18075  0.42987 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)      -5.69125    0.37701 -15.096 2.14e-12 ***
feromona$T_media  0.06172    0.01729   3.569  0.00192 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.2658 on 20 degrees of freedom
Multiple R-squared:  0.3891,    Adjusted R-squared:  0.3585 
F-statistic: 12.74 on 1 and 20 DF,  p-value: 0.001923
Note💻 Actividad

Considerando todos los análisis realizados anteriormente ¿se puede decir que el modelo se ajusta a los datos?