Este análisis
de varianza involucra 4 localidades. “La localidad tiene un propósito,
por lo que no se puede considerar una variable aleatoria”, ver pp. 184 – 190.
 |
|
Instalamos
el paquete “dae” en R para graficar la interacción entre las variables localidad, época y
variedad. También, empleamos los códigos ‘sp.plot’ y ‘ssp.plot’ de “agricolae”
para analizar datos como parcelas divididas.
Códigos para R
#***********************************************************
# EXPERIMENTACIÓN EN AGRICULTURA
# Fernández Escobar, R.; Trapero, A.; Domínguez, J. 2010
# CAPÍTULO 13: ANÁLISIS DE LA VARIANZA COMBINADO
# Análisis de diferentes localidades (páginas 184 - 190)
#************************************************************
#### Importar Tabla 13.5
local <- read.csv("C:/Users/Administrator/Desktop/tabla 13.5.csv", header=TRUE)
attach(local)
print(local)## LOCALIDAD EPOCA VARIEDAD BLOQUE RENDIMIENTO
## 1 L1 invierno V1 1 0
## 2 L1 invierno V1 2 214
## 3 L1 invierno V1 3 425
## 4 L1 invierno V1 4 40
## 5 L1 invierno V2 1 1252
## 6 L1 invierno V2 2 627
## 7 L1 invierno V2 3 716
## 8 L1 invierno V2 4 1068
## 9 L1 invierno V3 1 2163
## 10 L1 invierno V3 2 2714
## 11 L1 invierno V3 3 2521
## 12 L1 invierno V3 4 2240
## 13 L1 primavera V1 1 1036
## 14 L1 primavera V1 2 697
## 15 L1 primavera V1 3 849
## 16 L1 primavera V1 4 1258
## 17 L1 primavera V2 1 1524
## 18 L1 primavera V2 2 1861
## 19 L1 primavera V2 3 2220
## 20 L1 primavera V2 4 1744
## 21 L1 primavera V3 1 1312
## 22 L1 primavera V3 2 874
## 23 L1 primavera V3 3 695
## 24 L1 primavera V3 4 1133
## 25 L2 invierno V1 1 26
## 26 L2 invierno V1 2 125
## 27 L2 invierno V1 3 0
## 28 L2 invierno V1 4 358
## 29 L2 invierno V2 1 813
## 30 L2 invierno V2 2 986
## 31 L2 invierno V2 3 657
## 32 L2 invierno V2 4 432
## 33 L2 invierno V3 1 2475
## 34 L2 invierno V3 2 2671
## 35 L2 invierno V3 3 3044
## 36 L2 invierno V3 4 2850
## 37 L2 primavera V1 1 488
## 38 L2 primavera V1 2 369
## 39 L2 primavera V1 3 714
## 40 L2 primavera V1 4 983
## 41 L2 primavera V2 1 1815
## 42 L2 primavera V2 2 1662
## 43 L2 primavera V2 3 1233
## 44 L2 primavera V2 4 1486
## 45 L2 primavera V3 1 1958
## 46 L2 primavera V3 2 1604
## 47 L2 primavera V3 3 1760
## 48 L2 primavera V3 4 1325
## 49 L3 invierno V1 1 257
## 50 L3 invierno V1 2 64
## 51 L3 invierno V1 3 756
## 52 L3 invierno V1 4 452
## 53 L3 invierno V2 1 789
## 54 L3 invierno V2 2 491
## 55 L3 invierno V2 3 608
## 56 L3 invierno V2 4 1066
## 57 L3 invierno V3 1 2801
## 58 L3 invierno V3 2 3125
## 59 L3 invierno V3 3 2977
## 60 L3 invierno V3 4 3453
## 61 L3 primavera V1 1 412
## 62 L3 primavera V1 2 1118
## 63 L3 primavera V1 3 833
## 64 L3 primavera V1 4 621
## 65 L3 primavera V2 1 1325
## 66 L3 primavera V2 2 1950
## 67 L3 primavera V2 3 1507
## 68 L3 primavera V2 4 1773
## 69 L3 primavera V3 1 1582
## 70 L3 primavera V3 2 1236
## 71 L3 primavera V3 3 1844
## 72 L3 primavera V3 4 1371
## 73 L4 invierno V1 1 1060
## 74 L4 invierno V1 2 1383
## 75 L4 invierno V1 3 889
## 76 L4 invierno V1 4 723
## 77 L4 invierno V2 1 816
## 78 L4 invierno V2 2 1125
## 79 L4 invierno V2 3 1463
## 80 L4 invierno V2 4 1277
## 81 L4 invierno V3 1 892
## 82 L4 invierno V3 2 1112
## 83 L4 invierno V3 3 1521
## 84 L4 invierno V3 4 1366
## 85 L4 primavera V1 1 960
## 86 L4 primavera V1 2 474
## 87 L4 primavera V1 3 899
## 88 L4 primavera V1 4 763
## 89 L4 primavera V2 1 725
## 90 L4 primavera V2 2 583
## 91 L4 primavera V2 3 1179
## 92 L4 primavera V2 4 914
## 93 L4 primavera V3 1 897
## 94 L4 primavera V3 2 985
## 95 L4 primavera V3 3 544
## 96 L4 primavera V3 4 718#### Creando variables independientes como factores VARIEDAD <- factor(VARIEDAD , levels = c('V1','V2','V3'), labels = c("Variedad 1", "Variedad 2", "Variedad 3")) LOCALIDAD <- factor(LOCALIDAD, levels = c('L1','L2','L3','L4'), labels = c("Localidad 1", "Localidad 2", "Localidad 3", "Localidad 4")) BLOQUE <- factor(BLOQUE) EPOCA<- factor(EPOCA) LOCAL <- data.frame(LOCALIDAD, EPOCA,VARIEDAD, BLOQUE, RENDIMIENTO)
attach(LOCAL); str(LOCAL)#### Agregando función para calcular estadísticas descriptivas de las interacciones
library(plyr)source("C:/Users/Administrator/Desktop/summarySE.r")
MEDIAS <- summarySE(LOCAL, measurevar="RENDIMIENTO", groupvars=c("LOCALIDAD", "VARIEDAD", "EPOCA")); MEDIAS## LOCALIDAD VARIEDAD EPOCA N RENDIMIENTO sd se
## 1 Localidad 1 Variedad 1 invierno 4 169.75 193.8735 96.93673
## 2 Localidad 1 Variedad 1 primavera 4 960.00 242.2602 121.13010
## 3 Localidad 1 Variedad 2 invierno 4 915.75 294.1206 147.06029
## 4 Localidad 1 Variedad 2 primavera 4 1837.25 290.9082 145.45410
## 5 Localidad 1 Variedad 3 invierno 4 2409.50 254.7188 127.35940
## 6 Localidad 1 Variedad 3 primavera 4 1003.50 273.1819 136.59093
## 7 Localidad 2 Variedad 1 invierno 4 127.25 162.9875 81.49374
## 8 Localidad 2 Variedad 1 primavera 4 638.50 270.5926 135.29628
## 9 Localidad 2 Variedad 2 invierno 4 722.00 235.4443 117.72213
## 10 Localidad 2 Variedad 2 primavera 4 1549.00 249.9000 124.94999
## 11 Localidad 2 Variedad 3 invierno 4 2760.00 243.5173 121.75864
## 12 Localidad 2 Variedad 3 primavera 4 1661.75 267.1783 133.58916
## 13 Localidad 3 Variedad 1 invierno 4 382.25 295.2540 147.62699
## 14 Localidad 3 Variedad 1 primavera 4 746.00 301.7361 150.86804
## 15 Localidad 3 Variedad 2 invierno 4 738.50 250.3950 125.19751
## 16 Localidad 3 Variedad 2 primavera 4 1638.75 277.3065 138.65327
## 17 Localidad 3 Variedad 3 invierno 4 3089.00 276.4537 138.22687
## 18 Localidad 3 Variedad 3 primavera 4 1508.25 265.2827 132.64136
## 19 Localidad 4 Variedad 1 invierno 4 1013.75 282.0064 141.00318
## 20 Localidad 4 Variedad 1 primavera 4 774.00 216.2884 108.14419
## 21 Localidad 4 Variedad 2 invierno 4 1170.25 273.6413 136.82067
## 22 Localidad 4 Variedad 2 primavera 4 850.25 257.7148 128.85740
## 23 Localidad 4 Variedad 3 invierno 4 1222.75 277.5697 138.78483
## 24 Localidad 4 Variedad 3 primavera 4 786.00 195.8826 97.94131# Otras medias
library(psych)describeBy(RENDIMIENTO, LOCALIDAD:EPOCA:VARIEDAD, mat = TRUE, digits = 2)## item group1 vars n mean sd median
## X11 1 Localidad 1:invierno:Variedad 1 1 4 169.75 193.87 127.0
## X12 2 Localidad 1:invierno:Variedad 2 1 4 915.75 294.12 892.0
## X13 3 Localidad 1:invierno:Variedad 3 1 4 2409.50 254.72 2380.5
## X14 4 Localidad 1:primavera:Variedad 1 1 4 960.00 242.26 942.5
## X15 5 Localidad 1:primavera:Variedad 2 1 4 1837.25 290.91 1802.5
## X16 6 Localidad 1:primavera:Variedad 3 1 4 1003.50 273.18 1003.5
## X17 7 Localidad 2:invierno:Variedad 1 1 4 127.25 162.99 75.5
## X18 8 Localidad 2:invierno:Variedad 2 1 4 722.00 235.44 735.0
## X19 9 Localidad 2:invierno:Variedad 3 1 4 2760.00 243.52 2760.5
## X110 10 Localidad 2:primavera:Variedad 1 1 4 638.50 270.59 601.0
## X111 11 Localidad 2:primavera:Variedad 2 1 4 1549.00 249.90 1574.0
## X112 12 Localidad 2:primavera:Variedad 3 1 4 1661.75 267.18 1682.0
## X113 13 Localidad 3:invierno:Variedad 1 1 4 382.25 295.25 354.5
## X114 14 Localidad 3:invierno:Variedad 2 1 4 738.50 250.40 698.5
## X115 15 Localidad 3:invierno:Variedad 3 1 4 3089.00 276.45 3051.0
## X116 16 Localidad 3:primavera:Variedad 1 1 4 746.00 301.74 727.0
## X117 17 Localidad 3:primavera:Variedad 2 1 4 1638.75 277.31 1640.0
## X118 18 Localidad 3:primavera:Variedad 3 1 4 1508.25 265.28 1476.5
## X119 19 Localidad 4:invierno:Variedad 1 1 4 1013.75 282.01 974.5
## X120 20 Localidad 4:invierno:Variedad 2 1 4 1170.25 273.64 1201.0
## X121 21 Localidad 4:invierno:Variedad 3 1 4 1222.75 277.57 1239.0
## X122 22 Localidad 4:primavera:Variedad 1 1 4 774.00 216.29 831.0
## X123 23 Localidad 4:primavera:Variedad 2 1 4 850.25 257.71 819.5
## X124 24 Localidad 4:primavera:Variedad 3 1 4 786.00 195.88 807.5# Gráficos
library(dae)interaction.ABC.plot(RENDIMIENTO, x.factor = LOCALIDAD,
groups.factor = EPOCA,
trace.factor = VARIEDAD,
data = LOCAL,
title = "Rendimiento de Cultivares (kg/ha) en diferentes LOCALIDAD:EPOCA")
interaction.ABC.plot(RENDIMIENTO, x.factor = EPOCA,
groups.factor = LOCALIDAD,
trace.factor = VARIEDAD,
data = LOCAL,
title = "Rendimiento de Cultivares (kg/ha) en diferentes ÉPOCA:LOCALIDAD")
interaction.ABC.plot(RENDIMIENTO, x.factor = VARIEDAD,
groups.factor = LOCALIDAD,
trace.factor = EPOCA,
data = LOCAL,
title = "Rendimiento de Cultivares (kg/ha) por ÉPOCA")
interaction.ABC.plot(RENDIMIENTO, x.factor = VARIEDAD,
groups.factor = EPOCA,
trace.factor = LOCALIDAD,
data = LOCAL,
title = "Rendimiento de Cultivares (kg/ha) por LOCALIDAD")
#### ANOVA para LOCALIDAD 1
# Separando datos
LOCALuno <- LOCAL[1:24, ]; LOCALuno## LOCALIDAD EPOCA BLOQUE VARIEDAD RENDIMIENTO
## 1 Localidad 1 invierno 1 Variedad 1 0
## 2 Localidad 1 invierno 2 Variedad 1 214
## 3 Localidad 1 invierno 3 Variedad 1 425
## 4 Localidad 1 invierno 4 Variedad 1 40
## 5 Localidad 1 invierno 1 Variedad 2 1252
## 6 Localidad 1 invierno 2 Variedad 2 627
## 7 Localidad 1 invierno 3 Variedad 2 716
## 8 Localidad 1 invierno 4 Variedad 2 1068
## 9 Localidad 1 invierno 1 Variedad 3 2163
## 10 Localidad 1 invierno 2 Variedad 3 2714
## 11 Localidad 1 invierno 3 Variedad 3 2521
## 12 Localidad 1 invierno 4 Variedad 3 2240
## 13 Localidad 1 primavera 1 Variedad 1 1036
## 14 Localidad 1 primavera 2 Variedad 1 697
## 15 Localidad 1 primavera 3 Variedad 1 849
## 16 Localidad 1 primavera 4 Variedad 1 1258
## 17 Localidad 1 primavera 1 Variedad 2 1524
## 18 Localidad 1 primavera 2 Variedad 2 1861
## 19 Localidad 1 primavera 3 Variedad 2 2220
## 20 Localidad 1 primavera 4 Variedad 2 1744
## 21 Localidad 1 primavera 1 Variedad 3 1312
## 22 Localidad 1 primavera 2 Variedad 3 874
## 23 Localidad 1 primavera 3 Variedad 3 695
## 24 Localidad 1 primavera 4 Variedad 3 1133# Estadísticas descriptivas para EPOCA:VARIEDAD
source("C:/Users/Administrator/Desktop/schwarz.functions.r") # Función Schwarz
MEDIAuno <- ddply(LOCALuno, c("EPOCA","VARIEDAD"), sf.simple.summary, variable="RENDIMIENTO"); MEDIAuno## EPOCA VARIEDAD n nmiss mean sd
## 1 invierno Variedad 1 4 0 169.75 193.8735
## 2 invierno Variedad 2 4 0 915.75 294.1206
## 3 invierno Variedad 3 4 0 2409.50 254.7188
## 4 primavera Variedad 1 4 0 960.00 242.2602
## 5 primavera Variedad 2 4 0 1837.25 290.9082
## 6 primavera Variedad 3 4 0 1003.50 273.1819# Gráfico de interacción simple
library(dae)
interaction.ABC.plot(RENDIMIENTO, x.factor = VARIEDAD,
groups.factor = EPOCA,
trace.factor = VARIEDAD,
data = LOCALuno, title = "Effecto Epoca en Rendimiento de variedades")
interaction.plot(LOCALuno$VARIEDAD, LOCALuno$EPOCA, LOCALuno$RENDIMIENTO, fixed = TRUE,
col = 2:3, leg.bty = "o", xlab = "Cultivares", ylab = "Rendimento (kg/ha)")
# ANOVA Diseño en Parcelas Divididas
library(agricolae)modelo1 <- with(LOCALuno, sp.plot(BLOQUE,EPOCA,VARIEDAD,RENDIMIENTO)) # ANOVA para diseño en parcelas divididas 'Agricolae'##
## ANALYSIS SPLIT PLOT: RENDIMIENTO
## Class level information
##
## EPOCA : invierno primavera
## VARIEDAD : Variedad 1 Variedad 2 Variedad 3
## BLOQUE : 1 2 3 4
##
## Number of observations: 24
##
## Analysis of Variance Table
##
## Response: RENDIMIENTO
## Df Sum Sq Mean Sq F value Pr(>F)
## BLOQUE 3 24572 8191 0.3073 0.8208
## EPOCA 1 62322 62322 2.3380 0.2237
## Ea 3 79970 26657
## VARIEDAD 2 5522514 2761257 29.6850 2.259e-05 ***
## EPOCA:VARIEDAD 2 6838665 3419332 36.7597 7.633e-06 ***
## Eb 12 1116223 93019
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## cv(a) = 13.4 %, cv(b) = 25.1 %, Mean = 1215.958modelo2 <- aov(RENDIMIENTO ~ BLOQUE + BLOQUE*EPOCA + EPOCA*VARIEDAD, data = LOCALuno)
summary(modelo2)## Df Sum Sq Mean Sq F value Pr(>F)
## BLOQUE 3 24572 8191 0.088 0.965
## EPOCA 1 62322 62322 0.670 0.429
## VARIEDAD 2 5522514 2761257 29.685 2.26e-05 ***
## BLOQUE:EPOCA 3 79970 26657 0.287 0.834
## EPOCA:VARIEDAD 2 6838665 3419332 36.760 7.63e-06 ***
## Residuals 12 1116223 93019
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Gráfico de diagnósticos
par(mfrow=c(1,2))
plot(modelo2)
par(mfrow=c(1,1))
#### ANOVA para LOCALIDAD 2
# Separando datos
LOCALdos <- LOCAL[25:48, ]; LOCALdos## LOCALIDAD EPOCA BLOQUE VARIEDAD RENDIMIENTO
## 25 Localidad 2 invierno 1 Variedad 1 26
## 26 Localidad 2 invierno 2 Variedad 1 125
## 27 Localidad 2 invierno 3 Variedad 1 0
## 28 Localidad 2 invierno 4 Variedad 1 358
## 29 Localidad 2 invierno 1 Variedad 2 813
## 30 Localidad 2 invierno 2 Variedad 2 986
## 31 Localidad 2 invierno 3 Variedad 2 657
## 32 Localidad 2 invierno 4 Variedad 2 432
## 33 Localidad 2 invierno 1 Variedad 3 2475
## 34 Localidad 2 invierno 2 Variedad 3 2671
## 35 Localidad 2 invierno 3 Variedad 3 3044
## 36 Localidad 2 invierno 4 Variedad 3 2850
## 37 Localidad 2 primavera 1 Variedad 1 488
## 38 Localidad 2 primavera 2 Variedad 1 369
## 39 Localidad 2 primavera 3 Variedad 1 714
## 40 Localidad 2 primavera 4 Variedad 1 983
## 41 Localidad 2 primavera 1 Variedad 2 1815
## 42 Localidad 2 primavera 2 Variedad 2 1662
## 43 Localidad 2 primavera 3 Variedad 2 1233
## 44 Localidad 2 primavera 4 Variedad 2 1486
## 45 Localidad 2 primavera 1 Variedad 3 1958
## 46 Localidad 2 primavera 2 Variedad 3 1604
## 47 Localidad 2 primavera 3 Variedad 3 1760
## 48 Localidad 2 primavera 4 Variedad 3 1325# Estadísticas descriptivas
MEDIAdos <- ddply(LOCALdos, c("EPOCA","VARIEDAD"), sf.simple.summary, variable="RENDIMIENTO"); MEDIAdos## EPOCA VARIEDAD n nmiss mean sd
## 1 invierno Variedad 1 4 0 127.25 162.9875
## 2 invierno Variedad 2 4 0 722.00 235.4443
## 3 invierno Variedad 3 4 0 2760.00 243.5173
## 4 primavera Variedad 1 4 0 638.50 270.5926
## 5 primavera Variedad 2 4 0 1549.00 249.9000
## 6 primavera Variedad 3 4 0 1661.75 267.1783# Interacción
interaction.ABC.plot(RENDIMIENTO, x.factor = VARIEDAD,
groups.factor = EPOCA,
trace.factor = VARIEDAD,
data = LOCALdos, title = "Effecto Epoca en Rendimiento de variedades")
interaction.plot(LOCALdos$VARIEDAD, LOCALdos$EPOCA, LOCALdos$RENDIMIENTO, fixed = TRUE,
col = 2:3, leg.bty = "o", xlab = "Cultivares", ylab = "Rendimento (kg/ha)")
# ANOVA Diseño en Parcelas Divididas
modelo1 <- with(LOCALdos, sp.plot(BLOQUE,EPOCA,VARIEDAD,RENDIMIENTO)) # ANOVA para diseño en parcelas divididas 'Agricolae'##
## ANALYSIS SPLIT PLOT: RENDIMIENTO
## Class level information
##
## EPOCA : invierno primavera
## VARIEDAD : Variedad 1 Variedad 2 Variedad 3
## BLOQUE : 1 2 3 4
##
## Number of observations: 24
##
## Analysis of Variance Table
##
## Response: RENDIMIENTO
## Df Sum Sq Mean Sq F value Pr(>F)
## BLOQUE 3 3074 1025 0.0259 0.9932
## EPOCA 1 38400 38400 0.9711 0.3971
## Ea 3 118628 39543
## VARIEDAD 2 13505226 6752613 87.7571 6.869e-08 ***
## EPOCA:VARIEDAD 2 4264517 2132259 27.7109 3.179e-05 ***
## Eb 12 923360 76947
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## cv(a) = 16 %, cv(b) = 22.3 %, Mean = 1243.083modelo2 <- aov(RENDIMIENTO ~ BLOQUE + BLOQUE*EPOCA + EPOCA*VARIEDAD, data = LOCALdos)
summary(modelo2)## Df Sum Sq Mean Sq F value Pr(>F)
## BLOQUE 3 3074 1025 0.013 0.998
## EPOCA 1 38400 38400 0.499 0.493
## VARIEDAD 2 13505226 6752613 87.757 6.87e-08 ***
## BLOQUE:EPOCA 3 118628 39543 0.514 0.680
## EPOCA:VARIEDAD 2 4264517 2132259 27.711 3.18e-05 ***
## Residuals 12 923360 76947
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Gráfico de diagnósticos
par(mfrow=c(1,2))
plot(modelo2)
par(mfrow=c(1,1))
#### ANOVA para LOCALIDAD 3
# Separando datos
LOCALtres <- LOCAL[49:72, ]; LOCALtres## LOCALIDAD EPOCA BLOQUE VARIEDAD RENDIMIENTO
## 49 Localidad 3 invierno 1 Variedad 1 257
## 50 Localidad 3 invierno 2 Variedad 1 64
## 51 Localidad 3 invierno 3 Variedad 1 756
## 52 Localidad 3 invierno 4 Variedad 1 452
## 53 Localidad 3 invierno 1 Variedad 2 789
## 54 Localidad 3 invierno 2 Variedad 2 491
## 55 Localidad 3 invierno 3 Variedad 2 608
## 56 Localidad 3 invierno 4 Variedad 2 1066
## 57 Localidad 3 invierno 1 Variedad 3 2801
## 58 Localidad 3 invierno 2 Variedad 3 3125
## 59 Localidad 3 invierno 3 Variedad 3 2977
## 60 Localidad 3 invierno 4 Variedad 3 3453
## 61 Localidad 3 primavera 1 Variedad 1 412
## 62 Localidad 3 primavera 2 Variedad 1 1118
## 63 Localidad 3 primavera 3 Variedad 1 833
## 64 Localidad 3 primavera 4 Variedad 1 621
## 65 Localidad 3 primavera 1 Variedad 2 1325
## 66 Localidad 3 primavera 2 Variedad 2 1950
## 67 Localidad 3 primavera 3 Variedad 2 1507
## 68 Localidad 3 primavera 4 Variedad 2 1773
## 69 Localidad 3 primavera 1 Variedad 3 1582
## 70 Localidad 3 primavera 2 Variedad 3 1236
## 71 Localidad 3 primavera 3 Variedad 3 1844
## 72 Localidad 3 primavera 4 Variedad 3 1371# Estadísticas descriptivas
MEDIAtres <- ddply(LOCALtres, c("EPOCA","VARIEDAD"), sf.simple.summary, variable="RENDIMIENTO"); MEDIAtres## EPOCA VARIEDAD n nmiss mean sd
## 1 invierno Variedad 1 4 0 382.25 295.2540
## 2 invierno Variedad 2 4 0 738.50 250.3950
## 3 invierno Variedad 3 4 0 3089.00 276.4537
## 4 primavera Variedad 1 4 0 746.00 301.7361
## 5 primavera Variedad 2 4 0 1638.75 277.3065
## 6 primavera Variedad 3 4 0 1508.25 265.2827# Interacción
interaction.ABC.plot(RENDIMIENTO, x.factor = VARIEDAD,
groups.factor = EPOCA,
trace.factor = VARIEDAD,
data = LOCALtres, title = "Effecto Epoca en Rendimiento de variedades")
interaction.plot(LOCALtres$VARIEDAD, LOCALtres$EPOCA, LOCALtres$RENDIMIENTO, fixed = TRUE,
col = 2:3, leg.bty = "o", xlab = "Cultivares", ylab = "Rendimento (kg/ha)")
# ANOVA Diseño en Parcelas Divididas
modelo1 <- with(LOCALtres, sp.plot(BLOQUE,EPOCA,VARIEDAD,RENDIMIENTO)) # ANOVA para diseño en parcelas divididas 'Agricolae'##
## ANALYSIS SPLIT PLOT: RENDIMIENTO
## Class level information
##
## EPOCA : invierno primavera
## VARIEDAD : Variedad 1 Variedad 2 Variedad 3
## BLOQUE : 1 2 3 4
##
## Number of observations: 24
##
## Analysis of Variance Table
##
## Response: RENDIMIENTO
## Df Sum Sq Mean Sq F value Pr(>F)
## BLOQUE 3 245150 81717 0.8425 0.5544
## EPOCA 1 66887 66887 0.6896 0.4672
## Ea 3 290987 96996
## VARIEDAD 2 12348241 6174121 86.3800 7.507e-08 ***
## EPOCA:VARIEDAD 2 6816182 3408091 47.6814 1.950e-06 ***
## Eb 12 857716 71476
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## cv(a) = 23.1 %, cv(b) = 19.8 %, Mean = 1350.458modelo2 <- aov(RENDIMIENTO ~ BLOQUE + BLOQUE*EPOCA + EPOCA*VARIEDAD, data = LOCALtres)
summary(modelo2)## Df Sum Sq Mean Sq F value Pr(>F)
## BLOQUE 3 245150 81717 1.143 0.371
## EPOCA 1 66887 66887 0.936 0.352
## VARIEDAD 2 12348241 6174121 86.380 7.51e-08 ***
## BLOQUE:EPOCA 3 290987 96996 1.357 0.303
## EPOCA:VARIEDAD 2 6816182 3408091 47.681 1.95e-06 ***
## Residuals 12 857716 71476
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#### Gráfico de diagnósticos
par(mfrow=c(1,2))
plot(modelo2)
par(mfrow=c(1,1))
#### ANOVA para LOCALIDAD 4
# Separando datos
LOCALcuatro <- LOCAL[73:96, ]; LOCALcuatro## LOCALIDAD EPOCA BLOQUE VARIEDAD RENDIMIENTO
## 73 Localidad 4 invierno 1 Variedad 1 1060
## 74 Localidad 4 invierno 2 Variedad 1 1383
## 75 Localidad 4 invierno 3 Variedad 1 889
## 76 Localidad 4 invierno 4 Variedad 1 723
## 77 Localidad 4 invierno 1 Variedad 2 816
## 78 Localidad 4 invierno 2 Variedad 2 1125
## 79 Localidad 4 invierno 3 Variedad 2 1463
## 80 Localidad 4 invierno 4 Variedad 2 1277
## 81 Localidad 4 invierno 1 Variedad 3 892
## 82 Localidad 4 invierno 2 Variedad 3 1112
## 83 Localidad 4 invierno 3 Variedad 3 1521
## 84 Localidad 4 invierno 4 Variedad 3 1366
## 85 Localidad 4 primavera 1 Variedad 1 960
## 86 Localidad 4 primavera 2 Variedad 1 474
## 87 Localidad 4 primavera 3 Variedad 1 899
## 88 Localidad 4 primavera 4 Variedad 1 763
## 89 Localidad 4 primavera 1 Variedad 2 725
## 90 Localidad 4 primavera 2 Variedad 2 583
## 91 Localidad 4 primavera 3 Variedad 2 1179
## 92 Localidad 4 primavera 4 Variedad 2 914
## 93 Localidad 4 primavera 1 Variedad 3 897
## 94 Localidad 4 primavera 2 Variedad 3 985
## 95 Localidad 4 primavera 3 Variedad 3 544
## 96 Localidad 4 primavera 4 Variedad 3 718## Estadísticas descriptivas
MEDIAcuatro <- ddply(LOCALcuatro, c("EPOCA","VARIEDAD"), sf.simple.summary, variable="RENDIMIENTO"); MEDIAcuatro## EPOCA VARIEDAD n nmiss mean sd
## 1 invierno Variedad 1 4 0 1013.75 282.0064
## 2 invierno Variedad 2 4 0 1170.25 273.6413
## 3 invierno Variedad 3 4 0 1222.75 277.5697
## 4 primavera Variedad 1 4 0 774.00 216.2884
## 5 primavera Variedad 2 4 0 850.25 257.7148
## 6 primavera Variedad 3 4 0 786.00 195.8826# Interacción
interaction.ABC.plot(RENDIMIENTO, x.factor = VARIEDAD,
groups.factor = EPOCA,
trace.factor = VARIEDAD,
data = LOCALcuatro, title = "Effecto Epoca en Rendimiento de variedades")
interaction.plot(LOCALcuatro$VARIEDAD, LOCALcuatro$EPOCA, LOCALcuatro$RENDIMIENTO, fixed = TRUE,
col = 2:3, leg.bty = "o", xlab = "Cultivares", ylab = "Rendimento (kg/ha)")
# ANOVA Diseño en Parcelas Divididas
modelo1 <- with(LOCALcuatro,sp.plot(BLOQUE,EPOCA,VARIEDAD,RENDIMIENTO)) # ANOVA para diseño en parcelas divididas 'Agricolae'##
## ANALYSIS SPLIT PLOT: RENDIMIENTO
## Class level information
##
## EPOCA : invierno primavera
## VARIEDAD : Variedad 1 Variedad 2 Variedad 3
## BLOQUE : 1 2 3 4
##
## Number of observations: 24
##
## Analysis of Variance Table
##
## Response: RENDIMIENTO
## Df Sum Sq Mean Sq F value Pr(>F)
## BLOQUE 3 117489 39163 0.6647 0.62734
## EPOCA 1 662008 662008 11.2366 0.04399 *
## Ea 3 176745 58915
## VARIEDAD 2 68768 34384 0.4827 0.62861
## EPOCA:VARIEDAD 2 39253 19627 0.2755 0.76385
## Eb 12 854825 71235
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## cv(a) = 25 %, cv(b) = 27.5 %, Mean = 969.5modelo2 <- aov(RENDIMIENTO ~ BLOQUE + BLOQUE*EPOCA + EPOCA*VARIEDAD, data = LOCALcuatro)
summary(modelo2)## Df Sum Sq Mean Sq F value Pr(>F)
## BLOQUE 3 117489 39163 0.550 0.6578
## EPOCA 1 662008 662008 9.293 0.0101 *
## VARIEDAD 2 68768 34384 0.483 0.6286
## BLOQUE:EPOCA 3 176745 58915 0.827 0.5040
## EPOCA:VARIEDAD 2 39253 19627 0.276 0.7639
## Residuals 12 854825 71235
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Gráfico de diagnósticos
par(mfrow=c(1,2))
plot(modelo2)
par(mfrow=c(1,1))
dev.off()#### ANOVA combinado para los datos del ejemplo sería (pp. 187):
# Opción 1
modelo1 <- with(LOCAL, ssp.plot(LOCALIDAD, BLOQUE, EPOCA, VARIEDAD, RENDIMIENTO)) # Agricolae: Split-Split-Plot Analysis## ANALYSIS SPLIT-SPLIT PLOT: RENDIMIENTO
## Class level information
##
## BLOQUE : 1 2 3 4
## EPOCA : invierno primavera
## VARIEDAD : Variedad 1 Variedad 2 Variedad 3
## LOCALIDAD : Localidad 1 Localidad 2 Localidad 3 Localidad 4
##
## Number of observations: 96
##
## Analysis of Variance Table
##
## Response: RENDIMIENTO
## Df Sum Sq Mean Sq F value Pr(>F)
## LOCALIDAD 3 1866445 622148 25.0817 0.0001052 ***
## BLOQUE 3 167041 55680 2.2447 0.1522928
## Ea 9 223244 24805
## EPOCA 1 98176 98176 0.9207 0.3562321
## BLOQUE:EPOCA 3 118156 39385 0.3693 0.7765430
## Eb 12 1279616 106635
## VARIEDAD 2 23194537 11597268 36.2394 2.558e-10 ***
## VARIEDAD:BLOQUE 6 29262 4877 0.0152 0.9999827
## VARIEDAD:EPOCA 2 13852499 6926249 21.6433 1.995e-07 ***
## VARIEDAD:BLOQUE:EPOCA 6 718295 119716 0.3741 0.8918695
## Ec 48 15360896 320019
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## cv(a) = 13.2 %, cv(b) = 27.3 %, cv(c) = 47.3 %, Mean = 1194.75# Opción 2
modelo2 <- aov(RENDIMIENTO ~ LOCALIDAD + BLOQUE + EPOCA*LOCALIDAD +
BLOQUE*EPOCA + LOCALIDAD*EPOCA*VARIEDAD, data = LOCAL); summary(modelo2)## Df Sum Sq Mean Sq F value Pr(>F)
## LOCALIDAD 3 1866445 622148 9.077 4.11e-05 ***
## BLOQUE 3 167041 55680 0.812 0.4915
## EPOCA 1 98176 98176 1.432 0.2357
## VARIEDAD 2 23194537 11597268 169.208 < 2e-16 ***
## LOCALIDAD:EPOCA 3 731441 243814 3.557 0.0189 *
## BLOQUE:EPOCA 3 118156 39385 0.575 0.6337
## LOCALIDAD:VARIEDAD 6 8250212 1375035 20.062 3.40e-13 ***
## EPOCA:VARIEDAD 2 13852499 6926249 101.056 < 2e-16 ***
## LOCALIDAD:EPOCA:VARIEDAD 6 4106118 684353 9.985 7.95e-08 ***
## Residuals 66 4523542 68539
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Opción 3
modelo3 <- aov(RENDIMIENTO ~ LOCALIDAD + BLOQUE/LOCALIDAD + EPOCA*LOCALIDAD +
LOCALIDAD*EPOCA*VARIEDAD, data = LOCAL); summary(modelo3)## Df Sum Sq Mean Sq F value Pr(>F)
## LOCALIDAD 3 1866445 622148 8.448 9.07e-05 ***
## BLOQUE 3 167041 55680 0.756 0.523
## EPOCA 1 98176 98176 1.333 0.253
## VARIEDAD 2 23194537 11597268 157.484 < 2e-16 ***
## LOCALIDAD:BLOQUE 9 223244 24805 0.337 0.959
## LOCALIDAD:EPOCA 3 731441 243814 3.311 0.026 *
## LOCALIDAD:VARIEDAD 6 8250212 1375035 18.672 4.11e-12 ***
## EPOCA:VARIEDAD 2 13852499 6926249 94.054 < 2e-16 ***
## LOCALIDAD:EPOCA:VARIEDAD 6 4106118 684353 9.293 3.38e-07 ***
## Residuals 60 4418454 73641
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Contrastes ortogonales no incluidosCódigos para SAS
/***********************************************************
EXPERIMENTACIÓN EN AGRICULTURA
Fernández
Escobar, R.; Trapero, A.; Domínguez, J. 2010
CAPÍTULO 13:
ANÁLISIS DE LA VARIANZA COMBINADO
Análisis de
diferentes localidades (páginas 184 - 190)
************************************************************/
DATA local;
INPUT LOCALIDAD $ EPOCA $ VARIEDAD $ BLOQUE RENDIMIENTO;
LABEL
LOCALIDAD = 'Localidades del Experimento'
EPOCA = 'Épocas de cultivo'
VARIEDAD = 'Variedades evaluadas'
BLOQUE = 'Repeticiones'
RENDIMIENTO = 'Rendimiento (kg/ha)';
DATALINES;
L1 invierno V1 1 0
L1 invierno V1 2 214
L1 invierno V1 3 425
L1 invierno V1 4 40
L1 invierno V2 1 1252
L1 invierno V2 2 627
L1 invierno V2 3 716
L1 invierno V2 4 1068
L1 invierno V3 1 2163
L1 invierno V3 2 2714
L1 invierno V3 3 2521
L1 invierno V3 4 2240
L1 primavera V1 1 1036
L1 primavera V1 2 697
L1 primavera V1 3 849
L1 primavera V1 4 1258
L1 primavera V2 1 1524
L1 primavera V2 2 1861
L1 primavera V2 3 2220
L1 primavera V2 4 1744
L1 primavera V3 1 1312
L1 primavera V3 2 874
L1 primavera V3 3 695
L1 primavera V3 4 1133
L2 invierno V1 1 26
L2 invierno V1 2 125
L2 invierno V1 3 0
L2 invierno V1 4 358
L2 invierno V2 1 813
L2 invierno V2 2 986
L2 invierno V2 3 657
L2 invierno V2 4 432
L2 invierno V3 1 2475
L2 invierno V3 2 2671
L2 invierno V3 3 3044
L2 invierno V3 4 2850
L2 primavera V1 1 488
L2 primavera V1 2 369
L2 primavera V1 3 714
L2 primavera V1 4 983
L2 primavera V2 1 1815
L2 primavera V2 2 1662
L2 primavera V2 3 1233
L2 primavera V2 4 1486
L2 primavera V3 1 1958
L2 primavera V3 2 1604
L2 primavera V3 3 1760
L2 primavera V3 4 1325
L3 invierno V1 1 257
L3 invierno V1 2 64
L3 invierno V1 3 756
L3 invierno V1 4 452
L3 invierno V2 1 789
L3 invierno V2 2 491
L3 invierno V2 3 608
L3 invierno V2 4 1066
L3 invierno V3 1 2801
L3 invierno V3 2 3125
L3 invierno V3 3 2977
L3 invierno V3 4 3453
L3 primavera V1 1 412
L3 primavera V1 2 1118
L3 primavera V1 3 833
L3 primavera V1 4 621
L3 primavera V2 1 1325
L3 primavera V2 2 1950
L3 primavera V2 3 1507
L3 primavera V2 4 1773
L3 primavera V3 1 1582
L3 primavera V3 2 1236
L3 primavera V3 3 1844
L3 primavera V3 4 1371
L4 invierno V1 1 1060
L4 invierno V1 2 1383
L4 invierno V1 3 889
L4 invierno V1 4 723
L4 invierno V2 1 816
L4 invierno V2 2 1125
L4 invierno V2 3 1463
L4 invierno V2 4 1277
L4 invierno V3 1 892
L4 invierno V3 2 1112
L4 invierno V3 3 1521
L4 invierno V3 4 1366
L4 primavera V1 1 960
L4 primavera V1 2 474
L4 primavera V1 3 899
L4 primavera V1 4 763
L4 primavera V2 1 725
L4 primavera V2 2 583
L4 primavera V2 3 1179
L4 primavera V2 4 914
L4 primavera V3 1 897
L4 primavera V3 2 985
L4 primavera V3 3 544
L4 primavera V3 4 718
;
ODS HTML;
* Imprimir datos originales;
PROC PRINT DATA =
local;
RUN;
* Estadísticas
descriptivas;
PROC MEANS MIN MAX MEDIAN MEAN RANGE Q1 Q3 QRANGE CV VAR
STD STDERR SUM KURT SKEW MAXDEC=2;
VAR
RENDIMIENTO;
CLASS LOCALIDAD EPOCA VARIEDAD;
RUN;
* Separando y analizando datos de 'Localidad 1';
title "ANOVA para
Localidad 1";
DATA local1;
SET
local(obs = 24); RUN;
PROC PRINT DATA =
local1;
RUN;
*Opción 1;
PROC GLM DATA =
local1;
CLASS BLOQUE EPOCA VARIEDAD;
MODEL RENDIMIENTO = BLOQUE EPOCA BLOQUE*EPOCA VARIEDAD
EPOCA*VARIEDAD / SS3;
RUN;
*Opción 2;
PROC MIXED METHOD =
TYPE3 DATA = local1;
CLASS BLOQUE EPOCA VARIEDAD;
MODEL RENDIMIENTO = BLOQUE EPOCA BLOQUE*EPOCA VARIEDAD
EPOCA*VARIEDAD/DDFM = SATTERTH;
RUN;
* Separando y analizando datos de 'Localidad 2';
title "ANOVA para
Localidad 2";
DATA local2;
SET
local(firstobs = 25 obs = 48); RUN;
PROC PRINT DATA =
local2;
RUN;
*Opción 1;
PROC GLM DATA =
local2;
CLASS BLOQUE EPOCA VARIEDAD;
MODEL RENDIMIENTO = BLOQUE EPOCA BLOQUE*EPOCA VARIEDAD
EPOCA*VARIEDAD / SS3;
RUN;
*Opción 2;
PROC MIXED METHOD =
TYPE3 DATA = local2;
CLASS BLOQUE EPOCA VARIEDAD;
MODEL RENDIMIENTO = BLOQUE EPOCA BLOQUE*EPOCA VARIEDAD
EPOCA*VARIEDAD/DDFM = SATTERTH;
RUN;
* Separando y analizando datos de 'Localidad 3';
title "ANOVA para
Localidad 3";
DATA local3;
SET
local(firstobs = 49 obs = 72); RUN;
PROC PRINT DATA =
local3;
RUN;
*Opción 1;
PROC GLM DATA =
local3;
CLASS BLOQUE EPOCA VARIEDAD;
MODEL RENDIMIENTO = BLOQUE EPOCA BLOQUE*EPOCA VARIEDAD
EPOCA*VARIEDAD / SS3;
RUN;
*Opción 2;
PROC MIXED METHOD =
TYPE3 DATA = local3;
CLASS BLOQUE EPOCA VARIEDAD;
MODEL RENDIMIENTO = BLOQUE EPOCA BLOQUE*EPOCA VARIEDAD
EPOCA*VARIEDAD/DDFM = SATTERTH;
RUN;
* Separando y analizando datos de 'Localidad 4';
title "ANOVA para
Localidad 4";
DATA local4;
SET
local(firstobs = 73); RUN;
PROC PRINT DATA =
local4;
RUN;
*Opción 1;
PROC GLM DATA =
local4;
CLASS BLOQUE EPOCA VARIEDAD;
MODEL RENDIMIENTO = BLOQUE EPOCA BLOQUE*EPOCA VARIEDAD
EPOCA*VARIEDAD / SS3;
RUN;
*Opción 2;
PROC MIXED METHOD =
TYPE3 DATA = local4;
CLASS BLOQUE EPOCA VARIEDAD;
MODEL RENDIMIENTO = BLOQUE EPOCA BLOQUE*EPOCA VARIEDAD
EPOCA*VARIEDAD/DDFM = SATTERTH;
RUN;
*ANOVA combinado;
title "ANOVA
combinado";
*Opción 1;
PROC GLM DATA =
local;
CLASS LOCALIDAD BLOQUE EPOCA VARIEDAD;
MODEL RENDIMIENTO = LOCALIDAD BLOQUE EPOCA
LOCALIDAD*BLOQUE LOCALIDAD*EPOCA VARIEDAD LOCALIDAD*VARIEDAD EPOCA*VARIEDAD
LOCALIDAD*EPOCA*VARIEDAD / SS3;
RUN;
*Opción 1;
PROC MIXED METHOD =
TYPE3 DATA = local;
CLASS LOCALIDAD BLOQUE EPOCA VARIEDAD;
MODEL RENDIMIENTO = LOCALIDAD BLOQUE EPOCA
LOCALIDAD*EPOCA VARIEDAD LOCALIDAD*VARIEDAD EPOCA*VARIEDAD
LOCALIDAD*EPOCA*VARIEDAD/DDFM = SATTERTH;
RUN;
*Opción 3;
PROC GLM DATA =
local;
CLASS LOCALIDAD BLOQUE EPOCA VARIEDAD;
MODEL RENDIMIENTO = LOCALIDAD | BLOQUE(LOCALIDAD) | EPOCA
| VARIEDAD /SS3 ;
TEST H=LOCALIDAD E = BLOQUE(LOCALIDAD);
TEST H=EPOCA E =
EPOCA*BLOQUE(LOCALIDAD);
TEST H=LOCALIDAD*EPOCA E =
EPOCA*BLOQUE(LOCALIDAD);
TEST H=VARIEDAD E = VARIEDAD*BLOQUE(LOCALIDAD);
TEST H=LOCALIDAD*VARIEDAD E = VARIEDAD*BLOQUE(LOCALIDAD);
TEST H=EPOCA*VARIEDAD E =
EPOCA*VARIEDAD*BLOQUE(LOCALIDAD);
TEST H=LOCALIDAD*EPOCA*VARIEDAD E =
EPOCA*VARIEDAD*BLOQUE(LOCALIDAD);
MEANS LOCALIDAD EPOCA VARIEDAD EPOCA*LOCALIDAD;
*Contrastes ortogonales no
incluidos;
RUN;
QUIT;
ODS HTML CLOSE;
-----------------------------------------------------
Archivos
adjuntos
Datos en
Excel : tabla 13.5.csv
Códigos
en R: capítulo 13e.R, capítulo 13eMg.r, summarySE.r, schwarz.functions.r
Códigos en SAS: capitulo 13e.sas
No hay comentarios:
Publicar un comentario