R/arch_varmat.R

In esteful/ArchFlow: Archaeometric Analysis of Ceramics

#' Compositional Variation Matrix.
#' Returns the variation matrix based on the chemical evenness of the dataset. 
#' Prints the less varying variable to use as divisor in alr transformation
#' From Buxeda i Garrigós evariation.matrix2 R function

#' @param df_chem chemical data

#' @return Variation matrix based on \code{df_chem} dataset.

#' @export


"arch_varmat"<-  function(df_chem)
  
{
  #create the matrix 
  p <- dim(df_chem)[2]
  varmat <- matrix(0, p, p)
  varmat2<-matrix(0,p+4,p)
  
  #add values to the matrix (the log of the df_chem/)
  for(i in 1:p) {
    varmat[, i] <- diag(var(log(df_chem/df_chem[, i])))
  }
  varsum <- apply(varmat, 2, sum) #sumatorio
  totvar <- sum(varmat)/(2 * p)
  varprop <- totvar/varsum
  varcor<-vector(mode="numeric",length=p)
  for(i in 1:p) {
    varcor[i] <- compositions::cor(varmat[-c(i),i],varsum[-i])
  }

  par(mar=c(5,5,4,2)+0.1,mgp=c(3,1,0))
  
  #add the total variation, etc
  for(i in 1:p) varmat2[i,]<- varmat[i,]
  varmat2[p+1,]<- varsum
  varmat2[p+2,]<- varprop
  varmat2[p+3,]<- varcor
  varmat2[p+4,1]<- totvar
  dimnames(varmat2)<-list(c(dimnames(df_chem)[[2]],"t.i","vt/t.i","r v,t","vt"),c(dimnames(df_chem)[[2]]))
  cat("%Values of vt/t.i (in Buxeda's notation).\n\n
  %Maximum determines variable to use as divisor in subsequent log-ratio transformation.\n\n\n", "%", varprop, "\n\n")
  ord <- order(varprop)
  cat("%Maximum value : ", max(varprop), "\n")
  cat("%Variable : ", ord[p], "\n")
  

  
  #return a the variation matrix
  return(varmat2)
}