Ejercicio 1

Se presenta una base de datos de mediciones de NDVI realizadas en las cuatro estaciones del año 2020. En cada época del año, se registró el valor del índice en diferentes sitios dentro de la región Pampeana clasificados como cobertura agrícola (ndvi_4_estaciones.txt), por lo que cada observación puede considerarse independiente. Se desea realizar un ANAVA teniendo en cuenta el efecto de la época.

1.1.1. Actividades

  1. Construya tablas y/o gráficos para realizar una estadística descriptiva de la variable NDVI en función de la época del año.

  2. Escribir un modelo para estudiar el efecto de la estación sobre el NDVI.

  3. Controle si se cumplen los supuestos del análisis ajustado.

  4. ¿Existen diferencias de NDVI entre las épocas del año? Realizar una comparación de medias y concluir.


Respuestas:
1) Construya tablas y/o gráficos para realizar una estadística descriptiva de la variable NDVI en función de la época del año.

1.1.1 - Medidas de posicion de la variable (NDVI)

kable(ndvi_estaciones %>%
        summarise(
          n = n(),
          min = min(ndvi),
          max = max(ndvi),
          media = mean(ndvi),
          mediana = median(ndvi),
        ))
n min max media mediana
112 -0.0472835 0.8996883 0.3181682 0.2756814

1.1.2 - Medidas de dispersión de la variable (NDVI)

kable(ndvi_estaciones %>%
        summarise(
          n = n(),
          rango = max(ndvi) - min(ndvi),
          var = var(ndvi),
          de = sd(ndvi),
          cv = sd(ndvi) / mean(ndvi) * 100,
          EE = de / sqrt(n()),
          iqr = IQR(ndvi)
        ))
n rango var de cv EE iqr
112 0.9469717 0.0382696 0.1956261 61.48514 0.0184849 0.2389991

1.1.3 - Percentiles de la variable (NDVI) - (Son también un tipo de medidas de posición)

kable(ndvi_estaciones %>%
        summarise(
          n = n(),
          p05 = quantile(ndvi, 0.05),
          p25 = quantile(ndvi, 0.25),
          p50 = quantile(ndvi, 0.50),
          p75 = quantile(ndvi, 0.75),
          p95 = quantile(ndvi, 0.95)
        ))
n p05 p25 p50 p75 p95
112 0.0876447 0.1806242 0.2756814 0.4196233 0.6575365

El “n” de la muestra no es una medida de posición ni dispersión, pero es información muy útil por lo tanto la agregué en todas las tablas.

1.1.4 - Gráfico de la variable (NDVI)

ggplot(ndvi_estaciones, aes(y = ndvi)) +
  geom_boxplot() +
  ggtitle("Gráfico Boxplot") + 
  ylab("NDVI") +
  theme(plot.title = element_text(hjust = 0.5, size = 14))

1.2.1 - Medidas de posicion de la variable (NDVI) particionadas por los niveles del factor (Estación)

kable(ndvi_estaciones %>%
        group_by(estacion) %>%
        summarise(
          n = n(),
          min = min(ndvi),
          max = max(ndvi),
          media = mean(ndvi),
          mediana = median(ndvi),
        ))
estacion n min max media mediana
invierno 28 0.0384022 0.3492455 0.2033106 0.2044973
otono 28 -0.0472835 0.3168594 0.1721308 0.1646397
primavera 28 0.1291092 0.4391698 0.2891186 0.2968141
verano 28 0.3845793 0.8996883 0.6081126 0.6093405

1.2.2 - Medidas de dispersión de la variable (NDVI) particionadas por los niveles del factor (Estación)

kable(ndvi_estaciones %>%
        group_by(estacion) %>%
       summarise(
          n = n(),
          rango = max(ndvi) - min(ndvi),
          var = var(ndvi),
          de = sd(ndvi),
          cv = sd(ndvi) / mean(ndvi) * 100,
          EE = de / sqrt(n()),
          iqr = IQR(ndvi)
        ))
estacion n rango var de cv EE iqr
invierno 28 0.3108433 0.0081843 0.0904670 44.49692 0.0170966 0.1597686
otono 28 0.3641429 0.0084200 0.0917604 53.30852 0.0173411 0.1276084
primavera 28 0.3100606 0.0061769 0.0785932 27.18373 0.0148527 0.0755484
verano 28 0.5151089 0.0106952 0.1034177 17.00634 0.0195441 0.1218786

1.2.3 -Percentiles de la variable (NDVI) particionadas por los niveles del factor (Estación)

kable(ndvi_estaciones %>%
 group_by(estacion) %>%
        summarise(
          n = n(),
          p05 = quantile(ndvi, 0.05),
          p25 = quantile(ndvi, 0.25),
          p50 = quantile(ndvi, 0.50),
          p75 = quantile(ndvi, 0.75),
          p95 = quantile(ndvi, 0.95)
        ))
estacion n p05 p25 p50 p75 p95
invierno 28 0.0793948 0.1113486 0.2044973 0.2711172 0.3452533
otono 28 0.0278970 0.1161130 0.1646397 0.2437214 0.3023113
primavera 28 0.1541201 0.2449843 0.2968141 0.3205327 0.4130891
verano 28 0.4625638 0.5314781 0.6093405 0.6533566 0.7653440

1.2.4 Gráfico de Cajas

ggplot(ndvi_estaciones, aes(y = ndvi, x = estacion, fill = estacion)) +
  geom_boxplot() +
  ggtitle("Gráfico de cajas por estación") + 
  ylab("NDVI") +
  xlab("Estación") +
  theme(plot.title = element_text(hjust = 0.5, size = 14))

1.2.5 Gráfico de Frecuencias Acumuladas

ggplot(ndvi_estaciones, aes(x = ndvi, color = estacion)) +
  stat_ecdf() +
  ylab("Frecuencia relativa acumulada") +
  xlab("NDVI") +
  ggtitle("Frecuencia relativa acumulada - NDVI por estación") +
  theme(plot.title = element_text(hjust = 0.5, size = 14))

Comentarios

A partir de las tablas y gráficos obtenidos se puede señalar que la media muestral de NDVI más alta sucede en verano y la más baja en otoño. La mayor variabilidad (varianza) sucede en verano mientras que la menor sucede en primavera.


2) Escribir un modelo para estudiar el efecto de la estación sobre el NDVI.

Primero definimos el rol de cada variable que interviene en el modelo.
Variable Respuesta (VR): NVDI
Factor: Estación

Ahora si, definimos el modelo:
\(y_{ij} = \mu + \tau_{i} + e_{ij}\)

En donde…
\(i\): es el nivel del factor.
\(j\): es la repetición dentro del nivel del factor \(i\).
\(y\): es la variable respuesta.
\(y_{ij}\): es la observación \(j\) del nivel del factor \(i\) de la variable respuesta.
\(\mu\): es el parametro poblacional de la variable respuesta sin discriminar los niveles del factor.
\(\tau_{i}\): es el corrimiento de la media del nivel del factor \(i\) respecto a \(\mu\).
\(e_{ij}\): es el error asociado a la observación \(y_{ij}\).

Es importante resaltar que: \(\mu + \tau_{i} = \mu_{i}\)
\(\tau_{i}\): es la media poblacional del nivel del factor.
El modelo entonces también puede entenderse como:
\(y_{ij} = \mu_{i} + e_{ij}\)

El modelo es con parámetros poblacionales. Nosotros haremos las estimaciones de cada término a partir de estimadores muestrales. En tal caso, para los datos del ejercicios las estimaciones son:

El estimador de la media poblacional de cada nivel del factor es la media muestral de cada nivel del factor:

# Medias muestrales de cada nivel del factor
medias_nivel <- tapply(ndvi_estaciones$ndvi, ndvi_estaciones$estacion, mean)
medias_nivel
##  invierno     otono primavera    verano 
## 0.2033106 0.1721308 0.2891186 0.6081126

El estimador de la media poblacional, es la media de las medias muestrales:

# Media de las medias, muestral
media_de_las_medias <- mean(medias_nivel)
media_de_las_medias
## [1] 0.3181682

Un error común el usar como estimador de la media poblacional a la media de todos los datos de la muestra. En el modelo ANOVA La media de todos los datos es igual a la media de las medias solo cuando el “n” de todos los niveles del factor es el mismo, pero en cualquier otro caso obtendremos valores diferentes. Siempre hay que usar la media de las medias como estimador de \(\mu\).

El estimador de los efectos Tau, es la diferencia entre la media muestral de cada nivel del factor respecto a la media de las medias:

# Efectos de corrimiento para cada nivel del factor
corrimiento_nivel <- medias_nivel - media_de_las_medias
corrimiento_nivel
##    invierno       otono   primavera      verano 
## -0.11485757 -0.14603731 -0.02904958  0.28994446

Aquí también pueden verse algunos detalles que están definidos por demostración matemática, y es que la suma de todos los efectos de corrimiento es cero. Para el cálculo informático podríamos obtener resultado cero o muy valor muy muy chico.

# Suma efectos de corrimiento para cada nivel del factor(debe dar cero o muy muy chico)
sum(corrimiento_nivel)
## [1] 0

ANOVA

Juego de Hipótesis
H0: Todas las medias de NDVI por estación son iguales a nivel α=0.05.
H1: Por lo menos dos de las medias de NDVI por estación son diferentes a nivel α=0.05.

fit_ndvi <- aov(ndvi ~ estacion, data = ndvi_estaciones)
summary(fit_ndvi)
##              Df Sum Sq Mean Sq F value Pr(>F)    
## estacion      3  3.344  1.1147   133.2 <2e-16 ***
## Residuals   108  0.904  0.0084                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Los valores obtenidos de la prueba de ANOVA se deben validar mediante la verificación de los requisitos del modelo. Los requisitos pueden expresarse matemáticamente como:
\(e_{ij} \sim N,I,(0,\sigma)\)
Esto quiere decir que los errores del modelo ANOVA sean independientes, y tengan distribución normal con media cero y homogeneidad de varianzas. 


3) Controle si se cumplen los supuestos del análisis ajustado.

Requisito de independencia

El diseño de la expeciencia debe garantizar la independencia. El enunciado indica que las observaciones pueden considerarse como independientes. No tenemos por que dudar de la información proporcionada. Manifestamos que los datos son independientes, y de esa manera lo son los residuos del modelo también.


Requisito de Normalidad de los residuos

Se puede hacer de manera formal e informal. La manera formal es con un test estadístico que otorgue un valor p, en tal caso la herramienta propuesta es el test de Shapiro-Wilk. La prueba informal es un gráfico denominado QQplot que se realiza con los residuos estandarizados respecto a los cuantiles de una distribución normal y se los compara con una recta y=x para comprobar la normalidad.


Test de Shapiro-Wilks

H0: Los residuos estandarizados tienen una distribución normal a nivel α=0.05.

H1: Los residuos estandarizados no tienen una distribución normal a nivel α=0.05.

shapiro.test(residuals(fit_ndvi))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(fit_ndvi)
## W = 0.99132, p-value = 0.7025

El valor p es 0.7025. El valor alfa es 0.05.
El valor p es mayor que el valor de alfa, por lo tanto no se rechaza Ho.
Los residuos presentan una distribución estadísticamente normal.


Gráfico QQplot
qqnorm(rstandard(fit_ndvi)); qqline(rstandard(fit_ndvi))
qqline(rstandard(fit_ndvi), col = "orange", lwd = 3)

Se realiza un gráfico de los residuos estandarizados vs valores predichos.
A simple vista, a mi ojímetro, los residuos se ajustan bien a lo esperado para una distribución normal.
Lo correcto sería tomar decisión con el test de Shapiro-Wilk ó con el QQplot, y no mirar ambos y elegir el que a uno le conviene.



Requisito media cero

Esta característica se da si o si. Hay un teorema que demuestra que la suma de los residuos es igual a cero. Por lo tanto no eay un test estadístico que realizar. Este punto se cumple si o si. En términos informáticos, la suma de los errores dará cero o un valor muy muy chico.

# La suma de los residuos debe dar cero o un valor muy muy pequeño
sum(residuals(fit_ndvi))
## [1] -5.70724e-16



Requisito de Homogeneidad de varianzas de los residuos

Se puede hacer de manera formal e informal. La manera formal es con un test estadístico que otorgue un valor p, en tal caso la herramienta propuesta en el curso fue el test de Levenne que utiliza los residuos absolutos. La prueba informal es un gráfico de Residuos estandarizados vs predichos del modelo.


Prueba de Levene

H0: Las varianzas de todos los niveles del factor son homogéneas.
H1: Al menos dos varianzas de los niveles del factor son heterogéneas.

summary(aov(abs(residuals(fit_ndvi))~ndvi_estaciones$estacion))
##                           Df  Sum Sq  Mean Sq F value Pr(>F)
## ndvi_estaciones$estacion   3 0.00496 0.001652    0.57  0.636
## Residuals                108 0.31312 0.002899

El valor p es 0.636. El valor alfa es 0.05. El valor p es mayor que el valor alfa, por lo tanto no se rechaza Ho de Levenne. Las varianzas de los niveles del factor son estadísticamente homogéneas.


Gráfico de Residuos Estandarizados vs Predichos del Modelo
# Generaremos el gráfico con el mismo rango de valores en el eje Y, centrado en 0.
minimo <- min(rstandard(fit_ndvi))
maximo <- max(rstandard(fit_ndvi))
minimo_absoluto <- min(abs(minimo), abs(maximo))
maximo_absoluto <- max(abs(minimo), abs(maximo))
maximo_total <- max(minimo_absoluto, maximo_absoluto)

# Gráfico de residuos estudentizados vs predichos del modelo
plot(fitted.values(fit_ndvi), rstandard(fit_ndvi),
     xlab = "Valores predichos", ylab = "Residuos estandarizados", main = "Gráfico de Valores Predichos vs Residuos Estandarizados", ylim = c(-maximo_total, maximo_total))

# Agregamos una linea horizontal en cero.
abline(h=0, col = "blue", lwd = 3)

A mi ojímetro… Me parece que ocurre la homogeneidad de varianzas de los residuos.
Nuevamente, esta decisión hay que tomarla sin mirar el test y el gráfico y entonces quejarse con la que a uno le viene bien. Hay que decidir con qué criterio se tomará la decisión sin ver nuestros datos en ellos.


Visión Global de los requisitos de ANOVA

Se cumplen simultáneamente todos los requisitos del test de ANOVA a 1 Factor, por lo tanto es válido sacar conclusiones del valor p de la tabla de ANOVA.


4) ¿Existen diferencias de NDVI entre las épocas del año? Realizar una comparación de medias y concluir.
Hemos ya definido previamente que es válido interpretar y sacar conclusiones del valor p de ANOVA.
Para el test de ANOVA, el valor p <2e-16. El valor alfa es 0.05. El valor p es menor que el valor de alfa, por lo tanto se rechaza Ho. Existen diferencias estadísticamente significativas en al menos dos medias del factor estación.

Método de Tukey para verificar qué medias son significativamente diferentes

TukeyHSD(fit_ndvi, conf.level = 0.95)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = ndvi ~ estacion, data = ndvi_estaciones)
## 
## $estacion
##                           diff         lwr       upr     p adj
## otono-invierno     -0.03117974 -0.09498108 0.0326216 0.5806052
## primavera-invierno  0.08580799  0.02200665 0.1496093 0.0036165
## verano-invierno     0.40480202  0.34100068 0.4686034 0.0000000
## primavera-otono     0.11698773  0.05318639 0.1807891 0.0000321
## verano-otono        0.43598176  0.37218042 0.4997831 0.0000000
## verano-primavera    0.31899404  0.25519270 0.3827954 0.0000000
plot(TukeyHSD(fit_ndvi, conf.level = 0.95), las = 1)

Si bien el eje Y no se aprecia bien, la tabla anterior presenta el mismo orden de combinaciones de estación que el eje Y de este gráfico. Las diferencias de las medias de otoño e invierno incluye a cero, el resto de las diferencias de medias no incluye a cero. Esto indica que otoño e invierno tienen medias estadísticamente iguales entre si y diferentes al resto de las estaciones. Sin embargo, aún falta contextualizar, ya que sabemos que son iguales entre si y diferentes al resto, pero no sabemos si son iguales y superiores al resto o iguales e inferiores. Este detalle será presentado más adelante en otra tabla.

plot(HSD.test(fit_ndvi, trt = "estacion", alpha = 0.05), main = "Medias de NDVI agrupados por estación", las = 1)

Agrego a continuación la tabla con la cual es generado el gráfico anterior. Allí podrá visualizarse el valor de cada media de cada nivel del factor directamente, y también ver el grupo estadístico asignado por Tukey.

TUKEY_COMPLETO <- HSD.test(fit_ndvi, "estacion", alpha = 0.05)


# Tabla de Tukey
TABLA_TUKEY <- TUKEY_COMPLETO$groups
TABLA_TUKEY
##                ndvi groups
## verano    0.6081126      a
## primavera 0.2891186      b
## invierno  0.2033106      c
## otono     0.1721308      c

Conclusión final ejercicio 1

Ahora si, entre este último gráfico y tabla, queda de manifiesto que las medias de invierno y otoño son estadísticamente iguales entre si, y que son las menores medias. En el contexto de NDVI, los menores valores de NDVI se observan en otoño e invierno (sin poder discriminar entre ambos). El valor NDVI se observa en verano, y un nivel intermedio de NDVI se observa en primavera.  Tanto en el gráfico como en la tabla, las letras que se asignan a los niveles del factor están en realación a su semejanza en términos estadísticos. Letras iguales implican niveles del factor estadísticamente iguales. Letras diferentes implican niveles del factor estadísticamente diferentes.


Ejercicio 2

La base de datos (evaporacion.txt) contiene datos de la evaporación diaria de suelo (agua_eva), de la temperatura máxima del suelo (tpsuelma) y del aire (tairemax) registrada en ese día (en fahrenheit). El objetivo general es encontrar un modelo funcional para predecir la evaporación diaria.

1.2.1. Actividades

  1. Realizar una estadística descriptiva de la base de datos, seleccionando las medidas resumen y/o gráficos que prefiera.

  2. Plantee un modelo especificando como regresoras a la temperatura máxima del suelo y la del aire.

  3. ¿Considera que utilizar ambas variables en el modelo podría ser adecuado? ¿Por qué?

  4. Para cada una de las potenciales variables regresoras, ajuste un nuevo modelo de regresión lineal simple.

  5. ¿Que proporción de la variabilidad total es explicada por cada modelo?

  6. ¿Cuál modelo considera más adecuado?

  7. Especifique la ecuación ajustada para el modelo que considere más adecuado.

  8. Prediga cuánto se espera que se evapore si la temperatura máxima del suelo es de 88°F.




Respuestas:

1) Realizar una estadística descriptiva de la base de datos, seleccionando las medidas resumen y/o gráficos que prefiera.

1- Tabla Resumen

1.1.1 - Medidas de posicion de las 3 variables.

kable(evaporacion %>%
        summarise(
          variables = colnames(evaporacion),
          n = c(n(),n(), n()),
          min = c(min(tpsuelma), min(tairemax), min(agua_eva)),
          max = c(max(tpsuelma), max(tairemax), max(agua_eva)),
          media = c(mean(tpsuelma), mean(tairemax), mean(agua_eva)),
          mediana = c(median(tpsuelma), median(tairemax), median(agua_eva))
        ))
variables n min max media mediana
tpsuelma 25 73 93 83.92 84
tairemax 25 77 94 87.72 88
agua_eva 25 4 54 31.28 33

1.1.2 - Medidas de dispersión de las 3 variables.

kable(evaporacion %>%
        summarise(
          variables = colnames(evaporacion),
          n = c(n(),n(), n()),
          rango = c((max(tpsuelma) - min(tpsuelma)), (max(tairemax) - min(tairemax)), (max(agua_eva) - min(agua_eva))),
          var = c(var(tpsuelma), var(tairemax), var(agua_eva)),
          de = c(sd(tpsuelma), sd(tairemax), sd(agua_eva)),
          cv = c((sd(tpsuelma) / mean(tpsuelma) * 100), (sd(tairemax) / mean((tairemax) * 100)), (sd(agua_eva) / mean(agua_eva) * 100)),
          EE = de / sqrt(n()),
          iqr = c(IQR(tpsuelma), IQR(tairemax), IQR(agua_eva))
        ))
variables n rango var de cv EE iqr
tpsuelma 25 20 30.32667 5.506965 6.5621607 1.1013931 7
tairemax 25 17 23.62667 4.860727 0.0005541 0.9721454 8
agua_eva 25 50 198.29333 14.081667 45.0181156 2.8163333 20

1.1.2 - Percentiles de las 3 variables. - (Son también un tipo de medidas de posición)

kable(evaporacion %>%
        summarise(
          variables = colnames(evaporacion),
          n = c(n(),n(), n()),
          p05 = c(quantile(tpsuelma, 0.05), quantile(tairemax, 0.05), quantile(agua_eva, 0.05)),
          p25 = c(quantile(tpsuelma, 0.25), quantile(tairemax, 0.25), quantile(agua_eva, 0.25)),
          p50 = c(quantile(tpsuelma, 0.50), quantile(tairemax, 0.50), quantile(agua_eva, 0.50)),
          p75 = c(quantile(tpsuelma, 0.75), quantile(tairemax, 0.75), quantile(agua_eva, 0.75)),
          p95 = c(quantile(tpsuelma, 0.95), quantile(tairemax, 0.95), quantile(agua_eva, 0.95))
        ))
variables n p05 p25 p50 p75 p95
tpsuelma 25 74.2 81 84 88 92.4
tairemax 25 79.0 84 88 92 94.0
agua_eva 25 6.0 23 33 43 49.4

Gráfico de puntos T° Máxima de suelo vs Agua Evaporada

ggplot(evaporacion, aes(tpsuelma, agua_eva)) +
  geom_point() +
  ylab("Agua evaporada") +
  xlab("T° Máxima del suelo") +
  ggtitle("T° Máxima del suelo vs Agua Evaporada") + 
  theme(plot.title = element_text(hjust = 0.5, size = 14))

Gráfico de puntos T° Máxima del Aire vs Agua Evaporada

ggplot(evaporacion, aes(tairemax, agua_eva)) +
  geom_point() +
  ylab("Agua evaporada") +
  xlab("T° Máxima del aire") +
  ggtitle("T° Máxima del aire vs Agua Evaporada") + 
  theme(plot.title = element_text(hjust = 0.5, size = 14))

Gráfico 3D

plot3d( 
  x = evaporacion$tpsuelma, 
  y = evaporacion$tairemax,
  z = evaporacion$agua_eva,
  type = 's', 
  radius = 1,
  xlab = "Temp Max Suelo", 
  ylab = "Temp Max Aire", 
  zlab = "Agua Evaporada",
  col = "blue")
rglwidget(elementId = "plot3drgl")

Es posible interactuar con este gráfico. Se puede rotar y amplicar con el mouse.


El detalle importante del gráfico 3D es que el software al generar si o si ejes cartesianos ortogonales ha asumido que las variables regresoras no están correlacionadas. Este gráfico es útil, es lindo, pero puede ser una linda mentira 3D. Más adelante veremos si las regresoras están correlacionadas o no.


Comentarios de los gráficos XY y 3D.

Gráficamente se puede observar que hay una relación positiva entre las variables regresoras y la variable respuesta.

2. Plantee un modelo especificando como regresoras a la temperatura máxima del suelo y la del aire.

Variable Respuesta o dependiente (Y): Agua evaporada.

Variables Independientes o regresoras 1 (X1): T° Máxima de suelo.

Variables Independientes o regresoras 2 (X2): T° Máxima del aire.

Modelo de regresión lineal que considera temperatura del suelo y del aire como variables explicadoras.

rlm_eva <- lm(agua_eva ~ tpsuelma + tairemax, data = evaporacion)
summary(rlm_eva)
## 
## Call:
## lm(formula = agua_eva ~ tpsuelma + tairemax, data = evaporacion)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -14.4001  -2.1488   0.0814   3.3825  13.2875 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -169.4393    24.7330  -6.851    7e-07 ***
## tpsuelma       1.8969     0.6466   2.934  0.00768 ** 
## tairemax       0.4734     0.7325   0.646  0.52479    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.549 on 22 degrees of freedom
## Multiple R-squared:  0.8018, Adjusted R-squared:  0.7837 
## F-statistic: 44.49 on 2 and 22 DF,  p-value: 1.859e-08


El modelo ha sido planteado según lo requerido en el orden de los items del práctico. A esta altura, si bien el modelo está pantleado, no está definido que sea válido sacar conclusiones de ese modelo ya que no han sido puestos a preuba aún los requisitos.


3. ¿Considera que utilizar ambas variables en el modelo podría ser adecuado? ¿Por qué?

Para considerar la utilización de ambas variables, hay toda una serie de requisitos que deben cumplirse para tomar la decisión. Hay que corroborar primeramente la no correlación entre las variables regresoras, que se lleva a cabo con un test de correlación de Pearson. La propia correlación de Pearson tiene requisitos de normalidad de ambas variables ye de homogeneidad de varianzas.

Análisis de normalidad de la variable de temperatura del suelo:

normalidad_x1 <- shapiro.test(evaporacion$tpsuelma)
normalidad_x1
## 
##  Shapiro-Wilk normality test
## 
## data:  evaporacion$tpsuelma
## W = 0.9544, p-value = 0.3143

El valor p es 0.3143. El valor alfa es 0.05.
El valor p es mayor que el valor alfa. No se rechaza Ho.
La variable “tpsueloma” tiene una distribucion estadisticamente normal.


Análisis de normalidad de la variable de temperatura del aire:

normalidad_x2 <- shapiro.test(evaporacion$tairemax)
normalidad_x2
## 
##  Shapiro-Wilk normality test
## 
## data:  evaporacion$tairemax
## W = 0.93164, p-value = 0.09481

El valor p es 0.09481. El valor alfa es 0.05.
El valor p es mayor que el valor alfa. No se rechaza Ho.
La variable “tairemax” tiene una distribucion estadisticamente normal.


Análisis de homogeneidad de varianzas de Bartlett entre las variables X1 y X2:

homogeneidad_regresoras <- bartlett.test(list(evaporacion$tpsuelma, evaporacion$tairemax))
homogeneidad_regresoras
## 
##  Bartlett test of homogeneity of variances
## 
## data:  list(evaporacion$tpsuelma, evaporacion$tairemax)
## Bartlett's K-squared = 0.36538, df = 1, p-value = 0.5455

El valor p es 0.5455. El valor alfa es 0.05.
El valor p es mayor que el valor alfa. No se rechaza Ho.
La variables “tpsuelma” y “tairemax” son estadisticamente homogeneas.


Ahora si podemos realizar el test de correlación de Pearson, ya que cumplimos sus requisitos: # Ahora si, es valido aplicar el test de correlacion de Pearson.

correlacion_pearson <- cor.test(evaporacion$tpsuelma, evaporacion$tairemax)
correlacion_pearson
## 
##  Pearson's product-moment correlation
## 
## data:  evaporacion$tpsuelma and evaporacion$tairemax
## t = 11.84, df = 23, p-value = 2.886e-11
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.8389909 0.9676179
## sample estimates:
##      cor 
## 0.926858

El valor p es 2.886e-11. El valor alfa es 0.05.
El valor p es menor que el valor de alfa. Se rechaza Ho.
La variables “tpsuelma” y “tairemax” están estadísticamente correlacionadas.


Esto implica que no es valida la generacion de un modelo de regresion lineal que posea simultaneamente ambas variables regresoras.
El modelo de regresión lineal que incluye ambas variables regresoras en el modelo, no es válido para sacar conclusiones. No es factible hacer ningún tipo de interpretación de los valores p obtenidos de esa regresión.


Los caminos que quedan son de generar dos modelos de regrsión lineal simple, uno para cada variable regresoras, y corroborar los requisitos de cada modelo. Si ambos modelos son válidos, se define la preferencia por la utilización de uno u otro por aquel modelo que tenga el mayor valor de \(R^2\).

4. Para cada una de las potenciales variables regresoras, ajuste un nuevo modelo de regresión lineal simple.

Modelo 1

Modelo de regresión lineal que considera temperatura del suelo como variable X:

# Modelo 1: Aporacion de Agua vs. Temperatura de suelo
modelo01 <- lm(agua_eva ~ tpsuelma, data = evaporacion)

reg_lineal01_completa <- summary(modelo01)

residuos01 <- modelo01$residuals
predichos01 <- modelo01$fitted.values

modelo01
## 
## Call:
## lm(formula = agua_eva ~ tpsuelma, data = evaporacion)
## 
## Coefficients:
## (Intercept)     tpsuelma  
##    -160.413        2.284



Test de Normalidad de los residuos del modelo 1

normalidad01 <- shapiro.test(residuos01)
normalidad01
## 
##  Shapiro-Wilk normality test
## 
## data:  residuos01
## W = 0.98013, p-value = 0.8876

El valor p del test de normalidad de shapiro-wilk es 0.8876.
El valor alfa es 0.05.
El valor p es mayor que el valor alfa, por lo tanto no se rechaza Ho.
Los residuos del modelo 1 presentan distribucion estadisticamente normal.


Homogeneidad de varianzas de los residuos del modelo 1.

# Grafico de dispersion para evaluar la homogeneidad de los errores del modelo 1
plot(predichos01, rstandard(modelo01), main = "Modelo 1",
     xlab = "Valores predichos", ylab = "Residuos estandarizados")
abline(h=0, col="blue")

Visualizando el grafico de residuos estandarizados y predichos, a mi ojo, se cumple el requisito de homogeneidad de los errores.

Decisión general sobre el modelo 1

Se cumplen ambos requisitos del modelo, por lo tanto el modelo 1 es valido.


Modelo 2

Modelo de regresión lineal que considera temperatura del aire como variable X:

# Modelo 2: Aporacion de Agua vs. Temperatura de aire
modelo02 <- lm(agua_eva ~ tairemax, data = evaporacion)

reg_lineal02_completa <- summary(modelo02)

residuos02 <- modelo02$residuals
predichos02 <- modelo02$fitted.values

modelo02
## 
## Call:
## lm(formula = agua_eva ~ tairemax, data = evaporacion)
## 
## Coefficients:
## (Intercept)     tairemax  
##    -184.982        2.465



Test de Normalidad de los residuos del modelo 2

normalidad02 <- shapiro.test(residuos02)
normalidad02
## 
##  Shapiro-Wilk normality test
## 
## data:  residuos02
## W = 0.97364, p-value = 0.7377

El valor p del test de normalidad de shapiro-wilk es 0.7377.
El valor alfa es 0.05.
El valor p es mayor que el valor alfa, por lo tanto no se rechaza Ho.
Los residuos del modelo 2 presentan distribucion estadisticamente normal.


Homogeneidad de varianzas de los residuos del modelo 2.

# Grafico de dispersion para evaluar la homogeneidad de los errores del modelo 2
maximo_general <- max(abs(rstandard(modelo02)))

plot(predichos02, rstandard(modelo02), main = "Modelo 2",
     xlab = "Valores predichos", ylab = "Residuos estandarizados",
     ylim = c(-maximo_general, maximo_general))
abline(h=0, col="blue")


Visualizando el grafico de residuos estandarizados y predichos, a mi ojo, se cumple el requisito de homogeneidad de los errores del modelo 2.


Decisión general sobre el modelo 2.

Se cumplen ambos requisitos del modelo, por lo tanto el modelo 2 es valido.

5. ¿Que proporción de la variabilidad total es explicada por cada modelo? Para el modelo 1, el \(R^2\) es:

r1 <- reg_lineal01_completa$r.squared
r1
## [1] 0.7979935

El modelo 1 explica un 0.7979935 de la variabilidad total.
Podemos redondearlo a 0.80 para el modelo 1.


Para el modelo 2, el \(R^2\) es:

r2 <- reg_lineal02_completa$r.squared
r2
## [1] 0.7241965

El modelo 1 explica un 0.7241965 de la variabilidad total.
Podemos redondearlo a 0.72 para el modelo 2.


6. ¿Cuál modelo considera más adecuado? Los criterios son:
1) Modelos que cumplen todos los requisitos
2) Aquellos que cumplen el punto1, aquel cuyo \(R^2\) sea mayor.

En tal caso, nos quedamos con el modelo 1, que utilizó como variable regresona a la temperatura del suelo.


7. Especifique la ecuación ajustada para el modelo que considere más adecuado. Para el modelo 1, la ecuación de la recta es:
\(y = b_0 + b_1*x + e_{ij}\)

Las estimaciones de ordenada y pendiente que usaremos son las obtenidas del modelo 1:

# Modelo 1: Aporacion de Agua vs. Temperatura de aire
modelo01
## 
## Call:
## lm(formula = agua_eva ~ tpsuelma, data = evaporacion)
## 
## Coefficients:
## (Intercept)     tpsuelma  
##    -160.413        2.284

Y con esa información nuestra recta es:

\(y = -160.413 +2.284*x + e_{ij}\)

8. Prediga cuánto se espera que se evapore si la temperatura máxima del suelo es de 88°F.

Lo primero que hay que decir es que la recta predice solo dentro del rango de mínimo y máximo de la variable regresora.

# Modelo 1: Aporacion de Agua vs. Temperatura de aire
c(min(evaporacion$tpsuelma), max(evaporacion$tpsuelma))
## [1] 73 93

El valor x = 88 está dentro del rango, por lo tanto la predicción es válida. El valor predicho para x = 88 es:

x <- 88
prediccion_88 <- -160.413 + 2.284*x
prediccion_88
## [1] 40.579

El valor predicho para x=88 es 40.579 de agua evaporada.


Ejercicio 3

En 1986, debido a un incendio en una de las piezas de sus propulsores el transbordador espacial Challenger tuvo un accidente catastrófico. En lanzamientos de transbordadores espaciales lanzados con anterioridad, se habían inspeccionado los propulsores de las naves, y en algunos de ellos se habían encontrando defectos. El archivo orings.txt contiene observaciones de las siguientes variables: mission, es un indicador de número de misión del transbordador espacial lanzado; defecto, que toma valores 1 ó 0, en función de si se encontraron defectos o no en los propulsores; y temp, la temperatura (en grados Celcius) al momento del lanzamiento.

1.3.1. Actividades

  1. Ajuste un modelo con el cual sea posible predecir la probabilidad de encontrar defectos en función de la temperatura.

  2. ¿Se puede afirmar con un nivel de significación α = 0,05 que la temperatura influye en la probabilidad de que los propulsores tengan defectos?

  3. Interprete el valor del coeficientes estimados para la variable temperatura. A medida que aumenta la temperatura, ¿Cómo se comporta la probabilidad de encontrar defectos?

  4. ¿Cuál es la probabilidad estimada de falla si la temperatura en el momento del lanzamiento es de 24°?




Respuestas:

1. Ajuste un modelo con el cual sea posible predecir la probabilidad de encontrar defectos en función de la temperatura.

El modelo ha utilizar es un modelo de regresión logística.

fit_transbordador <- glm(defecto ~ temp, data = transbordador, family = binomial(link = "logit") )
summary(fit_transbordador)
## 
## Call:
## glm(formula = defecto ~ temp, family = binomial(link = "logit"), 
##     data = transbordador)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.0611  -0.7613  -0.3783   0.4524   2.2175  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)   7.6137     3.9334   1.936   0.0529 .
## temp         -0.4179     0.1948  -2.145   0.0320 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 28.267  on 22  degrees of freedom
## Residual deviance: 20.315  on 21  degrees of freedom
## AIC: 24.315
## 
## Number of Fisher Scoring iterations: 5

Validación del Modelo

Se realiza un test de Hosmer-Lemeshaw para validar el modelo.
Ho) Hay bondad y ajuste.
Hi) No hay bondad y ajuste.


hoslem.test(transbordador$defecto,fitted.values(fit_transbordador), g = 8)
## 
##  Hosmer and Lemeshow goodness of fit (GOF) test
## 
## data:  transbordador$defecto, fitted.values(fit_transbordador)
## X-squared = 6.8321, df = 6, p-value = 0.3367
#Nota: Se modificó la cantidad de grupos para obtener resultados satisfactorios.

El valor p es 0.2723. El valor alfa es 0.05.
El valor p es mayor que el valor de alfa. No se rechaza Ho. Es modelo presenta un buen ajuste en términos estadíticos.
Es válido utilizar el modelo de regresión lineal simple sobre estos datos.


Gráfico del Modelo de Regresión Logística

x <- data.frame(temp = seq(min(transbordador$temp),max(transbordador$temp),by = 0.005))
plot(transbordador$temp, transbordador$defecto, main = "Gráfico de Regresión Logística Simple", xlab = "Temperatura", ylab = "Prob. de defecto")
prob <- predict(fit_transbordador, newdata = x, type = "response")
lines(seq(min(transbordador$temp),max(transbordador$temp),by = 0.005), prob, lwd=2)

2. ¿Se puede afirmar con un nivel de significación α = 0,05 que la temperatura influye en la probabilidad de que los propulsores tengan defectos?
Para la pendiente hay un juego de hipótesis:
Ho) La pendiente es igual a cero.
Hi) La pendiente es distinta de cero.

La estimación de la pendiente es -0.4179, con un valor p de 0.0320. El valor alfa es 0.05. Se puede afirmar que la temperatura tiene efectos sobre la probabilidad de defectos, y que una disminución de la temperatura se relaciona con un aumento de la probabilidad de defecto.

3. Interprete el valor del coeficientes estimados para la variable temperatura. A medida que aumenta la temperatura, ¿Cómo se comporta la probabilidad de encontrar defectos? El ser el modelo válido, la estimación de la pendiente de la regresión logística ser estadísticamente significativa y tener signo negativo y haber determinado como “0” al estado inicial (no defectuoso) y “1” al cambio de estado (defectuoso) esto implica que una disminución de la temperatura se relaciona con un aumento de la probabilidad de que los propulsores tengan fallas.

4. ¿Cuál es la probabilidad estimada de falla si la temperatura en el momento del lanzamiento es de 24°?

ordenada <- coef(fit_transbordador)[1]
pendiente <- coef(fit_transbordador)[2]

x_temperatura <- 24

estimado <- exp(ordenada + pendiente*x_temperatura)/(1+exp(ordenada+pendiente*x_temperatura))
estimado
## (Intercept) 
##  0.08198054

La probabilidad de tener una falla con 24 grados farenheit es de 0.08198054.
Lo redondeamos a 0.082.


Hubicamos la estimación en el gráfico.

x <- data.frame(temp = seq(min(transbordador$temp),max(transbordador$temp),by = 0.005))
plot(transbordador$temp, transbordador$defecto, main = "Gráfico de Regresión Logística Simple", xlab = "Temperatura", ylab = "Prob. de defecto")
prob <- predict(fit_transbordador, newdata = x, type = "response")
lines(seq(min(transbordador$temp),max(transbordador$temp),by = 0.005), prob, lwd=2)
points(x_temperatura, estimado, col = "red", lwd = 3)
matriz_valores <- matrix(NA, 3, 2)
matriz_valores[1,] <- c(x_temperatura, 0)
matriz_valores[2,] <- c(x_temperatura, estimado)
matriz_valores[3,] <- c(0, estimado)
lines(matriz_valores,lwd=2, col = "blue")
points(x_temperatura, estimado, col = "red", lwd = 3)