Randomized Block design
Datos
Descargar la base de datos aquÃ
Cargar la base de datos
# Install.pakages("readxl")
library(readxl)
datos <- read_excel("C:/Users/LVARGAS/Desktop/materialCi2016 inferencia/baseDatos/datos.xlsx", sheet = 4)
datos <- transform(datos, Trt = factor(Trt))
datos <- transform(datos, Bloque = factor(Bloque))
head(datos, 2)## Bloque Parc Trt Wdry50S Wgrain50S WidthArea LengthArea WGrain WHumid
## 1 1 5 1 96.8 38.14 1.5 9 2370 200.05
## 2 1 7 11 106.9 38.40 1.5 9 2670 200.88
## WDry W200g Area AreaCount Count Tasseling Lodging Flowering Maturity
## 1 180.43 6.99 13.5 0.187 26.7 56 0 60 107
## 2 182.06 7.19 13.5 0.187 25.8 53 0 56 106
## Height HI X.Humidity H2O_Grain Yield_Dry Yield_12_H2O Biomass
## 1 78.5 39.40083 0.09807548 232.4389 1583.379 1799.294 4018.643
## 2 74.7 35.92142 0.09368777 250.1464 1792.484 2036.914 4990.015
## Straw Spikes.m2 Thou Grains.m² Grains.spike WStem WGrain.spike
## 1 2435.265 207.5745 34.95 4530.411 21.82546 1.936 0.7628
## 2 3197.530 233.3964 35.95 4986.048 21.36300 2.138 0.7680
## Plants.m² DaysGrain BioRate GrainRateA GrainRateB VegRate KernelRate
## 1 142.7807 47 37.55741 16.81583 38.28285 40.58774 0.743617
## 2 137.9679 50 47.07561 19.21617 40.73828 57.09876 0.719000
## P.H.
## 1 727.8
## 2 711.2str(datos)## 'data.frame': 14 obs. of 40 variables:
## $ Bloque : Factor w/ 2 levels "1","2": 1 1 1 1 1 1 1 2 2 2 ...
## $ Parc : num 5 7 8 10 11 12 14 16 17 19 ...
## $ Trt : Factor w/ 7 levels "1","3","5","7",..: 1 6 4 7 3 5 2 4 1 3 ...
## $ Wdry50S : num 96.8 106.9 129.9 107.3 128.3 ...
## $ Wgrain50S : num 38.1 38.4 52.5 38.6 53.5 ...
## $ WidthArea : num 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5 ...
## $ LengthArea : num 9 9 9 9 9 9 9 9 9 9 ...
## $ WGrain : num 2370 2670 1732 3556 4277 ...
## $ WHumid : num 200 201 200 200 200 ...
## $ WDry : num 180 182 182 182 182 ...
## $ W200g : num 6.99 7.19 7.19 6.87 7.78 6.97 7.77 7.11 7.25 7.68 ...
## $ Area : num 13.5 13.5 13.5 13.5 13.5 13.5 13.5 13.5 13.5 13.5 ...
## $ AreaCount : num 0.187 0.187 0.187 0.187 0.125 0.187 0.187 0.187 0.187 0.125 ...
## $ Count : num 26.7 25.8 23.5 22.8 20.7 21.7 24.5 17.3 19.5 22.3 ...
## $ Tasseling : num 56 53 53 54 53 53 55 53 56 54 ...
## $ Lodging : num 0 0 0 0 0 0 0 0 0 0 ...
## $ Flowering : num 60 56 57 57 57 56 59 57 59 57 ...
## $ Maturity : num 107 106 107 108 107 106 107 106 106 107 ...
## $ Height : num 78.5 74.7 64.3 88.3 86.5 69.3 78 64.8 77.8 88.8 ...
## $ HI : num 39.4 35.9 40.5 36 41.7 ...
## $ X.Humidity : num 0.0981 0.0937 0.0925 0.0944 0.0939 ...
## $ H2O_Grain : num 232 250 160 336 401 ...
## $ Yield_Dry : num 1583 1792 1164 2385 2871 ...
## $ Yield_12_H2O: num 1799 2037 1323 2711 3262 ...
## $ Biomass : num 4019 4990 2878 6629 6888 ...
## $ Straw : num 2435 3198 1714 4244 4018 ...
## $ Spikes.m2 : num 208 233 111 309 268 ...
## $ Thou : num 35 36 36 34.4 38.9 ...
## $ Grains.m² : num 4530 4986 3239 6945 7380 ...
## $ Grains.spike: num 21.8 21.4 29.2 22.5 27.5 ...
## $ WStem : num 1.94 2.14 2.6 2.15 2.57 ...
## $ WGrain.spike: num 0.763 0.768 1.051 0.772 1.069 ...
## $ Plants.m² : num 143 138 126 122 166 ...
## $ DaysGrain : num 47 50 50 51 50 50 48 49 47 50 ...
## $ BioRate : num 37.6 47.1 26.9 61.4 64.4 ...
## $ GrainRateA : num 16.8 19.2 12.4 25.1 30.5 ...
## $ GrainRateB : num 38.3 40.7 26.5 53.2 65.2 ...
## $ VegRate : num 40.6 57.1 30.1 74.5 70.5 ...
## $ KernelRate : num 0.744 0.719 0.719 0.674 0.778 ...
## $ P.H. : num 728 711 726 709 729 ...names(datos)## [1] "Bloque" "Parc" "Trt" "Wdry50S"
## [5] "Wgrain50S" "WidthArea" "LengthArea" "WGrain"
## [9] "WHumid" "WDry" "W200g" "Area"
## [13] "AreaCount" "Count" "Tasseling" "Lodging"
## [17] "Flowering" "Maturity" "Height" "HI"
## [21] "X.Humidity" "H2O_Grain" "Yield_Dry" "Yield_12_H2O"
## [25] "Biomass" "Straw" "Spikes.m2" "Thou"
## [29] "Grains.m²" "Grains.spike" "WStem" "WGrain.spike"
## [33] "Plants.m²" "DaysGrain" "BioRate" "GrainRateA"
## [37] "GrainRateB" "VegRate" "KernelRate" "P.H."Con la función agregate() se pueden hacer resumenes por grupo. Las opciones son: sum, mean, median, sd, se, min, max.
aggregate(Yield_12_H2O ~ Trt, datos, mean)## Trt Yield_12_H2O
## 1 1 1750.456
## 2 3 2874.998
## 3 5 3364.312
## 4 7 1037.812
## 5 9 1681.568
## 6 11 2309.546
## 7 13 2693.390Utilizar la función boxplot para construir un diagrama
boxplot(datos$Yield_12_H2O ~ datos$Trt, xlab="Tratamiento",ylab="Rendimiento 12% H2O (kg ha-1)")
Construir un diagrama de interacción
interaction.plot(datos$Trt, datos$Bloque, datos$Yield_12_H2O, xlab="Tratamiento",ylab="Rendimiento 12% H2O (kg ha-1)")
interaction.plot(datos$Bloque, datos$Trt, datos$Yield_12_H2O, xlab="Bloque",ylab="Rendimiento 12% H2O (kg ha-1)")
Realizar un ANOVA de acuerdo al modelo para un Randomized Block design
anova(lm(datos$Yield_12_H2O ~ datos$Bloque + datos$Trt))## Analysis of Variance Table
##
## Response: datos$Yield_12_H2O
## Df Sum Sq Mean Sq F value Pr(>F)
## datos$Bloque 1 0 0 0.000 0.9995682
## datos$Trt 6 7748613 1291436 22.192 0.0007507 ***
## Residuals 6 349161 58193
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1La misma acción en dos pasos
ajuste <- lm(datos$Yield_12_H2O ~ datos$Bloque + datos$Trt)
anova(ajuste)## Analysis of Variance Table
##
## Response: datos$Yield_12_H2O
## Df Sum Sq Mean Sq F value Pr(>F)
## datos$Bloque 1 0 0 0.000 0.9995682
## datos$Trt 6 7748613 1291436 22.192 0.0007507 ***
## Residuals 6 349161 58193
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Prueba de medias con el paquete TukeyHSD()
prueba <- aov(datos$Yield_12_H2O ~ datos$Bloque + datos$Trt)
TukeyHSD(prueba, "datos$Trt", ordered = TRUE)## Tukey multiple comparisons of means
## 95% family-wise confidence level
## factor levels have been ordered
##
## Fit: aov(formula = datos$Yield_12_H2O ~ datos$Bloque + datos$Trt)
##
## $`datos$Trt`
## diff lwr upr p adj
## 9-7 643.75649 -361.851512 1649.364 0.2499356
## 1-7 712.64420 -292.963805 1718.252 0.1833526
## 11-7 1271.73420 266.126197 2277.342 0.0171362
## 13-7 1655.57810 649.970096 2661.186 0.0045229
## 3-7 1837.18615 831.578146 2842.794 0.0025999
## 5-7 2326.49976 1320.891755 3332.108 0.0007094
## 1-9 68.88771 -936.720294 1074.496 0.9999073
## 11-9 627.97771 -377.630292 1633.586 0.2680896
## 13-9 1011.82161 6.213607 2017.430 0.0486973
## 3-9 1193.42966 187.821657 2199.038 0.0231802
## 5-9 1682.74327 677.135266 2688.351 0.0041515
## 11-1 559.09000 -446.517999 1564.698 0.3615919
## 13-1 942.93390 -62.674100 1948.542 0.0654782
## 3-1 1124.54195 118.933950 2130.150 0.0305106
## 5-1 1613.85556 608.247559 2619.464 0.0051696
## 13-11 383.84390 -621.764101 1389.452 0.6938596
## 3-11 565.45195 -440.156052 1571.060 0.3519511
## 5-11 1054.76556 49.157558 2060.374 0.0406471
## 3-13 181.60805 -823.999951 1187.216 0.9821915
## 5-13 670.92166 -334.686342 1676.530 0.2213177
## 5-3 489.31361 -516.294391 1494.922 0.4805923plot(TukeyHSD(prueba, "datos$Trt", ordered = TRUE))
Prueba de medias con el paquete agricolae
#install.packages("agricolae")
library(agricolae)
comparaciones <- HSD.test(ajuste, "datos$Trt")
comparaciones## $statistics
## Mean CV MSerror HSD
## 2244.583 10.74735 58193.47 1005.608
##
## $parameters
## Df ntr StudentizedRange alpha test name.t
## 6 7 5.895309 0.05 Tukey datos$Trt
##
## $means
## datos$Yield_12_H2O std r Min Max
## 1 1750.456 69.06707 2 1701.6183 1799.294
## 11 2309.546 385.56019 2 2036.9138 2582.178
## 13 2693.390 24.61455 2 2675.9849 2710.795
## 3 2874.998 58.10076 2 2833.9146 2916.081
## 5 3364.312 144.38696 2 3262.2146 3466.409
## 7 1037.812 403.34192 2 752.6061 1323.018
## 9 1681.568 90.66437 2 1617.4590 1745.678
##
## $comparison
## NULL
##
## $groups
## trt means M
## 1 5 3364.312 a
## 2 3 2874.998 ab
## 3 13 2693.390 abc
## 4 11 2309.546 bcd
## 5 1 1750.456 cde
## 6 9 1681.568 de
## 7 7 1037.812 e