Modelos Lineales: Diseños experimentales
E. Martinez y E. Maure
2022-10-17
Introducción
Anteriormente hemos trabajado con modelos lineales donde interesa la relación funcional entre variables cuantitativas. En reiteradas ocasiones el objetivo que se sigue no es ver la relación sino cómo una variable o varias a¿fectan la variabilidad de una característica en estudio. Generalmente los datos para este tipo de análisis son experimentales, donde se controlan las variables externas y las condiciones del experimento, pero hay veces que esos datos son observacionales y el análisis es similar. El análisis de estos datos se realiza con el ADEVA (ANOVA en inglés), a fin de realizar la comparación de las respuestas medias correspondiente a los tratamientos aplicados, ya sean simples o compuestos. En ocasiones, interesa examinar el efecto que se produce en la respuesta al aplicar simultáneamente más de un factor sobre una misma unidad experimental. Cuando se investiga el efecto de dos o más factores, no es apropiado analizar la variación que provoca en una respuesta cada uno de ellos por separado, dado que puede suceder que los efectos medidos en la variable respuesta se deban a la acción de los factores que actúan como variables independientes entre sí o que interactúen de modo no independiente. En consecuencia, para medir su eficacia hay que experimentar con ellos en forma simultánea, aplicándolos sobre una misma unidad experimental. Hay que llevar a cabo lo que se conoce como un experimento polifactorial o multifactorial o simplemente factorial. Genéricamente se trata de experimentos en los cuales se investiga la acción de dos o más factores, y donde los tratamientos responden a arreglos combinados de tratamientos simples (tratamientos combinados).
Ejemplo con un factor
N. glauca es endémica de la cordillera austral y está sujeta a intensa presión antrópica. Se estudió el efecto del ácido indolbutílico (AIB) su capacidad de arraigamiento.
Se colectaron 18 estacas de tallo de árboles de más de 20 años de edad que se pusieron a enraizar en macetas individuales con un mismo sustrato. Cada maceta reicibó de menera aleatoria e independiente una de las 3 dosis de AIB estudiadas: 0, 0.5, 1,5%. Así cada dosis fue repetida 6 veces.
Al cabo de tres meses se midió la longitud de la raíz (en cm).
Para el análisis y presentación de datos utilizaremos los paquetes “agricolae” y “pacman”
library(agricolae)
library(pacman)
library(tidyverse)## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.3.6 ✔ purrr 0.3.4
## ✔ tibble 3.1.8 ✔ dplyr 1.0.10
## ✔ tidyr 1.2.1 ✔ stringr 1.4.1
## ✔ readr 2.1.2 ✔ forcats 0.5.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()aib = read.csv2("AIB.csv", row.names=NULL)El paquete ‘pacman’ debe instalarse con `install.packages(“pacman”) Sus componentes son: pacman::p_load( tidyverse, # manejo general de datos y graficos rio, # importacion de datos agricolae, # algunas funciones para diseño/analisis de experimento nlme, # modelos mixtos emmeans, # estimacion y comparacion de medias DescTools, # funciones para comparacione test car # funciones para comparacione test )
Verificar y/o convertir tipos de datos
head(aib)## Dosis... Long.mm. X X.1
## 1 0 0.40 NA NA
## 2 0 1.32 NA NA
## 3 0 1.02 NA NA
## 4 0 2.47 NA NA
## 5 0 2.12 NA NA
## 6 0 1.66 NA NAstr(aib)## 'data.frame': 18 obs. of 4 variables:
## $ Dosis...: num 0 0 0 0 0 0 0.5 0.5 0.5 0.5 ...
## $ Long.mm.: num 0.4 1.32 1.02 2.47 2.12 1.66 5.4 5.59 4.82 3.67 ...
## $ X : logi NA NA NA NA NA NA ...
## $ X.1 : logi NA NA NA NA NA NA ...Para eliminar las columnas que están vacías podemos utilizar la selección de columnas y nombrar el objeto de la misma manera:
aib = aib[,1:2]
head(aib)## Dosis... Long.mm.
## 1 0 0.40
## 2 0 1.32
## 3 0 1.02
## 4 0 2.47
## 5 0 2.12
## 6 0 1.66Podemos agregar una columna para definir a dosis como factor.
aib$dosis = aib[,1]
aib$dosis = as.factor(aib$dosis)
str(aib)## 'data.frame': 18 obs. of 3 variables:
## $ Dosis...: num 0 0 0 0 0 0 0.5 0.5 0.5 0.5 ...
## $ Long.mm.: num 0.4 1.32 1.02 2.47 2.12 1.66 5.4 5.59 4.82 3.67 ...
## $ dosis : Factor w/ 3 levels "0","0.5","1.5": 1 1 1 1 1 1 2 2 2 2 ...aib = aib[,2:3]Explorar
Descriptiva
En la descripción numérica, por factor, armamos un resumen de las medidas que nos interesan:
aib %>%
group_by(dosis) %>%
select(- dosis) %>%
summarise_all(.funs = c(
n = length,
media = mean,
DE = sd,
min = min,
max = max
)
)## Adding missing grouping variables: `dosis`## # A tibble: 3 × 6
## dosis n media DE min max
## <fct> <int> <dbl> <dbl> <dbl> <dbl>
## 1 0 6 1.50 0.751 0.4 2.47
## 2 0.5 6 4.92 0.743 3.67 5.59
## 3 1.5 6 5.83 0.554 5.27 6.46Gráficos de puntos
Para la parte de representación gráfica debemos tener en cuenta que las características pueden ir cambiandose de acuerdo a como cada uno lo considere necesario, en cuanto a detalles de forma, color, etc. Para mayor detalle pueden consultarse por Google, por ejemplo.
ggplot(aib) +
aes(x = dosis, y = Long.mm.) +
geom_jitter(width = 0.2) +
stat_summary(fun.data = "mean_se", geom = "pointrange", color = "blue")
Gráficos boxplot
Hemos visto también que los diagramas de caja nos dan mayor información respecto a las medidas de posición y distribución de los datos.
ggplot(aib) +
aes(x = dosis, y = Long.mm.) +
geom_boxplot()
Modelo lineal
El modelo que utilizaremos es de efectos fijos ya que hemos fijado los niveles del factor, en este caso es uno: “Dosis de AIB”, por lo tanto el modelo a comprobar es: \[Y_{ij} =\mu + \alpha_j +\varepsilon_{ij}\] Teniendo en cuenta que: \(Y_{ij}\) : es la respuesta de la i-ésima estaca que se le aplicó la j-ésima dosis de AIB.
\(\mu\): es la longitud promedio (mm) de las estacas \(\alpha_j\): es el efecto producido por la j-ésima dosis de AIB sobre la longitud promedio (mm) \(\varepsilon_{ij}\): es el efecto del error aleatorio sobre la i-ésima estaca a la que se le aplicó la j-ésima dosis de AIB. Entonces, el modelo a aplicar es:
m_lm <- lm(Long.mm. ~ dosis, aib)
m_lm##
## Call:
## lm(formula = Long.mm. ~ dosis, data = aib)
##
## Coefficients:
## (Intercept) dosis0.5 dosis1.5
## 1.498 3.420 4.330Para la prueba de hipótesis del modelo, se pone a prueba la variabilidad entre las dosis (el efecto).
\[H_0:\alpha_j = 0\forall j\] \[H_1:\alpha_j \neq 0 algun j\]
m_aov <- aov(Long.mm.~ dosis, aib)
m_aov## Call:
## aov(formula = Long.mm. ~ dosis, data = aib)
##
## Terms:
## dosis Residuals
## Sum of Squares 62.54680 7.11645
## Deg. of Freedom 2 15
##
## Residual standard error: 0.6887888
## Estimated effects may be unbalancedResumen del modelo
summary(m_lm) ; summary.lm(m_aov)##
## Call:
## lm(formula = Long.mm. ~ dosis, data = aib)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.24833 -0.45833 0.03167 0.58417 0.97167
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.4983 0.2812 5.328 8.44e-05 ***
## dosis0.5 3.4200 0.3977 8.600 3.49e-07 ***
## dosis1.5 4.3300 0.3977 10.888 1.61e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6888 on 15 degrees of freedom
## Multiple R-squared: 0.8978, Adjusted R-squared: 0.8842
## F-statistic: 65.92 on 2 and 15 DF, p-value: 3.711e-08##
## Call:
## aov(formula = Long.mm. ~ dosis, data = aib)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.24833 -0.45833 0.03167 0.58417 0.97167
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.4983 0.2812 5.328 8.44e-05 ***
## dosis0.5 3.4200 0.3977 8.600 3.49e-07 ***
## dosis1.5 4.3300 0.3977 10.888 1.61e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6888 on 15 degrees of freedom
## Multiple R-squared: 0.8978, Adjusted R-squared: 0.8842
## F-statistic: 65.92 on 2 and 15 DF, p-value: 3.711e-08ADEVA
anova(m_lm) ; anova(m_aov)## Analysis of Variance Table
##
## Response: Long.mm.
## Df Sum Sq Mean Sq F value Pr(>F)
## dosis 2 62.547 31.2734 65.918 3.711e-08 ***
## Residuals 15 7.116 0.4744
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Analysis of Variance Table
##
## Response: Long.mm.
## Df Sum Sq Mean Sq F value Pr(>F)
## dosis 2 62.547 31.2734 65.918 3.711e-08 ***
## Residuals 15 7.116 0.4744
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Medias de tratamientos
En la lista del modelo y del ADEVA tenemos diferentes elementos para el análisis.
library(emmeans)
lsm <- emmeans(m_aov, ~ dosis)
lsm## dosis emmean SE df lower.CL upper.CL
## 0 1.50 0.281 15 0.899 2.10
## 0.5 4.92 0.281 15 4.319 5.52
## 1.5 5.83 0.281 15 5.229 6.43
##
## Confidence level used: 0.95plot(lsm)
# Supuestos
Residuos y predichos
fitted(m_aov)## 1 2 3 4 5 6 7 8
## 1.498333 1.498333 1.498333 1.498333 1.498333 1.498333 4.918333 4.918333
## 9 10 11 12 13 14 15 16
## 4.918333 4.918333 4.918333 4.918333 5.828333 5.828333 5.828333 5.828333
## 17 18
## 5.828333 5.828333resid(m_aov)## 1 2 3 4 5 6
## -1.09833333 -0.17833333 -0.47833333 0.97166667 0.62166667 0.16166667
## 7 8 9 10 11 12
## 0.48166667 0.67166667 -0.09833333 -1.24833333 0.59166667 -0.39833333
## 13 14 15 16 17 18
## 0.63166667 -0.39833333 -0.55833333 -0.51833333 0.28166667 0.56166667rstandard(m_aov)## 1 2 3 4 5 6 7
## -1.7467820 -0.2836201 -0.7607382 1.5453322 0.9886945 0.2571136 0.7660395
## 8 9 10 11 12 13 14
## 1.0682142 -0.1563887 -1.9853410 0.9409827 -0.6335067 1.0045985 -0.6335067
## 15 16 17 18
## -0.8879696 -0.8243539 0.4479608 0.8932709Combinar en el set de datos
aib <- aib %>%
mutate(est = fitted(m_aov), res_ord = resid(m_aov), res_std = rstandard(m_aov))
aib## Long.mm. dosis est res_ord res_std
## 1 0.40 0 1.498333 -1.09833333 -1.7467820
## 2 1.32 0 1.498333 -0.17833333 -0.2836201
## 3 1.02 0 1.498333 -0.47833333 -0.7607382
## 4 2.47 0 1.498333 0.97166667 1.5453322
## 5 2.12 0 1.498333 0.62166667 0.9886945
## 6 1.66 0 1.498333 0.16166667 0.2571136
## 7 5.40 0.5 4.918333 0.48166667 0.7660395
## 8 5.59 0.5 4.918333 0.67166667 1.0682142
## 9 4.82 0.5 4.918333 -0.09833333 -0.1563887
## 10 3.67 0.5 4.918333 -1.24833333 -1.9853410
## 11 5.51 0.5 4.918333 0.59166667 0.9409827
## 12 4.52 0.5 4.918333 -0.39833333 -0.6335067
## 13 6.46 1.5 5.828333 0.63166667 1.0045985
## 14 5.43 1.5 5.828333 -0.39833333 -0.6335067
## 15 5.27 1.5 5.828333 -0.55833333 -0.8879696
## 16 5.31 1.5 5.828333 -0.51833333 -0.8243539
## 17 6.11 1.5 5.828333 0.28166667 0.4479608
## 18 6.39 1.5 5.828333 0.56166667 0.8932709Análisis de suspuestos
Gráficos de diagnóstico
Si queremos obtener los cuatro gráficos del análisis para el diagnóstico en una sola salida, utilizamos del paquete base:
par(mfcol = c(2,2)) # para obtener todo en un grafico
plot(m_aov)
Comprobación de la Normalidad
plot(m_lm, which = 2)
ggplot(aib) +
aes(sample = res_ord) +
geom_qq() +
geom_qq_line() +
labs(x = "Cuantiles de residuos ordinarios", y = "Cuantiles teóricos normales")
Recordemos la prueba de hipótesis de Shapiro Wilks de bondad de ajuste a
la normalidad:
shapiro.test(resid(m_lm))##
## Shapiro-Wilk normality test
##
## data: resid(m_lm)
## W = 0.93609, p-value = 0.2481Homogeneidad
Apliquemos lo aprendido anteriormente:
plot(m_aov, which = 1)
ggplot(aib) +
aes(x = est, y = res_ord) +
geom_point() +
labs(x = "Predichos", y = "Residuos")
library(car)## Loading required package: carData##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recode## The following object is masked from 'package:purrr':
##
## someleveneTest(m_lm)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 2 0.1383 0.8719
## 15Comparaciones múltiples
Tukey
TukeyHSD
tukey1 <- TukeyHSD(m_aov, which = "dosis")
tukey1## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Long.mm. ~ dosis, data = aib)
##
## $dosis
## diff lwr upr p adj
## 0.5-0 3.42 2.3870578 4.452942 0.0000010
## 1.5-0 4.33 3.2970578 5.362942 0.0000000
## 1.5-0.5 0.91 -0.1229422 1.942942 0.0885165plot(tukey1, las = 1)
agricolae
HSD.test(m_lm, trt = "dosis", console = T)##
## Study: m_lm ~ "dosis"
##
## HSD Test for Long.mm.
##
## Mean Square Error: 0.47443
##
## dosis, means
##
## Long.mm. std r Min Max
## 0 1.498333 0.7514896 6 0.40 2.47
## 0.5 4.918333 0.7429513 6 3.67 5.59
## 1.5 5.828333 0.5536937 6 5.27 6.46
##
## Alpha: 0.05 ; DF Error: 15
## Critical Value of Studentized Range: 3.673378
##
## Minimun Significant Difference: 1.032942
##
## Treatments with the same letter are not significantly different.
##
## Long.mm. groups
## 1.5 5.828333 a
## 0.5 4.918333 a
## 0 1.498333 bEjemplo con más de un factor
Se desea conocer si existe un efecto conjunto entre el origen del tipo de pastura y la zona (en el centro oeste del país), en la producción de miel (en kg) por colmena. Se tomaron al azar 24 colmenas y se colocó reinas hermanas a todas por igual. Las colmenas fueron colocadas en 4 zonas diferentes y dentro de cada zona los terrenos linderos poseían tipos de cultivo nativo, introducido y mixto.
Datos
miel <- read_delim("miel.csv", delim = ";",
escape_double = FALSE, trim_ws = TRUE)## Rows: 23 Columns: 3
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ";"
## chr (2): zona, cultivo
## dbl (1): Miel_kg
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.head(miel)## # A tibble: 6 × 3
## Miel_kg zona cultivo
## <dbl> <chr> <chr>
## 1 75 A Nativo
## 2 72 A Nativo
## 3 62 B Nativo
## 4 64 B Nativo
## 5 85 C Nativo
## 6 88 C Nativomiel$zona = as.factor(miel$zona)
miel$cultivo =as.factor(miel$cultivo)
str(miel)## spec_tbl_df [23 × 3] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
## $ Miel_kg: num [1:23] 75 72 62 64 85 88 62 65 78 80 ...
## $ zona : Factor w/ 4 levels "A","B","C","D": 1 1 2 2 3 3 4 4 1 1 ...
## $ cultivo: Factor w/ 3 levels "Introducido",..: 3 3 3 3 3 3 3 3 1 1 ...
## - attr(*, "spec")=
## .. cols(
## .. Miel_kg = col_double(),
## .. zona = col_character(),
## .. cultivo = col_character()
## .. )
## - attr(*, "problems")=<externalptr>Descriptiva —-
miel %>%
group_by(zona,cultivo) %>%
summarise_all(.funs = c(
n = length,
media = mean,
sd = sd,
min = min,
max = max
)
)## # A tibble: 12 × 7
## # Groups: zona [4]
## zona cultivo n media sd min max
## <fct> <fct> <int> <dbl> <dbl> <dbl> <dbl>
## 1 A Introducido 2 79 1.41 78 80
## 2 A Mixto 2 86 1.41 85 87
## 3 A Nativo 2 73.5 2.12 72 75
## 4 B Introducido 2 71 1.41 70 72
## 5 B Mixto 2 79.5 2.12 78 81
## 6 B Nativo 2 63 1.41 62 64
## 7 C Introducido 2 92.5 0.707 92 93
## 8 C Mixto 2 99 1.41 98 100
## 9 C Nativo 2 86.5 2.12 85 88
## 10 D Introducido 2 76.5 2.12 75 78
## 11 D Mixto 1 92 NA 92 92
## 12 D Nativo 2 63.5 2.12 62 65Graficos interacción
ggplot(miel) +
aes(x = zona, y = Miel_kg,
color = cultivo,
group = cultivo) +
geom_jitter(width = 0.2) +
stat_summary(fun.data = "mean_se",
geom = "pointrange") +
geom_line(stat = "summary", fun = mean)## Warning: Removed 1 rows containing missing values (geom_segment).
ggplot(miel) +
aes(x = cultivo, y = Miel_kg,
color = zona,
group = zona) +
geom_jitter(width = 0.2) +
stat_summary(fun.data = "mean_se",
geom = "pointrange") +
geom_line(stat = "summary", fun = mean)## Warning: Removed 1 rows containing missing values (geom_segment).
Modelo lineal
El modelo que utilizaremos es de efectos fijos ya que hemos fijado los niveles de los factores, por lo tanto el modelo a comprobar es: \[Y_{ijl} =\mu + \alpha_j +\beta_l+(\alpha\beta)_{jl}+\varepsilon_{ijl}\] Para el análisis:
m_aov2 <- aov(Miel_kg ~ zona*cultivo,miel)
m_aov2## Call:
## aov(formula = Miel_kg ~ zona * cultivo, data = miel)
##
## Terms:
## zona cultivo zona:cultivo Residuals
## Sum of Squares 1586.3507 1012.3968 129.4699 33.0000
## Deg. of Freedom 3 2 6 11
##
## Residual standard error: 1.732051
## Estimated effects may be unbalancedADEVA
anova(m_aov2)## Analysis of Variance Table
##
## Response: Miel_kg
## Df Sum Sq Mean Sq F value Pr(>F)
## zona 3 1586.35 528.78 176.2612 1.404e-09 ***
## cultivo 2 1012.40 506.20 168.7328 5.569e-09 ***
## zona:cultivo 6 129.47 21.58 7.1928 0.002607 **
## Residuals 11 33.00 3.00
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Supuestos
Normalidad
m_aov3 <- aov(Miel_kg~zona:cultivo,miel)
shapiro.test(resid(m_aov3))##
## Shapiro-Wilk normality test
##
## data: resid(m_aov3)
## W = 0.82683, p-value = 0.001064Homoscedasticidad
leveneTest(m_aov2)## Warning in anova.lm(lm(resp ~ group)): ANOVA F-tests on an essentially perfect
## fit are unreliable## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 11 3.4373e+31 < 2.2e-16 ***
## 11
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1leveneTest(m_aov3)## Warning in anova.lm(lm(resp ~ group)): ANOVA F-tests on an essentially perfect
## fit are unreliable## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 11 3.4373e+31 < 2.2e-16 ***
## 11
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Medias —–
HSD.test(m_aov2, trt = "zona", console = T)##
## Study: m_aov2 ~ "zona"
##
## HSD Test for Miel_kg
##
## Mean Square Error: 3
##
## zona, means
##
## Miel_kg std r Min Max
## A 79.50000 5.753260 6 72 87
## B 71.16667 7.494442 6 62 81
## C 92.66667 5.715476 6 85 100
## D 74.40000 11.886968 5 62 92
##
## Alpha: 0.05 ; DF Error: 11
## Critical Value of Studentized Range: 4.256143
##
## Groups according to probability of means differences and alpha level( 0.05 )
##
## Treatments with the same letter are not significantly different.
##
## Miel_kg groups
## C 92.66667 a
## A 79.50000 b
## D 74.40000 c
## B 71.16667 dHSD.test(m_aov2, trt = "cultivo", console = T)##
## Study: m_aov2 ~ "cultivo"
##
## HSD Test for Miel_kg
##
## Mean Square Error: 3
##
## cultivo, means
##
## Miel_kg std r Min Max
## Introducido 79.75000 8.531454 8 70 93
## Mixto 88.71429 8.320943 7 78 100
## Nativo 71.62500 10.322479 8 62 88
##
## Alpha: 0.05 ; DF Error: 11
## Critical Value of Studentized Range: 3.819588
##
## Groups according to probability of means differences and alpha level( 0.05 )
##
## Treatments with the same letter are not significantly different.
##
## Miel_kg groups
## Mixto 88.71429 a
## Introducido 79.75000 b
## Nativo 71.62500 cHSD.test(m_aov2, trt = c("zona","cultivo"), console = T)##
## Study: m_aov2 ~ c("zona", "cultivo")
##
## HSD Test for Miel_kg
##
## Mean Square Error: 3
##
## zona:cultivo, means
##
## Miel_kg std r Min Max
## A:Introducido 79.0 1.4142136 2 78 80
## A:Mixto 86.0 1.4142136 2 85 87
## A:Nativo 73.5 2.1213203 2 72 75
## B:Introducido 71.0 1.4142136 2 70 72
## B:Mixto 79.5 2.1213203 2 78 81
## B:Nativo 63.0 1.4142136 2 62 64
## C:Introducido 92.5 0.7071068 2 92 93
## C:Mixto 99.0 1.4142136 2 98 100
## C:Nativo 86.5 2.1213203 2 85 88
## D:Introducido 76.5 2.1213203 2 75 78
## D:Mixto 92.0 NA 1 92 92
## D:Nativo 63.5 2.1213203 2 62 65
##
## Alpha: 0.05 ; DF Error: 11
## Critical Value of Studentized Range: 5.713009
##
## Groups according to probability of means differences and alpha level( 0.05 )
##
## Treatments with the same letter are not significantly different.
##
## Miel_kg groups
## C:Mixto 99.0 a
## C:Introducido 92.5 ab
## D:Mixto 92.0 ab
## C:Nativo 86.5 b
## A:Mixto 86.0 bc
## B:Mixto 79.5 cd
## A:Introducido 79.0 d
## D:Introducido 76.5 de
## A:Nativo 73.5 de
## B:Introducido 71.0 e
## D:Nativo 63.5 f
## B:Nativo 63.0 f