Problema 1

n <- 10000 # número de coordenadas a generar

# Se genera n coordenadas de distribución uniforme entre 0 y 1
x <- runif(n, min = 0,max = 1)
y <- runif(n, min = 0,max = 1)

# Se genera un dataFrame con las coordenadas
coordenadas <- data.frame(x=x,y=y)
head(coordenadas)
##           x          y
## 1 0.1329068 0.01774489
## 2 0.3861688 0.04087019
## 3 0.9621572 0.54587758
## 4 0.1039092 0.72781064
## 5 0.4416577 0.50726031
## 6 0.7675104 0.97995288
# Se grafican las coordenadas 
plot(coordenadas$x, coordenadas$y, 
     main = "Gráfico de dispersión de coordenadas", 
     xlab = "X", 
     ylab = "Y", 
     pch = 19, 
     col = "blue",
     asp=1,
     xlim = c(0, 1), 
     ylim = c(0, 1))
# Se dibuja un circulo de radio 0.5 centrado en (0.5,0.5)
symbols(0.5, 0.5, circles = 0.5, inches = FALSE, add = TRUE)

#Se añade cuadricula para facilidad de lectura
grid()

# Se calculan las distancias de cada punto al centro del circulo
distancias <- sqrt((coordenadas$x - 0.5)^2 + (coordenadas$y - 0.5)^2)

# Se cuentan los puntos que estan dentro del círculo
puntos_dentro <- sum(distancias <= 0.5)
puntos_dentro
## [1] 7897
# Se calcula la probabilidad de que cada punto esté dentro del círculo
probabilidad <- puntos_dentro / n
probabilidad
## [1] 0.7897
# Se estima el valor de pi sabiendo que el área del circulo es de pi/4
estimacion_pi <- probabilidad*4
estimacion_pi
## [1] 3.1588

Problema 2

# Cargar las librerías necesarias
library(ggplot2)
library(reshape2)

# Función para simular y evaluar los estimadores
evaluar_estimadores <- function(n_muestras, n, theta) {
  # Crear vectores para almacenar los valores de cada estimador en las simulaciones
  estimador1 <- numeric(n_muestras)
  estimador2 <- numeric(n_muestras)
  estimador3 <- numeric(n_muestras)
  estimador4 <- numeric(n_muestras)
  
  # Simular las muestras y calcular los estimadores
  for (i in 1:n_muestras) {
    muestra <- rexp(n, rate = 1/theta)  # Muestra de una distribución exponencial con parámetro theta
    
    # Calcular los estimadores
    estimador1[i] <- (muestra[1] + muestra[2])/6 + (muestra[3] + muestra[4])/3
    estimador2[i] <- (muestra[1] + 2*muestra[2] + 3*muestra[3] + 4*muestra[4])/5
    estimador3[i] <- mean(muestra)  # La media muestral
    estimador4[i] <- (min(muestra) + max(muestra)) / 2
  }
  
  # Calcular la esperanza (media) y la varianza de cada estimador
  mean_estimador1 <- mean(estimador1)
  mean_estimador2 <- mean(estimador2)
  mean_estimador3 <- mean(estimador3)
  mean_estimador4 <- mean(estimador4)
  
  var_estimador1 <- var(estimador1)
  var_estimador2 <- var(estimador2)
  var_estimador3 <- var(estimador3)
  var_estimador4 <- var(estimador4)
  
  # Guardar los resultados en una lista
  resultados <- list(
    medias = c(mean_estimador1, mean_estimador2, mean_estimador3, mean_estimador4),
    varianzas = c(var_estimador1, var_estimador2, var_estimador3, var_estimador4),
    estimadores = data.frame(
      Estimador1 = estimador1,
      Estimador2 = estimador2,
      Estimador3 = estimador3,
      Estimador4 = estimador4
    )
  )
  
  return(resultados)
}

# Parámetros
n_muestras <- 10000  # Número de muestras simuladas
theta <- 2  # Parámetro verdadero de la distribución exponencial
tamano_muestras <- c(20, 50, 100, 1000)  # Tamaños de muestras a evaluar

# Evaluar para diferentes tamaños de muestra
resultados_lista <- list()

for (n in tamano_muestras) {
  resultados_lista[[as.character(n)]] <- evaluar_estimadores(n_muestras, n, theta)
}

# Visualización de los resultados: Comparación de la media (insesgadez)
for (n in tamano_muestras) {
  cat("\nResultados para n =", n, ":\n")
  cat("Media de los estimadores:\n")
  print(resultados_lista[[as.character(n)]]$medias)
  cat("Varianza de los estimadores:\n")
  print(resultados_lista[[as.character(n)]]$varianzas)
}
## 
## Resultados para n = 20 :
## Media de los estimadores:
## [1] 2.007335 4.016012 2.012611 3.672154
## Varianza de los estimadores:
## [1] 1.1477793 4.9595578 0.2028293 1.6317815
## 
## Resultados para n = 50 :
## Media de los estimadores:
## [1] 1.993721 3.988572 2.001853 4.509132
## Varianza de los estimadores:
## [1] 1.07107205 4.61595961 0.07856664 1.59799883
## 
## Resultados para n = 100 :
## Media de los estimadores:
## [1] 1.980706 3.961987 1.999606 5.198852
## Varianza de los estimadores:
## [1] 1.09365917 4.74059289 0.03931599 1.58820434
## 
## Resultados para n = 1000 :
## Media de los estimadores:
## [1] 2.009266 4.015224 2.000539 7.490045
## Varianza de los estimadores:
## [1] 1.106545971 4.768575169 0.004020642 1.620953122
# Crear boxplots para cada tamaño de muestra

for (n in tamano_muestras) {
  # Añadir una columna ficticia "ID" para melt()
  resultados <- resultados_lista[[as.character(n)]]$estimadores
  resultados$ID <- 1:n_muestras  # Crear columna ID
  
  # Transformar el data frame para el gráfico
  estimadores_long <- melt(resultados, id.vars = "ID")  # Usar "ID" como identificador
  
  # Crear boxplot comparativo
  p <- ggplot(estimadores_long, aes(x = variable, y = value)) +
    geom_boxplot(fill = "lightblue", color = "black") +
    geom_hline(yintercept = theta, linetype = "dashed", color = "red") +  # Línea del valor verdadero de theta
    labs(title = paste("Comparación de los Estimadores para n =", n),
         x = "Estimadores",
         y = "Valor") +
    theme_minimal()
  
  print(p)
}

# Los estimadores 1 y 3 son insesgados ya que los valores que toman a traves de las simulaciones se acercan al valor real del parámetro
# El tercer estimador es el mas eficiente ya que es el que tiene menor varianza en su resultado
# Los estimadores 1 y 3 son los mas consistentes ya que, a medida que el tamaño de la muestra aumenta, estos estimadores convergen al dato real del parámetro (2)

Problema 3

# Función para generar una muestra y calcular el estimador de la proporción muestral
calcular_proporcion_muestral <- function(poblacion, n_muestra) {
  # Obtener una muestra aleatoria de tamaño n_muestra
  muestra <- sample(poblacion, size = n_muestra, replace = FALSE)
  
  # Calcular el estimador de la proporción muestral (pˆ)
  p_estimado <- mean(muestra)
  
  # Retornar el resultado
  return(p_estimado)
}

# Simulación de la población
set.seed(123)  # Fijar semilla para reproducibilidad
n_poblacion <- 1000  # Tamaño de la población
p_enfermas <- 0.50  # Proporción de plantas enfermas

# Generar la población donde 1 = enferma, 0 = sana
poblacion <- rbinom(n_poblacion, size = 1, prob = p_enfermas)

# Ejemplo: Calcular el estimador de la proporción muestral para diferentes tamaños de muestra
n_muestra <- 100  # Tamaño de la muestra
p_estimado <- calcular_proporcion_muestral(poblacion, n_muestra)

# Mostrar el resultado
cat("El estimador de la proporción muestral pˆ para una muestra de tamaño", n_muestra, "es:", p_estimado, "\n")
## El estimador de la proporción muestral pˆ para una muestra de tamaño 100 es: 0.53
# Cargar las librerías necesarias
library(ggplot2)

# Función para generar una muestra y calcular el estimador de la proporción muestral
calcular_proporcion_muestral <- function(poblacion, n_muestra) {
  # Obtener una muestra aleatoria de tamaño n_muestra
  muestra <- sample(poblacion, size = n_muestra, replace = FALSE)
  
  # Calcular el estimador de la proporción muestral (pˆ)
  p_estimado <- mean(muestra)
  
  # Retornar el resultado
  return(p_estimado)
}

# Simulación de la población
set.seed(123)  # Fijar semilla para reproducibilidad
n_poblacion <- 1000  # Tamaño de la población
p_enfermas <- 0.50  # Proporción de plantas enfermas

# Generar la población donde 1 = enferma, 0 = sana
poblacion <- rbinom(n_poblacion, size = 1, prob = p_enfermas)

# Parámetros para la simulación
n_muestra <- 100  # Tamaño de cada muestra
n_simulaciones <- 500  # Número de simulaciones

# Almacenar los resultados de los 500 estimadores pˆ
p_estimados <- numeric(n_simulaciones)

# Repetir el muestreo 500 veces
for (i in 1:n_simulaciones) {
  p_estimados[i] <- calcular_proporcion_muestral(poblacion, n_muestra)
}

# Calcular estadísticas descriptivas
mean_p <- mean(p_estimados)
sd_p <- sd(p_estimados)
skewness_p <- sum((p_estimados - mean_p)^3) / (n_simulaciones * sd_p^3)
kurtosis_p <- sum((p_estimados - mean_p)^4) / (n_simulaciones * sd_p^4) - 3  # Curtosis normalizada

# Mostrar los resultados descriptivos
cat("Media de los estimadores pˆ:", mean_p, "\n")
## Media de los estimadores pˆ: 0.4939
cat("Desviación estándar de los estimadores pˆ:", sd_p, "\n")
## Desviación estándar de los estimadores pˆ: 0.04647699
cat("Coeficiente de asimetría (skewness):", skewness_p, "\n")
## Coeficiente de asimetría (skewness): -0.1065394
cat("Curtosis (kurtosis):", kurtosis_p, "\n")
## Curtosis (kurtosis): 0.01873123
# Graficar el histograma de los resultados de pˆ
ggplot(data.frame(p_estimados), aes(x = p_estimados)) +
  geom_histogram(aes(y = ..density..), bins = 30, fill = "lightblue", color = "black") +
  geom_density(color = "red", size = 1) +
  geom_vline(aes(xintercept = mean_p), color = "blue", linetype = "dashed", size = 1) +
  labs(title = "Distribución de los 500 estimadores pˆ",
       x = "Estimador pˆ",
       y = "Densidad") +
  theme_minimal()
## 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.
## 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.

# Graficar el boxplot de los resultados
ggplot(data.frame(p_estimados), aes(y = p_estimados)) +
  geom_boxplot(fill = "lightblue", color = "black") +
  geom_hline(aes(yintercept = mean_p), color = "blue", linetype = "dashed", size = 1) +
  labs(title = "Boxplot de los 500 estimadores pˆ",
       y = "Estimador pˆ") +
  theme_minimal()

# Cargar las librerías necesarias
library(ggplot2)

# Función para generar una muestra y calcular el estimador de la proporción muestral
calcular_proporcion_muestral <- function(poblacion, n_muestra) {
  # Obtener una muestra aleatoria de tamaño n_muestra
  muestra <- sample(poblacion, size = n_muestra, replace = FALSE)
  
  # Calcular el estimador de la proporción muestral (pˆ)
  p_estimado <- mean(muestra)
  
  # Retornar el resultado
  return(p_estimado)
}

# Simulación de la población
set.seed(123)  # Fijar semilla para reproducibilidad
n_poblacion <- 1000  # Tamaño de la población
p_enfermas <- 0.50  # Proporción de plantas enfermas

# Generar la población donde 1 = enferma, 0 = sana
poblacion <- rbinom(n_poblacion, size = 1, prob = p_enfermas)

# Definir los tamaños de muestra
tamanos_muestra <- c(5, 10, 15, 20, 30, 50, 60, 100, 200, 500)
n_simulaciones <- 500  # Número de simulaciones

# Función para realizar las simulaciones y análisis de normalidad
analizar_normalidad <- function(poblacion, n_muestra, n_simulaciones) {
  # Almacenar los resultados de los estimadores pˆ
  p_estimados <- numeric(n_simulaciones)
  
  # Repetir el muestreo n_simulaciones veces
  for (i in 1:n_simulaciones) {
    p_estimados[i] <- calcular_proporcion_muestral(poblacion, n_muestra)
  }
  
  # Prueba de Shapiro-Wilk
  shapiro_result <- shapiro.test(p_estimados)
  
  # Graficar Q-Q plot para verificar la normalidad
  qqnorm(p_estimados, main = paste("Gráfico Q-Q para n =", n_muestra))
  qqline(p_estimados, col = "red")
  
  # Retornar los resultados de la prueba y los estimadores
  return(list(
    p_estimados = p_estimados,
    shapiro_p_value = shapiro_result$p.value
  ))
}

# Realizar las simulaciones para cada tamaño de muestra
resultados <- list()

for (n_muestra in tamanos_muestra) {
  cat("Análisis para tamaño de muestra n =", n_muestra, "\n")
  resultado <- analizar_normalidad(poblacion, n_muestra, n_simulaciones)
  
  # Mostrar el valor p de la prueba de Shapiro-Wilk
  cat("Valor p de la prueba de Shapiro-Wilk:", resultado$shapiro_p_value, "\n\n")
  
  # Guardar los resultados
  resultados[[as.character(n_muestra)]] <- resultado
}
## Análisis para tamaño de muestra n = 5

## Valor p de la prueba de Shapiro-Wilk: 3.333753e-15 
## 
## Análisis para tamaño de muestra n = 10

## Valor p de la prueba de Shapiro-Wilk: 7.877291e-10 
## 
## Análisis para tamaño de muestra n = 15

## Valor p de la prueba de Shapiro-Wilk: 4.472199e-08 
## 
## Análisis para tamaño de muestra n = 20

## Valor p de la prueba de Shapiro-Wilk: 7.172682e-07 
## 
## Análisis para tamaño de muestra n = 30

## Valor p de la prueba de Shapiro-Wilk: 2.618341e-05 
## 
## Análisis para tamaño de muestra n = 50

## Valor p de la prueba de Shapiro-Wilk: 0.003260273 
## 
## Análisis para tamaño de muestra n = 60

## Valor p de la prueba de Shapiro-Wilk: 0.005412488 
## 
## Análisis para tamaño de muestra n = 100

## Valor p de la prueba de Shapiro-Wilk: 0.04345206 
## 
## Análisis para tamaño de muestra n = 200

## Valor p de la prueba de Shapiro-Wilk: 0.08366892 
## 
## Análisis para tamaño de muestra n = 500

## Valor p de la prueba de Shapiro-Wilk: 0.1720724
# Cargar las librerías necesarias
library(ggplot2)

# Función para generar una muestra y calcular el estimador de la proporción muestral
calcular_proporcion_muestral <- function(poblacion, n_muestra) {
  # Obtener una muestra aleatoria de tamaño n_muestra
  muestra <- sample(poblacion, size = n_muestra, replace = FALSE)
  
  # Calcular el estimador de la proporción muestral (pˆ)
  p_estimado <- mean(muestra)
  
  # Retornar el resultado
  return(p_estimado)
}

# Definir los tamaños de muestra y número de simulaciones
tamanos_muestra <- c(5, 10, 15, 20, 30, 50, 60, 100, 200, 500)
n_simulaciones <- 500  # Número de simulaciones

# Función para realizar las simulaciones y análisis de normalidad
analizar_normalidad <- function(poblacion, n_muestra, n_simulaciones) {
  # Almacenar los resultados de los estimadores pˆ
  p_estimados <- numeric(n_simulaciones)
  
  # Repetir el muestreo n_simulaciones veces
  for (i in 1:n_simulaciones) {
    p_estimados[i] <- calcular_proporcion_muestral(poblacion, n_muestra)
  }
  
  # Prueba de Shapiro-Wilk
  shapiro_result <- shapiro.test(p_estimados)
  
  # Graficar Q-Q plot para verificar la normalidad
  qqnorm(p_estimados, main = paste("Gráfico Q-Q para n =", n_muestra))
  qqline(p_estimados, col = "red")
  
  # Retornar los resultados de la prueba y los estimadores
  return(list(
    p_estimados = p_estimados,
    shapiro_p_value = shapiro_result$p.value
  ))
}

# Función para ejecutar todo el análisis para un lote con una proporción específica
ejecutar_simulacion <- function(p_enfermas) {
  # Simulación de la población
  set.seed(123)  # Fijar semilla para reproducibilidad
  n_poblacion <- 1000  # Tamaño de la población
  
  # Generar la población donde 1 = enferma, 0 = sana
  poblacion <- rbinom(n_poblacion, size = 1, prob = p_enfermas)
  
  # Realizar las simulaciones para cada tamaño de muestra
  resultados <- list()
  
  for (n_muestra in tamanos_muestra) {
    cat("Análisis para tamaño de muestra n =", n_muestra, "\n")
    resultado <- analizar_normalidad(poblacion, n_muestra, n_simulaciones)
    
    # Mostrar el valor p de la prueba de Shapiro-Wilk
    cat("Valor p de la prueba de Shapiro-Wilk:", resultado$shapiro_p_value, "\n\n")
    
    # Guardar los resultados
    resultados[[as.character(n_muestra)]] <- resultado
  }
  
  return(resultados)
}

# Escenario 1: Lote con 10% de plantas enfermas
cat("=== Análisis para lote con 10% de plantas enfermas ===\n")
## === Análisis para lote con 10% de plantas enfermas ===
resultados_10 <- ejecutar_simulacion(0.10)
## Análisis para tamaño de muestra n = 5

## Valor p de la prueba de Shapiro-Wilk: 1.296402e-29 
## 
## Análisis para tamaño de muestra n = 10

## Valor p de la prueba de Shapiro-Wilk: 6.926359e-24 
## 
## Análisis para tamaño de muestra n = 15

## Valor p de la prueba de Shapiro-Wilk: 3.778667e-18 
## 
## Análisis para tamaño de muestra n = 20

## Valor p de la prueba de Shapiro-Wilk: 2.048936e-16 
## 
## Análisis para tamaño de muestra n = 30

## Valor p de la prueba de Shapiro-Wilk: 3.249148e-11 
## 
## Análisis para tamaño de muestra n = 50

## Valor p de la prueba de Shapiro-Wilk: 2.307379e-09 
## 
## Análisis para tamaño de muestra n = 60

## Valor p de la prueba de Shapiro-Wilk: 1.599034e-07 
## 
## Análisis para tamaño de muestra n = 100

## Valor p de la prueba de Shapiro-Wilk: 0.0001187274 
## 
## Análisis para tamaño de muestra n = 200

## Valor p de la prueba de Shapiro-Wilk: 0.003415887 
## 
## Análisis para tamaño de muestra n = 500

## Valor p de la prueba de Shapiro-Wilk: 0.08313573
# Escenario 2: Lote con 90% de plantas enfermas
cat("=== Análisis para lote con 90% de plantas enfermas ===\n")
## === Análisis para lote con 90% de plantas enfermas ===
resultados_90 <- ejecutar_simulacion(0.90)
## Análisis para tamaño de muestra n = 5

## Valor p de la prueba de Shapiro-Wilk: 1.296402e-29 
## 
## Análisis para tamaño de muestra n = 10

## Valor p de la prueba de Shapiro-Wilk: 6.926359e-24 
## 
## Análisis para tamaño de muestra n = 15

## Valor p de la prueba de Shapiro-Wilk: 3.778667e-18 
## 
## Análisis para tamaño de muestra n = 20

## Valor p de la prueba de Shapiro-Wilk: 2.048936e-16 
## 
## Análisis para tamaño de muestra n = 30

## Valor p de la prueba de Shapiro-Wilk: 3.249148e-11 
## 
## Análisis para tamaño de muestra n = 50

## Valor p de la prueba de Shapiro-Wilk: 2.307379e-09 
## 
## Análisis para tamaño de muestra n = 60

## Valor p de la prueba de Shapiro-Wilk: 1.599034e-07 
## 
## Análisis para tamaño de muestra n = 100

## Valor p de la prueba de Shapiro-Wilk: 0.0001187274 
## 
## Análisis para tamaño de muestra n = 200

## Valor p de la prueba de Shapiro-Wilk: 0.003415887 
## 
## Análisis para tamaño de muestra n = 500

## Valor p de la prueba de Shapiro-Wilk: 0.08313573

Conclusión

  • Estimación precisa de la proporción: La media de los estimadores de la proporción muestral, p^, calculada a partir de las 500 simulaciones, es cercana al valor verdadero de la proporción de plantas enfermas en la población (0.50). Esto indica que el método utilizado para estimar la proporción muestral es insesgado, ya que la media de los estimadores se aproxima al valor verdadero a lo largo de las simulaciones.

  • Dispersión de los estimadores: La desviación estándar de los estimadores muestra que, aunque hay cierta variabilidad en los resultados, los estimadores son razonablemente precisos. A medida que aumenta el tamaño de la muestra, esta variabilidad disminuye, lo que indica que las estimaciones se vuelven más precisas con muestras más grandes, lo que es consistente con la teoría del muestreo.

  • Asimetría y curtosis: os valores calculados de asimetría y curtosis muestran que la distribución de los estimadores no es perfectamente simétrica ni tiene la curtosis de una distribución normal estándar. Sin embargo, los valores son relativamente pequeños, lo que indica que las desviaciones de la normalidad no son extremas.

  • Evaluación de normalidad: La prueba de Shapiro-Wilk indica que la normalidad de los estimadores mejora a medida que el tamaño de la muestra aumenta. Para tamaños de muestra pequeños (como 5, 10 o 15), la prueba rechaza la hipótesis de normalidad, mientras que para muestras más grandes (especialmente de 200 o 500), la distribución de los estimadores se aproxima más a la normalidad, con valores de p mayores en la prueba de Shapiro-Wilk.

  • Los gráficos Q-Q corroboran estos resultados, mostrando cómo los estimadores se alinean más con los cuantiles teóricos de una distribución normal conforme el tamaño de la muestra aumenta.

  • En resumen, el estimador de proporción muestral es una herramienta robusta para estimar la proporción real de una población, con mayor precisión y normalidad a medida que se incrementa el tamaño de la muestra.

Problema 4

# Datos de la muestra
muestra <- c(7.69, 4.97, 4.56, 6.49, 4.34, 6.24, 4.45)

# Tamaño de la muestra y número de réplicas bootstrap
n <- length(muestra)
k <- 1000

# Inicializar vector para almacenar las medias bootstrap
medias_bootstrap <- numeric(k)

# Semilla para reproducibilidad
set.seed(123)

# Generar muestras bootstrap y calcular medias
for (i in 1:k) {
  muestra_bootstrap <- sample(muestra, size = n, replace = TRUE)
  medias_bootstrap[i] <- mean(muestra_bootstrap)
}

# Calcular los percentiles 2.5 y 97.5
p2.5 <- quantile(medias_bootstrap, 0.025)
p97.5 <- quantile(medias_bootstrap, 0.975)

# Intervalo de confianza con el Método 1
intervalo_metodo_1 <- c(p2.5, p97.5)

# Media muestral original
media_muestral <- mean(muestra)

# Intervalo de confianza ajustado con el Método 2
intervalo_metodo_2 <- c(2 * media_muestral - p97.5, 2 * media_muestral - p2.5)

# Mostrar resultados
cat("Intervalo de confianza (Método 1):", intervalo_metodo_1, "\n")
## Intervalo de confianza (Método 1): 4.748393 6.508643
cat("Intervalo de confianza (Método 2):", intervalo_metodo_2, "\n")
## Intervalo de confianza (Método 2): 4.559929 6.320179

Interpretación

  • Método 1: Este intervalo de confianza se construye tomando los percentiles de las medias de las muestras bootstrap. Proporciona un intervalo que refleja la variabilidad observada a partir del muestreo con reemplazo, sin hacer supuestos sobre la distribución subyacente de los datos.

  • Método 2: Este método ajusta el intervalo calculado en el Método 1 utilizando la media muestral original. Este enfoque intenta corregir cualquier sesgo que pudiera existir en la estimación del intervalo de confianza.

  • Confiabilidad: Dado que los métodos de bootstrap no dependen de suposiciones sobre la distribución subyacente de los datos, pueden ser particularmente útiles cuando se trabaja con muestras pequeñas y desconocemos la distribución de los datos.

Problema 5

Caso 1 original:

#Gráfica 1
#install.packages("plotly")

library(reshape2)

#Se necesita el paquete pwr 
if(!require(pwr)){install.packages("pwr");library("pwr")}
## Cargando paquete requerido: pwr
# t-TEST
# Se aplicará power.t.test del paquete stats (ya en R). Calcula la potencia de la prueba t de una o dos muestras, o determina los parámetros para obtener un valor particular de la potencia.

d<-seq(.1,2,by=.1) # 20 tamaños de los efectos
n<-1:150 # Tamaños muestrales

t.test.power.effect <-as.data.frame(do.call("cbind",lapply(1:length(d),function(i)
{
  sapply(1:length(n),function(j)
  {
    power.t.test(n=n[j],d=d[i],sig.level=0.05,power=NULL,type= "two.sample")$power
  })
})))

# Si algunas potencias no se pueden calcular, se ajustan a cero:
t.test.power.effect[is.na(t.test.power.effect)] <- 0 
colnames(t.test.power.effect)<-paste (d,"effect size")

#Graficando los resultados

prueba <-t.test.power.effect #data frame de 150 X 20 (para graficar)
cuts_num<-c(2,5,8) # cortes

#Cortes basados en: Cohen, J. (1988). Statistical Power Analysis for the Behavioral Sciences (2nd ed.). Hillsdale, NJ: Lawrence Erlbaum Associates, Publishers.
cuts_cat<-c("pequenio","medio","grande") 

columnas <- 1:ncol(prueba) #Lista de los valores 1:20
color_linea<-rainbow(length(columnas), alpha=.5) # Lista de 20 colores
grosor_linea=3 # Grosor de la línea

#Para el tipo de línea: (“blank”, “solid”, “dashed”, “dotted”, “dotdash”, “longdash”, “twodash”) ó  (0, 1, 2, 3, 4, 5, 6). 
#Note que lty = “solid” is idéntica a lty=1.

tipo_linea <- rep(1,length(color_linea))        #Repetir length(color)=20 veces el 1
tipo_linea[cuts_num]<-c(2:(length(cuts_num)+1)) #Asignar 2, 3, 4 en las posiciones 2, 5, 8 de tipo_linea

#Resaltar posiciones importantes
cuts_num<-c(2,5,8) # Cortes

#Cortes basados en: Cohen, J. (1988). Statistical Power Analysis for the Behavioral Sciences (2nd ed.). Hillsdale, NJ: Lawrence Erlbaum Associates, Publishers.
cuts_cat<-c("pequenio","medio","grande") 
color_linea[cuts_num]<-c("black")

efecto <- d # Listado de los 20 valores de 20
efecto[cuts_num] <- cuts_cat  #Reemplazar en "efecto" las posiciones cuts_num (2, 5, 8) por las categorías de cuts_cat

par(fig=c(0,.8,0,1),new=TRUE)
## Warning in par(fig = c(0, 0.8, 0, 1), new = TRUE): llamada par(new=TRUE) sin
## gráfico
#Gráfica
plot(1, type="n", #no produce puntos ni líneas
     frame.plot=FALSE, 
     xlab="Tamanio muestral", ylab="Potencia", 
     xlim=c(1,150),  ylim=c(0,1), 
     main="t-Test", axes = FALSE)

#Editando los ejes, grid, etc.
abline(v=seq(0,150,by=10), col = "lightgray", lty = "dotted") # Grid vertical
abline(h=seq(0,1,by=.05), col = "lightgray", lty = "dotted")  # Grid horizontal
axis(1,seq(0,150,by=10)) # Números en eje X
axis(2,seq(0,1,by=.05))  # Números en eje Y

#Plot de las lineas 
#columnas <- 1:ncol(prueba) # lista de los valores 1:20
for(i in 1:length(columnas)) #length(columnas)=20
{
  lines(1:150,
        #prueba (data frame de 150 X 20, para graficar)
        #columna <- 1:ncol(prueba) listado de valores 1:20 
        prueba[,columnas[i]], #filtrar "prueba" para valor de columna
        col=color_linea[i],   #color_linea[cuts_num]<-c("black")
        lwd=grosor_linea,     #grosor de cada linea
        lty=tipo_linea[i]     #tipo_linea[cuts_num]<-c(2:(length(cuts_num)+1))
  )
}

#Leyendas
par(fig=c(.65,1,0,1),new=TRUE)
plot.new()
legend("top",legend=efecto, col=color_linea, lwd=3, lty=tipo_linea, title="Tamanio efecto", 
       bty="n" #Opciones: o (complete box), n (no box), 7, L, C, U
)

#Gráfica 2

#plot using ggplot2

#library(ggplot2)
#library(reshape)
#library(plotly)

obj <- cbind(size=1:150, prueba) #Agregando el tamanio al data frame "prueba" 

# Usar melt y unir con "effect" para el mapeo
#El data frame "obj" se reconstruye con respecto al parámetro id="size". 
melted <- cbind(reshape::melt(obj, id="size"), effect=rep(d,each=150)) 

p<- ggplot(data=melted, aes(x=size, y=value, color=as.factor(effect))) + 
  geom_line(size=0.7,alpha=.5) +
  ylab("Potencia") + 
  xlab("Tamaño muestral") + 
  ggtitle("t-Test")+
  theme_bw() +
  #guides(fill=guide_legend(title="Efecto"))
  #scale_fill_discrete(name = "Efecto")
  #labs(fill='Efecto') 
  #scale_fill_manual("Efecto"#,values=c("orange","red")
  scale_color_discrete(name = "Tamaño del efecto")    

# Interactive plot
plotly::ggplotly(p)

Caso 1: Variando los tamaños de los efectos (d)

# Cargar el paquete necesario
if (!require(pwr)) {
  install.packages("pwr")
  library(pwr)
}

# Definir los parámetros
d <- c(0.1, 0.2, 0.5, 0.8)  # Tamaños de efecto: muy pequeño, pequeño, mediano y grande
n <- c(20, 60, 100, 140)    # Tamaños muestrales a analizar
alpha_values <- c(0.01, 0.025, 0.05, 0.1, 0.2)  # Niveles de significancia

# Inicializar un data frame para almacenar los resultados
resultados <- data.frame()

# Iterar sobre los niveles de significancia
for (alpha in alpha_values) {
  # Iterar sobre los tamaños del efecto
  for (efecto in d) {
    # Iterar sobre los tamaños muestrales
    for (tamano in n) {
      # Calcular la potencia con el valor de alpha, el tamaño del efecto y la muestra
      potencia <- power.t.test(n=tamano, d=efecto, sig.level=alpha, type="two.sample")$power
      
      # Guardar los resultados en el data frame
      resultados <- rbind(resultados, 
                          data.frame(Tamano_Muestral=tamano, 
                                     Tamano_Efecto=efecto, 
                                     Alpha=alpha, 
                                     Potencia=potencia))
    }
  }
}

# Visualizar los resultados
library(ggplot2)

# Gráfica para comparar los niveles de significancia
ggplot(resultados, aes(x=Tamano_Muestral, y=Potencia, color=factor(Alpha))) +
  geom_line(size=1) +
  facet_wrap(~Tamano_Efecto, scales = "free_y") +
  labs(title="Relacion entre tamanio muestral y potencia para diferentes niveles de significancia",
       x="Tamanio muestral",
       y="Potencia",
       color="Nivel de significancia") +
  theme_minimal()

Análisis:

  • Con un nivel de significancia más alto: La potencia aumentará, es decir, será más fácil rechazar la hipótesis nula. Sin embargo, esto aumenta el riesgo de errores de tipo I.

  • Con un nivel de significancia más bajo (ej. 0.01): La potencia será menor para el mismo tamaño muestral, ya que es más difícil rechazar la hipótesis nula con un criterio más estricto.

*Efecto del tamaño muestral: Para tamaños de muestra más grandes (100 y 140), la potencia será mayor en general, lo que reduce la probabilidad de errores de tipo II.

  • Para muestras pequeñas como 20, la potencia será baja a menos que el tamaño del efecto sea grande.

  • Para muestras de 100 y 140, es probable que la potencia sea suficientemente alta para detectar tamaños de efecto incluso pequeños, especialmente si el nivel de significancia es más alto (ej. 0.05 o 0.1).

Caso 2 original:

library(dplyr)    
## 
## Adjuntando el paquete: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(tidyr)    #Para manipulación de datos: separate, gather, spread
## 
## Adjuntando el paquete: 'tidyr'
## The following object is masked from 'package:reshape2':
## 
##     smiths
library(ggplot2)  
library(plotly)   #Para curvas de potencias interactivas
## 
## Adjuntando el paquete: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
library(pwr)      #Para cálculo de las potencias

#Generar cálculos de las potencias con la funcion pwr.t2n.test.
#Es un t-test para 2 muestras con tamaños diferentes 
#Aquí: d es el tamaño del efecto, Power= potencia de la prueba= 1-beta): 

#pwr.t2n.test(n1 = NULL, n2= NULL, d = NULL, sig.level = 0.05, power = NULL,  alternative = c("two.sided",  "less","greater"))

ptab <- cbind(NULL, NULL)       

for (i in seq(0,1, length.out = 200)){
  pwrt1 <- pwr.t2n.test(n1 = 28, n2 = 1406, 
                        sig.level = 0.05, power = NULL, 
                        d = i, alternative="two.sided")
  pwrt2 <- pwr.t2n.test(n1 = 144, n2 = 1290, 
                        sig.level = 0.05, power = NULL, 
                        d = i, alternative="two.sided")
  pwrt3 <- pwr.t2n.test(n1 = 287, n2 = 1147, 
                        sig.level = 0.05, power = NULL, 
                        d = i, alternative="two.sided")
  pwrt4 <- pwr.t2n.test(n1 = 430, n2 = 1004, 
                        sig.level = 0.05, power = NULL, 
                        d = i, alternative="two.sided")
  pwrt5 <- pwr.t2n.test(n1 = 574, n2 = 860, 
                        sig.level = 0.05, power = NULL, 
                        d = i, alternative="two.sided")
  pwrt6 <- pwr.t2n.test(n1 = 717, n2 = 717, 
                        sig.level = 0.05, power = NULL, 
                        d = i, alternative="two.sided")
  
  #Es un data frame de tamaño 200 por 12: 
  ptab <- rbind(ptab, cbind(pwrt1$d, pwrt1$power,
                            pwrt2$d, pwrt2$power,
                            pwrt3$d, pwrt3$power,
                            pwrt4$d, pwrt4$power,
                            pwrt5$d, pwrt5$power,
                            pwrt6$d, pwrt6$power))
}


#Es un data frame de tamaño 200 por 13 (la 1ra columna es ID)
ptab <- cbind(seq_len(nrow(ptab)), ptab)

colnames(ptab) <- c("id","n1=28, n2=1406;effect size","n1=28, n2=1406;power",
                    "n1=144, n2=1290;effect size","n1=144, n2=1290;power",
                    "n1=287, n2=1147;effect size","n1=287, n2=1147;power",
                    "n1=430, n2=1004;effect size","n1=430, n2=1004;power",
                    "n1=574, n2=860;effect size","n1=574, n2=860;power",
                    "n1=717, n2=717;effect size","n1=717, n2=717;power")

#gather se usa  para "reunir" un par key-value. En este caso, en 3 columnas: ID, variables y respuestas numericas
temp1 <- ptab %>% as.data.frame() %>%   gather(key = name, value = val, 2:13)

#Separar celdas en columnas, de acuerdo a una condición (sep=). En este caso, se separó "name" en dos columnas: samples y pruebas 
temp2 <- temp1 %>%   separate(col = name, into = c("samples", "pruebas"), sep = ";")


#La función spread hace lo opuesto a gather. Son funciones complementarias. 
#Es decir, si al resultado de aplicar la función spread le aplicamos la función gather llegamos al dataset original.
temp3 <- temp2 %>%   spread(key = pruebas, value = val)

#Convertir la variable "samples" a factor.
temp3$samples <- factor(temp3$samples, 
                        levels = c("n1=28, n2=1406", "n1=144, n2=1290", 
                                   "n1=287, n2=1147", "n1=430, n2=1004",
                                   "n1=574, n2=860", "n1=717, n2=717")
)

#Gráfica
p<- ggplot(temp3, aes(x = `effect size`, y = power, color = samples)) +
  geom_line(size=1) + 
  
  theme_bw() + 
  theme(axis.text=element_text(size=10), 
        axis.title=element_text(size=10), 
        legend.text=element_text(size=10)) +
  
  geom_vline(xintercept = .54, linetype = 2) +
  geom_hline(yintercept = 0.80, linetype = 2)+
  
  labs(x="Effect size", y="Power") +
  scale_color_discrete(name = "Sampling size") 

# so simple to make interactive plots
plotly::ggplotly(p)

Caso 2: Variando los tamaños muestrales

# Cargar los paquetes necesarios
if (!require(ggplot2)) {
  install.packages("ggplot2")
  library(ggplot2)
}
if (!require(pwr)) {
  install.packages("pwr")
  library(pwr)
}

# Definir los tamaños de muestra y los niveles de significancia
n1_values <- c(28, 144, 287, 430, 574, 717)  # n1 valores
n2_values <- c(1406, 1290, 1147, 1004, 860, 717)  # n2 valores
alpha_values <- c(0.01, 0.025, 0.05, 0.1, 0.2)  # Niveles de significancia
d_values <- seq(0.1, 1.5, length.out = 100)  # Tamaños de efecto

# Inicializar un data frame para almacenar los resultados
resultados <- data.frame()

# Realizar los cálculos para cada combinación de n1, n2 y alpha
for (alpha in alpha_values) {
  for (i in 1:length(n1_values)) {
    n1 <- n1_values[i]
    n2 <- n2_values[i]
    
    for (d in d_values) {
      # Calcular la potencia usando pwr.t2n.test
      potencia <- pwr.t2n.test(n1=n1, n2=n2, d=d, sig.level=alpha, power=NULL, alternative="two.sided")$power
      
      # Guardar los resultados en el data frame
      resultados <- rbind(resultados, 
                          data.frame(Tamano_Efecto=d, Tamano_Muestral=paste("n1=", n1, ", n2=", n2), Alpha=alpha, Potencia=potencia))
    }
  }
}

# Gráfica de comparacion
plot <- ggplot(resultados, aes(x=Tamano_Efecto, y=Potencia, color=Tamano_Muestral, linetype=factor(Alpha))) +
  geom_line(size=1) +
  labs(title="Relacion entre el tamanio del efecto y la potencia para diferentes niveles de significancia",
       x="Tamanio del efecto",
       y="Potencia",
       color="Tamanio de muestra",
       linetype="Nivel de significancia") +
  geom_hline(yintercept = 0.80, linetype = "dashed") +  # Línea de potencia 80%
  geom_vline(xintercept = 0.54, linetype = "dashed", color = "red") +  # Línea de tamaño del efecto 0.54
  theme_minimal()

# Forzar la impresión del gráfico
print(plot)

Análisis:

  • Impacto del nivel de significancia: Con un nivel de significancia más bajo), se espera que la potencia sea menor para los mismos tamaños del efecto en comparación con \(\alpha = 0.05\) o \(\alpha = 0.1\).Esto significa que será más difícil rechazar la hipótesis nula bajo un criterio más estricto.

  • Desbalance en los tamaños muestrales (n1 = 28 y n2 = 1406): En esta configuración, es probable que sea necesario un tamaño del efecto mucho mayor para alcanzar una potencia aceptable (80 %) debido al desbalance significativo entre los tamaños de los grupos. Según el análisis original, se necesita un tamaño del efecto mayor a 0.54 para alcanzar un poder del 80 % con estos tamaños de muestra desiguales.