## Polar Ordination (aka Bray-Curtis Ordination) ## Mark Gardener 2015 ## www.dataanalytics.org.uk ## Main polar.ord function polar.ord <- function(diss) { # diss = Distance matrix, e.g. from dist() or vegdist() if(inherits(diss, what = "dist") == FALSE) stop("\nInput must be a dist object\n") dMeth <- attributes(diss)\$method # dissimilarity metric D <- as.dist(diss) # Make sure diagonal and upper.triangle removed M <- as.matrix(D) # Make a matrix from dist object nrows <- nrow(M) # How many rows M0 <- M # Make a copy of the matrix M0[upper.tri(M0)] <- NA # Set upper triangle to NA diag(M0) <- NA # Set diagonals to NA Mi <- order(M0, na.last = NA, decreasing = TRUE) # Get index for lower triangle lD <- length(Mi) # The number of distances # Function to work out co-ordinates co.ord <- function(pos, rows) { if(pos %% rows == 0) { row.no <- rows col.no <- pos %/% rows } if(pos %% rows != 0) { row.no <- pos %% rows col.no <- pos %/% rows +1 } out <- c(row.no, col.no) } # Make matrix for co-ordinates coord <- matrix(nrow = lD, ncol = 2) # row for each dist, col for r-value, c-value for(a in 1:lD) { coord[a,1] <- co.ord(Mi[a], nrows)[1] # Row No. coord[a,2] <- co.ord(Mi[a], nrows)[2] # Col No. } # Make co-ordinates into data.frame and add labels for potential Axes coord <- data.frame(coord) colnames(coord) <- c("R", "C") coord\$End1 <- rownames(M0)[coord\$R] coord\$End2 <- rownames(M0)[coord\$C] coord\$Axis <- paste(coord\$End1, coord\$End2, sep = ".") # Axis labels coord\$dist <- numeric(nrow(coord)) for(i in 1:nrow(coord)) { # get original dist for each pair coord\$dist[i] <- M0[coord\$R[i], coord\$C[i]] } AxDist <- coord[,5:6] # Axis names and original distances AxNo <- nrow(coord) # How many potential Axes there are # All Axis scores Axes <- matrix(nrow = nrows, ncol = AxNo) # To hold results for (b in 1:AxNo) { Dm <- M[coord[b,1], coord[b,2]] Di <- coord\$R[b] Dj <- coord\$C[b] for (a in 1:nrows) { Axes[a,b] <- (Dm^2 + M[a,Di]^2 - M[a,Dj]^2) / (2*Dm) } } colnames(Axes) <- paste(coord\$End1, coord\$End2, sep = ".") rownames(Axes) <- rownames(M0) Axes <- data.frame(Axes) # Get correlations AxCor <- cor(Axes[1], Axes) AxCor[1] <- 0 # Set 1st to zero Axes <- Axes[, order(abs(AxCor))] # Re-order Axes according to correl to 1st AxDist <- AxDist[order(abs(AxCor)),] # Re-order axes and original distances AxD <- as.numeric(AxDist\$dist) names(AxD) <- AxDist\$Axis AxCor <- cor(Axes[1], Axes) # Get correlations again # Variance Explained axd = apply(Axes, 2, FUN = dist, method = "man") # Get all axis distances d = matrix(D, ncol = 1) # Original distances as matrix column dd = cbind(d, axd) colnames(dd)[1] <- "original" # All residuals cf Original options(warn = -1) # Turn off warnings for NaN resA <- matrix(ncol = ncol(dd)-1, nrow = nrow(dd)) colnames(resA) <- colnames(dd)[-1] for(i in 1:ncol(resA)) { resA[, i] <- sqrt(dd[, 1]^2 - dd[, i+1]^2) # With NaN } options(warn = 0) # Restore warnings options # All var explained cf Original varA <- vector(mode = "numeric", length = ncol(resA)) names(varA) <- colnames(resA) for(i in 1:ncol(resA)) { varA[i] <- 1 - sum(resA[, i]^2, na.rm = TRUE) / sum(dd[, 1]^2) } # Get cumulative distances disC <- matrix(ncol = ncol(Axes), nrow = ncol(Axes)) colnames(disC) <- colnames(Axes) for(i in 1:ncol(Axes)) { disC[,i] <- dist(Axes[,1:i]) } disC <- cbind(d, disC) colnames(disC)[1] <- "original" # Cumulative residuals options(warn = -1) resC <- matrix(ncol = ncol(disC)-1, nrow = nrow(disC)) colnames(resC) <- colnames(disC)[-1] for(i in 1:ncol(resC)) { resC[,i] <- sqrt(disC[,1]^2 - disC[,i+1]^2) } options(warn = 0) # Cumulative Variance Explained SSc <- colSums(resC^2, na.rm = TRUE) SStot <- sum(d[,1]^2) varC <- numeric(ncol(resC)) names(varC) <- colnames(resC) for(i in 1:ncol(resC)) { varC[i] <- 1-SSc[i]/SStot } # Variance accumulation varAccum <- numeric(length = length(varC)) names(varAccum) <- names(varC) varAccum[1] <- varC[1] for(i in 2:length(varC)) { varAccum[i] <- (varC[i] - varC[i-1]) } # Results result <- list(sites = Axes, # Site scores dist = AxD, # Original distances cor = AxCor, # Axis correlations var.exp = varA, # Var Expl by axis var.add = varAccum, # Addtl var exp by axis samples = rownames(Axes), # Sample names axes = colnames(Axes), # Axis names diss = D, # Original diss method = dMeth) class(result) <- "polar.ord" invisible(result) } ## END polar.ord Function ## Summary for class polar.ord summary.polar.ord <- function(ord, k = 2, cor = TRUE, var.exp = TRUE, digits = getOption("digits")) { # ord = polar.ord result # k = No. axes to show # cor = show correlations? # var.exp = show variance explained each axis? # digits = No. digits to display if(inherits(ord, what = "polar.ord") == FALSE) stop("\nResult must be from polar.ord command\n") Ax.Tot <- length(ord\$sites) if(k > Ax.Tot) k <- Ax.Tot site.scores <- ord\$sites[,1:k] if(cor == TRUE) { ax.cor <- ord\$cor[,1:k] ax.cor[1] <- NA } if(var.exp == TRUE) ax.var <- ord\$var.add[1:k] Sample.Tot <- length(ord\$samples) method <- ord\$method cat("\nSamples:", Sample.Tot) cat("\nAxes available:", Ax.Tot) cat("\nMethod: ", method) cat("\n\nSite Scores:\n\n") print(site.scores, digits = digits) if(cor == TRUE) { cat("\nAxis correlation:\n\n") print(ax.cor, digits = digits) } if(var.exp == TRUE) { cat("\nAdditional variance explained:\n\n") print(ax.var, digits = digits) } } ## END summary.polar.ord ## Plot method for class polar.ord plot.polar.ord <- function(ord, x.axis = 1, y.axis = 2, type = c("t", "p", "n"), ...) { # ord = result from polar.ord # x.axis = which axes to plot on x-axis # y.axis = which axis to plot on y-axis # ... = other args to plot() if(inherits(ord, what = "polar.ord") == FALSE) stop("\nResult must be from polar.ord command\n") if(any(c(x.axis, y.axis)> length(ord\$sites)) == TRUE) stop("\nThere are only ", length(ord\$sites), " axes available!\n") type <- match.arg(type) dat <- ord\$sites[, c(x.axis, y.axis)] # get axis scores to plot lab <- ord\$samples # labels if(type == "t") { plot(dat, type = "n", ...) text(dat, labels = lab, ...) } if(type == "p") { plot(dat, type = "p", ...) } if(type == "n") { plot(dat, type = "n", ...) } result <- list(x = dat[,1], y = dat[,2], labels = lab) class(result) <- "plot.polar.ord" invisible(result) } # End plot.polar.ord ## Identify method for class polar.ord identify.plot.polar.ord <- function(pord, labels = pord\$labels, ...) { # plot = plot from plot.polar.ord # labels = vector of labels # ... = other arguments if(inherits(pord, what = "plot.polar.ord") == FALSE) stop("\nPlot must be from polar.ord\n") x <- pord\$x y <- pord\$y labels <- pord\$labels identify(x, y, labels, ...) } # End identify.polar.ord # Print method for polar.ord print.polar.ord <- function(ord, digits = getOption("digits")) { res <- data.frame(Axis = names(ord\$dist), dist = ord\$dist, var.exp = ord\$var.exp, row.names = NULL) print(res, digits = digits) } # End print.polar.ord ## END