Cargar las librerías que usaremos

library(tidyverse)
library(knitr)

Problema 7: Análisis de Sismos (Quakes)

Tenemos en R el data frame quakes que registra sismos en cierta área alrededor de FIJI a partir de 1964, que consideraremos como población.

1. Definir nuestra población

# Usamos el dataset 'quakes' de R y le añadimos un ID
poblacion_quakes <- quakes %>%
  mutate(id = row_number())

# Definir N (Tamaño de la población)
N_quakes <- nrow(poblacion_quakes)

2. Extraer una muestra aleatoria simple sin reemplazo de tamaño 150 para estimar la media de la magnitud de los sismos.

Primero, definimos nuestra población (N = r N_quakes) y fijamos una semilla (set.seed(123)) para que el muestreo sea reproducible. Luego, extraemos la muestra de n = 150.

# Fijar semilla para reproducibilidad
set.seed(123)
n_quakes <- 150

# Extraer la muestra aleatoria simple
muestra_quakes <- poblacion_quakes %>%
  sample_n(size = n_quakes, replace = FALSE)

# Estimar la media poblacional usando la media muestral
media_muestral_mag <- mean(muestra_quakes$mag)

# Comparar con la media poblacional real (parámetro)
media_poblacional_mag <- mean(poblacion_quakes$mag)

cat(sprintf("La media muestral (nuestro estimador) es: %f . \n",round(media_muestral_mag, 4) ))
## La media muestral (nuestro estimador) es: 4.580000 .
cat(sprintf("El parámetro real de la población es: %f . \n",round(media_poblacional_mag, 4) ))
## El parámetro real de la población es: 4.620400 .

3. Calcular el error estándar del promedio muestral

Usamos la fórmula del Error Estándar (SE) para la media, \(SE(\bar{x}) = \frac{s}{\sqrt{n}}\), donde \(s\) es la desviación estándar de la muestra.

# Calcular la desviación estándar de la muestra
s_mag <- sd(muestra_quakes$mag)

# Calcular el Error Estándar
se_media_mag <- s_mag / sqrt(n_quakes)

cat(sprintf("El error estándar del promedio muestral es: %f . \n",round(se_media_mag, 4) ))
## El error estándar del promedio muestral es: 0.032400 .

4. Calcular el intervalo de confianza para el parametro poblacional considerando una confianza del 90%

Usamos la fórmula del Intervalo de Confianza: \(\bar{x} \pm z_{\alpha/2} \cdot SE(\bar{x})\)

# Nivel de confianza y valor Z
confianza_mag <- 0.90
alfa_mag <- 1 - confianza_mag
z_valor_mag <- qnorm(1 - (alfa_mag / 2))

# Calcular el margen de error
margen_error_mag <- z_valor_mag * se_media_mag

# Calcular el Intervalo de Confianza
ic_lim_inf_mag <- media_muestral_mag - margen_error_mag
ic_lim_sup_mag <- media_muestral_mag + margen_error_mag

cat(sprintf("El intervalo de confianza del 90%% para la media de la magnitud es: [%.4f, %.4f].", ic_lim_inf_mag, ic_lim_sup_mag))
## El intervalo de confianza del 90% para la media de la magnitud es: [4.5268, 4.6332].

5. Definiendo a un sismo como grave si la magnitud es mayor o igual a 5, calcular las proporciones poblacional y muestral de la gravedad

Ahora estimaremos una proporción (p), no una media.

# 1. Calcular Proporción POBLACIONAL (Parámetro 'p')
p_poblacional_grave <- poblacion_quakes %>%
  summarise(
    total_graves = sum(mag >= 5),
    p_grave = total_graves / n()
  ) %>%
  pull(p_grave)

# 2. Calcular Proporción MUESTRAL (Estimador 'p_gorro')
p_muestral_grave <- muestra_quakes %>%
  summarise(
    total_graves = sum(mag >= 5),
    p_grave = total_graves / n()
  ) %>%
  pull(p_grave)

cat(sprintf("La proporción poblacional real (parámetro) es: %.4f\n", p_poblacional_grave))
## La proporción poblacional real (parámetro) es: 0.1980
cat(sprintf("La proporción muestral (estimador) es: %.4f\n", p_muestral_grave))
## La proporción muestral (estimador) es: 0.2000

6. Calcular el error estándar de este último estadístico

Usamos la fórmula del Error Estándar (SE) para una proporción, \(SE(\hat{p}) = \sqrt{\frac{\hat{p}(1-\hat{p})}{n}}\).

# Usamos el estimador p_muestral (p_gorro) para la fórmula
p_hat <- p_muestral_grave
q_hat <- 1 - p_hat

se_proporcion_grave <- sqrt((p_hat * q_hat) / n_quakes)

cat(sprintf("El error estándar del estadístico de proporción es: %.4f.\n", se_proporcion_grave))
## El error estándar del estadístico de proporción es: 0.0327.

Problema 9: Muestreo Estratificado (ChickWeight)

Suponga que tenemos la siguiente población con información del peso y dieta para pollos df <- data.frame(peso = ChickWeight\(weight, dieta = ChickWeight\)Diet).

1. Preparamos

# Definir nuestra población
poblacion_chick <- data.frame(peso = ChickWeight$weight, dieta = ChickWeight$Diet)

# Calcular parámetros poblacionales generales
N_chick <- nrow(poblacion_chick)
media_poblacional_peso <- mean(poblacion_chick$peso)

# Calcular parámetros de los ESTRATOS (Dieta)
# Necesitamos Nh (tamaño) y Sh (DE) de cada estrato para la asignación
poblacion_estratos_info <- poblacion_chick %>%
  group_by(dieta) %>%
  summarise(
    Nh = n(),      # N_h: Tamaño del estrato
    Sh = sd(peso)  # S_h: Desviación Estándar del estrato
  ) %>%
  mutate(
    Wh = Nh / N_chick # W_h: Peso del estrato (Nh / N)
  )

2. Extraer una muestra aleatoria simple sin reemplazo de tamaño 120 y estimar la media poblacional.

Primero, realizamos un Muestreo Aleatorio Simple (MAS) para tener una estimación base.

# Fijar semilla para reproducibilidad
set.seed(456)
n_chick_mas <- 120

# Extraer la muestra aleatoria simple
muestra_chick_mas <- poblacion_chick %>%
  sample_n(size = n_chick_mas, replace = FALSE)

# Estimar la media poblacional usando la media muestral
media_muestral_mas <- mean(muestra_chick_mas$peso)

cat(sprintf("La media poblacional real (Parámetro) es: %.4f\n", media_poblacional_peso))
## La media poblacional real (Parámetro) es: 121.8183
cat(sprintf("La media muestral (Estimador MAS) es: %.4f\n", media_muestral_mas))
## La media muestral (Estimador MAS) es: 120.9250

3. Si consideramos que los pollos que siguieron una determinada dieta como un estrato, construir muestras con asignación proporcional y con asignación optima.

Ahora aplicaremos Muestreo Estratificado, usando la dieta como estrato. Usaremos el mismo tamaño total de muestra n = 120 para poder comparar los métodos.

cat("La tabla de parámetros de nuestros estratos (calculada de la población) es:\n")
## La tabla de parámetros de nuestros estratos (calculada de la población) es:
kable(poblacion_estratos_info, 
      digits = 4, 
      col.names = c("Estrato (Dieta)", "N_h (Tamaño)", "S_h (DE)", "W_h (Peso)"))
Estrato (Dieta) N_h (Tamaño) S_h (DE) W_h (Peso)
1 220 56.6566 0.3806
2 120 71.6075 0.2076
3 120 86.5418 0.2076
4 118 68.8287 0.2042

Cálculo de Asignación (Proporcional y Óptima)

Asignación Proporcional: \(n_h = n \cdot W_h\)

Asignación Óptima (Neyman): \(n_h = n \cdot \frac{N_h S_h}{\sum N_k S_k}\)

n_total_estratificado <- 120

calculos_asignacion <- poblacion_estratos_info %>%
  mutate(
    # 1. Asignación Proporcional
    nh_prop = n_total_estratificado * Wh,
    
    # 2. Asignación Óptima (Neyman)
    Nh_x_Sh = Nh * Sh,
    sum_Nh_Sh = sum(Nh_x_Sh),
    nh_opt = n_total_estratificado * (Nh_x_Sh / sum_Nh_Sh)
  )

# El redondeo puede causar que la suma no sea exactamente n_total.
# Lo ajustamos redondeando y asegurando que la suma sea 120.
calculos_asignacion <- calculos_asignacion %>%
  mutate(
    nh_prop = round(nh_prop),
    nh_opt = round(nh_opt)
  )

# Ajuste por redondeo (si n_prop no suma 120, se ajusta)
if (sum(calculos_asignacion$nh_prop) != n_total_estratificado) {
  diff <- n_total_estratificado - sum(calculos_asignacion$nh_prop)
  calculos_asignacion$nh_prop[which.max(calculos_asignacion$nh_prop)] <- 
    calculos_asignacion$nh_prop[which.max(calculos_asignacion$nh_prop)] + diff
}
if (sum(calculos_asignacion$nh_opt) != n_total_estratificado) {
  diff <- n_total_estratificado - sum(calculos_asignacion$nh_opt)
  calculos_asignacion$nh_opt[which.max(calculos_asignacion$nh_opt)] <- 
    calculos_asignacion$nh_opt[which.max(calculos_asignacion$nh_opt)] + diff
}


cat("Tabla de Asignación de Muestras (n=120):\n")
## Tabla de Asignación de Muestras (n=120):
kable(calculos_asignacion %>% select(dieta, Nh, Sh, nh_prop, nh_opt), 
      digits = 2,
      col.names = c("Estrato", "N_h", "S_h", "n_h (Prop)", "n_h (Óptima)"))
Estrato N_h S_h n_h (Prop) n_h (Óptima)
1 220 56.66 46 38
2 120 71.61 25 26
3 120 86.54 25 31
4 118 68.83 24 25

Extracción de las Muestras Estratificadas

Ahora extraemos las muestras de la población poblacion_chick basándonos en los tamaños nh_prop y nh_opt que calculamos.

# Fijamos una nueva semilla para este muestreo
set.seed(789)

# 1. Muestra Proporcional
# Agrupamos por dieta y tomamos el n_h correspondiente a ese grupo
muestra_prop <- poblacion_chick %>%
  group_by(dieta) %>%
  sample_n(
    # Busca el 'nh_prop' para el grupo (dieta) actual
    size = calculos_asignacion$nh_prop[match(cur_group()$dieta, calculos_asignacion$dieta)],
    replace = FALSE
  ) %>%
  ungroup() 

# 2. Muestra Óptima
muestra_opt <- poblacion_chick %>%
  group_by(dieta) %>%
  sample_n(
    # Busca el 'nh_opt' para el grupo (dieta) actual
    size = calculos_asignacion$nh_opt[match(cur_group()$dieta, calculos_asignacion$dieta)],
    replace = FALSE
  ) %>%
  ungroup()

cat(sprintf("Tamaño de la muestra proporcional: %d\n", nrow(muestra_prop)))
## Tamaño de la muestra proporcional: 120
cat(sprintf("Tamaño de la muestra óptima: %d\n", nrow(muestra_opt)))
## Tamaño de la muestra óptima: 120

4. Estimar las medias poblacionales de los estratos

Ahora calculamos la media muestral para cada estrato (\(\bar{y}_h\)) de ambas muestras.

# Medias de estratos de la muestra PROPORCIONAL
medias_estratos_prop <- muestra_prop %>%
  group_by(dieta) %>%
  summarise(media_estrato_h = mean(peso))

cat("Medias de Estratos (Muestra Proporcional):\n")
## Medias de Estratos (Muestra Proporcional):
kable(medias_estratos_prop, 
      digits = 4,
      col.names = c("Estrato", "Media Muestral (y_bar_h)"))
Estrato Media Muestral (y_bar_h)
1 98.0652
2 116.7200
3 141.8000
4 135.6250 
# Medias de estratos de la muestra ÓPTIMA
medias_estratos_opt <- muestra_opt %>%
  group_by(dieta) %>%
  summarise(media_estrato_h = mean(peso))
  
cat("\nMedias de Estratos (Muestra Óptima):\n")
## 
## Medias de Estratos (Muestra Óptima):
kable(medias_estratos_opt, 
      digits = 4,
      col.names = c("Estrato", "Media Muestral (y_bar_h)"))
Estrato Media Muestral (y_bar_h)
1 107.5000
2 106.1923
3 146.6452
4 127.9200

5. Estimar media poblacional general con los dos tipos de muestra (proporcional y optimo).

Finalmente, estimamos la media poblacional general (\(\bar{y}_{st}\)) usando la fórmula: \(\bar{y}_{st} = \sum_{h=1}^{L} W_h \cdot \bar{y}_h\)

Usamos los pesos \(W_h\) de la población (que calculamos en poblacion_estratos_info) y las medias \(\bar{y}_h\) de cada muestra (calculadas en el paso 4).

# Unir los pesos (Wh) con las medias de estratos (y_bar_h)

# --- Estimación Proporcional ---
est_final_prop <- medias_estratos_prop %>%
  # Unimos con la tabla de info de estratos para obtener Wh
  left_join(poblacion_estratos_info %>% select(dieta, Wh), by = "dieta") %>%
  # Calculamos el producto ponderado
  mutate(Wh_x_media_h = Wh * media_estrato_h)

# La estimación general es la suma de los productos ponderados
media_general_prop <- sum(est_final_prop$Wh_x_media_h)

# --- Estimación Óptima ---
est_final_opt <- medias_estratos_opt %>%
  left_join(poblacion_estratos_info %>% select(dieta, Wh), by = "dieta") %>%
  mutate(Wh_x_media_h = Wh * media_estrato_h)
  
media_general_opt <- sum(est_final_opt$Wh_x_media_h)

Resultados Finales de Estimación

cat(sprintf("La media poblacional real (Parámetro) es: %.4f\n", media_poblacional_peso))
## La media poblacional real (Parámetro) es: 121.8183
cat(sprintf("La estimación general con Muestreo Proporcional es: %.4f\n", media_general_prop))
## La estimación general con Muestreo Proporcional es: 118.6860
cat(sprintf("La estimación general con Muestreo Óptimo es: %.4f\n", media_general_opt))
## La estimación general con Muestreo Óptimo es: 119.5243
cat(sprintf("(Nota: En este caso, la Muestra Óptima dio una estimación ligeramente más cercana al parámetro real que la Proporcional. El Muestreo Aleatorio Simple (MAS) del punto 2 dio %.4f).\n", media_muestral_mas))
## (Nota: En este caso, la Muestra Óptima dio una estimación ligeramente más cercana al parámetro real que la Proporcional. El Muestreo Aleatorio Simple (MAS) del punto 2 dio 120.9250).