Análisis de covarianza
Isabel Lozano
2023-05-02
Análisis de covarianza
Para este ejercicio se quiere saber cual de los 5 tratamientos de fertilización en un cultivo dado tiene mejores efectos en el rendimiento del mismo. Cabe resaltar que las parcelas experimentales en donde se encuentran los cultivos tienen diferentes valores de materia orgánica (covariable)
# Simulación de datos
set.seed(123)
data = expand.grid(x=1:10, y=1:10)
# Rendimiento
data$rto = rnorm(100, 3, 0.3)
data$rto = sort(data$rto) + runif(100, 0, 0.1)
# Materia orgánica aplicada al suelo
data$mo = rnorm(100, 2.5, 0.1)
data$mo = sort(data$mo) + runif(100, 0, 0.1)library(ggplot2)
g1 = ggplot(data)+
aes(x,y, fill=rto)+
geom_tile()
g2 = ggplot(data)+
aes(x,y, fill=mo)+
geom_tile()
gridExtra::grid.arrange(g1, g2, nrow=1)
Análisis de Regresión Lineal Simple
plot(x = data$mo, y = data$rto, pch=16)
Modelo estadistico:
\[\hat{y} = \beta_0 + \beta_1x\] ### Supuesto de relación lineal
mod1 = lm(rto ~ mo, data)
summary(mod1)##
## Call:
## lm(formula = rto ~ mo, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.34934 -0.05963 -0.00802 0.06565 0.21694
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.47622 0.24748 -14.05 <2e-16 ***
## mo 2.56595 0.09685 26.49 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.09624 on 98 degrees of freedom
## Multiple R-squared: 0.8775, Adjusted R-squared: 0.8762
## F-statistic: 702 on 1 and 98 DF, p-value: < 2.2e-16plot(data$mo, data$rto, pch=16)
abline(mod1, col='red', lwd=3)
\[\widehat{RTO} = -3.48 + 2.57MO\] De manera visual, se puede observar que si hay una relación lineal entre la covariable, que corresponde a la materia orgánica disponible en el suelo y el rendimiento, esto debido a que cuando aumenta el valor de la covariable, de igual modo el rendimiento.
Intercepto: Normalmente no tienen importancia. Si la materia orgánica tuviera un valor= 0, no tendría sentido porque el rendimiento sería negativo; o sea no habria siquiera producción
Pendiente: Indica cuanto mejora el rendimiento cuando aumenta la materia orgánica en el suelo.
\[R^2ajustado = 87.62%\] Si este es mayor del 70% o más, es un buen modelo para predecir rendimiento.
Revisión de supuestos
# Normalidad de residuos
# # Se espera que los residuos tengan comportamiento normal
shapiro.test(mod1$residuals)##
## Shapiro-Wilk normality test
##
## data: mod1$residuals
## W = 0.98357, p-value = 0.2495# Homosedasticidad
plot(mod1$residuals, pch=16)
Por otro lado, respecto a la normalidad de los residuos, se tiene un p-valor de 0.2495, mayor al 5%, indicando así que no se rechaza la hipótesis nula y los residuos tienen un comportamiento normal
Análisis de Covarianza
FSCA - Factorial Simple Completamente al Azar
# Simulación de los datos
set.seed(123)
data = expand.grid(x=1:10, y=1:10)
# Rendimiento
data$rto = rnorm(100, 3, 0.3)
data$rto = sort(data$rto) + runif(100, 0, 0.1)
# Covariable → Materia orgánica
data$mo = rnorm(100, 2.5, 0.1)
data$mo = sort(data$mo) + runif(100, 0, 0.1)
# Tratamientos
data$trt = gl(4, 25, 100,
c('F0', 'FB', 'FM', 'FA'))
data## x y rto mo trt
## 1 1 1 2.331122 2.399377 F0
## 2 2 1 2.506251 2.404943 F0
## 3 3 1 2.554129 2.397630 F0
## 4 4 1 2.586877 2.379526 F0
## 5 5 1 2.660638 2.410530 F0
## 6 6 1 2.708506 2.392674 F0
## 7 7 1 2.670194 2.429978 F0
## 8 8 1 2.687383 2.431847 F0
## 9 9 1 2.680132 2.480773 F0
## 10 10 1 2.695680 2.424533 F0
## 11 1 2 2.727857 2.440139 F0
## 12 2 2 2.717370 2.403335 F0
## 13 3 2 2.713824 2.484451 F0
## 14 4 2 2.761865 2.438728 F0
## 15 5 2 2.786099 2.474339 F0
## 16 6 2 2.857325 2.427638 F0
## 17 7 2 2.826777 2.458925 F0
## 18 8 2 2.834492 2.498158 F0
## 19 9 2 2.876039 2.427405 F0
## 20 10 2 2.903514 2.503355 F0
## 21 1 3 2.840741 2.455981 F0
## 22 2 3 2.916033 2.462614 F0
## 23 3 3 2.904697 2.458905 F0
## 24 4 3 2.901885 2.432015 F0
## 25 5 3 2.854587 2.481421 F0
## 26 6 3 2.892213 2.501424 FB
## 27 7 3 2.905947 2.529207 FB
## 28 8 3 2.916029 2.443091 FB
## 29 9 3 2.936128 2.541469 FB
## 30 10 3 2.970703 2.472891 FB
## 31 1 4 2.947694 2.542307 FB
## 32 2 4 2.931644 2.518541 FB
## 33 3 4 2.954246 2.494774 FB
## 34 4 4 2.908068 2.533270 FB
## 35 5 4 2.934297 2.456325 FB
## 36 6 4 2.951194 2.528067 FB
## 37 7 4 2.934343 2.490756 FB
## 38 8 4 3.012483 2.541239 FB
## 39 9 4 2.946235 2.512120 FB
## 40 10 4 3.012611 2.489996 FB
## 41 1 5 2.988537 2.542851 FB
## 42 2 5 3.000839 2.472799 FB
## 43 3 5 2.954795 2.548246 FB
## 44 4 5 3.021638 2.495217 FB
## 45 5 5 3.006176 2.547271 FB
## 46 6 5 3.053882 2.569901 FB
## 47 7 5 3.027033 2.484519 FB
## 48 8 5 3.088372 2.563898 FB
## 49 9 5 3.098469 2.560506 FB
## 50 10 5 3.088572 2.520225 FB
## 51 1 6 3.046874 2.514855 FM
## 52 2 6 3.055384 2.508549 FM
## 53 3 6 3.096461 2.546841 FM
## 54 4 6 3.065538 2.579774 FM
## 55 5 6 3.099119 2.547569 FM
## 56 6 6 3.132920 2.555289 FM
## 57 7 6 3.081589 2.570352 FM
## 58 8 6 3.112059 2.558964 FM
## 59 9 6 3.123153 2.557514 FM
## 60 10 6 3.177869 2.510090 FM
## 61 1 7 3.192105 2.514229 FM
## 62 2 7 3.196142 2.604860 FM
## 63 3 7 3.181311 2.598373 FM
## 64 4 7 3.210601 2.541562 FM
## 65 5 7 3.171876 2.589395 FM
## 66 6 7 3.185591 2.595199 FM
## 67 7 7 3.164188 2.570573 FM
## 68 8 7 3.169195 2.613661 FM
## 69 9 7 3.140277 2.587923 FM
## 70 10 7 3.199636 2.625410 FM
## 71 1 8 3.251623 2.597230 FM
## 72 2 8 3.166805 2.634257 FM
## 73 3 8 3.182590 2.607292 FM
## 74 4 8 3.209734 2.594033 FM
## 75 5 8 3.283625 2.649278 FM
## 76 6 8 3.283925 2.627203 FA
## 77 7 8 3.331177 2.601401 FA
## 78 8 8 3.293122 2.563635 FA
## 79 9 8 3.258775 2.620330 FA
## 80 10 8 3.328322 2.619765 FA
## 81 1 9 3.344397 2.663273 FA
## 82 2 9 3.290391 2.658189 FA
## 83 3 9 3.337710 2.664184 FA
## 84 4 9 3.324220 2.625255 FA
## 85 5 9 3.313467 2.685092 FA
## 86 6 9 3.368641 2.620405 FA
## 87 7 9 3.351135 2.639166 FA
## 88 8 9 3.384977 2.677663 FA
## 89 9 9 3.372687 2.720986 FA
## 90 10 9 3.443173 2.684592 FA
## 91 1 10 3.437970 2.683900 FA
## 92 2 10 3.420653 2.655684 FA
## 93 3 10 3.462132 2.761995 FA
## 94 4 10 3.547827 2.679593 FA
## 95 5 10 3.543037 2.699073 FA
## 96 6 10 3.596180 2.783606 FA
## 97 7 10 3.634288 2.750558 FA
## 98 8 10 3.625385 2.746595 FA
## 99 9 10 3.660591 2.809314 FA
## 100 10 10 3.736083 2.827301 FAModelo
\[y_{ij} = \mu + \tau_i + \theta(x_{ij}-\bar{x}) + \epsilon_{ij}\] \[y_{ij}: \text{Respuesta}\] \[\mu: \text{Media general}\] \[\tau_{i}: \text{Efecto de los tratamientos}\] \[\theta(x_{ij}-\bar{x}): \text{Efecto de la covariable}\] \[\epsilon_{ij}: \text{Error experimental}\]
Hipotesis 1
\[H_0: \theta = 0\]
Hipotesis 2
\[H_0: \mu_{F0} = \mu_{FB} = \mu_{FM} = \mu_{FA}\]
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 ' ' 1Teniendo en cuenta que el p-valor para los tratamientos dió menor al 5% si se rechazaría la hipótesis nula, en ese sentido si hay diferencias significativas entre los tratamientos, aspecto que se puede visulizar mejor en el siguiente boxplot
boxplot(rto ~ trt, data)
Revisión de supuestos
# Normalidad
shapiro.test(mod2$residuals)##
## Shapiro-Wilk normality test
##
## data: mod2$residuals
## W = 0.96113, p-value = 0.004843hist(mod2$residuals)
boxplot(mod2$residuals, pch=16)
Se puede observar que no hay normalidad en los residuos, puesto que el
p-valor es menor al 5%, rechazando la normalidad. Al realizar un
boxplot, se observa un dato atípico en los residuos, el cual podría ser
el causante del rechazo de la hipótesis. Para poder hallar el dato
atípico se hace lo siguiente:
# Revisando posible dato atípico
boxplot(rto ~ trt, data, pch=16)
which.min(data$rto)## [1] 1data[which.min(data$rto), ]## x y rto mo trt
## 1 1 1 2.331122 2.399377 F0# install.packages("outliers")
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.346216964518421 is an outlier¿Cómo arreglar el dato atípico? Se pueden tomar tres caminos: * Imputar * Eliminar * Repetir experimento para esa parcela ¿Qué es la imputación? Es reemplazar la media del dato atípico por la media global. Para llevar a cabo la imputación se hace lo siguiente:
# Imputación
# Media por tratamiento
med_trt = tapply(data$rto,
data$trt,
mean)
med_trt## F0 FB FM FA
## 2.740161 2.979286 3.155851 3.427611data2 = data
data2$rto[1] = med_trt['F0']
head(data2)## x y rto mo trt
## 1 1 1 2.740161 2.399377 F0
## 2 2 1 2.506251 2.404943 F0
## 3 3 1 2.554129 2.397630 F0
## 4 4 1 2.586877 2.379526 F0
## 5 5 1 2.660638 2.410530 F0
## 6 6 1 2.708506 2.392674 F0Se vuelve a llevar a cabo el ANCOVA junto con el dato atípico imputado
# De nuevo ANCOVA en el dato atípico 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.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# Normalidad residuales para el modelo con dato imputado
shapiro.test(mod3$residuals)##
## Shapiro-Wilk normality test
##
## data: mod3$residuals
## W = 0.98906, p-value = 0.5891# Igualdad de varianzas
bartlett.test(mod3$residuals, data2$trt)##
## Bartlett test of homogeneity of variances
##
## data: mod3$residuals and data2$trt
## Bartlett's K-squared = 7.859, df = 3, p-value = 0.04902#
library(ggplot2)
ggplot(data2)+
aes(rto, mo)+
geom_point(aes(color=trt),
size=3)+
geom_smooth(aes(color=trt),
linewidth=2,
method = 'lm',
formula = 'y~x',
se=F)+
geom_smooth(method = 'lm',
formula = 'y~x',
se=F,
col='black')
Nuevamente se observa un p-valor menor al 5% para los tratamientos, lo cual rechazaría la hipótesis nula indicando que si hay diferencias significativas entre tratamientos. Por otro lado, se observa un p-valor para el supuesto de normalidad mayor al 5%, lo cual indica que los datos si son normales. Finalmente para el supuesto de homocedasticidad, no se cumple dicho supuesto, puesto que, si se tiene un umbral del 5% de significancia, si se rechazaría la hipótesis nula de que las varianzas son iguales entre tratamientos. No obstante, si se tiene en cuenta un valor umbral de significancia del 1%, no se rechazaría dicha hipótesis y por tanto, si habría homocedasticidad, aspecto que es imprescindible para que el análisis ANCOVA tenga validez.
Como conclusión se puede decir entonces que el tratamiento que mejores resultados ofreció en rendimiento corresponde con la fertilización alta.