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

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

\[f_{jl}^{(i)} = \begin{cases} 1 & \text{ si el tratamiento esta presente}\\ 0 & \text{ si el tratamiento no esta presente} \end{cases}\]

\[y_{jl} = \mu + f_{jl}^{(1)} \alpha_1+ f_{jl}^{(2)} \alpha_2+ f_{jl}^{(3)} \alpha_3+ f_{jl}^{(4)} \alpha_4+ \beta_j+\gamma_l+\epsilon_{jl}\]

\[y_{3,4} = \mu + f_{3,4}^{(1)} \alpha_1+ \beta_3+\gamma_4+\epsilon_{3,4}\]

library(collapsibleTree)
b<-8#columnas
k<-7#filas
library(readxl)
Actividad5 <- read_excel("C:/Users/juanc/OneDrive/Escritorio/2021-1/Suelos/disenio/Actividad5.xlsx", 
    col_types = c("text", "text", "text", 
        "numeric"))
row1= Actividad5$`row 1`
col1= Actividad5$`col 1`
treat=Actividad5$treat
y= Actividad5$y

dat2<-data.frame(row1,col1,treat,y)
dat2
##    row1 col1 treat    y
## 1     1    1     D 16.6
## 2     2    1     H 16.9
## 3     3    1     C 17.4
## 4     4    1     B 17.4
## 5     5    1     E 15.8
## 6     6    1     A 18.2
## 7     7    1     G 15.7
## 8     8    1     F 15.8
## 9     1    2     F 15.9
## 10    2    2     E 16.4
## 11    3    2     G 15.8
## 12    4    2     A 19.0
## 13    5    2     H 17.6
## 14    6    2     B 17.8
## 15    7    2     C 18.9
## 16    8    2     D 17.1
## 17    1    3     B 17.1
## 18    2    3     C 16.8
## 19    3    3     H 19.2
## 20    4    3     D 16.6
## 21    5    3     G 15.8
## 22    6    3     F 17.8
## 23    7    3     E 18.4
## 24    8    3     A 18.3
## 25    1    4     A 17.7
## 26    2    4     G 15.9
## 27    3    4     E 16.3
## 28    4    4     F 16.0
## 29    5    4     C 17.6
## 30    6    4     D 17.8
## 31    7    4     H 18.1
## 32    8    4     B 18.3
## 33    1    5     C 17.4
## 34    2    5     B 17.0
## 35    3    5     D 16.8
## 36    4    5     H 19.2
## 37    5    5     A 20.3
## 38    6    5     E 18.4
## 39    7    5     F 15.9
## 40    8    5     G 15.7
## 41    1    6     E 16.5
## 42    2    6     F 16.0
## 43    3    6     A 16.9
## 44    4    6     G 15.9
## 45    5    6     D 17.1
## 46    6    6     C 17.5
## 47    7    6     B 17.4
## 48    8    6     H 19.6
## 49    1    7     G 15.8
## 50    2    7     A 16.9
## 51    3    7     F 15.9
## 52    4    7     E 16.5
## 53    5    7     B 17.6
## 54    6    7     H 19.4
## 55    7    7     D 17.1
## 56    8    7     C 18.3
collapsibleTreeSummary(dat2, hierarchy = c('col1', 'row1', 'treat'))
N<-xtabs(~treat + col1, data <- dat2);N
##      col1
## treat 1 2 3 4 5 6 7
##     A 1 1 1 1 1 1 1
##     B 1 1 1 1 1 1 1
##     C 1 1 1 1 1 1 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(dat2$y, by<-list(treat<-dat2$treat),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(dat2$y, by<-list(row1<-dat2$row1),FUN<-sum);
Tr<-matrix(T2$x,k,1);# totales filas
## Warning in matrix(T2$x, k, 1): la longitud de los datos [8] no es un submúltiplo
## o múltiplo del número de filas [7] en la matriz
T3<-aggregate(dat2$y, by<-list(col1<-dat2$col1),FUN<-sum);
Tc<-matrix(T3$x,b,1)#totales Columna
## Warning in matrix(T3$x, b, 1): la longitud de los datos [7] no es un submúltiplo
## o múltiplo del número de filas [8] en la matriz
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(y~treat+row1+col1, dat2);
                anova(myfit);
                  }}}
med_trat = tapply(dat2$y, dat2$treat, mean)
barplot(med_trat, ylim = c(0,50),
        main='Distribución de la respuesta por tratamiento')
abline(h=mean(dat2$y))

med_trat = tapply(dat2$y, dat2$row1, mean)
barplot(med_trat, ylim = c(0,70),
        main='Distribución de la respuesta por Bloque fila')
abline(h=mean(dat2$y))

med_trat = tapply(dat2$y, dat2$col1, mean)
barplot(med_trat, ylim = c(0,50),
        main='Distribución de la respuesta por Bloque columna')
abline(h=mean(dat2$y))

Asignación

  • 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