Diseño Completamente Aleatorizado (DCA)

Características principales de un DCA

Se investiga el efecto de un factor con a niveles (o a tratamientos) sobre una variable respuesta.
Las unidades experimentales son homogéneas.
Los factores perturbadores han sido satisfactoriamente controlados en el experimento.
Se realizó un proceso de aleatorización en la asignación de unidades experimentales.

Caso de Estudio

Se quiere ver cómo un factor, Nutriente, afecta el crecimiento de una especie de planta, donde el factor tiene 5 niveles (Sin Nutriente, N, NK, NP, NPK).

Un DCA se diseñaría de tal manera que tuvieramos varias unidades experimentales independientes en cada tratamiento.

La aleatorización se realizaría espacialmente para la selección de las plantas o la ubicación de las parcelas de muestreo, pero también se realizaría en la asignación de los tratamientos a cada una de las plantas o parcelas.

0. Configuración inicial-Librerias requeridas

#wd="C:/Users/MAO/Desktop/Aplicada III"       
#setwd(wd)                               # Establecer el directorio de trabajo  

#library("easypackages")                  # Libreria especial para hacer carga automática de librerias

#lib_req<-c("ggplot2", "car", "nortest", "lmtest", "agricolae", "lsmeans", "emmeans", "multcomp", "PMCMRplus", "nlme", "coin", "rcompanion")# Listado de librerias requeridas por el script
#easypackages::packages(lib_req)         # Verificación, instalación y carga de librerias.

1. Lectura y transformación de datos

data <- read.table("Plantas.txt", header=TRUE, dec=".", stringsAsFactors = F,encoding ="UTF-8")

names(data)

## [1] "Finca"     "Nutriente" "Altura"

str(data)

## 'data.frame':    25 obs. of  3 variables:
##  $ Finca    : int  1 2 3 4 5 1 2 3 4 5 ...
##  $ Nutriente: chr  "1Sin" "1Sin" "1Sin" "1Sin" ...
##  $ Altura   : int  7 7 9 11 15 12 12 17 18 18 ...

table(data$Nutriente)

## 
## 1Sin   2N  3NP 4NPK  5NK 
##    5    5    5    5    5

data$Nutriente = factor(data$Nutriente)
str(data)

## 'data.frame':    25 obs. of  3 variables:
##  $ Finca    : int  1 2 3 4 5 1 2 3 4 5 ...
##  $ Nutriente: Factor w/ 5 levels "1Sin","2N","3NP",..: 1 1 1 1 1 2 2 2 2 2 ...
##  $ Altura   : int  7 7 9 11 15 12 12 17 18 18 ...

2. Visualización de datos

2.1 Resumen descriptivo

coef.var=function(x){
  sd(x)/mean(x)
}

media <- tapply(data$Altura, data$Nutriente, mean)
desvest <- tapply(data$Altura, data$Nutriente, sd)
CV <- tapply(data$Altura, data$Nutriente, coef.var)

Resumen <- cbind(media, desvest, CV)
round(Resumen,2)

##      media desvest   CV
## 1Sin   9.8    3.35 0.34
## 2N    15.4    3.13 0.20
## 3NP   17.6    2.07 0.12
## 4NPK  21.6    2.61 0.12
## 5NK   10.8    2.86 0.27

2.2. Representación Gráfica

Miremos el comportamiento de la Altura de las plantas por nutriente:

stripchart(Altura~Nutriente, data, vertical=T,pch=18,xlab="Nutriente")
abline(h=mean(data$Altura),col="blue", lty=2)
points(media, pch=19, col="red")

¿Que podemos decir de la grafica?

3. Estimación de ANOVA para un DCA

Contraste de hipotesis

Con el fin de evaluar si el factor Nutriente, afecta el crecimiento de la planta en estudio, debemos comprobar el siguiente contraste de hipotesis:

H0: No hay un efecto del factor de tratamiento sobre el crecimiento de las plantas H1: Al menos uno de los tratamientos tiene un efecto significativo sobre el crecimiento de las plantas

3.1. Ajuste modelo DCA:

mod.1<-lm(Altura~Nutriente, data)

Veamos el resultado (Recuerda que antes de interpretar, es necesario validar los supuestos del modelo)

anova(mod.1)

## Analysis of Variance Table
## 
## Response: Altura
##           Df Sum Sq Mean Sq F value    Pr(>F)    
## Nutriente  4 475.76  118.94  14.757 9.128e-06 ***
## Residuals 20 161.20    8.06                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Al evaluar la significacia del estadístico de prueba F, se observa un p-valor=0.000, indicando que al menos uno uno de los tratamientos tiene un efecto significativo sobre el crecimiento de las plantas.

Recordemos que el ANOVA solo nos indica si existe o no al menos un tratamiento que presente un efecto significativo sobre el crecimiento de las plantas, pero: ¿cuales tratamientos difieren o presentan mejores resultados?

3.2. Inspección de los residuos

Estimando los riduos

resid <- residuals(mod.1)
pred <- fitted(mod.1)

Media Cero:

plot(pred,resid)
abline(h=mean(resid))

t.test(resid, mu = 0, alternative = c("two.sided"))

## 
##  One Sample t-test
## 
## data:  resid
## t = -1.6704e-16, df = 24, p-value = 1
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -1.069782  1.069782
## sample estimates:
##     mean of x 
## -8.658005e-17

Normalidad

qqPlot(resid, pch=20)

## [1]  5 25

#Cuando n<30
shapiro.test(resid)

## 
##  Shapiro-Wilk normality test
## 
## data:  resid
## W = 0.94387, p-value = 0.1818

#Cuando n>=30
#nortest::ad.test(resid)

Homogeneidad de Varianzas:

#Bajo cumplimiento de normalidad.
bartlett.test(resid~Nutriente, data)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  resid by Nutriente
## Bartlett's K-squared = 0.93309, df = 4, p-value = 0.9198

#Bajo no normalidad
#car::leveneTest(resid~Nutriente, data)

Independencia:

dwtest(mod.1)

## 
##  Durbin-Watson test
## 
## data:  mod.1
## DW = 1.7819, p-value = 0.07288
## alternative hypothesis: true autocorrelation is greater than 0

4. Pruebas de Comparaciones Múltiple (Pruebas Postanova)

Aunque existen diferentes pruebas, Se eligirá una u otra dependiendo del riesgo que estemos dispuestos a correr al aceptar más o menos diferencias significativas.

Utilizando la libreria “agricolae”:

# Prueba de Tukey (consevardor) - Alpha Global
agricolae::HSD.test(mod.1,"Nutriente", alpha=0.05, console=TRUE, group=TRUE)

## 
## Study: mod.1 ~ "Nutriente"
## 
## HSD Test for Altura 
## 
## Mean Square Error:  8.06 
## 
## Nutriente,  means
## 
##      Altura      std r Min Max
## 1Sin    9.8 3.346640 5   7  15
## 2N     15.4 3.130495 5  12  18
## 3NP    17.6 2.073644 5  14  19
## 4NPK   21.6 2.607681 5  19  25
## 5NK    10.8 2.863564 5   7  15
## 
## Alpha: 0.05 ; DF Error: 20 
## Critical Value of Studentized Range: 4.231857 
## 
## Minimun Significant Difference: 5.372958 
## 
## Treatments with the same letter are not significantly different.
## 
##      Altura groups
## 4NPK   21.6      a
## 3NP    17.6     ab
## 2N     15.4     bc
## 5NK    10.8     cd
## 1Sin    9.8      d

# Prueba de Fisher (no consevardor) -Alpha individual para cada comparación.
#agricolae::LSD.test(mod.1,"Nutriente", alpha=0.05, console=TRUE, group=TRUE)

# Prueba de Newman-Kewls (Intermedio) -Alpha Intermedio entre Fisher y Tukey
#agricolae::SNK.test(mod.1,"Nutriente", alpha=0.05, console=TRUE, group=TRUE)

Utilizando la libreria “lsmeans” y “multcomp”:

leastsquare1 <- lsmeans(mod.1,~Nutriente)
pairs(leastsquare1, adjust="Tukey")

##  contrast    estimate  SE df t.ratio p.value
##  1Sin - 2N       -5.6 1.8 20  -3.119  0.0385
##  1Sin - 3NP      -7.8 1.8 20  -4.344  0.0026
##  1Sin - 4NPK    -11.8 1.8 20  -6.572  <.0001
##  1Sin - 5NK      -1.0 1.8 20  -0.557  0.9798
##  2N - 3NP        -2.2 1.8 20  -1.225  0.7372
##  2N - 4NPK       -6.2 1.8 20  -3.453  0.0189
##  2N - 5NK         4.6 1.8 20   2.562  0.1163
##  3NP - 4NPK      -4.0 1.8 20  -2.228  0.2101
##  3NP - 5NK        6.8 1.8 20   3.787  0.0091
##  4NPK - 5NK      10.8 1.8 20   6.015  0.0001
## 
## P value adjustment: tukey method for comparing a family of 5 estimates

multcomp::cld(leastsquare1, alpha=.05, Letters=letters)

##  Nutriente lsmean   SE df lower.CL upper.CL .group
##  1Sin         9.8 1.27 20     7.15     12.4  a    
##  5NK         10.8 1.27 20     8.15     13.4  ab   
##  2N          15.4 1.27 20    12.75     18.0   bc  
##  3NP         17.6 1.27 20    14.95     20.2    cd 
##  4NPK        21.6 1.27 20    18.95     24.2     d 
## 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 5 estimates 
## significance level used: alpha = 0.05 
## NOTE: Compact letter displays can be misleading
##       because they show NON-findings rather than findings.
##       Consider using 'pairs()', 'pwpp()', or 'pwpm()' instead.

5. No Normalidad en los residuales

¿Y cuando no se cumple el supuesto de Normalidad?

Hacer uso de la prueba No Paramétrica de Kruskal-Wallis, intuitivamente es igual que un ANOVA con los datos reemplazados por categorías, pero en este caso bajo el incumplimiento del supuesto de Normalidad. Sin embargo, si asume, que los datos vienen de la misma distribución (asume homogeneidad de varianzas).

Prueba No Paramétrica de Kruskal-Wallis

kruskal.test(Altura~Nutriente, data)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  Altura by Nutriente
## Kruskal-Wallis chi-squared = 19.064, df = 4, p-value = 0.0007636

Podemos observar que los nutrientes presentan un efecto significativo.

¿Que pruebas de comparación utilizar en este caso?

Para estos casos, hacemos uso de la función: posthoc.kruskal.nemenyi.test(), sin embargo, se debe tener en cuenta que esta prueba se ve afectada por los empates. De acuerdo a esto, se debe hacer un ajuste en la varianza de los datos, cambiando el metodo para determinar los p-valor.

Casos en los que no hay empates en los datos: se hace uso de dist=“Tukey”. Casos en que los se presenten empates en los datos: se hace uso de dist=“Chisq”

Prueba de comparación multiple

#Cuando No hay empates.
#PMCMRplus::kwAllPairsNemenyiTest(x=data$Altura, g=data$Nutriente, dist="Tukey")

#Si hay empates (algunos valores de la variables respuesta son iguales).
PMCMRplus::kwAllPairsNemenyiTest(x=data$Altura, g=data$Nutriente, dist="Chisquare")

## Warning in kwAllPairsNemenyiTest.default(x = data$Altura, g = data$Nutriente, :
## Ties are present. Chi-sq was corrected for ties.

## 
##  Pairwise comparisons using Nemenyi's all-pairs test with chi-square approximation

## data: data$Altura and data$Nutriente

##      1Sin   2N     3NP    4NPK  
## 2N   0.5973 -      -      -     
## 3NP  0.1864 0.9545 -      -     
## 4NPK 0.0085 0.3893 0.8329 -     
## 5NK  0.9994 0.7408 0.2921 0.0189

## 
## P value adjustment method: single-step

## alternative hypothesis: two.sided

6. No homogeneidad de varianzas en los residuales

Prueba alternativa Paramétrica cuando no hay homogeneidad de varianzas

Para este caso utilizamos la función gls(). Los minimos cuadrados generalizados (en ingles, generalized least squares (GLS)) es una técnica para la estimación de los parametros desconocidos en un modelo de regresión lineal.

El GLS es adecuado cuando las varianzas de las observaciones son desiguales, o cuando existe un cierto grado de correlación entre las observaciones.

6.1. Estimación mediante Mínimos Cuadrados Generalizados (asume normalidad).

Incorporamos en el modelo la matriz de varianzas desiguales, donde cada varianza depende de su tratamiento:

mod_heterog<-nlme::gls(Altura~1+Nutriente,
                       weight=varComb(varIdent(form=~1|Nutriente)),
                       data)

summary(mod_heterog)

## Generalized least squares fit by REML
##   Model: Altura ~ 1 + Nutriente 
##   Data: data 
##        AIC      BIC   logLik
##   125.5166 135.4739 -52.7583
## 
## Combination of variance functions: 
##  Structure: Different standard deviations per stratum
##  Formula: ~1 | Nutriente 
##  Parameter estimates:
##      1Sin        2N       3NP      4NPK       5NK 
## 1.0000000 0.9354152 0.6196199 0.7791983 0.8556557 
## 
## Coefficients:
##               Value Std.Error  t-value p-value
## (Intercept)     9.8  1.496660 6.547913  0.0000
## Nutriente2N     5.6  2.049387 2.732524  0.0128
## Nutriente3NP    7.8  1.760678 4.430111  0.0003
## Nutriente4NPK  11.8  1.897367 6.219144  0.0000
## Nutriente5NK    1.0  1.969770 0.507674  0.6172
## 
##  Correlation: 
##               (Intr) Ntrn2N Ntr3NP Nt4NPK
## Nutriente2N   -0.730                     
## Nutriente3NP  -0.850  0.621              
## Nutriente4NPK -0.789  0.576  0.671       
## Nutriente5NK  -0.760  0.555  0.646  0.599
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -1.7360769 -0.8366617  0.1533924  0.6751410  1.5538002 
## 
## Residual standard error: 3.346634 
## Degrees of freedom: 25 total; 20 residual

El weight nos permite determinar que tanto crecen o decrecen las variaznas dentro de cada tratamiento respecto a la primera.

¿Como determinar si los tratamiento tienen un efecto sobre la variable respuesta?

Estimamos el ANOVA para el modelo construido:

anova(mod_heterog)

## Denom. DF: 20 
##             numDF  F-value p-value
## (Intercept)     1 878.3935  <.0001
## Nutriente       4  14.9618  <.0001

Podemos observar que los nutrientes presentan un efecto significativo.

6.2. Prueba de comparación

leastsquare1 <- lsmeans(mod_heterog,~Nutriente)
pairs(leastsquare1, adjust="tukey")

##  contrast    estimate   SE   df t.ratio p.value
##  1Sin - 2N       -5.6 2.05 7.99  -2.733  0.1341
##  1Sin - 3NP      -7.8 1.76 6.86  -4.430  0.0187
##  1Sin - 4NPK    -11.8 1.90 7.72  -6.219  0.0019
##  1Sin - 5NK      -1.0 1.97 7.93  -0.508  0.9842
##  2N - 3NP        -2.2 1.68 7.07  -1.310  0.6945
##  2N - 4NPK       -6.2 1.82 7.85  -3.403  0.0548
##  2N - 5NK         4.6 1.90 8.08   2.424  0.2015
##  3NP - 4NPK      -4.0 1.49 7.62  -2.685  0.1469
##  3NP - 5NK        6.8 1.58 7.30   4.301  0.0193
##  4NPK - 5NK      10.8 1.73 7.95   6.235  0.0017
## 
## Degrees-of-freedom method: satterthwaite 
## P value adjustment: tukey method for comparing a family of 5 estimates

multcomp::cld(leastsquare1, alpha=0.05, Letters=letters)

##  Nutriente lsmean    SE   df lower.CL upper.CL .group
##  1Sin         9.8 1.497 4.13     5.69     13.9  a    
##  5NK         10.8 1.281 4.01     7.25     14.4  a    
##  2N          15.4 1.400 4.08    11.54     19.3  ab   
##  3NP         17.6 0.927 4.00    15.03     20.2   b   
##  4NPK        21.6 1.166 4.00    18.36     24.8   b   
## 
## Degrees-of-freedom method: satterthwaite 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 5 estimates 
## significance level used: alpha = 0.05 
## NOTE: Compact letter displays can be misleading
##       because they show NON-findings rather than findings.
##       Consider using 'pairs()', 'pwpp()', or 'pwpm()' instead.

Para este caso se observa que las diferencias entre tratamientos cambian un poco respecto al modelo con varianzas homogeneas.

Si desearamos comparar entre el modelo con Varianzas Homogeneas vs Varianzas Heterogeneas

Debemos estimar el modelo con varianzas homogeneas, pero para que sean comparables se debe utilizar el metodo de Mestimación de maxima Verosimilitud en ambos casos.

mod_hom<-gls(Altura~1+Nutriente, data, method="ML")

mod_het<-gls(Altura~1+Nutriente,
             weight=varComb(varIdent(form=~1|Nutriente)),
             data, method="ML")

anova(mod_hom, mod_het)

##         Model df      AIC      BIC    logLik   Test  L.Ratio p-value
## mod_hom     1  6 129.5412 136.8544 -58.77059                        
## mod_het     2 10 136.2582 148.4469 -58.12909 1 vs 2 1.282999  0.8643

¿Que podemos concluir de estos resultados?

7. No normalidad y no homogeneidad de varianzas en los residuales

7.1. Prueba de indpendencia: Para determinar si existen diferencias significativas entre los tratamiento.

Para determinar si existen diferencias significativas entre los tratamiento, haremos una prueba de indpendencia. Si los tratamientos son independientes, es porque cada tratamiento proviene de una distribución de probabilidad diferente. Indicando que hay diferencias entre los tratamiento.

coin::independence_test(Altura ~ Nutriente, data)

## 
##  Asymptotic General Independence Test
## 
## data:  Altura by Nutriente (1Sin, 2N, 3NP, 4NPK, 5NK)
## maxT = 3.1834, p-value = 0.007017
## alternative hypothesis: two.sided

7.2. Test de Permutaciones: Para determinar entre que tratamientos existe diferencias.

PT <- rcompanion::pairwisePermutationTest(Altura ~ Nutriente, data, method="fdr")
PT

##         Comparison   Stat  p.value p.adjust
## 1    1Sin - 2N = 0 -2.084  0.03712  0.05557
## 2   1Sin - 3NP = 0 -2.529  0.01145  0.03048
## 3  1Sin - 4NPK = 0 -2.731 0.006317  0.03048
## 4   1Sin - 5NK = 0  -0.53   0.5961  0.59610
## 5     2N - 3NP = 0 -1.261   0.2074  0.23040
## 6    2N - 4NPK = 0 -2.307  0.02105  0.04210
## 7     2N - 5NK = 0  1.952  0.05089  0.06361
## 8   3NP - 4NPK = 0 -2.065   0.0389  0.05557
## 9    3NP - 5NK = 0  2.507  0.01219  0.03048
## 10  4NPK - 5NK = 0  2.732 0.006294  0.03048

Presentando las diferecnias en letras a partir de los p.adjust:

cldList(p.adjust ~ Comparison, data = PT, threshold  = 0.05)

##   Group Letter MonoLetter
## 1  1Sin      a        a  
## 2    2N     ab        ab 
## 3   3NP     bc         bc
## 4  4NPK      c          c
## 5   5NK      a        a

Diseño Completamente Aleatorizado (DCA)

Ivan Mauricio Bermudez Vera

Estadística Aplicada III

Características principales de un DCA

Caso de Estudio

0. Configuración inicial-Librerias requeridas

1. Lectura y transformación de datos

2. Visualización de datos

2.1 Resumen descriptivo

2.2. Representación Gráfica

Miremos el comportamiento de la Altura de las plantas por nutriente:

¿Que podemos decir de la grafica?

3. Estimación de ANOVA para un DCA

Contraste de hipotesis

3.1. Ajuste modelo DCA:

3.2. Inspección de los residuos

Estimando los riduos

Media Cero:

Normalidad

Homogeneidad de Varianzas:

Independencia:

4. Pruebas de Comparaciones Múltiple (Pruebas Postanova)

Utilizando la libreria “agricolae”:

Utilizando la libreria “lsmeans” y “multcomp”:

5. No Normalidad en los residuales

¿Y cuando no se cumple el supuesto de Normalidad?

Prueba No Paramétrica de Kruskal-Wallis

¿Que pruebas de comparación utilizar en este caso?

Prueba de comparación multiple

6. No homogeneidad de varianzas en los residuales

Prueba alternativa Paramétrica cuando no hay homogeneidad de varianzas

6.1. Estimación mediante Mínimos Cuadrados Generalizados (asume normalidad).

¿Como determinar si los tratamiento tienen un efecto sobre la variable respuesta?

6.2. Prueba de comparación

Si desearamos comparar entre el modelo con Varianzas Homogeneas vs Varianzas Heterogeneas

¿Que podemos concluir de estos resultados?

7. No normalidad y no homogeneidad de varianzas en los residuales

7.1. Prueba de indpendencia: Para determinar si existen diferencias significativas entre los tratamiento.

7.2. Test de Permutaciones: Para determinar entre que tratamientos existe diferencias.