Introducción

library(readxl)
library(dplyr)
library(skimr)
datos <- read_excel("2- datos.xlsx")
datos %>% mutate(tratamiento = as.factor(tratamiento)) %>% select(2,3) -> datos
datos %>% skim()

Data summary
Name	Piped data
Number of rows	25
Number of columns	2
_______________________
Column type frequency:
factor	1
numeric	1
________________________
Group variables	None

Variable type: factor

skim_variable	n_missing	complete_rate	ordered	n_unique	top_counts
tratamiento	0	1	FALSE	5	T1: 5, T2: 5, T3: 5, T4: 5

Variable type: numeric

skim_variable	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
rendimiento	0	1	66.72	21.69	29	51	62	80	105	▃▇▇▅▆

Inferencia

datos %>% attach()
modelo <- lm(rendimiento~tratamiento, datos)

Análisis de varianza

modelo %>% aov %>% summary

##             Df Sum Sq Mean Sq F value   Pr(>F)    
## tratamiento  4  10277  2569.4    50.7 3.43e-10 ***
## Residuals   20   1014    50.7                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Comparación Múltiple

Sabemos que en nuestro análisis de varianza nos salió que al menos un efecto del tratamiento era distinto de cero. ´¿Cuál es el mejor tratamiento ahora?, para ello empleamos comparaciones múltiples para saber ello.

Prueba t

Es una prueba para comparaciones planeadas con anterioridad, es decir antes de recolectar los datos. El nivel de significación se toma como un error individual para cada prueba o también para cada tamaño de muestra . Los supuestos para la realización de esta prueba son:

La prueba F del análisis de variancia debe ser significativa.
Los \(\varepsilon_{ij}\) son variables aleatorias independientes y \(\varepsilon_{ij}\) ~ N(0,\(\sigma^2\))

Se busca comparar algunas de estas hipótesis

\(H_{0}: \mu_{i} \geq \mu_{j}\)

\(H_{1}: \mu_{i} < \mu_{j}\)

\(H_{0}: \mu_{i} \leq \mu_{j}\)

\(H_{1}: \mu_{i} > \mu_{j}\)

\(H_{0}: \mu_{i} = \mu_{j}\)

\(H_{1}: \mu_{i} \neq \mu_{j}\)

El estadístico de prueba es:

\(t_{calc} = \frac{y\bar{}_{i\bullet} - y\bar{}_{j\bullet}}{\sqrt{CME(\frac{1}{n_{i}} + \frac{1}{n_{j}})}}\) ~ \(t_{(N-k)}\)

Pero ya que estamos comparando por pares, es mejor emplear esta prueba:

\(H_{0}: \mu_{i} = \mu_{j}\)

\(H_{1}: \mu_{i} \neq \mu_{j}\)

Las medias del rendimiento de las mandarinas por cada tratamiento son:

(datos %>% group_by(tratamiento) %>% 
  summarise(Media = mean(rendimiento)) %>% select(Media) %>% as.matrix() %>% as.vector() -> Medias)

## [1] 101.0  59.6  76.6  53.0  43.4

Los cuadrados medios son:

(modelo %>% anova() %>% select(3) %>% as.matrix() %>%  as.vector()-> CM)

## [1] 2569.36   50.68

El estadístico de prueba si queremos comparar:

\(H_{0}: \mu_{1} = \mu_{2}\)

\(H_{1}: \mu_{1} \neq \mu_{2}\)

((Medias[1] - Medias[2])/sqrt(CM[2]*(2/5)) -> tc)

## [1] 9.195007

Valores Críticos

(qt(c(0.025,0.975),20) -> tcri)

## [1] -2.085963  2.085963

El p-valor

TENER EN CUENTA:

(2 * (1 - pt(abs(tc), 20)) -> pvalor)

## [1] 1.272565e-08

Decisión: Se rechaza \(H_{0}\)

Conclusión: A un nivel de significancia de 5% se puede concluir que las medias del tratamiento 1 y 2 son diferentes

Intervalos de confianza

LI <- (Medias[1] - Medias[2]) - qt(0.975,20)*sqrt(CM[2]*(1/5+1/5))
LS <- (Medias[1] - Medias[2]) + qt(0.975,20)*sqrt(CM[2]*(1/5+1/5))
c(LI,LS)

## [1] 32.00807 50.79193

Tener en cuenta que: Si el intervalo contiene al 0, se puede decir que no habría diferencia significativa entre las medias evaluadas.

Esta manera realizada es una opción y repetiríamos media por media de cada tratamiento, es decir, \(\mu_{1}\) con \(\mu_{2}\),… así \(\binom{2}{k}\), donde k = Número de tratamientos.

Una manera más sencilla es usar la librería agricolae de la siguiente manera:

library(agricolae)
LSD.test(modelo, "tratamiento", alpha = 0.05, p.adj = "none", console = T, group = F)

## 
## Study: modelo ~ "tratamiento"
## 
## LSD t Test for rendimiento 
## 
## Mean Square Error:  50.68 
## 
## tratamiento,  means and individual ( 95 %) CI
## 
##    rendimiento       std r       se     LCL      UCL Min Max Q25 Q50 Q75
## T1       101.0  2.915476 5 3.183709 94.3589 107.6411  97 105 100 101 102
## T2        59.6  4.277850 5 3.183709 52.9589  66.2411  54  65  57  60  62
## T3        76.6  5.683309 5 3.183709 69.9589  83.2411  68  83  75  77  80
## T4        53.0 10.770330 5 3.183709 46.3589  59.6411  40  67  46  52  60
## T5        43.4  8.848729 5 3.183709 36.7589  50.0411  29  51  42  45  50
## 
## Alpha: 0.05 ; DF Error: 20
## Critical Value of t: 2.085963 
## 
## Comparison between treatments means
## 
##         difference pvalue signif.         LCL       UCL
## T1 - T2       41.4 0.0000     ***  32.0080668 50.791933
## T1 - T3       24.4 0.0000     ***  15.0080668 33.791933
## T1 - T4       48.0 0.0000     ***  38.6080668 57.391933
## T1 - T5       57.6 0.0000     ***  48.2080668 66.991933
## T2 - T3      -17.0 0.0012      ** -26.3919332 -7.608067
## T2 - T4        6.6 0.1582          -2.7919332 15.991933
## T2 - T5       16.2 0.0018      **   6.8080668 25.591933
## T3 - T4       23.6 0.0000     ***  14.2080668 32.991933
## T3 - T5       33.2 0.0000     ***  23.8080668 42.591933
## T4 - T5        9.6 0.0456       *   0.2080668 18.991933

LSD.test(modelo, "tratamiento", alpha = 0.05, p.adj = "none", console = T) %>% plot()

## 
## Study: modelo ~ "tratamiento"
## 
## LSD t Test for rendimiento 
## 
## Mean Square Error:  50.68 
## 
## tratamiento,  means and individual ( 95 %) CI
## 
##    rendimiento       std r       se     LCL      UCL Min Max Q25 Q50 Q75
## T1       101.0  2.915476 5 3.183709 94.3589 107.6411  97 105 100 101 102
## T2        59.6  4.277850 5 3.183709 52.9589  66.2411  54  65  57  60  62
## T3        76.6  5.683309 5 3.183709 69.9589  83.2411  68  83  75  77  80
## T4        53.0 10.770330 5 3.183709 46.3589  59.6411  40  67  46  52  60
## T5        43.4  8.848729 5 3.183709 36.7589  50.0411  29  51  42  45  50
## 
## Alpha: 0.05 ; DF Error: 20
## Critical Value of t: 2.085963 
## 
## least Significant Difference: 9.391933 
## 
## Treatments with the same letter are not significantly different.
## 
##    rendimiento groups
## T1       101.0      a
## T3        76.6      b
## T2        59.6      c
## T4        53.0      c
## T5        43.4      d

Sabemos que en la prueba “t” la P(error tipo I) = \(\alpha\) es por cada prueba o cada par, por ello la probabilidad de que la prueba esté bien es 1-\(\alpha\), siguiendo esa idea, entonces como hemos realizado 10 comparaciones, la probabilidad de que todos nuestros resultados sean correctos son:

\((1-\alpha)^\binom{5}{2}\)

(1-0.05)**choose(5,2)

## [1] 0.5987369

Aproximadamente la probabilidad de que todos nuestros resultados en conjunto esten correctos es 0.6 = 60%, ello es muy bajo, por ello usaremos la corrección de Bonferroni. Pero antes de ello tengamos en cuenta lo siguiente:

Corrección de bonferroni

\(H_{0}: \mu_{1} = \mu_{2}\)

\(H_{1}: \mu_{1} \neq \mu_{2}\)

El \(t_{cal}\) es el mismo

tc

## [1] 9.195007

El “t” crítico es:

\(t_{cri} = t_{(1-\frac{\alpha}{2m},N-t)}\)

m <- choose(5,2)
(tcri_bonf <- qt(1-0.05/(2*m), 25-5))

## [1] 3.153401

El pvalor es:

m x 2 x \(P(\left | t_{calc}\right| > t_{(N-t)})\)

min(m * 2 * (1 - pt(abs(tc), 20)))

## [1] 1.272565e-07

Si consideramos el primer par que evaluamos entre t1 y t2, sería:

p.adjust(pvalor, method = "bonferroni", n = m)

## [1] 1.272565e-07

El intervalo de confianza a un 95% es:

LI_bonf <- (Medias[1] - Medias[2]) - qt(1-0.05/(2*m),20)*sqrt(CM[2]*(1/5+1/5))
LS_bonf <- (Medias[1] - Medias[2]) + qt(1-0.05/(2*m),20)*sqrt(CM[2]*(1/5+1/5))
c(LI_bonf,LS_bonf)

## [1] 27.20199 55.59801

Podemos observar que bonferroni es una corrección más justa.

Utilizando el paquete agricolae

LSD.test(modelo, "tratamiento", alpha = 0.05, p.adj = "bonferroni", group = F, console = T)

## 
## Study: modelo ~ "tratamiento"
## 
## LSD t Test for rendimiento 
## P value adjustment method: bonferroni 
## 
## Mean Square Error:  50.68 
## 
## tratamiento,  means and individual ( 95 %) CI
## 
##    rendimiento       std r       se     LCL      UCL Min Max Q25 Q50 Q75
## T1       101.0  2.915476 5 3.183709 94.3589 107.6411  97 105 100 101 102
## T2        59.6  4.277850 5 3.183709 52.9589  66.2411  54  65  57  60  62
## T3        76.6  5.683309 5 3.183709 69.9589  83.2411  68  83  75  77  80
## T4        53.0 10.770330 5 3.183709 46.3589  59.6411  40  67  46  52  60
## T5        43.4  8.848729 5 3.183709 36.7589  50.0411  29  51  42  45  50
## 
## Alpha: 0.05 ; DF Error: 20
## Critical Value of t: 3.153401 
## 
## Comparison between treatments means
## 
##         difference pvalue signif.        LCL       UCL
## T1 - T2       41.4 0.0000     ***  27.201991 55.598009
## T1 - T3       24.4 0.0003     ***  10.201991 38.598009
## T1 - T4       48.0 0.0000     ***  33.801991 62.198009
## T1 - T5       57.6 0.0000     ***  43.401991 71.798009
## T2 - T3      -17.0 0.0119       * -31.198009 -2.801991
## T2 - T4        6.6 1.0000          -7.598009 20.798009
## T2 - T5       16.2 0.0180       *   2.001991 30.398009
## T3 - T4       23.6 0.0004     ***   9.401991 37.798009
## T3 - T5       33.2 0.0000     ***  19.001991 47.398009
## T4 - T5        9.6 0.4558          -4.598009 23.798009

LSD.test(modelo, "tratamiento", alpha = 0.05, p.adj = "bonferroni", console = T) %>% plot()

## 
## Study: modelo ~ "tratamiento"
## 
## LSD t Test for rendimiento 
## P value adjustment method: bonferroni 
## 
## Mean Square Error:  50.68 
## 
## tratamiento,  means and individual ( 95 %) CI
## 
##    rendimiento       std r       se     LCL      UCL Min Max Q25 Q50 Q75
## T1       101.0  2.915476 5 3.183709 94.3589 107.6411  97 105 100 101 102
## T2        59.6  4.277850 5 3.183709 52.9589  66.2411  54  65  57  60  62
## T3        76.6  5.683309 5 3.183709 69.9589  83.2411  68  83  75  77  80
## T4        53.0 10.770330 5 3.183709 46.3589  59.6411  40  67  46  52  60
## T5        43.4  8.848729 5 3.183709 36.7589  50.0411  29  51  42  45  50
## 
## Alpha: 0.05 ; DF Error: 20
## Critical Value of t: 3.153401 
## 
## Minimum Significant Difference: 14.19801 
## 
## Treatments with the same letter are not significantly different.
## 
##    rendimiento groups
## T1       101.0      a
## T3        76.6      b
## T2        59.6      c
## T4        53.0     cd
## T5        43.4      d

Con bonferroni la tasa de error de todos en conjunto es de 0.05 = 5%, es un resultado mucho más confiable, respecto a la tradicional.

Método de diferencia de contrastes LSD o DLS

Es una forma abreviada de la prueba “t” del caso bilateral donde LSD (Least Significant Difference) es el criterio de decisión significativa de fisher, de modo que, cualquier diferencia entre las medias de dos tratamientos mayor a dicho límite sea significativa, esta prueba debe ser planeada con anterioridad.

Diferencia Límite Significativa:

\(LSD = t_{(1-\frac{\alpha}{2})}\sqrt{CME(\frac{1}{n_{i}}+\frac{1}{nj})}\)

Regla de decisión:

Se rechaza la \(H_{o}\) si:

\(\left |\bar{Y}_{i\bullet} - \bar{Y}_{j\bullet}\right| > LSD\)

LSD <- qt(0.975,20)*sqrt(CM[2]*(2/5))
Dif <- abs(Medias[1]-Medias[2])

if (Dif > LSD){
  print("Se rechaza H0")
}else{
  print("No se rechaza H0")
}

## [1] "Se rechaza H0"

Decisión: Se rechaza H0

Conclusión: A un nivel de significancia del 5%, se puede afirmar que las diferencias de las medias de tratamiento 1 con el tratamiento 2 es diferente de 0

De la misma manera se puede crear el IC.

Si queremos hallar el LSD con la corrección de bonferroni sería lo siguiente

Diferencia Límite Significativa:

\(LSD = t_{(1-\frac{\alpha}{2m})}\sqrt{CME(\frac{1}{n_{i}}+\frac{1}{nj})}\)

LSD_bonf <- qt(1-0.05/(2*m),20)*sqrt(CM[2]*(2/5))

if (Dif > LSD_bonf){
  print("Se rechaza H0")
}else{
  print("No se rechaza H0")
}

## [1] "Se rechaza H0"

En conclusión : Según la prueba “t” con la corrección de bonferroni el mejor tratamiento es el tratamiento 1, ya que se obtiene el mayor rendimiento de las mandarinas.

Prueba de Tukey

Usa la distribución de rangos estudentizados; compara todas las medias en pares considerando un tamaño de prueba \(\alpha\) si se cuenta en un experimento con “t” tratamientos el número total de hipótesis a analizar es obtenido mediante la combinatoria de “t” en 2, para aplicar esta prueba se debe comprobar la homocedasticidad de los errores, a diferencia de la prueba “t” NO necesita planear que grupo en específico se quiere comparar, porque tukey compara todas las pruebas sin excepción y mantiene el error global de tipo I (\(\alpha\) = 0.05), a diferencia de la “t”, que mientras más comparaciones hagas más se incrementa el error tipo I, para ello es útil la corrección de bonferroni. Por ello se dice que se considera la prueba tukey un error por familia (también llamado error por comparación múltiple o error Tipo I acumulado).

Hipótesis:

\(H_{0}: \mu_{i} = \mu_{j}\) \(\forall i \neq j\) j = 1,2,3,..,t

\(H_{1}: \mu_{i} \neq \mu_{j}\)

Estadístico de prueba:

\(Q_{calc} = \frac{\bar{y}_{max} - \bar{y}_{min}}{\sqrt{CME/2(\frac{1}{n_{i}} + \frac{1}{n_{j}})}}\) ~ \(q_{(t,N-t)}\)

Amplitud significativa de Tukey:

\(ALS(T) = q_{(1-\alpha,t,N-t)}\sqrt{\frac{CME}{2}(\frac{1}{n_{i}} + \frac{1}{n_{j}})}\)

Regla de decisión:

Se rechaza la \(H_{0}\) si:

\(\left |\bar{Y}_{i\bullet} - \bar{Y}_{j\bullet}\right| > ALS(T)\)

Qt <- abs(Medias[1] - Medias[2])/sqrt(CM[2]/2*(1/5+1/5))

qtukey(0.95,5,20)

## [1] 4.231857

A simple vista como Qt > cuartil tukey, Se rechaza la hipótesis nula

Según regla de decisión:

abs(Medias[1] - Medias[2])

## [1] 41.4

(HSD <- qtukey(0.95,5,20)*sqrt(CM[2]/2*(1/5 + 1/5)))

## [1] 13.473

Como 41.4 > 13.473, se rechaza la hipótesis nula, a un nivel de significancia del 5%, se puede afirmar que las diferencias de las medias de tratamiento 1 con el tratamiento 2 es diferente de 0, es decir es significativo

Usando la librería agricolae

HSD.test(modelo, "tratamiento", alpha = 0.05, group = F, console = T)

## 
## Study: modelo ~ "tratamiento"
## 
## HSD Test for rendimiento 
## 
## Mean Square Error:  50.68 
## 
## tratamiento,  means
## 
##    rendimiento       std r       se Min Max Q25 Q50 Q75
## T1       101.0  2.915476 5 3.183709  97 105 100 101 102
## T2        59.6  4.277850 5 3.183709  54  65  57  60  62
## T3        76.6  5.683309 5 3.183709  68  83  75  77  80
## T4        53.0 10.770330 5 3.183709  40  67  46  52  60
## T5        43.4  8.848729 5 3.183709  29  51  42  45  50
## 
## Alpha: 0.05 ; DF Error: 20 
## Critical Value of Studentized Range: 4.231857 
## 
## Comparison between treatments means
## 
##         difference pvalue signif.        LCL       UCL
## T1 - T2       41.4 0.0000     ***  27.927002 54.872998
## T1 - T3       24.4 0.0002     ***  10.927002 37.872998
## T1 - T4       48.0 0.0000     ***  34.527002 61.472998
## T1 - T5       57.6 0.0000     ***  44.127002 71.072998
## T2 - T3      -17.0 0.0093      ** -30.472998 -3.527002
## T2 - T4        6.6 0.5948          -6.872998 20.072998
## T2 - T5       16.2 0.0138       *   2.727002 29.672998
## T3 - T4       23.6 0.0003     ***  10.127002 37.072998
## T3 - T5       33.2 0.0000     ***  19.727002 46.672998
## T4 - T5        9.6 0.2456          -3.872998 23.072998

HSD.test(modelo, "tratamiento", alpha = 0.05, group = T, console = T)

## 
## Study: modelo ~ "tratamiento"
## 
## HSD Test for rendimiento 
## 
## Mean Square Error:  50.68 
## 
## tratamiento,  means
## 
##    rendimiento       std r       se Min Max Q25 Q50 Q75
## T1       101.0  2.915476 5 3.183709  97 105 100 101 102
## T2        59.6  4.277850 5 3.183709  54  65  57  60  62
## T3        76.6  5.683309 5 3.183709  68  83  75  77  80
## T4        53.0 10.770330 5 3.183709  40  67  46  52  60
## T5        43.4  8.848729 5 3.183709  29  51  42  45  50
## 
## Alpha: 0.05 ; DF Error: 20 
## Critical Value of Studentized Range: 4.231857 
## 
## Minimun Significant Difference: 13.473 
## 
## Treatments with the same letter are not significantly different.
## 
##    rendimiento groups
## T1       101.0      a
## T3        76.6      b
## T2        59.6      c
## T4        53.0     cd
## T5        43.4      d

En conclusión: El mejor tratamiento según “t” y Tukey sigue siendo el tratamiento 1, entendiendo como funciona cada prueba y los métodos que se emplean para determinar que tratamiento es el mejor.

Prueba de Student Newman Keuls

Se requiere que las medias estén ordenadas de forma ascendente para así evaluar las medias por posiciones y saber si son iguales o diferentes, la prueba NKS es una prueba tukey por posiciones, ello es importante, porque podemos ahorrar pruebas, usando la NKS, a comparación de tukey que aplica la prueba de manera global.

Buscamos contrastar las hipótesis de comparación de medias en pares, para \(i\neq j\)

\(H_{0} : \mu_{(i)} = \mu_{(j)}\)

\(H_{1} : \mu_{(i)} \neq \mu_{(j)}\)

Criterio de decisión:

\(SNK = q_{(1-\alpha,h,N-t)}\sqrt{\frac{CME}{t}(\frac{1}{n_{1}} + \frac{1}{n_{2}} + ... \frac{1}{n_{t}})}\)

Medias

## [1] 101.0  59.6  76.6  53.0  43.4

u5 < u4 < u2 < u3 < u1 (ordenado de forma ascendente)

Vamos a comparar la media 1 con 2, por ello el h = 5-3+1 = 3

qtukey(0.95,3,20)*sqrt(CM[2]/5*(1/5+1/5+1/5+1/5+1/5))

## [1] 11.3911

Medias[1] - Medias[2]

## [1] 41.4

Como:

41.4 > 11.39 -> se rechaza la \(H_{0}\), por lo tanto existe diferencia significativas entre las media 1 y media 2.

Valor crítico:

(ValorCritico <- abs(Medias[1] - Medias[2])/sqrt(CM[2]/5*(1/5+1/5+1/5+1/5+1/5)))

## [1] 13.0037

ptukey(ValorCritico, 3, 20, lower.tail = F)

## [1] 3.720337e-08

1 - ptukey(ValorCritico, 3, 20)

## [1] 3.720337e-08

Ahora supongamos que queramos comparar la media 1 con media 5; el h = 5-1+1 = 5

qtukey(0.95,5,20)*sqrt(CM[2]/5*(1/5+1/5+1/5+1/5+1/5))

## [1] 13.473

Medias[1] - Medias[5]

## [1] 57.6

Igual se rechaza \(H_{0}\) ello es suponible, ya que hallamos una diferencia más angosta que es significativa (u1 con u2) y al evaluar una media más ancha es evidente que también saldrá que es significativo, pero aquí entremos en debate, hemos cojido la media máxima que es (u1) con la mínima (u5), ¿Ello no es la prueba tukey?; pues sí, entonces concluímos que la prueba tukey es más particular y evita estar evaluando, pero la Newman hace menos comparaciones, las dos son viables.

Se puede realizar el mismo procedimiento con las demás medias, y verificar con los resultados que arroja R con el comando SNK.test

SNK.test(modelo, "tratamiento", group = F, alpha = 0.05, console = T)

## 
## Study: modelo ~ "tratamiento"
## 
## Student Newman Keuls Test
## for rendimiento 
## 
## Mean Square Error:  50.68 
## 
## tratamiento,  means
## 
##    rendimiento       std r       se Min Max Q25 Q50 Q75
## T1       101.0  2.915476 5 3.183709  97 105 100 101 102
## T2        59.6  4.277850 5 3.183709  54  65  57  60  62
## T3        76.6  5.683309 5 3.183709  68  83  75  77  80
## T4        53.0 10.770330 5 3.183709  40  67  46  52  60
## T5        43.4  8.848729 5 3.183709  29  51  42  45  50
## 
## Comparison between treatments means
## 
##         difference pvalue signif.         LCL       UCL
## T1 - T2       41.4 0.0000     ***  30.0088992 52.791101
## T1 - T3       24.4 0.0000     ***  15.0080673 33.791933
## T1 - T4       48.0 0.0000     ***  35.3979473 60.602053
## T1 - T5       57.6 0.0000     ***  44.1270018 71.072998
## T2 - T3      -17.0 0.0012      ** -26.3919327 -7.608067
## T2 - T4        6.6 0.1582          -2.7919327 15.991933
## T2 - T5       16.2 0.0049      **   4.8088992 27.591101
## T3 - T4       23.6 0.0001     ***  12.2088992 34.991101
## T3 - T5       33.2 0.0000     ***  20.5979473 45.802053
## T4 - T5        9.6 0.0456       *   0.2080673 18.991933

SNK.test(modelo, "tratamiento", group = T, alpha = 0.05, console = T) %>% plot(variation = "IQR")

## 
## Study: modelo ~ "tratamiento"
## 
## Student Newman Keuls Test
## for rendimiento 
## 
## Mean Square Error:  50.68 
## 
## tratamiento,  means
## 
##    rendimiento       std r       se Min Max Q25 Q50 Q75
## T1       101.0  2.915476 5 3.183709  97 105 100 101 102
## T2        59.6  4.277850 5 3.183709  54  65  57  60  62
## T3        76.6  5.683309 5 3.183709  68  83  75  77  80
## T4        53.0 10.770330 5 3.183709  40  67  46  52  60
## T5        43.4  8.848729 5 3.183709  29  51  42  45  50
## 
## Alpha: 0.05 ; DF Error: 20 
## 
## Critical Range
##         2         3         4         5 
##  9.391933 11.391101 12.602053 13.472998 
## 
## Means with the same letter are not significantly different.
## 
##    rendimiento groups
## T1       101.0      a
## T3        76.6      b
## T2        59.6      c
## T4        53.0      c
## T5        43.4      d

Haciendo la prueba de SNK, también llegamos a la conclusión que el mejor tratamiento que hace que nuestro rendimiento crezca más es el tratamiento 1, hasta ahora parece ser la ganadora.

datos %>% group_by(tratamiento) %>% skim()

Data summary
Name	Piped data
Number of rows	25
Number of columns	2
_______________________
Column type frequency:
numeric	1
________________________
Group variables	tratamiento

Variable type: numeric

skim_variable	tratamiento	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
rendimiento	T1	1	101.0	2.92	97	100	101	102	105	▇▇▇▇▇
rendimiento	T2	1	59.6	4.28	54	57	60	62	65	▇▇▇▇▇
rendimiento	T3	1	76.6	5.68	68	75	77	80	83	▃▁▇▃▃
rendimiento	T4	1	53.0	10.77	40	46	52	60	67	▇▇▇▇▇
rendimiento	T5	1	43.4	8.85	29	42	45	50	51	▃▁▃▃▇

Prueba de Duncan

Buscamos contrastar las hipótesis de comparación de medias en pares, para \(i\neq j\)

\(H_{0} : \mu_{(i)} = \mu_{(j)}\)

\(H_{1} : \mu_{(i)} \neq \mu_{(j)}\)

Criterio de decisión:

\(DU = q_{(\alpha´,h,N-t)}\sqrt{\frac{CME}{t}(\frac{1}{n_{1}} + \frac{1}{n_{2}} + ... \frac{1}{n_{t}})}\)

Para determinar el \(\alpha^´\) = (1-\(\alpha^{h-1}\)), donde “h” es lo mismo que habíamos hecho para la prueba de SNK.

Vamos a comparar la media 1 con 2, por ello el h = 5-3+1 = 3, entonces \(\alpha^{'}\) = 1 - (1 - \(\alpha)^{h-1}\); \(\alpha^{'}\) = 1 - \(0.95^2\) = 0.0975 = 0.0975

qtukey(1-0.0975,3,20)*sqrt(CM[2]/5*(1/5+1/5+1/5+1/5+1/5))

## [1] 9.858374

Medias[1] - Medias[2]

## [1] 41.4

Como:

41.4 > 8.97 -> se rechaza la \(H_{0}\), por lo tanto existe diferencia significativas entre las media 1 y media 2.

Valor crítico:

(ValorCritico <- abs(Medias[1] - Medias[2])/sqrt(CM[2]/5*(1/5+1/5+1/5+1/5+1/5)))

## [1] 13.0037

ptukey(ValorCritico, 3, 20, lower.tail = F)

## [1] 3.720337e-08

1 - ptukey(ValorCritico, 3, 20)

## [1] 3.720337e-08

Ahora supongamos que queramos comparar la media 1 con media 5; el h = 5-1+1 = 5, entonces \(\alpha^{'}\) = 1 - (1 - \(\alpha)^{h-1}\); \(\alpha^{'}\) = 1 - \(0.95^4\) = 0.1854938 = 0.19

qtukey(1-0.1854938,5,20)*sqrt(CM[2]/5*(1/5+1/5+1/5+1/5+1/5))

## [1] 10.36185

Medias[1] - Medias[5]

## [1] 57.6

Como:

57.6 > 10.361 -> se rechaza la \(H_{0}\), por lo tanto existe diferencia significativas entre las media 1 y media 5.

Se puede realizar el mismo procedimiento con las demás medias, y verificar con los resultados que arroja R con el comando duncan.test

duncan.test(modelo, "tratamiento", alpha = 0.05, group = F, console = T)

## 
## Study: modelo ~ "tratamiento"
## 
## Duncan's new multiple range test
## for rendimiento 
## 
## Mean Square Error:  50.68 
## 
## tratamiento,  means
## 
##    rendimiento       std r       se Min Max Q25 Q50 Q75
## T1       101.0  2.915476 5 3.183709  97 105 100 101 102
## T2        59.6  4.277850 5 3.183709  54  65  57  60  62
## T3        76.6  5.683309 5 3.183709  68  83  75  77  80
## T4        53.0 10.770330 5 3.183709  40  67  46  52  60
## T5        43.4  8.848729 5 3.183709  29  51  42  45  50
## 
## Comparison between treatments means
## 
##         difference pvalue signif.         LCL       UCL
## T1 - T2       41.4 0.0000     ***  31.5416261 51.258374
## T1 - T3       24.4 0.0000     ***  15.0080673 33.791933
## T1 - T4       48.0 0.0000     ***  37.8451925 58.154807
## T1 - T5       57.6 0.0000     ***  47.2381501 67.961850
## T2 - T3      -17.0 0.0012      ** -26.3919327 -7.608067
## T2 - T4        6.6 0.1582          -2.7919327 15.991933
## T2 - T5       16.2 0.0024      **   6.3416261 26.058374
## T3 - T4       23.6 0.0001     ***  13.7416261 33.458374
## T3 - T5       33.2 0.0000     ***  23.0451925 43.354807
## T4 - T5        9.6 0.0456       *   0.2080673 18.991933

duncan.test(modelo, "tratamiento", alpha = 0.05, group = T, console = T) %>% plot(variation = "SD")

## 
## Study: modelo ~ "tratamiento"
## 
## Duncan's new multiple range test
## for rendimiento 
## 
## Mean Square Error:  50.68 
## 
## tratamiento,  means
## 
##    rendimiento       std r       se Min Max Q25 Q50 Q75
## T1       101.0  2.915476 5 3.183709  97 105 100 101 102
## T2        59.6  4.277850 5 3.183709  54  65  57  60  62
## T3        76.6  5.683309 5 3.183709  68  83  75  77  80
## T4        53.0 10.770330 5 3.183709  40  67  46  52  60
## T5        43.4  8.848729 5 3.183709  29  51  42  45  50
## 
## Alpha: 0.05 ; DF Error: 20 
## 
## Critical Range
##         2         3         4         5 
##  9.391933  9.858374 10.154807 10.361850 
## 
## Means with the same letter are not significantly different.
## 
##    rendimiento groups
## T1       101.0      a
## T3        76.6      b
## T2        59.6      c
## T4        53.0      c
## T5        43.4      d

Conclusión: Todas las pruebas hasta ahora coinciden en que el tratamiento 1 es el mejor.

Prueba de Dunnet

Permite la comparación de las medias en pares respecto a un tratamiento control,testigo o placebo, correspondiente a un tratamiento cuyo efecto se conoce de antemano o corresponde a un método o procedimiento estándar que se utiliza de manera tradicional.

uscamos contrastar las hipótesis de comparación de medias en pares, para \(i\neq j\)

\(H_{0} : \mu_{i} = \mu_{0}\)

\(H_{1} : \mu_{i} \neq \mu_{0}\)

para todo i diferente 0

Criterio de decisión:

\(Dn = q_{(1-\alpha/2,t-1,N-t)}\sqrt{\frac{CME}{t}(\frac{1}{n_{i}} + \frac{1}{n_{0}}}\)

donde D1−α/2,t−1,GLE es el cuantil 1−α/2 de la tabla Dunnett para las t−1 comparaciones de tratamientos y GLE grados de libertad. Así, se tiene:

Con paquete multcomp; pongamos de ejemplo que nuestro tratamiento testigo va ser el tratamiento 2

datos1 <- data.frame(tratamiento = tratamiento, rendimiento = rendimiento
)
datos1$tratamiento <- factor(datos1$tratamiento, levels = c("T2","T1","T3","T4","T5"))
modelo_<- lm(rendimiento~tratamiento, datos1)
library(multcomp)
dunnet = glht(modelo_, linfct = mcp(tratamiento = "Dunnett"))
dunnet

## 
##   General Linear Hypotheses
## 
## Multiple Comparisons of Means: Dunnett Contrasts
## 
## 
## Linear Hypotheses:
##              Estimate
## T1 - T2 == 0     41.4
## T3 - T2 == 0     17.0
## T4 - T2 == 0     -6.6
## T5 - T2 == 0    -16.2

dunnet %>% summary

## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Dunnett Contrasts
## 
## 
## Fit: lm(formula = rendimiento ~ tratamiento, data = datos1)
## 
## Linear Hypotheses:
##              Estimate Std. Error t value Pr(>|t|)    
## T1 - T2 == 0   41.400      4.502   9.195  < 0.001 ***
## T3 - T2 == 0   17.000      4.502   3.776  0.00418 ** 
## T4 - T2 == 0   -6.600      4.502  -1.466  0.41229    
## T5 - T2 == 0  -16.200      4.502  -3.598  0.00635 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

dunnet %>% plot()

Con paquete DescTools:

library(DescTools)
DunnettTest(datos$rendimiento, datos$tratamiento, control = "T2")

## 
##   Dunnett's test for comparing several treatments with a control :  
##     95% family-wise confidence level
## 
## $T2
##        diff     lwr.ci    upr.ci    pval    
## T1-T2  41.4  29.457341 53.342659 8.8e-09 ***
## T3-T2  17.0   5.057341 28.942659  0.0043 ** 
## T4-T2  -6.6 -18.542659  5.342659  0.4122    
## T5-T2 -16.2 -28.142659 -4.257341  0.0064 ** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Conclusión: Tanto todas las pruebas coincide que el mejor tratamiento es el tratamiento 1.

Prueba de contrastes (Scheffé)

Este método está diseñado para probar todos los contrastes de medias que pudieran interesar al experimentador. Esto quiere decir, que se puede realizar infinitas pruebas simultáneas, aunque en la práctica se realice un número finito de pruebas simultáneas.

\(H_{0}: \sum_{i = 1}^{t} c_{i}\mu _{i} = 0\)

\(H_{1}: \sum_{i = 1}^{t} c_{i}\mu _{i} \neq 0\)

Nivel de significación \(\alpha\)

Valor Crítico de la prueba

Pongamos de ejemplo que queremos hacer la siguiente hipótesis:

\(H_{0}: \mu_{1} + \mu_{2} = \mu_{3} + \mu_{4} + \mu_{5}\)

\(H_{1}: \mu_{1} + \mu_{2} \neq \mu_{3} + \mu_{4} + \mu_{5}\)

(ybar <- tapply(datos$rendimiento,datos$tratamiento, mean))

##    T1    T2    T3    T4    T5 
## 101.0  59.6  76.6  53.0  43.4

ci <- c(1,1,-1,-1,-1)

ni<-table(datos$tratamiento)

ce<-abs(sum(ci*ybar))

anva <- modelo %>% anova
cme<-anva[2,3]

VCS<-sqrt(cme*sum(ci^2/ni))*sqrt(anva[1,1]*qf(0.95,anva[1,1],anva[2,1]))

ce > VCS

## [1] FALSE

Por lo tanto No se rechaza \(H_{0}\): El tratamiento 1 en conjunto con el tratamiento 2, es significativamente igual a las medias de los demas tratamiento en conjunto.

scheffe.test(modelo, "tratamiento", alpha = 0.05, group = F, console = T)

## 
## Study: modelo ~ "tratamiento"
## 
## Scheffe Test for rendimiento 
## 
## Mean Square Error  : 50.68 
## 
## tratamiento,  means
## 
##    rendimiento       std r       se Min Max Q25 Q50 Q75
## T1       101.0  2.915476 5 3.183709  97 105 100 101 102
## T2        59.6  4.277850 5 3.183709  54  65  57  60  62
## T3        76.6  5.683309 5 3.183709  68  83  75  77  80
## T4        53.0 10.770330 5 3.183709  40  67  46  52  60
## T5        43.4  8.848729 5 3.183709  29  51  42  45  50
## 
## Alpha: 0.05 ; DF Error: 20 
## Critical Value of F: 2.866081 
## 
## Comparison between treatments means
## 
##         Difference pvalue sig         LCL       UCL
## T1 - T2       41.4 0.0000 ***  26.1551711 56.644829
## T1 - T3       24.4 0.0008 ***   9.1551711 39.644829
## T1 - T4       48.0 0.0000 ***  32.7551711 63.244829
## T1 - T5       57.6 0.0000 ***  42.3551711 72.844829
## T2 - T3      -17.0 0.0238   * -32.2448289 -1.755171
## T2 - T4        6.6 0.7101      -8.6448289 21.844829
## T2 - T5       16.2 0.0335   *   0.9551711 31.444829
## T3 - T4       23.6 0.0012  **   8.3551711 38.844829
## T3 - T5       33.2 0.0000 ***  17.9551711 48.444829
## T4 - T5        9.6 0.3678      -5.6448289 24.844829

scheffe.test(modelo, "tratamiento", alpha = 0.05, group = T, console = T) %>% plot()

## 
## Study: modelo ~ "tratamiento"
## 
## Scheffe Test for rendimiento 
## 
## Mean Square Error  : 50.68 
## 
## tratamiento,  means
## 
##    rendimiento       std r       se Min Max Q25 Q50 Q75
## T1       101.0  2.915476 5 3.183709  97 105 100 101 102
## T2        59.6  4.277850 5 3.183709  54  65  57  60  62
## T3        76.6  5.683309 5 3.183709  68  83  75  77  80
## T4        53.0 10.770330 5 3.183709  40  67  46  52  60
## T5        43.4  8.848729 5 3.183709  29  51  42  45  50
## 
## Alpha: 0.05 ; DF Error: 20 
## Critical Value of F: 2.866081 
## 
## Minimum Significant Difference: 15.24483 
## 
## Means with the same letter are not significantly different.
## 
##    rendimiento groups
## T1       101.0      a
## T3        76.6      b
## T2        59.6      c
## T4        53.0     cd
## T5        43.4      d