Parcelas Divididas

Ejemplo

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

Para acceder a los datos de este ejemplo vamos al repositorio del docente donde reposan varios conjuntos de datos.

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

Sobre Repeticiones

Es importante entender las repeticiones que se disponen por tratamiento y parcelas, para eso generamos algunas tablas cruzadas.

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.

En lo que respecta a las repeticiones por parcela, generamos una nueva tabla:

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

Interacción

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.

Ajute del modelo

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

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.”.

Cuadrados latinos

Definición

Cuadrado Latino

Definición: Un diseño cuadrado latino es un diseño en el que se involucran tres factores (normalmente un tratamiento y dos de bloqueo), cada uno de los cuales es de naturaleza categórica y poseen el mismo número de niveles, de modo que cada par de variables factor tienen el mismo número de observaciones para cada posible par de niveles. Los tratamientos se colocan en un cuadrado de modo que cada fila y columna contiene cada tratamiento una vez.

Un diseño de cuadrado latino es un ejemplo de un diseño de bloques incompletos cuando hay un solo tratamiento y dos variables de bloqueo, cada una con el mismo número de niveles. Se aplica un solo tratamiento dentro de cada combinación de variables de bloqueo.

El término “cuadrado latino” es tomado directamente de las matemáticas. Antes de pasar a la discusión del diseño cuadrado latino, vamos a describir lo que significa un “cuadrado latino” en matemáticas. De acuerdo con la matemáticos, “Un cuadrado latino (puede abreviarse como”LS”) es una tabla cuadrada \((n\times n)\) llena de \(n\) símbolos diferentes (por ejemplo, alfabetos griegos α, β, γ, δ, etc.), de tal manera que cada símbolo aparece exactamente una vez en cada fila y exactamente una vez en cada columna.” En este sentido, el estudio de los cuadrados latinos es una cuestión de matemáticas en general y de combinatoria en particular. Aquí presentamos una breve descripción de Latin Square solo en referencia a Latin Square Design.

Con cualquier número de letras, podemos ver que se puede formar más de un número de cuadrados latinos. Por ejemplo, con 4 letras griegas α, β, γ, δ podemos construir los siguientes cuadrados latinos:

donde R1, R2, R3, R4 son los números de fila y C1, C2, C3, C4 los números de columna. En todos los tres cuadrados latinos anteriores, se puede ver que las tablas son cuadradas, cada uno con cuatro filas y cuatro columnas, de modo que cada letra griega se produce una vez en fila y una vez en una columna. Sin embargo, los tres cuadrados latinos anteriores no son los únicos cuadrados que podemos construir a partir de 4 letras. Es una cuestión de Combinatoria y Matemáticas para listar todos los cuadrados posibles con números \(p\) símbolos. Cuadrados estándar: la totalidad de los cuadrados latinos obtenidos de un solo cuadrado latino permutando las filas, columnas y las letras se llama un “conjunto transformado”. Un cuadrado latino \(p\times p\) con \(p\) letras en el orden natural que se produce en la primera fila y en la primera columna se llama _“cuadrado latino estándar_”, así, LS-1 es un cuadrado estándar mientras que LS-2 y LS3 son Conjuntos Transformados de LS-1.

Se puede demostrar que el número de cuadrados que pueden generarse a partir de un cuadrado estándar de orden \(p\times p\), mediante la permutación de filas, columnas y letras es \(p!\) . Todos estos cuadrados no son necesariamente diferentes. También, obtenemos \(p!(p-1)!\) cuadrados latinos diferentes permutando todas las \(p\) columnas y las \((p-1)\) filas, excepto la primera fila. Así, con p = 3, podemos generar 12 cuadrados Latinos; con p = 4, obtenemos 144 y así sucesivamente.

Uso

  • Los diseños de cuadrados latinos son opciones razonables cuando es imposible utilizar cada nivel de tratamiento para la misma combinación de niveles de bloqueo. Por ejemplo, considere un experimento con cuatro dietas, cada una de las cuales debe darse a cuatro vacas de forma consecutivas. Si cada vaca recibió las dietas en el mismo orden, el efecto del tratamiento se confundiría con el efecto debido al orden en que se dieron las dietas. Cada vaca puede recibir una sola dieta durante un único período de tiempo. Es común utilizar múltiples cuadrados latinos en un solo experimento.

Creando Cuadrados en R

latinos= function (n, nrand = 20){
x = matrix(LETTERS[1:n], n, n)
x = t(x)
for (i in 2:n) x[i, ] = x[i, c(i:n, 1:(i - 1))]
if (nrand > 0) {
for (i in 1:nrand) {
x = x[sample(n), ]
x = x[, sample(n)]
 }
}
x
}

A continuación ejecutamos la función para crear un cuadrado lationo \(5\times 5\)

latinos(5)
##      [,1] [,2] [,3] [,4] [,5]
## [1,] "C"  "A"  "E"  "B"  "D" 
## [2,] "E"  "C"  "B"  "D"  "A" 
## [3,] "A"  "D"  "C"  "E"  "B" 
## [4,] "B"  "E"  "D"  "A"  "C" 
## [5,] "D"  "B"  "A"  "C"  "E"

Ejemplo

A continuación se ilustra el diseño considerando los datos relacionados con el peso de las raíces de la acelga. Utilizaron un diseño cuadrado latino de orden 5. El campo rectangular en el que se realizó el experimento se dividió en cinco filas y cinco columnas. Se crearon 25 parcelas, dispuestas en un cuadrado, sobre el cual se aplicaron los tratamientos asociados a caudal de riego (A: sin riego, B: riego limitado fuertemente, C: riego limitado moderadamente, D: riego poco limitado y E: riego usual). Cada fila del cuadrado fue considerada como un bloque, por lo que cada tratamiento se aplicó en cada fila. La característica única de los diseños cuadrados latinos es que hay un segundo conjunto de bloques. Cada columna se consideraba también un bloque, por lo que cada tratamiento se aplicaba también a todas las columnas. Los datos se dan en el cuadro que se muestra a continuación, ordenados por filas y columnas con el tratamiento dado en el lugar apropiado y el peso radicular observado entre paréntesis.

Las razones de bloqueo se asociaron a la pendiente y a la proximidad a una zona con alta pedregosidad.

Modelo

El modelo para un cuadrado latino de orden \(p\) es el del análisis de varianza de tres vías:

\[y_{ijk}=\mu+\rho_i+\gamma_j+ \tau_{k}+\epsilon_{ijk} \]

donde \(\mu\) es la media poblacional, \(\rho_i\) representa al efectó de la i-ésima fila,\(\gamma_j\) es ele efcto de la j-ésima columna, \(\tau_{k}\) es el efecto del k-ésimo tratamiento ( i =j = k = 1,p) . Además, se tienen las condicones laterales de que \(\sum_{i=1}^p{\rho_i}=0\),\(\sum_{j=1}^p{\gamma_j}=0\) y \(\sum_{k=1}^p{\tau_k}=0\) . Además, se tiene el supuesto de que \(\epsilon_{ijk}\sim {N(0, \sigma^2)}\), con independencia e idéntica distribución.

Datos en R

peso=c(376,371,355,356,335,316,338,336,356,332,326,326,335,343,330,317,343,330,327,336,321,332,317,318,306)
fila=rep(c("R5","R4","R3","R2","R1"),each=5)
columna=rep(c("C1","C2","C3","C4","C5"),5)
tratamiento=(c(matrix(c("D","E","C","B","A","B","D","E","A","C","C","A","B","D","E","E","B","A","C","D","A","C","D","E","B"),byrow = T)))
dfcl=data.frame(peso,fila, columna, tratamiento)
datatable(dfcl)

Descripción

Para estudiar el comportamiento de la respuesta en los bloques veremos como lucen los siguientes diagramas:

ggplot(dfcl,aes(x = columna,y = fila,fill = peso)) + 
  geom_tile(color = "white",
            lwd = 1.5,
            linetype = 1) + 
  geom_text(aes(label = tratamiento), color = "white", size = 4)+
  scale_fill_gradient(low = "green", high = "red")+
  scale_y_discrete(labels = c("R5","R4","R3","R2","R1"))

Con la imagen es más claro donde se concentran los mayores rendimientos y qué tratamientos se asocian. Imagine que de R1 a R5 la pendiente aumenta que de C1 a C5 la proximidad a la región de alta pedregosidad disminuye. Con esta información es claro que los mayores rendimientos se observan en los tratamientos E y D que son las mejores condiciones de riego, pero en las menores pendientes en las condiciones de menor pedregosidad.

Análisis de Varianza

acelga= lm(peso ~ fila + columna + tratamiento)
anova(acelga)
## Analysis of Variance Table
## 
## Response: peso
##             Df Sum Sq Mean Sq F value   Pr(>F)   
## fila         4 4240.2 1060.06  7.2511 0.003294 **
## columna      4  701.8  175.46  1.2002 0.360412   
## tratamiento  4  330.2   82.56  0.5647 0.692978   
## Residuals   12 1754.3  146.19                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Claramente los datos proveen evidencia a favor de la hipótesis nula de efecto nulo de los tratamientos del riego. Al menos en estos valores probados no se ve afectado claramente el rendimiento. Los valores del estadístico F claramente muestran las ventajas que tuvo el bloqueo, siendo más eficiente el bloque por filas que el de columnas.

Cuadrados de Youden

Surgimiento

Cuadrados de Youden

No siempre los bloques se pueden construir de modo que las dos razones de bloqueo tengan los mismos niveles. Si se trata de bloqueo por aspectos espaciales, es posible generar los niveles de modo que se conforme facilmente un cuadrado, por ejemplo, en pedregodidad, las distancias se pueden generar en 5 niveles como los tratamientos particionando las distancias. Lo mismo puede hacerse con el riego, sin embargo , si la razón de blqoueo se generara por el tipo de bandejas con las que se hace el ensayo en la etapa de germinación, es probable que se disponga de menos tipos de materiales que de niveles de riego. Otro ejemplo puede asociarse a la procedencia del material o del fertilizantes, donde podemos tener tres casos comerciales y cinco niveles de riego. Esta estructura ya no permite generar un cuadrado, por lo que allí entran los cuadrados de Youden.

Datos

Considere los datos de las raíces de acelgas del caso del diseño del cuadrado latino, pero asuma que hay cinco filas, cuatro columnas y cinco tratamientos. Si ignoramos las columnas, las filas y los tratamientos forman un diseño de bloques incompletos balanceado, en el que cada par de tratamientos ocurre juntos tres veces. La característica clave de los cuadrados de Youden es que, además, los tratamientos, también se configuran de tal manera que cada tratamiento ocurre una vez en cada columna. Dado que cada fila también se produce una vez en cada columna, el análisis de las columnas puede realizarse independientemente del análisis de las filas y los tratamientos. Las columnas están balanceadas en relación con los tratamientos y las filas.

El cuadro anterior contiene los datos sin la última columna del problema de las acelgas. Las filas deben ajustarse antes de los tratamientos. Mientras se mantenga el equilibrio, no importa dónde se montan las columnas. Si los datos se desequilibran, los tratamientos deben ser ajustados por último.

pesoy=c(376,371,355,356,316,338,336,356,326,326,335,343,317,343,330,327,321,332,317,318)
filay=rep(c("R5","R4","R3","R2","R1"),each=4)
columnay=rep(c("C1","C2","C3","C4"),5)
tratamientoy=(c(matrix(c("D","E","C","B","B","D","E","A","C","A","B","D","E","B","A","C","A","C","D","E"),byrow = T)))
dfcly=data.frame(pesoy,filay, columnay, tratamientoy)

ANOVA Youden

anova_youden <- aov(pesoy ~ filay + columnay + tratamientoy, data = dfcly)
summary(anova_youden)
##              Df Sum Sq Mean Sq F value Pr(>F)  
## filay         4   4247  1061.8   6.869 0.0106 *
## columnay      3    367   122.3   0.791 0.5320  
## tratamientoy  4    224    56.0   0.362 0.8289  
## Residuals     8   1237   154.6                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Los datos proveen evidencia a facor de la hipotesis nula de efecto nulo de los tratamientos.

Vericación del Cuadrado de Youden

Para verificar si tus datos forman un cuadrado de Youden, puedes comprobar si cada par de tratamientos aparece juntos el mismo número de veces. En un cuadrado de Youden, cada par de tratamientos debería aparecer juntos exactamente una vez en cada fila y columna.Aquí proporcionamos el código en R para verificar esto:

library(stringr)
## Warning: package 'stringr' was built under R version 4.2.3
library(stringi)
library(gtools)
# Datos proporcionados
pesoy <- c(376, 371, 355, 356, 316, 338, 336, 356, 326, 326, 335, 343, 317, 343, 330, 327, 321, 332, 317, 318)
filay <- rep(c("R5", "R4", "R3", "R2", "R1"), each = 4)
columnay <- rep(c("C1", "C2", "C3", "C4"), 5)
tratamientoy <- c("D", "E", "C", "B", "B", "D", "E", "A", "C", "A", "B", "D", "E", "B", "A", "C", "A", "C", "D", "E")
dfcly <- data.frame(pesoy, filay, columnay, tratamientoy)

matriz_tratamientos <- matrix(tratamientoy, nrow = 5, byrow = TRUE);matriz_tratamientos
##      [,1] [,2] [,3] [,4]
## [1,] "D"  "E"  "C"  "B" 
## [2,] "B"  "D"  "E"  "A" 
## [3,] "C"  "A"  "B"  "D" 
## [4,] "E"  "B"  "A"  "C" 
## [5,] "A"  "C"  "D"  "E"
co = list()

for(i in 1:dim(matriz_tratamientos)[1]){
  m = matriz_tratamientos[i,]
  c = combn(x = m,m = 2, simplify = F)
  c_l = lapply(c, paste, collapse = '-')
  co[[i]] = unlist(c_l)
  tco=c(unlist(co))
  ttco=table(tco)
  ordenar_cadena <- function(cadena) {
  partes <- unlist(strsplit(cadena, "-"))
  partes_ordenadas <- sort(partes)
  paste(partes_ordenadas, collapse = "-")
  }
  cadenas_ordenadas <- c(sapply(tco, ordenar_cadena))
}
table(cadenas_ordenadas)
## cadenas_ordenadas
## A-B A-C A-D A-E B-C B-D B-E C-D C-E D-E 
##   3   3   3   3   3   3   3   3   3   3

Como todos los pares de tratamientos aparecen una misma cantidad de veces, el diseño desarrollado en un cuadrado de Youden.

Diseño Anidado: Factorial Incompleto

Definición

Diseño Anidado en dos etapas

Un diseño anidado es un tipo de diseño experimental en el que los niveles de un factor están jerárquicamente anidados dentro de los niveles de otro factor. Por ejemplo, imaginemos que estamos interesados en el efecto de un tipo de medicamento sobre la expresión de un gen específico. Para ello diseñamos un experimento anidado con ratones como unidades experimentales (en la figura dos ratones por tratamiento), donde incluimos un control. De cada ratón se tomaron tres células y en cada una se evaluó la expresión del gen dos veces (repeticiones técnicas).

Como se muestra en la figura, nuestro diseño tiene una apariencia jerárquica. En este tipo de diseño distinguimos entre dos tipos de factores: fijos y aleatorios. Un factor fijo es uno que tiene valores discretos o finitos, mientras que un factor aleatorio puede tomar muchos valores. En nuestro ejemplo, el factor asociado al medicamento se consideraría un factor fijo y los factores asociados al ratón, la célula y de las mediciones repetidas se considerarían factores aleatorios pues son muestras aleatorias de una población mayor de ratones, células y medidas. Nótese que los factores aleatorios son similares, como el ratón y la célula, pero no idénticos entre sí y se anidan sucesivamente hasta llegar al factor fijo, es decir, las células se anidan en los ratones , son sus propias células y en esencia representan al ratón, además, las medidas son propias de las células selecionadas. Es claro que la estructura no es completa pues los niveles del cualquier factor posterior al ratón solo son inherentes a los ratones usados con el medicamente o en el control.

Ahora veamos cómo podemos analizar los resultados de este tipo de diseño experimental con la ayuda del código R.

Ejemplo

Habiendo ilustrado el diseño, ahora imaginemos un problema más amplio pero respetando la naturaleza del ejemplo ilustrado. En el marco de datos tenemos las réplicas técnicas (rept), de las que se hicieron tres en cada célula, El factor A indica los niveles del medicamento, donde se uso el control y otros dos medicamentos diferentes, el factor B los niveles del factor asociado al ratón, de los cuales se usaron cinco grupos y el factor C representó al factor célula con cinco niveles tambíen y la última columna contiene los datos de expresión génica (expg).

GitHub

Ahora importaremos el archivo del GitHub

url <- "https://github.com/darghan/unal/raw/refs/heads/main/rat.xlsx"
temp_file <- tempfile(fileext = ".xlsx")
download.file(url, temp_file, mode = "wb")
rat <- read_excel(temp_file)
colnames(rat)=c("rept","cel","rat","med", "expg")
rat$rept=as.factor(rat$rept)
rat$rept
##   [1] 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1
##  [38] 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2
##  [75] 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3
## [112] 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1
## [149] 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2
## [186] 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3
## [223] 1 2 3
## Levels: 1 2 3

Visualización

Convirtiendo variables en factores

rat2 <- rat %>%
  mutate(
    rept = as.factor(rept),
    cel = as.factor(cel),
    rat = as.factor(rat),
    med= as.factor(med)
  )

Generando los gráficos que se anidan

rat2$med=factor(rat2$med)
expr_plot <- ggplot(rat, aes(x = med, y = expg)) +
  geom_point(aes(color = rept), size = 1) +
  stat_summary(
    fun = mean, geom = "crossbar", width = 0.5, color = "darkblue", linewidth = 0.3
  ) +
  facet_grid(rat ~ cel, labeller = label_both)+
  scale_x_continuous(breaks = seq(min(as.numeric(rat2$med)), max(as.numeric(rat2$med)), by = 1))

expr_plot

ANOVA

Organizando los factores segun su naturaleza

rat2$med <- as.fixed(rat2$med)
rat2$rat <- as.random(rat2$rat)
rat2$cel <- as.random(rat2$cel)

Para ajustar el modelo lineal debemos tener en cuenta la relación entre nuestra respuesta y los factores especificados anteriormente:

data_aov <- aov(expg ~ med + rat:med + cel:rat:med, data = rat2)

El término rat:med denota la variabilidad de los ratones dentro de cada tratamiento, y el término cel:rat:med denota la variabilidad de las células dentro de cada ratón y a su vez dentro de cada tratamiento.

Para mostrar la tabla ANOVA utilizamos la función gad():

gad_result <- gad(data_aov)
print(gad_result)
## $anova
## Analysis of Variance Table
## 
## Response: expg
##              Df Sum Sq Mean Sq F value    Pr(>F)    
## med           2 743.44  371.72  12.796  0.001058 ** 
## med:rat      12 348.59   29.05   5.163 7.413e-06 ***
## med:rat:cel  60 337.59    5.63  11.153 < 2.2e-16 ***
## Residuals   150  75.67    0.50                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

La función gad() distingue entre efectos fijos y aleatorios, así como la estructura anidada entre estos factores, por lo que realiza correcciones para calcular las estadísticas F. La razón de esto es que, dependiendo de si estamos considerando efectos fijos o aleatorios, los valores esperados para los cuadrados medios de error (EM) cambian.

Estimando Variabilidad

rat3 <- rat2 %>% 
  mutate(med = as.factor(med), rat = as.factor(rat), cel = as.factor(cel))
data_lme <- lmer(expg ~ 1 + med + (1|rat:med) + (1|cel:rat:med), data = rat3)

Tenga en cuenta que el ajuste del modelo lineal con la función lmer() requiere una sintaxis algo diferente a la utilizada con aov() y gad(). Para mostrar la contribución a la variabilidad de cada factor utilizamos la función summary() con el objeto data_lme como argumento:

summary(data_lme)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: expg ~ 1 + med + (1 | rat:med) + (1 | cel:rat:med)
##    Data: rat3
## 
## REML criterion at convergence: 684.4
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.76545 -0.53144 -0.02996  0.59495  2.20017 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  cel:rat:med (Intercept) 1.7073   1.3066  
##  rat:med     (Intercept) 1.5615   1.2496  
##  Residual                0.5045   0.7103  
## Number of obs: 225, groups:  cel:rat:med, 75; rat:med, 15
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)  10.0519     0.6224 12.0000  16.151 1.66e-09 ***
## med2          1.4654     0.8801 12.0000   1.665  0.12180    
## med3         -2.9085     0.8801 12.0000  -3.305  0.00629 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##      (Intr) med2  
## med2 -0.707       
## med3 -0.707  0.500

En efectos aleatorios, el término residual se refiere a réplicas técnicas, el término rat:med a ratones y el término cel:rat:med a células.

Comparaciones Múltiples

mult_comp <- glht(data_lme, linfct = mcp(med = "Tukey"))
summary(mult_comp )
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: lmer(formula = expg ~ 1 + med + (1 | rat:med) + (1 | cel:rat:med), 
##     data = rat3)
## 
## Linear Hypotheses:
##            Estimate Std. Error z value Pr(>|z|)    
## 2 - 1 == 0   1.4654     0.8801   1.665  0.21884    
## 3 - 1 == 0  -2.9085     0.8801  -3.305  0.00269 ** 
## 3 - 2 == 0  -4.3739     0.8801  -4.970  < 1e-04 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

Notas

  • Un diseño anidado es un tipo de diseño experimental en el que los niveles de un factor están jerárquicamente anidados dentro de los niveles de otro factor.

*En este tipo de diseño se distingue entre factores fijos (que tienen valores discretos o finitos) y factores aleatorios (que pueden tener muchos valores). En el ejemplo, se consideró que el factor del medicamento era fijo, mientras que los factores de las mediciones del ratón, la célula y las repeticiones fueron aleatorios.

  • Es posible utilizar el paquete ggplot2 para visualizar los resultados de este tipo de diseño.

  • Se utilizó el conjunto de datos GAD para realizar el análisis de la varianza de los resultados de este diseño.

*El paquete lme4 se utilizó para estimar la contribución a la variabilidad de los factores aleatorios en ratones, células y réplicas técnicas.