GGE<-function(variety,envir,block,yield,PC=2){

  ## GENOTYPE PLUS GENOTYPE X INTERACTION, AS APPLIED TO AGRICULTURE VARIETY TRIALS
  ## REFERENCE: Yan et al., 2000. Cultivar evaluation and mega-environment investaigation
  ## based on the GGE biplot
  ## Update: 30/10/2007

  ## 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)

  ## 1 - 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)
  
  ## 2 - ANOVA (additive model) and table of GE effects
  variety<-factor(variety)
  envir<-factor(envir)
  block<-factor(block)  
  add.anova<-aov(yield~envir+block%in%envir+variety+envir:variety)
  add.anova.residual.SS<-deviance(add.anova)
  add.anova.residual.DF<-add.anova$df.residual
  add.anova.residual.MS<-add.anova.residual.SS/add.anova.residual.DF
  anova.table<-summary(add.anova)
  row.names(anova.table[[1]])<-c("Environments","Genotypes","Blocks(Environments)","Environments x Genotypes", "Residuals")
  ge.anova<-aov(yield~envir+envir:variety)
  ge.eff<-model.tables(ge.anova,type="effects",cterms="envir:variety")$tables$"envir:variety"
  
  ## 3 - SVD
  dec <- svd(ge.eff, nu = PC, nv = PC)
  if (PC > 1){ 
    D <- diag(dec$d[1:PC]) 
  } else {
    D <- dec$d[1:PC]
  }
  E<-dec$u%*%sqrt(D)
  G<-dec$v%*%sqrt(D)
  Ecolnumb<-c(1:PC)
  Ecolnames<-paste("PC",Ecolnumb,sep="")
  dimnames(E)<-list(levels(envir),Ecolnames)
  dimnames(G)<-list(levels(variety),Ecolnames)
 
  ## 4 - Singular values and significance of PCs  
  svalues<-dec$d
  PC.n<-c(1:length(svalues))
  PC.SS<-svalues^2
  percSS<-PC.SS/sum(PC.SS)*100
  GGE.table<-data.frame("PC"=PC.n,"Singular_value"=svalues,"PC_SS"=PC.SS, "Perc_of_Total_SS"=percSS)
  numblock<-length(levels(block))
  GGE.SS<-(t(as.vector(ge.eff))%*%as.vector(ge.eff))*numblock
  
  ## 5 - GGE Biplot
  plot(1,type='n',xlim=range(c(E[,1],G[,1])),ylim=range(c(E[,2],G[,2])),xlab="PC 1",ylab="PC 2")
  points(E[,1],E[,2],"n",col="red",lwd=5)
  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")

  ## 6 - Results
  list(Genotype_means=var.mean, Environment_means=envir.mean, Interaction_means=int.mean, 
  GE_effect=ge.eff, Additive_ANOVA=anova.table, "GGE_Sum of Squares"=GGE.SS, 
  GGE_summary = GGE.table, Environment_scores=E, "Genotype_scores"=G)

}