PROBLEMA 1

Introducción

En este ejercicio, se estima el valor de \(\pi\). Se genera un número aleatorio de puntos en un cuadrado y se calcula la proporción de puntos que caen dentro de un círculo inscrito en el cuadrado. A partir de esta proporción, se estima el valor de \(\pi\).

Código y Resultados

knitr::opts_chunk$set(echo = TRUE)

# Número de puntos a generar
n <- 10000  # Cambia a 100000 para una estimación más precisa

# Generar coordenadas X e Y aleatorias
set.seed(42)  # Para reproducibilidad
X <- runif(n)
Y <- runif(n)

# Centro del círculo y radio al cuadrado
center_x <- 0.5
center_y <- 0.5
radius_squared <- 0.25

# Determinar si los puntos están dentro del círculo
inside_circle <- (X - center_x)^2 + (Y - center_y)^2 < radius_squared
num_inside_circle <- sum(inside_circle)

# Estimación de π
pi_estimate <- 4 * (num_inside_circle / n)

# Resultados
list(
  num_inside_circle = num_inside_circle,
  pi_estimate = pi_estimate
)
## $num_inside_circle
## [1] 7818
## 
## $pi_estimate
## [1] 3.1272
# Crear un data frame para visualización
data <- data.frame(X = X, Y = Y, InsideCircle = inside_circle)

# Gráfico de los puntos
library(ggplot2)
ggplot(data, aes(x = X, y = Y, color = InsideCircle)) +
  geom_point(alpha = 0.5) +
  scale_color_manual(values = c("red", "blue"), labels = c("Fuera del círculo", "Dentro del círculo")) +
  labs(title = "Puntos Aleatorios en el Cuadrado",
       subtitle = paste("Estimación de π:", round(pi_estimate, 4)),
       x = "X",
       y = "Y",
       color = "Categoría") +
  theme_minimal() +
  theme(legend.position = "bottom")

Explicación Detallada del Ejercicio

  1. Generar Puntos Aleatorios:
    • x <- runif(n) y y <- runif(n) generan n coordenadas aleatorias uniformemente distribuidas en el intervalo [0, 1] para las coordenadas x y y, respectivamente.
  2. Cálculo de Distancias:
    • La distancia al centro del círculo (0.5, 0.5) se calcula usando la fórmula (x - 0.5)^2 + (y - 0.5)^2. Esta es la distancia cuadrada desde el punto (x, y) al centro. La variable inside_circle almacena una serie de valores lógicos (TRUE o FALSE) indicando si el punto está dentro del círculo.
  3. Estimación de \(\pi\):
    • La proporción de puntos dentro del círculo sobre el total n es una estimación de \(\frac{\pi}{4}\). Multiplicando esta proporción por 4 da una estimación de \(\pi\).
  4. Visualización:
    • Utilizamos ggplot2 para crear un gráfico que muestra los puntos dentro y fuera del círculo. Los puntos dentro del círculo se colorean de manera diferente a los puntos fuera del círculo para una visualización clara.

Ajustar el número de puntos (n) en el código permitirá observar cómo varía la estimación con diferentes tamaños de muestra.

Este código permite generar una estimación de \(\pi\) con una cantidad significativa de puntos y visualizar la distribución de estos puntos en el cuadrado y el círculo.

PROBLEMA 2

Introducción

En este ejercicio, se analizarán las propiedades de varios estimadores para una muestra de una distribución exponencial con un parámetro \(\theta\) desconocido. Evaluaremos la insesgadez, eficiencia y consistencia de cada estimador.

Código y Resultados

knitr::opts_chunk$set(echo = TRUE)

# Parámetro theta
theta <- 2  # Supongamos un valor para el parámetro theta

# Función para calcular los estimadores
calc_estimators <- function(x) {
  theta1 <- (x[1] + x[2]) / 6 + (x[3] + x[4]) / 3
  theta2 <- (x[1] + 2*x[2] + 3*x[3] + 4*x[4]) / 5
  theta3 <- mean(x)
  theta4 <- (min(x) + max(x)) / 2
  return(c(theta1 = theta1, theta2 = theta2, theta3 = theta3, theta4 = theta4))
}

# Función para realizar la simulación
simulate_estimators <- function(n, theta, num_simulations = 1000) {
  estimators <- replicate(num_simulations, {
    sample <- rexp(4, rate = 1/theta)
    calc_estimators(sample)
  })
  return(t(estimators))
}

# Tamaños de muestra
sample_sizes <- c(20, 50, 100, 1000)


# Ejecutar simulaciones y almacenar resultados
results <- lapply(sample_sizes, function(n) simulate_estimators(n, theta))

# Función para calcular medias y varianzas
analyze_estimators <- function(estimates, theta) {
  means <- colMeans(estimates)
  variances <- apply(estimates, 2, var)
  bias <- means - theta
  list(means = means, variances = variances, bias = bias)
}

# Análisis de cada tamaño de muestra
analysis_results <- lapply(results, analyze_estimators, theta = theta)

# Resultados para n = 20
analysis_results[[1]]
## $means
##   theta1   theta2   theta3   theta4 
## 2.025120 4.025994 2.025619 2.391544 
## 
## $variances
##   theta1   theta2   theta3   theta4 
## 1.148173 4.838096 1.027036 1.673902 
## 
## $bias
##     theta1     theta2     theta3     theta4 
## 0.02512016 2.02599355 0.02561868 0.39154414

Explicación del Código

  1. Generación de Muestras:
    • rexp(4, rate = 1/theta) genera una muestra de tamaño 4 de una distribución exponencial con parámetro \(\theta\).
  2. Cálculo de Estimadores:
    • calc_estimators calcula los estimadores para una muestra dada.
  3. Simulación:
    • simulate_estimators genera múltiples muestras y calcula los estimadores para cada una. Se repite el proceso para cada tamaño de muestra especificado.
  4. Análisis:
    • analyze_estimators calcula la media, varianza y sesgo de los estimadores para evaluar insesgadez, eficiencia y consistencia.

Este código permite analizar y comparar las propiedades de los estimadores propuestos en función de diferentes tamaños de muestra, por lo cual ayuda a entender mejor sus características estadísticas.

PROBLEMA 3

Introducción

El Teorema del Límite Central (TLC) establece que la distribución de la proporción muestral se aproxima a una distribución normal a medida que el tamaño de la muestra aumenta. En este ejercicio, se analizará cómo se comporta esta aproximación para diferentes tamaños de muestra y porcentajes de individuos enfermos en la población.

a.Realice una simulación en la cual genere una población de n=1000 (Lote), donde el porcentaje de individuos (supongamos plantas) enfermas sea del 50%.

set.seed(123)
poblacion <- rbinom(n = 1000, size = 1, prob = 0.5)
numero_enfermas <- sum(poblacion)
numero_sanas <- 1000 - numero_enfermas

cat("Número de plantas enfermas:", numero_enfermas, "\n")
## Número de plantas enfermas: 493
cat("Número de plantas sanas:", numero_sanas, "\n")
## Número de plantas sanas: 507
  1. Genere una función que permita obtener una muestra aleatoria de la población y calcule el estimador de la proporción muestral para un tamaño de muestra dado n.
set.seed(123)
obtener_proporcion_muestral <- function(tamano_muestra, poblacion) {
  muestra <- sample(poblacion, size = tamano_muestra, replace = TRUE)
  p_hat <- mean(muestra)
  return(p_hat)
}
# Tamaño de la muestra n
tamano_muestra <- 300

#Proporción muestral
proporcion_muestral <- obtener_proporcion_muestral(tamano_muestra, poblacion)
print(paste("Para una muestra de tamaño:", tamano_muestra, "el estimador de la proporción es =", proporcion_muestral))
## [1] "Para una muestra de tamaño: 300 el estimador de la proporción es = 0.48"
  1. Repita el escenario anterior (b) n=500 veces y analice los resultados en cuanto al comportamiento de los 500 resultados del estimador pˆ¿Qué tan simétricos o sesgados son los resultados obtenidos? y ¿qué se puede observar en cuanto a la variabilidad?. Realice en su informe un comentario sobre los resultados obtenidos.
library("DescTools")
set.seed(123)
library(e1071)
#Repetimos el escenario 500 veces
proporciones <- replicate(500, obtener_proporcion_muestral(tamano_muestra = 300, poblacion))

# Análisis de los resultados
hist(proporciones, main = "Densidad de la proporción muestral", xlab = "Proporción muestral",freq=FALSE)
densidad <- density(proporciones)
lines(densidad, col = "blue")
abline(v = mean(proporciones), col = "red", lwd = 2)

Tabla=data.frame("ID"=0,"Tamaño_muestra"=tamano_muestra, "Media"=mean(proporciones),
                 "Mediana"=median(proporciones),"Moda"=Mode(proporciones),
                 "Desvest"=sd(proporciones),"Varianza"=var(proporciones),
                 "Asimetría"=skewness(proporciones))
Tabla
shapiro_test <- shapiro.test(proporciones)
shapiro_test
## 
##  Shapiro-Wilk normality test
## 
## data:  proporciones
## W = 0.99608, p-value = 0.2539
qqnorm(proporciones, main = "QQ-plot")
qqline(proporciones, col = "red")

Análisis punto C:

En primer lugar, con los resultados del data frame generado, se observó que la media, la moda y la mediana son muy cercanas entre sí (casi iguales), lo que nos indica que existe una gran probabilidad de que la distribución de los datos es simetrica.

Con respecto al Coeficiente de Asimetría, el valor es muy cercano a 0, lo que sugiere que la distribución no se encuentra muy lejos de ser simetrica

De las pruebas de bondad y ajuste, se resalta que el estadistico de prueba es cercano a 1, lo que nos indica que los datos tienen una distribución que se asemeja estrechamente a una distribución normal. Por su parte, al analizar el valor P (donde la Hipotesis Nula (H0) establece que los datos se distribuyen normalmente), el resultado nos arroja un valor P > 0.05; lo que nos indica que no hay suficientes evidencias para rechazar la hipotesis nula. Por lo anterior, hay una gran probabilidad que los datos se encuentran normalmente distribuidos.

Ahora bien, del grafico Q-Q, observamos que los puntos se alinean estrechamente con la línea de referencia, proporcionando evidencia visual de que los datos pueden ser normalmente distribuidos.

La desviación estándar de 0.029 es relativamente pequeña en comparación con la media de 0.49. Esto indica que la mayoría de los valores en el conjunto de datos están bastante cerca de la media. En otras palabras, hay una variabilidad baja en los datos y los valores no están muy dispersos.

  1. Realice los ejercicios completos b y c para tamaños de muestra n=5, 10, 15, 20, 30, 50, 60, 100, 200, 500. Y compare los resultados de los estimadores en cuanto a la normalidad. Investigue y utilice pruebas de bondad y ajuste (shapiro wilks) y métodos gráficos (grafico qq de normalidad)
set.seed(123)
# Tamaños de muestra
tamanos_muestra <- c(5, 10, 15, 20, 30, 50, 60, 100, 200, 500)

# Realizar simulaciones y análisis para cada tamaño de muestra
for (n in tamanos_muestra) {
  proporciones_n <- replicate(1000, obtener_proporcion_muestral(n, poblacion))
  
  # Configurar el área de gráficos para mostrar tres gráficos en una matriz de 2x2
  par(mfrow=c(2, 2))
  
  # Histograma de las proporciones muestrales
  hist(proporciones_n, main = paste("Histograma para n =", n), xlab = "Proporción muestral", freq = FALSE)
  lines(density(proporciones_n), col = "blue")
  
  # Gráfico QQ-Norm para las proporciones muestrales
  qqnorm(proporciones_n, main = paste("QQ-plot para n =", n))
  qqline(proporciones_n, col = "red")
  
  # Dejar un espacio en blanco para los resultados de la prueba de Shapiro-Wilk
  plot.new()
  
  # Prueba de Shapiro-Wilk para las proporciones muestrales
  shapiro_test <- shapiro.test(proporciones_n)
  
  # Mostrar los resultados de la prueba de Shapiro-Wilk en el espacio en blanco
  text(0.5, 0.6, paste("Shapiro-Wilk Test para n =", n, 
                        "\nW = ", round(shapiro_test$statistic, 5),
                        "\np-value = ", round(shapiro_test$p.value, 5)), cex = 1.2)
  
  # Restablecer la configuración de gráficos a la configuración por defecto
  par(mfrow=c(1, 1))
  
  # Pausar en entornos interactivos para permitir la visualización de los gráficos
  if (interactive()) {
    readline(prompt = "Presiona [Enter] para continuar...")
  }
}

set.seed(123)
# Suponiendo que las funciones y variables necesarias ya están definidas

# Tamaños de muestra
tamanos_muestra <- c(5, 10, 15, 20, 30, 50, 60, 100, 200, 500)

# Inicializar un data.frame para almacenar los resultados
resultados <- data.frame(
  Tamano_Muestra = integer(),
  Media = numeric(),
  Mediana = numeric(),
  Desviacion_Estandar = numeric(),
  Varianza = numeric(),
  Asimetria = numeric(),
  stringsAsFactors = FALSE
)

for (n in tamanos_muestra) {
  proporciones_n <- replicate(500, obtener_proporcion_muestral(n, poblacion))
  
  # Calcular estadísticas descriptivas
  media <- mean(proporciones_n)
  mediana <- median(proporciones_n)
  desv_estandar <- sd(proporciones_n)
  varianza <- var(proporciones_n)
  asimetria <- skewness(proporciones_n)
  
  # Agregar los resultados al data.frame
  resultados <- rbind(resultados, data.frame(
    Tamano_Muestra = n,
    Media = media,
    Mediana = mediana,
    Desviacion_Estandar = desv_estandar,
    Varianza = varianza,
    Asimetria = asimetria
  ))
}

# Mostrar la tabla de resultados
print(resultados)
##    Tamano_Muestra     Media   Mediana Desviacion_Estandar     Varianza
## 1               5 0.4996000 0.4000000          0.22667732 0.0513826052
## 2              10 0.4938000 0.5000000          0.15947547 0.0254324248
## 3              15 0.4876000 0.4666667          0.12883881 0.0165994389
## 4              20 0.4996000 0.5000000          0.11245057 0.0126451303
## 5              30 0.4968667 0.5000000          0.09312342 0.0086719706
## 6              50 0.4980400 0.5000000          0.07211768 0.0052009603
## 7              60 0.4978667 0.5000000          0.06511132 0.0042394834
## 8             100 0.4960200 0.4900000          0.04651329 0.0021634866
## 9             200 0.4919300 0.4900000          0.03454144 0.0011931113
## 10            500 0.4949920 0.4940000          0.02112032 0.0004460681
##      Asimetria
## 1   0.16154307
## 2   0.04791485
## 3   0.02325124
## 4  -0.06600513
## 5   0.02381993
## 6  -0.08987917
## 7   0.02171997
## 8   0.15185606
## 9  -0.09186055
## 10 -0.04466524

Analisis punto D:

De las pruebas de bondad y ajustes, evidenciamos que para el estadistico de prueba (valor W) a mayor número de muestras, se obtiene un valor más cercano a 1; lo que nos indica que los datos tienen una distribución que se asemeja estrechamente a una distribución normal.

Por otro lado, observamos que a mayor número de muestras, el valor P incrementa hasta superar el nivel de significancia elegido , identificando que para una muestra n=500 el valor P > 0.05, lo que nos permitió concluir que no existen suficientes evidencias para rechazar la Hipoteis Nula, por lo anterior, se considera que los datos se encuentran normalmente distribuidos.

Al analizar la gráfica de dispersión QQ, evidenciamos que a mayor número de muestras, los puntos se comienzan a alinear estrechamente con la línea de referencia; lo que sugiere que los datos pueden ser normalmente distribuidos.

Al calcular los estimadores ^P con diferentes tamaños de muestra, notamos que ^P tiende a aproximarse al valor del parámetro P = 0.5. Por otro lado, se observa que la varianza disminuye a medida que el tamaño de la muestra se incrementa, lo cual respalda la validez del Teorema del Limite Central

e.1. Repita toda la simulación (puntos a – d), pero ahora para lotes con 10% de plantas enfermas

a.Realice una simulación en la cual genere una población de n=1000 (Lote), donde el porcentaje de individuos (supongamos plantas) enfermas sea del 10%.

poblacion <- rbinom(n = 1000, size = 1, prob = 0.1)
numero_enfermas <- sum(poblacion)
numero_sanas <- 1000 - numero_enfermas

cat("Número de plantas enfermas:", numero_enfermas, "\n")
## Número de plantas enfermas: 87
cat("Número de plantas sanas:", numero_sanas, "\n")
## Número de plantas sanas: 913
  1. Genere una función que permita:Obtener una muestra aleatoria de la población y Calcule el estimador de la proporción muestral pˆpara un tamaño de muestra dado n
obtener_proporcion_muestral <- function(tamano_muestra, poblacion) {
  muestra <- sample(poblacion, size = tamano_muestra, replace = TRUE)
  p_hat <- mean(muestra)
  return(p_hat)
}
# Tamaño de la muestra n
tamano_muestra <- 300

#Proporción muestral
proporcion_muestral <- obtener_proporcion_muestral(tamano_muestra, poblacion)
print(paste("Para una muestra de tamaño:", tamano_muestra, "el estimador de la proporción es =", proporcion_muestral))
## [1] "Para una muestra de tamaño: 300 el estimador de la proporción es = 0.0833333333333333"
  1. Repita el escenario anterior (b) n=500 veces y analice los resultados en cuanto al comportamiento de los 500 resultados del estimador p. ¿Qué tan simétricos o sesgados son los resultados obtenidos? y ¿qué se puede observar en cuanto a la variabilidad?. Realice en su informe un comentario sobre los resultados obtenidos.
library("DescTools")
set.seed(123)
library(e1071)
#Repetimos el escenario 500 veces
proporciones <- replicate(500, obtener_proporcion_muestral(tamano_muestra = 300, poblacion))

# Análisis de los resultados
hist(proporciones, main = "Densidad de la proporción muestral", xlab = "Proporción muestral",freq=FALSE)
densidad <- density(proporciones)
lines(densidad, col = "blue")
abline(v = mean(proporciones), col = "red", lwd = 2)

Tabla=data.frame("ID"=0,"Tamaño_muestra"=tamano_muestra, "Media"=mean(proporciones),
                 "Mediana"=median(proporciones),"Moda"=Mode(proporciones),
                 "Desvest"=sd(proporciones),"Varianza"=var(proporciones),
                 "Asimetría"=skewness(proporciones))
Tabla
shapiro_test <- shapiro.test(proporciones)
shapiro_test
## 
##  Shapiro-Wilk normality test
## 
## data:  proporciones
## W = 0.99393, p-value = 0.04322
qqnorm(proporciones, main = "QQ-plot")
qqline(proporciones, col = "red")

  1. Realice los ejercicios completos b y c para tamaños de muestra n=5, 10, 15, 20, 30, 50, 60, 100, 200, 500. Y compare los resultados de los estimadores en cuanto a la normalidad. Investigue y utilice pruebas de bondad y ajuste (shapiro wilks) y métodos gráficos (grafico qq de normalidad)
# Tamaños de muestra
tamanos_muestra <- c(5, 10, 15, 20, 30, 50, 60, 100, 200, 500)

# Realizar simulaciones y análisis para cada tamaño de muestra
for (n in tamanos_muestra) {
  proporciones_n <- replicate(1000, obtener_proporcion_muestral(n, poblacion))
  
  # Configurar el área de gráficos para mostrar tres gráficos en una matriz de 2x2
  par(mfrow=c(2, 2))
  
  # Histograma de las proporciones muestrales
  hist(proporciones_n, main = paste("Histograma para n =", n), xlab = "Proporción muestral", freq = FALSE)
  lines(density(proporciones_n), col = "blue")
  
  # Gráfico QQ-Norm para las proporciones muestrales
  qqnorm(proporciones_n, main = paste("QQ-plot para n =", n))
  qqline(proporciones_n, col = "red")
  
  # Dejar un espacio en blanco para los resultados de la prueba de Shapiro-Wilk
  plot.new()
  
  # Prueba de Shapiro-Wilk para las proporciones muestrales
  shapiro_test <- shapiro.test(proporciones_n)
  
  # Mostrar los resultados de la prueba de Shapiro-Wilk en el espacio en blanco
  text(0.5, 0.6, paste("Shapiro-Wilk Test para n =", n, 
                        "\nW = ", round(shapiro_test$statistic, 5),
                        "\np-value = ", round(shapiro_test$p.value, 5)), cex = 1.2)
  
  # Restablecer la configuración de gráficos a la configuración por defecto
  par(mfrow=c(1, 1))
  
  # Pausar en entornos interactivos para permitir la visualización de los gráficos
  if (interactive()) {
    readline(prompt = "Presiona [Enter] para continuar...")
  }
}

# Suponiendo que las funciones y variables necesarias ya están definidas

# Tamaños de muestra
tamanos_muestra <- c(5, 10, 15, 20, 30, 50, 60, 100, 200, 500)

# Inicializar un data.frame para almacenar los resultados
resultados <- data.frame(
  Tamano_Muestra = integer(),
  Media = numeric(),
  Mediana = numeric(),
  Desviacion_Estandar = numeric(),
  Varianza = numeric(),
  Asimetria = numeric(),
  stringsAsFactors = FALSE
)

for (n in tamanos_muestra) {
  proporciones_n <- replicate(500, obtener_proporcion_muestral(n, poblacion))
  
  # Calcular estadísticas descriptivas
  media <- mean(proporciones_n)
  mediana <- median(proporciones_n)
  desv_estandar <- sd(proporciones_n)
  varianza <- var(proporciones_n)
  asimetria <- skewness(proporciones_n)
  
  # Agregar los resultados al data.frame
  resultados <- rbind(resultados, data.frame(
    Tamano_Muestra = n,
    Media = media,
    Mediana = mediana,
    Desviacion_Estandar = desv_estandar,
    Varianza = varianza,
    Asimetria = asimetria
  ))
}

# Mostrar la tabla de resultados
print(resultados)
##    Tamano_Muestra      Media    Mediana Desviacion_Estandar     Varianza
## 1               5 0.08040000 0.00000000          0.13285760 0.0176511423
## 2              10 0.08940000 0.10000000          0.09101763 0.0082842084
## 3              15 0.08373333 0.06666667          0.06943049 0.0048205923
## 4              20 0.08250000 0.05000000          0.06210043 0.0038564629
## 5              30 0.08666667 0.06666667          0.05254596 0.0027610777
## 6              50 0.08440000 0.08000000          0.03702136 0.0013705812
## 7              60 0.08876667 0.08333333          0.03473876 0.0012067813
## 8             100 0.08816000 0.09000000          0.02902239 0.0008422990
## 9             200 0.08707000 0.08500000          0.01964528 0.0003859370
## 10            500 0.08874800 0.08800000          0.01278810 0.0001635356
##     Asimetria
## 1  1.99879069
## 2  1.02134322
## 3  0.83559764
## 4  0.72099127
## 5  0.54834333
## 6  0.38854996
## 7  0.42925928
## 8  0.34110210
## 9  0.13413697
## 10 0.06499163

e.2. Repita toda la simulación (puntos a – d), pero ahora para lotes con 90% de plantas enfermas

a.Realice una simulación en la cual genere una población de n=1000 (Lote), donde el porcentaje de individuos (supongamos plantas) enfermas sea del 90%.

poblacion <- rbinom(n = 1000, size = 1, prob = 0.9)
numero_enfermas <- sum(poblacion)
numero_sanas <- 1000 - numero_enfermas

cat("Número de plantas enfermas:", numero_enfermas, "\n")
## Número de plantas enfermas: 910
cat("Número de plantas sanas:", numero_sanas, "\n")
## Número de plantas sanas: 90
  1. Genere una función que permita:Obtener una muestra aleatoria de la población y Calcule el estimador de la proporción muestral pˆpara un tamaño de muestra dado n
obtener_proporcion_muestral <- function(tamano_muestra, poblacion) {
  muestra <- sample(poblacion, size = tamano_muestra, replace = TRUE)
  p_hat <- mean(muestra)
  return(p_hat)
}
# Tamaño de la muestra n
tamano_muestra <- 300

#Proporción muestral
proporcion_muestral <- obtener_proporcion_muestral(tamano_muestra, poblacion)
print(paste("Para una muestra de tamaño:", tamano_muestra, "el estimador de la proporción es =", proporcion_muestral))
## [1] "Para una muestra de tamaño: 300 el estimador de la proporción es = 0.903333333333333"
  1. Repita el escenario anterior (b) n=500 veces y analice los resultados en cuanto al comportamiento de los 500 resultados del estimador p. ¿Qué tan simétricos o sesgados son los resultados obtenidos? y ¿qué se puede observar en cuanto a la variabilidad?. Realice en su informe un comentario sobre los resultados obtenidos.
library("DescTools")
library(e1071)
#Repetimos el escenario 500 veces
proporciones <- replicate(500, obtener_proporcion_muestral(tamano_muestra = 300, poblacion))

# Análisis de los resultados
hist(proporciones, main = "Densidad de la proporción muestral", xlab = "Proporción muestral",freq=FALSE)
densidad <- density(proporciones)
lines(densidad, col = "blue")
abline(v = mean(proporciones), col = "red", lwd = 2)

Tabla=data.frame("ID"=0,"Tamaño_muestra"=tamano_muestra, "Media"=mean(proporciones),
                 "Mediana"=median(proporciones),"Moda"=Mode(proporciones),
                 "Desvest"=sd(proporciones),"Varianza"=var(proporciones),
                 "Asimetría"=skewness(proporciones))
Tabla
shapiro_test <- shapiro.test(proporciones)
shapiro_test
## 
##  Shapiro-Wilk normality test
## 
## data:  proporciones
## W = 0.99519, p-value = 0.1241
qqnorm(proporciones, main = "QQ-plot")
qqline(proporciones, col = "red")

  1. Realice los ejercicios completos b y c para tamaños de muestra n=5, 10, 15, 20, 30, 50, 60, 100, 200, 500. Y compare los resultados de los estimadores en cuanto a la normalidad. Investigue y utilice pruebas de bondad y ajuste (shapiro wilks) y métodos gráficos (grafico qq de normalidad)
# Tamaños de muestra
tamanos_muestra <- c(5, 10, 15, 20, 30, 50, 60, 100, 200, 500)

# Realizar simulaciones y análisis para cada tamaño de muestra
for (n in tamanos_muestra) {
  proporciones_n <- replicate(1000, obtener_proporcion_muestral(n, poblacion))
  
  # Configurar el área de gráficos para mostrar tres gráficos en una matriz de 2x2
  par(mfrow=c(2, 2))
  
  # Histograma de las proporciones muestrales
  hist(proporciones_n, main = paste("Histograma para n =", n), xlab = "Proporción muestral", freq = FALSE)
  lines(density(proporciones_n), col = "blue")
  
  # Gráfico QQ-Norm para las proporciones muestrales
  qqnorm(proporciones_n, main = paste("QQ-plot para n =", n))
  qqline(proporciones_n, col = "red")
  
  # Dejar un espacio en blanco para los resultados de la prueba de Shapiro-Wilk
  plot.new()
  
  # Prueba de Shapiro-Wilk para las proporciones muestrales
  shapiro_test <- shapiro.test(proporciones_n)
  
  # Mostrar los resultados de la prueba de Shapiro-Wilk en el espacio en blanco
  text(0.5, 0.6, paste("Shapiro-Wilk Test para n =", n, 
                        "\nW = ", round(shapiro_test$statistic, 5),
                        "\np-value = ", round(shapiro_test$p.value, 5)), cex = 1.2)
  
  # Restablecer la configuración de gráficos a la configuración por defecto
  par(mfrow=c(1, 1))
  
  # Pausar en entornos interactivos para permitir la visualización de los gráficos
  if (interactive()) {
    readline(prompt = "Presiona [Enter] para continuar...")
  }
}

# Suponiendo que las funciones y variables necesarias ya están definidas

# Tamaños de muestra
tamanos_muestra <- c(5, 10, 15, 20, 30, 50, 60, 100, 200, 500)

# Inicializar un data.frame para almacenar los resultados
resultados <- data.frame(
  Tamano_Muestra = integer(),
  Media = numeric(),
  Mediana = numeric(),
  Desviacion_Estandar = numeric(),
  Varianza = numeric(),
  Asimetria = numeric(),
  stringsAsFactors = FALSE
)

for (n in tamanos_muestra) {
  proporciones_n <- replicate(500, obtener_proporcion_muestral(n, poblacion))
  
  # Calcular estadísticas descriptivas
  media <- mean(proporciones_n)
  mediana <- median(proporciones_n)
  desv_estandar <- sd(proporciones_n)
  varianza <- var(proporciones_n)
  asimetria <- skewness(proporciones_n)
  
  # Agregar los resultados al data.frame
  resultados <- rbind(resultados, data.frame(
    Tamano_Muestra = n,
    Media = media,
    Mediana = mediana,
    Desviacion_Estandar = desv_estandar,
    Varianza = varianza,
    Asimetria = asimetria
  ))
}

# Mostrar la tabla de resultados
print(resultados)
##    Tamano_Muestra     Media   Mediana Desviacion_Estandar     Varianza
## 1               5 0.9124000 1.0000000          0.12180276 0.0148359118
## 2              10 0.9040000 0.9000000          0.08831534 0.0077995992
## 3              15 0.9089333 0.9333333          0.07502821 0.0056292318
## 4              20 0.9101000 0.9000000          0.06672251 0.0044518938
## 5              30 0.9101333 0.9000000          0.05262624 0.0027695213
## 6              50 0.9066000 0.9200000          0.03855662 0.0014866132
## 7              60 0.9098000 0.9166667          0.03458833 0.0011963527
## 8             100 0.9091800 0.9100000          0.03007552 0.0009045367
## 9             200 0.9098900 0.9100000          0.02000596 0.0004002384
## 10            500 0.9112360 0.9120000          0.01366652 0.0001867739
##     Asimetria
## 1  -1.1669256
## 2  -0.8789934
## 3  -0.7219211
## 4  -0.6494139
## 5  -0.4423549
## 6  -0.3341238
## 7  -0.3447533
## 8  -0.2790996
## 9  -0.3673508
## 10 -0.0686704

Análisis punto D:

En primer lugar, observamos que el valor de estimadores ^P estan muy cercanos del valor del parametro P = 0.1 y P=0.9, respectivamente.

Al analizar el comportamiento de los datos al repetir el experimento 500 veces, se evidencia que para lotes con el 10% de enfermas la media , la mediana y la moda distan un poco entre si, en especial para la última medida lo que nos podria indicar que los datos probablemente puede que no se encuentren distribuidos normalmente. Pues es importante mencionar que, para que se encuentren distribudios simetricamente, deberian ser iguales los 3 valores. Para el caso, con el 90% de enfermas, los valores de la media, la mediana y la moda, son casi iguales, lo que sugiere que la distribución podria ser simetrica.

De igual manera, para el primer caso las pruebas de bondad, nos arrojaron un valor P < 0.05; rechazando de esta manera la hipotesis nula, por lo cual se podria concluir que los datos no siguen una distribución normal. En el segundo caso, el valor P > 0.05, concluyendo que los datos podrian seguir una distribución normal

Del analisis del grafico Q-Q, podemos observar que los puntos se desvian de manera significativa de la linea de referencia; en comparación a cuando la probabilidad de enfermas es del 50.

Con respecto al análisis para la aleatoredad en las muestras, se evidencia que a pesar de que el valor W y P son mas alto (el W acercandose a 1 y el P a 0.05 o más), al incrementar el número de muestras, en los 2 casos se rechazaron las hipotesis nula, dado que utilizando el mayor valor n=500, los valores P <0.05, sugiriendo que los datos no se distribuyen normalmente.

Explicación del Código

  1. Generación de la Población:
    • rbinom(population_size, 1, percent_sick) genera una población binomial con una probabilidad de enfermedad del 50%.
  2. Cálculo de Proporción Muestral:
    • calculate_proportion calcula la proporción de individuos enfermos en una muestra aleatoria de la población.
  3. Simulación:
    • simulate_proportions realiza 500 simulaciones para cada tamaño de muestra y almacena los resultados.
  4. Análisis:
    • analyze_simulation calcula estadísticas descriptivas como la media, varianza, asimetría y curtosis de las proporciones muestrales.
  5. Visualización y Pruebas de Normalidad:
    • plot_and_test crea histogramas, gráficos QQ y realiza la prueba de Shapiro-Wilk para evaluar la normalidad de la distribución de las proporciones muestrales.
  6. Variación en el Porcentaje de Enfermos:
    • Se repite la simulación para poblaciones con diferentes porcentajes de individuos enfermos y se comparan los resultados.

Este código ayuda a verificar el Teorema del Límite Central en diferentes condiciones y tamaños de muestra, y permite analizar la aproximación a la normalidad en función de estas variables.

PROBLEMA 4

Introducción

library(ggplot2)
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ lubridate 1.9.3     ✔ tibble    3.2.1
## ✔ purrr     1.0.2     ✔ tidyr     1.3.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
combustible <- c(7.69, 4.97, 4.56, 6.49, 4.34, 6.24,4.45)
print(combustible)
## [1] 7.69 4.97 4.56 6.49 4.34 6.24 4.45

La variable xbarra contiene la media de la muestra de datos (5.53) de consumo de combustible, que es la mejor estimación de la media poblacional basada en la muestra proporcionada.

xbarra<-mean(combustible)
xbarra
## [1] 5.534286

En primer lugar, se tomó la muestra original “combustible” y se generaron 1000 muestras remuestreadas con reemplazo de la muestra original. Despues, calculamos las dimensiones de la matriz, de donde se obtuvieron 7 filas y 10000 columnas; lo que nos indica que las muestras resultantes se distribuyeron por columna.

boostrap<-replicate(n=1000,sample(combustible,replace=TRUE))
dim(boostrap)
## [1]    7 1000

Posteriormente calculamos las medias de cada muestra remuestreada y ordenamos las medias de mayor a menor.

medias<-apply(boostrap, MARGIN=2, FUN=mean)
# Ordenamiento de las medias bootstrap de menor a mayor 
medias_ordenadas <- sort(medias)

Método 1

Este método utiliza la técnica de remuestreo bootstrap para generar muchas muestras de la distribución de la media. Luego, se calculan los percentiles 2.5% y 97.5% de las medias de estas muestras remuestreadas para obtener el intervalo de confianza.

Resultados del Método 1:

Percentil 2.5%: 4.705429 Percentil 97.5%: 6.444536

###Método 1###
# Encontrando el percentil 2.5% más bajo
percentil_2.5<- quantile(medias_ordenadas, 0.025)
# Encontrando el percentil 97.5% más alto
percentil_97.5<- quantile(medias_ordenadas, 0.975)
percentil_2.5
##     2.5% 
## 4.721429
percentil_97.5
##    97.5% 
## 6.418964
#(P2.5)=4.705429 y (P97.5)=6.444536

Método 2 Este método también utiliza las medias de las muestras bootstrap, pero en lugar de tomar directamente los percentiles, calcula el intervalo de confianza utilizando la media de todas las medias bootstrap, la cual es multiplicada por 2 y a la que se le restan los percentiles encontrados en el Método 1.

###Método 2###
limite_inferior<- 2*mean(medias_ordenadas)-percentil_97.5
limite_superior<- 2*mean(medias_ordenadas)-percentil_2.5
limite_inferior
##   97.5% 
## 4.60371
limite_superior
##     2.5% 
## 6.301246
#Histograma resultante
hist(medias_ordenadas, main="Histograma de valores medios", 
     xlab = "Medias de las muestras", 
     ylab = "Frecuencia",
     col = "blue")
color_met_1 = "red"
color_met_2 = "green"
abline(v = percentil_2.5, col = color_met_1, lwd = 2)
abline(v = percentil_97.5, col = color_met_1, lwd = 2)
abline(v = limite_inferior, col = color_met_2, lwd = 2)
abline(v = limite_superior, col = color_met_2, lwd = 2)

legend("topright", legend = c("Método 1", "Método 2"), 
       col = c(color_met_1, color_met_2), lwd = c(1, 2))

De acuerdo a los resultados obtenidos, el promedio real(con un cierto nivel de confianza) se encuentra entre el rango de valores; para el Método 1 [4.768464-6.534536] y para el Método 2 [4.58147-6.347541]; existiendo una variación poco representativa (aprox. de 0,19) para los limites obtenidos. De igual manera, se identificó

Ambos métodos se basan en el remuestreo bootstrap y proporcionan intervalos de confianza para la media. El Método 1 es más directo y utiliza los percentiles de las medias remuestreadas. El Método 2 ajusta estos percentiles basándose en la media de las medias bootstrap.

Los intervalos de confianza resultantes pueden ser ligeramente diferentes debido a la forma en que se calculan. El Método 1 da un intervalo basado directamente en la variabilidad de las medias remuestreadas, mientras que el Método 2 intenta centrar el intervalo alrededor de la media de las medias bootstrap.

En cuanto a la confianza en estas estimaciones, ambas son técnicas válidas y comúnmente utilizadas en estadística. La elección entre ellas puede depender de la preferencia del analista o de consideraciones específicas del contexto del estudio. En general, el bootstrap es una herramienta poderosa cuando no se pueden hacer suposiciones sobre la distribución subyacente de la población, y ambos métodos deberían proporcionar estimaciones razonables del intervalo de confianza para la media.

El código también incluye la creación de un histograma con las medias de las muestras bootstrap y las líneas que representan los intervalos de confianza de ambos métodos, lo que puede ayudar a visualizar la distribución de las medias y la comparación de los intervalos.

Para concluir, se puede confiar en las estimaciones de ambos métodos, pero es importante recordar que el intervalo de confianza es una estimación probabilística y no una garantía de que la media real de la población caiga dentro de ese rango en todos los casos.

PROBLEMA 5

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

#Gráfica 1

#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("pequeño","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("pequeño","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="Tamaño 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="Tamaño 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 tamaño 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")    
## 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.
# Interactive plot
plotly::ggplotly(p)
if (!require("pwr")) install.packages("pwr")
if (!require("ggplot2")) install.packages("ggplot2")
library(pwr)
library(ggplot2)

# Definimos los tamaños de efecto y los niveles de significancia a analizar
effect_sizes <- c(0.1, 0.2, 0.5, 0.8)
significance_levels <- c(0.01, 0.05, 0.10, 0.15, 0.20)
sample_sizes <- c(20, 60, 100, 140)

# Creamos un data frame vacío almacenar los resultados
results <- data.frame()

# Bucle para calcular la potencia para cada combinación de tamaño de efecto, nivel de significancia y tamaño de muestra
for (d in effect_sizes) {
  for (alpha in significance_levels) {
    for (n in sample_sizes) {
      # Calculamos la potencia de la prueba t
      power <- pwr.t.test(d = d, n = n, sig.level = alpha, type = "two.sample", alternative = "two.sided")$power
      
      # Almacenamos los resultados en el data frame
      results <- rbind(results, data.frame(effect_size = d, significance_level = alpha, sample_size = n, power = power))
    }
  }
}

print(results)
##    effect_size significance_level sample_size      power
## 1          0.1               0.01          20 0.01350839
## 2          0.1               0.01          60 0.02180466
## 3          0.1               0.01         100 0.03094690
## 4          0.1               0.01         140 0.04089447
## 5          0.1               0.05          20 0.06095912
## 6          0.1               0.05          60 0.08444035
## 7          0.1               0.05         100 0.10837184
## 8          0.1               0.05         140 0.13264737
## 9          0.1               0.10          20 0.11632725
## 10         0.1               0.10          60 0.14989019
## 11         0.1               0.10         100 0.18296555
## 12         0.1               0.10         140 0.21548307
## 13         0.1               0.15          20 0.16967584
## 14         0.1               0.15          60 0.20918436
## 15         0.1               0.15         100 0.24737543
## 16         0.1               0.15         140 0.28424260
## 17         0.1               0.20          20 0.22176378
## 18         0.1               0.20          60 0.26476788
## 19         0.1               0.20         100 0.30578795
## 20         0.1               0.20         140 0.34488342
## 21         0.2               0.01          20 0.02513156
## 22         0.2               0.01          60 0.06749459
## 23         0.2               0.01         100 0.12034615
## 24         0.2               0.01         140 0.18076972
## 25         0.2               0.05          20 0.09456733
## 26         0.2               0.05          60 0.19237465
## 27         0.2               0.05         100 0.29064587
## 28         0.2               0.05         140 0.38512417
## 29         0.2               0.10          20 0.16472928
## 30         0.2               0.10          60 0.29234051
## 31         0.2               0.10         100 0.40804897
## 32         0.2               0.10         140 0.51019013
## 33         0.2               0.15          20 0.22687425
## 34         0.2               0.15          60 0.36928271
## 35         0.2               0.15         100 0.49060948
## 36         0.2               0.15         140 0.59214983
## 37         0.2               0.20          20 0.28417898
## 38         0.2               0.20          60 0.43350233
## 39         0.2               0.20         100 0.55515082
## 40         0.2               0.20         140 0.65304748
## 41         0.5               0.01          20 0.14395508
## 42         0.5               0.01          60 0.54945334
## 43         0.5               0.01         100 0.82382249
## 44         0.5               0.01         140 0.94322683
## 45         0.5               0.05          20 0.33793903
## 46         0.5               0.05          60 0.77526589
## 47         0.5               0.05         100 0.94042720
## 48         0.5               0.05         140 0.98640721
## 49         0.5               0.10          20 0.46406530
## 50         0.5               0.10          60 0.85949024
## 51         0.5               0.10         100 0.96984801
## 52         0.5               0.10         140 0.99426817
## 53         0.5               0.15          20 0.54908911
## 54         0.5               0.15          60 0.90097004
## 55         0.5               0.15         100 0.98154365
## 56         0.5               0.15         140 0.99688999
## 57         0.5               0.20          20 0.61339503
## 58         0.5               0.20          60 0.92615073
## 59         0.5               0.20         100 0.98766879
## 60         0.5               0.20         140 0.99810756
## 61         0.8               0.01          20 0.43797260
## 62         0.8               0.01          60 0.95941915
## 63         0.8               0.01         100 0.99879065
## 64         0.8               0.01         140 0.99997722
## 65         0.8               0.05          20 0.69340420
## 66         0.8               0.05          60 0.99148114
## 67         0.8               0.05         100 0.99987838
## 68         0.8               0.05         140 0.99999876
## 69         0.8               0.10          20 0.79942625
## 70         0.8               0.10          60 0.99665226
## 71         0.8               0.10         100 0.99996732
## 72         0.8               0.10         140 0.99999976
## 73         0.8               0.15          20 0.85442599
## 74         0.8               0.15          60 0.99826590
## 75         0.8               0.15         100 0.99998680
## 76         0.8               0.15         140 0.99999992
## 77         0.8               0.20          20 0.88896170
## 78         0.8               0.20          60 0.99898181
## 79         0.8               0.20         100 0.99999360
## 80         0.8               0.20         140 0.99999997

Análisis:

-> Variando los niveles de significancia:

Al mantener estaticos los valores del tamaño de efecto y muestra y al iterar sobre el nivel de significancia, se puede observar una relación directa: a medida que el nivel de significancia aumenta, el tamaño de la potencia aumenta. Tomando como ejemplo, las iteraciones realizadas con tamaño de efecto 0.1 de muestra 20, se logró evidenciar un impacto en la potencia de la siguiente manera:

  • Para nivel de significancia (0.01), el valor de la potencia fue: 0.0135
  • Para nivel de significancia (0.05), el valor de la potencia fue: 0.0609
  • Para nivel de significancia (0.10), el valor de la potencia fue: 0.1163
  • Para nivel de significancia (0.15), el valor de la potencia fue: 0.1696
  • Para nivel de significancia (0.20), el valor de la potencia fue: 0.2217

En conclusión, la potencia de la prueba aumenta a medida que el nivel de significancia se incrementa. Esto se debe a que un nivel de significancia más alto implica un umbral menos estricto para rechazar la hipótesis nula, lo que a su vez aumenta la probabilidad de detectar un efecto si realmente existe (es decir, de rechazar la hipótesis nula cuando la hipótesis alternativa es verdadera).

-> Variando los tamaño de muestra:

Al mantener estaticos los valores del tamaño de efecto y de nivel de significancia y al iterar sobre el tamaño de la muestra, se puede observar una relación directa: a medida que el tamaño de la muestra, el tamaño de la potencia aumenta. Tomando como ejemplo, las iteraciones realizadas con tamaño de efecto 0.1 nivel de significancia 0.01, se logró evidenciar un impacto en la potencia de la siguiente manera:

  • Para tamaño de muestra (20), el valor de la potencia fue: 0.0135
  • Para tamaño de muestra (60), el valor de la potencia fue: 0.0218
  • Para tamaño de muestra (100), el valor de la potencia fue: 0.0309
  • Para tamaño de muestra (140), el valor de la potencia fue: 0.0408

En conclusión, la potencia también aumenta con el tamaño de la muestra. Esto es esperado, ya que muestras más grandes proporcionan más información y reducen la variabilidad del estimador, lo que mejora la capacidad de la prueba para detectar diferencias significativas entre grupos.

-> Variando el tamaño del efecto:

Por otro lado, se logra evidenciar que a medida que aumenta el tamaño el efecto, la potencia aumenta independientemente del tamaño de nivel de significancia y de muestras, tal como se evidenció en la tabla generada, donde manteniendo un nivel de significancia de 0.01 y de muestra de 20, se logran observar lo siguiente:

  • Para tamaño de efecto (0.1), el valor de la potencia fue: 0.0135
  • Para tamaño de efecto (0.2), el valor de la potencia fue: 0.0251
  • Para tamaño de efecto (0.5), el valor de la potencia fue: 0.1439
  • Para tamaño de efecto (0.8), el valor de la potencia fue: 0.4379

Como es lógico, la potencia es mayor para tamaños de efecto más grandes. Un tamaño de efecto más grande indica una diferencia más pronunciada entre los grupos que se están comparando, lo que facilita su detección por parte de la prueba estadística.

Caso 2: Variando los tamaños muestrales

Para poder realizar el análisis con 5 niveles diferentes de significancia (α), se debe ajustar el código siguinete, para iterar sobre esos valores. Se agrega un bucle para calcular la potencia con distintos niveles de α y generar gráficos comparativos.

Código inicial:

library(dplyr)    
library(tidyr)    #Para manipulación de datos: separate, gather, spread
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)

Ajustes:

  1. Incluir 5 valores de α: por ejemplo, 0.01, 0.025, 0.05, 0.075, y 0.1.
  2. Ajustar el código para iterar sobre esos niveles y generar tablas de resultados para cada uno.
library(dplyr)     
library(tidyr)    
library(ggplot2)   
library(plotly)   
library(pwr)      

# Valores de significancia
sig_levels <- c(0.01, 0.025, 0.05, 0.075, 0.1)

# Inicializar tabla para almacenar resultados
ptab_total <- NULL       

for (alpha in sig_levels) {
  ptab <- cbind(NULL, NULL) 
  
  for (i in seq(0, 1, length.out = 200)) {
    pwrt1 <- pwr.t2n.test(n1 = 28, n2 = 1406, 
                          sig.level = alpha, power = NULL, 
                          d = i, alternative = "two.sided")
    pwrt2 <- pwr.t2n.test(n1 = 144, n2 = 1290, 
                          sig.level = alpha, power = NULL, 
                          d = i, alternative = "two.sided")
    pwrt3 <- pwr.t2n.test(n1 = 287, n2 = 1147, 
                          sig.level = alpha, power = NULL, 
                          d = i, alternative = "two.sided")
    pwrt4 <- pwr.t2n.test(n1 = 430, n2 = 1004, 
                          sig.level = alpha, power = NULL, 
                          d = i, alternative = "two.sided")
    pwrt5 <- pwr.t2n.test(n1 = 574, n2 = 860, 
                          sig.level = alpha, power = NULL, 
                          d = i, alternative = "two.sided")
    pwrt6 <- pwr.t2n.test(n1 = 717, n2 = 717, 
                          sig.level = alpha, power = NULL, 
                          d = i, alternative = "two.sided")
    
    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))
  }
  
  # Añadir columna de nivel de significancia
  ptab <- cbind(sig.level = alpha, seq_len(nrow(ptab)), ptab)
  
  # Unir todos los resultados
  ptab_total <- rbind(ptab_total, ptab)
}

# Agregar nombres a las columnas
colnames(ptab_total) <- c("sig.level", "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")

# Transformar el data frame y preparar para graficar
temp1 <- ptab_total %>% as.data.frame() %>% gather(key = name, value = val, 3:14)
temp2 <- temp1 %>% separate(col = name, into = c("samples", "pruebas"), sep = ";")
temp3 <- temp2 %>% spread(key = pruebas, value = val)
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 comparativa con diferentes valores de significancia
p <- ggplot(temp3, aes(x = `effect size`, y = power, color = samples)) +
  geom_line(size = 1) +
  facet_wrap(~sig.level, nrow = 3) +  # Comparar en distintos niveles de significancia
  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")

# Hacer el gráfico interactivo
plotly::ggplotly(p)