setwd("~/Google Drive/Agrosavia/Colaboraciones/Laura/BD Estadística SGR Santander/data")
datos4<-read.table("sanjoabando.csv", header=T, sep=',')
datos4$gen<-as.factor(datos4$gen)
datos4$forestal<-as.factor(datos4$forestal)
datos4$bloque<-as.factor(datos4$bloque)
attach(datos4)
library(ggplot2)
library(emmeans)
#Gráfica diámetro
ggplot(datos4, aes(semana, diam, group = gen, colour = gen)) +
facet_grid(~forestal) +
geom_smooth(method="lm", se=F) +
theme_classic() +
xlab ("Semana") +
ylab ("Diámetro") +
labs(colour = "Genotipo") +
theme_linedraw() +
theme(
plot.title = element_text(hjust = 0.5, size = 16),
strip.text = element_text(size = 16),
axis.title.y = element_text(size = 16),
axis.title.x = element_text(size = 16),
axis.text.x = element_text(size = 14),
axis.text.y = element_text(size = 14)
)
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 794 rows containing non-finite values (stat_smooth).
# Gráfica altura
ggplot(datos4, aes(semana, alt, group = gen, colour = gen)) +
facet_grid(~forestal) +
geom_smooth(method="lm", se=F) +
theme_classic() +
xlab ("Semana") +
ylab ("Altura") +
labs(colour = "Genotipo") +
theme_linedraw() +
theme(
plot.title = element_text(hjust = 0.5, size = 16),
strip.text = element_text(size = 16),
axis.title.y = element_text(size = 16),
axis.title.x = element_text(size = 16),
axis.text.x = element_text(size = 14),
axis.text.y = element_text(size = 14)
)
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 793 rows containing non-finite values (stat_smooth).
# Anova general
aov.diam<-aov(diam~semana+forestal+gen+bloque+forestal*gen, na.action=na.exclude)
aov.alt<-aov(alt~semana+forestal+gen+bloque+forestal*gen, na.action=na.exclude)
#Análisis para diámetro
library(nlme)
fit.compsym.diam <- gls(diam ~ semana+forestal+gen+bloque+forestal*gen, data=datos4, corr=corCompSymm(, form= ~ 1 | gen),na.action=na.exclude)
fit.ar1.diam <- gls(diam ~ semana+forestal+gen+bloque+forestal*gen, data=datos4, corr=corAR1(, form= ~ 1 | gen), na.action=na.exclude)
anova(fit.compsym.diam, fit.ar1.diam) #compares the models
## Model df AIC BIC logLik
## fit.compsym.diam 1 20 16644.01 16777.79 -8302.005
## fit.ar1.diam 2 20 15766.18 15899.97 -7863.092
anova(fit.ar1.diam)
## Denom. DF: 5938
## numDF F-value p-value
## (Intercept) 1 39831.31 <.0001
## semana 1 3839.21 <.0001
## forestal 2 10.25 <.0001
## gen 4 4.03 0.0029
## bloque 2 119.03 <.0001
## forestal:gen 8 12.32 <.0001
# Análisis para altura
fit.compsym.alt <- gls(alt ~ semana+forestal+gen+bloque+forestal*gen, data=datos4, corr=corCompSymm(, form= ~ 1 | gen),na.action=na.exclude)
fit.ar1.alt <- gls(alt ~ semana+forestal+gen+bloque+forestal*gen, data=datos4, corr=corAR1(, form= ~ 1 | gen), na.action=na.exclude)
anova(fit.compsym.alt, fit.ar1.alt) #compares the models
## Model df AIC BIC logLik
## fit.compsym.alt 1 20 59850.83 59984.61 -29905.41
## fit.ar1.alt 2 20 59070.41 59204.20 -29515.21
anova(fit.ar1.alt)
## Denom. DF: 5939
## numDF F-value p-value
## (Intercept) 1 37408.85 <.0001
## semana 1 3901.63 <.0001
## forestal 2 4.68 0.0093
## gen 4 10.61 <.0001
## bloque 2 98.61 <.0001
## forestal:gen 8 8.60 <.0001
#Tukey diámetro
library(multcompView)
interac.tuk.diam<-TukeyHSD(aov.diam, "forestal:gen", ordered = TRUE)
## Warning in replications(paste("~", xx), data = mf): non-factors ignored: semana
fores.tuk.diam<-TukeyHSD(aov.diam, "forestal", ordered = TRUE)
## Warning in replications(paste("~", xx), data = mf): non-factors ignored: semana
gen.tuk.diam<-TukeyHSD(aov.diam, "gen", ordered = TRUE)
## Warning in replications(paste("~", xx), data = mf): non-factors ignored: semana
#Tukey altura
interac.tuk.alt<-TukeyHSD(aov.alt, "forestal:gen", ordered = TRUE)
## Warning in replications(paste("~", xx), data = mf): non-factors ignored: semana
fores.tuk.alt<-TukeyHSD(aov.alt, "forestal", ordered = TRUE)
## Warning in replications(paste("~", xx), data = mf): non-factors ignored: semana
gen.tuk.alt<-TukeyHSD(aov.alt, "gen", ordered = TRUE)
## Warning in replications(paste("~", xx), data = mf): non-factors ignored: semana
#Etiquetas Tukey diámetro
#Genotipos
generate_label_df_gen_diam <- function(gen.tuk.diam, variable){
Tukey.levels <- gen.tuk.diam[[variable]][,4]
Tukey.labels <- data.frame(multcompLetters(Tukey.levels)['Letters'])
Tukey.labels$treatment=rownames(Tukey.labels)
Tukey.labels=Tukey.labels[order(Tukey.labels$treatment) , ]
return(Tukey.labels)
}
labels.gen.diam <- generate_label_df_gen_diam(gen.tuk.diam, "gen")
labels.gen.diam
## Letters treatment
## CCN51 a CCN51
## TCS01 a TCS01
## TCS06 a TCS06
## TCS13 b TCS13
## TCS19 a TCS19
# Forestal
generate_label_df_forestal_diam <- function(fores.tuk.diam, variable){
Tukey.levels <- fores.tuk.diam[[variable]][,2]
Tukey.labels <- data.frame(multcompLetters(Tukey.levels)['Letters'])
Tukey.labels$treatment=rownames(Tukey.labels)
Tukey.labels=Tukey.labels[order(Tukey.labels$treatment) , ]
return(Tukey.labels)
}
labels.forestal.diam <- generate_label_df_forestal_diam(fores.tuk.diam, "forestal")
labels.forestal.diam
## Letters treatment
## Abarco b Abarco
## Roble a Roble
## Terminalia ab Terminalia
# Interacción Forestal:Genotipo
generate_label_df_interac_diam <- function(interac.tuk.diam, variable){
Tukey.levels <- interac.tuk.diam[[variable]][,4]
Tukey.labels <- data.frame(multcompLetters(Tukey.levels)['Letters'])
Tukey.labels$treatment=rownames(Tukey.labels)
Tukey.labels=Tukey.labels[order(Tukey.labels$treatment) , ]
return(Tukey.labels)
}
labels.interac.diam <- generate_label_df_interac_diam(interac.tuk.diam, "forestal:gen")
labels.interac.diam
## Letters treatment
## Abarco:CCN51 ab Abarco:CCN51
## Abarco:TCS01 d Abarco:TCS01
## Abarco:TCS06 ab Abarco:TCS06
## Abarco:TCS13 d Abarco:TCS13
## Abarco:TCS19 ab Abarco:TCS19
## Roble:CCN51 ac Roble:CCN51
## Roble:TCS01 b Roble:TCS01
## Roble:TCS06 abc Roble:TCS06
## Roble:TCS13 cd Roble:TCS13
## Roble:TCS19 cd Roble:TCS19
## Terminalia:CCN51 ab Terminalia:CCN51
## Terminalia:TCS01 ab Terminalia:TCS01
## Terminalia:TCS06 abc Terminalia:TCS06
## Terminalia:TCS13 ab Terminalia:TCS13
## Terminalia:TCS19 ab Terminalia:TCS19
#Etiquetas Tukey altura
#Genotipos
generate_label_df_gen_alt <- function(gen.tuk.alt, variable){
Tukey.levels <- gen.tuk.alt[[variable]][,4]
Tukey.labels <- data.frame(multcompLetters(Tukey.levels)['Letters'])
Tukey.labels$treatment=rownames(Tukey.labels)
Tukey.labels=Tukey.labels[order(Tukey.labels$treatment) , ]
return(Tukey.labels)
}
labels.gen.alt <- generate_label_df_gen_alt(gen.tuk.alt, "gen")
labels.gen.alt
## Letters treatment
## CCN51 a CCN51
## TCS01 a TCS01
## TCS06 a TCS06
## TCS13 a TCS13
## TCS19 b TCS19
# Forestal
generate_label_df_forestal_alt <- function(fores.tuk.alt, variable){
Tukey.levels <- fores.tuk.alt[[variable]][,2]
Tukey.labels <- data.frame(multcompLetters(Tukey.levels)['Letters'])
Tukey.labels$treatment=rownames(Tukey.labels)
Tukey.labels=Tukey.labels[order(Tukey.labels$treatment) , ]
return(Tukey.labels)
}
labels.forestal.alt <- generate_label_df_forestal_alt(fores.tuk.alt, "forestal")
labels.forestal.alt
## Letters treatment
## Abarco b Abarco
## Roble a Roble
## Terminalia ab Terminalia
# Interacción Forestal:Genotipo
generate_label_df_interac_alt <- function(interac.tuk.alt, variable){
Tukey.levels <- interac.tuk.alt[[variable]][,4]
Tukey.labels <- data.frame(multcompLetters(Tukey.levels)['Letters'])
Tukey.labels$treatment=rownames(Tukey.labels)
Tukey.labels=Tukey.labels[order(Tukey.labels$treatment) , ]
return(Tukey.labels)
}
labels.interac.alt <- generate_label_df_interac_alt(interac.tuk.alt, "forestal:gen")
labels.interac.alt
## Letters treatment
## Abarco:CCN51 def Abarco:CCN51
## Abarco:TCS01 f Abarco:TCS01
## Abarco:TCS06 cd Abarco:TCS06
## Abarco:TCS13 ef Abarco:TCS13
## Abarco:TCS19 ab Abarco:TCS19
## Roble:CCN51 def Roble:CCN51
## Roble:TCS01 abc Roble:TCS01
## Roble:TCS06 cd Roble:TCS06
## Roble:TCS13 def Roble:TCS13
## Roble:TCS19 cd Roble:TCS19
## Terminalia:CCN51 def Terminalia:CCN51
## Terminalia:TCS01 acd Terminalia:TCS01
## Terminalia:TCS06 cde Terminalia:TCS06
## Terminalia:TCS13 acd Terminalia:TCS13
## Terminalia:TCS19 b Terminalia:TCS19
## Gráficas contrastes de medias diametro
#Gen
contrast <- emmeans(aov.diam, ~gen)
## NOTE: Results may be misleading due to involvement in interactions
plot(contrast, comparisons = TRUE, xlab ="Diámetro")
medias.gen <- emmeans(aov.diam, pairwise ~ gen)
## NOTE: Results may be misleading due to involvement in interactions
medias.gen
## $emmeans
## gen emmean SE df lower.CL upper.CL
## CCN51 3.64 0.0281 5938 3.59 3.70
## TCS01 3.72 0.0282 5938 3.67 3.78
## TCS06 3.69 0.0281 5938 3.63 3.74
## TCS13 3.85 0.0281 5938 3.80 3.91
## TCS19 3.66 0.0281 5938 3.60 3.71
##
## Results are averaged over the levels of: forestal, bloque
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## CCN51 - TCS01 -0.0809 0.0398 5938 -2.032 0.2505
## CCN51 - TCS06 -0.0443 0.0397 5938 -1.115 0.7986
## CCN51 - TCS13 -0.2117 0.0397 5938 -5.332 <.0001
## CCN51 - TCS19 -0.0135 0.0397 5938 -0.339 0.9972
## TCS01 - TCS06 0.0366 0.0398 5938 0.920 0.8892
## TCS01 - TCS13 -0.1308 0.0398 5938 -3.285 0.0091
## TCS01 - TCS19 0.0674 0.0398 5938 1.693 0.4382
## TCS06 - TCS13 -0.1675 0.0397 5938 -4.216 0.0002
## TCS06 - TCS19 0.0308 0.0397 5938 0.775 0.9378
## TCS13 - TCS19 0.1983 0.0397 5938 4.989 <.0001
##
## Results are averaged over the levels of: forestal, bloque
## P value adjustment: tukey method for comparing a family of 5 estimates
#Forestal
contrast <- emmeans(aov.diam, ~forestal)
## NOTE: Results may be misleading due to involvement in interactions
plot(contrast, comparisons = TRUE, xlab ="Diámetro")
medias.forestal <- emmeans(aov.diam, pairwise ~ forestal)
## NOTE: Results may be misleading due to involvement in interactions
medias.forestal
## $emmeans
## forestal emmean SE df lower.CL upper.CL
## Abarco 3.77 0.0214 5938 3.73 3.82
## Roble 3.75 0.0220 5938 3.71 3.79
## Terminalia 3.61 0.0220 5938 3.57 3.66
##
## Results are averaged over the levels of: gen, bloque
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## Abarco - Roble 0.0228 0.0306 5938 0.744 0.7371
## Abarco - Terminalia 0.1593 0.0307 5938 5.196 <.0001
## Roble - Terminalia 0.1365 0.0311 5938 4.390 <.0001
##
## Results are averaged over the levels of: gen, bloque
## P value adjustment: tukey method for comparing a family of 3 estimates
#Forestal*Gen
contrast <- emmeans(aov.diam, ~gen|forestal)
plot(contrast, comparisons = TRUE, xlab ="Diámetro")
medias.forestal.gen <- emmeans(aov.diam, pairwise ~ gen|forestal)
medias.forestal.gen
## $emmeans
## forestal = Abarco:
## gen emmean SE df lower.CL upper.CL
## CCN51 3.58 0.0480 5938 3.49 3.68
## TCS01 4.07 0.0469 5938 3.98 4.16
## TCS06 3.66 0.0480 5938 3.57 3.76
## TCS13 4.00 0.0481 5938 3.90 4.09
## TCS19 3.55 0.0478 5938 3.46 3.65
##
## forestal = Roble:
## gen emmean SE df lower.CL upper.CL
## CCN51 3.74 0.0485 5938 3.64 3.83
## TCS01 3.50 0.0508 5938 3.40 3.60
## TCS06 3.71 0.0488 5938 3.61 3.80
## TCS13 3.91 0.0488 5938 3.81 4.00
## TCS19 3.90 0.0487 5938 3.81 4.00
##
## forestal = Terminalia:
## gen emmean SE df lower.CL upper.CL
## CCN51 3.60 0.0494 5938 3.51 3.70
## TCS01 3.60 0.0489 5938 3.50 3.70
## TCS06 3.69 0.0491 5938 3.59 3.79
## TCS13 3.66 0.0490 5938 3.56 3.76
## TCS19 3.51 0.0496 5938 3.42 3.61
##
## Results are averaged over the levels of: bloque
## Confidence level used: 0.95
##
## $contrasts
## forestal = Abarco:
## contrast estimate SE df t.ratio p.value
## CCN51 - TCS01 -0.48624 0.0671 5938 -7.248 <.0001
## CCN51 - TCS06 -0.07889 0.0678 5938 -1.163 0.7726
## CCN51 - TCS13 -0.41217 0.0679 5938 -6.067 <.0001
## CCN51 - TCS19 0.03256 0.0677 5938 0.481 0.9891
## TCS01 - TCS06 0.40735 0.0671 5938 6.069 <.0001
## TCS01 - TCS13 0.07407 0.0672 5938 1.102 0.8055
## TCS01 - TCS19 0.51880 0.0670 5938 7.747 <.0001
## TCS06 - TCS13 -0.33328 0.0680 5938 -4.903 <.0001
## TCS06 - TCS19 0.11146 0.0677 5938 1.646 0.4681
## TCS13 - TCS19 0.44473 0.0678 5938 6.557 <.0001
##
## forestal = Roble:
## contrast estimate SE df t.ratio p.value
## CCN51 - TCS01 0.23952 0.0702 5938 3.411 0.0059
## CCN51 - TCS06 0.03287 0.0688 5938 0.478 0.9894
## CCN51 - TCS13 -0.16607 0.0688 5938 -2.413 0.1119
## CCN51 - TCS19 -0.16215 0.0687 5938 -2.361 0.1264
## TCS01 - TCS06 -0.20665 0.0704 5938 -2.934 0.0277
## TCS01 - TCS13 -0.40559 0.0705 5938 -5.757 <.0001
## TCS01 - TCS19 -0.40167 0.0704 5938 -5.709 <.0001
## TCS06 - TCS13 -0.19894 0.0690 5938 -2.882 0.0324
## TCS06 - TCS19 -0.19502 0.0689 5938 -2.830 0.0376
## TCS13 - TCS19 0.00392 0.0690 5938 0.057 1.0000
##
## forestal = Terminalia:
## contrast estimate SE df t.ratio p.value
## CCN51 - TCS01 0.00400 0.0695 5938 0.058 1.0000
## CCN51 - TCS06 -0.08678 0.0696 5938 -1.247 0.7239
## CCN51 - TCS13 -0.05698 0.0696 5938 -0.819 0.9249
## CCN51 - TCS19 0.08918 0.0699 5938 1.275 0.7065
## TCS01 - TCS06 -0.09078 0.0693 5938 -1.310 0.6851
## TCS01 - TCS13 -0.06098 0.0693 5938 -0.880 0.9042
## TCS01 - TCS19 0.08518 0.0696 5938 1.224 0.7376
## TCS06 - TCS13 0.02980 0.0694 5938 0.429 0.9929
## TCS06 - TCS19 0.17596 0.0698 5938 2.523 0.0858
## TCS13 - TCS19 0.14616 0.0697 5938 2.097 0.2215
##
## Results are averaged over the levels of: bloque
## P value adjustment: tukey method for comparing a family of 5 estimates
## Gráficas contrastes de medias altura
#Gen
contrast <- emmeans(aov.alt, ~gen)
## NOTE: Results may be misleading due to involvement in interactions
plot(contrast, comparisons = TRUE, xlab ="Altura")
medias.gen <- emmeans(aov.alt, pairwise ~ gen)
## NOTE: Results may be misleading due to involvement in interactions
medias.gen
## $emmeans
## gen emmean SE df lower.CL upper.CL
## CCN51 137 1.07 5939 135 139
## TCS01 134 1.07 5939 132 136
## TCS06 133 1.07 5939 131 136
## TCS13 137 1.07 5939 135 139
## TCS19 125 1.07 5939 123 127
##
## Results are averaged over the levels of: forestal, bloque
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## CCN51 - TCS01 3.1765 1.51 5939 2.100 0.2199
## CCN51 - TCS06 3.7671 1.51 5939 2.498 0.0911
## CCN51 - TCS13 -0.0294 1.51 5939 -0.020 1.0000
## CCN51 - TCS19 12.5018 1.51 5939 8.288 <.0001
## TCS01 - TCS06 0.5906 1.51 5939 0.390 0.9951
## TCS01 - TCS13 -3.2059 1.51 5939 -2.119 0.2120
## TCS01 - TCS19 9.3253 1.51 5939 6.163 <.0001
## TCS06 - TCS13 -3.7965 1.51 5939 -2.516 0.0871
## TCS06 - TCS19 8.7347 1.51 5939 5.788 <.0001
## TCS13 - TCS19 12.5312 1.51 5939 8.303 <.0001
##
## Results are averaged over the levels of: forestal, bloque
## P value adjustment: tukey method for comparing a family of 5 estimates
#Forestal
contrast <- emmeans(aov.alt, ~forestal)
## NOTE: Results may be misleading due to involvement in interactions
plot(contrast, comparisons = TRUE, xlab ="Altura")
medias.forestal <- emmeans(aov.alt, pairwise ~ forestal)
## NOTE: Results may be misleading due to involvement in interactions
medias.forestal
## $emmeans
## forestal emmean SE df lower.CL upper.CL
## Abarco 136 0.812 5939 135 138
## Roble 134 0.835 5939 132 135
## Terminalia 130 0.837 5939 129 132
##
## Results are averaged over the levels of: gen, bloque
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## Abarco - Roble 2.70 1.16 5939 2.321 0.0530
## Abarco - Terminalia 6.04 1.16 5939 5.182 <.0001
## Roble - Terminalia 3.34 1.18 5939 2.823 0.0132
##
## Results are averaged over the levels of: gen, bloque
## P value adjustment: tukey method for comparing a family of 3 estimates
#Forestal*Gen
contrast <- emmeans(aov.alt, ~gen|forestal)
plot(contrast, comparisons = TRUE, xlab ="Altura")
medias.forestal.gen <- emmeans(aov.alt, pairwise ~ gen|forestal)
medias.forestal.gen
## $emmeans
## forestal = Abarco:
## gen emmean SE df lower.CL upper.CL
## CCN51 138 1.82 5939 134 141
## TCS01 145 1.78 5939 141 148
## TCS06 133 1.82 5939 129 137
## TCS13 143 1.83 5939 139 146
## TCS19 123 1.81 5939 119 126
##
## forestal = Roble:
## gen emmean SE df lower.CL upper.CL
## CCN51 136 1.84 5939 133 140
## TCS01 127 1.93 5939 123 130
## TCS06 134 1.85 5939 130 137
## TCS13 138 1.86 5939 135 142
## TCS19 133 1.85 5939 129 137
##
## forestal = Terminalia:
## gen emmean SE df lower.CL upper.CL
## CCN51 138 1.88 5939 134 141
## TCS01 131 1.86 5939 127 134
## TCS06 134 1.87 5939 130 137
## TCS13 131 1.86 5939 127 134
## TCS19 118 1.88 5939 114 122
##
## Results are averaged over the levels of: bloque
## Confidence level used: 0.95
##
## $contrasts
## forestal = Abarco:
## contrast estimate SE df t.ratio p.value
## CCN51 - TCS01 -7.289 2.55 5939 -2.860 0.0345
## CCN51 - TCS06 4.547 2.58 5939 1.764 0.3947
## CCN51 - TCS13 -5.091 2.58 5939 -1.973 0.2796
## CCN51 - TCS19 14.709 2.57 5939 5.720 <.0001
## TCS01 - TCS06 11.836 2.55 5939 4.642 <.0001
## TCS01 - TCS13 2.198 2.55 5939 0.861 0.9111
## TCS01 - TCS19 21.998 2.54 5939 8.647 <.0001
## TCS06 - TCS13 -9.638 2.58 5939 -3.733 0.0018
## TCS06 - TCS19 10.162 2.57 5939 3.950 0.0008
## TCS13 - TCS19 19.800 2.58 5939 7.685 <.0001
##
## forestal = Roble:
## contrast estimate SE df t.ratio p.value
## CCN51 - TCS01 9.802 2.67 5939 3.675 0.0022
## CCN51 - TCS06 2.792 2.61 5939 1.069 0.8227
## CCN51 - TCS13 -1.876 2.61 5939 -0.718 0.9525
## CCN51 - TCS19 3.219 2.61 5939 1.234 0.7317
## TCS01 - TCS06 -7.009 2.68 5939 -2.620 0.0669
## TCS01 - TCS13 -11.678 2.68 5939 -4.363 0.0001
## TCS01 - TCS19 -6.582 2.67 5939 -2.463 0.0993
## TCS06 - TCS13 -4.668 2.62 5939 -1.780 0.3854
## TCS06 - TCS19 0.427 2.62 5939 0.163 0.9998
## TCS13 - TCS19 5.095 2.62 5939 1.945 0.2936
##
## forestal = Terminalia:
## contrast estimate SE df t.ratio p.value
## CCN51 - TCS01 7.017 2.64 5939 2.658 0.0605
## CCN51 - TCS06 3.962 2.64 5939 1.498 0.5638
## CCN51 - TCS13 6.879 2.64 5939 2.604 0.0697
## CCN51 - TCS19 19.577 2.66 5939 7.369 <.0001
## TCS01 - TCS06 -3.055 2.63 5939 -1.160 0.7740
## TCS01 - TCS13 -0.138 2.63 5939 -0.053 1.0000
## TCS01 - TCS19 12.560 2.64 5939 4.749 <.0001
## TCS06 - TCS13 2.917 2.63 5939 1.107 0.8028
## TCS06 - TCS19 15.615 2.65 5939 5.893 <.0001
## TCS13 - TCS19 12.699 2.65 5939 4.798 <.0001
##
## Results are averaged over the levels of: bloque
## P value adjustment: tukey method for comparing a family of 5 estimates
detach(datos4)