Al finalizar esta clase, los estudiantes serán capaces de:
Somos un laboratorio de control de calidad de una empresa de jugos y queremos evaluar cómo diferentes métodos de conservación afectan el contenido de vitamina C (ácido ascórbico) en jugo de naranja comercial durante 7 días de almacenamiento.
Población: Todos los jugos de naranja pasteurizados comercializados en la empresa
Muestra: Envases de jugo de naranja de una marca específica, mismo lote de producción
Unidad experimental: Cada envase individual de 200 mL de jugo de naranja
Unidad de observación: La medición del contenido de vitamina C (mg/100 mL) en cada envase
Tratamientos (factor: método de conservación, 3 niveles):
Variable dependiente: Contenido de vitamina C en mg/100 mL (cuantitativa continua)
Variables de control:
Antes de comenzar cualquier experimento, debemos determinar cuántas repeticiones necesitamos por tratamiento para tener suficiente poder estadístico para detectar diferencias si existen.
Opción A - Basado en literatura: Buscamos estudios similares y calculamos f a partir de las medias y desviaciones reportadas.
Opción B - Estudio piloto: Realizamos un pequeño ensayo con 3-5 repeticiones por tratamiento y estimamos f.
Opción C - Convención de Cohen: Si no tenemos información previa, usamos f = 0.25 (efecto mediano) como punto de partida conservador.
Opción D - Cálculo a partir de diferencias esperadas: Si esperamos que T1 tenga media 50 mg/100mL, T2 tenga 45 mg/100mL, T3 tenga 40 mg/100mL, con desviación estándar esperada de 5 mg/100mL, podemos calcular f.
\[ f = \sqrt{\frac{\sum_{i=1}^{k}(\mu_i - \mu)^2 / k}{\sigma^2}} \]
Donde:
Abrimos RStudio y creamos un nuevo script llamado 01_calculo_poder.R:
# ============================================================================
# CÁLCULO DE TAMAÑO MUESTRAL CON ANÁLISIS DE POTENCIA
# Diseño: Completamente Aleatorizado (DCA)
# Experimento: Efecto de conservación en vitamina C de jugo de naranja
# ============================================================================
# Instalar y cargar paquete pwr (solo la primera vez)
# install.packages("pwr")
library(pwr)
# ----------------------------------------------------------------------------
# OPCIÓN 1: Cálculo basado en convención de Cohen (efecto mediano)
# ----------------------------------------------------------------------------
# Parámetros
k <- 3 # Número de tratamientos (T1, T2, T3)
f <- 0.25 # Tamaño del efecto mediano (convención de Cohen para ANOVA)
sig.level <- 0.05 # Nivel de significación (5%)
power <- 0.80 # Potencia deseada (80%)
# Calcular tamaño muestral por grupo
resultado <- pwr.anova.test(k = k, f = f, sig.level = sig.level, power = power)
# Mostrar resultados
print("=== ANÁLISIS DE POTENCIA - EFECTO MEDIANO ===")
print(resultado)
print(paste("Repeticiones por tratamiento (redondeado hacia arriba):",
ceiling(resultado$n)))
print(paste("Total de unidades experimentales:",
ceiling(resultado$n) * k))
# ----------------------------------------------------------------------------
# OPCIÓN 2: Cálculo basado en diferencias esperadas (más preciso)
# ----------------------------------------------------------------------------
# Supongamos que, basado en literatura o estudio piloto, esperamos:
# T1 (25°C): media = 45 mg/100mL
# T2 (4°C): media = 50 mg/100mL
# T3 (-5°C): media = 52 mg/100mL
# Desviación estándar dentro de grupos: σ = 4 mg/100mL
# Calcular medias
mu1 <- 45
mu2 <- 50
mu3 <- 52
mu_general <- mean(c(mu1, mu2, mu3)) # Media general
sigma <- 4 # Desviación estándar esperada
# Calcular tamaño del efecto f
varianza_entre <- sum((c(mu1, mu2, mu3) - mu_general)^2) / k
f_calculado <- sqrt(varianza_entre / sigma^2)
print("=== TAMAÑO DEL EFECTO CALCULADO ===")
print(paste("Media general:", round(mu_general, 2), "mg/100mL"))
print(paste("Desviación estándar esperada:", sigma, "mg/100mL"))
print(paste("Tamaño del efecto f calculado:", round(f_calculado, 3)))
# Recalcular tamaño muestral con f calculado
resultado_realista <- pwr.anova.test(k = k, f = f_calculado,
sig.level = sig.level, power = power)
print("=== ANÁLISIS DE POTENCIA - BASADO EN DIFERENCIAS ESPERADAS ===")
print(resultado_realista)
print(paste("Repeticiones por tratamiento (redondeado hacia arriba):",
ceiling(resultado_realista$n)))
print(paste("Total de unidades experimentales:",
ceiling(resultado_realista$n) * k))
# ----------------------------------------------------------------------------
# OPCIÓN 3: Análisis de sensibilidad - ¿Qué pasa si cambiamos la potencia?
# ----------------------------------------------------------------------------
# Crear secuencia de potencias desde 0.70 hasta 0.95
potencias <- seq(0.70, 0.95, by = 0.05)
n_por_potencia <- numeric(length(potencias))
for(i in 1:length(potencias)){
n_por_potencia[i] <- pwr.anova.test(k = k, f = f_calculado,
sig.level = sig.level,
power = potencias[i])$n
}
# Crear tabla resumen
tabla_sensibilidad <- data.frame(
Potencia = potencias,
n_por_tratamiento = round(n_por_potencia, 1),
Total_unidades = round(n_por_potencia * k, 1)
)
print("=== ANÁLISIS DE SENSIBILIDAD ===")
print(tabla_sensibilidad)
# Gráfico de sensibilidad
library(ggplot2)
ggplot(tabla_sensibilidad, aes(x = Potencia, y = Total_unidades)) +
geom_line(color = "steelblue", linewidth = 1.5) +
geom_point(color = "steelblue", size = 3) +
geom_hline(yintercept = ceiling(resultado_realista$n) * k,
linetype = "dashed", color = "red") +
annotate("text", x = 0.85, y = ceiling(resultado_realista$n) * k + 2,
label = paste("Punto elegido:\nTotal =",
ceiling(resultado_realista$n) * k, "unidades"),
color = "red", hjust = 0) +
labs(title = "Relación entre Potencia Estadística y Tamaño Muestral Total",
subtitle = "DCA con 3 tratamientos, f = 0.391",
x = "Potencia Estadística (1 - β)",
y = "Total de Unidades Experimentales",
caption = "Línea roja punteada: punto seleccionado (potencia = 0.80)") +
theme_minimal() +
theme(plot.title = element_text(face = "bold", size = 14),
plot.subtitle = element_text(size = 11))
# Guardar gráfico
ggsave("grafico_potencia_vs_n.png", width = 8, height = 5, dpi = 300)
# ----------------------------------------------------------------------------
# DECISIÓN FINAL
# ----------------------------------------------------------------------------
n_final <- ceiling(resultado_realista$n) # Redondear hacia arriba siempre
total_unidades <- n_final * k
print("=== DECISIÓN DE DISEÑO ===")
print(paste("Número de tratamientos (k):", k))
print(paste("Repeticiones por tratamiento (n):", n_final))
print(paste("Total de unidades experimentales (N):", total_unidades))
print(paste("Potencia estadística alcanzada:",
round(resultado_realista$power, 3)))
print(paste("Nivel de significación:", sig.level))
print(paste("Tamaño del efecto detectable (f):", round(f_calculado, 3)))
# Guardar parámetros para uso posterior
parametros <- list(
k = k,
n_por_tratamiento = n_final,
total_unidades = total_unidades,
f = f_calculado,
sig.level = sig.level,
power = resultado_realista$power
)
save(parametros, file = "parametros_diseno.RData")
print("Parámetros guardados en 'parametros_diseno.RData'")Al ejecutar el script anterior, obtendremos resultados como:
=== ANÁLISIS DE POTENCIA - BASADO EN DIFERENCIAS ESPERADAS ===
Balanced one-way analysis of variance power calculation
k = 3
n = 7.672891 ← ¡IMPORTANTE! Necesitamos 7.67 por grupo
f = 0.390625 ← Tamaño del efecto calculado
sig.level = 0.05
power = 0.8 ← Potencia deseada
NOTE: n is number in each group
Repeticiones por tratamiento (redondeado hacia arriba): 8
Total de unidades experimentales: 24Interpretación:
¿Por qué es importante este cálculo?
Basado en el cálculo de potencia, necesitamos:
Importante: La aleatorización es importante para evitar sesgos y cumplir el supuesto de independencia.
Creamos un segundo script en RStudio:
02_aleatorizacion.R
# ============================================================================
# ALEATORIZACIÓN DEL DISEÑO COMPLETAMENTE ALEATORIZADO
# ============================================================================
# Cargar parámetros guardados
load("parametros_diseno.RData")
# Verificar parámetros
print(parametros)
# ----------------------------------------------------------------------------
# Crear lista de unidades experimentales
# ----------------------------------------------------------------------------
# Identificadores únicos para cada unidad experimental
unidades <- paste0("UE", sprintf("%02d", 1:parametros$total_unidades))
# Crear vector de tratamientos (8 de cada uno)
tratamientos <- rep(c("T1_25C", "T2_4C", "T3_menos5C"),
each = parametros$n_por_tratamiento)
# Verificar balance
table(tratamientos)
# ----------------------------------------------------------------------------
# Aleatorizar asignación de tratamientos a unidades
# ----------------------------------------------------------------------------
# Fijar semilla para reproducibilidad
set.seed(2026) # Cualquier número entero
# Aleatorizar el orden de los tratamientos
tratamientos_aleatorizados <- sample(tratamientos)
# Crear tabla de asignación
asignacion <- data.frame(
Unidad_Experimental = unidades,
Tratamiento_Asignado = tratamientos_aleatorizados,
Orden_Ejecucion = 1:parametros$total_unidades
)
# Ver distribución después de aleatorizar
print("=== DISTRIPUCIÓN DE TRATAMIENTOS DESPUÉS DE ALEATORIZAR ===")
print(table(asignacion$Tratamiento_Asignado))
# ----------------------------------------------------------------------------
# Guardar planilla de aleatorización para imprimir y usar en campo
# ----------------------------------------------------------------------------
# Ordenar por unidad experimental para facilitar lectura
asignacion <- asignacion[order(asignacion$Unidad_Experimental), ]
# Agregar columnas para registro manual
asignacion$Fecha_Inicio <- ""
asignacion$Fecha_Fin <- ""
asignacion$Temperatura_Real <- ""
asignacion$Observaciones <- ""
asignacion$Responsable <- ""
# Exportar a Excel para imprimir y llevar al laboratorio
# install.packages("writexl")
library(writexl)
write_xlsx(asignacion, "planilla_aleatorizacion.xlsx")
print("Planilla de aleatorización guardada como 'planilla_aleatorizacion.xlsx'")
print("IMPRIMIR esta planilla y llevar al laboratorio")
# ----------------------------------------------------------------------------
# Visualizar aleatorización
# ----------------------------------------------------------------------------
library(ggplot2)
ggplot(asignacion, aes(x = Orden_Ejecucion, y = Unidad_Experimental,
fill = Tratamiento_Asignado)) +
geom_tile(color = "white", size = 0.5) +
scale_fill_brewer(palette = "Set2",
labels = c("T1: 25°C", "T2: 4°C", "T3: -18°C")) +
labs(title = "Esquema de aleatorización - Diseño Completamente Aleatorizado",
subtitle = paste("24 unidades experimentales, 8 repeticiones por tratamiento"),
x = "Orden de Ejecución",
y = "Unidad Experimental",
fill = "Tratamiento") +
theme_minimal() +
theme(axis.text.y = element_text(size = 8),
plot.title = element_text(face = "bold", size = 14),
legend.title = element_text(face = "bold"))
ggsave("esquema_aleatorizacion.png", width = 10, height = 6, dpi = 300)
# ----------------------------------------------------------------------------
# Generar lista de trabajo por tratamiento (para organizar el laboratorio)
# ----------------------------------------------------------------------------
# Lista por tratamiento
for(trat in unique(asignacion$Tratamiento_Asignado)){
cat("\n=== ", trat, " ===\n")
unidades_trat <- asignacion$Unidad_Experimental[asignacion$Tratamiento_Asignado == trat]
cat("Unidades:", paste(unidades_trat, collapse = ", "), "\n")
cat("Total:", length(unidades_trat), "unidades\n")
}Importante: Los datos deben registrarse inmediatamente durante la ejecución del experimento, no de memoria al final.
================================================================================
PLANILLA DE CAMPO - REGISTRO DE DATOS
Vitamina C en Jugo de Naranja
================================================================================
Fecha de inicio: _______________ Responsable: _________________________
Condiciones ambientales:
- Temperatura laboratorio: ______ °C
- Humedad relativa: ______ %
--------------------------------------------------------------------------------
UNIDAD | TRATAMIENTO | FECHA INI | FECHA FIN | VOL. MUESTRA | OBSERVACIONES
--------------------------------------------------------------------------------
UE01 | | | | |
UE02 | | | | |
UE03 | | | | |
... (continuar hasta UE24)
--------------------------------------------------------------------------------
================================================================================
REGISTRO DE MEDICIONES DE VITAMINA C
Método: Titulación con DCPIP
================================================================================
Fecha de análisis: _______________ Analista: _________________________
Solución estándar de ácido ascórbico: ______ mg/mL
Volumen de muestra titulado: ______ mL
Factor de la bureta: ______
--------------------------------------------------------------------------------
UE | Trat. | Vol. DCPIP (mL) | Cálculo | Vit. C (mg/100mL) | Iniciales
--------------------------------------------------------------------------------
UE01 | | | | |
UE02 | | | | |
UE03 | | | | |
... (24 filas)
--------------------------------------------------------------------------------
Cálculo:
Vitamina C (mg/100mL) = (Vol. DCPIP × Factor × 100) / Vol. muestra
================================================================================Día 0 - Recepción y aleatorización:
Días 1-6 - Monitoreo:
Día 7 - Análisis de vitamina C:
\[ \text{Vit. C (mg/100mL)} = \frac{V_{DCPIP} \times F \times 100}{V_{muestra}} \]
Donde:
- \(V_{DCPIP}\) = volumen de DCPIP
gastado (mL)
- \(F\) = factor de la solución de
DCPIP (mg ácido ascórbico/mL DCPIP)
- \(V_{muestra}\) = volumen de muestra
(mL)
4. Registrar todos los valores en la planilla de campo con iniciales del
analista
Después de completar el ensayo, creamos una planilla Excel estructurada para importar a R.
Nombre del archivo:
datos_vitamina_c.xlsx
Estructura de la hoja “Datos_Raw”:
| UE | Tratamiento | Temp_Conservacion | Vol_DCPIP_mL | VitaminaC_mg100mL | Fecha_Analisis | Analista | Observaciones |
|---|---|---|---|---|---|---|---|
| UE01 | T1_25C | 25 | 12.5 | 45.2 | 15/08/2026 | AB | Sin novedades |
| UE02 | T2_4C | 4 | 14.1 | 51.0 | 15/08/2026 | AB | Sin novedades |
| UE03 | T3_menos18C | -18 | 14.8 | 53.5 | 15/08/2026 | AB | Sin novedades |
| … | … | … | … | … | … | … | … |
Reglas de oro para la carga de datos:
VitaminaC_mg100mL)Hoja adicional “Diccionario”:
| Variable | Tipo | Descripción | Unidad | Rango válido |
|---|---|---|---|---|
| UE | Caracter | Identificador único de unidad experimental | - | UE01-UE24 |
| Tratamiento | Factor | Método de conservación | - | T1_25C, T2_4C, T3_menos5C |
| Temp_Conservacion | Numérica | Temperatura de conservación | °C | -20 a 30 |
| Vol_DCPIP_mL | Numérica | Volumen de DCPIP gastado en titulación | mL | 0-50 |
| VitaminaC_mg100mL | Numérica | Contenido de vitamina C calculado | mg/100mL | 0-100 |
| Fecha_Analisis | Fecha | Fecha de realización del análisis | dd/mm/yyyy | - |
| Analista | Caracter | Iniciales del analista | - | - |
| Observaciones | Caracter | Notas relevantes | - | - |
Creamos el script 03_importar_datos.R:
# ============================================================================
# IMPORTACIÓN DE DATOS
# ============================================================================
# Instalar paquetes necesarios (solo la primera vez)
# install.packages("readxl")
# install.packages("tidyverse")
# install.packages("janitor")
# Cargar librerías
library(readxl)
library(tidyverse)
library(janitor)
# ----------------------------------------------------------------------------
# Importar datos desde Excel
# ----------------------------------------------------------------------------
# Establecer directorio de trabajo (cambiar según tu ruta)
# setwd("C:/Users/TuUsuario/Documentos/Diseño_Experimental")
# Leer archivo Excel
datos_raw <- read_excel("datos_vitamina_c.xlsx", sheet = "Datos_Raw")
# Ver estructura inicial
str(datos_raw)
dim(datos_raw)
head(datos_raw)
# ----------------------------------------------------------------------------
# Limpieza de datos
# ----------------------------------------------------------------------------
# Renombrar columnas para mayor claridad
datos <- datos_raw %>%
rename(
ue = UE,
tratamiento = Tratamiento,
temp_Conservacion = Temp_Conservacion,
vol_dcpip_ml = Vol_DCPIP_mL,
vitamina_c = VitaminaC_mg100mL,
fecha_analisis = Fecha_Analisis,
analista = Analista,
observaciones = Observaciones
)
# ----------------------------------------------------------------------------
# Verificar tipos de variables
# ----------------------------------------------------------------------------
# Convertir tratamiento a factor (variable categórica)
datos$tratamiento <- as.factor(datos$tratamiento)
# Convertir ue a factor (identificador, no variable numérica)
datos$ue <- as.factor(datos$ue)
# Convertir analista a factor
datos$analista <- as.factor(datos$analista)
# Verificar niveles del factor tratamiento
levels(datos$tratamiento)
# ----------------------------------------------------------------------------
# Verificación de calidad de datos
# ----------------------------------------------------------------------------
# Verificar valores faltantes
print("=== VALORES FALTANTES POR COLUMNA ===")
colSums(is.na(datos))
# Verificar valores duplicados
print("=== ¿HAY UNIDADES EXPERIMENTALES DUPLICADAS? ===")
duplicados <- datos[duplicated(datos$ue), ]
if(nrow(duplicados) > 0){
print("ALERTA: Hay unidades duplicadas:")
print(duplicados)
} else {
print("OK: No hay duplicados")
}
# Verificar rangos de valores (control de calidad)
summary(datos)
# Verificar que todas las unidades estén presentes
print(paste("Unidades esperadas:", 24))
print(paste("Unidades en datos:", nrow(datos)))
print(paste("Unidades únicas:", length(unique(datos$ue))))
# ----------------------------------------------------------------------------
# Guardar datos limpios para análisis posterior
# ----------------------------------------------------------------------------
# Guardar como archivo .RData (formato nativo de R)
save(datos, file = "datos_limpios.RData")
# Guardar como CSV (formato universal)
write.csv(datos, "datos_limpios.csv", row.names = FALSE)
# ----------------------------------------------------------------------------
# Resumen final
# ----------------------------------------------------------------------------
print("=== RESUMEN DE IMPORTACIÓN ===")
print(paste("Total de observaciones:", nrow(datos)))
print(paste("Total de variables:", ncol(datos)))
print("Variables:")
print(names(datos))
print("Estructura final:")
str(datos)Creamos el script principal 04_analisis_completo.R:
# ============================================================================
# ANÁLISIS ESTADÍSTICO - DISEÑO COMPLETAMENTE ALEATORIZADO
# Efecto de métodos de conservación en vitamina C de jugo de naranja
# ============================================================================
# Cargar librerías
# install.packages("tidyverse")
# install.packages("ggplot2")
# install.packages("car")
# install.packages("agricolae")
# install.packages("DescTools")
# install.packages("patchwork")
# install.packages("broom")
library(tidyverse)
library(ggplot2)
library(car)
library(agricolae)
library(DescTools)
library(patchwork)
library(broom)
# Cargar datos limpios
# load("datos_limpios.RData")
# O si prefieres desde CSV:
datos <- read.csv("datos_limpios.csv", header = TRUE)
# ============================================================================
# SECCIÓN 1: ESTADÍSTICA DESCRIPTIVA
# ============================================================================
# ----------------------------------------------------------------------------
# 1.1 Tabla de estadísticos descriptivos por tratamiento
# ----------------------------------------------------------------------------
desc_por_tratamiento <- datos %>%
group_by(tratamiento) %>%
summarise(
n = n(),
media = mean(vitamina_c, na.rm = TRUE),
mediana = median(vitamina_c, na.rm = TRUE),
sd = sd(vitamina_c, na.rm = TRUE),
var = var(vitamina_c, na.rm = TRUE),
se = sd / sqrt(n), # Error estándar
cv = (sd / media) * 100, # Coeficiente de variación (%)
min = min(vitamina_c, na.rm = TRUE),
max = max(vitamina_c, na.rm = TRUE),
q1 = quantile(vitamina_c, 0.25, na.rm = TRUE),
q3 = quantile(vitamina_c, 0.75, na.rm = TRUE)
) %>%
mutate(
media = round(media, 2),
mediana = round(mediana, 2),
sd = round(sd, 2),
var = round(var, 2),
se = round(se, 2),
cv = round(cv, 2),
min = round(min, 2),
max = round(max, 2),
q1 = round(as.numeric(q1), 2),
q3 = round(as.numeric(q3), 2)
)
print("=== ESTADÍSTICOS DESCRIPTIVOS POR TRATAMIENTO ===")
print(desc_por_tratamiento)
# Guardar tabla descriptiva
library(writexl)
# write.csv(desc_por_tratamiento, "tabla_descriptiva.csv", row.names = FALSE)
write_xlsx(desc_por_tratamiento, "tabla_descriptiva.xlsx")
# ============================================================================
# SECCIÓN 2: GRÁFICOS EXPLORATORIOS
# ============================================================================
# ----------------------------------------------------------------------------
# 2.1 Boxplot con puntos individuales
# ----------------------------------------------------------------------------
p1 <- ggplot(datos, aes(x = tratamiento, y = vitamina_c, fill = tratamiento)) +
geom_boxplot(alpha = 0.6, outlier.shape = NA, width = 0.6) +
geom_jitter(width = 0.15, size = 2.5, alpha = 0.7, color = "black") +
stat_summary(fun = "mean", geom = "point", shape = 23, size = 4,
fill = "white", color = "black", stroke = 1.5) +
labs(title = "Contenido de Vitamina C por metodo de conservacion",
subtitle = "Los diamantes blancos indican la media de cada tratamiento",
x = "Tratamiento (Metodo de conservacion)",
y = "Vitamina C (mg/100mL)",
fill = "Tratamiento") +
theme_minimal(base_size = 12) +
theme(legend.position = "none",
plot.title = element_text(face = "bold", size = 14, hjust = 0.5),
plot.subtitle = element_text(size = 11, hjust = 0.5, color = "gray40"))
print(p1)
ggsave("grafico_boxplot.png", p1, width = 8, height = 6, dpi = 300)
# ----------------------------------------------------------------------------
# 2.2 Gráfico de medias con barras de error (SE)
# ----------------------------------------------------------------------------
p2 <- ggplot(desc_por_tratamiento,
aes(x = tratamiento, y = media, fill = tratamiento)) +
geom_bar(stat = "identity", width = 0.6, alpha = 0.7) +
geom_errorbar(aes(ymin = media - se, ymax = media + se),
width = 0.2, size = 0.8, color = "black") +
geom_text(aes(label = sprintf("%.2f \u00B1 %.2f", media, se)),
vjust = -1.5, size = 4, fontface = "bold") +
labs(title = "Medias y error estandar por tratamiento",
subtitle = "Barras de error: ± 1 SE",
x = "Tratamiento",
y = "Vitamina C (mg/100mL)",
fill = "Tratamiento") +
theme_minimal(base_size = 12) +
theme(legend.position = "none",
plot.title = element_text(face = "bold", size = 14, hjust = 0.5),
plot.subtitle = element_text(size = 11, hjust = 0.5, color = "gray40"),
axis.text.x = element_text(angle = 15, hjust = 1)) +
ylim(0, max(desc_por_tratamiento$media + desc_por_tratamiento$se) * 1.2)
print(p2)
ggsave("grafico_medias_se.png", p2, width = 8, height = 6, dpi = 300)
# ============================================================================
# SECCIÓN 3: MODELO MATEMÁTICO
# ============================================================================
cat("=== MODELO ESTADÍSTICO PARA DISEÑO COMPLETAMENTE ALEATORIZADO ===\n\n")
cat("El modelo lineal para un DCA con un factor se expresa como:\n\n")
cat("Y_ij = μ + τ_i + ε_ij\n\n")
cat("Donde:\n")
cat(" Y_ij = Contenido de vitamina C en la repetición j del tratamiento i\n")
cat(" μ = Media general del experimento\n")
cat(" τ_i = Efecto del tratamiento i (método de conservación)\n")
cat(" ε_ij = Error experimental ~ N(0, σ²)\n\n")
cat("Índices:\n")
cat(" i = 1, 2, 3 (tratamientos: T1, T2, T3)\n")
cat(" j = 1, 2, ..., 8 (repeticiones por tratamiento)\n\n")
cat("Hipótesis:\n")
cat(" H₀: τ₁ = τ₂ = τ₃ = 0 (no hay efecto de los tratamientos)\n")
cat(" H₁: al menos un τ_i ≠ 0 (al menos un tratamiento difiere)\n\n")
cat("Supuestos del modelo:\n")
cat(" 1. Independencia de las observaciones (garantizada por aleatorización)\n")
cat(" 2. Normalidad de los errores (ε_ij ~ N(0, σ²))\n")
cat(" 3. Homocedasticidad (varianza constante entre tratamientos)\n\n")
# ============================================================================
# SECCIÓN 4: VERIFICACIÓN DE SUPUESTOS
# ============================================================================
# ----------------------------------------------------------------------------
# Ajustar modelo ANOVA
# ----------------------------------------------------------------------------
modelo <- aov(vitamina_c ~ tratamiento, data = datos)
# Extraer residuos
residuos <- residuals(modelo)
valores_ajustados <- fitted(modelo)
# ----------------------------------------------------------------------------
# 4.1 Supuesto de Normalidad
# ----------------------------------------------------------------------------
cat("=== SUPUESTO 1: NORMALIDAD DE RESIDUOS ===\n\n")
# Prueba de Shapiro-Wilk
shapiro <- shapiro.test(residuos)
print(shapiro)
if(shapiro$p.value > 0.05){
cat("✓ INTERPRETACIÓN: No hay evidencia para rechazar normalidad (p > 0.05)\n")
cat(" Los residuos siguen una distribución normal.\n\n")
} else {
cat("✗ INTERPRETACIÓN: Los residuos NO siguen distribución normal (p ≤ 0.05)\n")
cat(" Considerar transformación de datos o prueba no paramétrica (Kruskal-Wallis).\n\n")
}
# Gráfico Q-Q
qq_plot <- ggplot(data.frame(residuos = residuos), aes(sample = residuos)) +
stat_qq(size = 3, color = "steelblue") +
stat_qq_line(color = "red", linewidth = 1.2) +
labs(title = "Grafico Q-Q de Residuos",
subtitle = paste("Shapiro-Wilk: p =", round(shapiro$p.value, 4)),
x = "Cuantiles teoricos",
y = "Cuantiles observados") +
theme_minimal(base_size = 12) +
theme(plot.title = element_text(face = "bold", size = 13),
plot.subtitle = element_text(size = 11, color = "gray40"))
ggsave("grafico_qq.png", qq_plot, width = 7, height = 6, dpi = 300)
print(qq_plot)
# Histograma de residuos con curva normal
hist_plot <- ggplot(data.frame(residuos = residuos), aes(x = residuos)) +
geom_histogram(aes(y = after_stat(density)), bins = 8,
fill = "steelblue", alpha = 0.6, color = "black") +
stat_function(fun = dnorm, args = list(mean = mean(residuos), sd = sd(residuos)),
color = "red", linewidth = 1.5) +
labs(title = "Histograma de residuos con curva normal superpuesta",
x = "Residuos",
y = "Densidad") +
theme_minimal(base_size = 12) +
theme(plot.title = element_text(face = "bold", size = 13))
ggsave("grafico_histograma_residuos.png", hist_plot, width = 7, height = 5, dpi = 300)
print(hist_plot)
# ----------------------------------------------------------------------------
# 4.2 Supuesto de Homocedasticidad (Homogeneidad de Varianzas)
# ----------------------------------------------------------------------------
# Prueba de Levene
levene <- leveneTest(vitamina_c ~ tratamiento, data = datos)
print(levene)
# Gráfico de residuos vs. valores ajustados
resid_plot <- ggplot(data.frame(ajustados = valores_ajustados, residuos = residuos),
aes(x = ajustados, y = residuos)) +
geom_point(size = 3, color = "steelblue", alpha = 0.7) +
geom_hline(yintercept = 0, linetype = "dashed", color = "red", linewidth = 1.2) +
geom_smooth(method = "lm", se = FALSE, color = "darkgreen", linewidth = 1) +
labs(title = "Residuos vs. valores ajustados",
subtitle = paste("Prueba de Levene: p =", round(levene[1, "Pr(>F)"], 4))) +
theme_minimal(base_size = 12)
ggsave("grafico_residuos_ajustados.png", resid_plot, width = 7, height = 6, dpi = 300)
print(resid_plot)
# ----------------------------------------------------------------------------
# 4.3 Supuesto de Independencia
# ----------------------------------------------------------------------------
cat("El supuesto de independencia se garantiza mediante:\n")
cat(" ✓ Aleatorización en la asignación de tratamientos\n")
cat(" ✓ Unidades experimentales independientes (envases separados)\n")
cat(" ✓ Mediciones realizadas en orden aleatorio\n")
cat(" ✓ Sin relación entre observaciones (cada envase es único)\n\n")
cat("Verificación gráfica: Gráfico de residuos vs. orden de ejecución\n")
# Agregar orden de ejecución a los datos
datos$orden <- 1:nrow(datos)
orden_plot <- ggplot(datos, aes(x = orden, y = residuos)) +
geom_point(size = 3, color = "steelblue", alpha = 0.7) +
geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
geom_smooth(method = "lm", se = TRUE, color = "darkgreen", alpha = 0.3) +
labs(title = "Residuos vs. orden de ejecucion",
subtitle = "Verificación de independencia temporal",
x = "Orden de ejecucion",
y = "Residuos") +
theme_minimal(base_size = 12) +
theme(plot.title = element_text(face = "bold", size = 13),
plot.subtitle = element_text(size = 11, color = "gray40"))
ggsave("grafico_residuos_orden.png", orden_plot, width = 8, height = 5, dpi = 300)
print(orden_plot)
# ============================================================================
# SECCIÓN 5: ANÁLISIS DE VARIANZA (ANOVA)
# ============================================================================
# ----------------------------------------------------------------------------
# 5.1 Tabla ANOVA
# ----------------------------------------------------------------------------
tabla_anova <- anova(modelo)
print(tabla_anova)
# Guardar tabla ANOVA
write_xlsx(as.data.frame(tabla_anova), "tabla_anova.xlsx")
p_valor <- anova(modelo)[1, "Pr(>F)"]
if(p_valor < 0.05){
cat("✓ CONCLUSION: Hay diferencias significativas entre tratamientos (p < 0.05)\n")
cat(" Rechazamos H₀. Al menos un tratamiento difiere de los demas.\n")
cat(" Procedemos a pruebas de comparacion multiple de medias.\n\n")
} else {
cat("✗ CONCLUSION: No hay diferencias significativas entre tratamientos (p ≥ 0.05)\n")
cat(" No rechazamos H₀. Los tratamientos no difieren estadisticamente.\n")
cat(" Las pruebas post-hoc no son necesarias pero pueden realizarse de forma exploratoria.\n\n")
}
# ============================================================================
# SECCIÓN 6: PRUEBA DE COMPARACIÓN MÚLTIPLE DE MEDIAS
# ============================================================================
# ----------------------------------------------------------------------------
# 6.1 Prueba de Tukey (HSD - Honest Significant Difference)
# ----------------------------------------------------------------------------
cat("=== PRUEBA DE TUKEY (HSD) ===\n")
cat("Propósito: Comparaciones por pares entre TODOS los tratamientos\n")
cat("Controla la tasa de error familiar (FWER) para múltiples comparaciones\n\n")
tukey <- TukeyHSD(modelo, "tratamiento", conf.level = 0.95)
print(tukey)
# Interpretación
cat("\n=== INTERPRETACIÓN DE TUKEY ===\n\n")
tukey_df <- as.data.frame(tukey$tratamiento)
tukey_df$significativo <- ifelse(tukey_df$`p adj` < 0.05, "Sí *", "No")
print(tukey_df)
cat("\nRegla de decisión:\n")
cat(" - Si p-valor ajustado < 0.05 → Diferencia significativa\n")
cat(" - Si p-valor ajustado ≥ 0.05 → Diferencia no significativa\n\n")Después de ejecutar todos los scripts, tendrás:
01_calculo_poder.R - Análisis de potencia y tamaño
muestral02_aleatorizacion.R - Aleatorización del diseño03_importar_datos.R - Importación y limpieza de
datos04_analisis_completo.R - Análisis estadísticoplanilla_aleatorizacion.xlsx - Para imprimir y usar en
campodatos_vitamina_c.xlsx - Datos cargados desde Exceldatos_limpios.RData y datos_limpios.csv -
Datos procesadosresultados_analisis.RData - Objetos de resultadostabla_descriptiva.csv, tabla_anova.csv,
etc.esquema_aleatorizacion.png - Esquema visual de
asignacióngrafico_boxplot.png - Boxplot con puntosgrafico_medias_se.png - Medias con error estándargrafico_qq.png - Q-Q plot de normalidadgrafico_histograma_residuos.png - Histograma con curva
normalgrafico_residuos_ajustados.png - Residuos
vs. ajustadosgrafico_residuos_orden.png - Residuos vs. ordenset.seed() para reproducibilidadDurante esta clase, pueden usar IA para:
✅ Buscar fórmulas de cálculo de potencia ✅ Entender la interpretación de pruebas estadísticas ✅ Depurar errores en código R ✅ Mejorar redacción del informe
❌ NO usar IA para:
Pregunta reflexiva para el informe: ¿Cómo validaron que el código generado o sugerido por IA era correcto? ¿Qué errores encontraron y cómo los resolvieron?
¡Éxitos en su análisis! Recuerden que la estadística es una herramienta para responder preguntas científicas, no un fin en sí misma. 📊🔬
https://search.r-project.org/CRAN/refmans/pwr/html/pwr.anova.test.html↩︎
http://rstudio-pubs-static.s3.amazonaws.com/235504_3b3b13081ea24c06ad7642f96dba5242.html↩︎
https://stevenggoni.github.io/tutoriales_R/potencia_modelos.html↩︎
https://r-statistics.co/Statistical-Power-Analysis-in-R.html↩︎
https://jdleongomez.info/en/publication/leongomez2020power/↩︎
https://cran.r-project.org/web/packages/pwranova/pwranova.pdf↩︎
https://cran.r-project.org/web/packages/pwr/vignettes/pwr-vignette.html↩︎
https://cran.r-project.org/web/packages/powertools/powertools.pdf↩︎
https://cran.rstudio.com/web/packages/pwrss/vignettes/examples.html↩︎
https://zenodo.org/records/3988777/files/Análisis de Poder en R.pdf↩︎