Ejercicio Cuadrados de Youden

Diseño Factorial Simple en Bloques Incompletos (Cuadrados de Youden)

• Modelo

\[y_{jl} = \mu + \sum_{i=1}^v{f_{jl}^{(i)} \alpha_i}+ \beta_j+\gamma_l+\epsilon_{jl}\]

library(readxl)
Ejercicio_Cuadrados_de_Youden <- read_excel("C:/Users/ERICK/Desktop/Semestre 2021-I/Disenio de Experimentos/Ejercicios/Ejercicio Cuadrados de Youden.xlsx")
df=Ejercicio_Cuadrados_de_Youden
Trat=Ejercicio_Cuadrados_de_Youden$Tratamiento
Long=Ejercicio_Cuadrados_de_Youden$Longitud
Latit=Ejercicio_Cuadrados_de_Youden$Latitud
Rend=Ejercicio_Cuadrados_de_Youden$Rendimiento
b<-8#columnas
k<-7#filas

library(collapsibleTree)
collapsibleTree::collapsibleTreeSummary(df, hierarchy = c('Longitud', 'Latitud', 'Tratamiento','Rendimiento' ))

N<-xtabs(~Trat + Long, data <- df);N

##     Long
## Trat 74 O 75 O 76 O 77 O 78 O 79 O 80 O 81 O
##    A    1    0    1    1    1    1    1    1
##    B    1    1    1    1    0    1    1    1
##    C    1    1    1    1    1    1    1    0
##    D    1    1    1    1    1    1    0    1
##    E    1    1    1    0    1    1    1    1
##    F    1    1    0    1    1    1    1    1
##    G    0    1    1    1    1    1    1    1
##    H    1    1    1    1    1    0    1    1

v<-nlevels(Trat);
n<-b*k;
T1<-aggregate(df$Rendimiento, by<-list(Trat<-df$Tratamiento),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
## matriz

T2<-aggregate(df$Rendimiento, by<-list(Latit<-df$Latitud),FUN<-sum);
Tr<-matrix(T2$x,k,1);# totales filas
T3<-aggregate(df$Rendimiento, by<-list(Long<-df$Longitud),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
Youden_Erick <- lm(Rend~Trat+Latit+Long, df)
anova(Youden_Erick) }}}

med_trat = tapply(df$Rendimiento, df$Tratamiento, mean)
barplot(med_trat, ylim = c(0,50),
        main='Distribución de la respuesta por tratamiento')
abline(h=mean(df$Rendimiento))

med_trat2 = tapply(df$Rendimiento, df$Latitud, mean)
barplot(med_trat2, ylim = c(0,70),
        main='Distribución de la respuesta por Bloque fila')
abline(h=mean(df$Rendimiento))