supuestos - ANOVA con diseño de bloque y medidas repetidas

supuestos de anova en r

Question

Estoy intentando ejecutar algunos análisis estadísticos en una prueba de campo que se construyó en 2 sitios durante la misma temporada de crecimiento.

En ambos sitios ( Site , niveles: HF | NW) el diseño experimental fue un RCBD con 4 (n = 4) bloques ( Block , niveles: 1 | 2 | 3 | 4 dentro de cada Site ). Hubo 4 tratamientos - 3 formas diferentes de fertilizante nitrogenado y un control (sin fertilizante nitrogenado) ( Treatment , niveles: AN, U, IU, C). Durante la prueba de campo hubo 3 períodos distintos que comenzaron con la adición de fertilizantes y terminaron con la cosecha de la hierba. Estos períodos han recibido los niveles 1 | 2 | 3 bajo el factor N_app .

Hay un rango de medidas que me gustaría probar la siguiente hipótesis nula H0 en:

Treatment (H0) no tuvo efecto en la medición

Dos de las medidas en las que estoy particularmente interesado son: el rendimiento de la hierba y las emisiones de amoníaco.

Comenzando con el rendimiento de pasto ( Dry_tonnes_ha ) como se muestra aquí, un conjunto de datos equilibrado agradable

Los datos se pueden descargar en R usando el siguiente código:

library(tidyverse)
download.file(''https://www.dropbox.com/s/w5ramntwdgpn0e3/HF_NW_grass_yield_data.csv?raw=1'', destfile = "HF_NW_grass_yield_data.csv", method = "auto")
raw_data <- read.csv("HF_NW_grass_yield_data.csv", stringsAsFactors = FALSE)
HF_NW_grass <- raw_data %>% mutate_at(vars(Site, N_app, Block, Plot, Treatment), as.factor) %>%
mutate(Date = as.Date(Date, format = "%d/%m/%Y"),
Treatment = factor(Treatment, levels = c("AN", "U", "IU", "C")))

He intentado ejecutar un ANOVA con el siguiente enfoque:

model_1 <- aov(formula = Dry_tonnes_ha ~ Treatment * N_app + Site/Block, data = HF_NW_grass, projections = TRUE)

Tengo algunas preocupaciones con esto.

En primer lugar, ¿cuál es la mejor manera de probar las suposiciones? Para un ANOVA simple de una vía, usaría shapiro.test() y bartlett.test() en la variable dependiente ( Dry_tonnes_ha ) para evaluar la normalidad y la heterogeneidad de la varianza. ¿Puedo usar el mismo enfoque aquí?

En segundo lugar, me preocupa que N_app sea ​​una medida repetida, ya que la misma medición se toma de la misma gráfica en 3 períodos diferentes: ¿cuál es la mejor manera de incorporar estas medidas repetidas en el modelo?

En tercer lugar, no estoy seguro de la mejor manera de anidar Block dentro del Site . En ambos sitios, los niveles de Block son 1: 4. ¿Necesito tener niveles de Block únicos para cada sitio?

Tengo otro conjunto de datos para las emisiones de NH3 aquí . Código R para descargar:

download.file(''https://www.dropbox.com/s/0ax16x95m2z3fb5/HF_NW_NH3_emissions.csv?raw=1'', destfile = "HF_NW_NH3_emissions.csv", method = "auto")
raw_data_1 <- read.csv("HF_NW_NH3_emissions.csv", stringsAsFactors = FALSE)
HF_NW_NH3 <- raw_data_1 %>% mutate_at(vars(Site, N_app, Block, Plot, Treatment), as.factor) %>%
mutate(Treatment = factor(Treatment, levels = c("AN", "U", "IU", "C")))

Para esto tengo todas las preocupaciones anteriores con la adición de que el conjunto de datos está desequilibrado. En HF para N_app 1 n = 3, pero para N_app 2 y 3 n = 4 En NW n = 4 para todos los niveles de N_app . En NF mediciones solo se realizaron en los niveles de Treatment U e IU En NW mediciones se realizaron en los niveles de Treatment AN , U e IU

No estoy seguro de cómo lidiar con este nivel adicional de complejidad. Estoy tentado de analizar solo como 2 sitios separados (el hecho de que los períodos de N_app no sean los mismos en cada sitio puede alentar este enfoque). ¿Puedo usar un ANOVA de suma de cuadrados tipo iii aquí?

Se me ha sugerido que un enfoque lineal de modelado mixto puede ser el camino a seguir, pero no estoy familiarizado con su uso.

Me gustaría recibir sus opiniones sobre cualquiera de los anteriores. Gracias por tu tiempo.

Rory

Answer 1

Estoy de acuerdo con @David C en el uso de diagnósticos visuales. Parcelas QQ simples deberían funcionar

# dependent variable.
par(mfrow=c(1,2))
qqnorm(dt[,dry_tonnes_ha]); qqline(dt[,dry_tonnes_ha], probs= c(0.15, 0.85))
qqnorm(log(dt[,dry_tonnes_ha])); qqline(log(dt[,dry_tonnes_ha]), probs= c(0.15, 0.85))

La transformación de registro parece razonable para mí. También puedes ver esto desde el diagrama de densidad, que es de larga cola y algo bimodal

par(mfrow=c(1,1))
plot(density(dt[,dry_tonnes_ha]))

Alternativamente, podría usar gráficos de alineación (Buja et al, 2009) si lo desea. No estoy seguro de que sean necesarios en este caso. Viñeta proporcionada

library(nullabor)
# this may not be the best X variable. I''m not familiar with your data
dt_l <- lineup(null_permute("dry_tonnes_ha"), dt)
qplot(dry_tonnes_ha, treatment, data = dt_l) + facet_wrap(~ .sample)

Para las otras suposiciones, puede usar los diagramas de diagnóstico estándar de la lm

lm2 <- lm(log(dry_tonnes_ha) ~ treatment * n_app + site/block, data = dt)
plot(lm2)

No veo nada demasiado problemático en estas tramas.

Answer 2

Para responder a su primera pregunta sobre la mejor manera de probar las suposiciones. Si bien su intento de utilizar otra prueba estadística, implementada en R, es razonable, realmente solo visualizaría la distribución y vería si los datos cumplen con los supuestos de ANOVA. Este enfoque puede parecer algo subjetivo, pero funciona en la mayoría de los casos.

• datos distribuidos de manera idéntica e independiente (iid): esta es una pregunta que ya puede tener una respuesta basada en cuánto sabe sobre sus datos. Es posible usar una prueba de chi-cuadrado para determinar la independencia (o no).
• datos distribuidos normalmente: utilice un gráfico de histograma / QQ para verificar. Basado en la distribución, creo que es razonable usar aov pesar de la distribución ligeramente bimodal.

(Parece que la transformación de log ayuda a cumplir con la suposición de normalidad. Esto es algo que puede considerar, especialmente para análisis posteriores).

par(mfrow=c(2,2))
plot(density(HF_NW_grass$Dry_tonnes_ha), col="red", main="Density")
qqnorm(HF_NW_grass$Dry_tonnes_ha, col="red", main="qqplot")
qqline(HF_NW_grass$Dry_tonnes_ha)
DTH_trans <- log10(HF_NW_grass$Dry_tonnes_ha)
plot(density(DTH_trans), col="blue", main="transformed density")
qqnorm(DTH_trans, col="blue", main="transformed density")
qqline(DTH_trans)

Con respecto a su segunda pregunta sobre cuál es la mejor manera de construir medidas repetidas en el modelo: Desafortunadamente, es difícil determinar cuál es el "mejor" modelo, pero basado en mi conocimiento (principalmente a través de big data de genómica), es posible que desee utilizar un modelo lineal de efectos mixtos. Esto se puede implementar a través del paquete lme4 R, por ejemplo. Como parece que ya sabes cómo construir un modelo lineal en R, no deberías tener problemas para aplicar funciones lme4 .

Su tercera pregunta sobre si anidar dos variables es complicada. Si yo fuera usted, comenzaría con Site y Block como si fueran factores independientes. Sin embargo, si sabes que no son independientes, probablemente deberías anidarlos.

Creo que sus preguntas y preocupaciones son bastante abiertas. Mi recomendación es que, mientras tengas una justificación plausible, adelante y proceda.