A continuación se muestra el caso y los datos con los que se trabajará el análisis estadístico para el examen.

Se realizó una investigación para conocer el Índice de Masa Corporal (IMC) de cuatro poblaciones distintas ubicadas en el Sur de Jalisco, una vez creado el estudio y el diseño, el tamaño de muestra arrojo la cantidad de 30 personas. A continuación se describen las variables:

Población=(c("Sayula","Sayula","Sayula","Sayula","Sayula","Sayula","Sayula","Sayula","Sayula","Sayula","Sayula","Sayula","Sayula","Sayula","Sayula","Sayula","Sayula","Sayula","Sayula","Sayula","Sayula","Sayula","Sayula","Sayula","Sayula","Sayula","Sayula","Sayula","Sayula","Sayula","GomezFarias","GomezFarias","GomezFarias","GomezFarias","GomezFarias","GomezFarias","GomezFarias","GomezFarias","GomezFarias","GomezFarias","GomezFarias","GomezFarias","GomezFarias","GomezFarias","GomezFarias","GomezFarias","GomezFarias","GomezFarias","GomezFarias","GomezFarias","GomezFarias","GomezFarias","GomezFarias","GomezFarias","GomezFarias","GomezFarias","GomezFarias","GomezFarias","GomezFarias","GomezFarias","Zacoalco","Zacoalco","Zacoalco","Zacoalco","Zacoalco","Zacoalco","Zacoalco","Zacoalco","Zacoalco","Zacoalco","Zacoalco","Zacoalco","Zacoalco","Zacoalco","Zacoalco","Zacoalco","Zacoalco","Zacoalco","Zacoalco","Zacoalco","Zacoalco","Zacoalco","Zacoalco","Zacoalco","Zacoalco","Zacoalco","Zacoalco","Zacoalco","Zacoalco","Zacoalco","Techaluta","Techaluta","Techaluta","Techaluta","Techaluta","Techaluta","Techaluta","Techaluta","Techaluta","Techaluta","Techaluta","Techaluta","Techaluta","Techaluta","Techaluta","Techaluta","Techaluta","Techaluta","Techaluta","Techaluta","Techaluta","Techaluta","Techaluta","Techaluta","Techaluta","Techaluta","Techaluta","Techaluta","Techaluta","Techaluta"))
IMC=(c(27,25,23,23,30,24,20,25,26,27,32,24,22,28,25,22,23,26,26,30,25,21,25,25,25,23,29,26,26,24,25,18,25,19,24,26,25,25,24,26,24,22,26,21,25,23,25,26,20,23,20,26,23,30,31,24,27,27,26,23,25,26,27,30,21,29,31,27,22,26,19,21,23,26,26,25,27,28,21,26,23,28,20,25,28,26,19,23,29,26,24,23,30,28,27,25,24,20,24,21,29,23,26,21,21,17,23,20,25,21,24,21,30,24,23,20,25,29,27,27))
datos <-data.frame (Población= Población, IMC= IMC)
head(datos)

Estadística descriptiva:

library(ggplot2)
library(gridExtra)

table(datos$Población)
## 
## GomezFarias      Sayula   Techaluta    Zacoalco 
##          30          30          30          30
aggregate(IMC ~ Población, data = datos, FUN = mean)
aggregate(IMC ~ Población, data = datos, FUN = sd)

Pruebas de normalidad

Gráficas

A continuación se muestran los gráficos QQ en los cuales visualemente observaremos la normalidad de los datos:

library(ggplot2)
par(mfrow = c(2,2))
qqnorm(datos[datos$Población == "Sayula","IMC"], main = "Sayula")
qqline(datos[datos$Población == "Sayula","IMC"], col="aquamarine4")
qqnorm(datos[datos$Población == "GomezFarias","IMC"], main = "GomezFarias")
qqline(datos[datos$Población == "GomezFarias","IMC"], col="aquamarine4")
qqnorm(datos[datos$Población == "Techaluta","IMC"], main = "Techaluta")
qqline(datos[datos$Población == "Techaluta","IMC"], col="aquamarine4")
qqnorm(datos[datos$Población == "Zacoalco","IMC"], main = "IMC")
qqline(datos[datos$Población == "Zacoalco","IMC"], col="aquamarine4")

par(mfrow = c(1,1))

Como se observa, pareciera que los datos siguen una distribución normal.

Formales

Para comprobar se realizarán las pruebas formales de normalidad con Shapiro Wilk test:

by(data = datos,INDICES = datos$Población,FUN = function(x){ shapiro.test(x$IMC)})
## datos$Población: GomezFarias
## 
##  Shapiro-Wilk normality test
## 
## data:  x$IMC
## W = 0.95488, p-value = 0.228
## 
## ------------------------------------------------------------ 
## datos$Población: Sayula
## 
##  Shapiro-Wilk normality test
## 
## data:  x$IMC
## W = 0.96395, p-value = 0.3892
## 
## ------------------------------------------------------------ 
## datos$Población: Techaluta
## 
##  Shapiro-Wilk normality test
## 
## data:  x$IMC
## W = 0.9644, p-value = 0.3992
## 
## ------------------------------------------------------------ 
## datos$Población: Zacoalco
## 
##  Shapiro-Wilk normality test
## 
## data:  x$IMC
## W = 0.95285, p-value = 0.2014

Como los resultados muestran un valor de p mayor a 0.05, nos indica que los datos son normales.

Pruebas de Homocedasticidad

Se realizan las pruevas de homocedasticidad para determinar si es un modelo equilibrado y como se comportan los datos, homogeneidad de varianzas.

Gráficas

A manera de graficas se pueden observar la homocedasticidad:

library(ggplot2)
ggplot(data = datos, aes(x = Población, y = IMC, color = Población)) +    geom_boxplot() +    theme_bw()

Observamos que las varianzas son similares en los grupos, siguiendo una distribución simétrica.

Formales

Para la prueba formal se realizará test de Levene o test de Fligner-Killeen.y , para ver la homocedasticidad.

fligner.test(IMC ~ Población,datos)
## 
##  Fligner-Killeen test of homogeneity of variances
## 
## data:  IMC by Población
## Fligner-Killeen:med chi-squared = 1.8859, df = 3, p-value = 0.5964
library(car)
leveneTest(IMC ~ Población,datos,center = "median")

El resultado es p>0.05 por lo que no hay indicios de que no haya homocedasticidad.

Prueba ANOVA o Kruskal Wallis

Se procede a comparar las diferencias entre medias con la prueba ANOVA o estudiar los posibles efectos de los factores sobre la varianza de una variable.

anova <- aov(datos$IMC ~ datos$Población)
summary(anova)
##                  Df Sum Sq Mean Sq F value Pr(>F)
## datos$Población   3   30.1  10.031   1.081   0.36
## Residuals       116 1076.2   9.278
plot(anova)

El resultado es p>0.05 significa que no hay evidencias suficientes para considerar que al menos dos medias son distintas. Además, la representación gráfica de los residuos no muestra falta de homocedasticidad.

Tamaño del efecto del ANOVA con ETA Cuadrado

Describe la proporción de variabilidad total atribuible a un factor:

eta_cuadrado <- 30.1/(30.1 + 1076.2)
eta_cuadrado
## [1] 0.02720781

Lo que me indica que es significativo pero debil, por lo que se comparan los diferentes grupos:

pairwise.t.test(x = datos$IMC, g = datos$Población, p.adjust.method = "holm",
                pool.sd = TRUE, paired = FALSE, alternative = "two.sided")
## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  datos$IMC and datos$Población 
## 
##           GomezFarias Sayula Techaluta
## Sayula    0.96        -      -        
## Techaluta 1.00        0.84   -        
## Zacoalco  0.96        1.00   0.96     
## 
## P value adjustment method: holm
TukeyHSD(anova)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = datos$IMC ~ datos$Población)
## 
## $`datos$Población`
##                             diff       lwr       upr     p adj
## Sayula-GomezFarias     0.9333333 -1.116714 2.9833803 0.6364219
## Techaluta-GomezFarias -0.2333333 -2.283380 1.8167136 0.9908727
## Zacoalco-GomezFarias   0.8000000 -1.250047 2.8500470 0.7397154
## Techaluta-Sayula      -1.1666667 -3.216714 0.8833803 0.4507644
## Zacoalco-Sayula       -0.1333333 -2.183380 1.9167136 0.9982579
## Zacoalco-Techaluta     1.0333333 -1.016714 3.0833803 0.5560019
plot(TukeyHSD(anova))

No se encuentrán diferencias significativas entre las medias de los grupos, por lo tanto no se hacen pruebas post-hoc.

Pruebas de correlación

Para observar la relación de las variables, en este caso si se correlacionan los valores de IMC entre cada uno de los grupos.

A continuación, opté por volver a subir la matríz de datos para separar los datos de las poblaciones como variables:

Sayula=c(27,25,23,23,30,24,20,25,26,27,32,24,22,28,25,22,23,26,26,30,25,21,25,25,25,23,29,26,26,24)
GomezFarias=c(25,18,25,19,24,26,25,25,24,26,24,22,26,21,25,23,25,26,20,23,20,26,23,30,31,24,27,27,26,23)
Zacoalco=c(25,26,27,30,21,29,31,27,22,26,19,21,23,26,26,25,27,28,21,26,23,28,20,25,28,26,19,23,29,26)
Techaluta=c(24,23,30,28,27,25,24,20,24,21,29,23,26,21,21,17,23,20,25,21,24,21,30,24,23,20,25,29,27,27)
Lugar <-data.frame (Sayula= Sayula, GomezFarias=GomezFarias, Zacoalco=Zacoalco, Techaluta=Techaluta)
head(Lugar)

Dado que son múltiples variables, se realizará una patriz de correlación.

Se generó gráficos de dispersión dos a dos:

pairs(x = Lugar, lower.panel = NULL)

A simple vista no se observan correlaciones entre variables.

Posteriormente se realiza la pruba de correlación de Pearson, dado que los datos son paramétricos.

cor(x = Lugar, method = "pearson")
##                  Sayula GomezFarias   Zacoalco  Techaluta
## Sayula       1.00000000 -0.01810688 -0.5322067  0.1633216
## GomezFarias -0.01810688  1.00000000  0.1334118 -0.0202623
## Zacoalco    -0.53220670  0.13341175  1.0000000 -0.2962594
## Techaluta    0.16332161 -0.02026230 -0.2962594  1.0000000

Como se puede observar hay correlaciones positivas y negativas, además que en algunas variables no existe correlación:

Observémoslo de manera gráfica:

library(psych)
## 
## Attaching package: 'psych'
## The following object is masked from 'package:car':
## 
##     logit
## The following objects are masked from 'package:ggplot2':
## 
##     %+%, alpha
multi.hist(x = Lugar, dcol = c("aquamarine4", "darkslategray"), dlty = c("dotted", "solid"), 
           main = "")

library(corrplot)
## corrplot 0.92 loaded
corrplot(corr = cor(x = Lugar, method = "pearson"), method = "number",
         tl.cex = 0.7,number.cex = 0.8, cl.pos = "n")

En la última gráfica se puede observar en el Heat map que solo hay una relación inversa moderada entre Zacoalco y Sayula (-0.5) y relación directa debil entre Techaluta y Zacoalco (0.3), muy debil entre Goméz Farías y Techaluta (0.13) además de Goméz y Zacoalco (0.13).

library(psych)
pairs.panels(x = Lugar, ellipses = FALSE, lm = TRUE, method = "pearson")

library(GGally)
## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2
library(psych)
ggpairs(Lugar, lower = list(continuous = "smooth"), 
        diag = list(continuous = "bar"), axisLabels = "none")
## Warning in check_and_set_ggpairs_defaults("diag", diag, continuous =
## "densityDiag", : Changing diag$continuous from 'bar' to 'barDiag'
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

En los gráficos anteriones podemos observar diagramas de dispersión, en los que se ve que la relación no es clara entre variables.

Ecuación de regresión, con análisis de residuos.

La regresión lineal simple consiste en generar un modelo de regresión (ecuación de una recta) que permita explicar la relación lineal que existe entre dos variables, en este caso la población y el IMC. Se usarán solo las variables de Sayula y Zacoalco para el modelo, dado que fueron las que tuvieron mayor correlación.

Se inicia con la representación gráfica de intuir si existe una relación y cuantificar dicha relación mediante un coeficiente de correlación.

library(ggplot2)
ggplot(data = Lugar, mapping = aes(x = Sayula, y = Zacoalco)) +
  geom_point(color = "darkslategray", size = 2) +
  labs(title  =  'Diagrama de dispersión', x  =  'Sayula') +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5))

A simple vista no se ve que exista una relación, pero por motivos del examen se seguirá el proceso. Además, como se observó en el apartado anterior la correlación entre las variables no es fuerte.

modelo_lineal <- lm(Sayula ~ Zacoalco, Lugar)

summary(modelo_lineal)
## 
## Call:
## lm(formula = Sayula ~ Zacoalco, data = Lugar)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.1686 -1.3872 -0.1098  1.5561  5.1675 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  36.4118     3.3875  10.749 1.91e-11 ***
## Zacoalco     -0.4454     0.1339  -3.326  0.00247 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.337 on 28 degrees of freedom
## Multiple R-squared:  0.2832, Adjusted R-squared:  0.2576 
## F-statistic: 11.06 on 1 and 28 DF,  p-value: 0.002468

El valor de R2 indica que el modelo calculado explica el 28.32% de la variabilidad presente en la variable Sayula mediante la variable Zacoalco.

El p-value obtenido en el test F (0.002468) determina que sí es significativamente superior la varianza explicada por el modelo en comparación a la varianza total.

Formula de regresión lineal: Sayula=36.4118-0.4454 Zacoalco

confint(modelo_lineal)
##                  2.5 %    97.5 %
## (Intercept) 29.4727768 43.350765
## Zacoalco    -0.7196082 -0.171104
ggplot(data = Lugar, mapping = aes(x = Sayula, y = Zacoalco)) +
  geom_point(color = "aquamarine4", size = 2) +
  labs(title  =  'Sayula ~ Zacoalco', x  =  'Zacoalco') +
  geom_smooth(method = "lm", se = FALSE, color = "firebrick") +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5))
## `geom_smooth()` using formula = 'y ~ x'

puntos <- seq(from = min(Lugar$Sayula),
              to = max(Lugar$Sayula),
              length.out = 100)
limites_intervalo <- predict(object = modelo_lineal,
                             newdata = data.frame(Sayula = puntos),
                             interval = "confidence", level = 0.95)
## Warning: 'newdata' had 100 rows but variables found have 30 rows
head(limites_intervalo, 3) 
##        fit      lwr      upr
## 1 25.27787 24.40341 26.15233
## 2 24.83251 23.92430 25.74073
## 3 24.38716 23.36959 25.40473
library(ggplot2)
plot(Lugar$Sayula, Lugar$Zacoalco, col = "firebrick4", pch = 19, ylab = "Zacoalco",
     xlab = "Sayula", main = 'Sayula ~ Zacoalco')

library(ggplot2)
ggplot(data = Lugar, mapping = aes(x = Sayula, y = Zacoalco)) +
  geom_point(color = "aquamarine4", size = 2) +
  labs(title  =  'Diagrama de dispersión', x  =  'Sayula') +
  geom_smooth(method = "lm", se = TRUE, color = "black") +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5))
## `geom_smooth()` using formula = 'y ~ x'

Se predice el valor de la variable y junto con su intervalo de confianza, genera la regresión y su intervalo de forma directa.

Verificar condiciones para poder aceptar un modelo lineal:

Lugar$prediccion <- modelo_lineal$fitted.values
Lugar$residuos   <- modelo_lineal$residuals
head(Lugar)

A continuación se hace el análisis de residuos:

ggplot(data = Lugar, aes(x = prediccion, y = residuos)) +
  geom_point(aes(color = residuos)) +
  scale_color_gradient2(low = "aquamarine4", mid = "darkslategray", high = "firebrick4") +
  geom_hline(yintercept = 0) +
  geom_segment(aes(xend = prediccion, yend = 0), alpha = 0.2) +
  labs(title = "Distribución de los residuos", x = "predicción modelo",
       y = "residuo") +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5), legend.position = "none")

Los residuos se distribuyen de forma aleatoria entorno al 0 por lo que se acepta la linealidad.

ggplot(data = Lugar, aes(x = residuos)) +
  geom_histogram(aes(y = ..density..)) +
  labs(title = "Histograma de los residuos") +
  theme_light()
## Warning: The dot-dot notation (`..density..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(density)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Los residuos siguen una distribución normal:

qqnorm(modelo_lineal$residuals)
qqline(modelo_lineal$residuals)

shapiro.test(modelo_lineal$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  modelo_lineal$residuals
## W = 0.98506, p-value = 0.9381

En efecto los residuos son normales.

ggplot(data = Lugar, aes(x = prediccion, y = residuos)) +
  geom_point(aes(color = residuos)) +
  scale_color_gradient2(low = "aquamarine4", mid = "darkslategray", high = "firebrick4") +
  geom_segment(aes(xend = prediccion, yend = 0), alpha = 0.2) +
  geom_smooth(se = FALSE, color = "firebrick") +
  labs(title = "Distribución de los residuos", x = "predicción modelo",
       y = "residuo") +
  geom_hline(yintercept = 0) +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5), legend.position = "none")
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

El contraste de hipótesis muestra que hay homocedasticidad.

library(ggplot2)
ggplot(data = Lugar, aes(x = seq_along(residuos), y = residuos)) +
  geom_point(aes(color = residuos)) +
  scale_color_gradient2(low = "aquamarine4", mid = "darkslategray", high = "firebrick4") +
  geom_line(size = 0.3) +
  labs(title = "Distribución de los residuos", x = "index", y = "residuo") +
  geom_hline(yintercept = 0) +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5), legend.position = "none")
## 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.

La representación de los residuos no muestra ninguna tendencia.

library(car)
summary(influence.measures(model = modelo_lineal))
## Potentially influential observations of
##   lm(formula = Sayula ~ Zacoalco, data = Lugar) :
## 
##    dfb.1_ dfb.Zclc dffit   cov.r   cook.d hat  
## 11  0.80  -0.75     0.85_*  0.97    0.33   0.16
## 20 -0.07   0.13     0.47    0.75_*  0.09   0.04
## 27  0.19  -0.18     0.21    1.25_*  0.02   0.16
influencePlot(model = modelo_lineal)

La observación 11 como atípica.

library(ggplot2)
ggplot(data = Lugar, mapping = aes(x = Sayula, y = Zacoalco)) +
  geom_point(color = "aquamarine4", size = 2) +
  geom_smooth(method = "lm", se = FALSE, color = "darkslategray") +
  #se resalta el valor excluido
  geom_point(data = Lugar[2,], color = "firebrick4", size = 2) +
  #se añade la nueva recta de mínimos cuadrados
  geom_smooth(data = Lugar[-2,], method = "lm", se = FALSE, color = "skyblue4") +
  labs(title  =  'Diagrama de dispersión', x  =  'Sayula') +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5), legend.position = "none")
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'

La exclusión del valor atípoco no cambió la distribución.

lm(formula = Sayula ~ Zacoalco, data = Lugar)$coefficients
## (Intercept)    Zacoalco 
##  36.4117711  -0.4453561
lm(formula = Sayula ~ Zacoalco, data = Lugar[-2,])$coefficients
## (Intercept)    Zacoalco 
##  36.4188606  -0.4458693

4 gráficos de evaluación de recursos:

par(mfrow = c(1,2))
plot(modelo_lineal)

Preguntas

  1. ¿Existe diferencia significativa entre el IMC de las muestras? No.

  2. ¿Cuáles son las localidades que mejor se correlacionan? Zacoalco y Sayula se correlacionan negativa y moderadamente.

  3. ¿Cuál es la ecuación que modela el comportamiento del IMC de esas dos localidades? Sayula=36.4188606 -0.4458693 Zacoalco

  4. A qué conclusión llega.

    -No hay diferencias significativas entre medias (F=1.081).

    -No hay correlación fuerte entre la variable de población e IMC. Solo el IMC de Zacoalco y Sayula se correlacionan negativa y moderadamente(-0.5).

    -La variable de Población no es predictiva del IMC.

    -El modelo lineal entre las dos variables con mayor correlación fue Sayula=36.4188606 -0.4458693 Zacoalco.