# Mixture 0.0, Jul 21, 2003 # Mixture 0.1, Sep 26, 2003 # Mixture 0.2, Jul 29, 2005 # MixeR 0.2b, Aug 14, 2005 # MixeR 0.2c, Dec 28, 2005 # # # checks the mix object for the rows sums and determines the status m$sta of the mix object # mix.Check <- function(m, eps=1e-6){ m$sum <- NA if(any(m$mat < 0)) m$sta <- -2 # matrix contains negative elements else{ s <- drop(m$mat %*% rep(1, ncol(m$mat))) if( any(s<=0)) m$sta <- -1 # zero sum row exists else{ if(any(m$mat = 0, output is the new mixture with the m$sum = c # default value for c=1 # mix.Normalize <- function(m, c=1){ if (m$sta<0){cat("Negative elements or zero row sum in the data", '\n')} else { m$mat <- c*(m$mat / drop(m$mat %*% rep(1, ncol(m$mat)))) mix.Check(m) } } # # concatenate matrix 'a' and the title 't' into mix object # mix.Matrix <- function(a,t) { mix.Check(structure(list(tit=t,sum=NA,sta=NA, mat=a,class=c("mixture")))) } # # constructs the random mix object with 'nr' rows and 'nc' columns and constant row sum s # mix.Random <- function(nr, nc, s=1) { m <- mix.Normalize(structure(list(tit="random",sum=NA,sta=0, mat=matrix(runif(nr*nc),nr,nc)),class=c("mixture"))) m$mat <- s*m$mat; m$sta <- NA mix.Check(m) } # ==========================================================================# # subcompositions of mixture objects # =========================================================================# # the mix object computed as subcomposition of 'm' without the k=(k_1,...,k_r) columns # normalized if Normalize=TRUE # mix.Sub <- function(m, k, Normalize=TRUE) { if(m$sta<0) {cat("Negative elements or zero row sum in the data", '\n')} else {m$mat <- m$mat[,-k] if(Normalize==TRUE) {mix.Normalize(m)} else { mix.Check(m)} } } # # the mix object computed as subcomposition of 'm' with only the k=(k_1,...,k_r) columns # normalized if Normalize=TRUE # mix.Extract <- function(m, k, Normalize=TRUE) { if(m$sta<0) {cat("Negative elements or zero row sum in the data", '\n')} else {m$mat <- m$mat[, k] if(Normalize==TRUE) {mix.Normalize(m)} else { mix.Check(m)} } } # compute the subcomposition with the k=(k_1,...,k_r) columns # all the rest is amalgamated in the residual # output is the normalized mixture object with the r+1 columns, # r columns k_1,...,k_r extracted and all the rest is summarized in the last as residual # mix.ExtractRes <- function(m, k){ if(m$sta>=0){ residual <- apply(m$mat[,-k], 1, sum) m$mat <- cbind(m$mat[,k], residual) mix.Normalize(m) } } # comput the subcomposition without the k=(k_1,...,k_r) columns amalgamated in the residual # output is the normalized mixture object with the k columns summarized in the residual # mix.SubRes <- function(m, k){ if(m$sta>=0){ residual <- 0 for( i in k) residual <- residual+m$mat[,i] m$mat <- cbind(m$mat[,-k], residual) mix.Normalize(m) } } # ========================================================================== # Aitchison geometry on mixture objects # ======================================================================== # is.mixObject <- function(m){ if (m$sta<0){is.mixObject <- FALSE ; is.mixObject } else { is.mixObject <- TRUE ; is.mixObject} } # Barycenter of adequate dimension mix.Unit <- function(m){ e <- mix.Matrix(matrix(unit.dim(ncol(m$mat)),nrow=1),"Barycenter") e } # Barycenter of dimension n mix.Unitn <- function(n){ e <- mix.Matrix(matrix(unit.dim(n),nrow=1),"Barycenter of dimension n") e } # inverse data mix.InvData<- function(m){ mix.Matrix(inv(m$mat), "Inverse of the data") } # inverse mixture mix.InvMix<- function(m){ if(m$sta < 1) cat("Data contains negative or/and zero elements", '\n') else{ mix.Normalize(mix.Matrix(inv(m$mat), "Mixtures Inverses of the data")) }} # Geometric mean of the data in the mix object mix.Gmean <- function(m){ G <- as.matrix(m$mat[1,]) G <- t(G) {for (i in 1:ncol(m$mat)) G[1,i] <- Gmean(m$mat[,i])} mix.Matrix(as.comp(G),"Geometric mean of the data") } # inverse Geometric mean mix.InvGmean <- function(m){ if(m$sta < 1) cat("Data contains negative or/and zero elements", '\n') else{ mix.Matrix(as.comp(inv(mix.Gmean(m)$mat)), "Inverse geometric mean of the data")} } # Perurbation with inverse Geometric Mean -- Centered Data mix.Center <- function(m){ if(m$sta < 1) cat("Data contains negative or/and zero elements", '\n') else{ G <- mix.InvGmean(m) m1 <- mix.Normalize(m) G1 <- mix.Normalize(mix.Matrix(t(t(m1$mat)*drop(G$mat)),"Centered Data")) G1 }} # ========================================================================== # ploting mixture objects into ternary diagram # ========================================================================= # # plots mix object 'm' into ternary diagram # ternary() is from the on-line answers. # mix.Ter <- function(m, lcex = 1, add = FALSE, ord = 1:3, ...) { if(m$sta>=0){ X <- mix.Normalize(m)$mat cn <- dimnames(X)[[2]][ord] X <- matrix(data=X[,ord],ncol=3,dimnames=list(dimnames(X)[[1]],cn)) s3 <- sqrt(1/3) if(!add){ oldpty <- par("pty"); on.exit(par(pty=oldpty)); par(pty="s") plot(c(-s3, s3), c(0.5-s3, 0.5+s3), type="n", axes=FALSE, xlab="", ylab="") polygon(c(0, -s3, s3), c(1, 0, 0), density=0) if(length(cn) < 3) cn <- as.character(1:3) eps <- 0.05 * lcex text(c(0, s3+eps*0.7, -s3-eps*0.7), c(1+0.1*eps, -0.9*eps, -0.9*eps), cn, cex=lcex, pos=c(3,2,4)) } points((X[,2] - X[,3])*s3, X[,1], ...) } } # # ternary diagram and No/Centered No/Borders No/Gmean # mix.Ternary <- function(m, dist=c(.05,.05,.05), distG=c(.05,.05,.05), ..., cls=c("green","green3","yellowgreen"), Center=FALSE, Borders=FALSE, Gmean=FALSE) { s3 <- 1/sqrt(3) if(Center == FALSE) { mix.Ter(m,...) if(Borders == TRUE) { y <- c(min(m$mat[,1]),max(m$mat[,1])) for (p in 1:2){ lines(c(s3*(y[p]-1),s3*(1-y[p])),c(y[p],y[p]), lty = 4, col=cls[1]) text(s3*(1-y[p])+ dist[1],y[p], round(y[p]*100,1), col=cls[1])} y <- c(min(m$mat[,2]),max(m$mat[,2])) for (p in 1:2){ lines(c(s3*(2*y[p]-1),s3*y[p]),c(0,1-y[p]),lty=3, col=cls[2]) text(s3*(2*y[p]-1),0-dist[2], round(y[p]*100,1), col=cls[2])} y <- c(min(m$mat[,3]),max(m$mat[,3])) for (p in 1:2){ lines(c(s3*(1-2*y[p]),s3*(-y[p])),c(0,1-y[p]),lty=2, col=cls[3]) text(s3*(-y[p])-dist[3],1-y[p]+dist[3], round((y[p])*100,1), col=cls[3])} } if(Gmean == TRUE) { G <- mix.Gmean(m)$mat points(s3*(G[1,2]-G[1,3]),G[1,1], col="red", pch=13, cex=2) lines(c(s3*(G[1,1] -1),s3*(1-G[1,1])),c(G[1,1],G[1,1]), lty = 4, col="red") text(s3*(1-G[1,1])+ distG[1],G[1,1], round(G[1,1]*100,1), col="red") lines(c(s3*(2*G[1,2]-1),s3*G[1,2]),c(0,1-G[1,2]),lty=3, col="red") text(s3*(2*G[1,2]-1),-distG[2], round(G[1,2]*100,1), col="red") lines(c(s3*(1-2*G[1,3]),s3*(-G[1,3])),c(0,1-G[1,3]),lty=2, col="red") text(s3*(-G[1,3])-distG[3],1-G[1,3]+distG[3], round(G[1,3]*100,1), col="red") } } else { mix.Ter(mix.Center(m),...) if(Borders == TRUE) { y <- c(min(m$mat[,1]),max(m$mat[,1])) y1 <- c(min(mix.Center(m)$mat[,1]),max(mix.Center(m)$mat[,1])) for (p in 1:2){ lines(c(s3*(y1[p]-1),s3*(1-y1[p])),c(y1[p],y1[p]), lty = 4, col=cls[1]) text(s3*(1-y1[p])+ dist[1],y1[p], round(y[p]*100,1), col=cls[1])} y <- c(min(m$mat[,2]),max(m$mat[,2])) y1 <- c(min(mix.Center(m)$mat[,2]),max(mix.Center(m)$mat[,2])) for (p in 1:2){ lines(c(s3*(2*y1[p]-1),s3*y1[p]),c(0,1-y1[p]),lty=3, col=cls[2]) text(s3*(2*y1[p]-1),0-dist[2], round(y[p]*100,1), col=cls[2])} y <- c(min(m$mat[,3]),max(m$mat[,3])) y1 <- c(min(mix.Center(m)$mat[,3]),max(mix.Center(m)$mat[,3])) for (p in 1:2){ lines(c(s3*(1-2*y1[p]),s3*(-y1[p])),c(0,1-y1[p]),lty=2, col=cls[3]) text(s3*(-y1[p])-dist[3],1-y1[p]+dist[3], round((y[p])*100,1), col=cls[3]) } if(Gmean == TRUE) { G <- mix.Gmean(m)$mat points(0,.33, col="red", pch=13, cex=2) lines(c(s3*(.33-1),s3*(1-.33)),c(.33,.33), lty = 4, col="red") text(s3*(1-.33)+ distG[1],.33, round(G[1,1]*100,1), col="red") lines(c(s3*(2*.33-1),s3*.33),c(0,1-.33),lty=3, col="red") text(s3*(2*.33-1),-distG[2], round(G[1,2]*100,1), col="red") lines(c(s3*(1-2*.33),s3*(-.33)),c(0,1-.33),lty=2, col="red") text(s3*(-.33)-distG[3],1-.33+distG[3], round(G[1,3]*100,1), col="red") } } } } # # Percentile lines # perc.lines <- function(y, direction = 1:3 , cls = yellow, dist =c(.05,.05,.05), lt=c(4,3,2)){ s3 <- 1/sqrt(3) yellow <- c("yellow" , "yellow2", "yellow3") {if(y<1) y1 <- y else {y1 <- y/100} } perc.linesM <- function(y, i){ if (i == 1){ for (p in 1:length(y)){ lines(c(s3*(y1[p]-1),s3*(1-y1[p])),c(y1[p],y1[p]), lty=lt[1], col=cls[1]) text(s3*(1-y1[p])+ dist[1],y1[p], round(y1[p]*100,1), col=cls[1])} } if (i == 2){ for (p in 1:length(y)){ lines(c(s3*(2*y1[p]-1),s3*y1[p]),c(0,1-y1[p]), lty=lt[2], col=cls[2]) text(s3*(2*y1[p]-1),0-dist[2], round(y1[p]*100,1), col=cls[2])} } if (i == 3){ for (p in 1:length(y)){ lines(c(s3*(1-2*y1[p]),s3*(-y1[p])),c(0,1-y1[p]), lty=lt[3], col=cls[3]) text(s3*(-y1[p])-dist[3],1-y1[p]+dist[3], round((y1[p])*100,1), col=cls[3])} } } for (i in direction) perc.linesM(y,i) } # # Ratio lines # ratio.lines <- function(y, direction = 1:3 , cls =green , dist =c(.05,.05,.05), lt=c(4,3,2)){ s3 <- 1/sqrt(3) green <- c("green","green3","yellowgreen" ) ratio.linesM <- function(y, i){ if (i == 1){ for (p in 1:length(y)){ lines(c(s3*(y[p]-1)/(y[p]+1),0),c(0,1), lty = lt[1], col=cls[1]) text( s3*(y[p]-1)/(y[p]+1),-dist[1] , round(y[p],2), col=cls[1])} } if (i == 2){ for (p in 1:length(y)){ lines(c(-s3/(1+y[p]),s3),c(y[p]/(1+y[p]),0),lty=lt[2], col=cls[2]) text(-s3/(1+y[p])-dist[2],y[p]/(1+y[p]) , round(y[p],2), col=cls[2])} } if (i == 3){ for (p in 1:length(y)){ lines(c(s3/(1+y[p]),-s3),c(y[p]/(1+y[p]),0),lty=lt[3], col=cls[3]) text(s3/(1+y[p])+dist[3],y[p]/(1+y[p]), round(y[p],2), col=cls[3])} } } for (i in direction) ratio.linesM(y,i) } # ========================================================================== # ploting mixture objects into tetrahedron # ========================================================================= # # plots mix object 'm' into tetrahedron # mix.Quad2Net <- function(fnet,m){ net <- file(fnet,"w") n <- nrow(m$mat) cat("*vertices",n+4,"\n",file=net) rn <- c(dimnames(m$mat)[[1]],"A","B","C","D") a <- 10/9+c(m$mat[,1],1,0,0,0) b <- c(m$mat[,2],0,1,0,0) c <- c(m$mat[,3],0,0,1,0) d <- c(m$mat[,4],0,0,0,1) x <- (a - b - c + d)*0.45 y <- (a - b + c - d)*0.45 z <- (a + b - c - d)*0.45 for (i in 1:length(rn)){ cat(i," \"",rn[i],"\" ",x[i]," ",y[i]," ",z[i],"\n",file=net,sep="") } cat("*edges\n",n+1,n+2,"\n",n+1,n+3,"\n",n+1,n+4,"\n", n+2,n+3,"\n",n+2,n+4,"\n",n+3,n+4,"\n",file=net) close(net) # cbind(x,y,z) } # ========================================================================== # transforms rows of 4-parts mixture objects into XYZ coordinates # ========================================================================= # # transforms rows of 4-parts mixture m quadrays into 3-dimensional # XYZ coordinates and writes them as file.kin. # to be displayed as points in a tetrahedron with KiNG viewer a free software # clu - partition -> color of points # vec - vector of values -> size of points # col - color of points if clu=NULL # scale - relative size of points # king - FALSE (for Mage); TRUE (for King) mix.Quad2kin <- function(kinfile,m,clu=NULL,vec=NULL,king=TRUE,scale=0.2,col=1){ kinColors <- c('white', 'red', 'blue', 'yellow', 'green', 'cyan', 'magenta', 'pink', 'lime', 'orange', 'peach', 'gold', 'purple', 'sea', 'gray', 'sky', 'brown', 'lilac', 'hotpink', 'yellowtint', 'pinktint', 'peachtint', 'greentint', 'bluetint', 'lilactint', 'deadwhite', 'deadblack', 'invisible') warn <- "" if (king) warn <- "\n*** works only with KiNG viewer: http://kinemage.biochem.duke.edu/" head <- paste("@text\n", kinfile," / ",date(), "\nDataset: ", m$tit,warn, "\nLayout obtained using MixeR http://vlado.fmf.uni-lj.si/pub/MixeR/ Vladimir Batagelj & Matevz Bren Institute of Mathematics, Physics and Mechanics Ljubljana, Slovenia @kinemage 1 @caption\n", m$tit, "\n@fontsizeword 18 @zoom 1.00 @thinline @zclipoff @group {Complete} dominant animate movieview = 1 @spherelist 0 radius= 0.20\n",sep="") foot <- "@vectorlist {} color= blue P 9.500, 0.500, 9.500 0.500, 9.500, 9.500 P 9.500, 0.500, 9.500 0.500, 0.500, 0.500 P 9.500, 0.500, 9.500 9.500, 9.500, 0.500 P 0.500, 9.500, 9.500 0.500, 0.500, 0.500 P 0.500, 9.500, 9.500 9.500, 9.500, 0.500 P 0.500, 0.500, 0.500 9.500, 9.500, 0.500\n" kin <- file(kinfile,"w") n <- nrow(m$mat) if (is.null(clu)) {clu <- rep(col,n)} clu <- c(clu,0,0,0,0) if (is.null(vec)) {vec <- rep(1,n)} vec <- c(vec,1,1,1,1) size <- scale/max(vec) rn <- c(dimnames(m$mat)[[1]],"A","B","C","D") a <- 10/9+c(m$mat[,1],1,0,0,0) b <- c(m$mat[,2],0,1,0,0) c <- c(m$mat[,3],0,0,1,0) d <- c(m$mat[,4],0,0,0,1) x <- (a - b - c + d)*0.45 y <- (a - b + c - d)*0.45 z <- (a + b - c - d)*0.45 cat(head,file=kin) for (i in 1:length(rn)) { color <- "white" if (clu[i]>0) color <- kinColors[2+(clu[i]-1)%%18] cat("{",rn[i],"} ", color," ",file=kin,sep="") if (king) cat(" r=", vec[i]*size," ",file=kin,sep=" ") cat(10*x[i],10*(1-y[i]),10*z[i],"\n",file=kin,sep=" ") } cat(foot,file=kin) close(kin) } # ========================================================================== # transformations from 'compositions' classes to mixtures # ========================================================================= # transform 'compositions' class to mixture mix.InvComp <- function (X) { t <- c("Mixture from 'compositions' class") m <- X[ , ] mix.Matrix(m,t) } # transform rmult to mixture mix.Inv2rmult <- function (X) { t <- c("Mixture from 'rmult' class") m <- X[ , ] mix.Matrix(m,t) } # transform rplus to mixture mix.Inv2rplus <- function (X) { t <- c("Mixture from 'rplus' class") m <- X[ , ] mix.Matrix(m,t) } # transform rcomp to mixture mix.Inv2rcomp <- function (X) { t <- c("Mixture from 'rcomp' class") m <- X[ , ] mix.Matrix(m,t) } # transform acomp to mixture mix.Inv2acomp <- function (X) { t <- c("Mixture from 'acomp' class") m <- X[ , ] mix.Matrix(m,t) } # transform aplus to mixture mix.Inv2aplus <- function (X) { t <- c("Mixture from 'aplus' class") m <- X[ , ] mix.Matrix(m,t) } # ========================================================================== # transformations from mixtures to 'compositions' classes 'aplus', 'acomp', 'rcomp', 'rplus' # ======================================================================== # transform mixture to aplus mix.2aplus <- function (X) { Y <- X$mat class(Y) <- "aplus" Y } # transform mixture to acomp mix.2acomp <- function (X) { Y <- X$mat class(Y) <- "acomp" Y } # transform mixture to rcomp mix.2rcomp <- function (X) { Y <- X$mat class(Y) <- "rcomp" Y } # transform mixture to rplus mix.2rplus <- function (X) { Y <- X$mat class(Y) <- "rplus" Y } # transform mixture to rmult mix.2rmult <- function (X) { Y <- X$mat class(Y) <- "rmult" Y } # ========================================================================== # reads the CODA data files -- same as MATLAB # ========================================================================= # reads the CODA data file and returns a mixture object, # or if Matrix=TRUE only the matrix with the data # mix.ReadCODA <- function(f, Matrix=FALSE) { at <- read.table(file=f, header=T, sep='\n') nc <- as.integer(as.vector(at[[1]][1])) nr <- as.integer(as.vector(at[[1]][2])) n <- as.vector(at[[1]][4:(nc+3)]) m <- matrix(as.numeric(as.vector(at[[1]][(nc+4):((nr+1)*(nc+1)+2)])), ncol=nc+1, byrow=T)[,-1] dimnames(m)[[2]] <- n dimnames(m)[[1]] <- 1:nr if (!Matrix ) m <- mix.Matrix(m,colnames(at)) m } # ===================================================================== # reads the MATLAB data files -- same as CODA # ===================================================================== # reads the MATLAB data file and returns a mixture object with a column of cases labels, # or if Matrix=TRUE only the matrix with the data # mix.ReadML <- function(f, Matrix=FALSE) { at <- read.table(file=f, header=T, sep='\n') nc <- as.integer(as.vector(at[[1]][1])) nr <- as.integer(as.vector(at[[1]][2])) n <- as.vector(at[[1]][4:(nc+3)]) m <- matrix(as.numeric(as.vector(at[[1]][(nc+4):((nr+1)*(nc+1)+2)])), ncol=nc+1, byrow=T)[,-1] dimnames(m)[[2]] <- n dimnames(m)[[1]] <- 1:nr if (!Matrix ) m <- mix.Matrix(m,colnames(at)) m } # ========================================================================== # checking for zero (i.e. less than eps) and replacement of the zero values in the data # Simple replacement strategy, Math. Geology Vol. 35, No. 3, pp.261 # Multiplicative replacement strategy, Math. Geology Vol. 35, No. 3, pp.262 # ========================================================================= # computs rows or columns sums of the data # if row=TRUE output are the rows sums # if row=FALSE output are the columns sums # mix.Sum <- function(m, row=TRUE){ if (row==TRUE) { s <- apply(m$mat, 1, sum); s} else { s <- apply(m$mat, 2, sum); s} } # # computing rows sums of the data # output are the rows sums # mix.RowSum <- function(m){ s <- drop(m$mat %*% rep(1, ncol(m$mat))); s } # # computing the columns sums of the data # output are the columns sums # mix.ColSum <- function(m){ s <- drop(t(rep(1, nrow(m$mat)) %*% m$mat )); s } # # checking for zero and negative elements in the data # output is column, row, value and No. of negative and/or zero cases # or "No negative, no zero values in the data" # if Detect = TRUE # output is also the mix object with matrix elements equal # -1 at negative cases, equal 0 at zero cases and equal 1 at positive cases # mix.CheckData <- function(m, eps=1e-6, Detect=FALSE){ if (m$sta>0) {cat("No negative, no zero values in the data", '\n')} else { n1 <- 0; n2 <- 0 for( j in 1:ncol(m$mat)) for( i in 1:nrow(m$mat)) { if( m$mat[i,j] < 0) { n1 <- n1+1; cat("negative values in the data ", '\n', '\t',"negative value ", '\t', "column ", j, '\t', "row ", i, '\t', "value ", m$mat[i,j],'\t', "No. ", n1, '\n')} else if( abs(m$mat[i,j]) < eps) { n2 <- n2+1; cat("zero values in the data ", '\n', '\t',"zero value ",'\t', "column ", j,'\t', "row ", i, '\t', "value ", m$mat[i,j],'\t', "No. ", n2, '\n')} } if( n1 == 0) cat("No negative values in the data", '\n') if( n2 == 0) cat("No zero values in the data", '\n') if (Detect==TRUE){ { for( j in 1:ncol(m$mat)) for( i in 1:nrow(m$mat)) if( abs(m$mat[i,j]) < eps) m$mat[i,j] <- 0 else if( m$mat[i,j] < 0) m$mat[i,j] <- -1 else m$mat[i,j] <- 1 } m <- mix.Matrix(m$mat,"Noting -1 for neg., 0 for zero and 1 for pozitive values ") m } } } # # replacement of the zero (i.e. less than eps) components # Simple replacement strategy, Math. Geology Vol. 35, No. 3, pp.261 # all the zero components are replaced, the row sums are preserved # de are inputet values, # if Col=TRUE, one component for each column # if Col=FALSE, one component for each row with zero components # mix.ZeroReplaceSimp <- function(m, de, Col=TRUE, eps=1e-6){ if (m$sta <= -1) {cat(" Negative values or zero row sum in the data ", '\n')} else { if (Col==TRUE) { de<- rep(de,times=ceiling(ncol(m$mat)/length(de))) for( i in 1:nrow(m$mat)) {s <- sum(m$mat[i,]) {n <- 0; sde <- 0 {for( j in 1:ncol(m$mat)) if( abs(m$mat[i,j]) < eps) {n <- n+1; sde <- sde+de[j]}}} if(n > 0) {for( k in 1:ncol(m$mat)) {if( abs(m$mat[i,k]) < eps) m$mat[i,k] <- s*de[k]/(s+sde) else m$mat[i,k] <- s*m$mat[i,k]/(s+sde)}}} } else { r <- 0 for( i in 1:nrow(m$mat)) {s <- sum(m$mat[i,]) {n <- 0 for( j in 1:ncol(m$mat)) if( abs(m$mat[i,j]) < eps) n <- n+1} if(n > 0) {r <- r+1 for( k in 1:ncol(m$mat)) if( abs(m$mat[i,k]) < eps) m$mat[i,k] <- s*de[r]/(s+n*de[r]) else m$mat[i,k] <- s*m$mat[i,k]/(s+n*de[r])}} } mix.Check(m) } } # # replacement of the zero (i.e. less than eps) data # Multiplicative replacement strategy, Math. Geology Vol. 35, No. 3, pp.262 # all the zero data are replaced, the row sums are preserved # de is the vector of inputet values, one component for each column # mix.ZeroReplaceMult <- function(m, de, eps=1e-6){ if (m$sta <= -1) {cat(" Negative values or zero row sum in the data ", '\n')} else { de<- rep(de,times=ceiling(ncol(m$mat)/length(de))) for( i in 1:nrow(m$mat)) {s <- sum(m$mat[i,]) {n <- 0; sde <- 0 {for( j in 1:ncol(m$mat)) if( abs(m$mat[i,j]) < eps) {n <- n+1; sde <- sde+de[j]}}} if(n > 0) {for( k in 1:ncol(m$mat)) {if( abs(m$mat[i,k]) < eps) m$mat[i,k] <- de[k] else m$mat[i,k] <- (1-sde/s)*m$mat[i,k]}} } mix.Check(m) } } # ========================================================================== # distances and dissimilarityes on mix object, with the R dist output # ========================================================================= # # Aitchison distance between cases i.e. rows of the matrix of mixture object # mix.AitDist <- function(m){ if(m$sta < 1) cat("Data contains negative or/and zero elements", '\n') else{ n <- nrow(m$mat) mE <- lower.tri(array(0,dim=c(n,n))) {for( i in 2:n) for( j in 1:(i-1) ) {mE[i,j] <- sqrt(sum((log(m$mat[i,]/Gmean(m$mat[i,]))-log(m$mat[j,]/Gmean(m$mat[j,])))^2)) }} as.dist(mE)} } # # Martin dissim -version 1- between the rows of the mixture object # mix.MartDist <- function(m){ if(m$sta < 1) cat("Data contains negative or/and zero elements", '\n') else{ n <- nrow(m$mat) nc <- ncol(m$mat) mE <- lower.tri(array(0,dim=c(n,n))) {for( i in 2:n) for( j in 1:(i-1) ) {mE[i,j] <- sqrt((nc/2)*log(mean(m$mat[i,]/m$mat[j,])*mean(m$mat[j,]/m$mat[i,]))) }} as.dist(mE)} } # # Kullbach - Leibler information number between the rows of the mixture object # mix.KuLeDist <- function(m){ if(m$sta < 1) cat("Data contains negative or/and zero elements", '\n') else{ n <- nrow(m$mat) nc <- ncol(m$mat) mE <- lower.tri(array(0,dim=c(n,n))) {for( i in 2:n) for( j in 1:(i-1) ) {mE[i,j] <- sum(m$mat[i,]*(log(m$mat[i,]/m$mat[j,]))) }} as.dist(mE)} } # # Bhattarcharyya distance between the rows of the mixture object # mix.BhatDist <- function(m){ if(m$sta < 1) cat("Data contains negative or/and zero elements", '\n') else{ n <- nrow(m$mat) nc <- ncol(m$mat) mE <- lower.tri(array(0,dim=c(n,n))) {for( i in 2:n) for( j in 1:(i-1) ) {mE[i,j] <- acos(sum(sqrt(m$mat[i,]*m$mat[j,]))) }} as.dist(mE)} } # # Euclidean distance between the rows of the mixture object # mix.EucDist <- function(m){ if(m$sta < 1) cat("Data contains negative or/and zero elements", '\n') else{ n <- nrow(m$mat) nc <- ncol(m$mat) mE <- lower.tri(array(0,dim=c(n,n))) {for( i in 2:n) for( j in 1:(i-1) ) {mE[i,j] <- sqrt(sum((m$mat[i,]-m$mat[j,])^2)) }} as.dist(mE)} } # # Minkovski p - distance between the rows of the mixture object # mix.MinkDist <- function(m,p){ if(m$sta < 1) cat("Data contains negative or/and zero elements", '\n') else{ n <- nrow(m$mat) nc <- ncol(m$mat) mE <- lower.tri(array(0,dim=c(n,n))) {for( i in 2:n) for( j in 1:(i-1) ) {mE[i,j] <- (sum(abs(m$mat[i,]-m$mat[j,])^p))^(1/p) }} as.dist(mE)} } # ========================================================================== # also necessary routines # ======================================================================== # inverse values inv <- function(x) 1/x # unit vector of dim n unit.dim <- function(n) rep(1/n, times=n) # unit vector of dim is length of comp x unit.comp <- function(x) rep(1/length(x), times=length(x)) # is TRUE if the sum of x components equal to 1 -- logical is.composition <- function(x, eps=1e-6 ){ is.numeric(x) & !is.complex(x) & all(x>=0) & abs(sum(x)-1.0) < eps } # definition of the binary perturbation operation, x,y are comp. "%+%" <- function(x, y){ x * y/sum(x*y) } # definition of the binary (power) multiplication, x - comp., k - power "%^%" <- function(x, k){ x^k/sum(x^k) } # definition of the binary perturbation operation with the inverse elt., x,y are comp. "%-%" <- function(x, y) { x %+% (y %^% -1) } # transformation of vector x into comp. as.comp <- function(x) x/sum(x) # inverse comp inv.comp <- function(x) as.comp(1/x) # geometric mean of a comp Gmean <- function(x) prod(x)^(1/length(x)) # harmonic mean of a comp Hmean <- function(x) length(x)/sum(1/x) #===================================================================== APPENDIX # percentile lines using polygon lines only # =========================== # percentile lines to the vertex No. 1 perc.linPoly1 <- function(y){ if(0