Diseño en bloques completos, completamente aleatorizados
Kenneth Roy Cabrera Torres
Jueves, 15 de marzo de 2018
Fuentes de nitrógeno para producción.
Un agrónomo desea determinar el efecto de diferentes fuentes de nitrógeno en la producción de una materia seca sobre cebada forrajera. Hay cinco fuentes a ser comparadas: \((NH_4)_2 SO_4\) , \(NH_4NO_3\), \(CO(NH_2)_2\), \(Ca(NO_3)_2\), \(NaNO_3\) y un tratamiento control sin nitrógeno. Se desea aplicar los resultados sobre un rango bastante amplio de condiciones, se hicieron ensayos sobre cuatro tipos de suelo.
Para el diseño experimental se eligió un diseño en bloques completamente aleatorizado con los tipos de suelo como factor de bloqueo, se localizaron seis parcelas en cada uno de los cuatro tipos de suelo, y se asignó aleatoriamente los tratamientos a las parcelas dentro de cada tipo de suelo. La variable de interés es la producción en (kg/parcela) de cebada bajo varias fuentes de nitrógeno.
Los datos obtenidos de realizar este experimento se presentan en la tabla 7.3.
library(readxl)
ferti1 <- read_excel("fertilidad.xlsx")
ferti1$tratamiento <- factor(ferti1$tratamiento,
levels = c("Control", "CaNO32", "CONH22",
"NaNO3", "NH42SO4", "NH4NO3"))
ferti1$nombre <- factor(ferti1$nombre,
levels = c( "Sin nitrógeno", "Nitrato de calcio",
"Urea", "Nitrato de sodio",
"Sulfato de amonio", "Nitrato de amonio"))
ferti1$tipo_suelo <- factor(ferti1$tipo_suelo,
levels = c("I", "II", "III", "IV"))Gráfica exploratoria
library(ggplot2)
g1 <- ggplot(ferti1, aes(nombre, produccion, col = tipo_suelo))
g1 + geom_point() +
geom_line(aes(group = tipo_suelo)) +
theme(axis.text.x = element_text(angle = 30, hjust = 1))
Modelo.
Se propone:
\[ y_{ij} = \mu + \alpha_i + \beta_j + \epsilon_{ij}, \qquad i = 1, 2, \cdots, 6, \quad j = 1, 2, \cdots, 4 \\ \epsilon_{ij} \sim \mathcal{N}(0,\sigma^2) \quad i.i.d. \] ## Hipótesis.
Como hipótesis principal se propone:
\[ H_0: \alpha_1 = \alpha_2 = \alpha_3 = \alpha_4 = \alpha_5 = \alpha_6 = 0 \\ H_A: \textrm{Alguna } \alpha_i \textrm{ diferente de } 0 \] También se tiene una hipótesis de verificiación en cuanto al bloque:
\[ H_0: \beta_1 = \beta_2 = \beta_3 = \beta_4 = 0 \\ H_A: \textrm{Alguna } \alpha_i \textrm{ diferente de } 0 \]
Ajuste del modelo.
modelo1 <- aov(produccion ~ nombre + tipo_suelo,
data = ferti1)
summary(modelo1)## Df Sum Sq Mean Sq F value Pr(>F)
## nombre 5 256.15 51.23 16.85 1.10e-05 ***
## tipo_suelo 3 192.75 64.25 21.13 1.22e-05 ***
## Residuals 15 45.62 3.04
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Estimación de los parámetros del modelo
Efectos
model.tables(modelo1, type = "effects", se = TRUE)## Tables of effects
##
## nombre
## nombre
## Sin nitrógeno Nitrato de calcio Urea Nitrato de sodio Sulfato de amonio
## -5.492 0.183 -1.492 -0.142 5.408
## Nitrato de amonio
## 1.533
##
## tipo_suelo
## tipo_suelo
## I II III IV
## -4.008 -0.342 3.975 0.375
##
## Standard errors of effects
## nombre tipo_suelo
## 0.8719 0.7119
## replic. 4 6Promedios
model.tables(modelo1, type = "means", se = TRUE)## Tables of means
## Grand mean
##
## 30.84167
##
## nombre
## nombre
## Sin nitrógeno Nitrato de calcio Urea Nitrato de sodio Sulfato de amonio
## 25.35 31.03 29.35 30.70 36.25
## Nitrato de amonio
## 32.38
##
## tipo_suelo
## tipo_suelo
## I II III IV
## 26.83 30.50 34.82 31.22
##
## Standard errors for differences of means
## nombre tipo_suelo
## 1.233 1.007
## replic. 4 6Estimación del error aleatorio \(\sigma^2\).
# Estimador de la varianza del error aleatorio
varErr <- summary(modelo1)[[1]]$"Mean Sq"[3]
# Estimación de la desviación estándar del error aleatorio
desvErr <- sqrt(varErr)La estimación de la varianza de los errores aleatorios es de: 3.04. Es decir estimación de la desviación estándar de los errores es de: 1.74.
Diagnósticos del modelo.
residuales <- residuals(modelo1)
residualesEstandarizados <- rstandard(modelo1) # Residuales estudentizados internos
residualesEstudentizado <- rstudent(modelo1) # Residuales estudentizados externosSupuesto de normalidad.
shapiro.test(residuales)##
## Shapiro-Wilk normality test
##
## data: residuales
## W = 0.96731, p-value = 0.6013shapiro.test(residualesEstandarizados)##
## Shapiro-Wilk normality test
##
## data: residualesEstandarizados
## W = 0.96731, p-value = 0.6013shapiro.test(residualesEstudentizado)##
## Shapiro-Wilk normality test
##
## data: residualesEstudentizado
## W = 0.96563, p-value = 0.5613Supesto de igualdad de varianzas.
bartlett.test(produccion ~ nombre,
data = ferti1)##
## Bartlett test of homogeneity of variances
##
## data: produccion by nombre
## Bartlett's K-squared = 2.8214, df = 5, p-value = 0.7275bartlett.test(produccion ~ tipo_suelo,
data = ferti1)##
## Bartlett test of homogeneity of variances
##
## data: produccion by tipo_suelo
## Bartlett's K-squared = 1.6051, df = 3, p-value = 0.6582Prueba de independiencia.
No se tiene el orden o la secuencia de experimentación, por lo tanto no se puede verificar el supuesto de independencia.
Supuestos de linealidad y supuesto de aditividad.
par.orig <- par(no.readonly = TRUE)
par(mfrow = c(2, 2))
plot(modelo1)
par(par.orig)ferti1$residuales <- residuales
ferti1$stdresiduales <- residualesEstandarizados
ferti1$studresiduales <- residualesEstudentizado
ferti1$predichos <- predict(modelo1)
library(ggplot2)
g1 <- ggplot(ferti1, aes(produccion, residuales))
g1 + geom_point() + geom_smooth(method = lm) +
geom_hline(yintercept = 0, col = "red")
library(ggplot2)
g1 <- ggplot(ferti1, aes(predichos, residuales))
g1 + geom_point() + geom_smooth(method = lm) +
geom_hline(yintercept = 0, col = "red")
Modelo final
predichosIC <- data.frame(predict(modelo1, interval = "confidence"))
ferti1$inf <- predichosIC$lwr
ferti1$sup <- predichosIC$upr
library(ggplot2)
g1 <- ggplot(ferti1, aes(nombre, produccion, col = tipo_suelo))
g1 + geom_point() +
geom_point(aes(nombre, predichos), shape = 4, size = 2) +
geom_errorbar(aes(ymax = sup, ymin = inf), width = 0.2) +
ylab("Producción (Kg/parcela)") +
xlab("Fertilizante") +
labs(title = "Producción vs fertilizante",
subtitle = "Tipo de suelo como bloque") +
geom_line(aes(nombre, predichos, group = tipo_suelo)) +
theme(axis.text.x = element_text(angle = 30, hjust = 1)) +
guides(col = guide_legend(title = "Tipo de suelo"))
library(multcomp)## Loading required package: mvtnorm## Loading required package: survival## Loading required package: TH.data## Loading required package: MASS##
## Attaching package: 'TH.data'## The following object is masked from 'package:MASS':
##
## geyserPruebaDunnet <- glht(modelo1, linfct = mcp(nombre = "Dunnett"))
summary(PruebaDunnet)##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Dunnett Contrasts
##
##
## Fit: aov(formula = produccion ~ nombre + tipo_suelo, data = ferti1)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## Nitrato de calcio - Sin nitrógeno == 0 5.675 1.233 4.602 0.00141 **
## Urea - Sin nitrógeno == 0 4.000 1.233 3.244 0.02195 *
## Nitrato de sodio - Sin nitrógeno == 0 5.350 1.233 4.339 0.00264 **
## Sulfato de amonio - Sin nitrógeno == 0 10.900 1.233 8.839 < 0.001 ***
## Nitrato de amonio - Sin nitrógeno == 0 7.025 1.233 5.697 < 0.001 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)(pruebaTukey <- TukeyHSD(modelo1))## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = produccion ~ nombre + tipo_suelo, data = ferti1)
##
## $nombre
## diff lwr upr p adj
## Nitrato de calcio-Sin nitrógeno 5.675 1.66867109 9.6813289 0.0038124
## Urea-Sin nitrógeno 4.000 -0.00632891 8.0063289 0.0504763
## Nitrato de sodio-Sin nitrógeno 5.350 1.34367109 9.3563289 0.0063062
## Sulfato de amonio-Sin nitrógeno 10.900 6.89367109 14.9063289 0.0000031
## Nitrato de amonio-Sin nitrógeno 7.025 3.01867109 11.0313289 0.0004979
## Urea-Nitrato de calcio -1.675 -5.68132891 2.3313289 0.7496184
## Nitrato de sodio-Nitrato de calcio -0.325 -4.33132891 3.6813289 0.9997870
## Sulfato de amonio-Nitrato de calcio 5.225 1.21867109 9.2313289 0.0076579
## Nitrato de amonio-Nitrato de calcio 1.350 -2.65632891 5.3563289 0.8760096
## Nitrato de sodio-Urea 1.350 -2.65632891 5.3563289 0.8760096
## Sulfato de amonio-Urea 6.900 2.89367109 10.9063289 0.0005982
## Nitrato de amonio-Urea 3.025 -0.98132891 7.0313289 0.1995232
## Sulfato de amonio-Nitrato de sodio 5.550 1.54367109 9.5563289 0.0046250
## Nitrato de amonio-Nitrato de sodio 1.675 -2.33132891 5.6813289 0.7496184
## Nitrato de amonio-Sulfato de amonio -3.875 -7.88132891 0.1313289 0.0608085
##
## $tipo_suelo
## diff lwr upr p adj
## II-I 3.6666667 0.7648371 6.5684962 0.0115167
## III-I 7.9833333 5.0815038 10.8851629 0.0000052
## IV-I 4.3833333 1.4815038 7.2851629 0.0028416
## III-II 4.3166667 1.4148371 7.2184962 0.0032348
## IV-II 0.7166667 -2.1851629 3.6184962 0.8908645
## IV-III -3.6000000 -6.5018295 -0.6981705 0.0131166Consecuencias de no hacer el bloqueo.
modelo2 <- aov(produccion ~ nombre, data = ferti1)
summary(modelo2)## Df Sum Sq Mean Sq F value Pr(>F)
## nombre 5 256.1 51.23 3.869 0.0148 *
## Residuals 18 238.4 13.24
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Estimador de la varianza del error aleatorio
(varErr2 <- summary(modelo2)[[1]]$"Mean Sq"[2])## [1] 13.2425# Estimación de la desviación estándar del error aleatorio
(desvErr2 <- sqrt(varErr2))## [1] 3.639025PruebaDunnet2 <- glht(modelo2, linfct = mcp(nombre = "Dunnett"))
summary(PruebaDunnet2)##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Dunnett Contrasts
##
##
## Fit: aov(formula = produccion ~ nombre, data = ferti1)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## Nitrato de calcio - Sin nitrógeno == 0 5.675 2.573 2.205 0.14468
## Urea - Sin nitrógeno == 0 4.000 2.573 1.554 0.41453
## Nitrato de sodio - Sin nitrógeno == 0 5.350 2.573 2.079 0.18146
## Sulfato de amonio - Sin nitrógeno == 0 10.900 2.573 4.236 0.00213 **
## Nitrato de amonio - Sin nitrógeno == 0 7.025 2.573 2.730 0.05327 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)(pruebaTukey <- TukeyHSD(modelo2))## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = produccion ~ nombre, data = ferti1)
##
## $nombre
## diff lwr upr p adj
## Nitrato de calcio-Sin nitrógeno 5.675 -2.502653 13.852653 0.2826336
## Urea-Sin nitrógeno 4.000 -4.177653 12.177653 0.6361203
## Nitrato de sodio-Sin nitrógeno 5.350 -2.827653 13.527653 0.3402462
## Sulfato de amonio-Sin nitrógeno 10.900 2.722347 19.077653 0.0056028
## Nitrato de amonio-Sin nitrógeno 7.025 -1.152653 15.202653 0.1172482
## Urea-Nitrato de calcio -1.675 -9.852653 6.502653 0.9851439
## Nitrato de sodio-Nitrato de calcio -0.325 -8.502653 7.852653 0.9999946
## Sulfato de amonio-Nitrato de calcio 5.225 -2.952653 13.402653 0.3642271
## Nitrato de amonio-Nitrato de calcio 1.350 -6.827653 9.527653 0.9944171
## Nitrato de sodio-Urea 1.350 -6.827653 9.527653 0.9944171
## Sulfato de amonio-Urea 6.900 -1.277653 15.077653 0.1279938
## Nitrato de amonio-Urea 3.025 -5.152653 11.202653 0.8425115
## Sulfato de amonio-Nitrato de sodio 5.550 -2.627653 13.727653 0.3039622
## Nitrato de amonio-Nitrato de sodio 1.675 -6.502653 9.852653 0.9851439
## Nitrato de amonio-Sulfato de amonio -3.875 -12.052653 4.302653 0.6650810