Cálculo de tamaño muestral

El cálculo del tamaño muestral es una etapa clave en el diseño experimental, ya que permite determinar cuántas réplicas son necesarias por tratamiento para garantizar que los resultados del estudio sean estadísticamente significativos y tengan suficiente potencia para detectar efectos relevantes. Dependiendo del enfoque adoptado, puede buscarse minimizar la varianza de un estimador, maximizar la eficiencia bajo restricciones presupuestarias o garantizar la detección de diferencias mínimas entre tratamientos con un nivel de confianza preestablecido.

A lo largo del desarrollo de esta aplicación, también se construyó una librería en R (dexp)con el objetivo de poner a disposición de la comunidad académica y profesional las funciones utilizadas en los cálculos. Esta librería permite aplicar los métodos implementados sin necesidad de interactuar con la interfaz gráfica de la app, siendo especialmente útil para usuarios que prefieren trabajar desde scripts reproducibles o integrarlos a otros proyectos estadísticos en R.

Cabe destacar que, en la aplicación, algunos parámetros avanzados como el número de iteraciones, la precisión de convergencia en algoritmos iterativos utilizados para el cálculo de la potencia y otros pármétros necesarios como el K en HHM no se muestran directamente en la interfaz para facilitar la usabilidad. Sin embargo, estos parámetros están disponibles en las funciones de la librería subyacente, lo cual ofrece a los usuarios más experimentados la posibilidad de ajustar con mayor precisión los métodos empleados, explorar su comportamiento y realizar estudios de sensibilidad.

A continuación se presentan distintos métodos para la determinación del tamaño muestral, adaptados a escenarios con o sin costos variables por tratamiento, así como modelos con componentes de varianza y métodos basados en pruebas post-hoc.

Determinación del Tamaño de Muestra sin Costo Variable por Tratamiento

Teoría

Asumiendo que los costos por tratamiento (\(c_i > 0\)) son variables, bajo la restricción:

\[ \sum_{i=1}^t c_i r_i = C \quad \text{(Restricción presupuestaria)} \]

y considerando que:

\[ Var(\text{MELI}(L)) = \sum_{i=1}^t \frac{\lambda_i^2 \sigma_i^2}{r_i} \quad \text{(Ecuación 5.25)} \]

donde \(L = \sum_{i=1}^t \lambda_i \mu_i\) y

\[ \text{MELI}(L) = \sum_{i=1}^t \lambda_i \bar{y}_i \quad \text{(Ecuación 5.26)} \]

Minimización Lagrangiana

Para minimizar \(Var(\text{MELI}(L))\) sujeto a la restricción presupuestaria, construimos el lagrangiano:

\[ Q = \sum_{i=1}^t \frac{\lambda_i^2 \sigma_i^2}{r_i} + \phi \left( \sum_{i=1}^t c_i r_i - C \right) \quad \text{(Ecuación 5.27)} \]

Al resolver lo anterior, obtenemos:

\[ r_i = \frac{|\lambda_i| \sigma_i}{\sqrt{\phi c_i}} \quad \text{(Ecuación 5.28)} \]

donde el multiplicador de Lagrange es:

\[ \phi = \frac{\left( \sum_{i=1}^t |\lambda_i| \sigma_i \sqrt{c_i} \right)^2}{C^2} \]

Recomendación: Los coeficientes \(\lambda_i\) deben expresarse como fracciones para simplificar el cálculo de \(r_i\).

Asignación Proporcional (Ecuación 5.29)

Si en lugar de optimizar con costos, asignamos las observaciones de forma proporcional a las desviaciones estándar (con \(n\) fijo):

\[ \tilde{r}_{i} = \frac{n \sigma_i}{\sum_{s=1}^t \sigma_s}; \quad i = 1, \dots, t \quad \text{(Ecuación 5.29)} \]

Interpretación: Los tratamientos con mayor variabilidad (\(\sigma_i\)) reciben más réplicas.

Ejemplo

Supongamos que un agricultor tiene los siguientes datos referentes a la producción en toneladas por hectárea de cuatro variedades de caña de azúcar:

Se observa que hay una proporcionalidad en las desviaciones estándar, entonces por (5.29) los tamaños de muestra adecuados para cada variedad serían:

\[ \tilde{r}_1=\frac{(5*4)*(6.27)}{6.27+ 9.57+ 12+ 3.32}= 4.0243\approx4 \]

\[ \tilde{r}_2=\frac{(5*4)*(9.57)}{6.27+ 9.57+ 12+ 3.32}= 6.1424\approx6 \]

\[ \tilde{r}_3=\frac{(5*4)*(12)}{6.27+ 9.57+ 12+ 3.32}= 7.7021\approx8 \]

\[ \tilde{r}_4=\frac{(5*4)*(3.32)}{6.27+ 9.57+ 12+ 3.32}= 2.1309\approx2 \]

Las réplicas asignadas por cada uno de los tratamientos son \(4,6,8,2\) respectivamente.

Determinación del Tamaño de Muestra con Costo Variable por Tratamiento

Teoría

Asumiendo que los costos por tratamiento (\(c_i > 0\)) son variables, bajo la restricción:

\[ \sum_{i=1}^t c_i r_i = C \quad \text{(Restricción presupuestaria)} \]

y considerando que:

\[ Var(\text{MELI}(L)) = \sum_{i=1}^t \frac{\lambda_i^2 \sigma_i^2}{r_i} \quad \text{(Ecuación 5.25)} \]

donde \(L = \sum_{i=1}^t \lambda_i \mu_i\) y

\[ \text{MELI}(L) = \sum_{i=1}^t \lambda_i \bar{y}_i \quad \text{(Ecuación 5.26)} \]

Minimización Lagrangiana

Para minimizar \(Var(\text{MELI}(L))\) sujeto a la restricción presupuestaria, construimos el lagrangiano:

\[ Q = \sum_{i=1}^t \frac{\lambda_i^2 \sigma_i^2}{r_i} + \phi \left( \sum_{i=1}^t c_i r_i - C \right) \quad \text{(Ecuación 5.27)} \]

Al resolver lo anterior, obtenemos:

\[ r_i = \frac{|\lambda_i| \sigma_i}{\sqrt{\phi c_i}} \quad \text{(Ecuación 5.28)} \]

donde el multiplicador de Lagrange es:

\[ \phi = \frac{\left( \sum_{i=1}^t |\lambda_i| \sigma_i \sqrt{c_i} \right)^2}{C^2} \]

Recomendación: Los coeficientes \(\lambda_i\) deben expresarse como fracciones para simplificar el cálculo de \(r_i\).

Ejemplo

Supongamos que un agricultor tiene los siguientes datos referentes a la producción en toneladas por hectárea de cuatro variedades de caña de azúcar y sus respectivos costos:

Se observa que hay tenemos un costo por cada tratamiento y un presupuesto total, entonces por (5.28) los tamaños de muestra adecuados para cada variedad serían:

Iniciando por el calculo de \(\Phi\)

\[ \Phi= \sqrt{\frac{\left(\left(\frac{1}{4}*6.27*\sqrt{1000}\right)+\left(\frac{1}{4}*9.57*\sqrt{200}\right)+\left(\frac{1}{4}*12*\sqrt{700}\right)+\left(\frac{1}{4}*3.32*\sqrt{1100}\right)\right)^2}{50000^2}}=1.448629\times10^{-05} \]

\[ r_1=\frac{\frac{1}{4}*6.27}{\sqrt{1.448629\times10^{-05}*1000}}=13.0235\approx13 \]

\[ r_2=\frac{\frac{1}{4}*9.57}{\sqrt{1.448629\times10^{-05}*200}}=44.4486\approx44 \]

\[ r_3=\frac{\frac{1}{4}*12}{\sqrt{1.448629\times10^{-05}*700}}=29.7915\approx30 \]

\[ r_4=\frac{\frac{1}{4}*3.32}{\sqrt{1.448629\times10^{-05}*1100}}=6.5751\approx7 \]

Las réplicas asignadas por cada uno de los tratamientos son \(13,44,30,7\) respectivamente.

Asignación de tratamientos y réplicas con Función de Costos y Varianza Máxima

Teoría

En el modelo de componentes de varianza tanto el número de tratamientos \(t\) como el número de réplicas \(r\) son variables y sus estimaciones están ligadas con el control de dichas varianzas. Un criterio usual para elegir los valores de \(r\) y \(t\) es el de minimizar los costos en la estimación de la media \(\mu\). Una medida de la cantidad de información disponible para estimar \(\mu\) es la varianza de la media muestral dada por:

\[ v(\bar{y}..) = \frac{\sigma^2}{rt} + \frac{\sigma_A^2}{t} \]

En el caso de un diseño C.A. (Completamente Aleatorizado).

El problema se reduce a encontrar los valores de \(r\) y \(t\) que minimicen la función de costos dada por:

\[ C = C_1 t + C_2 t r \]

para una varianza \(v(\bar{y}..)\) fija, donde \(C_1\) es el costo por unidad de tratamiento y \(C_2\) es el costo por unidad experimental. La solución matemática es, según Mendenhall (1968):

\[ t = \frac{1}{v(\bar{y}..)} \left( \hat{\sigma}_A^2 + \sqrt{\frac{\hat{\sigma}_A^2\hat\sigma^2C^2}{C^1}} \right) \]

y

\[ r = \sqrt{\frac{\hat\sigma^2 C_1}{\hat{\sigma}_A^2 C_2}} \]

Ejemplo

Un estudio genético en ganado, consistió en seleccionar aleatoriamente varios machos (toros) apareados con grupos separados de hembras. Cuando nacieron los terneros, se midieron los pesos iniciales como medida en un estudio de pesos hereditarios (estudios de progene). En la tabla a continuación se presentan los pesos al nacer de los terneros de cada uno de cinco grupos de apareamiento.

Obteniendo la tabla ANOVA

            Df Sum Sq Mean Sq F value Pr(>F)  
Toro         4   6070  1517.6   3.646 0.0155 *
Residuals   30  12486   416.2                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

De la anterior tabal pondemos obtener los valores de \(\hat\sigma_e^2\) y de \(\hat\sigma_A^2\)

\[ \hat\sigma_e^2=\text{CME}=416.2 \]

\[ \hat\sigma_A^2=\frac{\text{CMA}-\text{CME}}{r_0}=\frac{(1517.6-416.2)}{6.97}=158.02 \]

Donde,

\[ r_0=\frac{\left(35-\frac{\left(8^2+6^2+6^2+7^2+8^2\right)}{35}\right)}{4}=6.97 \]

Suponiendo una varianza máxima \(V(\bar{y}_{..})=43.49\), \(C1 = 150000\) y \(C2 = 50000\), se encuentra que:

\[ t = \frac{1}{43.49} \left( 158.015 + \sqrt{\frac{(158.015)(416.21409)(50000)}{150000}} \right)=7.04 \]

y

\[ r = \sqrt{\frac{(416.21409) (150000)}{(158.015) (50000)}}=3.35 \]

Para una varianza de la media muestral no mayor de 43,49, deberían seleccionarse 7 toros y 3 terneros por cada toro, asumiendo que el costo experimental de cada toro es de \(150000\) y el de cada ternero es de \(50000\).

Estimación S1 y df1

Teoría

Ocasionalmente el experimentador no define estimadores de varianza en la forma \(S_1^2\), pero conoce “algo” acerca de los órdenes de magnitud para relacionar la información con los límites superior e inferior de la desviación estándar dentro de las cantidades máximas \(S_1\) y \(df_1\).

Harris, Horvitz & Mood (1948) proponen un procedimiento simple: Primero el experimentador se cuestiona por los límites inferior y superior (\(S_I\), \(S_S\)) de la desviación estándar, suponiendo que él desea estimar 7% y 12% de la media para los límites superior e inferior, y la desviación estándar de la media es 30, entonces \(S_I = (0{,}7)(30) = 2{,}1\) y \(S_S = (0{,}12)(30) = 3{,}6\).

La estimación de la desviación estándar es el promedio de los dos estimadores:

\[ S_1 = \frac{S_I + S_S}{2} = 2{,}85 \]

Para obtener \(df_1\) es necesario que el investigador tenga alguna “confianza” sobre los estimadores. Se calculan los valores de \(\sqrt{\chi^2_{(0{,}1)} / \chi^2_{(0{,}9)}}\) para varios grados de libertad, el valor más cercano al cociente \(S_S / S_I\) se considera como los grados de libertad asociados con \(S_1\) (\(df_1\)). Para el ejemplo \(S_S / S_I = 1{,}72\) y de la tabla \(\sqrt{\chi^2_{(12;0{,}1)} / \chi^2_{(12;0{,}9)}} = 1{,}72\), de donde se asume que \(S_1\) tiene asociados 12 grados de libertad, con este valor se estima \(r\).

Ejemplo

En este caso usaremos la función calcular_S1_df1 para hacer el cálculo del ejemplo anterior.

library(dexp)
set.seed(123)
calcular_S1_df1(desviacion_estandar = 30, Si = 0.07, Ss = 0.12)
$S1
[1] 2.85

$grados_libertad
[1] 12

$valor_x
[1] 1.715391

$error_relativo
[1] 0.000645021

Así, la desviacón estandar estimada es de 2.85 con 12 grados de libertad. El valor de 1.715391 corresponde al cociente de la raiz cuadrada de las chi-cuadrado de la parte superior y dado que esto fue estimado por medio de simulaciones el error relativo es de 0.0645021%

Tamaño r por potencia

Teoría

En arreglos de DCA, a menudo, el investigador está interesado en determinar el número de réplicas que le permitan detectar diferencias significativas entre los tratamientos, es decir, determinar si hay o no evidencia para que con la prueba \(F\) del análisis de varianza se rechace o no la hipótesis nula.

La técnica desarrollada en la teoría estadística para decidir sobre el número de réplicas necesarias en un experimento es el cálculo de la potencia de las pruebas estadísticas de interés. En la prueba \(F\) del ANOVA a una vía, el cálculo directo de la potencia es generalmente complejo, pero se han construido algunas gráficas, llamadas curvas características de operación, que permiten estimar un valor para la probabilidad \(\beta\) o error de tipo II. La potencia \(1 - \beta\) se deduce a partir de esta probabilidad.

Tang desarrolló un procedimiento y preparó las tablas necesarias para determinar el número de réplicas a ser usado en un experimento. Asume que la variación de la componente del error se distribuye en forma normal con media cero y varianza \(\sigma^2\). Para obtener el número adecuado de réplicas, se requiere un buen estimador de \(\sigma^2\) y que la diferencia verdadera del efecto medio de tratamiento \(\alpha_i = \mu_i - \mu\) se especifique.

La sensibilidad de la prueba \(F\), o poder de la prueba, denotado por \(1 - \beta\), donde \(\beta\) es la probabilidad del error tipo II, depende de:

  1. El tamaño del test, es decir, de la probabilidad del error tipo I (\(\alpha\)).
  2. De los grados de libertad \((t - 1)\) y \(t(r - 1)\) (en el caso de una vía de clasificación).
  3. Del parámetro de no centralidad:

\[ \phi^2_{\beta,\alpha} = \frac{r \sum \alpha_i^2}{t \sigma^2} \tag{5.19} \]

de la distribución \(F\) no central, donde los \(\alpha_i\) son los valores verdaderos del efecto de tratamientos, especificado bajo la hipótesis alterna. En el procedimiento general se especifica \(\alpha\), \(1 - \beta\) y \(\phi/r\) y se formula la pregunta: ¿Cuántas réplicas, \(r\), son necesarias para detectar, con probabilidad \(1 - \beta\), diferencias de tratamientos especificados por \(\phi/r\) si se usa una prueba de tamaño \(\alpha\)?

El número de réplicas en el modelo I se obtiene haciendo uso de las gráficas construidas por Pearson y Hartley, las cuales se encuentran reproducidas en varios textos (aunque con ligeras modificaciones), por ejemplo en Kempthorne (1952), Scheffé (1959), Neter, Wasserman & Kutner (1990) y Montgomery (2003). En este texto se presentan las tablas de este último autor. Las curvas fueron construidas para dar el valor \(\beta\) en la ordenada cuando se proponen valores de un parámetro \(\phi\) sobre la abscisa y se asumen valores conocidos de \(\alpha\), \(\nu_1\) y \(\nu_2\). El parámetro \(\phi_{\beta,\alpha}\) se llama parámetro de no centralidad de la distribución \(F\) y es una medida del grado de desigualdad de los \(\alpha_i\).

En este caso, no se dificulta la especificación de \(\alpha\) y \(\beta\). Usualmente se toma \(\alpha = 0.05\) o \(\alpha = 0.10\) como valores razonables para evaluar el riesgo de cometer el error tipo I.

Un paso más difícil es la escogencia de \(\beta\) o \(1 - \beta\), es decir la probabilidad de detectar diferencias cuando existen. Es razonable tomar \(1 - \beta = 0.80\), aunque depende del problema. La parte más difícil es definir \(\phi / r\) porque representa el verdadero estado de la naturaleza.

Ejemplo 5.8

Se va a determinar el número de réplicas para un experimento similar al del ejemplo 5.1 suponiendo que la potencia no debe ser inferior a 0.80. Asignando los mismos valores del ejemplo citado a los parámetros, se tiene:

\[ \alpha = 0.05, \quad \nu_1 = t - 1 = 3, \quad \nu_2 = t(r - 1) = 16, \quad \hat{\sigma}^2 = 10.35 \]

\[ \sum_{i=1}^t \hat{D}_i^2 = 1.1342 + 16.2812 + 14.94 + 0.801 = 33.156 \]

Entonces:

\[ \phi^2 = \frac{5 \cdot 33.156}{4 \cdot 10.35} = 4.004 \quad \Rightarrow \quad \phi = 2.001 \]

Según la figura A.1 del apéndice, en la gráfica con \(\nu_1 = 3\) y \(\alpha = 0.05\), se localiza \(\phi = 2.001\), se sube hasta cortar la línea \(\nu_2 = 16\) y se lee la probabilidad \(\beta \approx 0.16\). Por lo tanto, la potencia es \(1 - \beta \approx 0.84\).

Si se quiere una potencia menor, se puede disminuir \(r\). Al suponer \(r = 4\), se calcula \(\phi = 1.79\) y \(\beta \approx 0.21\), por lo tanto la potencia será \(0.79\), que no cumple la condición inicial. Se recomienda tomar 5 réplicas para futuros estudios similares.

En la práctica, cuando no se conocen los parámetros, se debe estimar primero \(\sigma^2\). Si se conocen estudios previos, puede usarse como estimador, de lo contrario, se debe realizar una suposición razonable.

Si fuera posible proponer valores para los \(\alpha_i\), se calcularía como en el ejemplo anterior. Lo más frecuente es no tener esos valores, por lo tanto, se recurre a especificar la diferencia mínima entre los extremos de los tratamientos (el mejor y el peor), es decir, \(\Delta = \alpha_{\max} - \alpha_{\min}\).

Si se asume \(\alpha_{\max} = -\alpha_{\min}\) entonces \(\alpha_{\max} = \frac{\Delta}{2}\), \(\alpha_{\min} = -\frac{\Delta}{2}\) y \(\alpha_i = 0\) para los demás. Así, se tiene:

\[ \Phi^2 = \frac{2r \cdot \Delta^2 / 4}{t \sigma^2} = \frac{r \Delta^2}{2t \sigma^2} \tag{5.20} \]

Como la prueba \(F\) es creciente en función de \(\Phi\), se toma el menor valor de \(\Phi\) para todas las diferencias.

En los casos prácticos, se parte de un valor tentativo de \(r\), se calcula \(\Phi\), se consulta la tabla y se estima si la potencia es la deseada. Si es menor, se incrementa \(r\), si es mayor, se disminuye.

Ejemplo

Considere un estudio donde se detecta como significativa una diferencia entre medias igual o mayor de 3 kilogramos, es decir, \(\Delta = 3\). Supóngase que se estimó la varianza en \(\hat{\sigma}^2 = 10.35\). Hay \(t = 4\) tratamientos por comparar y se requiere una potencia de 0.80. Entonces:

\[ \Phi^2 = \frac{9r}{(8)(10.35)} = 0.1087r. \]

Iniciando con \(r = 15\), en la figura A.1 del apéndice con \(\nu_1 = 3\) y \(\alpha = 0{,}05\), se localiza \(\beta = 0{,}20\) y para \(\nu_2 = (4)(14) = 56\) se lee \(\Phi \simeq 1{,}72\). Entonces \(\Phi^2 = 2{,}96\) y el \(r\) despejado es 27.22. Se repite el proceso con \(r = 27\), de modo que:

\[ \Phi = \sqrt{(0.1087)(27)} = 1{,}71 \quad \text{y} \quad \nu_2 = (4)(26) = 104. \]

Estos valores dan una probabilidad \(\beta \simeq 0{,}20\), se concluye que se requieren 27 individuos como mínimo por tratamiento si se desea una potencia de 0.80 para la prueba \(F\) y, asumiendo diferencias de 3 kilogramos o mayores entre medias poblacionales como significativas.

\(\Delta\) \(\sigma^2\) \(1 - \beta\) \(r\) \(\Phi\) \(\beta\)
3 10.35 0.80 27 1.71 0.20
2 10.35 0.80 61 1.72 0.20
3 8.50 0.80 22 1.72 0.20
3 10.35 0.90 35 1.95 0.10
4 8.50 0.82 14 1.81 0.18

Tabla 5.13. Valores de \(r\) para diferentes valores de los parámetros \(\Delta\), \(\sigma\) y \(1 - \beta\).

Al crecer el número de réplicas \(r\), también crece la potencia. Para una potencia fija, se puede disminuir \(r\) si se aumenta \(\Phi\). Pero el parámetro \(\Phi\) depende básicamente del cociente \(\frac{\Delta}{\sigma}\), el cual puede aumentar ya sea porque la varianza es pequeña o porque la diferencia significativa se asume grande.

Diferencias grandes entre las \(\alpha_i\) son fáciles de detectar con pocos datos. Como no hay mucha precisión en cuanto a la varianza estimada, es aconsejable investigar varios tamaños de muestra dentro de un rango probable de valores de \(\sigma^2\), antes de decidirse por un tamaño definitivo.

Ahora implementaremos la función que hicimos.

library(dexp)

resultado <- encontrar_r_minimo_Potencia(t = 4, rho = 0.5, dif = 3, potencia_objetivo = 0.8, sigma2_ = 10.35)
resultado$r_optimo
[1] 22
resultado$grafico

Método de Harris–Hurvitz–Mood (HHM)

Teoría

Este método determina el número de réplicas requerido para obtener significancia en una proporción específica de experimentos, donde diferencias grandes o mayores que algunos valores de \(d\) preestablecidos existen. Se asume que los valores dentro de cada población se distribuyen en forma normal con varianza común para todas las poblaciones, las cuales tienen como estimador \(S_1^2\) y \(df_1\) grados de libertad.

Conociendo los valores de \(S_1^2\), \(df_1\) y \(d\), el número de réplicas requerido, con un nivel \(\alpha\) de significancia, se obtiene a partir de la expresión

\[ r \;=\; 2\,\bigl(df_2 + 1\bigr)\,\Bigl(\frac{K' \,S_1}{d}\Bigr)^2 \]

donde:

  • \(K'\) es el valor de la tabla A.9 del apéndice.
  • \(df_2\) son los grados de libertad estimados del segundo estimador de la varianza población \(S_2^2\).

Si \(r\) es demasiado pequeño para proporcionar una estimación confiable de \(df_2\), debe emplearse el menor valor posible de \(df_2\).

Ejemplo

Los datos de un experimento sobre cantidad de grasa absorbida por 64 buñuelos para diferentes tipos de grasa se muestran en la Tabla 5.15.

Cantidades de grasa Medias 4 3 2 6 1 5 8
7 161 24 21 17 15 11 4 1
8 162 23 20 16 14 10 3 —
5 165 20 17 13 11 7 — —
1 172 13 16 6 4 — — —
6 176 9 6 2 — — — —
2 178 7 4 — — — — —
3 182 3 — — — — — —

Suponga que se encontró \(S_1^2 = 141.6\) con \(df_1 = 40\) grados de libertad. Al asumir que \(d = 20\), \(1-\beta = 80\%\) e incluyendo solamente 6 tratamientos en el experimento, se obtiene

\[ r \;=\; \frac{2\,(141.6)\,(0.322)^2\,(60 + 1)}{400} \;=\; 4.48 \]

Observe que si \(df_2 = 60\), entonces se sobreestima los grados de libertad; en tanto que si \(df_2 = 25\), es decir, subestimando los grados de libertad, se obtiene

\[ r \;=\; \frac{2\,(141.6)\,(0.502)^2\,(25 + 1)}{400} \;=\; 4.64 \]

Por lo tanto, se necesitan 5 réplicas para alcanzar los resultados deseados en este experimento.

Si \(1 - \beta = 0.95\), entonces son necesarias más réplicas.

Ahora implementando nuestra librearía, tenemos:

library(dexp)
resultados <- calcular_r_minimo_HHM(
  t = 6,
  d = 20,
  r_inicial = 3,
  S1_sq = 141.6,
  df1 = 40,
  alpha_ = 0.05,
  beta_ = 0.8
)

resultados$r      # número estimado de réplicas
[1] 5
resultados$K_     # valor de K = a / s1 obtenido por simulación
[1] 0.5348283
resultados$df2    # grados de libertad del error final
[1] 22

Se implemento una función k_tabla_mc la cual nos permite calcular el K independientemente de los valores de df1 y df2, es decir, ya la tabla la cual fue calculada para ciertos valores específicos de grados de libertad ya no es una limitación.

Metodo de tukey

El método de Tukey para el cálculo del tamaño muestral de las réplicas por tratamiento tiene como objetivo asegurar que, al momento de realizar el diseño experimental, sea posible detectar una diferencia mínima significativa entre los tratamientos en la variable de interés al aplicar comparaciones pareadas mediante pruebas post-hoc.

En otras palabras, si al aplicar estas pruebas no se detecta una diferencia significativa, se debe a que la diferencia establecida previamente como mínima detectable no se logró identificar, lo que está directamente relacionado con el tamaño muestral calculado inicialmente.

Para aplicar esta prueba es necesario contar con información a priori sobre la varianza del error experimental. Si no se dispone de dicha información, y por lo tanto tampoco de los grados de libertad asociados, se procede a estimarla mediante la función calcular_S1_df1().

Parámetros:

Teoría

Un intervalo de confianza de \((1 - \alpha) \cdot 100\) de longitud \(\leq 2d\) para la diferencia de cualquier par de medias, cuando se tiene un conjunto de \(t\) tratamientos, se obtiene a partir de la expresión del siguiente proceso:

Si hay \(t\) medias (asociadas con los tratamientos), se hacen comparaciones dos a dos, es decir:

\[ \frac{\bar{Y}_{\max} - \bar{Y}_{\min}}{\sqrt{CME/r}} \sim q_{t;df_2} \]

Sea \(P_0 \leq P\) la longitud del intervalo que cubre todas las diferencias \(\leq 2d\), entonces:

\[ P_0 = P\left(2S q_{1 - \alpha}/\sqrt{r} \leq 2d\right) \]

donde \(S = \sqrt{CME}\) y \(q_{1 - \alpha}\) es el límite en la tabla A.12 del apéndice de rangos studentizados. Desde esta expresión se encuentra que:

\[ P_0 = P\left(S^2 \leq \frac{d^2 r}{q_{(1 - \alpha)}^2}\right) \]

En esta expresión, al dividir por un \(S^2\) previo, \(S_1^2\), para obtener una \(F\), se sigue:

\[ P_0 = P\left(\frac{S^2}{S_1^2} \leq \frac{d^2 r}{q_{(1 - \alpha)}^2 S_1^2}\right) \]

Con \(\frac{S^2}{S_1^2} \sim F(df_2, df_1)\), se obtiene finalmente:

\[ r = F_{(df_2, df_1; 1 - \alpha)} \cdot \frac{q^2_{t; df_2; 1 - \alpha/2}}{d^2} \]

Ejemplo

Retomando el ejemplo del libro 5.11, sea \(S_1^2 = 141{,}6\), \(df_1 = 40\) y \(d = 20\).
Los valores de \(gl_E\) y \(q_{t; df_2; 1 - \alpha}\) dependen del diseño y de los valores de \(1 - \beta\).
Si se supone \(S_1\) como en el diseño completamente aleatorizado con \(r = 6\), \(1 - \beta = 0{,}9\) y \(df_2 = 30\), entonces:

\[ \begin{align*} gl(\text{Total}) &= 35, \\ gl(\text{Tratamientos}) &= 5, \\ gl(\text{Error}) &= 30. \end{align*} \]

Así, \(F_{(30;40;0{,}10)} = 1{,}54\) y \(q_{(5;30;0{,}10)} = 4{,}30\), donde \(q\) corresponde al valor apropiado en la tabla A.12 de rangos studentizados.

Con estos resultados y de la ecuación (5.24), se obtiene:

\[ r = \frac{(141{,}6)(4{,}30)^2(1{,}54)}{400} = 10{,}08 \]

Para este caso, el estimador es sesgado. Si se desea garantizar que la longitud del intervalo sea \(2A\), entonces de (5.18) se encuentra:

\[ r = \frac{2\sigma^2 t_{1-\alpha}^2 \Gamma^2\left( \frac{r}{2} \right)}{A^2 (r - 1) \Gamma^2\left( \frac{r - 1}{2} \right)} \]

Como el valor de \(df_2\) fue subestimado, se toma \(r = 9\), entonces \(df_2 = 48\),
\(q_{(5;48;0{,}90)} = 4{,}2\), \(F_{(48;40;0{,}10)} = 1{,}48\) y \(r = 9{,}2\).

Como el valor de \(r > 9\), entonces con \(r = 10\), va a ser lo suficientemente grande para esta propuesta,
se puede obtener un intervalo de confianza con 10 réplicas y el 95 % de confianza, el cual es mayor que \(2d\) en longitud en el 90 % de los experimentos.

Ahora vamos a usar la librería que creamos:

library(dexp)
calcular_r_MT(
  T_ = 6,
  D = 20,
  ro = 6,
  S1 = sqrt(141.6),
  df1 = 40
)
$r
[1] 10

$A
[1] 9.488479

$r_i
[1] 9

$A_i
[1] 10.11104

$lista_r
 [1]  5  6  7  8  9 10 11 12 13 14 15

$valores_A
 [1] 15.081301 13.172866 11.857586 10.878156 10.111039  9.488479  8.969733
 [8]  8.528597  8.147330  7.813420  7.517754

$posicion_escogida
[1] 5

El valor r representa el tamaño de réplica calculado mediante el método de Tukey, siendo A la diferencia mínima que se desea detectar con dicho tamaño. Por otro lado, r_i es el tamaño de réplica cuyo valor se aproxima más a lograr la diferencia mínima requerida, denominada A_i. Las listas lista_r y valores_A contienen valores de tamaño de réplica y las correspondientes diferencias mínimas detectables (A), calculadas tomando como referencia el valor de r.

Según el método de Tukey, se concluye que para detectar una diferencia mínima de 20 unidades se necesita un tamaño de réplica de 10. No obstante, esto implica que la diferencia real podría estar subestimada, por ejemplo, en un valor como 9.48. Esto puede interpretarse también como que, al emplear un tamaño de 10 réplicas, el contraste será más estricto para rechazar la hipótesis nula \(H_0\). En cambio, utilizar el valor r_i proporciona un enfoque ligeramente más flexible, lo que puede facilitar la detección de la diferencia deseada con un balance adecuado entre error tipo I y tipo II.

Número de réplicas en el modelo de efectos aleatorios

El cálculo del número de réplicas en un diseño experimental con modelo de efectos aleatorios es fundamental cuando se desea asegurar una potencia estadística adecuada en pruebas ANOVA. En este contexto, el objetivo es determinar cuántas observaciones por tratamiento son necesarias para detectar diferencias significativas entre niveles de un factor aleatorio, considerando la variabilidad inherente del experimento.

Aunque la potencia de la prueba se basa en la distribución F, su determinación puede simplificarse mediante el uso de curvas características de operación (OC curves), que permiten visualizar la probabilidad de cometer errores tipo II (\(\beta\)) en función de diferentes configuraciones del diseño.

Este enfoque es especialmente útil en estudios donde los factores no son controlados directamente sino seleccionados aleatoriamente de una población más amplia (por ejemplo, lotes, sujetos, ubicaciones), y se requiere garantizar que el diseño tenga suficiente sensibilidad para detectar efectos reales entre tratamientos o unidades experimentales.

Teoría

En el caso del ANOVA de una vía de clasificación:

\[ Y_{ij} = \mu + T_i + \varepsilon_{ij} \]

  • \(\mu\): media general
  • \(T_i\): efecto diferencial del tratamiento
  • \(\varepsilon_{ij}\): error de la unidad observacional

Nosotros estimábamos a \(T_i\) como si fuera un valor fijo, pero cuando tenemos efectos aleatorios:

\[ \varepsilon_{ij} \sim \mathcal{N}(0, \sigma^2_\varepsilon) \]

Ahora \(T_i\) va a ser aleatorio de individuo a individuo:

\[ T_j \sim \mathcal{N}(0, \sigma^2_T) \]

Así tendremos dos efectos de varianza, entonces la idea se hace grande cuando vemos el diseño completamente al azar desde el punto de vista de los efectos aleatorios, ya no nos va a interesar evaluar la misma hipótesis.

Hipótesis en efectos fijos:

\[ H_0: \; T_1 = T_2 = \dots = T_K = 0 \]

Si no que:

Hipótesis en efectos aleatorios:

\[ H_0: \; \sigma^2_T = 0 \quad \text{(La varianza de los efectos aleatorios sea estadísticamente cero)} \]

El cálculo es el mismo, solo que la hipótesis cambia.

Aclaración:
Que \(T_k = 0 \Rightarrow \sigma^2_T = 0\), pero no al revés.
Es decir, si \(\sigma^2_T = 0\) implica que \(T_k\) son iguales, pero no a cero.

¿Cuándo decidimos usar cuál?

En efectos fijos, los \(T\) decimos que son parámetros fijos, que se espera que afecten a la variable dependiente de forma directa.
En cambio, los efectos aleatorios me ayudan para explicar la variabilidad no explicada de las explicativas (punto de vista aleatorio).

Ejemplo:
Si se tiene un estudio de diferentes fertilizantes y queremos evaluar el rendimiento en una cosecha, el investigador está pensando en elegir solo tres tipos de fertilizante.
Para él, es un valor fijo dado que lo decidió bajo experiencia y quiere evaluar entre esos tres.
→ Efecto fijo

Ahora, si en el mismo ejercicio se tiene un grupo más amplio de fertilizantes (20, 30, …, X) y elijo aleatoriamente 3 por medio de un muestreo, entonces estoy pensando en efectos aleatorios.
En ese caso, me interesa estimar la variabilidad del rendimiento de lo que estamos midiendo, causado por los diferentes tipos de fertilizantes.

Y en ese intento, la varianza asociada a los fertilizantes no se ve influenciada por la manera en que se eligen aleatoriamente los fertilizantes que se iban a usar en el estudio.

→ La diferencia de un modelo a otro reside en la concepción del experimento, pero en el contraste de las diferencias, o desde el punto de vista estadístico, es exactamente la misma tanto en efectos fijos como en efectos aleatorios.

Otra diferencia es que en E.F. a mí me interesa comparar niveles específicos, i.e., el fertilizante 1 vs. el 2, se comparan niveles específicos.

En cambio, en E.A. es controlar o evaluar qué tanto podrían cambiar la variabilidad en el rendimiento, o en la variable de respuesta por efecto del tratamiento.

→ Para contrastar el ejercicio se calcula el estadístico \(F\) ya sea visto en E.F. o E.A.:

\[ F = \frac{CM_{Trat}}{CM_{Res}} \]

Una cosa importante a tener en cuenta es que en el numerador de \(F\) hay una combinación entre los \(\sigma^2_T\) (variabilidad de tratamientos) y \(\sigma^2_\varepsilon\) (variabilidad del error experimental).
En cambio, en \(CM_{Res}\) solo se tiene \(\sigma^2_\varepsilon\).

→ Eso sería lo primero.

Ahora el ejercicio para calcular \(CM_{Trat}\) para E.A. es el siguiente:

Suponemos que:

\[ SC_{Tot} = SC_{Trat} + SC_{Res} \]

Sabemos que son chi-cuadrados, entonces:

\[ CM_{Trat} = \frac{SC_{Trat}}{K - 1} \sim \chi^2 \Rightarrow \mathbb{E}(CM_{Trat}) = \sigma^2_\varepsilon + r \cdot \sigma^2_T \]

Donde \(K\) es el número de tratamientos.
Y para \(CM_{Res}\):

\[ \mathbb{E}(CM_{Res}) = \sigma^2_\varepsilon \]

Al final, va a pasar que la estadística \(F\) bajo la hipótesis nula estará a tener distribución \(F\), pero cuando no estoy bajo elsupuesto de la hipótesis nula, i.e., si \(\sigma^2_T \neq 0\), el resultado de \(F\) va a depender obligatoriamente de \(\sigma^2_T\).

→ Las curvas características de operación (curvas OC) para el anterior escenario me dicen que a partir de la distribución exacta de la estadística \(F\), cuando el \(F\) no es centrado o no central, me dice que va a haber un cambio proporcional entre \(F\) y \(\sigma^2_T\).
Entonces, lo que uno hace es una suerte de cociente:

\[ \rho = \frac{\sigma^2_T}{\sigma^2_\varepsilon} \]

→ Para el ejercicio uno hace lo siguiente:

  1. Simulación de cómo va a ser la potencia.
    Hace el cálculo por el ANOVA:

    • Fijamos el número de tratamientos (\(t\))
    • Fijamos el número de réplicas (\(r\))
    • Fijamos \(\sigma^2_\varepsilon\)
    • Fijamos \(\beta\)

Supongamos que ya tenemos claro el modelo:

\[ y_{ij} = \mu + T_j + \varepsilon_{ij} \]

→ De \(\sigma^2_\varepsilon\) y \(\beta\) se puede recuperar \(\sigma^2_T\).
\(T_j\) tiene que ser aleatorio, entonces se simula vía rnorm() al igual que \(\varepsilon_{ij}\):

  • \(T_j \sim \mathcal{N}(0, \sigma^2_T)\)
  • \(\varepsilon_{ij} \sim \mathcal{N}(0, \sigma^2_\varepsilon)\)

Se fija el valor de \(\mu\).

Con estas asignaciones generamos un ANOVA:

anova(y ~ trt)
  • Luego se calcula \(F\)-stat
  • Se mira si se rechaza
  • Se supone que \(H_0\) es falso para saber

La potencia de la prueba según \(r\):

Referencia

Referencia Bibliográfica

Melo M., O. O., López P., L. A., & Melo M., S. E. (2007). Diseño de Experimentos: Métodos y Aplicaciones (1a ed.). Pro-Offset Editorial S.A.. Capítulo 5, pp. 189-207.