1 Diseño de Experimentos: Teoría Básica

1.1 1. ¿Qué es un Diseño de Experimentos?

Un Diseño de Experimentos (DE) es un conjunto de técnicas y metodologías que permiten planificar, conducir, analizar y extraer conclusiones de experimentos de forma eficiente y rigurosa. El objetivo principal es identificar y cuantificar el efecto de uno o varios factores (variables independientes) sobre una o varias respuestas (variables dependientes), minimizando la variabilidad no explicada y optimizando recursos.

1.1.1 Aspectos clave

  • Objetivo claro: Definir con precisión cuál es la pregunta o hipótesis que se desea responder.
  • Control de variabilidad: Reducir el ruido o variabilidad no deseada mediante bloqueos, réplicas o diseños adecuados.
  • Eficiencia: Diseñar la experimentación para obtener la máxima información con el menor número de ensayos posible.
  • Validez estadística: Asegurar que los resultados sean válidos (internamente y externamente) mediante aleatorización, réplica y control de sesgos.

1.2 2. Usos y Aplicaciones

Los DE se emplean en numerosas áreas donde se requiere optimizar procesos, validar hipótesis o mejorar la toma de decisiones: - Industria y manufactura: Ajustar parámetros de procesos de producción para maximizar rendimiento o minimizar defectos. - Agricultura: Evaluar tratamientos (fertilizantes, variedades de cultivo, riegos) para optimizar rendimiento de cosechas. - Farmacéutica y biomedicina: Determinar dosis óptimas, eficacia de tratamientos, condiciones de cultivo de células, etc. - Investigación científica: Probar hipótesis en ciencias naturales, psicología, sociología, etc. - Marketing y negocios: Experimentos A/B para lanzar campañas publicitarias, comparar versiones de una web o producto. - Ingeniería: Optimización de diseños de componentes, condiciones de prueba en sistemas mecánicos o electrónicos. - Tecnologías de la información: Experimentos en usabilidad, rendimiento de algoritmos bajo distintas configuraciones.

1.3 3. Elementos de un Diseño de Experimentos

Un DE se compone de varios elementos esenciales:

  1. Factor(es)
    • Definición: Variables independientes que se manipulan intencionalmente (p. ej., temperatura, dosis, método).
    • Niveles: Cada factor tiene niveles o valores específicos (p. ej., temperatura baja/media/alta).
  2. Respuesta(s)
    • Definición: Variables dependientes que se miden para evaluar el efecto de los factores (p. ej., rendimiento, tiempo de proceso, tasa de éxito).
    • Métricas: Debe definirse cómo se mide: escalas, precisión, instrumentos.
  3. Tratamientos
    • Combinaciones específicas de niveles de factores que se aplican en cada unidad experimental.
  4. Unidades experimentales
    • Elementos a los cuales se aplica cada tratamiento (p. ej., plantas en un campo, lotes de producción, usuarios en un A/B test).
    • Importancia del aleatorizado: asignar tratamientos de manera aleatoria para evitar sesgos sistemáticos.
  5. Réplicas
    • Repeticiones del mismo tratamiento en unidades experimentales distintas, para estimar la variabilidad inherente y obtener mayor precisión.
  6. Aleatorización
    • Proceso de asignar aleatoriamente tratamientos a las unidades experimentales, garantizando independencia y reduciendo sesgos.
  7. Bloqueo / Estratificación
    • Agrupar unidades experimentales en bloques con características homogéneas (p. ej., lotes de materia prima, días de prueba, operarios distintos) para controlar variabilidad conocida.
    • Dentro de cada bloque se aleatorizan los tratamientos.
  8. Control de factores extraños
    • Identificación de variables externas que pueden afectar la respuesta y, de ser posible, controlarlas o bloquearlas (p. ej., hora del día, lote de materia prima).
  9. Diseño de la muestra y tamaño muestral
    • Calcular cuántas réplicas son necesarias para detectar efectos de cierto tamaño con la potencia deseada.
    • Equilibrar presupuesto, tiempo y recursos con la precisión estadística.
  10. Análisis estadístico
    • Selección de modelos adecuados (ANOVA, regresión lineal, modelos mixtos si hay estructura de aleatoriedad adicional).
    • Verificación de supuestos (normalidad de residuos, homocedasticidad, independencia).
    • Estimación de efectos principales, interacciones, tamaños de efecto.
    • Pruebas post-hoc o comparaciones múltiples si aplica.
    • Visualización de resultados para interpretar tendencias y patrones.
  11. Interpretación y conclusiones
    • Relacionar los hallazgos con la pregunta original.
    • Considerar validez interna (diseño bien ejecutado) y externa (generalización de resultados).
    • Proponer recomendaciones o próximos pasos basados en evidencias.

1.4 4. Ejemplo Conceptual Práctico (Teórico)

Ejemplo sencillo en agricultura
- Pregunta: ¿Cómo afecta la dosis de fertilizante (Factor A) y el tipo de riego (Factor B) al rendimiento de un cultivo?
- Factores y niveles:
- Factor A (Fertilizante): Nivel 1 = baja dosis, Nivel 2 = dosis media, Nivel 3 = alta dosis.
- Factor B (Riego): Nivel 1 = riego diario, Nivel 2 = riego cada dos días.
- Unidades experimentales: Parcela de cultivo individual.
- Bloques: Si el terreno tiene ligera pendiente o variabilidad de suelo, dividir en bloques homogéneos (por ej., secciones con similar composición de suelo).
- Aleatorización: Dentro de cada bloque, asignar aleatoriamente cada combinación de dosis × riego a parcelas.
- Réplicas: Repetir cada combinación en varias parcelas (p. ej., 4 réplicas) para estimar variabilidad.
- Análisis:
- Modelo ANOVA de dos factores con interacción: Rendimiento ~ Fertilizante + Riego + Fertilizante:Riego + error.
- Verificar si hay efecto principal de dosis, de riego, o interacción (p. ej., quizá la dosis alta mejora rendimiento solo con riego cada dos días, pero no con riego diario).
- Interpretación:
- Si interacción significativa, investigar combinaciones específicas (ej. gráficos de interacción).
- Calcular tamaño de efecto (por ejemplo, η² parciales) para cuantificar magnitud práctica.
- Recomendaciones: dosis y frecuencia de riego óptimos para maximizar rendimiento.

Nota: Este ejemplo es conceptual: más adelante, se traduciría a código R que genere el plan experimental, realice ANOVA y visualice interacciones.

1.5 5. Tipos de Diseños de Experimentos: Usos, Ventajas y Desventajas

A continuación se listan los tipos más comunes de diseños, con una breve descripción, cuándo se usan, ventajas y desventajas:

  1. Diseño Completamente Aleatorizado (DCA / CRD)

    • Descripción: Todas las unidades experimentales se asignan aleatoriamente a los tratamientos sin estructura adicional.
    • Usos: Cuando las unidades experimentales son homogéneas o la variabilidad entre ellas es pequeña o irrelevante.
    • Ventajas:
      • Simple de planificar y analizar.
      • Máxima aleatorización, minimiza sesgos.
      • Flexible: permite añadir o quitar tratamientos fácilmente.
    • Desventajas:
      • Si hay variabilidad no controlada entre unidades, puede aumentar la varianza residual y reducir potencia.
      • No controla factores extraños conocidos (a menos que se incluyan en el modelo como covariables).

  2. Diseño en Bloques Aleatorizados (DBA / RBD)

    • Descripción: Unidades experimentales se agrupan en bloques homogéneos. Dentro de cada bloque, se asignan aleatoriamente los tratamientos.
    • Usos: Cuando existe una fuente de variación conocida (suelo, lotes, operarios, tiempo) que se quiere controlar.
    • Ventajas:
      • Reduce variabilidad residual al controlar la variabilidad entre bloques.
      • Aumenta potencia estadística con mismo número de réplicas.
    • Desventajas:
      • Requiere identificar y formar bloques homogéneos.
      • Puede ser más complejo de planificar si muchos bloques o si tratamientos no caben uniformemente en bloques.

  3. Diseño Factorial Completo

    • Descripción: Permite estudiar varios factores simultáneamente, probando todas las combinaciones de niveles (p. ej., 2×3 factorial).
    • Usos: Cuando se desea investigar efectos principales y posibles interacciones entre factores.
    • Ventajas:
      • Eficiencia: con un número moderado de ensayos, se obtiene información de múltiples factores.
      • Permite detectar interacciones, esencial en procesos complejos.
    • Desventajas:
      • El número de tratamientos crece exponencialmente con factores/niveles grandes.
      • Interpretación de interacciones de alto orden puede ser complicada.
    • Variantes:
      • Factorial fraccional: solo una fracción de combinaciones, útil cuando hay muchos factores y se busca screening para identificar factores relevantes.
        • Ventaja: Menos ensayos.
        • Desventaja: Puede confundir (confounding) efectos si no se diseña cuidadosamente.
      • Factorial completo con bloques: combina con bloque aleatorizado para controlar variabilidad adicional.

    Diseño de Cuadrado Latino

    • Descripción: Dos fuentes de variabilidad (dos factores de bloqueo) se controlan: se organiza en una matriz donde filas y columnas son bloques, y cada tratamiento aparece exactamente una vez en cada fila y columna.
    • Usos: Cuando hay dos fuentes de variabilidad conocidas (p. ej., posición en filas y columnas en campo agrícola, día y operario en laboratorio).
    • Ventajas:
      • Controla dos factores de variabilidad sin requerir tantas réplicas como en un diseño factorial completo de bloque doble.
    • Desventajas:
      • Número de tratamientos = tamaño de la cuadrícula; obliga a que el número de tratamientos concuerde con número de filas/columnas.
      • Si los bloques no son perfectamente homogéneos, puede quedar variabilidad no controlada.
      • No permite examinar interacción entre bloques.

    Diseño de Bloques Incompletos Balanceados (BIBD)

    • Descripción: No todos los tratamientos se prueban en cada bloque; los bloques son “incompletos” pero balanceados según ciertas condiciones.
    • Usos: Cuando el número de tratamientos es grande y no es factible incluir todos en cada bloque (p. ej., estudios clínicos con muchos tratamientos, ensayos con muchas variedades).
    • Ventajas:
      • Permite comparar muchos tratamientos controlando variabilidad por bloques sin exigir bloques grandes.
    • Desventajas:
      • Planificación matemática más compleja.
      • Requiere balanceo específico para obtener comparaciones precisas.

    Diseño en Parcelas Divididas (Split-Plot)

    • Descripción: Algunos factores (plot) se asignan a “parcelas principales” y otros factores (subplot) se asignan dentro de esas parcelas. Diseñado para cuando algunos factores son difíciles o costosos de variar a nivel pequeño.
    • Usos: Cuando manipular ciertos factores a nivel de unidad pequeña es costoso o poco práctico (p. ej., tratamiento principal en campo, subtratamiento en subparcelas).
    • Ventajas:
      • Flexibilidad para manejar factores con niveles de manipulación distintos (difícil vs. fácil).
    • Desventajas:
      • Análisis más complejo: errores experimentales de parcelas principales y subparcelas diferencian niveles de variabilidad.
      • Menor potencia para el factor aplicado a parcelas principales (menos réplicas efectivas).

  4. Diseño Secuencial o en Fases

    • Descripción: El experimento se realiza en etapas, usando resultados preliminares para ajustar el diseño o concentración de factores en etapas posteriores.
    • Usos: Cuando se desea explorar inicialmente un amplio rango de niveles, luego refinar la región de interés (experimentos de respuesta superficial, por ejemplo).
    • Ventajas:
      • Eficiencia: evita gastar muchos recursos en rangos no informativos.
      • Permite aprendizaje y adaptación entre fases.
    • Desventajas:
      • Planificación más dinámica y exigente.
      • Posible complejidad en análisis al combinar fases (diferentes errores, variabilidad entre fases).

    Diseños de Respuesta de Superficie (RSM)

    • Descripción: Serie de diseños factoriales (o fraccionados) seguidos de diseños adecuados (por ejemplo, diseños centrales compuestos o Box-Behnken) para modelar la respuesta como función continua de factores y ubicar óptimos.
    • Usos: Optimización de procesos donde la relación entre factores y respuesta es aproximadamente cuadrática o no lineal suave.
    • Ventajas:
      • Permite modelar superficies de respuesta y encontrar combinaciones óptimas.
    • Desventajas:
      • Requiere suposiciones sobre continuidad y comportamientos cuadráticos.
      • Puede necesitar más réplicas y análisis de diagnóstico de modelo.

    Diseño de Experimentos Adaptativo (Adaptive Designs)

    • Descripción: En campos como ensayos clínicos, el diseño se adapta durante la recolección de datos (por ejemplo, reasignar pacientes a tratamientos más prometedores).
    • Usos: Ensayos clínicos, ensayos en marketing online, A/B testing continuo.
    • Ventajas:
      • Ética y eficiencia: en clínica, menos pacientes reciben tratamientos menos efectivos; en marketing, despliegue más rápido del mejor contenido.
    • Desventajas:
      • Complejidad estadística y regulatoria (en clínica) para asegurar validez y control de error tipo I.
      • Necesita monitoreo interino y reglas predefinidas.

  5. Diseños de Experimentos con Modelos Mixtos

  • Descripción: Cuando existen efectos aleatorios (p. ej., variabilidad entre sujetos, lotes, etc.) junto con efectos fijos de interés. Utiliza modelos lineales mixtos para análisis.
  • Usos: Datos jerárquicos o anidados (por ejemplo, ensayos en bloques, mediciones repetidas en sujetos, experimentos en tiempo).
  • Ventajas:
    • Modela estructura de correlación y variabilidad entre unidades más complejas.
    • Mejora estimaciones al considerar componentes de varianza.
  • Desventajas:
    • Análisis más complejo, requiere entendimiento de modelos mixtos y supuestos.
    • Puede requerir software especializado o mayor cantidad de datos para estimar componentes.

1.6 6. Consideraciones Prácticas y Buenas Prácticas

  • Definir hipótesis claras: Antes de diseñar, formular hipótesis nulas y alternativas (por ejemplo, “la dosis no afecta el rendimiento” vs. “al menos una dosis difiere”).
  • Tamaño de muestra y potencia: Realizar análisis de potencia para determinar réplicas necesarias según el tamaño de efecto esperado y la variabilidad estimada.
  • Aleatorización rigurosa: Usar métodos aleatorios apropiados (generadores confiables) y documentar el proceso.
  • Control de sesgos: Ciego o doble ciego en experimentos clínicos o estudios de comportamiento cuando aplique.
  • Registro detallado: Documentar fecha, hora, condiciones, desviaciones o problemas durante el experimento.
  • Verificación de supuestos: Tras recabar datos, verificar normalidad de residuos, homogeneidad de varianzas, independencia. Si no se cumplen, considerar transformaciones o métodos no paramétricos.
  • Visualización exploratoria: Boxplots, gráficos de interacción, scatterplots, para entender patrones antes y después de ANOVA.
  • Análisis de efectos prácticos: Más allá de significancia estadística (p-valor), calcular tamaños de efecto (η², ω², d de Cohen, etc.) y analizar relevancia práctica.
  • Repetibilidad y reproducibilidad: Diseñar de modo que otros puedan replicar el experimento; compartir protocolos, datos y código.
  • Costos y recursos: Balancear alcance del experimento con costos (tiempo, material, personal). A veces un diseño fraccionado o secuencial es preferible.
  • Interpretación consciente: Reconocer limitaciones del diseño (alcance de inferencia) y evitar sobreinterpretar resultados.

1.7 7. Resumen de Tipos y Comparación

Diseño Uso principal Ventaja principal Desventaja principal
Completamente Aleatorizado Unidades homogéneas, sin fuente de variabilidad extra Simplicidad y máxima aleatorización No controla variabilidad conocida
Bloque Aleatorizado Controlar variabilidad conocida entre bloques Reduce variabilidad residual Requiere identificar bloques, más complejo
Factorial Completo Estudiar múltiples factores e interacciones Eficiencia en detección de efectos e interacciones Crece rápido con muchos factores
Factorial Fraccional Screening de muchos factores Menos ensayos iniciales Confounding de efectos si no se planifica bien
Cuadrado Latino Controlar dos fuentes de variabilidad Control eficiente de dos factores de bloqueo Número de tratamientos fijo, rigidez en diseño
Split-Plot Factores con diferentes escalas de manipulación Maneja factores difíciles de cambiar Análisis complejo, menor potencia para ciertos factores
Respuesta de Superficie (RSM) Optimización de respuestas continuas Modela superficie y encuentra óptimos Requiere supuestos de continuidad, más ensayos
Secuencial / Adaptativo Exploración y refinamiento por fases Eficiencia y adaptación Complejidad en planificación y análisis
Bloques Incompletos Balanceados Muchos tratamientos con bloques pequeños Controla variabilidad sin bloques grandes Planificación matemática compleja
Modelos Mixtos Datos jerárquicos o estructurados Modela correlaciones y variabilidad compleja Requiere experiencia en modelos avanzados

1.8 8. Ejemplos Claros y Sencillos (Conceptuales)

  1. Test de A/B en marketing digital
    • Factor: versión de la página web (A vs. B).
    • Unidades: usuarios asignados aleatoriamente.
    • Réplicas: muchos usuarios.
    • Métrica: tasa de conversión.
    • Diseño: completamente aleatorizado; si hay segmentos (bloques) como país o dispositivo, se puede bloquear o estratificar asignación.
    • Análisis: prueba de proporciones, ANOVA o regresión logística según enfoque.
  2. Ensayo clínico fase II
    • Factor: dosis de fármaco (baja, media, alta).
    • Bloqueo: centro hospitalario (varios hospitales).
    • Unidades: pacientes.
    • Réplicas: pacientes por dosis en cada centro.
    • Diseño: bloque aleatorizado dentro de cada centro.
    • Métrica: respuesta biométrica (p. ej., reducción de síntoma).
    • Análisis: ANOVA o modelos mixtos con efecto aleatorio de centro, verificación de supuestos, análisis de potencia y ajuste múltiple si se comparan varias dosis.
  3. Optimización de proceso industrial
    • Factores: temperatura (3 niveles), presión (2 niveles), velocidad de avance (2 niveles).
    • Diseño factorial completo 3×2×2 = 12 combinaciones.
    • Bloque: turnos de producción o lotes de materia prima.
    • Réplicas: repetir cada combinación varias veces.
    • Métrica: calidad del producto (valor continuo) o tasa de defectos.
    • Análisis: ANOVA factorial, gráficos de interacción, cálculo de tamaños de efecto, recomendación de parámetros óptimos.
  4. Experimento agrícola sencillo (CRD vs. RBD)
    • Si campo homogéneo: usar DCA.
    • Si heterogeneidad en terreno: usar bloques.
    • Factor: variedad de semilla (4 variedades).
    • Diseño:
      • DCA: 4 variedades asignadas aleatoriamente a parcelas.
      • RBD: dividir campo en bloques según calidad de suelo, asignar aleatoriamente variedades dentro de cada bloque.
    • Métrica: rendimiento.
    • Análisis: comparar varianza residual y poder detectivo en ambos diseños.

1.9 9. Próximos Pasos

  • Profundizar en análisis de modelos: cómo realizar ANOVA en R, verificar supuestos, interpretar salidas.
  • Visualización: gráficos de caja, interacción, diagnóstico de residuos.
  • Cálculo de potencia y tamaño de muestra en R: ejemplos prácticos.
  • Diseños factoriales y fraccionales en R: generar planes, analizar con paquetes como FrF2, AlgDesign.
  • Modelos mixtos: uso de lme4, interpretar componentes de varianza.
  • Diseños de superficie de respuesta: uso de rsm en R.
  • Ejemplos reales: desarrollar un caso completo desde la formulación de hipótesis hasta la interpretación final.

Nota: Esta sección está enfocada en la teoría y descripción de conceptos, empleando ejemplos conceptuales sin código R. En la siguiente etapa, traduciremos estos conceptos a scripts en R que generen planes experimentales, realicen análisis ANOVA, verifiquen supuestos y produzcan visualizaciones.

2 Ejercicios

2.1 Ejemplo de Diseño Completamente Aleatorizado (DCA) usando un dataset existente en R

https://www.youtube.com/watch?v=fWVacB2_JOs

2.2 Ejemplo Mejorado: Análisis de One-Way ANOVA para Tratamientos de Remojo de Frijoles

A continuación se presenta una versión refinada del enunciado y el análisis paso a paso en R, con formato claro en Markdown. Se incluye la descripción del experimento, la estructura de los datos, la carga y análisis en R, verificación de supuestos, comparaciones múltiples, tamaños de efecto y visualizaciones. El objetivo es comparar los cuatro tratamientos (Control, T2, T3, T4) sobre el tiempo de cocción de frijoles.

2.2.1 1. Enunciado Refinado

En un centro de investigación se realiza un experimento para comparar cuatro tratamientos aplicados previamente a frijoles crudos, con el fin de reducir su tiempo de cocción. Los tratamientos son:

  1. Control: No se aplica remojo ni aditivo.
  2. T2 (Bicarbonato): Remojo en agua con bicarbonato de sodio (NaHCO₃).
  3. T3 (Sal): Remojo en agua con cloruro de sodio (NaCl).
  4. T4 (Mezcla): Remojo en agua con una combinación de bicarbonato de sodio y sal común en proporciones iguales.
  • Variable de respuesta: Tiempo de cocción (en minutos) medido para cada lote de frijoles tras el respectivo tratamiento.
  • Unidades experimentales: Lotes homogéneos de frijoles, asignados de forma completamente aleatorizada a uno de los cuatro tratamientos.
  • Número de réplicas: 7 mediciones por tratamiento (total 28 observaciones).
  • Objetivo: Determinar si alguno de los tratamientos reduce significativamente el tiempo de cocción en comparación con los demás, y cuantificar la magnitud de las diferencias.

Nota: Asumimos que las condiciones de cocción (temperatura, presión, tipo de olla, volumen de agua, etc.) se mantienen constantes en todas las réplicas, de modo que la fuente principal de variabilidad proviene del tratamiento de remojo.

2.2.2 2. Datos Observados

A continuación se muestran los tiempos de cocción (en minutos), con 7 observaciones por cada tratamiento:

Tratamiento Observaciones de tiempo (min)
Control 213, 214, 204, 208, 212, 200, 207
T2 (Bicarbonato) 76, 85, 74, 78, 82, 75, 82
T3 (Sal) 57, 67, 55, 64, 61, 63, 63
T4 (Mezcla) 84, 82, 85, 92, 87, 79, 90

2.2.3 3. Preparación del Dataset en R

Se construye un data.frame con dos columnas:
- Tratamiento: factor con niveles Control, T2, T3, T4.
- Tiempo: numérico con el tiempo de cocción correspondiente.

# 3.1. Crear el data.frame manualmente
tiempos_control <- c(213, 214, 204, 208, 212, 200, 207)
tiempos_T2      <- c(76, 85, 74, 78, 82, 75, 82)
tiempos_T3      <- c(57, 67, 55, 64, 61, 63, 63)
tiempos_T4      <- c(84, 82, 85, 92, 87, 79, 90)

df <- data.frame(
  Tratamiento = factor(rep(c("Control", "T2", "T3", "T4"), each = 7),
                       levels = c("Control", "T2", "T3", "T4")),
  Tiempo = c(tiempos_control, tiempos_T2, tiempos_T3, tiempos_T4)
)

# 3.2. Inspección inicial
str(df)
## 'data.frame':    28 obs. of  2 variables:
##  $ Tratamiento: Factor w/ 4 levels "Control","T2",..: 1 1 1 1 1 1 1 2 2 2 ...
##  $ Tiempo     : num  213 214 204 208 212 200 207 76 85 74 ...
table(df$Tratamiento)  # Debe mostrar 7 observaciones de cada nivel
## 
## Control      T2      T3      T4 
##       7       7       7       7
summary(df)
##   Tratamiento     Tiempo      
##  Control:7    Min.   : 55.00  
##  T2     :7    1st Qu.: 72.25  
##  T3     :7    Median : 82.00  
##  T4     :7    Mean   :108.54  
##               3rd Qu.:119.00  
##               Max.   :214.00
# Estadísticas descriptivas por grupo:
res <- aggregate(Tiempo ~ Tratamiento, data = df, 
                 FUN = function(x) c(media = mean(x), sd = sd(x), n = length(x)))

res2 <- data.frame(
  Tratamiento = res$Tratamiento,
  media       = res$Tiempo[, "media"],
  sd          = res$Tiempo[, "sd"],
  n           = res$Tiempo[, "n"]
)
# Ahora res2 es un data.frame con columnas claras.
print(res2)
##   Tratamiento     media       sd n
## 1     Control 208.28571 5.122313 7
## 2          T2  78.85714 4.180453 7
## 3          T3  61.42857 4.157609 7
## 4          T4  85.57143 4.503967 7
library(ggplot2)
ggplot(df, aes(x = Tratamiento, y = Tiempo, fill = Tratamiento)) +
  geom_boxplot(alpha = 0.7) +
  ggtitle("Boxplot: Tiempo de Cocción según Tratamiento") +
  xlab("Tratamiento") +
  ylab("Tiempo de cocción (min)") +
  theme_minimal() +
  theme(legend.position = "none")

como podemos Evaluar en nuestra estadísticas descriptivas, las desviaciones estándar de los diferentes tratamientos, son sillares y podemos notar que el tratamiento T3 es el que menos tiempo toma de Cocción, pero sera sistemáticamente así, vamos a detallarlo, con nuestra tabla anova.

#Ajuste ANOVA
modelo <- aov(Tiempo ~ Tratamiento, data = df)
anova_res <- summary(modelo)
print(anova_res)
##             Df Sum Sq Mean Sq F value Pr(>F)    
## Tratamiento  3  95041   31680    1559 <2e-16 ***
## Residuals   24    488      20                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Verificación de supuestos
residuos <- residuals(modelo)
fitted_vals <- fitted(modelo)

# QQ-plot y Shapiro-Wilk
par(mfrow = c(1, 2))
qqnorm(residuos, main = "QQ-Plot de residuos")
qqline(residuos, col = "red")
shapiro_res <- shapiro.test(residuos)
print(shapiro_res)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuos
## W = 0.95991, p-value = 0.3469
# Residuals vs Fitted y Levene
plot(fitted_vals, residuos,
     xlab = "Valores ajustados", ylab = "Residuos",
     main = "Residuals vs Fitted")
abline(h = 0, col = "red")

if (!requireNamespace("car", quietly = TRUE)) {
  install.packages("car")
}
library(car)
## Cargando paquete requerido: carData
levene_res <- car::leveneTest(Tiempo ~ Tratamiento, data = df)
print(levene_res)
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  3  0.1631 0.9201
##       24
par(mfrow = c(1, 1))

En nuestra tabla anova notamos que existe una diferencia significativa, entre cada nivel del factor llamado tratamiento, en este experimento tenemos un solo factor y 4 niveles y 7 replicas por nivel y ningún bloqueo.

Los supuestos de la anova se cumple.

  • Normalidad de los residuos cumplida

  • Homogeneidad de varianzas cumplida

Vamos hacer prueba post hoc para validar que medias son diferentes entre los diferentes niveles.

Comparaciones Múltiples (Post-hoc)
Si ANOVA muestra efecto significativo de Tratamiento, se procede a identificar pares que difieren.

anova_tab <- anova_res[[1]]
p_val <- anova_tab["Tratamiento", "Pr(>F)"]
if (p_val < 0.05) {
  tukey_res <- TukeyHSD(modelo, "Tratamiento")
  print(tukey_res)
}
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Tiempo ~ Tratamiento, data = df)
## 
## $Tratamiento
##                   diff           lwr        upr     p adj
## T2-Control -129.428571 -136.07568671 -122.78146 0.0000000
## T3-Control -146.857143 -153.50425813 -140.21003 0.0000000
## T4-Control -122.714286 -129.36140099 -116.06717 0.0000000
## T3-T2       -17.428571  -24.07568671  -10.78146 0.0000010
## T4-T2         6.714286    0.06717044   13.36140 0.0471059
## T4-T3        24.142857   17.49574187   30.78997 0.0000000
plot(tukey_res)

Con esta prueba podemos validar que todos los tratamientos, son diferentes, pero entre T4-T2 existe una similitud, por 0.03 queda acteado que tiene medias iguales.

Vamos a validar con un test diferente a TukeyHSD.

duncan_res <- agricolae::duncan.test(modelo, "Tratamiento", 
                                     group = TRUE, 
                                     console = TRUE)
## 
## Study: modelo ~ "Tratamiento"
## 
## Duncan's new multiple range test
## for Tiempo 
## 
## Mean Square Error:  20.32143 
## 
## Tratamiento,  means
## 
##            Tiempo      std r       se Min Max   Q25 Q50   Q75
## Control 208.28571 5.122313 7 1.703837 200 214 205.5 208 212.5
## T2       78.85714 4.180453 7 1.703837  74  85  75.5  78  82.0
## T3       61.42857 4.157609 7 1.703837  55  67  59.0  63  63.5
## T4       85.57143 4.503967 7 1.703837  79  92  83.0  85  88.5
## 
## Alpha: 0.05 ; DF Error: 24 
## 
## Critical Range
##        2        3        4 
## 4.973148 5.223301 5.383910 
## 
## Means with the same letter are not significantly different.
## 
##            Tiempo groups
## Control 208.28571      a
## T4       85.57143      b
## T2       78.85714      c
## T3       61.42857      d
plot(duncan_res)

Podemos observar que todos los tratamientos, tienen grupos de clasificación diferentes, lo podemos observar en que tienen letra a,b,c,d y qyue el tiempo mas bajo es T3

Vamos a medir el efecto del tratamiento

# Tamaños de efecto
# η² parcial:
eta2_res <- effectsize::eta_squared(modelo, partial = TRUE)
## For one-way between subjects designs, partial eta squared is equivalent
##   to eta squared. Returning eta squared.
print(eta2_res)
## # Effect Size for ANOVA
## 
## Parameter   | Eta2 |       95% CI
## ---------------------------------
## Tratamiento | 0.99 | [0.99, 1.00]
## 
## - One-sided CIs: upper bound fixed at [1.00].
# Salida típica: muestra términos del modelo y su η² parcial.

# ω²:
omega2_res <- effectsize::omega_squared(modelo)
## For one-way between subjects designs, partial omega squared is
##   equivalent to omega squared. Returning omega squared.
print(omega2_res)
## # Effect Size for ANOVA
## 
## Parameter   | Omega2 |       95% CI
## -----------------------------------
## Tratamiento |   0.99 | [0.99, 1.00]
## 
## - One-sided CIs: upper bound fixed at [1.00].

Al solo existir un solo factor nos da igual

Eta2 (η² parcial):
- Valor: 0.99 (o 99%).
- Significado: indica que aproximadamente el 99% de la variabilidad total explicable (varianza explicada más varianza residual) en la variable respuesta se atribuye al factor “Tratamiento”. En otras palabras, el efecto del tratamiento en tu experimento es extremadamente grande, pues casi toda la variación observada en los datos se debe a las diferencias entre niveles de “Tratamiento”.
- Rango habitual: η² varía de 0 a 1. Valores cercanos a 0 indican efecto pequeño; valores cercanos a 1, efecto muy grande. Aquí, 0.99 es casi el máximo posible, lo que sugiere un efecto muy dominante.

  1. Omega2 (ω²):

    • Valor: 0.99 (o 99%).

    • Significado: ω² es una estimación más conservadora de la proporción de varianza en la población explicada por el factor, corrigiendo de cierto modo el sesgo de η² hacia sobreestimación. También varía de 0 a 1. Un ω² de 0.99 sugiere que, incluso tras corrección, el factor “Tratamiento” explica prácticamente toda la varianza explicable en la población.

library(ggplot2)
ggplot(df, aes(x = Tratamiento, y = Tiempo, fill = Tratamiento)) +
  geom_violin(alpha = 0.5) +
  geom_jitter(width = 0.1, size = 2) +
  ggtitle("Distribución de Tiempo de Cocción por Tratamiento") +
  xlab("Tratamiento") +
  ylab("Tiempo (min)") +
  theme_minimal() +
  theme(legend.position = "none")

2.3 Diseño de bloques completos al azar (DBCA)

https://rpubs.com/deliosalgado/1168360

https://www.youtube.com/watch?v=VKt50_XiORo

Un equipo de mejora investiga el efecto de cuatro métodos de ensamble A, B, C y D, sobre el tiempo de ensamble en minutos. Se va a controlar activamente en el experimento a los operadores que realizarán el ensamble .Los tiempos de ensamble obtenidos se muestran a continuación:


Operadores
Tratamientos 1
A 6
B 7
C 10
D 10

Generalización del experimento y notación usada.

Para el diseño de experimentos de bloques completos al azar tenemos la siguiente tabla general:

TratamientoTratamiento 11 22 33 …… bb TotalesTotales PromedioPromedio
11 y11y11 y12y12 y13y13 …… y1by1b y1.y1. y1.¯y1.¯
22 y21y21 y22y22 y23y23 …… y2by2b y2.y2. y2.¯y2.¯
33 y31y31 y32y32 y33y33 …… y3by3b y3.y3. y3.¯y3.¯
…… …… …… …… …… …… …… ……
aa ya1ya1 ya2ya2 ya3ya3 …… yabyab ya.ya. ya.¯ya.¯
TotalesTotales y.1y.1 y.2y.2 y.3y.3 …… y.by.b y..y.. NA
PromedioPromedio y.1¯y.1¯ y.2¯y.2¯ y.3¯y.3¯ …… y.b¯y.b¯ NA y..¯y..¯

Para la anterior tabla tenemos lo siguiente:

  • aa tratamientos, es decir, i=1,2,…,ai=1,2,…,a

  • bb bloques, es decir, j=1,2,…,bj=1,2,…,b

  • NN cantidad total de observaciones, N=a∗b

Aunque el modelo estadístico se puede expresar de distintas maneras, suele utilizarse el modelo de los efectos:

yij=u+τi+βj+ϵiji=1,2,..,aj=1,2,..,byij=u+τi+βj+ϵiji=1,2,..,aj=1,2,..,b

Donde:

  • μ: media globalμ: media global

  • τi:efecto del i−ésimo tratamientoτi:efecto del i−ésimo tratamiento

  • βj:efecto del j−ésimo bloqueβj:efecto del j−ésimo bloque

  • ϵij: error aleatorioϵij: error aleatorio

# Tratamientos A, B, C, D
# Operadores 1, 2, 3, 4
# Valores (fila A: 6,9,7,8; fila B:7,10,11,8; fila C:10,16,11,14; fila D:10,13,11,9)

# Definimos vectores para Tratamiento, Operador y Respuesta:
tratamiento <- rep(c("A","B","C","D"), each = 4)
operador    <- rep(1:4, times = 4)          # si prefieres nombres: factor(paste0("Op", rep(1:4, times=4)))
respuesta   <- c(
  6, 9, 7, 8,    # Tratamiento A para Operador 1,2,3,4
  7,10,11, 8,    # Tratamiento B
 10,16,11,14,    # Tratamiento C
 10,13,11, 9     # Tratamiento D
)

# Construir data.frame
datos <- data.frame(
  Tratamiento = factor(tratamiento),
  Operador    = factor(operador),
  Respuesta   = respuesta
)

# Vista rápida
print(datos)
##    Tratamiento Operador Respuesta
## 1            A        1         6
## 2            A        2         9
## 3            A        3         7
## 4            A        4         8
## 5            B        1         7
## 6            B        2        10
## 7            B        3        11
## 8            B        4         8
## 9            C        1        10
## 10           C        2        16
## 11           C        3        11
## 12           C        4        14
## 13           D        1        10
## 14           D        2        13
## 15           D        3        11
## 16           D        4         9
# 2. Estadísticas descriptivas básicas
library(dplyr)
## 
## Adjuntando el paquete: 'dplyr'
## The following object is masked from 'package:car':
## 
##     recode
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
resumen <- datos %>%
  group_by(Tratamiento) %>%
  summarise(
    n = n(),
    media = mean(Respuesta),
    sd = sd(Respuesta)
  )
print(resumen)
## # A tibble: 4 × 4
##   Tratamiento     n media    sd
##   <fct>       <int> <dbl> <dbl>
## 1 A               4   7.5  1.29
## 2 B               4   9    1.83
## 3 C               4  12.8  2.75
## 4 D               4  10.8  1.71
# También puedes ver medias por Operador:
medias_bloque <- datos %>%
  group_by(Operador) %>%
  summarise(
    n = n(),
    media = mean(Respuesta),
    sd = sd(Respuesta)
  )
print(medias_bloque)
## # A tibble: 4 × 4
##   Operador     n media    sd
##   <fct>    <int> <dbl> <dbl>
## 1 1            4  8.25  2.06
## 2 2            4 12     3.16
## 3 3            4 10     2   
## 4 4            4  9.75  2.87
library(ggplot2)
# Si no tienes patchwork instalado: install.packages("patchwork")
library(patchwork)

# Creamos los plots y los asignamos a variables
p1 <- ggplot(datos, aes(x = Tratamiento, y = Respuesta)) +
  geom_boxplot(fill = "lightblue") +
  ggtitle("Boxplot de Respuesta por Tratamiento") +
  theme_minimal()

p2 <- ggplot(datos, aes(x = Operador, y = Respuesta)) +
  geom_boxplot(fill = "lightgreen") +
  ggtitle("Boxplot de Respuesta por Operador") +
  theme_minimal()

# Los mostramos juntos: uno al lado del otro
print(p1 + p2)

Podemos observar en las gráficas de caja que el tratamiento A,B son similares al igual que el tratamiento C y D.

El operador 2 fue mas eficiente.

La idea no es hacer un analisis multi-factorial, ya que en teoría, el bloque no influye si estadísticamente influye debe hacer un experimento factorial, ya que no se trata de que el operador sea una variable a controlar sino que realmente es una variable que influye.

#--------------- Análisis de varianza ANOVA  ---------------
#ANOVA RCBD: Respuesta ~ Tratamiento + Operador
mod_anova <- aov(Respuesta ~ Tratamiento + Operador, data = datos)
cat("\n--- Resumen ANOVA (aov) ---\n")
## 
## --- Resumen ANOVA (aov) ---
print(summary(mod_anova))
##             Df Sum Sq Mean Sq F value  Pr(>F)   
## Tratamiento  3   61.5    20.5   10.25 0.00292 **
## Operador     3   28.5     9.5    4.75 0.02985 * 
## Residuals    9   18.0     2.0                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 4. Extraer residuales y parámetros de error para tests posteriores
residuales <- resid(mod_anova)
# Extraer de la tabla ANOVA los grados de libertad y MSerror (Residuals)
anova_tab <- summary(mod_anova)[[1]]
df_resid <- anova_tab["Residuals", "Df"]
ms_resid <- anova_tab["Residuals", "Mean Sq"]
cat("\nGrados de libertad residuales:", df_resid, "\n")
## 
## Grados de libertad residuales: 9
cat("Mean Sq residuales:", ms_resid, "\n")
## Mean Sq residuales: 2

Podemos observar que el operador estadísticamente tiene medias iguales, pero los tratamientos no es así.

# 6. Comparaciones múltiples de medias (post-hoc) si Tratamiento es significativo
# 6.1 Usando emmeans + Tukey
if (!requireNamespace("emmeans", quietly = TRUE)) {
  install.packages("emmeans")
}
library(emmeans)
## Welcome to emmeans.
## Caution: You lose important information if you filter this package's results.
## See '? untidy'
emm <- emmeans(mod_anova, ~ Tratamiento)
contrastes_tukey <- pairs(emm, adjust = "tukey")
cat("\n--- Comparaciones Tukey HSD (Tratamiento) ---\n")
## 
## --- Comparaciones Tukey HSD (Tratamiento) ---
print(contrastes_tukey)
##  contrast estimate SE df t.ratio p.value
##  A - B       -1.50  1  9  -1.500  0.4759
##  A - C       -5.25  1  9  -5.250  0.0024
##  A - D       -3.25  1  9  -3.250  0.0412
##  B - C       -3.75  1  9  -3.750  0.0196
##  B - D       -1.75  1  9  -1.750  0.3548
##  C - D        2.00  1  9   2.000  0.2567
## 
## Results are averaged over the levels of: Operador 
## P value adjustment: tukey method for comparing a family of 4 estimates
plot(emm)

# duncan.test requiere la respuesta, el factor, DFerror y MSerror
duncan_res <- agricolae::duncan.test(
  y = datos$Respuesta,
  trt = datos$Tratamiento,
  DFerror = df_resid,
  MSerror = ms_resid
)
cat("\n--- Resultados Duncan test ---\n")
## 
## --- Resultados Duncan test ---
print(duncan_res)
## $statistics
##   MSerror Df Mean       CV
##         2  9   10 14.14214
## 
## $parameters
##     test            name.t ntr alpha
##   Duncan datos$Tratamiento   4  0.05
## 
## $duncan
##      Table CriticalRange
## 2 3.199173      2.262157
## 3 3.339138      2.361127
## 4 3.419765      2.418139
## 
## $means
##   datos$Respuesta      std r        se Min Max   Q25  Q50   Q75
## A            7.50 1.290994 4 0.7071068   6   9  6.75  7.5  8.25
## B            9.00 1.825742 4 0.7071068   7  11  7.75  9.0 10.25
## C           12.75 2.753785 4 0.7071068  10  16 10.75 12.5 14.50
## D           10.75 1.707825 4 0.7071068   9  13  9.75 10.5 11.50
## 
## $comparison
## NULL
## 
## $groups
##   datos$Respuesta groups
## C           12.75      a
## D           10.75     ab
## B            9.00     bc
## A            7.50      c
## 
## attr(,"class")
## [1] "group"
# Gráfico de agrupaciones Duncan:
plot(duncan_res, main = "Agrupaciones según Duncan test")

Nos queda clao que el método mas eficiente es el el A; también concluimos que usar el método C-D son iguales, al igual que el D-B y el D-B.

Si se tratara de ahorrar y el método B, fuera mas económico de implementar lo haría, tiene solo una diferencia de 2 min; claro esto en producción puede significar dinero.

Validemos supuestos

qqnorm(residuos, main = "QQ-Plot de residuos")
qqline(residuos, col = "red")

shapiro_res <- shapiro.test(residuos)
print(shapiro_res)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuos
## W = 0.95991, p-value = 0.3469
# 5.2 Residuals vs Fitted y Levene
plot(fitted_vals, residuos,
     xlab = "Valores ajustados", ylab = "Residuos",
     main = "Residuals vs Fitted")
abline(h = 0, col = "red")

if (!requireNamespace("car", quietly = TRUE)) {
  install.packages("car")
}
library(car)
levene_res <- car::leveneTest(Tiempo ~ Tratamiento, data = df)
print(levene_res)
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  3  0.1631 0.9201
##       24
par(mfrow = c(1, 1))

2.4 Diseño de cuadrado latino

Un diseño de Cuadrado Latino es un tipo de diseño experimental bloqueado que permite controlar dos fuentes de variabilidad (dos factores de bloqueo) además del tratamiento de interés, utilizando una estructura cuadrada en la que cada tratamiento aparece exactamente una vez en cada fila y en cada columna. A continuación, se presenta un resumen de sus características, usos, ventajas y desventajas, seguido de un ejemplo práctico en R.

2.4.1 1. ¿Qué es un Diseño de Cuadrado Latino?

  • Es un diseño para experimentos con:
    • Un factor de interés (Tratamientos) con \(k\) niveles.
    • Dos factores de bloqueo, conocidos como “filas” y “columnas”, cada uno con \(k\) niveles.
  • Se organiza en una matriz \(k \times k\) tal que:
    • Cada tratamiento aparece exactamente una vez en cada fila.
    • Cada tratamiento aparece exactamente una vez en cada columna.
  • Así, la estructura bloquea simultáneamente la variabilidad asociada a filas y columnas, además de estimar el efecto de los tratamientos.

2.4.2 2. Cuándo usarlo

  • Hay dos fuentes conocidas de variabilidad que conviene controlar y cada una tiene \(k\) niveles (por ejemplo, posición en filas y columnas en un campo agrícola, día y operario en laboratorio, lote y máquina, etc.).
  • El número de tratamientos debe igualar el número de niveles de cada factor de bloqueo (\(k\)).
  • Se desea mantener el número de réplicas relativamente bajo: solo una observación por celda del cuadrado (total \(k^2\) unidades experimentales).

2.4.3 3. Ventajas

  • Controla dos fuentes de variabilidad sin necesidad de un número excesivo de réplicas.
  • Más eficiente que un diseño completamente aleatorizado o bloqueado simple cuando hay dos factores de bloqueo con muchas variaciones.
  • Reduce la varianza residual al ajustar la variabilidad debida a filas y columnas.

2.4.4 4. Desventajas y limitaciones

  • El número de tratamientos \(k\) debe coincidir con el número de niveles de filas y columnas; esto impone rigidez (no se puede tener diferente número de tratamientos).
  • No se pueden estimar interacciones entre tratamientos y bloqueos (fila o columna): se asume que no hay interacción o es despreciable.
  • Si alguno de los factores de bloqueo no es homogéneo o faltan observaciones, se complica el diseño o su análisis.
  • Menos flexible que diseños factoriales o de bloques incompletos; requiere que el experimento se ajuste exactamente a la estructura \(k \times k\).

2.4.5 5. Supuestos

  • Independencia: Las unidades experimentales (cada celda) son independientes.
  • Efectos aditivos: El modelo asume que los efectos de fila, columna y tratamiento se suman sin interacción significativa.
  • Normalidad y homogeneidad de varianzas: Para el análisis ANOVA, se asume que los residuos son aproximadamente normales y tienen varianza constante en todas las celdas.
  • Correcta aleatorización: La asignación de tratamientos a la matriz debe respetar la regla del cuadrado latino y la asignación de la disposición de filas/columnas debe ser aleatoria o balanceada según el contexto.

2.4.6 6. Estructura del modelo ANOVA

El modelo general es: \[ Y_{ij} = \mu + \alpha_i + \beta_j + \tau_{t(i,j)} + \epsilon_{ij}, \] donde: - \(i\) = índice de fila (1,…,k), - \(j\) = índice de columna (1,…,k), - \(\alpha_i\) = efecto de la i-ésima fila (bloque fila), - \(\beta_j\) = efecto de la j-ésima columna (bloque columna), - \(\tau_{t(i,j)}\) = efecto del tratamiento asignado a la celda \((i,j)\), - \(\epsilon_{ij}\) = error aleatorio ~ \(N(0, \sigma^2)\).

Se ajusta un ANOVA con tres factores: fila, columna y tratamiento, y se evalúa la significancia de \(\tau\).

2.4.7 7. Ejemplo clásico de aplicación

  • Agricultura: Parcela dividida en una cuadrícula \(k \times k\), donde filas y columnas representan dos fuentes de variabilidad del suelo (p.ej., pendiente y tipo de suelo). Los \(k\) tratamientos (variedades de cultivo o dosis) se asignan según un cuadrado latino.
  • Laboratorio: Experimento en microplacas \(k \times k\), donde filas pueden ser fecha de experimento y columnas operario o dispositivo; los tratamientos se asignan en la placa de modo latino.
  • Industria: Pruebas en una matriz de máquinas y operarios, asignando tratamientos de forma que cada tratamiento se pruebe una vez con cada operario y con cada máquina.

2.4.8 Diseño de Cuadrado Latino: Ejemplo con comparación de cuatro marcas de llantas

A continuación se presenta un ejemplo completo en R de análisis de un experimento en cuadrado latino para comparar cuatro marcas de llantas. Se incluyen:

  1. Planteamiento y datos observados (tomados del enunciado).
  2. Construcción del data.frame en R.
  3. Ajuste de ANOVA para diseño de cuadrado latino.
  4. Verificación de supuestos.
  5. Cálculo de tamaños de efecto.
  6. Interpretación de resultados.
  7. (Opcional) Comparaciones post-hoc y visualizaciones.

2.4.8.1 1. Enunciado resumido

Una compañía de mensajería quiere determinar cuál marca de llantas tiene mayor duración en términos de desgaste. Se planifica un experimento en cuadrado latino de orden 4:

  • Tratamientos (factor de interés): 4 marcas de llantas, denotadas A, B, C y D.
  • Bloque 1: Tipo de coche (4 tipos distintos), denotados Carro 1, Carro 2, Carro 3, Carro 4.
  • Bloque 2: Posición de la llanta en el coche (4 posiciones, por ejemplo: delanteras/traseras o zonas, aquí denotadas Posición 1, 2, 3 y 4).
  • Se asume que cada “celda” del cuadrado latino (combinación Carro i, Posición j) recibe exactamente una marca de llanta, asignada de forma tal que cada marca aparece una vez en cada fila (Posición) y una vez en cada columna (Carro).
  • Se realiza la prueba sometiendo cada llanta a una prueba de 32 000 km; se mide alguna cuantificación de desgaste o duración, por ejemplo, un puntaje donde mayor valor indica mayor duración (en el enunciado de ejemplo se muestran valores numéricos).
  • La tabla de observaciones (Tabla 4.6 del enunciado) es:
Posición  Carro Carro 1 Carro 2 Carro 3 Carro 4
Posición 1 C = 12 D = 11 A = 13 B = 8
Posición 2 B = 14 C = 12 D = 11 A = 13
Posición 3 A = 17 B = 14 C = 10 D = 9
Posición 4 D = 13 A = 14 B = 13 C = 9
  • Aquí, el número tras “=” es el resultado de la prueba de duración/desgaste (a mayor valor, mayor duración).
  • Objetivo: Analizar si hay diferencias significativas entre las cuatro marcas de llantas, controlando la variabilidad debida al tipo de coche y la posición de la llanta mediante el diseño en cuadrado latino.

2.4.8.2 2. Construcción del dataset en R

Creamos un data.frame con 16 observaciones y columnas: - Carro: factor con niveles "Carro1", "Carro2", "Carro3", "Carro4". - Posicion: factor con niveles "Pos1", "Pos2", "Pos3", "Pos4". - Marca: factor con niveles "A", "B", "C", "D". - Duracion: numérico con los valores observados.

df <- data.frame(
  Respuesta = c(
    # Carro 1, Posición 1 a 4
    12, 14, 17, 13,
    # Carro 2, Posición 1 a 4
    11, 12, 14, 14,
    # Carro 3, Posición 1 a 4
    13, 11, 10, 13,
    # Carro 4, Posición 1 a 4
    8,  13,  9,  9
  ),
  Carro = factor(rep(paste0("Carca", 1:4), each = 4),
                 levels = paste0("Carca", 1:4)),
  Posicion = factor(rep(paste0("Pos", 1:4), times = 4),
                     levels = paste0("Pos", 1:4)),
  Latina = factor(c(
    # Carro1: C, B, A, D
    "C","B","A","D",
    # Carro2: D, C, B, A
    "D","C","B","A",
    # Carro3: A, D, C, B
    "A","D","C","B",
    # Carro4: B, A, D, C
    "B","A","D","C"
  ), levels = c("A","B","C","D"))
)

# Inspección rápida
str(df)
## 'data.frame':    16 obs. of  4 variables:
##  $ Respuesta: num  12 14 17 13 11 12 14 14 13 11 ...
##  $ Carro    : Factor w/ 4 levels "Carca1","Carca2",..: 1 1 1 1 2 2 2 2 3 3 ...
##  $ Posicion : Factor w/ 4 levels "Pos1","Pos2",..: 1 2 3 4 1 2 3 4 1 2 ...
##  $ Latina   : Factor w/ 4 levels "A","B","C","D": 3 2 1 4 4 3 2 1 1 4 ...
print(df)
##    Respuesta  Carro Posicion Latina
## 1         12 Carca1     Pos1      C
## 2         14 Carca1     Pos2      B
## 3         17 Carca1     Pos3      A
## 4         13 Carca1     Pos4      D
## 5         11 Carca2     Pos1      D
## 6         12 Carca2     Pos2      C
## 7         14 Carca2     Pos3      B
## 8         14 Carca2     Pos4      A
## 9         13 Carca3     Pos1      A
## 10        11 Carca3     Pos2      D
## 11        10 Carca3     Pos3      C
## 12        13 Carca3     Pos4      B
## 13         8 Carca4     Pos1      B
## 14        13 Carca4     Pos2      A
## 15         9 Carca4     Pos3      D
## 16         9 Carca4     Pos4      C
# Usando aggregate
res_latina <- aggregate(Respuesta ~ Latina, data = df,
                       FUN = function(x) c(media = mean(x), sd = sd(x), n = length(x)))
# Convertir resultado a data.frame plano
res_latina2 <- data.frame(
  Latina = res_latina$Latina,
  media  = res_latina$Respuesta[, "media"],
  sd     = res_latina$Respuesta[, "sd"],
  n      = res_latina$Respuesta[, "n"]
)
print(res_latina2)
##   Latina media       sd n
## 1      A 14.25 1.892969 4
## 2      B 12.25 2.872281 4
## 3      C 10.75 1.500000 4
## 4      D 11.00 1.632993 4

Es importante entender que las letras latinas A,B,C,D se usan para las diferentes marcar de llantas

mod_anova <- aov(Respuesta ~ Carro + Posicion+Latina, data = df)
cat("\n--- Resumen ANOVA (aov) ---\n")
## 
## --- Resumen ANOVA (aov) ---
print(summary(mod_anova))
##             Df Sum Sq Mean Sq F value  Pr(>F)   
## Carro        3  38.69  12.896  14.395 0.00378 **
## Posicion     3   6.19   2.062   2.302 0.17695   
## Latina       3  30.69  10.229  11.419 0.00683 **
## Residuals    6   5.37   0.896                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 4. Extraer residuales y parámetros de error para tests posteriores
residuales <- resid(mod_anova)
# Extraer de la tabla ANOVA los grados de libertad y MSerror (Residuals)
anova_tab <- summary(mod_anova)[[1]]
df_resid <- anova_tab["Residuals", "Df"]
ms_resid <- anova_tab["Residuals", "Mean Sq"]
cat("\nGrados de libertad residuales:", df_resid, "\n")
## 
## Grados de libertad residuales: 6
cat("Mean Sq residuales:", ms_resid, "\n")
## Mean Sq residuales: 0.8958333

Como podemos observar el factor Latino que en realidad es la marca de llanta, es significativamente diferente y también el carro. según nuestra anova la posición de la llanta es poco significante, o da igual la posición de la llanta pero me absorbe varianza.

# 3. Post-hoc (Tukey)

# 6. Comparaciones múltiples de medias (post-hoc) si Tratamiento es significativo
# 6.1 Usando emmeans + Tukey
if (!requireNamespace("emmeans", quietly = TRUE)) {
  install.packages("emmeans")
}
library(emmeans)

TukeyHSD(mod_anova, "Latina")
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Respuesta ~ Carro + Posicion + Latina, data = df)
## 
## $Latina
##      diff       lwr        upr     p adj
## B-A -2.00 -4.316805  0.3168049 0.0872429
## C-A -3.50 -5.816805 -1.1831951 0.0078229
## D-A -3.25 -5.566805 -0.9331951 0.0112170
## C-B -1.50 -3.816805  0.8168049 0.2144635
## D-B -1.25 -3.566805  1.0668049 0.3315951
## D-C  0.25 -2.066805  2.5668049 0.9805997
library(agricolae)
emm <- emmeans(mod_anova, ~ Latina)
contrastes_tukey <- pairs(emm, adjust = "tukey")
cat("\n--- Comparaciones Tukey HSD (Tratamiento) ---\n")
## 
## --- Comparaciones Tukey HSD (Tratamiento) ---
print(contrastes_tukey)
##  contrast estimate    SE df t.ratio p.value
##  A - B        2.00 0.669  6   2.988  0.0872
##  A - C        3.50 0.669  6   5.230  0.0078
##  A - D        3.25 0.669  6   4.856  0.0112
##  B - C        1.50 0.669  6   2.241  0.2145
##  B - D        1.25 0.669  6   1.868  0.3316
##  C - D       -0.25 0.669  6  -0.374  0.9806
## 
## Results are averaged over the levels of: Carro, Posicion 
## P value adjustment: tukey method for comparing a family of 4 estimates
plot(contrastes_tukey)

# duncan.test requiere la respuesta, el factor, DFerror y MSerror
resultado_lsd <- LSD.test(mod_anova, "Latina", p.adj = "none")  # sin ajuste, típico en LSD
print(resultado_lsd)
## $statistics
##     MSerror Df    Mean       CV  t.value      LSD
##   0.8958333  6 12.0625 7.846505 2.446912 1.637634
## 
## $parameters
##         test p.ajusted name.t ntr alpha
##   Fisher-LSD      none Latina   4  0.05
## 
## $means
##   Respuesta      std r        se       LCL      UCL Min Max   Q25  Q50   Q75
## A     14.25 1.892969 4 0.4732424 13.092018 15.40798  13  17 13.00 13.5 14.75
## B     12.25 2.872281 4 0.4732424 11.092018 13.40798   8  14 11.75 13.5 14.00
## C     10.75 1.500000 4 0.4732424  9.592018 11.90798   9  12  9.75 11.0 12.00
## D     11.00 1.632993 4 0.4732424  9.842018 12.15798   9  13 10.50 11.0 11.50
## 
## $comparison
## NULL
## 
## $groups
##   Respuesta groups
## A     14.25      a
## B     12.25      b
## D     11.00      b
## C     10.75      b
## 
## attr(,"class")
## [1] "group"
plot(resultado_lsd)

Podemos concluir que algunas marcas de llantas no significativamente iguales B - C , C - D ,B - D, en la prueba de tukey; en la prueba LSD.test me da que solo el A, es diferente.

Podemos concluir que debería usarse la Marca A, debido que es el que maximiza el recorrido de kilómetros..

Tener muy en cuenta que el cuadrado latino de este ejemplo es solo un ejemplo se puede diseñar la matrix de cuadrado latino casi como se nos ocurra ejemplo

     [,1] [,2] [,3] [,4]
[1,] "A"  "B"  "C"  "D" 
[2,] "D"  "A"  "B"  "C" 
[3,] "C"  "D"  "A"  "B" 
[4,] "B"  "C"  "D"  "A" 
design.lsd(
  trt   = LETTERS[1:4], 
  r     = 4, 
   serie = 1)
## $parameters
## $parameters$design
## [1] "lsd"
## 
## $parameters$trt
## [1] "A" "B" "C" "D"
## 
## $parameters$r
## [1] 4
## 
## $parameters$serie
## [1] 1
## 
## $parameters$seed
## [1] 1643942451
## 
## $parameters$kinds
## [1] "Super-Duper"
## 
## $parameters[[7]]
## [1] 4
## 
## 
## $sketch
##      [,1] [,2] [,3] [,4]
## [1,] "C"  "D"  "A"  "B" 
## [2,] "A"  "B"  "C"  "D" 
## [3,] "B"  "C"  "D"  "A" 
## [4,] "D"  "A"  "B"  "C" 
## 
## $book
##    plots row col LETTERS[1:4]
## 1     11   1   1            C
## 2     12   1   2            D
## 3     13   1   3            A
## 4     14   1   4            B
## 5     21   2   1            A
## 6     22   2   2            B
## 7     23   2   3            C
## 8     24   2   4            D
## 9     31   3   1            B
## 10    32   3   2            C
## 11    33   3   3            D
## 12    34   3   4            A
## 13    41   4   1            D
## 14    42   4   2            A
## 15    43   4   3            B
## 16    44   4   4            C

Con este código podemos generar aleatoriamente cuadrados latinos.

Condiciones de un cudrado latino

Sí, ese cuadrado latino es válido. Un cuadrado latino debe cumplir dos condiciones:

  1. Cada símbolo aparece exactamente una vez en cada fila

  2. Cada símbolo aparece exactamente una vez en cada columna

2.5 Diseño cuadrado grecolatino

https://youtu.be/9AwREFvp6aY?si=iX6H7aETaYPIkN5p

https://rpubs.com/Gabo6381/768438

Con el Diseño en Cuadro Grecolatino (DCGL) se controlan tres factores de bloque, además del factor de tratamiento. Se llama cuadro grecolatino porque los cuatro factores involucrados se prueban en la misma cantidad de niveles, de aquí que se pueda escribir como un cuadro, como se muestra en la Tabla 1. Se utilizan letras para denotar los tratamientos, y letras griegas para nombrar los niveles del tercer factor de bloque, mientras que el primer y segundo factor de bloque se ubican en los renglones y en las columnas (Pulido & Vara Salazar, 2012):

Tabla 1. Estructura del Diseño en Cuadro Grecolatino
1 2 3 4
Aαα Bββ Cγγ Dδδ
Bδδ Aγγ Cββ Dαα
Cββ Dαα Aδδ Bγγ
Dγγ Cδδ Bαα Aβ

Tener encenta que no debe cumplir este orden de tabla y podemos generarla como deseamos con el siguiente código.

design.graeco(
    trt1 = LETTERS[1:4],
    trt2 = c("α","β","γ","δ"),
    serie = 1
)
## $parameters
## $parameters$design
## [1] "graeco"
## 
## $parameters$trt1
## [1] "A" "B" "C" "D"
## 
## $parameters$trt2
## [1] "α" "β" "γ" "δ"
## 
## $parameters$r
## [1] 4
## 
## $parameters$serie
## [1] 1
## 
## $parameters$seed
## [1] -2067387977
## 
## $parameters$kinds
## [1] "Super-Duper"
## 
## $parameters[[8]]
## [1] TRUE
## 
## 
## $sketch
##      [,1]  [,2]  [,3]  [,4] 
## [1,] "D δ" "A γ" "C α" "B β"
## [2,] "A α" "D β" "B δ" "C γ"
## [3,] "C β" "B α" "D γ" "A δ"
## [4,] "B γ" "C δ" "A β" "D α"
## 
## $book
##    plots row col LETTERS[1:4] c("α", "β", "γ", "δ")
## 1     11   1   1            D                     δ
## 2     12   1   2            A                     γ
## 3     13   1   3            C                     α
## 4     14   1   4            B                     β
## 5     21   2   1            A                     α
## 6     22   2   2            D                     β
## 7     23   2   3            B                     δ
## 8     24   2   4            C                     γ
## 9     31   3   1            C                     β
## 10    32   3   2            B                     α
## 11    33   3   3            D                     γ
## 12    34   3   4            A                     δ
## 13    41   4   1            B                     γ
## 14    42   4   2            C                     δ
## 15    43   4   3            A                     β
## 16    44   4   4            D                     α

El modelo estadístico que describe las mediciones en un cuadro grecolatino está dado por:

yijk=μ+τi+βj+δk+φl+εijklyijk=μ+τi+βj+δk+φl+εijkl

Donde yijklyijkl es la observación o respuesta que se encuentra en el tratamiento i (i-ésima letra latina), en el renglón j, en la columna k y en la l-ésima letra griega, que son los niveles del tercer factor de bloque; el término εijklεijkl representa el error aleatoio atribuible a la medición yijklyijkl.

##Caso Práctico

Para ejemplificar la metodología del DCGL tomaremos como base el ejercicio 25 de la página 111 del libro de Análsis y Diseño de Experimentos (Pulido & Vara Salazar, 2012), en el que se considera que un investigador está interesado en el efecto del porcentaje de lisina y del porcentaje de de proteína en la producción de vacas lecheras. Se consideran siete niveles de cada factor:

  • Porcentaje de lisina: 0.0(A),0.1(B),0.2(C),0.3(D),0.4(E),0.5(F),0.6(G).

  • Porcentaje de proteína: 2(αα),4(ββ),6(χχ),8(δδ),10(εε),12(φφ),14(γγ).

Para el estudio, se seleccionan siete vacas al azar, a las cuales se les da un seguimiento de siete periodos de tres meses. Los datos en galones de leche fueron los siguientes:

Tabla 2. Resultados del experimento.
Vaca/Periodo 1 2 3 4 5 6 7
1 304(Aαα) 436(Bεε) 350(Cββ) 504(Dφφ) 417(Eχχ) 519(Fγγ) 432(Gδδ)
2 381(Bββ) 505(Cφφ) 425(Dχχ) 564(Eγγ) 494(Fδδ) 350(Gαα) 413(GAεε)
3 432(Cχχ) 566(Dγγ) 479(Eδδ) 357(Fαα) 461(Gεε) 340(Aββ) 502(Bφφ)
4 442(Dδδ) 372(Eαα) 536(Fεε) 366(Gββ) 495(Aφφ) 425(Bχχ) 507(Cγγ)
5 496(Eεε) 449(Fββ) 493(Gφφ) 345(Aχχ) 509(Bγγ) 481(Cδδ) 380(Dαα)
6 534(Fφφ) 421(Gχχ) 352(Aγγ) 427(Bδδ) 346(Cαα) 478(Dεε) 397(Eββ)
7 543(Gγγ) 386(Aδδ) 435(Bαα) 485(Cεε) 406(Dββ) 554(Eφφ) 410(Fχχ)

  1. Analice este experimento. ¿Qué factores tienen efecto en la producción de leche?

  2. Interprete los resultados usando gráfico de medias.

  3. ¿Cómo puede explicar la falta de efectos en vacas y periodo?

  4. ¿Qué porcentajes de lisina y proteína dan los mejores resultados?

  5. Verifique los supuestos del modelo.

2.5.1 Analisis del experimento

# Data frame completo en el orden: Galones, Lisina, Vaca, Periodo, Proteina
df <- data.frame(
  Galones = c(
    # Periodo 1 (Vacas 1–7)
    304, 381, 432, 442, 496, 534, 543,
    # Periodo 2
    436, 505, 566, 372, 449, 421, 386,
    # Periodo 3
    350, 425, 479, 536, 493, 352, 435,
    # Periodo 4
    504, 564, 357, 366, 345, 427, 485,
    # Periodo 5
    417, 494, 461, 495, 509, 346, 406,
    # Periodo 6
    519, 350, 340, 425, 481, 478, 554,
    # Periodo 7
    432, 413, 502, 507, 380, 397, 410
  ),
  Lisina = c(
    # Periodo 1
    "A","B","C","D","E","F","G",
    # Periodo 2
    "B","C","D","E","F","G","A",
    # Periodo 3
    "C","D","E","F","G","A","B",
    # Periodo 4
    "D","E","F","G","A","B","C",
    # Periodo 5
    "E","F","G","A","B","C","D",
    # Periodo 6
    "F","G","A","B","C","D","E",
    # Periodo 7
    "G","A","B","C","D","E","F"
  ),
  Vaca =as.factor(rep(1:7, times = 7)),
  Periodo = as.factor(rep(1:7, each = 7)),
  Proteina = c(
    # Periodo 1
    "α","β","χ","δ","ε","φ","γ",
    # Periodo 2
    "ε","φ","γ","α","β","χ","δ",
    # Periodo 3
    "β","χ","δ","ε","φ","γ","α",
    # Periodo 4
    "φ","γ","α","β","χ","δ","ε",
    # Periodo 5
    "χ","δ","ε","φ","γ","α","β",
    # Periodo 6
    "γ","α","β","χ","δ","ε","φ",
    # Periodo 7
    "δ","ε","φ","γ","α","β","χ"
  ),
  stringsAsFactors = FALSE
)

# Mostrar el resultado
print(df)
##    Galones Lisina Vaca Periodo Proteina
## 1      304      A    1       1        α
## 2      381      B    2       1        β
## 3      432      C    3       1        χ
## 4      442      D    4       1        δ
## 5      496      E    5       1        ε
## 6      534      F    6       1        φ
## 7      543      G    7       1        γ
## 8      436      B    1       2        ε
## 9      505      C    2       2        φ
## 10     566      D    3       2        γ
## 11     372      E    4       2        α
## 12     449      F    5       2        β
## 13     421      G    6       2        χ
## 14     386      A    7       2        δ
## 15     350      C    1       3        β
## 16     425      D    2       3        χ
## 17     479      E    3       3        δ
## 18     536      F    4       3        ε
## 19     493      G    5       3        φ
## 20     352      A    6       3        γ
## 21     435      B    7       3        α
## 22     504      D    1       4        φ
## 23     564      E    2       4        γ
## 24     357      F    3       4        α
## 25     366      G    4       4        β
## 26     345      A    5       4        χ
## 27     427      B    6       4        δ
## 28     485      C    7       4        ε
## 29     417      E    1       5        χ
## 30     494      F    2       5        δ
## 31     461      G    3       5        ε
## 32     495      A    4       5        φ
## 33     509      B    5       5        γ
## 34     346      C    6       5        α
## 35     406      D    7       5        β
## 36     519      F    1       6        γ
## 37     350      G    2       6        α
## 38     340      A    3       6        β
## 39     425      B    4       6        χ
## 40     481      C    5       6        δ
## 41     478      D    6       6        ε
## 42     554      E    7       6        φ
## 43     432      G    1       7        δ
## 44     413      A    2       7        ε
## 45     502      B    3       7        φ
## 46     507      C    4       7        γ
## 47     380      D    5       7        α
## 48     397      E    6       7        β
## 49     410      F    7       7        χ
modelo=lm(Galones~(Lisina+Vaca+Periodo+Proteina),data=df)

summary(modelo)
## 
## Call:
## lm(formula = Galones ~ (Lisina + Vaca + Periodo + Proteina), 
##     data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -65.082 -13.796  -0.082  11.204  56.776 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 281.7959    22.4284  12.564 4.80e-12 ***
## LisinaB      68.5714    16.7839   4.086 0.000424 ***
## LisinaC      67.2857    16.7839   4.009 0.000515 ***
## LisinaD      80.8571    16.7839   4.818 6.60e-05 ***
## LisinaE      92.0000    16.7839   5.481 1.23e-05 ***
## LisinaF      94.8571    16.7839   5.652 8.07e-06 ***
## LisinaG      61.5714    16.7839   3.668 0.001212 ** 
## Vaca2        24.2857    16.7839   1.447 0.160843    
## Vaca3        25.0000    16.7839   1.490 0.149375    
## Vaca4        25.8571    16.7839   1.541 0.136499    
## Vaca5        27.2857    16.7839   1.626 0.117073    
## Vaca6        -1.0000    16.7839  -0.060 0.952983    
## Vaca7        36.7143    16.7839   2.187 0.038684 *  
## Periodo2      0.4286    16.7839   0.026 0.979840    
## Periodo3     -8.8571    16.7839  -0.528 0.602542    
## Periodo4    -12.0000    16.7839  -0.715 0.481525    
## Periodo5     -0.5714    16.7839  -0.034 0.973122    
## Periodo6      2.1429    16.7839   0.128 0.899471    
## Periodo7    -13.0000    16.7839  -0.775 0.446169    
## Proteinaβ    20.7143    16.7839   1.234 0.229087    
## Proteinaγ   145.1429    16.7839   8.648 7.74e-09 ***
## Proteinaδ    85.2857    16.7839   5.081 3.38e-05 ***
## Proteinaε   108.7143    16.7839   6.477 1.07e-06 ***
## Proteinaφ   149.0000    16.7839   8.878 4.77e-09 ***
## Proteinaχ    47.2857    16.7839   2.817 0.009537 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 31.4 on 24 degrees of freedom
## Multiple R-squared:  0.8938, Adjusted R-squared:  0.7876 
## F-statistic: 8.417 on 24 and 24 DF,  p-value: 8.969e-07

Podemos observar que nos da un R2Ajustado de 78.76% que explica nuestro modelo, podemos observar que tanto la variable vaca como variable periodo en algunos niveles, superan nuestro p_value de 0.05, indicándonos que la pendiente es 0, y que no no aporta a explicar el modelo, vamos hacer una prueba simple nuestro modelo lineal sin las variables proteína o vaca, para ver que tanto lo explica, recuerde que estas variables de bloqueo, la intención es que aportan variabilidad, de tal modo que me explique el modelo y que los residuos no absorban toda la variabilidad restante.

summary(lm(Galones~(Lisina+Proteina),data=df))
## 
## Call:
## lm(formula = Galones ~ (Lisina + Proteina), data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -90.122 -14.551  -1.694  15.163  69.449 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   296.98      15.87  18.712  < 2e-16 ***
## LisinaB        68.57      16.47   4.163 0.000187 ***
## LisinaC        67.29      16.47   4.085 0.000235 ***
## LisinaD        80.86      16.47   4.909 1.98e-05 ***
## LisinaE        92.00      16.47   5.586 2.49e-06 ***
## LisinaF        94.86      16.47   5.759 1.46e-06 ***
## LisinaG        61.57      16.47   3.738 0.000642 ***
## Proteinaβ      20.71      16.47   1.258 0.216594    
## Proteinaγ     145.14      16.47   8.813 1.62e-10 ***
## Proteinaδ      85.29      16.47   5.178 8.70e-06 ***
## Proteinaε     108.71      16.47   6.601 1.10e-07 ***
## Proteinaφ     149.00      16.47   9.047 8.41e-11 ***
## Proteinaχ      47.29      16.47   2.871 0.006812 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 30.81 on 36 degrees of freedom
## Multiple R-squared:  0.8466, Adjusted R-squared:  0.7955 
## F-statistic: 16.56 on 12 and 36 DF,  p-value: 3.388e-11

Como podemos observar proteína y lisina me explican el 79.55, antes me mejoro, si fuera un modelo predictivo quitaría vaca y periodo, pero en el diseño de experimento vamos a dejarla para que que absorba variabilidad.

boxplot(Galones~ Periodo, xlab="Identificador de Periodo", ylab = "Galones de leche",data=df)

boxplot(Galones~Vaca, col="skyblue", xlab = "Identificador de Vacas", ylab="Galones de leche",data=df)

Podemos observar que las variables de bloqueo, tienden a tener la misma media eso explica un modo el resultado del modelo lineal.

boxplot(Galones~Lisina, col="red", xlab= "Porcentaje de Lisina", ylab = "Galones de leche",data=df)

boxplot(Galones~Proteina, col="orange", xlab="Porcentaje de Proteína", ylab="Galones de leche",df)

Tanto la proteína y la lisina, que son variables de tratamiento, es notorio que tienen medias diferentes y que la proteína que maximiza mas la producción de leche es γ que es la proteína 12.

La lisina que maximiza mas la producción de leche es la E que es lisina 0.4

anova=aov(modelo)
summary(anova)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## Lisina       6  42784    7131   7.232 0.000171 ***
## Vaca         6   8754    1459   1.480 0.227194    
## Periodo      6   1761     293   0.298 0.931964    
## Proteina     6 145880   24313  24.660 3.77e-09 ***
## Residuals   24  23663     986                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
MSerror <-deviance(anova)/df.residual(anova)
DFerror <- df.residual(anova)

Podemos observar que tanto la lisina como la proteína tienen medias diferentes.

El periodo podemos concluir que definitivamente no influye.

Vamos a estudiar en la proteína y lisina que tiene cambios

library(agricolae)
rm_lisina=duncan.test(y=df$Galones,trt=df$Lisina,MSerror = MSerror, DFerror = DFerror, alpha = 0.05,group = TRUE)
print(rm_lisina)
## $statistics
##   MSerror Df     Mean       CV
##   985.949 24 442.8776 7.089956
## 
## $parameters
##     test    name.t ntr alpha
##   Duncan df$Lisina   7  0.05
## 
## $duncan
##      Table CriticalRange
## 2 2.918793      34.64029
## 3 3.065610      36.38272
## 4 3.159874      37.50144
## 5 3.226454      38.29162
## 6 3.276155      38.88146
## 7 3.314602      39.33776
## 
## $means
##   df$Galones      std r       se Min Max   Q25 Q50   Q75
## A   376.4286 62.77701 7 11.86802 304 495 342.5 352 399.5
## B   445.0000 45.36151 7 11.86802 381 509 426.0 435 469.0
## C   443.7143 69.90878 7 11.86802 346 507 391.0 481 495.0
## D   457.2857 63.65196 7 11.86802 380 566 415.5 442 491.0
## E   468.4286 75.68984 7 11.86802 372 564 407.0 479 525.0
## F   471.2857 68.58988 7 11.86802 357 536 429.5 494 526.5
## G   438.0000 68.10776 7 11.86802 350 543 393.5 432 477.0
## 
## $comparison
## NULL
## 
## $groups
##   df$Galones groups
## F   471.2857      a
## E   468.4286      a
## D   457.2857      a
## B   445.0000      a
## C   443.7143      a
## G   438.0000      a
## A   376.4286      b
## 
## attr(,"class")
## [1] "group"
plot(rm_lisina)

Podemos observar que de la lisina todos los niveles tenían medias iguales menos la A que es al 0.0

Los resultados de la prueba de rango múltiple arrojan como resultado que, en lo general, todos los niveles de lisina generan en promedio la misma producción de leche, salvo en el caso del nivel A, que corresponde al nivel en el que no se suministra lisina a las unidades experimentales.

library(agricolae)
rm_lisina=duncan.test(y=df$Galones,trt=df$Proteina,MSerror = MSerror, DFerror = DFerror, alpha = 0.05,group = TRUE)
print(rm_lisina)
## $statistics
##   MSerror Df     Mean       CV
##   985.949 24 442.8776 7.089956
## 
## $parameters
##     test      name.t ntr alpha
##   Duncan df$Proteina   7  0.05
## 
## $duncan
##      Table CriticalRange
## 2 2.918793      34.64029
## 3 3.065610      36.38272
## 4 3.159874      37.50144
## 5 3.226454      38.29162
## 6 3.276155      38.88146
## 7 3.314602      39.33776
## 
## $means
##   df$Galones      std r       se Min Max   Q25 Q50   Q75
## α   363.4286 39.84912 7 11.86802 304 435 348.0 357 376.0
## β   384.1429 37.19959 7 11.86802 340 449 358.0 381 401.5
## γ   508.5714 73.23673 7 11.86802 352 566 508.0 519 553.5
## δ   448.7143 38.16506 7 11.86802 386 494 429.5 442 480.0
## ε   472.1429 40.36264 7 11.86802 413 536 448.5 478 490.5
## φ   512.4286 22.76589 7 11.86802 493 554 498.5 504 519.5
## χ   410.7143 29.79214 7 11.86802 345 432 413.5 421 425.0
## 
## $comparison
## NULL
## 
## $groups
##   df$Galones groups
## φ   512.4286      a
## γ   508.5714      a
## ε   472.1429      b
## δ   448.7143      b
## χ   410.7143      c
## β   384.1429     cd
## α   363.4286      d
## 
## attr(,"class")
## [1] "group"
plot(rm_lisina)

Los resultados de la prueba de rango múltiple arrojan como resultado que, los niveles del factor de bloque denominado proteína que generan una producción promedio más alta son los niveles φφ y γγ, en orden descendente, los niveles que seguirían son εε y δδ, con la misma producción promedio entre ellos, para el caso del nivel χχ, pertenece a dos grupos, los cuales no están claramente diferenciados junto con el factor ββ, el misma factor ββ se incluye junto con el factor αα en el último grupo con la menor producción de leche.

Si queremos validar la mejor combinación de los factores vemos el siguiente ejemplo

https://claude.ai/share/00cde676-ab9c-43a3-a9d2-da0b5f02da6c

library(emmeans) 
combinaciones <- emmeans(modelo, ~ Lisina * Proteina) 
print(combinaciones)
##  Lisina Proteina emmean   SE df lower.CL upper.CL
##  A      α           297 16.2 24      264      330
##  B      α           366 16.2 24      332      399
##  C      α           364 16.2 24      331      398
##  D      α           378 16.2 24      344      411
##  E      α           389 16.2 24      356      422
##  F      α           392 16.2 24      358      425
##  G      α           359 16.2 24      325      392
##  A      β           318 16.2 24      284      351
##  B      β           386 16.2 24      353      420
##  C      β           385 16.2 24      352      418
##  D      β           399 16.2 24      365      432
##  E      β           410 16.2 24      376      443
##  F      β           413 16.2 24      379      446
##  G      β           379 16.2 24      346      413
##  A      γ           442 16.2 24      409      476
##  B      γ           511 16.2 24      477      544
##  C      γ           509 16.2 24      476      543
##  D      γ           523 16.2 24      490      556
##  E      γ           534 16.2 24      501      568
##  F      γ           537 16.2 24      504      570
##  G      γ           504 16.2 24      470      537
##  A      δ           382 16.2 24      349      416
##  B      δ           451 16.2 24      417      484
##  C      δ           450 16.2 24      416      483
##  D      δ           463 16.2 24      430      497
##  E      δ           474 16.2 24      441      508
##  F      δ           477 16.2 24      444      511
##  G      δ           444 16.2 24      410      477
##  A      ε           406 16.2 24      372      439
##  B      ε           474 16.2 24      441      508
##  C      ε           473 16.2 24      440      506
##  D      ε           487 16.2 24      453      520
##  E      ε           498 16.2 24      464      531
##  F      ε           501 16.2 24      467      534
##  G      ε           467 16.2 24      434      501
##  A      φ           446 16.2 24      413      479
##  B      φ           515 16.2 24      481      548
##  C      φ           513 16.2 24      480      547
##  D      φ           527 16.2 24      493      560
##  E      φ           538 16.2 24      505      571
##  F      φ           541 16.2 24      507      574
##  G      φ           508 16.2 24      474      541
##  A      χ           344 16.2 24      311      378
##  B      χ           413 16.2 24      379      446
##  C      χ           412 16.2 24      378      445
##  D      χ           425 16.2 24      392      459
##  E      χ           436 16.2 24      403      470
##  F      χ           439 16.2 24      406      473
##  G      χ           406 16.2 24      372      439
## 
## Results are averaged over the levels of: Vaca, Periodo 
## Confidence level used: 0.95
mejor <- summary(combinaciones)[which.max(summary(combinaciones)$emmean), ]
print(mejor)
##  Lisina Proteina emmean   SE df lower.CL upper.CL
##  F      φ           541 16.2 24      507      574
## 
## Results are averaged over the levels of: Vaca, Periodo 
## Confidence level used: 0.95
plot(mejor)

Deberíamos usar Lisina F que es al 0.5 y Proteína φ que es 12 para tener el mejor rendimiento de leche que seria aproximadamente 541 galones; es claro que en la medición real de estas dos combinaciones dio 534(Fφφ), pero el valor de 541 es un valor inferido de la regresión lineal.

Por eso es muy importante validar los supuestos, para validar que tan confiable es el modelo.

2.5.1.1 Normalidad de Residuos

normalidad=shapiro.test(resid(modelo))
print(normalidad)
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(modelo)
## W = 0.98663, p-value = 0.8467
#Gráfica de Probabiliad Normal
qqnorm(resid(modelo), main= "Gráfica de Probabilidad para los Residuales del Modelo", xlab="Cuantiles Teoricos", ylab = "Cuantiles de muestra")
qqline(resid(modelo))

Como es observarse, se acepta la hipótesis nula de la Prueba de Shapiro Wilks, y se concluye que los residuales provienen de una problación con Distribución de Probabilidad Normal con μ=0μ=0, conclusion que es confirmada por la Gráfica de Probabilidad Normal.

2.5.1.2 Prueba de Homocedasticidad

La Prueba de Homocedasticidad, o de Varianzas Iguales, se realiza mediante la Prueba de Bartlett, misma que para el caso particular, se aplicará al factor de tratamiento y al factor de bloque denominados proteína..

La Prueba de Bartlett tiene las siguientes hipótesis:

H0:σ2i=σ2j=ConstanteH0:σi2=σj2=Constante

H1:σ2i≠σ2j≠Constante

homocedasticidad_lisina=bartlett.test(resid(modelo)~Lisina,data=df)
print(homocedasticidad_lisina)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  resid(modelo) by Lisina
## Bartlett's K-squared = 6.4864, df = 6, p-value = 0.371
homocedasticidad_proteina=bartlett.test(resid(modelo)~Proteina,data=df)
print(homocedasticidad_proteina)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  resid(modelo) by Proteina
## Bartlett's K-squared = 5.7725, df = 6, p-value = 0.4492