Parcelas divididas

En esta sección vamos a aprender algo sobre diseños experimentales que contienen unidades experimentales de diferentes tamaños, con diferentes aleatorizaciones. Estos diseños de parcelas divididas son quizás los diseños más incomprendidos en la práctica; por lo tanto, a menudo se analizan de forma incorrecta. Desafortunadamente, los diseños de parcelas divididas son muy comunes, aunque no siempre se llevan a cabo con el propósito o conscientemente. El buen mensaje es que una vez que sabes cómo detectar estos diseños, el análisis es sencillo, solo tenemos que añadir los efectos aleatorios adecuados al modelo, la cual es una propuesta más conveniente a otras que han surgido para su análisis.

Introducción:

Comenzamos con un ejemplo agronómico. Un productor tiene ocho parcelas donde aleatoriza de forma balanceada la aplicación de dos esquemas de fertilización, un control y uno nuevo con Biochar. Además, cada parcela se divide en cuatro subparcelas donde siembra cuatro variedades diferentes de fresas asignándolas aleatoriamente a las subparcelas. El productor está interesado en el efecto del esquema de fertilización y la variedad de fresa sobre la biomasa de la fruta en la primera cosecha. Esto significa que tenemos un total de <8.*4 =32> observaciones si se recoge un dato por subparcela.

El diseño se resume en la tabla 1. Alli las ocho parcelas aparecen en las diferentes columnas. El esquema de fertilización se denota por control (ctrl) y Biochar (Bioc). Las variedades de fresa se etiquetan con las letras A,B,C y D.

La variable respuesta viene dada por la biomasa. Si consideramos los dos factores , fertilizante y variedad, el diseño se parece a un diseño factorial completo “clásico” a primera vista.

Importando datos del GitHub

1

url <- "https://github.com/darghan/unal/raw/refs/heads/main/dfa.xlsx"
temp_file <- tempfile(fileext = ".xlsx")
download.file(url, temp_file, mode = "wb")
dfa <- read_excel(temp_file)
head(dfa)
## # A tibble: 6 × 4
##   plot  fertilizer variety  mass
##   <chr> <chr>      <chr>   <dbl>
## 1 1     control    A         8.9
## 2 1     control    B         9.5
## 3 1     control    C        11.7
## 4 1     control    D        15  
## 5 2     control    A        10.8
## 6 2     control    B        11

Repeticiones por tratamiento 

xtabs(~fertilizer + variety, data = dfa)
##           variety
## fertilizer A B C D
##    control 4 4 4 4
##    new     4 4 4 4

Cuando se considera la parcela y el fertilizante, vemos que por ejemplo la parcela 1 solo contiene control de nivel de fertilizante, mientras que la parcela 3 solo contiene nivel nuevo(Biochar). Obtenemos el mismo patrón que en la tabla 1.

Repeticiones por Parcela 

xtabs(~ fertilizer + plot, data = dfa)
##           plot
## fertilizer 1 2 3 4 5 6 7 8
##    control 4 4 0 4 0 4 0 0
##    new     0 0 4 0 4 0 4 4

Tenga en cuenta que normalmente no podemos recuperar el procedimiento de aleatorización solo de un marco de datos. Realmente necesitamos toda la “historia”. Podemos visualizar los datos con un gráfico de interacción que muestra que la biomasa es mayor en promedio con el nuevo esquema de fertilización. Además, parece haber un efecto de variedad. La interacción no es muy pronunciada (el efecto varietal parece ser consistente entre los dos esquemas de fertilización).

with(dfa, interaction.plot(x.factor = variety,  ,col = 2:3, lty = 2:3,xtick = TRUE,
                            trace.factor = fertilizer, cex.axis = 0.8,type = "b",leg.bty = "o",
                            response = mass,xlab = "variedad de fresa", ylab="promedio de biomasa"))

¿Cómo podemos modelar tales datos? Para establecer un modelo correcto, tenemos que estudiar cuidadosamente el protocolo de aleatorización que se aplicó. Se han utilizado dos aleatorizaciones:

Los esquemas de fertilización se han aleatorizado y aplicado a parcelas. Se seleccionaron aleatoriamente variedades de fresa y se aplicaron a subparcelas, por lo tanto, una unidad experimental para el fertilizante es dada por una parcela de tierra, mientras que para la variedad de fresa, la unidad experimental es dada por una subparcela.Este no habría sido el caso si la aleatorización es completa y en todo el terreno.

El diseño

Este diseño es un ejemplo de diseño de parcela dividida. El esquema de fertilización es el llamado factor de la parcela principal y la variedad de fresa es el factor de la parcela dividida o subparcela. Seguramente la fertilización va en la principal pues es más dificil en la práctica seguramente aleatorizar este factor de modo que caiga “en cualquier parte”, por lo que en la parcela principal suele ir el factor que más complicaciones genera en campo (no es la razón única), sin embargo, otro productor pudo haber establecido en la parcela principal la variedad porque el manejo de cada variedad se hace más fácil cuando se ubican en una sola parcela y no en diferentes sitios del experimento.

Como tenemos dos tamaños diferentes de unidades experimentales, también necesitamos dos términos de error para modelar los errores residuales correspondientes, pues tenemos doble proceso de aleatorización. Necesitamos un término de error para la parcela principal y otro en el nivel de la subparcela. Sea \(y_{ijk}\) la biomasa observada de la \(k\)-ésima replica del \(i\)-ésimo nivel de fertilización y \(j\)-ésimo nivel de variedad. Utilizamos el modelo

\[y_{ijk}=\mu+\phi_i+\varepsilon_{k(i)}+ \upsilon_{j}+(\phi\upsilon)_{ij}+\epsilon_{ijk} \]

donde \(\phi_i\) es el efecto fijo del esquema de fertilización,\(\upsilon_{j}\) es el efecto fijo de la variedad de fresa y \((\phi\upsilon)_{ij}\) es el término de interacción correspondiente. Llamamos \(\varepsilon_{k(i)}\)
al error de la parcela principal. La notación de anidación (en paréntesis) asegura que obtengamos un error de la parcela principal (réplicas anidadas en las fertilizaciones). La notación matemática es un poco engorrosa aquí, la especificación del modelo que se muestra más adelante en R será más fácil. Al final tenemos lo que se llama error de la subparcela (el término de error “usual” que actúa en el nivel de una observación individual).

Utilizamos los supuestos: \[\varepsilon_{k(i)} \sim{N(0,\sigma_{\varepsilon}^2)} \hspace{1cm} \epsilon_{ijk}\sim{N(0,\sigma_{\epsilon}^2)} \]

además de la independencia e idéntica distribución.

La parte del modelo asociada con \(\phi_i+\varepsilon_{k(i)}\) puede considerarse como la “reacción” de una parcela individual en el esquema de fertilización. Todas las propiedades específicas de la parcela se incluyen en el error de la parcela completa. El hecho de que todas las subparcelas de la misma parcela compartan el mismo error de la parcela entera tiene el efecto secundario de que las observaciones de la misma parcela se modelen como datos correlacionados.

La parte \(\upsilon_{j}+(\phi\upsilon)_{ij}+\epsilon_{ijk}\) es la “reacción” de la subparcela en lo que respecta a la variedad,incluyendo una posible interacción con el esquema de fertilización. Todas las propiedades específicas de la subparcela se pueden encontrar en el error de la subparcela.

Si solo se considerara el esquema de fertilización, hacemos un diseño factorial simple, con parcelas como unidades experimentales. La primera parte del modelo es en realidad la ecuación modelo correspondiente del ANOVA unidireccional correspondiente a este diseño. Por otra parte, si sólo consideramos la variedad, podríamos tratar las parcelas como bloques y tendríamos un diseño factorial simple con bloques completos generalizados y al azar, incluyendo un término de interacción; esto es lo que vemos en la segunda parte del modelo.

Ajutando el modelo en R

Para ajustar un modelo de este tipo en R, utilizamos un enfoque de modelo mixto. El error de la parcela principal puede ser fácilmente incorporado con la notación (1 | plot) para indicar la aleatoriedad de este componente. El error de la parcela dividida o subparcela, se incluye automáticamente al igual que en las observaciones individuales. Por lo tanto, terminamos con la siguiente llamada de función.

modelo_pd <- lmer(mass ~ fertilizer * variety + (1 | plot), data = dfa)

y la tabla del anova la obtenemos de:

anova(modelo_pd)
## Type III Analysis of Variance Table with Satterthwaite's method
##                     Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
## fertilizer         137.413 137.413     1     6 68.2395 0.0001702 ***
## variety             96.431  32.144     3    18 15.9627 2.594e-05 ***
## fertilizer:variety   4.173   1.391     3    18  0.6907 0.5695061    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Los datos proveen evidencia a favor de la hipótes de interacción nula (ausencia de interacción), mientras que las dos hipótesis de los efectos principales son rechazadas. Examinemos más de cerca los grados de libertad del denominador para comprender mejor el modelo de parcelas divididas. Observamos que para la prueba del efecto principal de los fertilizantes sólo tenemos 6 grados de libertad (DenDF). ¿Por qué? Porque como se ha descrito anteriormente, básicamente realizamos un diseño completamente aleatorio con ocho unidades experimentales (ocho parcelas) y un factor de tratamiento que tiene dos niveles (control y nuevo-Biochar), por lo tanto, el error tiene 8-1-1=6= (2*(4-1)) grados de libertad (un grado menos para los niveles del factor en (4-1) parcelas). También podríamos argumentar que deberíamos comprobar si la variación entre los diferentes esquemas de fertilización es mayor que la variación entre las parcelas que reciben el mismo esquema de fertilización. Este último es dado por el error de todo-plot teniendo (solo!) grados de libertad, ya que está anidado en el fertilizante.

En contraste, la variedad (y la interacción fertilizante:variedad) se prueban contra el “término de error habitual”: Tenemos un total de 32 observaciones, lo que lleva a 31 grados de libertad. Utilizamos 1 grado de libertad para el fertilizante, 6 para el error de la parcela principal (ver arriba), 3 para la variedad y otros 3 para la interacción fertilizante:variedad. Por lo tanto, los grados de libertad que permanecen para el error de la subparcela son 31-1-6-3-3=18. Otra forma de pensar es que podemos interpretar el experimento en el nivel de la parcela dividida como un diseño de bloques completos aleatorios donde bloqueamos en las diferentes parcelas que utiliza 7 grados de libertad. Tanto el efecto principal de la variedad como el efecto interacción fertilizante:variedad utilizan otros 3 grados de libertad cada uno, de modo que llegamos a 18 grados de libertad para el término error.

Comentarios sobre el diseño en parcelas divididas

Los diseños de parcelas divididas implican el uso simultáneo de unidades experimentales de diferentes tamaños. Las correspondientes modelos implican más de un término de error. Según Casella (2008, p.171), “los experimentos de parcelas divididas son el caballo de batalla del diseño estadístico”. Existe un dicho que expresa que si la única herramienta que tienes es un martillo, entonces todo en el mundo se ve como un clavo. Podría ser justo decir que, a partir de ahora, casi todos los diseños que veas serán una especie de parcela dividida.”

En el ejemplo anterior, tenemos un experimento mucho más preciso para las variedadades que para los fertilizantes (ver los grados de libertad del denominador). Recuerden que a mayor grados de libertdad para el error, más potente es la estimación de las diferencias.

Quizás el aspecto más importante de este diseño es la interacción. Es fácil establecer un diseño de los fertilizantes y un buen experimento para las variedades; la dificultad está en conseguir mirar la interacción entre los dos elementos. El hecho más importante en el análisis es que la interacción entre fertilizantes y variedades está sujeta exactamente a la misma variabilidad que las comparaciones de variedades. Así hemos eliminado un importante fuente de variación, la variabilidad entre parcelas principales. Los efectos de interacción sólo están sujetos a la variabilidad de las subparcelas, es decir, la variabilidad dentro de las parcelas principales.

La idea básica detrás de los diseños de parcelas divididas es muy general. La idea clave es que una unidad (parcela entera o principal) se divide para permitir varias mediciones distintas en la unidad. Estos son a menudo llamados medidas repetidas. Cuando se realizan mediciones repetidas en la misma unidad de observación, es más probable que las mediciones sean similares , por lo tanto, las mediciones en la misma unidad están correlacionadas y esta correlación necesita ser modelada en el análisis.

Si existe interacción, necesitamos explorar las relaciones entre los tratamientos. Debido a su manejabilidad matemática, un enfoque útil es mirar las relaciones entre los tratamientos de subparcelas para cada tratamiento de parcela principal fija. Las interacciones implican ver si tales relaciones cambian de tratamiento a tratamiento en las parcelas principales. El siguiente código muestra el patrón que evidencia la ausencia de interacción (el patrón para las cuatro variedades se mantiene en los dos esquemas de fertilización)

ggplot(dfa, aes(x = variety, y = mass, fill = variety)) +
  geom_bar(stat = "identity", position = "dodge") +
  facet_wrap(~ fertilizer) +
  theme_minimal() +
  labs(title = "Biomasa por Fertilizante y Variedad",
       x = "Variedad",y = "Biomasa")

La fuerza de los diseños/modelos de parcelas divididas es su capacidad para analizar los efectos de las subparcelas y las interacciones entre los efectos de las subparcelas y los efectos de las parcelas principales. Normalmente, se dispone de menos información sobre los efectos de la parcela principal (incluidas sus interacciones entre sí). Si las subparcelas están desequilibradas o desbalanceadas, no es posible realizar un análisis limpio de los efectos de la parcela principal (independientemente de que las parcelas principales estén equilibradas). Los métodos ilustrados hasta ahora requieren subparcelas equilibradas, es decir, no faltan subparcelas (pero si falta toda una parcela principal no hay problema). Aunque existen propuestas para el desbalance, se recomienda a estudiantes en sus primeras etapas de formación en diseño, procurar conformar diseños de este tipo pero equilibrados.

¿Por qué debemos utilizar diseños de parcelas divididas? Normalmente, los diseños de parcelas divididas son adecuados para situaciones en las que uno de los factores sólo puede variar a gran escala. Por ejemplo, fertilizantes o riego en grandes parcelas. Mientras que “grande” era literalmente grande en el ejemplo anterior, este no es siempre el caso, pues no siempre tenemos problemas de contexto espacial. Consideremos un ejemplo en el que una asperjadora funciona con diferentes boquillas y utiliza soluciones de densidades. Mientras que es fácil cambiar las densidades es mucho más tedioso cambiar los ajustes de las boquillas, por lo tanto, no queremos cambiarlas con demasiada frecuencia. Podríamos pensar en un diseño experimental donde cambiáramos la configuración de la asperjadora y siguiéramos usando la misma configuración para diferentes soluciones y sus densidades o concentraciones. Esto significa que no estamos aleatorizando completamente la concentración, sino que lo hacemos en cada boquilla. Este sería otro ejemplo de un diseño de parcela dividida en el que los ajustes de la boquilla son el factor de parcela principal y el material o solución es el factor de parcela dividida. Utilizando esta terminología, el factor que es difícil de cambiar será el factor de la parcela principal.

El precio que pagamos por esta “pereza” en el nivel de la parcela completa es menos precisión, o menos potencia, para el efecto principal correspondiente porque tenemos muchas menos observaciones en este nivel. No aplicamos el tratamiento de la parcela principal con mucha frecuencia; por lo tanto, no podemos esperar que nuestros resultados sean muy precisos respecto al factor de la subparcela. Esto también puede observarse en la tabla ANOVA donde los grados de libertad del denominador para cada factor son diferentes. Nótese que el efecto principal del factor de la parcela dividida y la interacción entre el factor de la parcela dividida y el factor de la parcela principal no se ven afectados por esta pérdida de eficacia. De hecho, en el nivel de parcela dividida, estamos haciendo un experimento eficiente mientras bloqueamos en las parcelas principales.

Cuando se planifica un experimento: un pensamientos como, “Es más fácil si no cambiamos estos ajustes con demasiada frecuencia…”, puede ayudar a orientar la elección del factor a la parcela principal. Si no tenemos en cuenta la estructura especial de las parcelas divididas, los resultados a nivel de la parcela principal serán normalmente demasiado optimistas, lo que significa que los valores p son demasiado pequeños. De nuevo, no hay almuerzo gratis, este es el precio que pagamos por la “pereza.”.