diff --git a/DESCRIPTION b/DESCRIPTION index b506c99..6034cf4 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,14 +1,14 @@ Package: netassoc Type: Package -Title: Inference of Species Associations from Co-occurrence Data -Version: 0.4.1 -Date: 2015-01-29 +Title: Inference of Species Associations from Co-Occurrence Data +Version: 0.4.3 +Date: 2015-03-13 Author: Benjamin Blonder, Naia Morueta-Holme Maintainer: Benjamin Blonder -Description: Infers species associations from local and regional-scale co-occurrence data (species x site matrices) +Description: Infers species associations from local and regional-scale co-occurrence data (species x site matrices). License: GPL-3 Depends: igraph -Packaged: 2015-01-30 22:38:14 UTC; benjaminblonder +Packaged: 2015-03-13 17:01:47 UTC; benjaminblonder NeedsCompilation: no Repository: CRAN -Date/Publication: 2015-01-31 08:31:33 +Date/Publication: 2015-03-13 20:22:09 diff --git a/MD5 b/MD5 index e60e672..9f5f17a 100644 --- a/MD5 +++ b/MD5 @@ -1,8 +1,6 @@ -78461170b1f9f8456174c0e718ad3a87 *DESCRIPTION +648150b10eab05c75734ec845e3afd27 *DESCRIPTION 11ce28fad68ae62136749e9a93424939 *NAMESPACE -4b597216d0102f9a654b64b58d075ee1 *R/netassoc.R -b132e71152f9588edd175f90059adf9f *data/fia.rda -6517b6b86290db58559092a573032d24 *man/fia.Rd -0d1bedd3aec4bedf2e2ebe3a39b6a471 *man/makenetwork.Rd -9c4d412f0439a8cc43c0386aa0a381d7 *man/netassoc-package.Rd -ddcfc865933f30160ed3bc4ddf6a31de *man/plot_netassoc_network.Rd +c3806531e8b9051a4f8d142cf8d7184c *R/netassoc.R +9bc661c5f6f4054499e3872f585cd795 *man/makenetwork.Rd +6aa2c25359d4b486497bfa5ce19e4520 *man/netassoc-package.Rd +491d858df7cc04641612582e169b9d42 *man/plot_netassoc_network.Rd diff --git a/R/netassoc.R b/R/netassoc.R index 7bb25dd..ec6d1dd 100644 --- a/R/netassoc.R +++ b/R/netassoc.R @@ -1,115 +1,28 @@ -# plot network -plot_netassoc_network <- function(network, layout = layout.auto(network), - vertex.label = V(network)$name, - vertex.color = NA, - vertex.shape = "none", - vertex.label.color = "black", - vertex.label.family = "sans", - edge.width = NULL, - edge.color = NULL, - edge.arrow.size = 0.05, - vertex.label.cex = 0.5, - ...) -{ - if(is.null(edge.width)) - { - if(length(E(network)$weight)==0) - { - edge.width=1 - } - else - { - edge.width=sqrt(abs(E(network)$weight)) - } - } - - if(is.null(edge.color)) - { - if(length(E(network)$weight)==0) - { - edge.color <- 'black' - zlmin <- -1 - zlmax <- 1 - } - else - { - edge.color <- ifelse(E(network)$weight > 0, rgb(0,0,1,0.8),rgb(1,0,0,0.8)) - zlmax <- max(abs(E(network)$weight),na.rm=T) - zlmin = -1*zlmax - } - } - - plot(network, - layout=layout, - vertex.label=vertex.label, - edge.color=edge.color, - edge.width=edge.width, - vertex.color=vertex.color, - vertex.label.color=vertex.label.color, - vertex.shape=vertex.shape, - edge.arrow.size=edge.arrow.size, - vertex.label.cex=vertex.label.cex, - vertex.label.family=vertex.label.family, - ...) - - colors=colorRampPalette(c('red','white','blue'))(51) - legend('topleft',adj=c(0,0),legend=format(c(zlmin,zlmin/2+zlmax/2,zlmax),digits=2),fill=c(colors[1],colors[ceiling(length(colors)/2)],colors[length(colors)]),bg='white',cex=0.5) -} - -plot_netassoc_matrix <- function(data, colors, ylab='Species',xlab='Site',onesided=FALSE,main="") -{ - zlmax <- max(abs(as.numeric(data)),na.rm=T) - if (is.infinite(zlmax)) - { - zlmax <- 1 - } - if (onesided==TRUE) - { - zlmin = 0 - } - else - { - zlmin = -1*zlmax - } - - image(t(data),col=colors,axes=F,zlim=c(zlmin, zlmax),main=main) - - mtext(side=2,text=ylab,cex=0.5) - mtext(side=3,text=xlab,cex=0.5) - legend('topleft',adj=c(0,0),legend=format(c(zlmin,zlmin/2+zlmax/2,zlmax),digits=2),fill=c(colors[1],colors[ceiling(length(colors)/2)],colors[length(colors)]),bg='white',cex=0.5) - box() -} - - -spcor <- function(matrix, sp1id, sp2id, method) # problems ca -{ - # get pres/abs for each species across sites - vec1 <- matrix[sp1id,] - vec2 <- matrix[sp2id,] - - return(cor(vec1, vec2, method=method)) -} - # generate null sp x site matrices with same total abundance as observed and species probabilities drawn abundances within the base null matrix -generate_nul_resample <- function(nul, obs, verbose=TRUE) +generate_nul_resample <- function(nul, obs) { - nul_resample <- nul + nul_resample <- matrix(0, nrow=nrow(nul),ncol=ncol(nul),dimnames=dimnames(nul)) for (i in 1:ncol(nul_resample)) { - # zero out all values - nul_resample[,i] <- 0 - # fill in resamples with the same total abundance - abund_simulated <- table(sample(names(nul[,i]),size=sum(obs[,i]),replace=T,prob=nul[,i])) - nul_resample[names(abund_simulated),i] <- abund_simulated + samples <- sample(1:nrow(nul),size=sum(obs[,i]),replace=TRUE,prob=nul[,i]) + nul_resample[,i] <- tabulate(samples, nbins=nrow(nul)) } - if (verbose==TRUE) { cat('.') } return(nul_resample) } -makenetwork <- function(obs, nul, method='pearson', kappa=2, numnulls=1000, debug=FALSE,verbose=TRUE) -{ +makenetwork <- function(obs, nul, whichmethod='pearson', kappa=2, numnulls=1000, plot=FALSE,verbose=TRUE) +{ + if (is.matrix(obs) & is.null(dimnames(obs))) + { + dimnames(obs) <- list(1:nrow(obs),1:ncol(obs)) + } + if (is.matrix(nul) & is.null(dimnames(nul))) + { + dimnames(nul) <- list(1:nrow(nul),1:ncol(nul)) + } + obs <- as.matrix(obs,rownames.force=TRUE) nul <- as.matrix(nul,rownames.force=TRUE) @@ -118,13 +31,19 @@ makenetwork <- function(obs, nul, method='pearson', kappa=2, numnulls=1000, debu stop("obs and nul must be same dimensionalities.") } - problemspecies <- which(rowSums(obs)==0 | rowSums(nul)==0) - warning(sprintf("Some species do not occur in any sites: %s",paste(names(problemspecies),collapse=", "))) - - #obs <- obs[-problemspecies,] - #nul <- nul[-problemspecies,] + problemspecies <- which(rowSums(obs)==0 | rowSums(nul)==0 | is.na(rowSums(nul)) | is.na(rowSums(obs))) + if (length(problemspecies) > 0) + { + stop(sprintf("Some species do not occur in any sites: %s",paste(names(problemspecies),collapse=", "))) + } + + problemsites <- which(colSums(obs)==0 | colSums(nul)==0 | is.na(colSums(nul)) | is.na(colSums(obs))) + if (length(problemsites) > 0) + { + stop(sprintf("Some sites are empty: %s",paste(names(problemsites),collapse=", "))) + } - if (debug==TRUE) + if (plot==TRUE) { par(mfrow=c(2,4)) par(cex.lab=0.5) @@ -132,7 +51,7 @@ makenetwork <- function(obs, nul, method='pearson', kappa=2, numnulls=1000, debu par(mar=c(0.5,2,2,0.5)) } - if (debug==TRUE) + if (plot==TRUE) { plot_netassoc_matrix(obs, colors=colorRampPalette(c('white','green'))(51),onesided=TRUE,main="Observed sp x site") @@ -141,15 +60,17 @@ makenetwork <- function(obs, nul, method='pearson', kappa=2, numnulls=1000, debu if (verbose==TRUE) { cat('Generating null replicates...') } - nulls <- replicate(numnulls, generate_nul_resample(nul, obs, verbose=verbose)) + nulls <- vector(mode="list",length=numnulls) + + nulls <- replicate(numnulls, {if(verbose==TRUE) { cat('.') }; return(generate_nul_resample(nul, obs)) }) + if (verbose==TRUE) { cat('...done.\n') } - if (debug==TRUE) + if (plot==TRUE) { - plot_netassoc_matrix(nulls[,,sample(numnulls, 1)], colors=colorRampPalette(c('white','black'))(51),onesided=TRUE,main="Example null resample sp x site") + plot_netassoc_matrix(nulls[,,1], colors=colorRampPalette(c('white','black'))(51),onesided=TRUE,main="Example null resample sp x site") } - fm_obs <- matrix(NA,nrow=nrow(obs),ncol=nrow(obs)) fm_nul_mean <- matrix(NA,nrow=nrow(obs),ncol=nrow(obs)) fm_nul_sd <- matrix(NA,nrow=nrow(obs),ncol=nrow(obs)) @@ -165,26 +86,41 @@ makenetwork <- function(obs, nul, method='pearson', kappa=2, numnulls=1000, debu if (i!=j) { count <- count + 1 - if (verbose==TRUE) { cat (sprintf('%.3f ', count/(nrow(obs)*(nrow(obs)-1)/2))) } - - cor_obs <- spcor(obs, i, j, method=method) + if (verbose==TRUE) { cat (sprintf('%d %d %.3f ', i, j, count/(nrow(obs)*(nrow(obs)-1)/2))) } + if (verbose==TRUE) { cat('obs') } + veci_obs <- as.numeric(obs[i,]) + vecj_obs <- as.numeric(obs[j,]) + cor_obs <- cor(x=veci_obs, y=vecj_obs, method=whichmethod) fm_obs[i,j] <- cor_obs + if (verbose==TRUE) { cat('. ') } + if (!is.na(cor_obs)) { - cor_nul <- rep(NA, dim(nulls)[3]) + cor_nul <- rep(NA, numnulls) - for (k in 1:dim(nulls)[3]) + for (k in 1:numnulls) { - cor_nul[k] <- spcor(nulls[,,k], i, j, method=method) + if (verbose==TRUE) { cat('.') } + veci_nul <- as.numeric(nulls[i,,k]) + vecj_nul <- as.numeric(nulls[j,,k]) + cor_nul[k] <- cor(x=veci_nul, y=vecj_nul, method=whichmethod) + if (verbose==TRUE) { cat('. ') } } fm_nul_mean[i,j] <- mean(cor_nul,na.rm=T) fm_nul_sd[i,j] <- sd(cor_nul,na.rm=T) } + else + { + if (verbose==TRUE) { cat('nul NONE') } + } + if (verbose==TRUE) { cat('\n') } } } } + # save some memory + rm(nulls) if (verbose==TRUE) { cat('...done.\n') } # calculate standard effect size @@ -192,7 +128,7 @@ makenetwork <- function(obs, nul, method='pearson', kappa=2, numnulls=1000, debu finalmatrix <- (fm_obs - fm_nul_mean) / fm_nul_sd if (verbose==TRUE) { cat('...done.\n') } - if (debug==TRUE) + if (plot==TRUE) { plot_netassoc_matrix(fm_obs, xlab='Species',ylab='Species',colors=colorRampPalette(c('red','white','blue'))(51),main="Observed co-occurrence score for sp x sp") plot_netassoc_matrix(fm_nul_mean, xlab='Species',ylab='Species',colors=colorRampPalette(c('red','white','blue'))(51),main="Null mean co-occurrence score for sp x sp") @@ -209,7 +145,7 @@ makenetwork <- function(obs, nul, method='pearson', kappa=2, numnulls=1000, debu finalmatrix_trimmed[is.infinite(finalmatrix_trimmed)] <- 0 if (verbose==TRUE) { cat('...done.\n') } - if (debug==TRUE) + if (plot==TRUE) { plot_netassoc_matrix(finalmatrix_trimmed, xlab='Species',ylab='Species',colorRampPalette(c('red','white','blue'))(51),main="S.E.S. co-occurrence score for sp x sp") } @@ -217,36 +153,125 @@ makenetwork <- function(obs, nul, method='pearson', kappa=2, numnulls=1000, debu # convert to network representation if (verbose==TRUE) { cat('Building network...') } network_all <- graph.adjacency(finalmatrix_trimmed,mode='upper',weighted=T) + + network_pos <- delete.edges(network_all, E(network_all)[E(network_all)$weight < 0]) + network_neg <- delete.edges(network_all, E(network_all)[E(network_all)$weight > 0]) + E(network_neg)$weight <- -1*E(network_neg)$weight + if (verbose==TRUE) { cat('...done.\n') } - if (debug==TRUE) + if (plot==TRUE) { plot_netassoc_network(network_all) title('Association network') } - if (debug==TRUE) + if (plot==TRUE) { par(mfrow=c(1,1)) } - if (debug==TRUE) + + return(list( + matrix_spsite_obs=obs, + matrix_spsite_null=nul, + matrix_spsp_obs=fm_obs, + matrix_spsp_null_mean=fm_nul_mean, + matrix_spsp_null_sd=fm_nul_sd, + matrix_spsp_ses_all=finalmatrix, + matrix_spsp_ses_thresholded=finalmatrix_trimmed, + network_all=network_all, + network_pos=network_pos, + network_neg=network_neg + )) + +} + + + + +# plot network +plot_netassoc_network <- function(network, layout = layout.auto(network), + vertex.label = V(network)$name, + vertex.color = NA, + vertex.shape = "none", + vertex.label.color = "black", + vertex.label.family = "sans", + edge.width = NULL, + edge.color = NULL, + edge.arrow.size = 0.05, + vertex.label.cex = 0.5, + ...) +{ + if(is.null(edge.width)) { - return(list( - matrix_spsite_obs=obs, - matrix_spsite_null=nul, - matrix_spsite_nulls=nulls, - matrix_spsp_obs=fm_obs, - matrix_spsp_null_mean=fm_nul_mean, - matrix_spsp_null_sd=fm_nul_sd, - matrix_spsp_ses_all=finalmatrix, - matrix_spsp_ses_thresholded=finalmatrix_trimmed, - network_all=network_all - )) + if(length(E(network)$weight)==0) + { + edge.width=1 + } + else + { + edge.width=sqrt(abs(E(network)$weight)) + } + } + + if(is.null(edge.color)) + { + if(length(E(network)$weight)==0) + { + edge.color <- 'black' + zlmin <- -1 + zlmax <- 1 + } + else + { + edge.color <- ifelse(E(network)$weight > 0, rgb(0,0,1,0.8),rgb(1,0,0,0.8)) + zlmax <- max(abs(E(network)$weight),na.rm=T) + zlmin = -1*zlmax + } + } + + plot(network, + layout=layout, + vertex.label=vertex.label, + edge.color=edge.color, + edge.width=edge.width, + vertex.color=vertex.color, + vertex.label.color=vertex.label.color, + vertex.shape=vertex.shape, + edge.arrow.size=edge.arrow.size, + vertex.label.cex=vertex.label.cex, + vertex.label.family=vertex.label.family, + ...) + + colors=colorRampPalette(c('red','white','blue'))(51) + legend('topleft',adj=c(0,0),legend=format(c(zlmin,zlmin/2+zlmax/2,zlmax),digits=2),fill=c(colors[1],colors[ceiling(length(colors)/2)],colors[length(colors)]),bg='white',cex=0.5) +} + +plot_netassoc_matrix <- function(data, colors, ylab='Species',xlab='Site',onesided=FALSE,main="") +{ + zlmax <- max(abs(as.numeric(data)),na.rm=T) + if (is.infinite(zlmax)) + { + zlmax <- 1 + } + if (onesided==TRUE) + { + zlmin = 0 } else { - return(network_all) + zlmin = -1*zlmax } + + image(t(data),col=colors,axes=F,zlim=c(zlmin, zlmax),main=main) + + mtext(side=2,text=ylab,cex=0.5) + mtext(side=3,text=xlab,cex=0.5) + legend('topleft',adj=c(0,0),legend=format(c(zlmin,zlmin/2+zlmax/2,zlmax),digits=2),fill=c(colors[1],colors[ceiling(length(colors)/2)],colors[length(colors)]),bg='white',cex=0.5) + box() } + + + diff --git a/data/fia.rda b/data/fia.rda deleted file mode 100644 index 47832f6..0000000 Binary files a/data/fia.rda and /dev/null differ diff --git a/man/fia.Rd b/man/fia.Rd deleted file mode 100644 index aff2ca3..0000000 --- a/man/fia.Rd +++ /dev/null @@ -1,15 +0,0 @@ -\name{fia} -\alias{fia_obs} -\alias{fia_nul} -\docType{data} -\title{ -Demo input species x site matrices -} -\description{ -Data from referenced paper. -\code{fia_obs}{Observed data: Community matrix of abundance of tree species in Forest Inventory of America plots.} -\code{fia_nul}{Null data: Suitability score for tree species in locations of Forest Inventory of America plots, as predicted by MaxEnt regional models of species' distributions.} -} -\references{ -Morueta-Holme, N., Blonder, B., et al. in preparation. -} \ No newline at end of file diff --git a/man/makenetwork.Rd b/man/makenetwork.Rd index 2610d21..d65322e 100644 --- a/man/makenetwork.Rd +++ b/man/makenetwork.Rd @@ -18,8 +18,8 @@ The resulting network can be analyzed using functions from the \code{igraph} net } \usage{ makenetwork(obs, nul, - method = "pearson", kappa = 2, numnulls = 1000, - debug = FALSE, verbose = TRUE) + whichmethod = "pearson", kappa = 2, numnulls = 1000, + plot = FALSE, verbose = TRUE) } %- maybe also 'usage' for other objects documented here. \arguments{ @@ -29,7 +29,7 @@ A m x n community matrix describing the abundance of m species at n sites. Repre \item{nul}{ A m x n community matrix describing the abundance of m species at n sites. Represents the regional null expectation data. } - \item{method}{ + \item{whichmethod}{ Method used to calculate co-occurrence score between abundance vectors of species across sites. Can be either \code{'pearson'} or \code{'spearman'}, as per arguments in \code{\link{cor}} } \item{kappa}{ @@ -38,7 +38,7 @@ Threshold value determining the minimum absolute standardized effect size to be \item{numnulls}{ Number of resamples of the \code{nul} matrix used to assemble null communities. Larger values produce more accurate results. } - \item{debug}{ + \item{plot}{ If \code{TRUE}, plots all intermediate matrices calculated by the algorithm. Can be used to visualize input. } \item{verbose}{ @@ -49,18 +49,17 @@ If \code{TRUE}, prints status updates and progress bars during calculations. Full explanation of the algorithm is provided in the referenced paper. } \value{ -If \code{debug=FALSE}, an \code{igraph} object representing the association network, with edge weights between species corresponding to thresholded standardized effect sizes. - -If \code{debug=TRUE}, a list with the following components: +A list with the following components: \item{\code{matrix_spsite_obs}}{The trimmed \code{obs} matrix} \item{\code{matrix_spsite_null}}{The trimmed \code{nul} matrix} -\item{\code{matrix_spsite_nulls}}{The \code{numnulls} resampled null matrices} \item{\code{matrix_spsp_obs}}{The observed co-occurrence scores between species} \item{\code{matrix_spsp_null_mean}}{The mean null co-occurrence scores between species} \item{\code{matrix_spsp_null_sd}}{The s.d. null co-occurrence scores between species} \item{\code{matrix_spsp_ses_all}}{The raw standardized effect size co-occurrence scores between species} \item{\code{matrix_spsp_ses_thresholded}}{The thresholded standardized effect size co-occurrence scores between species} \item{\code{network_all}}{An \code{igraph} object representing the association network} +\item{\code{network_pos}}{An \code{igraph} object representing an association network including only positive associations} +\item{\code{network_pos}}{An \code{igraph} object representing an association network including only negative associations} } \note{ Warnings may be generated when calculating co-occurrence scores using the \code{cor} function in cases where a pair of species co-occur identically at all sites. These warnings can be ignored. @@ -70,14 +69,14 @@ Morueta-Holme, N., Blonder, B., et al. in preparation. } \examples{ -# example using New World forest plot data from referenced manuscript -data(fia) +# generate random data +set.seed(1) +nsp <- 10 +nsi <- 5 +m_obs <- floor(matrix(rgamma(nsp*nsi,shape=5),ncol=nsi,nrow=nsp)) +m_nul <- floor(matrix(rexp(nsp*nsi,rate=0.05),ncol=nsi,nrow=nsp)) -# run in debug mode for full output -# trim input to a few species to pass CRAN time checks -nsp <- 30 -nsites <- 50 -n <- makenetwork(fia_obs[1:nsp,1:nsites], fia_nul[1:nsp,1:nsites], numnulls=50, debug=TRUE) +n <- makenetwork(m_obs, m_nul, numnulls=50, plot=TRUE) # extract network n$network_all diff --git a/man/netassoc-package.Rd b/man/netassoc-package.Rd index 0b11a4b..1739402 100644 --- a/man/netassoc-package.Rd +++ b/man/netassoc-package.Rd @@ -8,15 +8,6 @@ Inference of species assocations from co-occurrence data \description{ Infers species associations from local and regional-scale co-occurrence data (species x site matrices). The resulting network can be analyzed using functions from the \code{igraph} network package. } -\details{ -\tabular{ll}{ -Package: \tab netassoc\cr -Type: \tab Package\cr -Version: \tab 0.4.1\cr -Date: \tab 2015-01-29\cr -License: \tab What license is it under?\cr -} -} \author{ Benjamin Blonder, Naia Morueta-Holme diff --git a/man/plot_netassoc_network.Rd b/man/plot_netassoc_network.Rd index f2b0c6f..4e64b9d 100644 --- a/man/plot_netassoc_network.Rd +++ b/man/plot_netassoc_network.Rd @@ -61,18 +61,18 @@ Other arguments to be passed to \code{\link{plot.igraph}}. } \examples{ -# example using New World forest plot data from referenced manuscript -data(fia) +# generate random data +set.seed(1) +nsp <- 10 +nsi <- 5 +m_obs <- floor(matrix(rgamma(nsp*nsi,shape=5),ncol=nsi,nrow=nsp)) +m_nul <- floor(matrix(rexp(nsp*nsi,rate=0.05),ncol=nsi,nrow=nsp)) -# run in debug mode for full output -# trim input to a few species to pass CRAN time checks -nsp <- 30 -nsites <- 50 -n <- makenetwork(fia_obs[1:nsp,1:nsites], fia_nul[1:nsp,1:nsites], numnulls=50) +n <- makenetwork(m_obs, m_nul, numnulls=50, plot=TRUE) # plot -plot_netassoc_network(n) +plot_netassoc_network(n$network_all) # plot using circular layout -plot_netassoc_network(n, layout=layout.circle(n)) +plot_netassoc_network(n$network_all, layout=layout.circle(n$network_all)) }