# Nested Design : Many thanks to Zi-Xuan Wang for the R code
# Purity example  
install.packages("EMSaov")
library(EMSaov)
y <-scan("purityEX141p607.txt")
purity <- data.frame(Supplier=gl(3,12,36),batch=gl(4,3,36),y)
purity.aov <- EMSanova(y~Supplier+batch, 
                     data=purity,
                      type=c("F","R"),
                      nested=c(NA,"Supplier"),
                      level=c(3,4))
purity.aov

# If treated as CRD
purity.CRDaov <- aov(y~Supplier*batch,purity)
summary(purity.CRDaov)  # Wrong F values....

# Gun-Loading example....  
#      Many thanks to Zi-Xuan Wang for the following R code
install.packages("EMSaov")
library(EMSaov)
method<-gl(2, 18, 36, labels = c("method1","method2"))
Group<-rep(gl(3, 3, 9, labels = c("I","II","III")),2)
Team <- rep(1:3,12)
Y<-c(20.2,26.2,23.8,22,22.6,22.9,23.1,22.9,21.8,24.1,26.9,24.9,23.5,24.6,25,22.9,23.7,23.5,14.2,18,12.5,14.1,14,13.7,14.1,12.2,12.7,16.2,19.1,15.4,16.1,18.1,16,16.1,13.8,15.1)
tot.data<-data.frame(Y=Y,method=method,Group=Group,Team=Team)
anova.result<-EMSanova(Y~Group+Team+method,
                data=tot.data,
                type=c("F","R","F"),
                nested=c(NA,"Group",NA),
               level=c(3,3,2))
anova.result  


#  Thanks to Shue-Poo Chen for the following R code specifically for calculating ANOVA table for GunLoading ex
testEMS<-function(data,A,b,c){
  m<-length(A)
  n<-length(A[[1]])
  a1<-A
  for(i in 1:m){
    for(j in 1:n){
      if(A[[i]][j]==0){
        a1[[i]][j]<-b[j]
      }
      if(A[[i]][j]==2){
        a1[[i]][j]<-1
      }
      if(A[[i]][j]==1){
        a1[[i]][j]<-c[j]
      }
    }
    
  }
  
  
  
  EMS<-A
  for(i in 1:m){
    s<-c()
    for(x in 1:m){
      s[x]<-1
    }
    for(j in 1:n){
      for(r in i:m){
        if(A[[i]][j]==1||A[[i]][j]==2){
          if(A[[r]][j]!=1&&A[[r]][j]!=2){
            s[r]<-0
          }
        }
      }
    }
    for(d in 1:i-1){
      s[d]<-0
    }
    EMS[[i]]<-s
  }
  
  
  EMSX<-EMS
  
  for(i in 1:m){
    for(j in 1:m){
      if(EMS[[i]][j]==1){
        x<-c()
        for(k in 1:n){
          if(A[[i]][k]==1){
            x[k]<-1
          }else{
            x[k]<-a1[[j]][k]
          }
          
        }
        z<-1
        for(k in 1:n){
          z<-z*x[k] 
        }
        EMSX[[i]][j]<-z
        
      }
    }
  }
  
  EMSY<-EMSX
  EMSZ<-c()
  for(i in 1:(m-1)){
    for(j in 1:m){
      if(EMSX[[i]][j]!=0){
        EMSY[[i]][j]=0
        for(r in (i+1):m){
          tf<-1
          for(l in 1:m){
            if(EMSY[[i]][l]!=EMSX[[r]][l]){
              tf<-0
              break
            }
            
          }
          if(tf==1){
            EMSZ[i]<-r
            break
          }
        }
        break
      }
    }
    
  }
  
  data2<-data
  d1=length(data)
  d2=length(data[[1]])
  d3=length(data[[1]][[1]])
  d4=length(data[[1]][[1]][[1]])
  
  for(i in 1:d1){
    for(j in 1:d2){
      for(k in 1:d3){
        x<-0
        for(l in 1:d4){
          x<-x+data[[i]][[j]][[k]][l]
        }
        data2[[i]][[j]][[k]]<-x
      }
    }
  }
  data3<-data2
  for(i in 1:d1){
    for(j in 1:d2){
      x<-0
      for(k in 1:d3){
        for(l in 1:d4){
          x<-x+data[[i]][[j]][[k]][l]
        }
      }
      data3[[i]][[j]]<-x
    }
  }
  data4<-list(group1=list(method1=0,method2=0),group2=list(method1=0,method2=0),group3=list(method1=0,method2=0))
  for(i in 1:d1){
    for(l in 1:d4){
      x<-0
      for(j in 1:d2){
        x<-x+data2[[i]][[j]][[l]]
      }
      data4[[i]][[l]]<-x
    }  
  }

  
  data5<-data3
  for(i in 1:d1){
    x<-0
    for(j in 1:d2){
      for(k in 1:d3){
        for(l in 1:d4){
          x<-x+data[[i]][[j]][[k]][l]
        }
      }
    }
    data5[[i]]<-x
  }
  data6<-list(method1=0,method2=0)
  for(l in 1:d4){
    x<-0
    for(i in 1:d1){
      x<-x+data4[[i]][[l]]
    }
    data6[[l]]<-x
  }
  data7<-0
  for(i in 1:d1){
    data7<-data7+data5[[i]]
  }
  #data2是每個group中的每個team中的每個method的total
  #data3是每個group中的每個team的total
  #data4是每個group中的每個method的total
  #data5是每個group的total
  #data6是每個method的total
  #data7是整個資料的total
  
  df<-c(d4-1,d1-1,(d1-1)*(d4-1),d1*(d2-1),d1*(d2-1)*(d4-1),(d1*d2*d3*d4-1)-(d4-1)-(d1-1)-(d1-1)*(d4-1)-d1*(d2-1)-d1*(d2-1)*(d4-1),d1*d2*d3*d4-1)
  ss<-c()
  ss[1]<-0
  for(k in 1:d4){
    ss[1]<-ss[1]+data6[[k]]^2/(d1*d2*d3)
  }
  ss[1]<-ss[1]-data7^2/(d1*d2*d3*d4)
  
  ss[2]<-0
  for(i in 1:d1){
    ss[2]<-ss[2]+data5[[i]]^2/(d2*d3*d4)
  }
  ss[2]<-ss[2]-data7^2/(d1*d2*d3*d4)
  
  ss[3]<-0
  for(i in 1:d1){
    for(k in 1:d3){
      ss[3]<-ss[3]+data4[[i]][[k]]^2/(d2*d4)
    }
  }
  for(k in 1:d4){
    ss[3]<-ss[3]-data6[[k]]^2/(d1*d2*d3)
  }
  for(i in 1:d1){
    ss[3]<-ss[3]-data5[[i]]^2/(d2*d3*d4)
  }
  ss[3]<-ss[3]+data7^2/(d1*d2*d3*d4)
  
  ss[4]<-0
  for(i in 1:d1){
    for(j in 1:d2){
      ss[4]<-ss[4]+data3[[i]][[j]]^2/(d3*d4)
    }
    ss[4]<-ss[4]-data5[[i]]^2/(d2*d3*d4)
  }

  
  ss[5]<-0
  for(i in 1:d1){
    for(j in 1:d2){
      for(k in 1:d3){
        ss[5]<-ss[5]+data2[[i]][[j]][[k]]^2/d4
      }
    }
    for(j in 1:d2){
      ss[5]<-ss[5]-data3[[i]][[j]]^2/(d3*d4)
    }
    for(k in 1:d3){
      ss[5]<-ss[5]-data4[[i]][[k]]^2/(d2*d4)
    }
    ss[5]<-ss[5]+data5[[i]]^2/(d2*d3*d4)
  }
  
  ss[6]<-0
  ss[7]<-0
  for(i in 1:d1){
    for(j in 1:d2){
      for(k in 1:d3){
        for(l in 1:d4){
          ss[7]<-ss[7]+data[[i]][[j]][[k]][l]^2
        }
      }
    }
  }
  ss[7]<-ss[7]-data7^2/(d1*d2*d3*d4)
  ss[6]<-ss[7]-ss[1]-ss[2]-ss[3]-ss[4]-ss[5]
  ms<-ss/df
  f<-c(0,0,0,0,0,NA,NA)
  for(i in 1:(m-1)){
    x<-ms[EMSZ[i]]
    f[i]<-ms[i]/x
  }
  p<-c(0,0,0,0,0,NA,NA)
  for(i in 1:(m-1)){
    x<-EMSZ[i]
    p[i]<-1-pf(f[i],df[i],df[x])
  }
  x<-c("M","G","MG","T","MT","Error","Total")
  anova<-data.frame(x,df,ss,ms,f,p)

  head(anova,7)
  
}
data<-list(group1=list(team1=list(method1=c(20.2,24.1),method2=c(14.2,16.2)),team2=list(method1=c(26.2,26.9),method2=c(18,19.1)),team3=list(method1=c(23.8,24.9),method2=c(12.5,15.4))),group2=list(team1=list(method1=c(22,23.5),method2=c(14.1,16.1)),team2=list(method1=c(22.6,24.6),method2=c(14,18.1)),team3=list(method1=c(22.9,25),method2=c(13.7,16))),group3=list(team1=list(method1=c(23.1,22.9),method2=c(14.1,16.1)),team2=list(method1=c(22.9,23.7),method2=c(12.2,13.8)),team3=list(method1=c(21.8,23.5),method2=c(12.7,15.1))))


c<-c(0,0,1,1)#按順序,fixed輸入0,random輸入1
b<-c(2,3,3,2)#按順序輸入變數的數量
A<-list(M=c(1,0,0,0),G=c(0,1,0,0),MG=c(1,1,0,0),T=c(0,2,1,0),MT=c(1,2,1,0),E=c(2,2,2,1))
#目標有哪些變數，按順序填入，如果沒有就填0，有的話填1，如果有被括號括起來的話改填2
testEMS(data,A,b,c)