setwd("~/Google Drive/Agrosavia/Colaboraciones/Laura/BD Estadística SGR Santander/data")
datos4<-read.table("cabaparcel.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 1208 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 1231 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 17512.73 17644.74 -8736.365
## fit.ar1.diam 2 20 16108.19 16240.19 -8034.093
anova(fit.ar1.diam)
## Denom. DF: 5434
## numDF F-value p-value
## (Intercept) 1 21937.577 <.0001
## semana 1 3806.706 <.0001
## forestal 2 30.815 <.0001
## gen 4 6.988 <.0001
## bloque 2 124.657 <.0001
## forestal:gen 8 14.813 <.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 56974.25 57106.18 -28467.13
## fit.ar1.alt 2 20 55053.18 55185.11 -27506.59
anova(fit.ar1.alt)
## Denom. DF: 5411
## numDF F-value p-value
## (Intercept) 1 17353.454 <.0001
## semana 1 2177.124 <.0001
## forestal 2 27.766 <.0001
## gen 4 8.630 <.0001
## bloque 2 109.474 <.0001
## forestal:gen 8 12.688 <.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 ab CCN51
## TCS01 a TCS01
## TCS06 ab TCS06
## TCS13 c TCS13
## TCS19 b 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 ab Roble
## Terminalia a 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 defg Abarco:TCS01
## Abarco:TCS06 def Abarco:TCS06
## Abarco:TCS13 efgh Abarco:TCS13
## Abarco:TCS19 b Abarco:TCS19
## Roble:CCN51 fgh Roble:CCN51
## Roble:TCS01 cde Roble:TCS01
## Roble:TCS06 defg Roble:TCS06
## Roble:TCS13 h Roble:TCS13
## Roble:TCS19 gh Roble:TCS19
## Terminalia:CCN51 defg Terminalia:CCN51
## Terminalia:TCS01 defg Terminalia:TCS01
## Terminalia:TCS06 ac Terminalia:TCS06
## Terminalia:TCS13 defg Terminalia:TCS13
## Terminalia:TCS19 cd 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 b CCN51
## TCS01 b TCS01
## TCS06 a TCS06
## TCS13 b TCS13
## TCS19 c 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 a Abarco
## Roble a Roble
## Terminalia a 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 a Abarco:CCN51
## Abarco:TCS01 de Abarco:TCS01
## Abarco:TCS06 bcd Abarco:TCS06
## Abarco:TCS13 bcd Abarco:TCS13
## Abarco:TCS19 g Abarco:TCS19
## Roble:CCN51 f Roble:CCN51
## Roble:TCS01 cd Roble:TCS01
## Roble:TCS06 bcd Roble:TCS06
## Roble:TCS13 def Roble:TCS13
## Roble:TCS19 de Roble:TCS19
## Terminalia:CCN51 ef Terminalia:CCN51
## Terminalia:TCS01 cd Terminalia:TCS01
## Terminalia:TCS06 abc Terminalia:TCS06
## Terminalia:TCS13 bcd Terminalia:TCS13
## Terminalia:TCS19 ab 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.83 0.0367 5434 3.76 3.90
## TCS01 3.98 0.0358 5434 3.91 4.05
## TCS06 3.85 0.0356 5434 3.78 3.92
## TCS13 4.16 0.0371 5434 4.09 4.23
## TCS19 3.70 0.0385 5434 3.63 3.78
##
## 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.1498 0.0506 5434 -2.960 0.0257
## CCN51 - TCS06 -0.0212 0.0505 5434 -0.419 0.9935
## CCN51 - TCS13 -0.3303 0.0514 5434 -6.427 <.0001
## CCN51 - TCS19 0.1296 0.0523 5434 2.480 0.0954
## TCS01 - TCS06 0.1286 0.0500 5434 2.572 0.0757
## TCS01 - TCS13 -0.1805 0.0509 5434 -3.543 0.0037
## TCS01 - TCS19 0.2794 0.0519 5434 5.388 <.0001
## TCS06 - TCS13 -0.3091 0.0508 5434 -6.082 <.0001
## TCS06 - TCS19 0.1508 0.0518 5434 2.914 0.0295
## TCS13 - TCS19 0.4599 0.0526 5434 8.746 <.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.72 0.0292 5434 3.66 3.77
## Roble 4.11 0.0290 5434 4.05 4.17
## Terminalia 3.89 0.0282 5434 3.84 3.95
##
## 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.396 0.0400 5434 -9.893 <.0001
## Abarco - Terminalia -0.175 0.0397 5434 -4.409 <.0001
## Roble - Terminalia 0.221 0.0395 5434 5.591 <.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.36 0.0652 5434 3.23 3.48
## TCS01 4.04 0.0627 5434 3.91 4.16
## TCS06 3.95 0.0612 5434 3.83 4.07
## TCS13 4.09 0.0635 5434 3.96 4.21
## TCS19 3.15 0.0670 5434 3.02 3.28
##
## forestal = Roble:
## gen emmean SE df lower.CL upper.CL
## CCN51 4.16 0.0628 5434 4.03 4.28
## TCS01 3.88 0.0618 5434 3.76 4.00
## TCS06 3.97 0.0612 5434 3.85 4.09
## TCS13 4.36 0.0637 5434 4.24 4.49
## TCS19 4.19 0.0666 5434 4.06 4.32
##
## forestal = Terminalia:
## gen emmean SE df lower.CL upper.CL
## CCN51 3.98 0.0607 5434 3.86 4.10
## TCS01 4.03 0.0603 5434 3.91 4.15
## TCS06 3.64 0.0616 5434 3.52 3.76
## TCS13 4.04 0.0638 5434 3.91 4.16
## TCS19 3.77 0.0640 5434 3.64 3.89
##
## 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.67883 0.0900 5434 -7.546 <.0001
## CCN51 - TCS06 -0.59391 0.0890 5434 -6.675 <.0001
## CCN51 - TCS13 -0.73163 0.0905 5434 -8.084 <.0001
## CCN51 - TCS19 0.20840 0.0928 5434 2.245 0.1636
## TCS01 - TCS06 0.08492 0.0873 5434 0.973 0.8675
## TCS01 - TCS13 -0.05279 0.0889 5434 -0.594 0.9760
## TCS01 - TCS19 0.88723 0.0913 5434 9.717 <.0001
## TCS06 - TCS13 -0.13771 0.0879 5434 -1.568 0.5184
## TCS06 - TCS19 0.80231 0.0903 5434 8.882 <.0001
## TCS13 - TCS19 0.94002 0.0918 5434 10.236 <.0001
##
## forestal = Roble:
## contrast estimate SE df t.ratio p.value
## CCN51 - TCS01 0.27430 0.0876 5434 3.130 0.0151
## CCN51 - TCS06 0.19045 0.0873 5434 2.183 0.1863
## CCN51 - TCS13 -0.20621 0.0888 5434 -2.322 0.1380
## CCN51 - TCS19 -0.03647 0.0909 5434 -0.401 0.9945
## TCS01 - TCS06 -0.08385 0.0866 5434 -0.968 0.8696
## TCS01 - TCS13 -0.48051 0.0883 5434 -5.444 <.0001
## TCS01 - TCS19 -0.31077 0.0904 5434 -3.439 0.0053
## TCS06 - TCS13 -0.39667 0.0879 5434 -4.514 0.0001
## TCS06 - TCS19 -0.22692 0.0900 5434 -2.522 0.0859
## TCS13 - TCS19 0.16975 0.0915 5434 1.856 0.3415
##
## forestal = Terminalia:
## contrast estimate SE df t.ratio p.value
## CCN51 - TCS01 -0.04487 0.0853 5434 -0.526 0.9847
## CCN51 - TCS06 0.33991 0.0862 5434 3.944 0.0008
## CCN51 - TCS13 -0.05306 0.0877 5434 -0.605 0.9744
## CCN51 - TCS19 0.21695 0.0879 5434 2.469 0.0978
## TCS01 - TCS06 0.38478 0.0859 5434 4.479 0.0001
## TCS01 - TCS13 -0.00819 0.0875 5434 -0.094 1.0000
## TCS01 - TCS19 0.26182 0.0876 5434 2.989 0.0235
## TCS06 - TCS13 -0.39298 0.0883 5434 -4.450 0.0001
## TCS06 - TCS19 -0.12296 0.0884 5434 -1.390 0.6339
## TCS13 - TCS19 0.27002 0.0899 5434 3.002 0.0226
##
## 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 154 1.42 5411 151 157
## TCS01 152 1.39 5411 149 154
## TCS06 146 1.38 5411 143 148
## TCS13 151 1.44 5411 148 154
## TCS19 134 1.49 5411 131 137
##
## 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 2.600 1.96 5411 1.329 0.6729
## CCN51 - TCS06 8.632 1.95 5411 4.417 0.0001
## CCN51 - TCS13 3.073 1.99 5411 1.547 0.5320
## CCN51 - TCS19 20.502 2.02 5411 10.160 <.0001
## TCS01 - TCS06 6.032 1.94 5411 3.114 0.0159
## TCS01 - TCS13 0.472 1.97 5411 0.240 0.9993
## TCS01 - TCS19 17.901 2.00 5411 8.932 <.0001
## TCS06 - TCS13 -5.560 1.97 5411 -2.823 0.0384
## TCS06 - TCS19 11.869 2.00 5411 5.929 <.0001
## TCS13 - TCS19 17.429 2.03 5411 8.575 <.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 139 1.13 5411 137 142
## Roble 154 1.13 5411 152 157
## Terminalia 148 1.09 5411 146 150
##
## 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 -15.09 1.55 5411 -9.745 <.0001
## Abarco - Terminalia -8.69 1.53 5411 -5.663 <.0001
## Roble - Terminalia 6.40 1.53 5411 4.190 0.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.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 134 2.52 5411 129 139
## TCS01 154 2.42 5411 150 159
## TCS06 148 2.38 5411 143 153
## TCS13 148 2.45 5411 143 153
## TCS19 112 2.59 5411 107 117
##
## forestal = Roble:
## gen emmean SE df lower.CL upper.CL
## CCN51 166 2.42 5411 161 171
## TCS01 150 2.41 5411 145 154
## TCS06 148 2.38 5411 143 153
## TCS13 156 2.47 5411 151 161
## TCS19 152 2.57 5411 147 157
##
## forestal = Terminalia:
## gen emmean SE df lower.CL upper.CL
## CCN51 162 2.34 5411 158 167
## TCS01 151 2.33 5411 146 156
## TCS06 141 2.38 5411 136 145
## TCS13 149 2.47 5411 144 154
## TCS19 137 2.47 5411 132 142
##
## 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 -20.1180 3.47 5411 -5.794 <.0001
## CCN51 - TCS06 -13.7327 3.45 5411 -3.986 0.0007
## CCN51 - TCS13 -13.8171 3.49 5411 -3.955 0.0007
## CCN51 - TCS19 22.5290 3.58 5411 6.287 <.0001
## TCS01 - TCS06 6.3852 3.38 5411 1.889 0.3232
## TCS01 - TCS13 6.3009 3.43 5411 1.837 0.3524
## TCS01 - TCS19 42.6470 3.52 5411 12.101 <.0001
## TCS06 - TCS13 -0.0844 3.40 5411 -0.025 1.0000
## TCS06 - TCS19 36.2618 3.50 5411 10.370 <.0001
## TCS13 - TCS19 36.3461 3.54 5411 10.254 <.0001
##
## forestal = Roble:
## contrast estimate SE df t.ratio p.value
## CCN51 - TCS01 16.6358 3.40 5411 4.899 <.0001
## CCN51 - TCS06 18.0299 3.38 5411 5.334 <.0001
## CCN51 - TCS13 9.9190 3.44 5411 2.885 0.0321
## CCN51 - TCS19 13.7047 3.51 5411 3.907 0.0009
## TCS01 - TCS06 1.3942 3.37 5411 0.414 0.9939
## TCS01 - TCS13 -6.7167 3.43 5411 -1.958 0.2868
## TCS01 - TCS19 -2.9311 3.50 5411 -0.837 0.9190
## TCS06 - TCS13 -8.1109 3.41 5411 -2.375 0.1222
## TCS06 - TCS19 -4.3252 3.48 5411 -1.241 0.7272
## TCS13 - TCS19 3.7856 3.54 5411 1.070 0.8222
##
## forestal = Terminalia:
## contrast estimate SE df t.ratio p.value
## CCN51 - TCS01 11.2830 3.29 5411 3.428 0.0055
## CCN51 - TCS06 21.5996 3.33 5411 6.493 <.0001
## CCN51 - TCS13 13.1158 3.39 5411 3.868 0.0011
## CCN51 - TCS19 25.2713 3.39 5411 7.452 <.0001
## TCS01 - TCS06 10.3165 3.32 5411 3.111 0.0161
## TCS01 - TCS13 1.8327 3.38 5411 0.542 0.9829
## TCS01 - TCS19 13.9882 3.38 5411 4.137 0.0003
## TCS06 - TCS13 -8.4838 3.41 5411 -2.485 0.0941
## TCS06 - TCS19 3.6717 3.41 5411 1.076 0.8192
## TCS13 - TCS19 12.1555 3.48 5411 3.497 0.0043
##
## Results are averaged over the levels of: bloque
## P value adjustment: tukey method for comparing a family of 5 estimates
detach(datos4)