# ROGER GUEVARA (roger.guevara@inecol.edu.mx)
# INSTITUTO DE ECOLOGIA A.C RED DE BIOLOGÍA EVOLUTIVA
# ESTADISTICA COMPUTACIONAL EN R
# ANOVAS DE PARCELAS DIVIDIDAS, split-plot
# Vamos a utilizar las bases de datos << splityield,csv >>
DATOS <- read.table("~/desktop/curso R 2012/splityield.csv", header = TRUE,
sep = ",")
attach(DATOS)
names(DATOS)
## [1] "YIELD" "BLOCK" "IRRIGATION" "DENSITY" "FERTILIZER"
# Los modelos de ANOVA de parcelas divididas son modelos factoriales donde
# se permiten todas las interacciones entre los factores que explican,
# pero el error está ordenado jerárquicamente, es decir para cada tamaño
# de parcela existe un estimado diferente de error.
# La sintaxis para definir estos modelos lo haremos con la función <<
# aov() >> y noten el término de <<
# +Error(BLOCK/IRRIGATION/DENSITY/FERTILIZER) >> lo que define el error
# jerárquico
M <- aov(YIELD ~ IRRIGATION * DENSITY * FERTILIZER + Error(BLOCK/IRRIGATION/DENSITY/FERTILIZER))
summary(M)
##
## Error: BLOCK
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 3 194 64.8
##
## Error: BLOCK:IRRIGATION
## Df Sum Sq Mean Sq F value Pr(>F)
## IRRIGATION 1 8278 8278 17.6 0.025 *
## Residuals 3 1412 471
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: BLOCK:IRRIGATION:DENSITY
## Df Sum Sq Mean Sq F value Pr(>F)
## DENSITY 2 1758 879 3.78 0.053 .
## IRRIGATION:DENSITY 2 2747 1374 5.91 0.016 *
## Residuals 12 2788 232
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: BLOCK:IRRIGATION:DENSITY:FERTILIZER
## Df Sum Sq Mean Sq F value Pr(>F)
## FERTILIZER 2 1977 989 11.45 0.00014 ***
## IRRIGATION:FERTILIZER 2 953 477 5.52 0.00811 **
## DENSITY:FERTILIZER 4 305 76 0.88 0.48405
## IRRIGATION:DENSITY:FERTILIZER 4 235 59 0.68 0.61067
## Residuals 36 3109 86
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# El resumen de estos modelos produce una tabla de innova para cada tamaño
# de parcela que exista. Y como se mencionó cada tamaño de parcela tiene
# su propio estimado error
# A continuación se muestran los cálculos paso por paso para calcular las
# diversas sumas de cuadrados en el modelo, y como nota se muestran los
# grados de libertad. Para obtener las sumas de cuadrados de los factores
# jerárquicos inferiores es necesario restar la suma de los cuadrados de
# los factores superiores los cuales ya fueron considerados, lo mismo se
# hace con los grados de libertad. Además recuerda que para el calculo de
# los grados de libertad siempre se resta uno que corresponde al promedio
# general.
MG <- mean(YIELD)
SST <- sum((YIELD - MG)^2) #GL N - 1
SSB <- sum((tapply(YIELD, BLOCK, mean) - MG)^2 * table(BLOCK)) # 4-1=3
SSB
## [1] 194.4
SSI <- sum((tapply(YIELD, IRRIGATION, mean) - MG)^2 * table(IRRIGATION)) #2-1=1
SSI
## [1] 8278
SSBI <- sum((tapply(YIELD, list(BLOCK, IRRIGATION), mean) - MG)^2 * table(BLOCK,
IRRIGATION)) - SSB - SSI #8 - 1 - 3 -1 = 3
SSBI
## [1] 1412
SSD <- sum((tapply(YIELD, DENSITY, mean) - MG)^2 * table(DENSITY)) #3-1=2
SSD
## [1] 1758
SSID <- sum((tapply(YIELD, list(IRRIGATION, DENSITY), mean) - MG)^2 * table(IRRIGATION,
DENSITY)) - SSI - SSD #6 - 1 - 1 -2 = 2
SSID
## [1] 2747
SSBID <- sum((tapply(YIELD, list(BLOCK, IRRIGATION, DENSITY), mean) - MG)^2 *
table(BLOCK, IRRIGATION, DENSITY)) - SSB - SSI - SSBI - SSD - SSID #24 - 1 - 3 -1 -3 - 2 -2 = 12
SSBID
## [1] 2788
SSF <- sum((tapply(YIELD, FERTILIZER, mean) - MG)^2 * table(FERTILIZER)) #3-1=2
SSF
## [1] 1977
SSIF <- sum((tapply(YIELD, list(IRRIGATION, FERTILIZER), mean) - MG)^2 * table(IRRIGATION,
FERTILIZER)) - SSI - SSF #6 - 1 - 1 -2 = 2
SSIF
## [1] 953.4
SSDF <- sum((tapply(YIELD, list(DENSITY, FERTILIZER), mean) - MG)^2 * table(DENSITY,
FERTILIZER)) - SSD - SSF #9 - 1 - 2 -2 = 4
SSDF
## [1] 304.9
SSIDF <- sum((tapply(YIELD, list(IRRIGATION, DENSITY, FERTILIZER), mean) - MG)^2 *
table(IRRIGATION, DENSITY, FERTILIZER)) - SSI - SSD - SSID - SSF - SSIF -
SSDF #18 - 1 - 1 -2 -2 -2 -2 -4 = 4
SSIDF
## [1] 234.7
SSBIDF <- sum((tapply(YIELD, list(BLOCK, IRRIGATION, DENSITY, FERTILIZER), mean) -
MG)^2 * table(BLOCK, IRRIGATION, DENSITY, FERTILIZER)) - SSB - SSI - SSBI -
SSD - SSID - SSBID - SSF - SSIF - SSDF - SSIDF #72 - 1 - 3 -1 -3 -2 -2 -12 -2 -2 -4 -4 = 36
SSBIDF
## [1] 3109
# Otro ejemplo mucho más complicado. Ahora al nivel de los bloques existen
# tratamientos con insecticida y con molusquicida. Estos tratamientos son
# factoriales y no forman parte del estructura jerárquica de error, y pero
# adentro de estos bloques hay parcelas de distintos tamaños
DATOS <- read.table("~/desktop/curso R 2012/splitplot.csv", header = TRUE, sep = ",")
attach(DATOS)
## The following object(s) are masked from 'DATOS (position 3)':
##
## BLOCK
names(DATOS)
## [1] "BLOCK" "INSECT" "MOLLUSC" "RABBIT" "LIME"
## [6] "COMPETITION" "NUTRIENT" "BIOMASS"
model <- aov(BIOMASS ~ INSECT * MOLLUSC * RABBIT * LIME * COMPETITION * NUTRIENT +
Error(BLOCK/RABBIT/LIME/COMPETITION/NUTRIENT))
summary(model)
##
## Error: BLOCK
## Df Sum Sq Mean Sq F value Pr(>F)
## INSECT 1 414 414 34.25 0.0043 **
## MOLLUSC 1 9 9 0.73 0.4423
## INSECT:MOLLUSC 1 11 11 0.92 0.3924
## Residuals 4 48 12
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: BLOCK:RABBIT
## Df Sum Sq Mean Sq F value Pr(>F)
## RABBIT 1 389 389 4490.46 3e-07 ***
## INSECT:RABBIT 1 0 0 4.60 0.099 .
## MOLLUSC:RABBIT 1 0 0 0.15 0.714
## INSECT:MOLLUSC:RABBIT 1 0 0 2.88 0.165
## Residuals 4 0 0
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: BLOCK:RABBIT:LIME
## Df Sum Sq Mean Sq F value Pr(>F)
## LIME 1 86.7 86.7 1889.54 8.7e-11 ***
## INSECT:LIME 1 0.0 0.0 0.73 0.417
## MOLLUSC:LIME 1 0.1 0.1 2.63 0.143
## RABBIT:LIME 1 0.1 0.1 3.22 0.111
## INSECT:MOLLUSC:LIME 1 0.1 0.1 1.15 0.314
## INSECT:RABBIT:LIME 1 0.0 0.0 0.08 0.791
## MOLLUSC:RABBIT:LIME 1 0.1 0.1 1.98 0.197
## INSECT:MOLLUSC:RABBIT:LIME 1 0.5 0.5 10.18 0.013 *
## Residuals 8 0.4 0.0
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: BLOCK:RABBIT:LIME:COMPETITION
## Df Sum Sq Mean Sq F value Pr(>F)
## COMPETITION 2 214.5 107.2 1187.96 <2e-16
## INSECT:COMPETITION 2 0.1 0.1 0.81 0.453
## MOLLUSC:COMPETITION 2 0.2 0.1 0.86 0.432
## RABBIT:COMPETITION 2 0.2 0.1 1.13 0.337
## LIME:COMPETITION 2 0.0 0.0 0.14 0.874
## INSECT:MOLLUSC:COMPETITION 2 0.4 0.2 2.28 0.119
## INSECT:RABBIT:COMPETITION 2 0.1 0.1 0.69 0.509
## MOLLUSC:RABBIT:COMPETITION 2 0.0 0.0 0.13 0.880
## INSECT:LIME:COMPETITION 2 0.1 0.0 0.28 0.757
## MOLLUSC:LIME:COMPETITION 2 0.0 0.0 0.16 0.853
## RABBIT:LIME:COMPETITION 2 0.0 0.0 0.08 0.927
## INSECT:MOLLUSC:RABBIT:COMPETITION 2 0.0 0.0 0.17 0.845
## INSECT:MOLLUSC:LIME:COMPETITION 2 0.1 0.0 0.37 0.696
## INSECT:RABBIT:LIME:COMPETITION 2 0.4 0.2 2.08 0.142
## MOLLUSC:RABBIT:LIME:COMPETITION 2 0.5 0.3 2.77 0.077
## INSECT:MOLLUSC:RABBIT:LIME:COMPETITION 2 0.0 0.0 0.07 0.936
## Residuals 32 2.9 0.1
##
## COMPETITION ***
## INSECT:COMPETITION
## MOLLUSC:COMPETITION
## RABBIT:COMPETITION
## LIME:COMPETITION
## INSECT:MOLLUSC:COMPETITION
## INSECT:RABBIT:COMPETITION
## MOLLUSC:RABBIT:COMPETITION
## INSECT:LIME:COMPETITION
## MOLLUSC:LIME:COMPETITION
## RABBIT:LIME:COMPETITION
## INSECT:MOLLUSC:RABBIT:COMPETITION
## INSECT:MOLLUSC:LIME:COMPETITION
## INSECT:RABBIT:LIME:COMPETITION
## MOLLUSC:RABBIT:LIME:COMPETITION .
## INSECT:MOLLUSC:RABBIT:LIME:COMPETITION
## Residuals
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: BLOCK:RABBIT:LIME:COMPETITION:NUTRIENT
## Df Sum Sq Mean Sq F value
## NUTRIENT 3 1426 475 6997.49
## INSECT:NUTRIENT 3 0 0 1.03
## MOLLUSC:NUTRIENT 3 0 0 0.60
## RABBIT:NUTRIENT 3 0 0 0.43
## LIME:NUTRIENT 3 0 0 0.17
## COMPETITION:NUTRIENT 6 1 0 1.81
## INSECT:MOLLUSC:NUTRIENT 3 0 0 0.56
## INSECT:RABBIT:NUTRIENT 3 1 0 3.89
## MOLLUSC:RABBIT:NUTRIENT 3 1 0 4.52
## INSECT:LIME:NUTRIENT 3 0 0 0.28
## MOLLUSC:LIME:NUTRIENT 3 0 0 2.34
## RABBIT:LIME:NUTRIENT 3 0 0 1.85
## INSECT:COMPETITION:NUTRIENT 6 1 0 1.36
## MOLLUSC:COMPETITION:NUTRIENT 6 0 0 0.84
## RABBIT:COMPETITION:NUTRIENT 6 0 0 1.07
## LIME:COMPETITION:NUTRIENT 6 0 0 1.11
## INSECT:MOLLUSC:RABBIT:NUTRIENT 3 0 0 2.02
## INSECT:MOLLUSC:LIME:NUTRIENT 3 1 0 3.41
## INSECT:RABBIT:LIME:NUTRIENT 3 0 0 2.15
## MOLLUSC:RABBIT:LIME:NUTRIENT 3 0 0 0.26
## INSECT:MOLLUSC:COMPETITION:NUTRIENT 6 1 0 1.38
## INSECT:RABBIT:COMPETITION:NUTRIENT 6 1 0 1.78
## MOLLUSC:RABBIT:COMPETITION:NUTRIENT 6 0 0 1.06
## INSECT:LIME:COMPETITION:NUTRIENT 6 0 0 0.64
## MOLLUSC:LIME:COMPETITION:NUTRIENT 6 0 0 1.02
## RABBIT:LIME:COMPETITION:NUTRIENT 6 0 0 0.99
## INSECT:MOLLUSC:RABBIT:LIME:NUTRIENT 3 0 0 1.72
## INSECT:MOLLUSC:RABBIT:COMPETITION:NUTRIENT 6 0 0 0.64
## INSECT:MOLLUSC:LIME:COMPETITION:NUTRIENT 6 0 0 1.01
## INSECT:RABBIT:LIME:COMPETITION:NUTRIENT 6 0 0 0.69
## MOLLUSC:RABBIT:LIME:COMPETITION:NUTRIENT 6 0 0 0.88
## INSECT:MOLLUSC:RABBIT:LIME:COMPETITION:NUTRIENT 6 1 0 2.41
## Residuals 144 10 0
## Pr(>F)
## NUTRIENT <2e-16 ***
## INSECT:NUTRIENT 0.3826
## MOLLUSC:NUTRIENT 0.6182
## RABBIT:NUTRIENT 0.7333
## LIME:NUTRIENT 0.9141
## COMPETITION:NUTRIENT 0.1011
## INSECT:MOLLUSC:NUTRIENT 0.6452
## INSECT:RABBIT:NUTRIENT 0.0104 *
## MOLLUSC:RABBIT:NUTRIENT 0.0046 **
## INSECT:LIME:NUTRIENT 0.8406
## MOLLUSC:LIME:NUTRIENT 0.0763 .
## RABBIT:LIME:NUTRIENT 0.1401
## INSECT:COMPETITION:NUTRIENT 0.2348
## MOLLUSC:COMPETITION:NUTRIENT 0.5385
## RABBIT:COMPETITION:NUTRIENT 0.3863
## LIME:COMPETITION:NUTRIENT 0.3621
## INSECT:MOLLUSC:RABBIT:NUTRIENT 0.1135
## INSECT:MOLLUSC:LIME:NUTRIENT 0.0193 *
## INSECT:RABBIT:LIME:NUTRIENT 0.0964 .
## MOLLUSC:RABBIT:LIME:NUTRIENT 0.8521
## INSECT:MOLLUSC:COMPETITION:NUTRIENT 0.2248
## INSECT:RABBIT:COMPETITION:NUTRIENT 0.1067
## MOLLUSC:RABBIT:COMPETITION:NUTRIENT 0.3896
## INSECT:LIME:COMPETITION:NUTRIENT 0.6991
## MOLLUSC:LIME:COMPETITION:NUTRIENT 0.4159
## RABBIT:LIME:COMPETITION:NUTRIENT 0.4338
## INSECT:MOLLUSC:RABBIT:LIME:NUTRIENT 0.1655
## INSECT:MOLLUSC:RABBIT:COMPETITION:NUTRIENT 0.6995
## INSECT:MOLLUSC:LIME:COMPETITION:NUTRIENT 0.4228
## INSECT:RABBIT:LIME:COMPETITION:NUTRIENT 0.6552
## MOLLUSC:RABBIT:LIME:COMPETITION:NUTRIENT 0.5119
## INSECT:MOLLUSC:RABBIT:LIME:COMPETITION:NUTRIENT 0.0298 *
## Residuals
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
MG <- mean(BIOMASS)
SSI <- sum((tapply(BIOMASS, INSECT, mean) - MG)^2 * table(INSECT)) #2 - 1 = 1
SSM <- sum((tapply(BIOMASS, MOLLUSC, mean) - MG)^2 * table(MOLLUSC)) #2 - 1 = 1
SSIM <- sum((tapply(BIOMASS, list(INSECT, MOLLUSC), mean) - MG)^2 * table(INSECT,
MOLLUSC)) - SSI - SSM # 4 - 1 -1 -1
SSB <- sum((tapply(BIOMASS, BLOCK, mean) - MG)^2 * table(BLOCK)) - SSI - SSM -
SSIM #8 -1 -1 -1 -1 = 4
SSR <- sum((tapply(BIOMASS, RABBIT, mean) - MG)^2 * table(RABBIT)) #2 - 1 = 1
SSIR <- sum((tapply(BIOMASS, list(INSECT, RABBIT), mean) - MG)^2 * table(INSECT,
RABBIT)) - SSI - SSR #4 - 1 -1 -1= 1
SSMR <- sum((tapply(BIOMASS, list(MOLLUSC, RABBIT), mean) - MG)^2 * table(MOLLUSC,
RABBIT)) - SSM - SSR # 4 - 1 -1 -1 = 1
SSIMR <- sum((tapply(BIOMASS, list(INSECT, MOLLUSC, RABBIT), mean) - MG)^2 *
table(INSECT, MOLLUSC, RABBIT)) - SSI - SSM - SSIM - SSR - SSIR - SSMR #8 -1 -1 -1 -1 - 1 -1 -1 = 1
SSBR <- sum((tapply(BIOMASS, list(BLOCK, RABBIT), mean) - MG)^2 * table(BLOCK,
RABBIT)) - SSI - SSM - SSIM - SSB - SSR - SSIR - SSMR - SSIMR #16 -1 -1 -1 -1 -4 -1 -1 -1 -1 = 4
SSL <- sum((tapply(BIOMASS, LIME, mean) - MG)^2 * table(LIME)) #2 - 1 = 1
SSL
## [1] 86.67
SSIL <- sum((tapply(BIOMASS, list(INSECT, LIME), mean) - MG)^2 * table(INSECT,
LIME)) - SSI - SSL #4 -1 -1 -1= 1
SSIL
## [1] 0.03356
SSML <- sum((tapply(BIOMASS, list(MOLLUSC, LIME), mean) - MG)^2 * table(MOLLUSC,
LIME)) - SSM - SSL #4 - 1 -1 -1 = 1
SSML
## [1] 0.1208
SSRL <- sum((tapply(BIOMASS, list(RABBIT, LIME), mean) - MG)^2 * table(RABBIT,
LIME)) - SSR - SSL #4 -1 -1 -1 =1
SSRL
## [1] 0.1477
SSIML <- sum((tapply(BIOMASS, list(INSECT, MOLLUSC, LIME), mean) - MG)^2 * table(INSECT,
MOLLUSC, LIME)) - SSI - SSM - SSIM - SSL - SSIL - SSML #8 -1 -1 -1 -1 -1 -1 -1 = 1
SSIML
## [1] 0.05297
SSIRL <- sum((tapply(BIOMASS, list(INSECT, RABBIT, LIME), mean) - MG)^2 * table(INSECT,
RABBIT, LIME)) - SSI - SSR - SSIR - SSL - SSIL - SSRL #8 -1 -1 -1 -1 -1 -1 -1 = 1
SSIRL
## [1] 0.003444
SSMRL <- sum((tapply(BIOMASS, list(MOLLUSC, RABBIT, LIME), mean) - MG)^2 * table(MOLLUSC,
RABBIT, LIME)) - SSM - SSR - SSMR - SSL - SSML - SSRL #8 -1 -1 -1 -1 -1 -1 -1 = 1
SSMRL
## [1] 0.09096
SSIMRL <- sum((tapply(BIOMASS, list(INSECT, MOLLUSC, RABBIT, LIME), mean) -
MG)^2 * table(INSECT, MOLLUSC, RABBIT, LIME)) - SSI - SSM - SSIM - SSR -
SSIR - SSMR - SSIMR - SSL - SSIL - SSML - SSRL - SSIML - SSIRL - SSMRL #16 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1
SSIMRL
## [1] 0.4669
SSBRL <- sum((tapply(BIOMASS, list(BLOCK, RABBIT, LIME), mean) - MG)^2 * table(BLOCK,
RABBIT, LIME)) - SSI - SSM - SSIM - SSB - SSR - SSIR - SSMR - SSIMR - SSBR -
SSL - SSIL - SSML - SSRL - SSIML - SSIRL - SSMRL - SSIMRL #32 -1 -1 -1 -1 -4 -1 -1 -1 -1 -4 -1 -1 -1 -1 -1 -1 -1 -1= 8
SSBRL
## [1] 0.3669
# Aún faltan dos tablas más correspondientes a los niveles jerárquicos mas
# internos, las parcelas mas pequeñas. Aquí vamos a dejar el ejercicio
# pero si tienen curiosidad y tiempovaldría la pena que completaran estas
# dos tablas faltates para que todo este procedimiento quede mas claro