youden
MANUEL VARGAS
17/8/2021
se realizo un ensayo observar la efectividad de 8 insecticidas en la reduccion del da;o foliar en determinado cultivo para lo cual se tomo como variable respuesta el porcentaje de da;o foliar observado 8 despues de la aplicacion de los tratamientos se considero la ubicacion espacial de cada unidad de muestreo como un factor que podia afectar la respuesta ; para lo cual se establecieron latitud y longitud como dos factores de bloqueo ##datos iniciales
resp= c(16.6,15.9,17.1,17.7,17.4,16.5,15.8,16.9,16.4,16.8,15.9,17.0,16.0,16.9,17.4,15.8,19.2,16.3,16.8,16.9,15.9,17.4,19.0,16.6,16.0,19.2,15.9,16.5,15.8,17.6,15.8,17.6,20.3,17.1,17.6,18.2,17.8,17.8,17.8,18.4,17.5,19.4,15.7,18.9,18.4,18.1,15.9,17.4,17.1,15.8,17.1,18.3,18.3,15.7,19.6,18.3)
trat= c('d','f','b','a','c','e','g','h','e','c','g','b','f','a','c','g','h','e','d','a','f','b','a','d','f','h','g','e','e','h','g','c','a','d','b','a','b','f','d','e','c','h','g','c','e','h','f','b','d','f','d','a','b','g','h','c')
lat= gl(7, 1, 56, c(1,2,3,4,5,6,7))
long= gl(8,7, 56, c(1,2,3,4,5,6,7,8))
DF= data.frame(resp,trat,lat,long)
DF## resp trat lat long
## 1 16.6 d 1 1
## 2 15.9 f 2 1
## 3 17.1 b 3 1
## 4 17.7 a 4 1
## 5 17.4 c 5 1
## 6 16.5 e 6 1
## 7 15.8 g 7 1
## 8 16.9 h 1 2
## 9 16.4 e 2 2
## 10 16.8 c 3 2
## 11 15.9 g 4 2
## 12 17.0 b 5 2
## 13 16.0 f 6 2
## 14 16.9 a 7 2
## 15 17.4 c 1 3
## 16 15.8 g 2 3
## 17 19.2 h 3 3
## 18 16.3 e 4 3
## 19 16.8 d 5 3
## 20 16.9 a 6 3
## 21 15.9 f 7 3
## 22 17.4 b 1 4
## 23 19.0 a 2 4
## 24 16.6 d 3 4
## 25 16.0 f 4 4
## 26 19.2 h 5 4
## 27 15.9 g 6 4
## 28 16.5 e 7 4
## 29 15.8 e 1 5
## 30 17.6 h 2 5
## 31 15.8 g 3 5
## 32 17.6 c 4 5
## 33 20.3 a 5 5
## 34 17.1 d 6 5
## 35 17.6 b 7 5
## 36 18.2 a 1 6
## 37 17.8 b 2 6
## 38 17.8 f 3 6
## 39 17.8 d 4 6
## 40 18.4 e 5 6
## 41 17.5 c 6 6
## 42 19.4 h 7 6
## 43 15.7 g 1 7
## 44 18.9 c 2 7
## 45 18.4 e 3 7
## 46 18.1 h 4 7
## 47 15.9 f 5 7
## 48 17.4 b 6 7
## 49 17.1 d 7 7
## 50 15.8 f 1 8
## 51 17.1 d 2 8
## 52 18.3 a 3 8
## 53 18.3 b 4 8
## 54 15.7 g 5 8
## 55 19.6 h 6 8
## 56 18.3 c 7 8b<-8#columnas
k<-7#filas##se establece entonces un arbol de jerarquia de los datos
library(collapsibleTree)
collapsibleTreeSummary(DF, hierarchy = c('long','lat','trat','resp'))##analisis descriptivo de los datos
medias = tapply(DF$resp, DF$trat, mean); medias## a b c d e f g h
## 18.18571 17.51429 17.70000 17.01429 16.90000 16.18571 15.80000 18.57143boxplot(DF$resp~DF$trat,main='varianza y tendencia central por tratamiento')
points(1:8, medias, pch = 15, col='blue')
en este grafico se puede evidenciar diferencias en la varianza de los tratamientos y en las tendencias centrales de los mismos
library(ggplot2)
ggplot(DF) + aes(trat,resp)+ geom_point()+ facet_grid(long~trat)
ggplot(DF) + aes(trat,resp)+ geom_point()+ facet_grid(lat~trat)
se puede evideciar en estos graficos el comportamiento de los tratamientos en cada factor de bloqueo, se evidencia tambien el caracter que define el cuadrado de youden ya que para cada nivel de logitud hace falta un tratamiento a evaluar
boxpbloq=boxplot(DF$resp~DF$long, main='Distribución de la varianza por niveles de long')
boxpbloq=boxplot(DF$resp~DF$lat, main='Distribución de la varianza por niveles de lat')
en estos dos graficos se puede evidenciar la varianza de la respuesta en los niveles de cada factor de bloqueo de los cuales se puede notar que la varianza en los niveles de latitud es homogenea lo cual puede indicar una baja eficiencia de bloqueo
med_trat = tapply(DF$resp, DF$trat, mean)
barplot(med_trat, ylim = c(0,20),
main='Distribución de la respuesta por tratamiento')
abline(h=mean(DF$resp))
med_lat = tapply(DF$resp, DF$lat, mean)
barplot(med_lat, ylim = c(0,20),
main='Distribución de la respuesta por bloque lat')
abline(h=mean(DF$resp))
med_long = tapply(DF$resp, DF$long, mean)
barplot(med_long, ylim = c(0,20),
main='Distribución de la respuesta por bloque long')
abline(h=mean(DF$resp))
nuevamente se evidencia homogeneidad en los niveles de latitud como factor de bloqueo se puede ver que los tratamientos f y g presentan la tendencia central mas baja que la media total lo cual los coloca como mejores tratamientos teniendo en cuenta que el objetivo es minimizar el da;o foliar
#MODELO ESTABLECIDO \[y_{jl} = \mu + \sum_{i=1}^v{f_{jl}^{(i)} \alpha_i}+ \beta_j+\gamma_l+\epsilon_{jl}\] \[y_{jl} = \mu + f_{jl}^{(1)} \alpha_1+ f_{jl}^{(2)} \alpha_2+ f_{jl}^{(3)} \alpha_3+ f_{jl}^{(4)} \alpha_4+ f_{jl}^{(5)} \alpha_5+ f_{jl}^{(6)} \alpha_6+ f_{jl}^{(7)} \alpha_7+ f_{jl}^{(8)} \alpha_8+ \beta_j+\gamma_l+\epsilon_{jl}\] #CORREMOS ENTONCES EL ANALISIS DE VARIANZA PERTINENTE
N<-xtabs(~trat + long, data <- DF);N## long
## trat 1 2 3 4 5 6 7 8
## a 1 1 1 1 1 1 0 1
## b 1 1 0 1 1 1 1 1
## c 1 1 1 0 1 1 1 1
## d 1 0 1 1 1 1 1 1
## e 1 1 1 1 1 1 1 0
## f 1 1 1 1 0 1 1 1
## g 1 1 1 1 1 0 1 1
## h 0 1 1 1 1 1 1 1v<-nlevels(trat);
n<-b*k;
T1<-aggregate(DF$resp, by<-list(trat<-DF$trat),FUN<-sum);
Tt<-matrix(T1$x,v,1);# totales de Trat## Warning in matrix(T1$x, v, 1): data length exceeds size of matrixT2<-aggregate(DF$resp, by<-list(lat<-DF$lat),FUN<-sum);
Tr<-matrix(T2$x,k,1);# totales filas
T3<-aggregate(DF$resp, by<-list(long<-DF$long),FUN<-sum);
Tc<-matrix(T3$x,b,1)#totales Columna
i<-1:k
j<-1:b
if((b==k)&&(N[i,j]==1)){
print("dado el diseño cuadrado latino");
#si es un DCL uno corre el código antes considerado.
}else{
x<-0;z<-0;
if(v>=2 && k<v)
{
A<-N%*%t(N);
for(i in 2:v){
x1<-A[1,1]-A[i,i];
x<-x+abs(x1);
}
sum1<-sum(x);sum1
for(i in 1:v){
for(j in 1:v){
if(i!=j){
x2<-A[1,2]-A[i,j];
z<-z+abs(x2)}
}}
sum2<-sum(z);
if(sum1==0 && sum2==0){
r<-A[1,1];
lambda<-A[1,2];
cat("dado el diseño cuadrado de Youden \n con r y lambda respectivamente")
print(r);
print(lambda);
#se define el vector columna de k entradas todas en 1
ek<-c(matrix(rep(1,k),k,1));
Dk<-diag(k)-(ek%*%t(ek))/v;
Q<-Tt-(N%*%Tc)/k;
alpha_hat<-(k/(lambda*v))*Q;
beta_hat<-(Tc/k)-(t(N)%*%Q)/(lambda*v);
gamma_hat<-Tr/v;
print("Las estimaciones son");
print("mu_hat");
print("0");
print("alpha_hat");
print(alpha_hat);
print("beta_hat")
print(beta_hat)
print("gamma_hat")
print(gamma_hat);
#probando hipótesis H_alpha, H_beta y H_gamma
myfit <- lm(resp~trat+lat+long, DF);
anova(myfit);
}}}ya que presenta un error en el tama;o de la matriz se corre solo el analisis de varianza
mod1 = lm(resp~trat+lat+long, DF)
anova(mod1)## Analysis of Variance Table
##
## Response: resp
## Df Sum Sq Mean Sq F value Pr(>F)
## trat 7 44.137 6.3053 13.2943 3.205e-08 ***
## lat 6 3.827 0.6378 1.3448 0.26402
## long 7 9.542 1.3631 2.8740 0.01763 *
## Residuals 35 16.600 0.4743
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1de acuerdo con el resultado obtenido en el analisis de varianza los datos dan evidencia con un nivel de significancia menor que 0.05 en contra de la hipotesis nula entonces se procede a hacer el test de duncan y de tukey con el objetivo de comparar los tratamientos por pares
tukey y barlett test
library(agricolae)
agricolae::duncan.test(mod1,'trat', console = T)##
## Study: mod1 ~ "trat"
##
## Duncan's new multiple range test
## for resp
##
## Mean Square Error: 0.4742857
##
## trat, means
##
## resp std r Min Max
## a 18.18571 1.20336829 7 16.9 20.3
## b 17.51429 0.44131837 7 17.0 18.3
## c 17.70000 0.68799225 7 16.8 18.9
## d 17.01429 0.41403934 7 16.6 17.8
## e 16.90000 1.05198226 7 15.8 18.4
## f 16.18571 0.71514234 7 15.8 17.8
## g 15.80000 0.08164966 7 15.7 15.9
## h 18.57143 1.04037539 7 16.9 19.6
##
## Alpha: 0.05 ; DF Error: 35
##
## Critical Range
## 2 3 4 5 6 7 8
## 0.7473175 0.7855956 0.8105396 0.8284457 0.8420478 0.8527707 0.8614466
##
## Means with the same letter are not significantly different.
##
## resp groups
## h 18.57143 a
## a 18.18571 ab
## c 17.70000 bc
## b 17.51429 bc
## d 17.01429 c
## e 16.90000 cd
## f 16.18571 de
## g 15.80000 ese establece entonces el tratamiento g como el mas efectivo en la reduccion del da;o foliar causado por la plaga
MOD= aov(resp~ trat+lat+long , DF)
summary(MOD)## Df Sum Sq Mean Sq F value Pr(>F)
## trat 7 44.14 6.305 13.294 3.21e-08 ***
## lat 6 3.83 0.638 1.345 0.2640
## long 7 9.54 1.363 2.874 0.0176 *
## Residuals 35 16.60 0.474
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1TukeyHSD(MOD,'trat')## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = resp ~ trat + lat + long, data = DF)
##
## $trat
## diff lwr upr p adj
## b-a -0.6714286 -1.8570921 0.51423497 0.6088571
## c-a -0.4857143 -1.6713778 0.69994925 0.8852514
## d-a -1.1714286 -2.3570921 0.01423497 0.0547690
## e-a -1.2857143 -2.4713778 -0.10005075 0.0257245
## f-a -2.0000000 -3.1856635 -0.81433646 0.0001081
## g-a -2.3857143 -3.5713778 -1.20005075 0.0000047
## h-a 0.3857143 -0.7999493 1.57137783 0.9632174
## c-b 0.1857143 -0.9999493 1.37137783 0.9995643
## d-b -0.5000000 -1.6856635 0.68566354 0.8694297
## e-b -0.6142857 -1.7999493 0.57137783 0.7063155
## f-b -1.3285714 -2.5142350 -0.14290789 0.0191256
## g-b -1.7142857 -2.8999493 -0.52862217 0.0010634
## h-b 1.0571429 -0.1285207 2.24280640 0.1097877
## d-c -0.6857143 -1.8713778 0.49994925 0.5839231
## e-c -0.8000000 -1.9856635 0.38566354 0.3914267
## f-c -1.5142857 -2.6999493 -0.32862217 0.0049599
## g-c -1.9000000 -3.0856635 -0.71433646 0.0002424
## h-c 0.8714286 -0.3142350 2.05709211 0.2883566
## e-d -0.1142857 -1.2999493 1.07137783 0.9999833
## f-d -0.8285714 -2.0142350 0.35709211 0.3480778
## g-d -1.2142857 -2.3999493 -0.02862217 0.0415179
## h-d 1.5571429 0.3714793 2.74280640 0.0035881
## f-e -0.7142857 -1.8999493 0.47137783 0.5340963
## g-e -1.1000000 -2.2856635 0.08566354 0.0852560
## h-e 1.6714286 0.4857650 2.85709211 0.0014874
## g-f -0.3857143 -1.5713778 0.79994925 0.9632174
## h-f 2.3857143 1.2000507 3.57137783 0.0000047
## h-g 2.7714286 1.5857650 3.95709211 0.0000002se considera nuevamente el tratamiento g como el de menor media y por lo tanto el mas efectivo en la reduccion del da;o
revision de los supuestos de normalidad e igualdad de varianzas en los residuales
#normalidad de los residuales
residuals=MOD$residuals
shapiro.test(residuals)##
## Shapiro-Wilk normality test
##
## data: residuals
## W = 0.98005, p-value = 0.478plot(residuals)
hist(residuals)
se cumple el primer supuesto de normalidad lo cual se puede observar en el histograma y grafico de dispersion #igualdad de varianzas
bartlett.test(residuals~ DF$trat)##
## Bartlett test of homogeneity of variances
##
## data: residuals by DF$trat
## Bartlett's K-squared = 6.8811, df = 7, p-value = 0.4414#conclusion una ves confirmados los dos supuestos necesarios para la validez del analisis de varianza se concluye; que el insecticida ‘g’ presenta entonces la mayor efectividad en la reduccion de da;o foliar para el cultivo tambien es importante concluir que teniendo en cuenta los valores F value para los factores de bloqueo, los dos factores de bloqueo tuvieron eficiencia positiva aunque como fue observado en las graficas la eficiencia del factor de bloqueo latitud fue muy cercana a 1 porlo tanto aunque positiva, fue muy leve