1 Introducción

Revisaremos la construcción y validación de los supuestos para modelos espaciales de temperaturas promedio en México, ver la Parte 1 de esta serie de ejemplos relacionados con Geoeestadística: https://rpubs.com/Humberto_Mtz_Bta/1157323. Utilizaremos la información de los Estados de Guanajuato (GTO) y Yucatán (YUC) quienes en el análisis distribucional con boxplots, exhiben similitud en variación y simetría, lo que favorece la aplicación correcta de técnicas estadísticas clásicas como ,por ejemplo, el análisis de varianza (ANOVA).

2 Instalación y activación de la paquetería R

list.of.packages <- c("tidyverse","ggplot2","writexl","here","glue","car","rstatix")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)


library(tidyverse)
library(ggplot2)
library(here)
library(glue)
library(writexl)
library(car)
library(rstatix)

3 Obtención y exploración de datos

rm(list = ls()) #Se limpia el ambiente de trabajo
datos <- read.csv('https://smn.conagua.gob.mx/tools/RESOURCES/com_archivo_datos_resumenes/202305010000TMed.csv')
datos %>%
  ggplot( aes(x=Edo, y=Tmed)) + 
  geom_boxplot() + 
  scale_x_discrete(name="Estado")

4 Delimitación de información

A continuación, mostramos la distribución de las temperaturas promedio en los Estados de Guanajuato y Yucatán, los cuales fueron seleccionados estratégicamente para esta segunda parte de análisis espacial. Note que en GTO se tiene posibles outliers, lo cuales nos darán oportunidad de estudiar este tipo de datos posteriormente.

datos1 <- datos %>% filter(Edo %in% c("GTO", "YUC") )
 str(datos1) 
#> 'data.frame':    73 obs. of  6 variables:
#>  $ Lon  : num  -102 -101 -101 -101 -101 ...
#>  $ Lat  : num  20.4 20 20.5 20.6 20.7 ...
#>  $ Clave: chr  "ABSGJ" "ACNGJ" "AEGGJ" "AMCGJ" ...
#>  $ Edo  : chr  "GTO" "GTO" "GTO" "GTO" ...
#>  $ Est  : chr  "Abasolo Gto." "Ac\xe1mbaro Gto." "Apaseo el Grande Gto." "Ameche Gto." ...
#>  $ Tmed : num  25.5 23.6 23 22 23.1 ...
#Boxplot
 datos1 %>%
  ggplot( aes(x=Edo, y=Tmed)) + 
  geom_boxplot() + 
  scale_x_discrete(name="Estado")

5 Estadísticos univariados para Tmed

######Resumen para GTO y YUC
#sort(datos$Tmed, decreasing = FALSE)
Tabla2 <- datos1 %>%
  group_by(Edo) %>% 
  summarize(
    Min = min(Tmed),
    Mediana = median(Tmed),
    Media = mean(Tmed),
    sd = sd(Tmed),
    Max = max(Tmed))
Tabla2
#> # A tibble: 2 × 6
#>   Edo     Min Mediana Media    sd   Max
#>   <chr> <dbl>   <dbl> <dbl> <dbl> <dbl>
#> 1 GTO    18.4    22.7  22.5  1.75  26.0
#> 2 YUC    27.4    29.7  29.6  1.17  32.3

6 Outliers

Mediante el uso de la función identify_outliers analizamos si existen valores atípicos en los datos.

datos1 %>%
  group_by(Edo) %>%
  identify_outliers(Tmed)
#> # A tibble: 1 × 8
#>   Edo     Lon   Lat Clave Est           Tmed is.outlier is.extreme
#>   <chr> <dbl> <dbl> <chr> <chr>        <dbl> <lgl>      <lgl>     
#> 1 GTO   -100.  20.2 CORGJ Coroneo Gto.  18.4 TRUE       FALSE

#Imputación por la mediana
datos1[which(datos1$Tmed<18.5),"Tmed"] <- median(datos1$Tmed)

7 Normalidad

7.1 Inspección visual/gráfica

h <- filter(datos1, Edo == "GTO")
m_h<-mean(h$Tmed)
std_h<-sqrt(var(h$Tmed))
hist(h$Tmed, xlab="Temperaturas promedio", ylab="Frecuencia", main="Distribución Temperaturas promedio para GTO ", prob=TRUE)
curve(dnorm(x, mean=m_h, sd=std_h), col="darkblue", add=TRUE)



c <- filter(datos1, Edo == "YUC")
m_c<-mean(c$Tmed)
std_c<-sqrt(var(c$Tmed))
hist(c$Tmed, xlab="Temperaturas promedio", ylab="Frecuencia", main="Distribución de Temperaturas promedio para YUC", prob=TRUE, ylim=c(0,.4))
curve(dnorm(x, mean=m_c, sd=std_c), col="darkblue", add=TRUE)

7.2 Inspección análitica



shapiro.test(h$Tmed)        #Prueba de normalidad para GTO
#> 
#>  Shapiro-Wilk normality test
#> 
#> data:  h$Tmed
#> W = 0.97606, p-value = 0.4135

shapiro.test(c$Tmed)        #Prueba de normalida para YUC
#> 
#>  Shapiro-Wilk normality test
#> 
#> data:  c$Tmed
#> W = 0.9868, p-value = 0.9824
  • NO SE RECHAZA \(H_0:\text{Hay normalidad}\), con la prueba de Shapiro-Wilk.

8 Homogeneidad de varianzas

datos1 %>%
  leveneTest(Tmed ~ Edo, data = .)      #Prueba de levene
#> Levene's Test for Homogeneity of Variance (center = median)
#>       Df F value Pr(>F)
#> group  1  1.8475 0.1784
#>       71

La prueba de levene NO RECHAZA la igualdad de varianzas, es decir tenemos homocedasticidad.

9 Comparación de medias usando ANOVA

La herramienta estadística de Análisis de la Varianza (ANOVA) fue desarrollada por Ronald Fisher a principios del siglo XX y su campo de aplicación inicial fue en la agricultura.

El ANOVA es una herramienta de análisis que particiona la variación total de los datos en un conjunto de fuentes de variación para determinar la igualdad (\(H_0\)) o diferencia (\(H_a\)) de medias de los grupos. Algunos casos de aplicación pueden ser:


Hipótesis estadísticas

Modelo simple en función del Estado : 

datos_red <- datos1[-73,]
Tmed_Modelo1 <- aov(formula = Tmed ~ Edo , data = datos_red)        #Prueba  clásica
summary(Tmed_Modelo1)
#>             Df Sum Sq Mean Sq F value Pr(>F)    
#> Edo          1  796.7   796.7   359.7 <2e-16 ***
#> Residuals   70  155.0     2.2                   
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Se aprecia que la temperatura promedio (Tmed) por Estado(Edo) es diferente para GTO y YUC, es decir se rechaza \(H_0\) con un nivel de significancia \(\alpha < 0.05\).

Modelo espacial en función del Estado y posición : 

Tmed_Modelo2 <- aov(formula = Tmed ~ Edo + Lon + Lat, data = datos_red)
summary(Tmed_Modelo2)
#>             Df Sum Sq Mean Sq F value   Pr(>F)    
#> Edo          1  796.7   796.7  492.27  < 2e-16 ***
#> Lon          1   20.8    20.8   12.82 0.000638 ***
#> Lat          1   24.2    24.2   14.97 0.000247 ***
#> Residuals   68  110.1     1.6                     
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Modelo espacial en función del Estado, posición y su interacción : 

Tmed_Modelo3 <- aov(formula = Tmed ~ Edo*Lon*Lat, data = datos_red)
summary(Tmed_Modelo3)
#>             Df Sum Sq Mean Sq F value   Pr(>F)    
#> Edo          1  796.7   796.7 575.623  < 2e-16 ***
#> Lon          1   20.8    20.8  14.995 0.000256 ***
#> Lat          1   24.2    24.2  17.509 8.89e-05 ***
#> Edo:Lon      1   16.9    16.9  12.208 0.000869 ***
#> Edo:Lat      1    3.3     3.3   2.372 0.128490    
#> Lon:Lat      1    0.3     0.3   0.182 0.670764    
#> Edo:Lon:Lat  1    1.0     1.0   0.753 0.388821    
#> Residuals   64   88.6     1.4                     
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

10 Selección del mejor modelo

library(AICcmodavg)
model.set <- list(Tmed_Modelo1, Tmed_Modelo2, Tmed_Modelo3)
model.names <- c("Por Estado", "Por estado y posición", "Con interacción edo y loc")
aictab(model.set, modnames = model.names)
#> 
#> Model selection based on AICc:
#> 
#>                           K   AICc Delta_AICc AICcWt Cum.Wt      LL
#> Con interacción edo y loc 9 240.15       0.00   0.94   0.94 -109.62
#> Por estado y posición     5 245.79       5.63   0.06   1.00 -117.44
#> Por Estado                3 265.91      25.75   0.00   1.00 -129.78

11 Verificación de supuestos

par(mfrow=c(2,2))
plot(Tmed_Modelo3)

par(mfrow=c(1,1))

12 Material de apoyo