1 Caso 1

  1. En una estación experimental se realizó un experimento en el que se evaluó el efecto del tiempo de cosecha sobre el rendimiento de grano de maíz (en kg/parcela). Se diseñó un experimento con cuatro tratamientos usando una distribución de bloques completos al azar. Los tratamientos fueron 30, 35, 40 y 45 días después de ocurrida la polinización (para el tiempo de cosecha). Se presume que el número de plantas por parcela útil (X) es una variable explicativa de la producción de grano seco (Y). La variedad usada fue “V1” y el cultivo se efectuó con riego. Los datos recolectados se presentan en la siguiente tabla.

2 Datos

y = c(4.08,4.26,4.72,4,2.78,4.23,4.92,4.53,2.79,5.6,4.5,4.83,4.24,6.36,5.62,4.30); y
##  [1] 4.08 4.26 4.72 4.00 2.78 4.23 4.92 4.53 2.79 5.60 4.50 4.83 4.24 6.36 5.62
## [16] 4.30
x = c(41,40,37,32,24,36,32,38,31,44,38,40,46,48,41,40); x
##  [1] 41 40 37 32 24 36 32 38 31 44 38 40 46 48 41 40
trata = factor(rep(1:4,4)) ; trata
##  [1] 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4
## Levels: 1 2 3 4
bloq = factor(rep(1:4, each = 4)) ; bloq
##  [1] 1 1 1 1 2 2 2 2 3 3 3 3 4 4 4 4
## Levels: 1 2 3 4
data.frame(y,x,trata,bloq) -> datos; datos
##       y  x trata bloq
## 1  4.08 41     1    1
## 2  4.26 40     2    1
## 3  4.72 37     3    1
## 4  4.00 32     4    1
## 5  2.78 24     1    2
## 6  4.23 36     2    2
## 7  4.92 32     3    2
## 8  4.53 38     4    2
## 9  2.79 31     1    3
## 10 5.60 44     2    3
## 11 4.50 38     3    3
## 12 4.83 40     4    3
## 13 4.24 46     1    4
## 14 6.36 48     2    4
## 15 5.62 41     3    4
## 16 4.30 40     4    4

2.1 Modelo aditivo lineal

\(Y_{ij} = \mu + \tau_{i} + \rho_{j}+ \beta(X_{ij} - \bar{X_{..}}) + \epsilon_{ij}\) \(i=1,...,t\) ; \(j = 1,...,b\)

\(Y_{ij}\) : Rendimiento de la producción de grano en el j-ésima parcela del i-ésimo días después de la polinización

\(\mu\) : Media global del rendimiento de la producción de grano

\(\tau_{i}\) : Efecto de los i-ésimo días despues de ocurrida la polinización

\(\rho_{j}\) : Efecto de de la j-ésima parcela

\(\beta\) : Coeficiente de relación lineal entre la producción de grano con el número de plantas por parcelas útil

\(X_{ij}\) : Número de plantas por parcelas útil en la j-ésima parcela del i-ésimo día después de ocurrida la polinización

\(\bar{X_{..}}\) : Media global del número de plantas por parcela útil.

\(\epsilon_{ij}\) : Error experimental de la producción de grano en la j-ésima parcela en los i-ésimos días de después ocurrida la polinización

modelo = lm(y~x+trata + bloq)

3 Supuestos

3.1 Normalidad

\(H_0\): Los errores se distribuyen normalmente

\(H_1\): Los errores no se distrubuyen normalmente

\(\alpha\) = 0.05

res = modelo$residuals
shapiro.test(res)
## 
##  Shapiro-Wilk normality test
## 
## data:  res
## W = 0.89364, p-value = 0.06363

p-value = 0.06363 > 0.05 (No se rechaza H0)

Conclusión : A un nivel de significancia del 5%, no se puede afirmar que los errores no se distribuyen normalmente

3.2 Homogeneidad

\(H_0 : \sigma^2_{1} = \sigma^2_{2} = \sigma^2_{3} = \sigma^2_{4}\)

\(H_1\) : Al menos una varianza es diferente a los demas i = 1,2,3,4,

bartlett.test(res,trata)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  res and trata
## Bartlett's K-squared = 1.2916, df = 3, p-value = 0.7311

p-value = 0.7311 > 0.05 (No se rechaza H0)

Conclusión : A un nivel de significancia del 5%, no se puede afirmar que los errores no son homogéneos

4 Prueba de hipótesis de la influencia de la covariable

\(H0: \beta = 0\)

\(H1 : \beta \neq 0\)

Estadístico de prueba

F = \(\frac{\frac{(E_{xy})^2}{Exx}}{CMError_{ajus}}\)

Anova(modelo, type = 3) -> ancova; ancova
## Anova Table (Type III tests)
## 
## Response: y
##             Sum Sq Df F value  Pr(>F)  
## (Intercept) 0.0055  1  0.0232 0.88283  
## x           1.9666  1  8.3037 0.02046 *
## trata       3.5102  3  4.9404 0.03150 *
## bloq        0.2780  3  0.3913 0.76267  
## Residuals   1.8947  8                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Este cuadro ancova, muestra medidas ajustadas, justamente usadas para corroborar la influencia o significancia de la covariable “Número de plantas por parcela útil”

p-value = 0.02046 < 0.05 (Se rechaza H0)

Conclusión : A un nivel de significancia del 5%, existe suficiente evidencia estadística para afirmar que existe influencia de la covariable número de plantas por parcela útil, es decir la producción de grano es explicada por el número de plantas por parcela útil.

Se justifica la inclusión de la covarible al experimento

5 Prueba de medias ajustadas

\(H0 : \mu_{1,aju} = \mu_{2,aju} = ...=\mu_{4,aju}\)

\(H1\) : Al menos un \(\mu_{i,aju}\) es distinta a los demás para un i = 1,2,3,4

\(F_{calc} = \frac{CMTrat_{ajus}}{CMError_{ajus}}\)

ancova[3,4]
## [1] 0.03149942

p-value = 0.03149942 < 0.05 (Se rechaza H0)

Conclusión : A un nivel de significancia del 5%, existe evidencia estadística para afirmar que al menos unos de los i-ésimos dias después ocurrida la polinización en la cosecha brinda una producción de grano diferente a las demás, ajustado por el número de plantas por parcela útil.

6 Pruebas de comparación

6.1 DLS

6.2 Forma tradicional

Para aplicar las pruebas de comparación de las medias ajustadas, primero debemos hallar el \(\hat{\beta}\), dado que hemos ajustado

\(\hat{\beta} = \frac{E_{xy}}{E_{xx}}\)

Pero, del estadístico de la prueba de hipótesis de la influencia de la covariable.

F = \(\frac{\frac{(E_{xy})^2}{Exx}}{CMError_{ajus}}\) (1.1)

F_coef_reg = ancova[2,3]; F_coef_reg 
## [1] 8.303673
CMEajus = ancova[5,1]/ancova[5,2] ; CMEajus
## [1] 0.2368331

Para hallar \(E_{xx}\) (suma de cuadrados del error explicado por X) lo hallamos haciendo lo siguiente

modelox = lm(x~trata + bloq)
modelox %>% anova() -> anvax; anvax
## Analysis of Variance Table
## 
## Response: x
##           Df Sum Sq Mean Sq F value  Pr(>F)  
## trata      3   94.0  31.333  1.3857 0.30886  
## bloq       3  254.5  84.833  3.7518 0.05355 .
## Residuals  9  203.5  22.611                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Exx = anvax[3,2]; Exx
## [1] 203.5

Despejando de la ecuación (1.1)

\(E_{xy} = \sqrt{E_{xx}*F*CME_{ajus}}\)

Exy = sqrt(Exx*F_coef_reg*CMEajus); Exy
## [1] 20.005

Por lo tanto:

\(\hat{\beta} = \frac{20.005}{203.5}\)

beta_estim = Exy/Exx; beta_estim
## [1] 0.09830467

Ahora que hemos obtenido \(\hat{\beta}\), podemos hacer las pruebas de comparación dado que, como se ha podido corroborar, al cumplirse la significancia de la variable predictora, debemos hacer nuestras comparaciones con las medias ajustadas por el coeficiente de regresión, por ello:

\(\bar{Y}_{iaj} = \bar{Y}_{i.}- \hat{\beta}(\bar{X_{ij}} - \bar{X_{..}})\)

\(\bar{Y}_{i.}\)

ybar = tapply(y, trata, mean); ybar
##      1      2      3      4 
## 3.4725 5.1125 4.9400 4.4150

\(X_{ij}\)

xbar = tapply(x, trata, mean); xbar
##    1    2    3    4 
## 35.5 42.0 37.0 37.5

\(\bar{X_{..}}\)

xbar_total = mean(x); xbar_total
## [1] 38

\(\bar{Y}_{iaj}\)

ybar_aju = ybar - beta_estim*(xbar - xbar_total); ybar_aju
##        1        2        3        4 
## 3.718262 4.719281 5.038305 4.464152

Usando la librería emmeans

emmeans(modelo,~trata) %>% as.data.frame() %>% dplyr::select(2) %>% as.vector()
## $emmean
## [1] 3.718262 4.719281 5.038305 4.464152

6.2.1 Prueba de hipótesis

\(H0: \mu_{i\bullet_{}aju} = \mu_{j\bullet_{}aju}\)

\(H1: \mu_{i\bullet_{}aju} \neq \mu_{j\bullet_{}aju}\) para un i,j = 1,2,3,4

Criterio de decisión

\[\left | \bar{Y}_{i\bullet aj} - \bar{Y}_{j\bullet aj} \right | > t_{(1-\alpha/2, GLE_{aj})} \sqrt{CME_{aju}\left[\frac{2}{b}+ \frac{(\bar{X}_{i\bullet} - \bar{X}_{j\bullet})^2}{E_{XX}}\right]}\]

modelo2 = emmeans(modelo, ~trata)
pairs(modelo2, adjust = "none")
##  contrast        estimate    SE df t.ratio p.value
##  trata1 - trata2   -1.001 0.409  8  -2.445  0.0402
##  trata1 - trata3   -1.320 0.348  8  -3.794  0.0053
##  trata1 - trata4   -0.746 0.351  8  -2.126  0.0662
##  trata2 - trata3   -0.319 0.384  8  -0.831  0.4303
##  trata2 - trata4    0.255 0.377  8   0.677  0.5175
##  trata3 - trata4    0.574 0.345  8   1.666  0.1342
## 
## Results are averaged over the levels of: bloq
dif <- c()
DLS <- c()
combi <- t(combn(4,2))
tcri <- qt(0.975,8) 
for (i in 1:6){
  dif[i] <- abs(ybar_aju[combi[i,1]] - ybar_aju[combi[i,2]])
  DLS[i] <- tcri*sqrt(CMEajus*(2/4 + (xbar[combi[i,1]]-xbar[combi[i,2]])**2/anvax[3,2]))
} 

dif > DLS
## [1]  TRUE  TRUE FALSE FALSE FALSE FALSE
ybar_aju
##        1        2        3        4 
## 3.718262 4.719281 5.038305 4.464152

Conclusión : A un nivel de significancia del 5%, el mejor día fue el 3,2 y 4 que representan; 40, 35 y 45 días después de ocurrida la polinización, ajustada por el número de plantas por parcela útil

6.3 Tukey

6.3.1 Prueba de hipótesis

\(H0 : \mu_{i.aj} = \mu_{j.aj}\)

\(H1 : \mu_{i.aj} \neq \mu_{j.aj}\)

pairs(modelo2, adjust = "tukey")
##  contrast        estimate    SE df t.ratio p.value
##  trata1 - trata2   -1.001 0.409  8  -2.445  0.1450
##  trata1 - trata3   -1.320 0.348  8  -3.794  0.0220
##  trata1 - trata4   -0.746 0.351  8  -2.126  0.2238
##  trata2 - trata3   -0.319 0.384  8  -0.831  0.8387
##  trata2 - trata4    0.255 0.377  8   0.677  0.9029
##  trata3 - trata4    0.574 0.345  8   1.666  0.3984
## 
## Results are averaged over the levels of: bloq 
## P value adjustment: tukey method for comparing a family of 4 estimates

6.4 Dunnett

6.4.1 Prueba de hipótesis

\(H0 : \mu_{i.aj} = \mu_{3.aj}\) para i = 1,3,4,

\(H1 : \mu_{i.aj} \neq \mu_{3.aj}\)

dunnet_emm = contrast(modelo2, method = "trt.vs.ctrl", ref = 3, adjust = "none")
summary(dunnet_emm, adjust = "dunnett")
##  contrast        estimate    SE df t.ratio p.value
##  trata1 - trata3   -1.320 0.348  8  -3.794  0.0139
##  trata2 - trata3   -0.319 0.384  8  -0.831  0.7355
##  trata4 - trata3   -0.574 0.345  8  -1.666  0.3004
## 
## Results are averaged over the levels of: bloq 
## P value adjustment: dunnettx method for 3 tests

Conclusión: A un nivel de significancia del 5%, se puede afirmar que teniendo en cuesta como día testigo “3”, el rendimiento promedio de la producción de grano difiere con el día 1.