Para este ejercicio se trata de un modelo lineal con mas de un factor y una varible respuesta que se van analizar. Para la resolución del ejercicio realizaremos un ANCOVA.

El experimento se trata de medir el rendimiento de un cultivo al aplicarle fungicida e inserctida o la ausencia de ellos, donde utilizaremos como covariable la materia organica (mo) y su modelo es:

\[y_{ij} = \mu + \tau_i + \theta(x_{ij}-\bar{x}) + \epsilon_{ij}\] Se van asumir los siguientes 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'))

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

Se platean 2 hipotesis, la primera es que no hay materia organica y la segundo es que la varianza entre tratamientos es nula.

#Hipotesis 1 \[H_0: \theta = 0\] #Hipotesis 2 \[H_0: \mu_{s0} = \mu_{sf} = \mu_{si} = \mu_{sfi}\] #Analisis descriptivo

boxplot(rto ~ trt, data)

Se observa que el rendimiento fue menor al no aplicarle ni insecticida ni fungicida, tuvo un rendimiento medio el cultivo si se le aplica fungicida o insecticida. Por ultimo, si le aplica el insecticida y el fungicida en simultaneo se obtuvo más rendimiento en comparación.

#ANCOVA

mod = aov(rto ~ mo + trt, data)
summary(mod)
##             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

Se obtuvo evidencia para rechazar las anteriores hipotesis nula, por ende se confirma que hay relación de la materia organica y el rendimiento de la planta; además se observa que si hay una diferencia en el rendimiento entre los tratamientos por ende el efecto no es igual entre los mismo.

#Revisión de supuestos

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

Se confrima que los datos no siguen una distribución normal (en forma de campana) al ser su p-value mayor al 5% como se envidencia acontinuación:

hist(mod$residuals)

$Datos Atipicos

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

El dato posible atipico se encuentra en el tratamiento testigo

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)
## Warning: package 'outliers' was built under R version 4.1.3
grubbs.test(mod$residuals)
## 
##  Grubbs test for one outlier
## 
## data:  mod$residuals
## G.1 = 4.36012, U = 0.80603, p-value = 0.0002265
## alternative hypothesis: lowest value -0.346216964518421 is an outlier

Hay tres manera de remediar un dato atipico: imputar, eliminar o repetir el experimento

#Imputación

#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 2.740161 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

Habiendo imputado el datos atipico realizamos el ANCOVA

mod2 =aov(rto ~ mo + trt, data2)
summary(mod2)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## mo           1  6.182   6.182 1174.12  < 2e-16 ***
## trt          3  0.283   0.094   17.93 2.65e-09 ***
## Residuals   95  0.500   0.005                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Apesar de que se llega la misma interpretación podemos ver que el p-value cambio.

Se vuelve a validar los supuestos

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

Encontramos evidencia para confirmar que hay normalidad en los datos