ANALISIS DE COVARIANZA

\[ Y_{ij} = \mu + \tau_i + \theta(x_{ij} - x) + \epsilon_{ij} \]

Desarrollo

Para este ejemplo vamos a usar MO como covariable.

Creación de datos

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"))

print(data)
##      x  y      rto       mo trt
## 1    1  1 2.331122 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
## 7    7  1 2.670194 2.429978  S0
## 8    8  1 2.687383 2.431847  S0
## 9    9  1 2.680132 2.480773  S0
## 10  10  1 2.695680 2.424533  S0
## 11   1  2 2.727857 2.440139  S0
## 12   2  2 2.717370 2.403335  S0
## 13   3  2 2.713824 2.484451  S0
## 14   4  2 2.761865 2.438728  S0
## 15   5  2 2.786099 2.474339  S0
## 16   6  2 2.857325 2.427638  S0
## 17   7  2 2.826777 2.458925  S0
## 18   8  2 2.834492 2.498158  S0
## 19   9  2 2.876039 2.427405  S0
## 20  10  2 2.903514 2.503355  S0
## 21   1  3 2.840741 2.455981  S0
## 22   2  3 2.916033 2.462614  S0
## 23   3  3 2.904697 2.458905  S0
## 24   4  3 2.901885 2.432015  S0
## 25   5  3 2.854587 2.481421  S0
## 26   6  3 2.892213 2.501424  sf
## 27   7  3 2.905947 2.529207  sf
## 28   8  3 2.916029 2.443091  sf
## 29   9  3 2.936128 2.541469  sf
## 30  10  3 2.970703 2.472891  sf
## 31   1  4 2.947694 2.542307  sf
## 32   2  4 2.931644 2.518541  sf
## 33   3  4 2.954246 2.494774  sf
## 34   4  4 2.908068 2.533270  sf
## 35   5  4 2.934297 2.456325  sf
## 36   6  4 2.951194 2.528067  sf
## 37   7  4 2.934343 2.490756  sf
## 38   8  4 3.012483 2.541239  sf
## 39   9  4 2.946235 2.512120  sf
## 40  10  4 3.012611 2.489996  sf
## 41   1  5 2.988537 2.542851  sf
## 42   2  5 3.000839 2.472799  sf
## 43   3  5 2.954795 2.548246  sf
## 44   4  5 3.021638 2.495217  sf
## 45   5  5 3.006176 2.547271  sf
## 46   6  5 3.053882 2.569901  sf
## 47   7  5 3.027033 2.484519  sf
## 48   8  5 3.088372 2.563898  sf
## 49   9  5 3.098469 2.560506  sf
## 50  10  5 3.088572 2.520225  sf
## 51   1  6 3.046874 2.514855  si
## 52   2  6 3.055384 2.508549  si
## 53   3  6 3.096461 2.546841  si
## 54   4  6 3.065538 2.579774  si
## 55   5  6 3.099119 2.547569  si
## 56   6  6 3.132920 2.555289  si
## 57   7  6 3.081589 2.570352  si
## 58   8  6 3.112059 2.558964  si
## 59   9  6 3.123153 2.557514  si
## 60  10  6 3.177869 2.510090  si
## 61   1  7 3.192105 2.514229  si
## 62   2  7 3.196142 2.604860  si
## 63   3  7 3.181311 2.598373  si
## 64   4  7 3.210601 2.541562  si
## 65   5  7 3.171876 2.589395  si
## 66   6  7 3.185591 2.595199  si
## 67   7  7 3.164188 2.570573  si
## 68   8  7 3.169195 2.613661  si
## 69   9  7 3.140277 2.587923  si
## 70  10  7 3.199636 2.625410  si
## 71   1  8 3.251623 2.597230  si
## 72   2  8 3.166805 2.634257  si
## 73   3  8 3.182590 2.607292  si
## 74   4  8 3.209734 2.594033  si
## 75   5  8 3.283625 2.649278  si
## 76   6  8 3.283925 2.627203 sfi
## 77   7  8 3.331177 2.601401 sfi
## 78   8  8 3.293122 2.563635 sfi
## 79   9  8 3.258775 2.620330 sfi
## 80  10  8 3.328322 2.619765 sfi
## 81   1  9 3.344397 2.663273 sfi
## 82   2  9 3.290391 2.658189 sfi
## 83   3  9 3.337710 2.664184 sfi
## 84   4  9 3.324220 2.625255 sfi
## 85   5  9 3.313467 2.685092 sfi
## 86   6  9 3.368641 2.620405 sfi
## 87   7  9 3.351135 2.639166 sfi
## 88   8  9 3.384977 2.677663 sfi
## 89   9  9 3.372687 2.720986 sfi
## 90  10  9 3.443173 2.684592 sfi
## 91   1 10 3.437970 2.683900 sfi
## 92   2 10 3.420653 2.655684 sfi
## 93   3 10 3.462132 2.761995 sfi
## 94   4 10 3.547827 2.679593 sfi
## 95   5 10 3.543037 2.699073 sfi
## 96   6 10 3.596180 2.783606 sfi
## 97   7 10 3.634288 2.750558 sfi
## 98   8 10 3.625385 2.746595 sfi
## 99   9 10 3.660591 2.809314 sfi
## 100 10 10 3.736083 2.827301 sfi

Hipotesis

Hipotesis de MO:

\[H_0: \theta = 0 \] Hipotesis de trt:

\[H_0: \mu_{s0} = \mu_{sf} = \mu_{f1} = \mu_{sf1} \]

ANCOVA

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

  • MO, vamos al pvalor 2x10-16 < 0.05 rechazo H0 Por tanto hay una relación entre la MO y el rendimiento

  • trt, vamos al pvalor 8.45x10 -8 < 0.05 rechazo H0, el efecto por tratamiento no es igual.

Visualización

boxplot(rto ~ trt, data)

  • Se observa que el rendimiento es mayor cuando se aplica el fungicida e insecticida a la vez. Es menos si no se usa ninguno de los dos.

Revision de supuestos

Normalidad de residuales.

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

  • Los residuales no son de distribución normal porque pvalor < 0.05.

  • Vamos a observar los datos en el histograma.

hist(mod2$residuals)

Datos atípicos.

boxplot(mod2$residuals)

  • Se observa que si hay un dato atipico, queremos observar en donde.
boxplot(rto ~ trt, data, pch=16)

  • En el SO (sin insecticida ni fungicida)
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 para remediar dato 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
#Calculamos la media de los datos para reemplazar este atipico. 


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
Solucionado
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
  • Se obtiene la misma conclusión pero el pavlor de trt cambió.
shapiro.test(mod3$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  mod3$residuals
## W = 0.98876, p-value = 0.5729
boxplot(rto~ trt, data2)

  • Ya no se observa dato atipico en el grafico.

  • pvalor de shapiro.test > 0.05, por tanto ya tienen distribución normal

Cuando se daña el analisis

  • Se decide que los resultados son cuestionables porque el modelo de una sola pendiente no se cumple.

Como interpretar analisis de regresión

  • Intercepto: estimado -3.47622, cuando MO= 0 el intercepto no tiene interpretación.

  • Pendiente: 2.57 significa que si aumento en 1 la unidad de MO, mejora el rendimiento 2.57. Cuanto mejora o empeora

  • p-value de la pendiente = 2.2e-16 < 0.05, rechazo la hipotesis nula, es decir que la pendiente sea 0. Por lo tanto, es diferente de 0.

  • R-squared (R^2-ajustado): 87.62, los puntos alrededor de una linea cayeron bien, existe un buen ajuste lineal.