 |
|
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.
 |
|
 |
|
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 factoresCULTIVO.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 interactionsFUNGICIDA.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 ' ' 1plot(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 -1CONTRASTAR <- 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 ' ' 1Có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;
------------------------------------------------------------------------
Archivos:
capítulo 11.R, capítulo 11.sas, tabla 11.1.csv, schwarz.functions.r