AMMI<-function(variety,envir,block,yield,PC=2,biplot=1){

  ## ADDITIVE MAIN EFFECTS MULTIPLICATIVE INTERACTION MODEL, AS APPLIED TO AGRICULTURE VARIETY TRIALS
  ## REFERENCE: H. F. Gollob, 1968. A statistical model which combine features of factor analytic and
  ## analysis of variance techniques. Psychometrika, 33, 1, 73-114.
  ## Update: 13/11/2006

  ## INPUTS:
  ## variety is a vector of variety codes (factor)
  ## envir is a vector of environment codes (factor)
  ## block is a vector of block codes (factor)
  ## yield is a vector of yields to be analysed (numerical)
  ## PC is the number of PC to be considered (default is 2)

  ## This part calculates some descriptive statistics
  overall.mean<-mean(yield)
  envir.mean<-tapply(yield,envir,mean)
  var.mean<-tapply(yield,variety,mean)  
  int.mean<-tapply(yield,list(variety,envir),mean)
  envir.num<-length(envir.mean)
  var.num<-length(var.mean)
  
  ## This part calculates the ANOVA (additive model)                                                                                                    
  variety<-factor(variety)
  envir<-factor(envir)
  block<-factor(block)  
  add.anova<-anova(lm(yield~envir+block%in%envir+variety+variety:envir))
  row.names(add.anova)[3]<-"Block(Environment)"
  row.names(add.anova)[4]<-"Environment x Variety"

  ## This part calculates the ANOVA effects.
  numdata<-length(int.mean)
  envir.tab<-matrix(c(as.vector(manova(int.mean~1)$coefficients)),var.num,envir.num,byrow=TRUE)
  variety.tab<-t(matrix(c(as.vector(manova(t(int.mean)~1)$coefficients)),envir.num,var.num,byrow=TRUE))
  variety.eff<-variety.tab-mean(int.mean)
  envir.eff<-envir.tab-mean(int.mean)
  int.eff<-int.mean-(mean(int.mean)+envir.eff+variety.eff)
  
  ## This part performs a factor analytic decomposition on interaction matrix (main symbols from Gollob, 1968)
  i<-array(c(1:PC,1:PC),dim=c(PC,2))
  D<-matrix(c(0),PC,PC)
  D[i]<-sqrt(eigen(t(int.eff)%*%int.eff)$values[1:PC])
  B<-eigen(t(int.eff)%*%int.eff)$vectors[,1:PC]
  A<-int.eff%*%B[,1:PC]%*%solve(D[1:PC,1:PC])
  E<-B%*%sqrt(D)
  G<-A%*%sqrt(D)
  Ecolnumb<-c(1:PC)
  Ecolnames<-paste("PC",Ecolnumb,sep="")
  dimnames(E)<-list(dimnames(int.mean)[[2]],Ecolnames)
  colnames(G)<-Ecolnames
  
  ## This part tests the significance of PCs (liberal F test)
  PCval<-eigen(t(int.eff)%*%int.eff)$values[1:PC]
  numblock<-length(levels(block))  
  SS<-PCval*numblock
  SS[PC+1]<-add.anova$Sum[4]-sum(SS)
  DF<-var.num+envir.num-1-2*Ecolnumb
  DF[PC+1]<-add.anova$Df[4]-sum(DF)
  MS<-SS/DF
  F<-MS/add.anova$Mean[5]
  probab<-pf(F,DF,add.anova$Df[5],lower.tail=FALSE)
  PercSS<-SS/add.anova$Sum[4]
  rowlab<-c(Ecolnames,"Residuals")
  mult.anova<-data.frame("Effect"=rowlab,"SS"=SS,"DF"=DF,"MS"=MS,"F"=F,"Prob."=probab)
  
  ## This part drows the biplots
  if(biplot==1) {
                plot(envir.mean,E[,1],"n",col="red",lwd=5,xlab="Yield",ylab="PC 1",ylim=c(-1,1))
                text(envir.mean,E[,1],labels=row.names(envir.mean),adj=c(0.5,0.5),col="red")
                points(var.mean,G[,1], "n", col="blue",lwd=5)
                text(var.mean,G[,1],labels=row.names(var.mean),adj=c(0.5,0.5),col="blue")
                abline(h=0, lty=5)
                abline(v=overall.mean, lty=5)
                }
  else {
                plot(E[,1],E[,2],"n",col="red",lwd=5,xlab="PC 1",ylab="PC 2",xlim=c(-1,1),ylim=c(-1,1))
                text(E[,1],E[,2],labels=row.names(E),adj=c(0.5,0.5),col="red")
                points(G[,1],G[,2], "n", col="blue",lwd=5)
                text(G[,1],G[,2],labels=row.names(G),adj=c(0.5,0.5),col="blue")
       }
       
  ## This part displays results  
  list("Variety_means"=var.mean, "Environment_means"=envir.mean, "Interaction_means"=int.mean, "Additive_ANOVA"=add.anova, "Multiplicative_Interaction" = mult.anova, "Environment_scores"=E, "Variety_scores"=G)
  
}