library(readxl)
datacuadrado <- read_excel("D:/DOCUMENTOS/datacuadrado.xlsx")
View(datacuadrado)
df=data.frame(datacuadrado)
b<-7
k<-8
row <- factor(rep(1:k, times <- b))
col <- factor(c(rep(1,k),rep(2,k),rep(3,k),rep(4,k),rep(5,k), rep(6,k),rep(7,k)))
treat<-factor(c("F", "E", "G", "A", "H", "B", "C", "D",
"B", "C", "H", "D", "G", "F", "E", "A",
"A", "G", "E", "F", "C", "D", "H", "B",
"C", "B", "D", "H", "A", "E", "F", "G",
"E", "F", "A", "G", "D", "C", "B", "H",
"G", "A", "F", "E", "B", "H", "D", "B",
"H", "D", "B", "C", "F", "G", "A", "E" ))
y= c(15.9, 16.4, 15.8, 19, 17.6, 17.8, 18.9, 17.1,
17.1, 16.8, 19.2, 16.6, 15.8, 17.8, 18.4, 18.3,
17.7, 15.9, 16.3, 16, 17.6, 17.8, 18.1, 18.3,
17.4, 17, 16.8, 19.2, 20.3, 18.4, 15.9, 15.7,
16.5, 16, 16.9, 15.9, 17.1, 17.5, 17.4, 19.6,
15.8, 16.9, 15.9, 16.5, 17.6, 19.4, 17.1, 18.3,
18.6, 17.4, 17.4, 19.2, 16.8, 15.7, 17.4, 18.4)
df2= data.frame(row,col,treat,y)
df2
## row col treat y
## 1 1 1 F 15.9
## 2 2 1 E 16.4
## 3 3 1 G 15.8
## 4 4 1 A 19.0
## 5 5 1 H 17.6
## 6 6 1 B 17.8
## 7 7 1 C 18.9
## 8 8 1 D 17.1
## 9 1 2 B 17.1
## 10 2 2 C 16.8
## 11 3 2 H 19.2
## 12 4 2 D 16.6
## 13 5 2 G 15.8
## 14 6 2 F 17.8
## 15 7 2 E 18.4
## 16 8 2 A 18.3
## 17 1 3 A 17.7
## 18 2 3 G 15.9
## 19 3 3 E 16.3
## 20 4 3 F 16.0
## 21 5 3 C 17.6
## 22 6 3 D 17.8
## 23 7 3 H 18.1
## 24 8 3 B 18.3
## 25 1 4 C 17.4
## 26 2 4 B 17.0
## 27 3 4 D 16.8
## 28 4 4 H 19.2
## 29 5 4 A 20.3
## 30 6 4 E 18.4
## 31 7 4 F 15.9
## 32 8 4 G 15.7
## 33 1 5 E 16.5
## 34 2 5 F 16.0
## 35 3 5 A 16.9
## 36 4 5 G 15.9
## 37 5 5 D 17.1
## 38 6 5 C 17.5
## 39 7 5 B 17.4
## 40 8 5 H 19.6
## 41 1 6 G 15.8
## 42 2 6 A 16.9
## 43 3 6 F 15.9
## 44 4 6 E 16.5
## 45 5 6 B 17.6
## 46 6 6 H 19.4
## 47 7 6 D 17.1
## 48 8 6 B 18.3
## 49 1 7 H 18.6
## 50 2 7 D 17.4
## 51 3 7 B 17.4
## 52 4 7 C 19.2
## 53 5 7 F 16.8
## 54 6 7 G 15.7
## 55 7 7 A 17.4
## 56 8 7 E 18.4
library(collapsibleTree)
collapsibleTreeSummary(df2, hierarchy = c('col', 'row', 'treat'))
N<-xtabs(~treat + col, df2 <- df2);N
## col
## treat 1 2 3 4 5 6 7
## A 1 1 1 1 1 1 1
## B 1 1 1 1 1 2 1
## C 1 1 1 1 1 0 1
## D 1 1 1 1 1 1 1
## E 1 1 1 1 1 1 1
## F 1 1 1 1 1 1 1
## G 1 1 1 1 1 1 1
## H 1 1 1 1 1 1 1
v<-nlevels(treat);
n<-b*k;
T1<-aggregate(df2$y, by<-list(treat<-df2$treat),FUN<-sum);
Tt<-matrix(T1$x,v,1);# totales de Trat
T2<-aggregate(df2$y, by<-list(row<-df2$row),FUN<-sum);
Tr<-matrix(T2$x,k,1);# totales filas
T3<-aggregate(df2$y, by<-list(col<-df2$col),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");
}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)
}}}
mod2= lm(y~treat+row+col, df2)
anova(mod2)
## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## treat 7 45.519 6.5027 11.7094 1.439e-07 ***
## row 7 8.099 1.1570 2.0833 0.07165 .
## col 6 1.899 0.3164 0.5698 0.75147
## Residuals 35 19.437 0.5553
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
med_trat = tapply(df2$y, df2$treat, mean)
barplot(med_trat, ylim = c(0,50),
main='Distribución de la respuesta por tratamiento')
abline(h=mean(df2$y))
med_trat = tapply(df2$y, df2$row, mean)
barplot(med_trat, ylim = c(0,70),
main='Distribución de la respuesta por Bloque fila')
abline(h=mean(df2$y))
med_trat = tapply(df2$y, df2$col, mean)
barplot(med_trat, ylim = c(0,50),
main='Distribución de la respuesta por Bloque columna')
abline(h=mean(df2$y))
Tratamiento H tiene el efecto positivo mas grande, y es el tratamiento que muestra mejor rendimiento
Tratamiento d tiene el efecto negativo mas grande