Tarea - D. Youden
Karen Abril Corredor
20/7/2021
TAREA: * Con la siguiente tabla realizar el anova como cuadrado de Youden * Realizar el diagrma de arbol * Las filas: Latitud * Las columnas: Longitud * Con 8 tratamientos: A-H * La respuesta: Rendimiento
Tabla de datos
De la tabla anterior elimine la ultima fila, para que quedara 7 filas y 8 columnas, quedando como lo muestra la siguente imagen
Tabla de datos modificados
library(readxl)## Warning: package 'readxl' was built under R version 4.0.4Youden <- read_excel("~/R/design/Youden.xlsx",
sheet = "Hoja2");Youden## # A tibble: 56 x 4
## rol col trat ren
## <dbl> <dbl> <chr> <dbl>
## 1 1 66 D 16.6
## 2 2 66 F 15.9
## 3 3 66 B 17.1
## 4 4 66 A 17.7
## 5 5 66 C 17.4
## 6 6 66 E 16.5
## 7 7 66 G 15.8
## 8 1 67 H 16.9
## 9 2 67 E 16.4
## 10 3 67 C 16.8
## # ... with 46 more rowsdf<- Youden
trat<-Youden$trat
long<-Youden$col
lati<-Youden$rol
rendi<-Youden$ren
b<-8 #Columnas
k<-7 #filas
df<-data.frame(lati,long,trat,rendi);df## lati long trat rendi
## 1 1 66 D 16.6
## 2 2 66 F 15.9
## 3 3 66 B 17.1
## 4 4 66 A 17.7
## 5 5 66 C 17.4
## 6 6 66 E 16.5
## 7 7 66 G 15.8
## 8 1 67 H 16.9
## 9 2 67 E 16.4
## 10 3 67 C 16.8
## 11 4 67 G 15.9
## 12 5 67 B 17.0
## 13 6 67 F 16.0
## 14 7 67 A 16.9
## 15 1 68 C 17.4
## 16 2 68 G 15.8
## 17 3 68 H 19.2
## 18 4 68 E 16.3
## 19 5 68 D 16.8
## 20 6 68 A 16.9
## 21 7 68 F 15.9
## 22 1 69 B 17.4
## 23 2 69 A 19.0
## 24 3 69 D 16.6
## 25 4 69 F 16.0
## 26 5 69 H 19.2
## 27 6 69 G 15.9
## 28 7 69 E 16.5
## 29 1 70 E 15.8
## 30 2 70 H 17.6
## 31 3 70 G 15.8
## 32 4 70 C 17.6
## 33 5 70 A 20.3
## 34 6 70 D 17.1
## 35 7 70 B 17.6
## 36 1 71 A 18.2
## 37 2 71 B 17.8
## 38 3 71 F 17.8
## 39 4 71 D 17.8
## 40 5 71 E 18.4
## 41 6 71 C 17.5
## 42 7 71 H 19.4
## 43 1 72 G 15.7
## 44 2 72 C 18.9
## 45 3 72 E 18.4
## 46 4 72 H 18.1
## 47 5 72 F 15.9
## 48 6 72 B 17.4
## 49 7 72 D 17.1
## 50 1 73 F 15.8
## 51 2 73 D 17.1
## 52 3 73 A 18.3
## 53 4 73 B 18.3
## 54 5 73 G 15.7
## 55 6 73 H 19.6
## 56 7 73 C 18.3library(collapsibleTree)## Warning: package 'collapsibleTree' was built under R version 4.0.4collapsibleTreeSummary(df, hierarchy = c('long', 'lati', 'trat','rendi'))N<-xtabs(~trat + long, data <- df);N## long
## trat 66 67 68 69 70 71 72 73
## 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$rendi, by<-list(trat<-df$trat),FUN<-sum);
Tt<-matrix(T1$x,v,1);# totales de Trat## Warning in matrix(T1$x, v, 1): la longitud de los datos excede el tamaño de la
## matrizT2<-aggregate(df$rendi, by<-list(lati<-df$lati),FUN<-sum);
Tr<-matrix(T2$x,k,1);# totales filas
T3<-aggregate(df$rendi, by<-list(long<-df$lati),FUN<-sum);
Tc<-matrix(T3$x,long,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,lat),lat,1));
Dk<-diag(k)-(ek%*%t(ek))/v;
Q<-Tt-(N%*%Tc)/lat;
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(rendi~trat+lati+long, df);
anova(myfit);
}}}Se pasaron los datos por que salia un error
rm(list=ls());
lon<-8;#Columnas
lat<-7;#Filas
rol <- factor(rep(1:lat, times <- lon));
col <-factor(c(rep(66,lat),rep(67,lat),rep(68,lat),rep(69,lat),rep(70,lat),rep(71,lat),rep(72,lat),rep(73,lat)));
trat<-factor(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"))
ren<-c(16.6,15.9,17.1,17.7,17.4,16.5,15.8,16.9,16.4,16.8,15.9,17,16,16.9,17.4,15.8,19.2,16.3,16.8,16.9,15.9,17.4,19,16.6,16,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)
data1<-data.frame(rol,col,trat,ren);data1## rol col trat ren
## 1 1 66 D 16.6
## 2 2 66 F 15.9
## 3 3 66 B 17.1
## 4 4 66 A 17.7
## 5 5 66 C 17.4
## 6 6 66 E 16.5
## 7 7 66 G 15.8
## 8 1 67 H 16.9
## 9 2 67 E 16.4
## 10 3 67 C 16.8
## 11 4 67 G 15.9
## 12 5 67 B 17.0
## 13 6 67 F 16.0
## 14 7 67 A 16.9
## 15 1 68 C 17.4
## 16 2 68 G 15.8
## 17 3 68 H 19.2
## 18 4 68 E 16.3
## 19 5 68 D 16.8
## 20 6 68 A 16.9
## 21 7 68 F 15.9
## 22 1 69 B 17.4
## 23 2 69 A 19.0
## 24 3 69 D 16.6
## 25 4 69 F 16.0
## 26 5 69 H 19.2
## 27 6 69 G 15.9
## 28 7 69 E 16.5
## 29 1 70 E 15.8
## 30 2 70 H 17.6
## 31 3 70 G 15.8
## 32 4 70 C 17.6
## 33 5 70 A 20.3
## 34 6 70 D 17.1
## 35 7 70 B 17.6
## 36 1 71 A 18.2
## 37 2 71 B 17.8
## 38 3 71 F 17.8
## 39 4 71 D 17.8
## 40 5 71 E 18.4
## 41 6 71 C 17.5
## 42 7 71 H 19.4
## 43 1 72 G 15.7
## 44 2 72 C 18.9
## 45 3 72 E 18.4
## 46 4 72 H 18.1
## 47 5 72 F 15.9
## 48 6 72 B 17.4
## 49 7 72 D 17.1
## 50 1 73 F 15.8
## 51 2 73 D 17.1
## 52 3 73 A 18.3
## 53 4 73 B 18.3
## 54 5 73 G 15.7
## 55 6 73 H 19.6
## 56 7 73 C 18.3library(collapsibleTree)
collapsibleTreeSummary(data1, hierarchy = c('col', 'rol', 'trat','ren'))N<-xtabs(~trat + col, data <- data1);N## col
## trat 66 67 68 69 70 71 72 73
## 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 1Codigo enviado por el profe para calcular el calculo del diseño cuadrado de Youden
v<-nlevels(trat);
n<-lon*lat;
T1<-aggregate(data1$ren, by<-list(trat<-data1$trat),FUN<-sum);
Tt<-matrix(T1$x,v,1);# totales de Trat
T2<-aggregate(data1$ren, by<-list(rol<-data1$rol),FUN<-sum);
Tr<-matrix(T2$x,lat,1);# totales filas
T3<-aggregate(data1$ren, by<-list(col<-data1$col),FUN<-sum);
Tc<-matrix(T3$x,lon,1)#totales Columna
i<-1:lat
j<-1:lon
if((lon==lat)&&(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 && lat<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,lat),lat,1));
Dk<-diag(lat)-(ek%*%t(ek))/v;
Q<-Tt-(N%*%Tc)/lat;
alpha_hat<-(lat/(lambda*v))*Q;
beta_hat<-(Tc/lat)-(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(ren~trat+rol+col, data1);
anova(myfit);
}}}## dado el diseño cuadrado de Youden
## con r y lambda respectivamente[1] 7
## [1] 6
## [1] "Las estimaciones son"
## [1] "mu_hat"
## [1] "0"
## [1] "alpha_hat"
##
## trat [,1]
## A 0.9895833
## B 0.2375000
## C 0.4750000
## D -0.3229167
## E -0.2895833
## F -1.0458333
## G -1.3333333
## H 1.2895833
## [1] "beta_hat"
##
## col [,1]
## 66 16.89851
## 67 16.51101
## 68 16.93393
## 69 17.29643
## 70 17.25060
## 71 17.93810
## 72 17.49851
## 73 17.54435
## [1] "gamma_hat"
## [,1]
## [1,] 16.7250
## [2,] 17.3125
## [3,] 17.5000
## [4,] 17.2125
## [5,] 17.5875
## [6,] 17.1125
## [7,] 17.1875## Analysis of Variance Table
##
## Response: ren
## Df Sum Sq Mean Sq F value Pr(>F)
## trat 7 44.137 6.3053 13.2943 3.205e-08 ***
## rol 6 3.827 0.6378 1.3448 0.26402
## col 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 ' ' 1med_trat = tapply(data1$ren, data1$trat, mean)
barplot(med_trat, ylim = c(0,20),
main='Distribución de la respuesta por tratamiento')
abline(h=mean(data$ren))
En el diagrama anterior se puede observar la respuesta de los tratamientos, se puede observar que el mejor tratamiento es H, seguido por A ya que presentan valores por encima de la media
med_trat = tapply(data1$ren, data1$rol, mean)
barplot(med_trat, ylim = c(0,20),
main='Distribución de la respuesta por Bloque fila')
abline(h=mean(data1$ren))
Con respecto a la respuesta por bloque de fila, lo valores por encima de la media solo se tiene la fila 5,tambien se puede observar que los valores de las filas fluctuan cerca de la media. el unico valor que se puede ver po9r debajo de la media es la fila 1
med_trat = tapply(data1$ren, data1$col, mean)
barplot(med_trat, ylim = c(0,20),
main='Distribución de la respuesta por Bloque columna')
abline(h=mean(data1$ren))
En el diagrama anterior se tiene la respueste por bloque de columna, se puede observar que donde tiene el valor mas alto es en la columna 71, mientras que las otras columnas tienen valores que oscilan cerca de la media y un poco por encima de esta. Tambien se pueden observar varios valores por debajo de la media como lo son la columna 66,67 y 68.