Clase Diseno 02/05
Wilder Harvey Clavijo Beltran
2023-05-02
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)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)
Analisis de Regresion Lineal Simple.
plot(x = data$mo,
y = data$rto,
pch = 16)
Modelo Objetivo
\[\hat{y} = \beta_0 + \beta_1x + \epsilon\] El epsilon se consideran todos los factores que influyen en el rendimiento, por fuera de la materia organica. Para ello se suelen emplear los minimos cuadrados para establecer los parametros de dicho modelo.
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\] En los modelos de regresion, el intercepto no suele tener importancia practica. (-3.48)
La pendiente expresa la relacion de la variables ((y) frente a la variable (x)) Por cada aumento de una unidad de Materia Organica, el Rendimiento aumenta 2.57. El aumento siempre tiene un limite en su relacion.
El P-value solo tiene importancia en relacion a la pendiente Al tener un valor menor al 5%, se puede rechazar la hipotesis nula, que dice que no existe relacion entre las variables, ya que SI existe la pendiente.
Para valorar la eficiencia de un modelo se revisa el valor de R^2 Ajustado, y se considera bueno con valores superiores al 70%. Este valor se considera bueno unicamente en la Regresion Lineal Simple
Revision de supuestos del modelo
# 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.2495Un p_value mayor que 5 se considera normal.
# Homosedasticidad
plot(mod1$residuals, pch=16)
Analisis de Covarianzas
Mezcla entre regresion y analisis de varianzas
FSCA - Factorial Simple con Arreglo al Azar
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 sfiModelo
Pasa de ser un ANOVA (Analisis de Varianzas) a un ANCOVA (Analisis de Covarianza) - En este caso la covariable es la MO, ya que varia de parcela a parcela.
- Theta es la regresion entre el Rto (rendimiento) y la MO (Materia Organica), LA PENDIENTE.
\[y_{ij} = \mu + \tau_i + \theta(x_{ij}-\bar{x}) + \epsilon_{ij} \] Epsilon es en esencia todo aquellos que se pudo medir y no se hizo.
Hipotesis 1 \[H_0: \theta = 0\] Hipotesis 2 \[H_0: \mu_{s0} = \mu_{sf} = \mu_{si} = \mu_{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 ' ' 1El P_value sigue siendo menor a 5% (0.05), por ende si existe relacion entre las variables ((Si existe la pendiente)) - De igual manera, el tratamiento tambien rechaza la hipotesis nula.
boxplot(rto ~ trt, data)
Revision 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)
# Revision de posible dato atipico
boxplot(rto ~ trt, data, pch=16)
which.min(data$rto)## [1] 1data[which.min(data$rto)]## x
## 1 1
## 2 2
## 3 3
## 4 4
## 5 5
## 6 6
## 7 7
## 8 8
## 9 9
## 10 10
## 11 1
## 12 2
## 13 3
## 14 4
## 15 5
## 16 6
## 17 7
## 18 8
## 19 9
## 20 10
## 21 1
## 22 2
## 23 3
## 24 4
## 25 5
## 26 6
## 27 7
## 28 8
## 29 9
## 30 10
## 31 1
## 32 2
## 33 3
## 34 4
## 35 5
## 36 6
## 37 7
## 38 8
## 39 9
## 40 10
## 41 1
## 42 2
## 43 3
## 44 4
## 45 5
## 46 6
## 47 7
## 48 8
## 49 9
## 50 10
## 51 1
## 52 2
## 53 3
## 54 4
## 55 5
## 56 6
## 57 7
## 58 8
## 59 9
## 60 10
## 61 1
## 62 2
## 63 3
## 64 4
## 65 5
## 66 6
## 67 7
## 68 8
## 69 9
## 70 10
## 71 1
## 72 2
## 73 3
## 74 4
## 75 5
## 76 6
## 77 7
## 78 8
## 79 9
## 80 10
## 81 1
## 82 2
## 83 3
## 84 4
## 85 5
## 86 6
## 87 7
## 88 8
## 89 9
## 90 10
## 91 1
## 92 2
## 93 3
## 94 4
## 95 5
## 96 6
## 97 7
## 98 8
## 99 9
## 100 10# 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# Medidas remediales para el dato atipico
# 1. Imputar
# 2. Eliminar
# 3. Repetir el experimento en esa parcela
# Medida de tratamiento
med_trt = tapply(data$rto,data$trt,mean)
med_trt## s0 sf si sfi
## 2.740161 2.979286 3.155851 3.427611data2 = 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- Se presenta un Cluster. Se puede juntar su dato con la media del tratamiento.
# De nuevo ANCOVA con el dato atipico 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.5891El p_value fue imputado, se concluye que con un solo valor se puede dañar todo el experimento.
# 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.04902Se puede usar un umbral del 1% para Rechazar o No Rechazar la Hipotesis Nula. - Se debe hacer la claridad frente a que umbral se emplea. - Puede aplicarse en el Analisis en si mismo, no unicamente en los supuestos.
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')