abril 04, 2017

ANOVA combinado (4) — diferentes años


Creamos múltiples gráficos para visualizar desde diferentes perspectivas la tendencia de los tratamientos en los años de muestreo.


#### Gráfico de interacción + errores típicos
pd <- position_dodge(0.1) # mover barras de errores .05 hacia la izquierda y la derecha
SEgraph <- ggplot(data = MEDIAS, aes(x = TRAT, y = AFECT,
                                     group    = FECHA, 
                                     color    = FECHA, 
                                     linetype = FECHA, 
                                     shape    = FECHA))                                  + 
  geom_errorbar(aes(ymin =  AFECT-se, ymax    = AFECT+se), 
                width    =  .2,      position = pd,  size     = 0.4, linetype = 'solid') +
  geom_line(position     =  pd,      size     = 1.0, linetype = 'solid')                 + 
  geom_point(position    =  pd,      size     = 3.0)                                     +     
  scale_color_brewer(palette = "Set2")                                                   + 
  theme_minimal()                                                                        + 
  ggtitle("Gráfico con errores típicos")                                                 +
  xlab("Tratamientos")                                                                   + 
  ylab("Tejido foliar afectado (Arcsen)"); SEgraph


“La variabilidad entre años es generalmente impredecible, por lo que los años se consideran como una variable aleatoria.” Análisis de diferentes años (pp. 180 – 184).



#### Representando "TRAT" en abscisas
# Gráfico de barras con errores típicos
ggplot(MEDIAS, aes(x        =  TRAT, 
                   y        =  AFECT, 
                   fill     =  FECHA))                                 + 
       geom_bar(position    =  position_dodge(), stat = "identity")    + 
  geom_errorbar(aes(ymin    =  AFECT-se, ymax = AFECT+se), width = .2, 
                position    =  position_dodge(.9))                     +   
  scale_fill_brewer(palette =  "Set2") + theme_minimal()


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 años (páginas 180 - 184)
#************************************************************

#### Importar Tabla 13.4
foliarDOBLE <- read.csv("C:/Users/Administrator/Desktop/tabla 13.4.csv", header=TRUE)
attach(foliarDOBLE)
print(foliarDOBLE)
##    FECHA REP      TRAT AFECT
## 1      1   I SinTratar  61.9
## 2      1  II SinTratar  55.9
## 3      1 III SinTratar  64.5
## 4      1  IV SinTratar  57.7
## 5      1   V SinTratar  50.9
## 6      1   I  Invierno  20.3
## 7      1  II  Invierno  16.5
## 8      1 III  Invierno  15.0
## 9      1  IV  Invierno   8.5
## 10     1   V  Invierno  23.8
## 11     1   I Primavera  54.0
## 12     1  II Primavera  60.4
## 13     1 III Primavera  49.7
## 14     1  IV Primavera  70.4
## 15     1   V Primavera  49.0
## 16     1   I  Inv+Prim  14.4
## 17     1  II  Inv+Prim  16.1
## 18     1 III  Inv+Prim  13.4
## 19     1  IV  Inv+Prim  10.9
## 20     1   V  Inv+Prim  22.8
## 21     2   I SinTratar  70.9
## 22     2  II SinTratar  61.2
## 23     2 III SinTratar  56.8
## 24     2  IV SinTratar  62.4
## 25     2   V SinTratar  64.5
## 26     2   I  Invierno  15.2
## 27     2  II  Invierno  29.9
## 28     2 III  Invierno  18.4
## 29     2  IV  Invierno  22.1
## 30     2   V  Invierno  17.7
## 31     2   I Primavera  46.3
## 32     2  II Primavera  32.0
## 33     2 III Primavera  40.1
## 34     2  IV Primavera  56.0
## 35     2   V Primavera  38.4
## 36     2   I  Inv+Prim   4.8
## 37     2  II  Inv+Prim  14.5
## 38     2 III  Inv+Prim  19.6
## 39     2  IV  Inv+Prim   9.8
## 40     2   V  Inv+Prim  12.2
#### Creando variables independientes como factores
FECHA   <- factor(FECHA, levels = c(1,2), labels = c("Año 1", "Año 2")) # "Años" en libro
REP    <- factor(REP) 
TRAT   <- factor(TRAT)
fDOBLE <- data.frame(FECHA, REP, TRAT, AFECT)
attach(fDOBLE); detach(foliarDOBLE); str(foliarDOBLE)
## 'data.frame':    40 obs. of  4 variables:
##  $ FECHA: int  1 1 1 1 1 1 1 1 1 1 ...
##  $ REP  : Factor w/ 5 levels "I","II","III",..: 1 2 3 4 5 1 2 3 4 5 ...
##  $ TRAT : Factor w/ 4 levels "Inv+Prim","Invierno",..: 4 4 4 4 4 2 2 2 2 2 ...
##  $ AFECT: num  61.9 55.9 64.5 57.7 50.9 20.3 16.5 15 8.5 23.8 ...
#### Agregando función para calcular estadísticas descriptivas de las interacciones
library(plyr)
source("C:/Users/Administrator/Desktop/summarySE.r") 
MEDIAS <- summarySE(fDOBLE, measurevar="AFECT", groupvars=c("FECHA", "TRAT"))
MEDIAS
##   FECHA      TRAT N AFECT       sd       se        ci
## 1 Año 1  Inv+Prim 5 15.52 4.484083 2.005343  5.567724
## 2 Año 1  Invierno 5 16.82 5.776418 2.583292  7.172370
## 3 Año 1 Primavera 5 56.70 8.901685 3.980955 11.052902
## 4 Año 1 SinTratar 5 58.18 5.296414 2.368628  6.576366
## 5 Año 2  Inv+Prim 5 12.18 5.490173 2.455280  6.816950
## 6 Año 2  Invierno 5 20.66 5.725644 2.560586  7.109326
## 7 Año 2 Primavera 5 42.56 9.075957 4.058891 11.269289
## 8 Año 2 SinTratar 5 63.16 5.161686 2.308376  6.409079
#### Gráficos
#### Representando "Años" en abscisas
# Gráfico de barras con errores típicos
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 3.3.3
ggplot(MEDIAS, aes(x =  FECHA, 
                   y =  AFECT, 
                   fill  = TRAT)) + 
  geom_bar(position      =  position_dodge(), stat = "identity") + 
  geom_errorbar(aes(ymin =  AFECT-se, 
                    ymax =  AFECT+se), 
                   width =  .2, 
                position =  position_dodge(.9))  +   
  scale_fill_brewer(palette=  "Set2") + theme_minimal()
#### Diagrama de dispersión
ggplot(fDOBLE, aes(x=FECHA, y=AFECT, color=TRAT)) + geom_point(shape=1) + 
  scale_color_brewer(palette = "Set2") + geom_point(size = 3.0)
#### Gráfico de interacción simple
interaction.plot(FECHA, TRAT, AFECT)
#### Gráfico de interacción coloreado
ggplot(MEDIAS, aes(x = factor(FECHA), y = AFECT, colour = TRAT, group = TRAT)) +
  scale_color_brewer(palette = "Set2") + geom_line(size = 1.0, linetype = 'solid')
#### Gráfico de interacción + errores típicos
pd <- position_dodge(0.1) # mover barras de errores .05 hacia la izquierda y la derecha
SEgraph <- ggplot(data = MEDIAS, aes(x = FECHA, y = AFECT,
                                     group    = TRAT, 
                                     color    = TRAT, 
                                     linetype = TRAT, 
                                     shape    = TRAT))                                   + 
  geom_errorbar(aes(ymin =  AFECT-se, ymax    = AFECT+se), 
                width    =  .2,      position = pd,  size     = 0.4, linetype = 'solid') +
  geom_line(position     =  pd,      size     = 1.0, linetype = 'solid')                 + 
  geom_point(position    =  pd,      size     = 3.0)                                     +     
  scale_color_brewer(palette = "Set2")                                                   + 
  theme_minimal()                                                                        + 
  ggtitle("Gráfico con errores típicos")                                                 +
  xlab("Años")                                                                           + 
  ylab("Tejido foliar afectado (Arcsen)"); SEgraph
#### Boxplot y líneas de interacción
library("plyr")
ipMEDIAS <- ddply(fDOBLE,.(FECHA,TRAT), summarise, val  = mean(AFECT)) # precalculando medias
ggplot(fDOBLE, aes(x = factor(FECHA), y = AFECT, colour = TRAT))                                     + 
  geom_boxplot()                                                                                     + 
  geom_point(data = ipMEDIAS, aes(y = val), position = position_dodge(width=0.75))                   +
  geom_line(data  = ipMEDIAS, aes(y = val, group     = TRAT), position = position_dodge(width=0.75)) + 
  scale_color_brewer(palette = "Set2")                                                               +
  scale_x_discrete("Años")                                                                           + 
  scale_y_continuous("Tejido foliar afectado (Arcsen)")                                + 
  theme_bw() + theme_bw()                                                                            +
  theme(axis.title.x = element_text(size = 12, hjust = 0.54, vjust = 0),
        axis.title.y = element_text(size = 12, angle = 90,  vjust = 0.25))
#### Boxplot para interacciones
library(plotly)
# FECHA:TRAT
gfDOBLE <- plot_ly(fDOBLE, y = ~AFECT, color = ~FECHA:TRAT,  type = "box", 
                    boxpoints = "all", jitter = 0.3, pointpos = -1.8); gfDOBLE
Año 1:Inv+PrimAño 1:InviernoAño 1:PrimaveraAño 1:SinTratarAño 2:Inv+PrimAño 2:InviernoAño 2:PrimaveraAño 2:SinTratar10203040506070
AFECTAño 1:Inv+PrimAño 1:InviernoAño 1:PrimaveraAño 1:SinTratarAño 2:Inv+PrimAño 2:InviernoAño 2:PrimaveraAño 2:SinTratar
# TRAT:FECHA
gfDOBLE  <- plot_ly(fDOBLE, y = ~AFECT, color = ~TRAT:FECHA,  type = "box", 
                    boxpoints = "all", jitter = 0.3, pointpos = -1.8); gfDOBLE
Inv+Prim:Año 1Inv+Prim:Año 2Invierno:Año 1Invierno:Año 2Primavera:Año 1Primavera:Año 2SinTratar:Año 1SinTratar:Año 210203040506070
AFECTInv+Prim:Año 1Inv+Prim:Año 2Invierno:Año 1Invierno:Año 2Primavera:Año 1Primavera:Año 2SinTratar:Año 1SinTratar:Año 2
#### Representando "TRAT" en abscisas
# Gráfico de barras con errores típicos
ggplot(MEDIAS, aes(x        =  TRAT, 
                   y        =  AFECT, 
                   fill     =  FECHA))                                 + 
       geom_bar(position    =  position_dodge(), stat = "identity")    + 
  geom_errorbar(aes(ymin    =  AFECT-se, ymax = AFECT+se), width = .2, 
                position    =  position_dodge(.9))                     +   
  scale_fill_brewer(palette =  "Set2") + theme_minimal()
#### Gráfico de interacción simple
interaction.plot(TRAT, FECHA, AFECT)
#### Gráfico de interacción coloreado
ggplot(MEDIAS, aes(x = factor(TRAT), y = AFECT, colour = FECHA, group = FECHA)) +
  scale_color_brewer(palette = "Set2") + geom_line(size = 1.0, linetype = 'solid')
#### Gráfico de interacción + errores típicos
pd <- position_dodge(0.1) # mover barras de errores .05 hacia la izquierda y la derecha
SEgraph <- ggplot(data = MEDIAS, aes(x = TRAT, y = AFECT,
                                     group    = FECHA, 
                                     color    = FECHA, 
                                     linetype = FECHA, 
                                     shape    = FECHA))                                  + 
  geom_errorbar(aes(ymin =  AFECT-se, ymax    = AFECT+se), 
                width    =  .2,      position = pd,  size     = 0.4, linetype = 'solid') +
  geom_line(position     =  pd,      size     = 1.0, linetype = 'solid')                 + 
  geom_point(position    =  pd,      size     = 3.0)                                     +     
  scale_color_brewer(palette = "Set2")                                                   + 
  theme_minimal()                                                                        + 
  ggtitle("Gráfico con errores típicos")                                                 +
  xlab("Tratamientos")                                                                   + 
  ylab("Tejido foliar afectado (Arcsen)"); SEgraph
#### Boxplot y líneas de interacción
library("plyr")
ipMEDIAS <- ddply(fDOBLE,.(FECHA,TRAT), summarise, val  = mean(AFECT)) # precalculando medias
ggplot(fDOBLE, aes(x = factor(TRAT), y = AFECT, colour = FECHA))                                     + 
  geom_boxplot()                                                                                     + 
  geom_point(data = ipMEDIAS, aes(y = val), position = position_dodge(width=0.75))                   +
  geom_line(data  = ipMEDIAS, aes(y = val, group     = FECHA), position = position_dodge(width=0.75))+ 
  scale_color_brewer(palette = "Set2")                                                               +
  scale_x_discrete("Tratamientos")                                                                   + 
  scale_y_continuous("Tejido foliar afectado (Arcsen)")                                              + 
  theme_bw() + theme_bw()                                                                            +
  theme(axis.title.x = element_text(size = 12, hjust = 0.54, vjust = 0),
        axis.title.y = element_text(size = 12, angle = 90,  vjust = 0.25))
#### Boxplot para interacciones
library(plotly)
# TRAT:FECHA 
gfDOBLE <- plot_ly(fDOBLE, y = ~AFECT, color = ~TRAT:FECHA,  type = "box", 
                   boxpoints = "all", jitter = 0.3, pointpos = -1.8); gfDOBLE
Inv+Prim:Año 1Inv+Prim:Año 2Invierno:Año 1Invierno:Año 2Primavera:Año 1Primavera:Año 2SinTratar:Año 1SinTratar:Año 210203040506070
AFECTInv+Prim:Año 1Inv+Prim:Año 2Invierno:Año 1Invierno:Año 2Primavera:Año 1Primavera:Año 2SinTratar:Año 1SinTratar:Año 2
# FECHA:TRAT
gfDOBLE  <- plot_ly(fDOBLE, y = ~AFECT, color = ~FECHA:TRAT,  type = "box", 
                    boxpoints = "all", jitter = 0.3, pointpos = -1.8); gfDOBLE
Año 1:Inv+PrimAño 1:InviernoAño 1:PrimaveraAño 1:SinTratarAño 2:Inv+PrimAño 2:InviernoAño 2:PrimaveraAño 2:SinTratar10203040506070
AFECTAño 1:Inv+PrimAño 1:InviernoAño 1:PrimaveraAño 1:SinTratarAño 2:Inv+PrimAño 2:InviernoAño 2:PrimaveraAño 2:SinTratar
#### ANOVA para el AFECT en Año 1
# Separando datos de Año 1
foliarUNO <- fDOBLE[1:20, ]; foliarUNO
##    FECHA REP      TRAT AFECT
## 1  Año 1   I SinTratar  61.9
## 2  Año 1  II SinTratar  55.9
## 3  Año 1 III SinTratar  64.5
## 4  Año 1  IV SinTratar  57.7
## 5  Año 1   V SinTratar  50.9
## 6  Año 1   I  Invierno  20.3
## 7  Año 1  II  Invierno  16.5
## 8  Año 1 III  Invierno  15.0
## 9  Año 1  IV  Invierno   8.5
## 10 Año 1   V  Invierno  23.8
## 11 Año 1   I Primavera  54.0
## 12 Año 1  II Primavera  60.4
## 13 Año 1 III Primavera  49.7
## 14 Año 1  IV Primavera  70.4
## 15 Año 1   V Primavera  49.0
## 16 Año 1   I  Inv+Prim  14.4
## 17 Año 1  II  Inv+Prim  16.1
## 18 Año 1 III  Inv+Prim  13.4
## 19 Año 1  IV  Inv+Prim  10.9
## 20 Año 1   V  Inv+Prim  22.8
# Estadísticas descriptivas para TRAT
library(doBy) 
library(pastecs)
summaryBy(AFECT ~ TRAT, data = foliarUNO, FUN = stat.desc)
##        TRAT AFECT.nbr.val AFECT.nbr.null AFECT.nbr.na AFECT.min AFECT.max
## 1  Inv+Prim             5              0            0      10.9      22.8
## 2  Invierno             5              0            0       8.5      23.8
## 3 Primavera             5              0            0      49.0      70.4
## 4 SinTratar             5              0            0      50.9      64.5
##   AFECT.range AFECT.sum AFECT.median AFECT.mean AFECT.SE.mean
## 1        11.9      77.6         14.4      15.52      2.005343
## 2        15.3      84.1         16.5      16.82      2.583292
## 3        21.4     283.5         54.0      56.70      3.980955
## 4        13.6     290.9         57.7      58.18      2.368628
##   AFECT.CI.mean.0.95 AFECT.var AFECT.std.dev AFECT.coef.var
## 1           5.567724    20.107      4.484083     0.28892287
## 2           7.172370    33.367      5.776418     0.34342554
## 3          11.052902    79.240      8.901685     0.15699621
## 4           6.576366    28.052      5.296414     0.09103496
# Gráfico de Medias con 95% de intervalo de confianza
library(gplots)
plotmeans(AFECT ~ TRAT,xlab="Tratamientos", data = foliarUNO,
          ylab="Hojas afectados (Arsen)", main="Gráfico de Medias\ncon 95% CI")
# Boxplot para AFECT de TRAT en Año 1
gfoliarUNO <- plot_ly(foliarUNO, y = ~AFECT, color = ~TRAT,  type = "box", 
                     boxpoints = "all", jitter = 0.3, pointpos = -1.8); gfoliarUNO
Inv+PrimInviernoPrimaveraSinTratar10203040506070
AFECTInv+PrimInviernoPrimaveraSinTratar
# ANOVA Diseño en Bloques al Azar
ANOVAfUNO <- aov(AFECT ~ REP + TRAT, data = foliarUNO) 
summary(ANOVAfUNO)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## REP          4      9     2.3   0.043    0.996    
## TRAT         3   8526  2841.9  53.789 3.13e-07 ***
## Residuals   12    634    52.8                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
confint(ANOVAfUNO)
##                    2.5 %    97.5 %
## (Intercept)     6.348644 26.381356
## REPII         -11.623626 10.773626
## REPIII        -13.198626  9.198626
## REPIV         -11.973626 10.423626
## REPV          -12.223626 10.173626
## TRATInvierno   -8.716356 11.316356
## TRATPrimavera  31.163644 51.196356
## TRATSinTratar  32.643644 52.676356
plot(ANOVAfUNO)
# Gráfico de residuos contra valores ajustados
ggplot(data.frame(Fitted = fitted(ANOVAfUNO), Residuals = resid(ANOVAfUNO), Treatment = foliarUNO$TRAT), 
       aes(Fitted, Residuals, colour = Treatment)) + geom_point()
#### ANOVA para el AFECT en Año 2
# Separando datos de Año 2
foliarDOS <- fDOBLE[21:40, ]; foliarDOS
##    FECHA REP      TRAT AFECT
## 21 Año 2   I SinTratar  70.9
## 22 Año 2  II SinTratar  61.2
## 23 Año 2 III SinTratar  56.8
## 24 Año 2  IV SinTratar  62.4
## 25 Año 2   V SinTratar  64.5
## 26 Año 2   I  Invierno  15.2
## 27 Año 2  II  Invierno  29.9
## 28 Año 2 III  Invierno  18.4
## 29 Año 2  IV  Invierno  22.1
## 30 Año 2   V  Invierno  17.7
## 31 Año 2   I Primavera  46.3
## 32 Año 2  II Primavera  32.0
## 33 Año 2 III Primavera  40.1
## 34 Año 2  IV Primavera  56.0
## 35 Año 2   V Primavera  38.4
## 36 Año 2   I  Inv+Prim   4.8
## 37 Año 2  II  Inv+Prim  14.5
## 38 Año 2 III  Inv+Prim  19.6
## 39 Año 2  IV  Inv+Prim   9.8
## 40 Año 2   V  Inv+Prim  12.2
# Estadísticas descriptivas para Año 2
summaryBy(AFECT ~ TRAT, data = foliarDOS, FUN = stat.desc)
##        TRAT AFECT.nbr.val AFECT.nbr.null AFECT.nbr.na AFECT.min AFECT.max
## 1  Inv+Prim             5              0            0       4.8      19.6
## 2  Invierno             5              0            0      15.2      29.9
## 3 Primavera             5              0            0      32.0      56.0
## 4 SinTratar             5              0            0      56.8      70.9
##   AFECT.range AFECT.sum AFECT.median AFECT.mean AFECT.SE.mean
## 1        14.8      60.9         12.2      12.18      2.455280
## 2        14.7     103.3         18.4      20.66      2.560586
## 3        24.0     212.8         40.1      42.56      4.058891
## 4        14.1     315.8         62.4      63.16      2.308376
##   AFECT.CI.mean.0.95 AFECT.var AFECT.std.dev AFECT.coef.var
## 1           6.816950    30.142      5.490173     0.45075312
## 2           7.109326    32.783      5.725644     0.27713669
## 3          11.269289    82.373      9.075957     0.21325088
## 4           6.409079    26.643      5.161686     0.08172397
# Gráfico de Medias con 95% de intervalo de confianza
plotmeans(AFECT ~ TRAT,xlab="Tratamientos", data = foliarDOS,
          ylab="Hojas afectados (Arcsen)", main="Gráfico de Medias\ncon 95% CI")
# Boxplot para AFECT de TRATes en Año 2
gfDOS <- plot_ly(foliarDOS,   y = ~AFECT, color = ~TRAT,    type = "box", 
                      boxpoints = "all", jitter = 0.3, pointpos  = -1.8); gfDOS
Inv+PrimInviernoPrimaveraSinTratar10203040506070
AFECTInv+PrimInviernoPrimaveraSinTratar
# ANOVA Diseño en Bloques al Azar
ANOVAfDOS <- aov(AFECT ~ REP + TRAT, data = foliarDOS) 
summary(ANOVAfDOS)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## REP          4     47    11.7   0.219    0.923    
## TRAT         3   7880  2626.7  49.176 5.14e-07 ***
## Residuals   12    641    53.4                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
confint(ANOVAfDOS)
##                    2.5 %   97.5 %
## (Intercept)     1.768863 21.91114
## REPII         -11.159873 11.35987
## REPIII        -11.834873 10.68487
## REPIV          -7.984873 14.53487
## REPV          -12.359873 10.15987
## TRATInvierno   -1.591137 18.55114
## TRATPrimavera  20.308863 40.45114
## TRATSinTratar  40.908863 61.05114
plot(ANOVAfDOS)
# Gráfico de residuos contra valores ajustados
ggplot(data.frame(Fitted = fitted(ANOVAfDOS), Residuals = resid(ANOVAfDOS), Treatment = foliarDOS$TRAT), 
       aes(Fitted, Residuals, colour = Treatment)) + geom_point()
#### Modelo FECHA FECHA*REP(E) TRAT TRAT*FECHA
# Opción 1
MODELO <- aov(AFECT ~ FECHA + Error(FECHA/REP) + FECHA*TRAT, fDOBLE)
summary(MODELO)
## 
## Error: FECHA
##       Df Sum Sq Mean Sq
## FECHA  1  46.87   46.87
## 
## Error: FECHA:REP
##           Df Sum Sq Mean Sq F value Pr(>F)
## Residuals  8  55.84    6.98               
## 
## Error: Within
##            Df Sum Sq Mean Sq F value   Pr(>F)    
## TRAT        3  15826    5275  99.302 1.15e-13 ***
## FECHA:TRAT  3    580     193   3.638   0.0271 *  
## Residuals  24   1275      53                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Opción 2
MODELO <- aov(AFECT ~ FECHA + FECHA/REP + FECHA*TRAT, fDOBLE)
summary(MODELO)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## FECHA        1     47      47   0.882   0.3569    
## TRAT         3  15826    5275  99.302 1.15e-13 ***
## FECHA:REP    8     56       7   0.131   0.9971    
## FECHA:TRAT   3    580     193   3.638   0.0271 *  
## Residuals   24   1275      53                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#### Descomposición de las sumas de cuadrados
#### de ambos factores en CONTRASTES ORTOGONALES,
#### con un grado de libertad

# Verificar niveles y orden de tratamientos
str(fDOBLE)
# Coeficientes de contraste
# Contraste      IxP   I    P     ST
# Invierno        1    1   -1    -1
# Primavera       1   -1    1    -1
# I x P           1   -1   -1     1


# Crear matrix temporal
contrastmatrix <- cbind(c(1,    1,   -1,    -1), 
                        c(1,   -1,    1,    -1),
                        c(1,   -1,   -1,     1)); contrastmatrix
##      [,1] [,2] [,3]
## [1,]    1    1    1
## [2,]    1   -1   -1
## [3,]   -1    1   -1
## [4,]   -1   -1    1
# Asignar matrix temporal a TRAT
contrasts(fDOBLE$TRAT)<-contrastmatrix
fDOBLE$TRAT
## attr(,"contrasts")
##           [,1] [,2] [,3]
## Inv+Prim     1    1    1
## Invierno     1   -1   -1
## Primavera   -1    1   -1
## SinTratar   -1   -1    1
## Levels: Inv+Prim Invierno Primavera SinTratar
# Correr nuevamente el ANOVA con interacciones
CONTRASTAR <- aov(AFECT ~ FECHA + FECHA/REP + FECHA*TRAT, fDOBLE)


# Separar suma de cuadrados para TRAT
summary(CONTRASTAR, split = list(TRAT = list("Invierno" = 1,"Primavera" = 2, "I x P" = 3)))
##                         Df Sum Sq Mean Sq F value   Pr(>F)    
## FECHA                    1     47      47   0.882  0.35693    
## TRAT                     3  15826    5275  99.302 1.15e-13 ***
##   TRAT: Invierno         1  15097   15097 284.184 8.32e-15 ***
##   TRAT: Primavera        1    634     634  11.942  0.00206 ** 
##   TRAT: I x P            1     95      95   1.780  0.19468    
## FECHA:REP                8     56       7   0.131  0.99709    
## FECHA:TRAT               3    580     193   3.638  0.02707 *  
##   FECHA:TRAT: Invierno   1     58      58   1.098  0.30518    
##   FECHA:TRAT: Primavera  1    432     432   8.138  0.00878 ** 
##   FECHA:TRAT: I x P      1     89      89   1.677  0.20761    
## Residuals               24   1275      53                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Nota: 
# La Suma de Cuadrados (SC) y los Cuadrados Medios (CM) son similares a los presentados 
# en el ANOVA de la página 184; no obstante, los valores de F y p son diferentes.


Codigos 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 años (páginas 180 - 184)
************************************************************/

DATA foliarDOBLE; 
INPUT FECHA BLO $ TRAT $ AFECT;
LABEL
       FECHA      = 'Año de Experimento'
       BLO        = 'Repeticiones'
       TRAT       = 'Momento de Aplicación'
       AFECT      = 'Tejido foliar afectado (Arcsen)';
DATALINES;
1      I      SinTratar    61.9
1      II     SinTratar    55.9
1      III    SinTratar    64.5
1      IV     SinTratar    57.7
1      V      SinTratar    50.9
1      I      Invierno     20.3
1      II     Invierno     16.5
1      III    Invierno     15
1      IV     Invierno     8.5
1      V      Invierno     23.8
1      I      Primavera    54
1      II     Primavera    60.4
1      III    Primavera    49.7
1      IV     Primavera    70.4
1      V      Primavera    49
1      I      Inv+Prim     14.4
1      II     Inv+Prim     16.1
1      III    Inv+Prim     13.4
1      IV     Inv+Prim     10.9
1      V      Inv+Prim     22.8
2      I      SinTratar    70.9
2      II     SinTratar    61.2
2      III    SinTratar    56.8
2      IV     SinTratar    62.4
2      V      SinTratar    64.5
2      I      Invierno     15.2
2      II     Invierno     29.9
2      III    Invierno     18.4
2      IV     Invierno     22.1
2      V      Invierno     17.7
2      I      Primavera    46.3
2      II     Primavera    32
2      III    Primavera    40.1
2      IV     Primavera    56
2      V      Primavera    38.4
2      I      Inv+Prim     4.8
2      II     Inv+Prim     14.5
2      III    Inv+Prim     19.6
2      IV     Inv+Prim     9.8
2      V      Inv+Prim     12.2
;
ODS HTML;
* Generando resultados en Rich Text Format (RTF) en estilo STATISTICAL;
* Otras opciones: LISTING, JOURNAL, JOURNAL2, RTF, ANALYSYS;
ODS RTF FILE  = 'C:\Users\Administrator\Desktop\SASoutput.rtf' STYLE= STATISTICAL;

* Imprimir datos originales;
PROC PRINT DATA = foliarDOBLE;
RUN;

* Estadísticas descriptivas;
PROC MEANS MEAN STD;
       VAR AFECT;
       CLASS FECHA TRAT;
RUN;

* Separando datos de Año 1;
DATA foliarUNO;
       SET foliarDOBLE(obs = 20); RUN; *Incluye las primeras 20 líneas;
PROC PRINT DATA = foliarUNO;
RUN;

* ANOVA Diseño Bloques Completos al Azar: Año 1;
PROC GLM DATA = foliarUNO;
       CLASS TRAT BLO;
       MODEL AFECT = BLO TRAT;
RUN;

* Separando datos de Año 2;
DATA foliarDOS;
       SET foliarDOBLE(firstobs = 21); RUN; *Lectura de datos inicia en línea 21;
PROC PRINT DATA = foliarDOS;
RUN;

PROC GLM DATA = foliarDOS;
       CLASS TRAT BLO;
       MODEL AFECT = BLO TRAT;
RUN;

* Opción 1: ANOVA Combinado AFECT (página 183);
PROC GLM DATA = foliarDOBLE;
       CLASS TRAT BLO FECHA;
       MODEL AFECT = FECHA BLO FECHA*BLO TRAT FECHA*TRAT;
RUN;

* Opción 2: ANOVA Combinado AFECT (página 183) Mejor!;
PROC MIXED DATA = foliarDOBLE METHOD = TYPE3;
       CLASS TRAT BLO FECHA;
       MODEL AFECT = FECHA TRAT FECHA*TRAT / DDFM=SATTERTH;
       RANDOM FECHA*BLO;
RUN;

* Opción 1: ANOVA Combinado con la descomposición de los tratamientos (página 184);
* SC y CM son similares a los de la tabla -- valores de F y p son diferentes;
PROC GLM DATA = foliarDOBLE;
       CLASS TRAT BLO FECHA;
       MODEL AFECT = FECHA BLO FECHA*BLO TRAT FECHA*TRAT / SOLUTION;
       RANDOM FECHA*BLO / TEST;
       *Contrastes          Tratamientos   IxP   I    P    ST;
       CONTRAST 'Invierno'  TRAT            1    1   -1    -1;
       CONTRAST 'Primavera' TRAT            1   -1    1    -1;
       CONTRAST 'I x P'     TRAT            1   -1   -1     1;
RUN;

* Opción 2: ANOVA Combinado con la descomposición de los tratamientos (página 184) Mejor!;
PROC MIXED DATA = foliarDOBLE METHOD = TYPE1;
       CLASS TRAT BLO FECHA;
       MODEL AFECT = FECHA TRAT FECHA*TRAT / DDFM=SATTERTH;
       RANDOM FECHA*BLO;
       *Contrastes          Tratamientos   IxP   I    P    ST;
       CONTRAST 'Invierno'  TRAT            1    1   -1    -1;
       CONTRAST 'Primavera' TRAT            1   -1    1    -1;
       CONTRAST 'I x P'     TRAT            1   -1   -1     1;
RUN;
QUIT;
ODS RTF CLOSE;
ODS HTML CLOSE;

-----------------------------------------------------------


Archivos con códigos completos:

R        : capítulo 13d.R, summarySE.R



No hay comentarios:

Publicar un comentario