ANALISIS DE COVARIANZA -Es la mezcla entre la regresión lineal y analisis de varianza.

set.seed(123)
data = expand.grid(x=1:10, y=1:10)
data$rto = rnorm(100, 3, 0.3)
data$rto = sort(data$rto) + runif(100, 0, 0.1)
data$mo = rnorm(100, 2.5, 0.1)
data$mo = sort(data$ mo) + runif(100, 0, 0.1)

data$trt = gl(4, 25, 100, c("S0", "sf", "si", "sfi"))

*TIP: Este es un diseño factorial simple completamente al azar, factorial simple sin bloqueo.

*En este caso vamos a tomar covariable MO, tenemos covariables, se usa ANCOVA.

MODELO ESTADISTICO \[y_{ij} = \mu + \tau_i + θ(x_{ij}−{x^-})+\epsilon_{ij}\] ¿De dónde sale este modelo?

Siendo theta la covariables que añadimos al modelo anterior Cambio de rendimiento (pendiente) El efecto que tiene la materia organica en el rendimiento. Los factores que afectan el rto son la MO, semilla y error.

Hipotesis 1 \[H_0:θ=0\]

Hipotesis 2 \[H_1:μ_{s0}=μ_{sf}=μ_{fi}=μ_{sfi}\]

mod2 = aov(rto ~ mo + trt, data)
summary(mod2)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## mo           1  6.502   6.502  989.51  < 2e-16 ***
## trt          3  0.283   0.094   14.38 8.45e-08 ***
## Residuals   95  0.624   0.007                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Interpretamos

*MO, vamos al pvalor 2x10-16 < 0.05 rechazo H01 Por tanto, la pendiente no es 0 y hay una relación entre la MO y el rendimiento

*trt, vamos al pvalor 8.45x10 -8 < 0.05 rechazo la H02

boxplot(rto ~ trt, data)

*Con este entendemos que el rendimiento es mejor con el cuarto tratamiento fungi + insecticida.

Revisión de supuestos

shapiro.test(mod2$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  mod2$residuals
## W = 0.96113, p-value = 0.004843

Normalidad: 0.0048< NO SON NORMALES, por esto vamos a usar el histograma para ver si tienen o no distribución normal

hist(mod2$residuals)

Revisión de posibles datos atípicos

boxplot(mod2$residuals)

*Vemos que si hay un dato atípico.

boxplot(rto ~ trt, data, pch=16)

which.min(data$rto)
## [1] 1
data[which.min(data$rto), ]
##   x y      rto       mo trt
## 1 1 1 2.331122 2.399377  S0
library(outliers)
grubbs.test(mod2$residuals)
## 
##  Grubbs test for one outlier
## 
## data:  mod2$residuals
## G.1 = 4.36012, U = 0.80603, p-value = 0.0002265
## alternative hypothesis: lowest value -0.346216964518422 is an outlier

Caminos remediales para el atipico

  1. Imputar: esto es sustituirlos por su media
  2. eliminar
  3. repetir experimento en la misma parcela

Vamos a imputar:

#Media por tratamiento 
med_trt = tapply(data$rto,
                data$trt,
                mean)
med_trt
##       S0       sf       si      sfi 
## 2.740161 2.979286 3.155851 3.427611
data2 = data
data2$ rto[1] = med_trt["s0"]
head(data2)
##   x y      rto       mo trt
## 1 1 1       NA 2.399377  S0
## 2 2 1 2.506251 2.404943  S0
## 3 3 1 2.554129 2.397630  S0
## 4 4 1 2.586877 2.379526  S0
## 5 5 1 2.660638 2.410530  S0
## 6 6 1 2.708506 2.392674  S0

Corriendo nuevamente el analisis de covarianza con el dato imputado

boxplot(rto~ trt, data2)

mod3 =aov(rto ~ mo + trt, data2)
summary(mod3)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## mo           1  6.068   6.068  1144.7  < 2e-16 ***
## trt          3  0.283   0.094    17.8 3.13e-09 ***
## Residuals   94  0.498   0.005                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 1 observation deleted due to missingness
##Normalida de residuales para el modelo con dato imputado
shapiro.test(mod3$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  mod3$residuals
## W = 0.98876, p-value = 0.5729

Este nos permite concluir que los residuales ya son normales pvalue>5% gracias a la imputación del dato atipico