DCA

En la siguiente data “2-datos.xlsx” se tienen los rendimientos en Kg de mandarinas por planta en un comparativo de cinco variedades de portainjertos, con cinco repeticiones. El experimento se ha conducido bajo un Diseño Completamente al Azar.

Análisis Exploratorio

library(readxl)
library(tidyverse)
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	▃▇▇▅▆

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

Dado que tenemos un diseño balanceado la suma de los efectos de los tratamientos por grupo respecto a la media debe dar 0, para ello hacemos lo siguiente:

mu = 66.72
mu1 = 101
mu2 = 59.6
mu3 = 76.6
mu4 = 53
mu5 = 43.4
tau1 = mu1 - mu
tau2 = mu2 - mu 
tau3 = mu3 - mu 
tau4 = mu4- mu 
tau5 = mu5 - mu
tau1 + tau2 + tau3 + tau4 + tau5

## [1] 0

datos %>% ggplot(aes(x = tratamiento, y = rendimiento)) + geom_boxplot(fill = "skyblue")

Estimación de parámetros

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

La matrix de datos es:

modelo %>% model.matrix()

##    (Intercept) tratamientoT2 tratamientoT3 tratamientoT4 tratamientoT5
## 1            1             0             0             0             0
## 2            1             0             0             0             0
## 3            1             0             0             0             0
## 4            1             0             0             0             0
## 5            1             0             0             0             0
## 6            1             1             0             0             0
## 7            1             1             0             0             0
## 8            1             1             0             0             0
## 9            1             1             0             0             0
## 10           1             1             0             0             0
## 11           1             0             1             0             0
## 12           1             0             1             0             0
## 13           1             0             1             0             0
## 14           1             0             1             0             0
## 15           1             0             1             0             0
## 16           1             0             0             1             0
## 17           1             0             0             1             0
## 18           1             0             0             1             0
## 19           1             0             0             1             0
## 20           1             0             0             1             0
## 21           1             0             0             0             1
## 22           1             0             0             0             1
## 23           1             0             0             0             1
## 24           1             0             0             0             1
## 25           1             0             0             0             1
## attr(,"assign")
## [1] 0 1 1 1 1
## attr(,"contrasts")
## attr(,"contrasts")$tratamiento
## [1] "contr.treatment"

Verificación de Supuestos

Normalidad de Errores

modelo %>% residuals() %>% shapiro.test()

## 
##  Shapiro-Wilk normality test
## 
## data:  .
## W = 0.97389, p-value = 0.7441

library(moments)
modelo %>% residuals() %>% skewness()

## [1] -0.3211397

modelo %>% residuals() %>% kurtosis()

## [1] 3.15697

library(gvlma)
modelo %>% gvlma

## 
## Call:
## lm(formula = rendimiento ~ tratamiento, data = datos)
## 
## Coefficients:
##   (Intercept)  tratamientoT2  tratamientoT3  tratamientoT4  tratamientoT5  
##         101.0          -41.4          -24.4          -48.0          -57.6  
## 
## 
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance =  0.05 
## 
## Call:
##  gvlma(x = .) 
## 
##                         Value p-value                Decision
## Global Stat         3.758e+00 0.43980 Assumptions acceptable.
## Skewness            4.297e-01 0.51213 Assumptions acceptable.
## Kurtosis            2.567e-02 0.87272 Assumptions acceptable.
## Link Function      -2.186e-17 1.00000 Assumptions acceptable.
## Heteroscedasticity  3.302e+00 0.06919 Assumptions acceptable.

library(broom)
modelo %>% augment() %>% ggplot(aes(x = .resid)) +
  geom_histogram(bins = round(1 + 3.3*log(nrow(datos)))) + theme_bw()

modelo %>% plot(which = 2)

Podemos observar que las observaciones 21, 17 y 16 están ocasionando que las colas sean ligeramente pesadas, y ello se puede observar justamente en los residuos estandarizados ya que tienen números no cercanos a cero y eso estaría pasando con las observaciones vista en el gráfico

modelo %>% augment() %>% .[,8] %>% print(n = 25)

## # A tibble: 25 × 1
##    .std.resid
##         <dbl>
##  1     0     
##  2    -0.628 
##  3    -0.157 
##  4     0.157 
##  5     0.628 
##  6    -0.879 
##  7     0.848 
##  8    -0.408 
##  9     0.377 
## 10     0.0628
## 11     1.01  
## 12    -1.35  
## 13     0.0628
## 14    -0.251 
## 15     0.534 
## 16     2.20  
## 17    -2.04  
## 18    -1.10  
## 19    -0.157 
## 20     1.10  
## 21    -2.26  
## 22     0.251 
## 23     1.19  
## 24    -0.220 
## 25     1.04

Homocedasticidad de errores

modelo %>% plot(which = 1)

modelo %>% augment() %>% 
  with(lowess(x = .fitted, y = .resid)) %>% 
    as.data.frame() -> smoothed

modelo %>% augment() %>%
  ggplot(aes(x = .fitted, y = .resid)) + 
  geom_point(size = 3) + geom_hline(yintercept = 0) +
  geom_path(data = smoothed, aes(x = x, y = y), color = "red") + 
  labs(title = "Gráfico de homocedasticidad de errores", 
       x = "Valores ajustados",
       y = "Residuales") + theme_minimal()

modelo %>% plot(which = 3)

Podemos observar que la mismas observaciones tanto en normalidad como homocedasticidad, coincide. ¿Será que están relacionadas?.

Test de Barlett

Requiere previamente que se comprobe normalidad

bartlett.test(x = modelo$residuals, g = datos$tratamiento)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  modelo$residuals and datos$tratamiento
## Bartlett's K-squared = 7.2418, df = 4, p-value = 0.1237

bartlett.test(rendimiento~tratamiento, datos)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  rendimiento by tratamiento
## Bartlett's K-squared = 7.2418, df = 4, p-value = 0.1237

Test de levene

Es una prueba no paramétrica, permite elegir el estadístico de centralidad, como la mediana (para datos no normales y/o asimétricos) o la media truncada (si la distribución de los errores tienen cola grande como Cauchy).

library(car)
leveneTest(y = modelo %>% residuals(), group = datos$tratamiento, center = median)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  4  1.6389 0.2036
##       20

leveneTest(rendimiento~tratamiento, datos, center = median)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  4  1.6389 0.2036
##       20

Test de Fligner - Killen

fligner.test(x = modelo %>% residuals() , g = datos$tratamiento)

## 
##  Fligner-Killeen test of homogeneity of variances
## 
## data:  modelo %>% residuals() and datos$tratamiento
## Fligner-Killeen:med chi-squared = 5.076, df = 4, p-value = 0.2796

fligner.test(rendimiento~tratamiento, datos)

## 
##  Fligner-Killeen test of homogeneity of variances
## 
## data:  rendimiento by tratamiento
## Fligner-Killeen:med chi-squared = 5.4928, df = 4, p-value = 0.2404

En este caso cumple homocedasticidad de errores, pero que pasaría si no cumple, en esos casos podemos aplicar la transformación de boxCox que sugiere algunos valores genéricos para mejorar dicho inconveniente

library(MASS)
(modelo %>% boxcox(lambda = seq(-4,4, 0.1)) -> bc)

## $x
##   [1] -4.00000000 -3.91919192 -3.83838384 -3.75757576 -3.67676768 -3.59595960
##   [7] -3.51515152 -3.43434343 -3.35353535 -3.27272727 -3.19191919 -3.11111111
##  [13] -3.03030303 -2.94949495 -2.86868687 -2.78787879 -2.70707071 -2.62626263
##  [19] -2.54545455 -2.46464646 -2.38383838 -2.30303030 -2.22222222 -2.14141414
##  [25] -2.06060606 -1.97979798 -1.89898990 -1.81818182 -1.73737374 -1.65656566
##  [31] -1.57575758 -1.49494949 -1.41414141 -1.33333333 -1.25252525 -1.17171717
##  [37] -1.09090909 -1.01010101 -0.92929293 -0.84848485 -0.76767677 -0.68686869
##  [43] -0.60606061 -0.52525253 -0.44444444 -0.36363636 -0.28282828 -0.20202020
##  [49] -0.12121212 -0.04040404  0.04040404  0.12121212  0.20202020  0.28282828
##  [55]  0.36363636  0.44444444  0.52525253  0.60606061  0.68686869  0.76767677
##  [61]  0.84848485  0.92929293  1.01010101  1.09090909  1.17171717  1.25252525
##  [67]  1.33333333  1.41414141  1.49494949  1.57575758  1.65656566  1.73737374
##  [73]  1.81818182  1.89898990  1.97979798  2.06060606  2.14141414  2.22222222
##  [79]  2.30303030  2.38383838  2.46464646  2.54545455  2.62626263  2.70707071
##  [85]  2.78787879  2.86868687  2.94949495  3.03030303  3.11111111  3.19191919
##  [91]  3.27272727  3.35353535  3.43434343  3.51515152  3.59595960  3.67676768
##  [97]  3.75757576  3.83838384  3.91919192  4.00000000
## 
## $y
##   [1] -37.4863376 -36.3001586 -35.1193739 -33.9441315 -32.7745832 -31.6108868
##   [7] -30.4532052 -29.3017070 -28.1565663 -27.0179631 -25.8860831 -24.7611180
##  [13] -23.6432656 -22.5327298 -21.4297209 -20.3344551 -19.2471551 -18.1680500
##  [19] -17.0973749 -16.0353715 -14.9822877 -13.9383775 -12.9039012 -11.8791248
##  [25] -10.8643207  -9.8597666  -8.8657461  -7.8825478  -6.9104659  -5.9497992
##  [31]  -5.0008512  -4.0639300  -3.1393473  -2.2274188  -1.3284635  -0.4428036
##  [37]   0.4292363   1.2873291   2.1311460   2.9603563   3.7746284   4.5736296
##  [43]   5.3570268   6.1244865   6.8756755   7.6102607   8.3279097   9.0282908
##  [49]   9.7110732  10.3759272  11.0225244  11.6505377  12.2596411  12.8495104
##  [55]  13.4198227  13.9702566  14.5004924  15.0102120  15.4990993  15.9668400
##  [61]  16.4131220  16.8376359  17.2400751  17.6201367  17.9775216  18.3119359
##  [67]  18.6230912  18.9107062  19.1745075  19.4142314  19.6296250  19.8204484
##  [73]  19.9864759  20.1274985  20.2433253  20.3337866  20.3987351  20.4380486
##  [79]  20.4516316  20.4394177  20.4013716  20.3374909  20.2478072  20.1323876
##  [85]  19.9913360  19.8247934  19.6329392  19.4159903  19.1742009  18.9078622
##  [91]  18.6173014  18.3028805  17.9649944  17.6040688  17.2205580  16.8149432
##  [97]  16.3877295  15.9394422  15.4706273  14.9818423

Dado que la transformación de BoxCox lo que busca es buscar cuál de todos los modelos creados alcanza la mayor verosimilitud (que de todos los modelos cuál más se asemeja a nuestros datos), el objetivo es encontrar el lambda que maximice esa verosimilitud y en consecuencia la log(verosimilitud).

(bc$x[which.max(bc$y)] -> lambda)

## [1] 2.30303

La función transform usa la transformación de boxCox mediante esta ecuación \(y^* = y^\lambda\)

transform(datos, trendimi = rendimiento ^ lambda) -> tdatos
aov(trendimi ~ tratamiento, tdatos) %>% summary

##             Df    Sum Sq   Mean Sq F value   Pr(>F)    
## tratamiento  4 3.994e+09 998590498   97.33 8.16e-13 ***
## Residuals   20 2.052e+08  10260149                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

La función BoxCox de la librería forecast usa la transformación \(y^* = \frac{y^\lambda - 1}{\lambda}\), si bien es cierto son distintas, pero las conclusiones serán iguales (Fc y p-valor).

library(forecast)
data.frame(trendimi = BoxCox(datos$rendimiento, lambda), datos) -> t1datos
aov(trendimi ~ tratamiento,t1datos) %>% summary

##             Df    Sum Sq   Mean Sq F value   Pr(>F)    
## tratamiento  4 753092142 188273035   97.33 8.16e-13 ***
## Residuals   20  38688719   1934436                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Si bien es cierto el lambda salio aprox 2, ello significa que se tendria que aplicar una transformación cuadrada, pero en este caso corroboramos que la homogeneidad existe; entonces ¿Por qué aplicando boxcox me recomienda usar una transformación cuadrática? Ello debe deberse a las posibles observacione 16,17 y 21

Independencia de Errores

modelo %>% 
  augment %>% 
  ggplot(aes(x=1:nrow(datos),y=.resid))+
  geom_line()+
  geom_point()+
    labs(x = "Orden",
       y = "Residual", 
       title = "Evaluación de independencia",
       subtitle = "")+
  theme_minimal()

modelo %>% augment %>% dplyr::select(.resid) %>% TSA::acf()

Se puede observar que las autocorrelaciones en los desfases son independientes o igual a 0

Prueba de rachas

Prueba de rachas evalúa si los signos de los residuos tienen un patrón, sí lo tiene no son independientes, y si no tienen un patrón, significa que los residuos son independientes

dwtest() puede ser más apropiado, ya que está específicamente diseñado para detectar autocorrelación. Sin embargo, si deseas analizar la aleatoriedad de los residuos, runs.test() puede ser una opción válida

library(tseries)
runs.test(factor(modelo %>% residuals() > 0))

## 
##  Runs Test
## 
## data:  factor(modelo %>% residuals() > 0)
## Standard Normal = 1.0314, p-value = 0.3024
## alternative hypothesis: two.sided

library(lmtest)

modelo %>% dwtest(alternative = "two.sided")

## 
##  Durbin-Watson test
## 
## data:  .
## DW = 2.5528, p-value = 0.5376
## alternative hypothesis: true autocorrelation is not 0

modelo %>% durbinWatsonTest(alternative = "two.sided",
                            max.lag=10,
                            reps=1e5)

##  lag Autocorrelation D-W Statistic p-value
##    1     -0.29786898     2.5527624 0.53218
##    2     -0.17012628     2.2795580 0.73578
##    3     -0.02225730     1.9258485 0.93004
##    4      0.14258090     1.5926598 0.65692
##    5     -0.10339384     1.8642463 0.51996
##    6     -0.05378848     1.6857537 0.65322
##    7      0.04356748     1.4612865 0.92162
##    8     -0.20970797     1.9128256 0.14264
##    9      0.29096290     0.7390687 0.10814
##   10     -0.16002368     1.4475138 0.48046
##  Alternative hypothesis: rho[lag] != 0

Análisis de varianza

El análisis de varianza (ANOVA) tiene como objetivo determinar si el tratamiento aplicado genera un efecto significativo, ya sea positivo o negativo. Para ello, primero debemos asegurarnos de que se cumplen todos los supuestos del modelo aditivo lineal. Una vez verificados estos requisitos, aplicamos el ANOVA para evaluar si las diferencias observadas entre los tratamientos son estadísticamente significativas.

\(H_{0} : μ1=μ2=μ3=μ4\)

\(H1\): Al menos un \(μ_{i}\) es distinto, para \(i=1,2,3\)

También podemos hacer las siguientes pruebas de hipótesis

\(H_{1}: τ_{i} \neq 0\) (Efecto del tratamiento es significativo)

\(H_{0}: τ_{i} = 0\) (El efecto del tratamiento no es significativo)

Estadístico de prueba:

\(F_{c} = \frac{CM_{trata}}{CM_{Error}}\) ~ \(F_{(1-\alpha,k-1,N-k)}\)

k : Número de tratamientos

N : Número total de observaciones en el experimento.

De manera intuitiva se busca que “N-k” sea pequeño para que así \(SC_{Error}\) sea pequeño y por ende, el \(CM_{Error}\) también lo sea, ello hará que se rechaze nuestra hipótesis que es lo que estaríamos buscando

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

En esta oportunidad tenemos las repeticiones por tratamiento que son 5, pero que sucedería si es que queremos saber, ¿Cuál es el tamaño de prueba adecuado para que nuestra potencia de prueba sea lo más óptimo posible?, hay que tener en cuenta que si los datos presentan mayor variabilidad lo recomendable es aumentar el número de repeticiones y por ende el IC será más angosto y preciso.

Coeficiente de variabilidad

library(agricolae)
modelo %>% cv.model()

## [1] 10.66995

Tamaño de prueba

Cantidad mínima de repeticiones

\(n = \frac{2t^2\sigma_{0}^2}{d^2}\)

t = \(t_{1-\alpha/2,N-k}\) según Gutierrez (2008) y según Melo t = \(t_{1-\alpha/2,\infty}\), pero sabemos que si una “t” con infinito grados de libertad se asemaja a la normal estándar, pero también es posible considerar la potencia de prueba, dado los valores para \(\beta\) (y en consecuencia (1 - \(\beta\)) dado los valores para \(\alpha\), \(GL_{Trat}\), \(GLE\) y \(\phi\), donde

\(\phi\) = \(\sqrt{\frac{\lambda}{k}} = \sqrt{\frac{n\sum_{i=1}^{k}\tau_{i}^2}{k\sigma^2}}\)

En el caso donde no tengamos todos los efectos del tratamiento podríamos usar lo siguiente:

\(\Delta\) = \(max(\tau) - min(\tau)\), donde:

\(\phi\) = \(\sqrt{\frac{n\Delta^2}{2k\sigma^2}}\)

Ejemplo:

Por ejemplo, suponga que se desea determinar el número de repeticiones que permita alcanzar una potencia de prueba de por lo menos 0.80, con un tamaño de prueba \(α =0.01\), en un diseño experimental completamente aleatorizado con 5 tratamientos en estudio. Además, se sabe que CME = 4.3 y que \(τ_{1}\) = 3.1, \(τ_{2}\) = −0.4, \(τ_{3}\) = −0.9, \(τ_{4}\)= −2 y \(τ_{5}\) =0.2. ¿Serán suficientes 6 repeticiones?

Podemos hacerlo manualmente, y usar el gráfico de potencia de prueba:

Pero usaremos la librería pwr

t1 = 3.1
t2 = -0.4
t3 = -0.9
t4 = -2 
t5 = 0.2
sumtau = sum(t1^2+t2^2+t3^2+t4^2+t5^2)

library(pwr)
pwr.anova.test(k = 5, sig.level = 0.01, power = 0.8, f = sqrt(sumtau/(5*4.3)))

## 
##      Balanced one-way analysis of variance power calculation 
## 
##               k = 5
##               n = 6.329472
##               f = 0.8246211
##       sig.level = 0.01
##           power = 0.8
## 
## NOTE: n is number in each group

Se recomienda tener 7 repeticiones por tratamiento para que la potencia de prueba sea 0.8, por lo tanto no será suficiente 6 repetciones, ahora hagamoslo para nuestro caso de estudio con un número de repeticiones de 5.

(modelo %>% anova %>% .[2,3]-> CME)

## [1] 50.68

(sumtaus <- tau1^2 + tau2^2 + tau3^2 + tau4^2 + tau5^2)

## [1] 2055.488

pwr.anova.test(k = 5, sig.level = 0.05, n = 5, f = sqrt(sumtaus/(5*CME)))

## 
##      Balanced one-way analysis of variance power calculation 
## 
##               k = 5
##               n = 5
##               f = 2.848093
##       sig.level = 0.05
##           power = 1
## 
## NOTE: n is number in each group

Se puede observar que con 5 repeticiones ya tenemos 100% de potencia de prueba, ello indica que tenemos una muy buena cantidad de unidades experimentales para nuestro diseño.

En Resumen:

Hemos hecho un DCA, asimismo hemos sacado algunas conclusiones y posibles soluciones si fuese el caso, ¿Ahora que sabemos que al menos un efecto del tratamiento es distinto a 0, cual será ese tratamiento?, ello lo veremos cuando apliquemos Pruebas de comparación múltiple más de dos medias.