febrero 28, 2017

Diseño en Parcelas Dividas


Instalamos tres librerías y una funcion para ejecutar este análisis con R: library(ggplot2), library(plyr), library (lsmeans), schwarz.functions.r.


# Código de gráfico de interacción
plot.interaction <- ggplot(data = REPORTE, aes(x = VARIEDAD, y = mean, 
                                            group    = FUNGICIDA, 
                                            color    = FUNGICIDA, 
                                            linetype = FUNGICIDA, 
                                            shape    = FUNGICIDA )) +
  ggtitle("Gráfico de interacción de tratamientos") +
  xlab("VARIEDAD") + ylab("kg/parcela") + geom_point(size = 3) + geom_line()
plot.interaction



El gráfico anterior y el siguiente tienen más sentido cuando es tiempo el segundo factor del experimento. Cuando las unidades experimentales son medidas en varios momentos.




# Gráfico de interacción: Cambiando valores en coordenadas 
plot.interaction <- ggplot(data = REPORTE, aes(x = FUNGICIDA, y = mean, 
                                               group    = VARIEDAD, 
                                               color    = VARIEDAD, 
                                               linetype = VARIEDAD, 
                                               shape    = VARIEDAD )) +
  ggtitle("Interacción: invirtiendo tratamientos en coordenadas") +
  xlab("FUNGICIDA") + ylab("kg/parcela") + geom_point(size = 3) + geom_line()
plot.interaction



El Diseño en Parcelas Dividas ofrece muchas posibilidades en experimentación (Capítulo 11, pp. 135–148). Preparar suelo, por ejemplo, requiere terreno con dimensiones amplias que faciliten maniobrar maquinaria agrícola (ver 11 Two-factor Split-plot designs).



# Código para guardar gráfico en escritorio "Desktop"
ggsave(plot = plot.dotplot, file = "C:/Users/Administrator/Desktop/Box-and-Whisker-Plot.png", height = 4, width = 6, units = "in")



Salvedad: Las tablas de ANOVAs, los gráficos, y las pruebas de homogeneidad unifactorial, resultaron similares. Los contrastes ortogonales son diferentes. Los datos analizados aparecen en la Tabla 11.1 de Experimentación en Agricultura.



Códigos de R


#***********************************************************
# EXPERIMENTACIÓN EN AGRICULTURA                              
# Fernández Escobar, R.; Trapero, A.; Domínguez, J. 2010       
# CAPÍTULO 11: DISEÑO EN PARCELAS DIVIDIDAS                     
#************************************************************

#### Opciones útiles
options(useFancyQuotes=FALSE) # presenta resumén de resultados correctamente
options(error = expression(NULL))  # permite  correr instrucciones for ficheros, incluso si existen errores detectados
options(width=200)

#### Datos
CULTIVO <- read.csv('C:/Users/Administrator/Desktop/tabla 11.1.csv', header=TRUE)

# Convirtiendo VARIEDAD y FUNGICIDA a factor
CULTIVO$VARIEDAD     <- as.factor(CULTIVO$VARIEDAD)
CULTIVO$FUNGICIDA    <- as.factor(CULTIVO$FUNGICIDA)
CULTIVO$TRATAMIENTOS <- with(CULTIVO, interaction(VARIEDAD,FUNGICIDA))
attach(CULTIVO); str(CULTIVO); CULTIVO # ver datos importados
##    VARIEDAD FUNGICIDA REPETICION PESO TRATAMIENTOS
## 1         1         1          1    5          1.1
## 2         1         1          2   18          1.1
## 3         1         1          3    7          1.1
## 4         1         1          4   10          1.1
## 5         1         2          1   22          1.2
## 6         1         2          2   13          1.2
## 7         1         2          3    6          1.2
## 8         1         2          4   16          1.2
## 9         1         3          1   29          1.3
## 10        1         3          2   35          1.3
## 11        1         3          3   40          1.3
## 12        1         3          4   24          1.3
## 13        2         1          1   33          2.1
## 14        2         1          2   43          2.1
## 15        2         1          3   28          2.1
## 16        2         1          4   30          2.1
## 17        2         2          1   38          2.2
## 18        2         2          2   47          2.2
## 19        2         2          3   51          2.2
## 20        2         2          4   57          2.2
## 21        2         3          1   62          2.3
## 22        2         3          2   52          2.3
## 23        2         3          3   54          2.3
## 24        2         3          4   44          2.3
## 25        3         1          1   47          3.1
## 26        3         1          2   43          3.1
## 27        3         1          3   51          3.1
## 28        3         1          4   60          3.1
## 29        3         2          1   52          3.2
## 30        3         2          2   56          3.2
## 31        3         2          3   49          3.2
## 32        3         2          4   38          3.2
## 33        3         3          1   54          3.3
## 34        3         3          2   34          3.3
## 35        3         3          3   53          3.3
## 36        3         3          4   48          3.3
####  Gráficos
# Box and Whisker Plot y gráfico de puntos atípicos
library(ggplot2)
plot.dotplot <- ggplot(data = CULTIVO, aes(x = TRATAMIENTOS, y = PESO))+
  ggtitle("Box-and-Whisker con puntos para inspeccionar valores atípicos") +
  xlab("Tratamientos (VARIEDAD:FUNGICIDA)") + ylab("kg/parcela") +
  geom_point(position = position_jitter(height = 0.2, width = 0.2), size = 2) +
  geom_boxplot(alpha = 0.1)
plot.dotplot
# Guardar gráfico en escritorio "Desktop"
ggsave(plot = plot.dotplot, file = "C:/Users/Administrator/Desktop/Box-and-Whisker-Plot.png", height = 4, width = 6, units = "in")
#### Estadísticas descriptivas para cada grupo
library(plyr)
source("C:/Users/Administrator/Desktop/schwarz.functions.r") # Función Schwarz
REPORTE <- ddply(CULTIVO, c("FUNGICIDA","VARIEDAD"), sf.simple.summary, variable="PESO")
REPORTE 
##   FUNGICIDA VARIEDAD n nmiss  mean       sd
## 1         1        1 4     0 10.00 5.715476
## 2         1        2 4     0 33.50 6.658328
## 3         1        3 4     0 50.25 7.274384
## 4         2        1 4     0 14.25 6.652067
## 5         2        2 4     0 48.25 7.973916
## 6         2        3 4     0 48.75 7.719024
## 7         3        1 4     0 32.00 6.976150
## 8         3        2 4     0 53.00 7.393691
## 9         3        3 4     0 47.25 9.215024
# Verificando si la desviación típica incrementa con la media
ggplot(data = REPORTE, aes(x = log(mean), y = log(sd))) + 
  ggtitle("¿Incrementa la desviación típica con la media?") + 
  geom_point()
# Código de gráfico de interacción
plot.interaction <- ggplot(data = REPORTE, aes(x = VARIEDAD, y = mean, 
                                            group    = FUNGICIDA, 
                                            color    = FUNGICIDA, 
                                            linetype = FUNGICIDA, 
                                            shape    = FUNGICIDA )) +
  ggtitle("Gráfico de interacción de tratamientos") +
  xlab("VARIEDAD") + ylab("kg/parcela") + geom_point(size = 3) + geom_line()
plot.interaction
# Guardar gráfico
ggsave(plot = plot.interaction, 
       file = "C:/Users/Administrator/Desktop/gráfico-interaccion.png", height = 4, width = 6, units = "in")
#### ANOVA para efectos fijos
str(CULTIVO) # Variables relevantes deben ser factores
CULTIVO.fit.lmerTest <- lmerTest::lmer(PESO ~ FUNGICIDA + VARIEDAD + FUNGICIDA:VARIEDAD + (1 | REPETICION), data = CULTIVO)
anova(CULTIVO.fit.lmerTest)  # aproximación Satterhwaite 
anova(CULTIVO.fit.lmerTest, ddf = "Kenward-Roger")
anova(CULTIVO.fit.lmerTest, ddf = "lme4")
invisible(summary(CULTIVO.fit.lmerTest)) #tabla resumén
plot(CULTIVO.fit.lmerTest)
#### LSMEANS 
#### Resultados de VARIEDAD ajustados con Tukey y graficados con Compact Letter Displays (CLD)
library(lsmeans)
VARIEDAD.lsmo   <- lsmeans::lsmeans(CULTIVO.fit.lmerTest, ~ VARIEDAD)
VARIEDAD.cld <- cld(VARIEDAD.lsmo, adjust = 'tukey')
VARIEDAD.cld
##  VARIEDAD   lsmean       SE df lower.CL upper.CL .group
##  1        18.75000 2.120411 27 13.35342 24.14658  1    
##  2        44.91667 2.120411 27 39.52009 50.31325   2   
##  3        48.75000 2.120411 27 43.35342 54.14658   2   
## 
## Results are averaged over the levels of: FUNGICIDA 
## Degrees-of-freedom method: satterthwaite 
## Confidence level used: 0.95 
## Conf-level adjustment: sidak method for 3 estimates 
## P value adjustment: tukey method for comparing a family of 3 estimates 
## significance level used: alpha = 0.05
# Gráfico de barras de error
VARIEDAD.cld.plot1 <- sf.cld.plot.bar(VARIEDAD.cld, "VARIEDAD") +
  xlab   ("VARIEDAD") +
  ylab   ("Medias para cultivo con 95% de IC") +
  ggtitle("Comparación de medias de VARIEDAD con CLD")
VARIEDAD.cld.plot1
# Gráfico de línea con 95% intervalos de confianza
VARIEDAD.cld.plot2 <- sf.cld.plot.line(VARIEDAD.cld, "VARIEDAD") +
  xlab   ("VARIEDAD") +
  ylab   ("Medias para cultivo con 95% de IC") +
  ggtitle("Comparación de medias de VARIEDAD con CLD")
VARIEDAD.cld.plot2
# Resultados
VARIEDAD.pairs <- pairs(VARIEDAD.lsmo, adjust='tukey')
summary(VARIEDAD.pairs, infer=TRUE)
#### Resultados de FUNGICIDA ajustados con prueba de Tukey y graficados con CLD
FUNGICIDA.lsmo <- lsmeans::lsmeans(CULTIVO.fit.lmerTest, ~ FUNGICIDA)
## NOTE: Results may be misleading due to involvement in interactions
FUNGICIDA.cld <- cld(FUNGICIDA.lsmo, adjust = 'tukey')
FUNGICIDA.cld
##  FUNGICIDA   lsmean       SE df lower.CL upper.CL .group
##  1         31.25000 2.120411 27 25.85342 36.64658  1    
##  2         37.08333 2.120411 27 31.68675 42.47991  12   
##  3         44.08333 2.120411 27 38.68675 49.47991   2   
## 
## Results are averaged over the levels of: VARIEDAD 
## Degrees-of-freedom method: satterthwaite 
## Confidence level used: 0.95 
## Conf-level adjustment: sidak method for 3 estimates 
## P value adjustment: tukey method for comparing a family of 3 estimates 
## significance level used: alpha = 0.05
# Gráfico de barras de error
FUNGICIDA.cld.plot1 <- sf.cld.plot.bar(FUNGICIDA.cld, "FUNGICIDA") +
  xlab("FUNGICIDA") +
  ylab("Medias para cultivo con 95% de IC") +
  ggtitle("Comparación de medias de FUNGICIDA con CLD")
FUNGICIDA.cld.plot1
# Gráfico de línea con 95% intervalos de confianza
FUNGICIDA.cld.plot2 <- sf.cld.plot.line(FUNGICIDA.cld, "FUNGICIDA") +
  xlab("FUNGICIDA") +
  ylab("Medias para cultivo con 95% de IC") +
  ggtitle("Comparación de medias de FUNGICIDA con CLD")
FUNGICIDA.cld.plot2
# Resultados
FUNGICIDA.pairs <- pairs(FUNGICIDA.lsmo, adjust = 'tukey')
summary(FUNGICIDA.pairs, infer = TRUE)
#### LSMEANS para la interacción VARIEDAD:FUNGICIDA
# Simultáneos intervalos de confianza (IC)
tm.lsmo <- lsmeans::lsmeans(CULTIVO.fit.lmerTest, ~ FUNGICIDA:VARIEDAD)
tm.cld <- cld(tm.lsmo, adjust='tukey')
tm.cld
##  FUNGICIDA VARIEDAD lsmean       SE df  lower.CL upper.CL .group
##  1         1         10.00 3.672659 27 -1.034401  21.0344  1    
##  2         1         14.25 3.672659 27  3.215599  25.2844  1    
##  3         1         32.00 3.672659 27 20.965599  43.0344   2   
##  1         2         33.50 3.672659 27 22.465599  44.5344   23  
##  3         3         47.25 3.672659 27 36.215599  58.2844   234 
##  2         2         48.25 3.672659 27 37.215599  59.2844   234 
##  2         3         48.75 3.672659 27 37.715599  59.7844   234 
##  1         3         50.25 3.672659 27 39.215599  61.2844    34 
##  3         2         53.00 3.672659 27 41.965599  64.0344     4 
## 
## Degrees-of-freedom method: satterthwaite 
## Confidence level used: 0.95 
## Conf-level adjustment: sidak method for 9 estimates 
## P value adjustment: tukey method for comparing a family of 9 estimates 
## significance level used: alpha = 0.05
# Gráfico de barras de error, agregando una nueva variable para crear el reporte
tm.cld$TM    <- interaction(tm.cld$FUNGICIDA, tm.cld$VARIEDAD)
tm.cld.plot1 <- sf.cld.plot.bar(tm.cld, "TM")  +
  xlab   ("VARIEDAD.FUNGICIDA")                +
  ylab   ("Medias para cultivo con 95% de IC") +
  ggtitle("Medias de VARIEDAD:FUNGICIDAD con comparaciones CLD")
tm.cld.plot1
# Gráfico de línea
tm.cld.plot2 <- sf.cld.plot.line(tm.cld, "TM") +
  xlab   ("FUNGICIDA.VARIEDAD")                +
  ylab   ("Medias para cultivo con 95% de IC") +
  ggtitle("Medias de VARIEDAD:FUNGICIDAD con comparaciones CLD")
tm.cld.plot2
# Resultados
tm.pairs <- pairs(tm.lsmo, adjust='tukey')
summary(tm.pairs, infer=TRUE)
#### ANOVA Unifactorial para VARIEDAD, Prueba de Bartlett, Shapiro-Wilk
# Diseño Completamente Aleatorio para Peso (VARIEDAD)
anova <- aov(PESO ~ VARIEDAD, data = CULTIVO); summary(anova)
model.tables(anova, type = "means", se = T)
bartlett.test(PESO ~ VARIEDAD, data = CULTIVO)
shapiro.test(anova$residuals)
# Diseño Completamente Aleatorio para Peso (FUNGICIDA)
anova <- aov(PESO ~ FUNGICIDA, data = CULTIVO); summary(anova)
model.tables(anova, type = "means", se = T)
bartlett.test(PESO ~ FUNGICIDA, data = CULTIVO)
shapiro.test(anova$residuals)
#### Gráfico de diagnósticos de ANOVA
par(mfrow=c(1,2))
DIAGNOSTICOS <- aov(PESO ~ VARIEDAD * FUNGICIDA + REPETICION:VARIEDAD, data = CULTIVO)
summary(DIAGNOSTICOS)
##                     Df Sum Sq Mean Sq F value   Pr(>F)    
## VARIEDAD             2   6398    3199  53.426 1.45e-09 ***
## FUNGICIDA            2    991     495   8.275  0.00185 ** 
## VARIEDAD:FUNGICIDA   4    944     236   3.941  0.01346 *  
## VARIEDAD:REPETICION  3     20       7   0.110  0.95334    
## Residuals           24   1437      60                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(DIAGNOSTICOS)
#### Contrastes ortogonales
# Asignar coeficientes de contraste
# Orden de tratamientos -----------V1F0 V2F0  V3F0  V1F1 V2F1 V3F1 V1F2 V2F2 V3F2;
# CONTRAST 'V1F0  vs V1F1+V1F2'      2    0    0    -1    0    0   -1    0    0; 
# CONTRAST 'V1F1  vs      V1F2'      0    0    0     1    0    0   -1    0    0;
# CONTRAST 'V2F0  vs V2F1+V2F2'      0    2    0     0   -1    0    0   -1    0;
# CONTRAST 'V2F1  vs      V2F2'      0    0    0     0    1    0    0   -1    0;
# CONTRAST 'V3F0  vs  V3F1+V3F2'     0    0    2     0    0   -1    0    0   -1;
# CONTRAST 'V3F1  vs      V3F2'      0    0    0     0    0    1    0    0   -1; 
CULTIVO$TRATAMIENTOS <- as.factor(CULTIVO$TRATAMIENTOS)
# Verificar con str() los niveles y orden de tratamientos
str(CULTIVO$TRATAMIENTOS)
# Creando Matrix
contrastmatrix <- cbind(c(2, 0, 0, -1, 0, 0, -1, 0, 0), c(0, 0, 0, 1, 0, 0, -1, 0, 0), c(0, 2, 0, 0, -1, 0, 0, -1, 0), c(0, 0, 0, 0, 1, 0, 0, -1, 0), c(0, 0, 2, 0, 0, -1, 0, 0, -1), c(0, 0, 0, 0, 0, 1, 0, 0, -1)); contrastmatrix
##       [,1] [,2] [,3] [,4] [,5] [,6]
##  [1,]    2    0    0    0    0    0
##  [2,]    0    0    2    0    0    0
##  [3,]    0    0    0    0    2    0
##  [4,]   -1    1    0    0    0    0
##  [5,]    0    0   -1    1    0    0
##  [6,]    0    0    0    0   -1    1
##  [7,]   -1   -1    0    0    0    0
##  [8,]    0    0   -1   -1    0    0
##  [9,]    0    0    0    0   -1   -1
CONTRASTAR <- aov(PESO ~ TRATAMIENTOS, CULTIVO)
summary(CONTRASTAR, split = list(TRATAMIENTOS = list("V1F0 vs V1F1+V1F2" = 1, "V1F1 vs V1F2" = 2, "V2F0 vs V2F1+V2F2" = 3, "V2F1 vs V2F2"  = 4, "V3F0 vs V3F1+V3F2" = 5, "V3F1 vs V3F2" = 6)))
##                                   Df Sum Sq Mean Sq F value   Pr(>F)    
## TRATAMIENTOS                       8   8332  1041.5  19.304 2.68e-09 ***
##   TRATAMIENTOS: V1F0 vs V1F1+V1F2  1     71    71.0   1.316  0.26137    
##   TRATAMIENTOS: V1F1 vs V1F2       1    690   689.5  12.780  0.00135 ** 
##   TRATAMIENTOS: V2F0 vs V2F1+V2F2  1   2251  2251.3  41.727 6.36e-07 ***
##   TRATAMIENTOS: V2F1 vs V2F2       1    337   336.7   6.240  0.01888 *  
##   TRATAMIENTOS: V3F0 vs V3F1+V3F2  1    557   556.5  10.315  0.00340 ** 
##   TRATAMIENTOS: V3F1 vs V3F2       1     68    67.7   1.255  0.27255    
## Residuals                         27   1457    54.0                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1



Códigos de SAS


/***********************************************************
 EXPERIMENTACIÓN EN AGRICULTURA                             
 Fernández Escobar, R.; Trapero, A.; Domínguez, J. 2010      
 CAPÍTULO 11: DISEÑO EN PARCELAS DIVIDIDAS                 
************************************************************/

DATA cultivo;
INPUT VAR FUN REP PESO;
CARDS;
1      1      1      5
1      1      2      18
1      1      3      7
1      1      4      10
1      2      1      22
1      2      2      13
1      2      3      6
1      2      4      16
1      3      1      29
1      3      2      35
1      3      3      40
1      3      4      24
2      1      1      33
2      1      2      43
2      1      3      28
2      1      4      30
2      2      1      38
2      2      2      47
2      2      3      51
2      2      4      57
2      3      1      62
2      3      2      52
2      3      3      54
2      3      4      44
3      1      1      47
3      1      2      43
3      1      3      51
3      1      4      60
3      2      1      52
3      2      2      56
3      2      3      49
3      2      4      38
3      3      1      54
3      3      2      34
3      3      3      53
3      3      4      48
;
ODS HTML;
PROC PRINT;
PROC MEANS MEAN STD STDERR;
       VAR VAR;
       CLASS FUN;
RUN;

* ANOVA parcelas divididas (primera opcion: valor-p no significativo en interaccion);
PROC GLM;
       CLASS REP VAR FUN;
       MODEL PESO = REP VAR REP*VAR FUN REP*FUN VAR*FUN / SS2; *Error type = 2;
       TEST H=REP VAR E=REP*VAR;
       LSMEANS VAR / STDERR E=REP*VAR;
       LSMEANS FUN / STDERR;
       LSMEANS VAR*FUN/ ETYPE=2;
RUN;

* ANOVA parcelas divididas (segunda opcion ¡Mejor!);
PROC MIXED METHOD = TYPE2;
       CLASS REP VAR FUN;
       MODEL PESO = REP VAR FUN VAR*FUN /DDFM =SATTERTH;
       RANDOM REP*VAR;
       * Contrastes;      
       CONTRAST 'V1 vs media(V2,V3)' VAR 2 -1 -1;
       CONTRAST 'V2 vs  V3         ' VAR 0  1 -1;
       CONTRAST 'FO vs media(F1,F2)' FUN 2 -1 -1;
       CONTRAST 'F1 vs  F2         ' FUN 0  1 -1;
       CONTRAST "VAR cuadratica * FUN lineal" VAR*FUN -1  0  1 -1  0  1  2  0 -2;
       CONTRAST 'F0 vs F1 en V1' FUN 1 -1  0
       FUN*VAR 1 -1 0 0 0 0 0 0 0;
       LSMEANS VAR FUN FUN*VAR / PDIFF;
RUN;

* Presentacion visual de interaccion de medias
* Medias para cada variedad*fungicida;
PROC MEANS DATA = cultivo NOPRINT;
       VAR PESO;
       BY VAR FUN;
       OUTPUT OUT = DOS MEAN = MEDIAS;
RUN; PROC PRINT; RUN;
* Grafico;
       SYMBOL1 INTERPOL = JOIN VALUE = DOT;
PROC GPLOT DATA = DOS;
       PLOT MEDIAS*VAR=FUN;
RUN;

* ANOVA Unifactorial y Prueba de Bartlett;
PROC ANOVA DATA = cultivo;
       CLASS VAR;
       MODEL PESO = VAR; *DCA para PESO (VARIEDAD);
       MEANS VAR / HOVTEST = BARTLETT;
RUN;

PROC ANOVA DATA = cultivo;
       CLASS FUN;
       MODEL PESO = FUN; *DCA para PESO (FUNGICIDA);
       MEANS FUN / HOVTEST = BARTLETT;
RUN; QUIT;

ODS HTML CLOSE;


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