Skip to content

Commit

Permalink
version 0.4.7
Browse files Browse the repository at this point in the history
  • Loading branch information
@bblonder @gaborcsardi
bblonder authored and gaborcsardi committed Jun 27, 2015
1 parent a878384 commit ffd2145
Show file tree
Hide file tree
Showing 5 changed files with 68 additions and 60 deletions.
8 changes: 4 additions & 4 deletions DESCRIPTION
View file
@@ -1,14 +1,14 @@
Package: netassoc
Type: Package
Title: Inference of Species Associations from Co-Occurrence Data
Version: 0.4.3
Date: 2015-03-13
Version: 0.4.7
Date: 2015-06-27
Author: Benjamin Blonder, Naia Morueta-Holme
Maintainer: Benjamin Blonder <bblonder@gmail.com>
Description: Infers species associations from local and regional-scale co-occurrence data (species x site matrices).
License: GPL-3
Depends: igraph
Packaged: 2015-03-13 17:01:47 UTC; benjaminblonder
Packaged: 2015-06-24 15:19:46 UTC; benjaminblonder
NeedsCompilation: no
Repository: CRAN
Date/Publication: 2015-03-13 20:22:09
Date/Publication: 2015-06-24 18:23:01
8 changes: 4 additions & 4 deletions MD5
View file
@@ -1,6 +1,6 @@
648150b10eab05c75734ec845e3afd27 *DESCRIPTION
72e45f04c3eb03fe020deb7a1419a699 *DESCRIPTION
11ce28fad68ae62136749e9a93424939 *NAMESPACE
c3806531e8b9051a4f8d142cf8d7184c *R/netassoc.R
9bc661c5f6f4054499e3872f585cd795 *man/makenetwork.Rd
15109ce762c0e73fdce2af7854de8986 *R/netassoc.R
ad3b4b828cffb83bc1bad5b07f10d861 *man/makenetwork.Rd
6aa2c25359d4b486497bfa5ce19e4520 *man/netassoc-package.Rd
491d858df7cc04641612582e169b9d42 *man/plot_netassoc_network.Rd
c4f15a61333a62adf788d4cf74d8232e *man/plot_netassoc_network.Rd
97 changes: 52 additions & 45 deletions R/netassoc.R
View file
@@ -1,7 +1,7 @@
# 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)
{
nul_resample <- matrix(0, nrow=nrow(nul),ncol=ncol(nul),dimnames=dimnames(nul))
nul_resample <- matrix(0, nrow=nrow(nul),ncol=ncol(nul),dimnames=dimnames(nul)) # auto-set sparseness depending on filling
for (i in 1:ncol(nul_resample))
{
# fill in resamples with the same total abundance
Expand All @@ -12,7 +12,7 @@ generate_nul_resample <- function(nul, obs)
return(nul_resample)
}

makenetwork <- function(obs, nul, whichmethod='pearson', kappa=2, numnulls=1000, plot=FALSE,verbose=TRUE)
makenetwork <- function(obs, nul, whichmethod='pearson', alpha=0.05, numnulls=1000, plot=FALSE,verbose=TRUE)
{
if (is.matrix(obs) & is.null(dimnames(obs)))
{
Expand Down Expand Up @@ -42,6 +42,13 @@ makenetwork <- function(obs, nul, whichmethod='pearson', kappa=2, numnulls=1000,
{
stop(sprintf("Some sites are empty: %s",paste(names(problemsites),collapse=", ")))
}

numtests <- nrow(obs)*(nrow(obs)-1)/2
alpha_bonferroni <- (alpha/numtests)

kappa <- qnorm(1-alpha_bonferroni/2)

cat(sprintf("Number of tests: %d\nBonferroni-corrected alpha value: %f\nKappa threshold for overall alpha=%f level: %f\n", numtests, alpha_bonferroni, alpha, kappa))

if (plot==TRUE)
{
Expand All @@ -58,26 +65,17 @@ makenetwork <- function(obs, nul, whichmethod='pearson', kappa=2, numnulls=1000,
plot_netassoc_matrix(nul, colors=colorRampPalette(c('white','black'))(51),onesided=TRUE,main="Null sp x site")
}


if (verbose==TRUE) { cat('Generating null replicates...') }
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 (plot==TRUE)
{
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))

finalmatrix <- matrix(NA,nrow=nrow(obs),ncol=nrow(obs))

if (verbose==TRUE) { cat('Calculating co-occurrence scores...') }


if (verbose==TRUE) { cat('Calculating observed co-occurrence scores...') }

count <- 0
for (i in 1:nrow(obs))
{
Expand All @@ -86,48 +84,57 @@ makenetwork <- function(obs, nul, whichmethod='pearson', kappa=2, numnulls=1000,
if (i!=j)
{
count <- count + 1

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)
suppressWarnings(
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 <- array(data=NA,dim=c(nrow(obs),nrow(obs),numnulls))
for (k in 1:numnulls)
{
# Generating null
if (verbose==TRUE) { cat(sprintf('Generating null replicate %d...', k)) }
thisnull <- generate_nul_resample(nul, obs)
if (verbose==TRUE) { cat('...done.\n') }

count <- 0
for (i in 1:nrow(obs))
{
for (j in i:nrow(obs))
{
cor_nul <- rep(NA, numnulls)

for (k in 1:numnulls)
if (i!=j)
{
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)
count <- count + 1
if (verbose==TRUE) { cat (sprintf('%d %d %.3f\n', i, j, count/(nrow(obs)*(nrow(obs)-1)/2))) }

veci_nul <- as.numeric(thisnull[i,])
vecj_nul <- as.numeric(thisnull[j,])
suppressWarnings(
cor_nul[i,j,k] <- cor(x=veci_nul, y=vecj_nul, method=whichmethod)
)
}
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
if (verbose==TRUE) { cat('Calculating standardized effect sizes...') }
# calculating null distributions
fm_nul_mean <- apply(cor_nul, c(1,2), mean, na.rm=T)
fm_nul_sd <- apply(cor_nul, c(1,2), sd, na.rm=T)

finalmatrix <- (fm_obs - fm_nul_mean) / fm_nul_sd
if (verbose==TRUE) { cat('...done.\n') }


if (verbose==TRUE) { cat('...done.\n') }

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")
Expand Down Expand Up @@ -191,7 +198,7 @@ makenetwork <- function(obs, nul, whichmethod='pearson', kappa=2, numnulls=1000,


# plot network
plot_netassoc_network <- function(network, layout = layout.auto(network),
plot_netassoc_network <- function(network, layout = layout.fruchterman.reingold(network),
vertex.label = V(network)$name,
vertex.color = NA,
vertex.shape = "none",
Expand Down
13 changes: 7 additions & 6 deletions man/makenetwork.Rd
View file
Expand Up @@ -18,7 +18,7 @@ The resulting network can be analyzed using functions from the \code{igraph} net
}
\usage{
makenetwork(obs, nul,
whichmethod = "pearson", kappa = 2, numnulls = 1000,
whichmethod = "pearson", alpha = 0.05, numnulls = 1000,
plot = FALSE, verbose = TRUE)
}
%- maybe also 'usage' for other objects documented here.
Expand All @@ -32,8 +32,8 @@ A m x n community matrix describing the abundance of m species at n sites. Repre
\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}{
Threshold value determining the minimum absolute standardized effect size to be retained in the network. Set to zero to retain all associations.
\item{alpha}{
Analysis-wide Type I error rate. Bonferroni correction for multiple species-species comparisons and a normal distribution cumulative distribution function are used to determine a kappa threshold value determining the minimum absolute standardized effect size to be retained in the network. Set alpha to one to retain all assocations.
}
\item{numnulls}{
Number of resamples of the \code{nul} matrix used to assemble null communities. Larger values produce more accurate results.
Expand Down Expand Up @@ -72,11 +72,12 @@ Morueta-Holme, N., Blonder, B., et al. in preparation.
# generate random data
set.seed(1)
nsp <- 10
nsi <- 5
nsi <- 50
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))
m_nul <- 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))
n <- makenetwork(m_obs, m_nul, numnulls=50, plot=TRUE)
n <- makenetwork(m_obs, m_nul, numnulls=100, plot=TRUE)
# extract network
n$network_all
Expand Down
2 changes: 1 addition & 1 deletion man/plot_netassoc_network.Rd
View file
Expand Up @@ -8,7 +8,7 @@ Plots species association network
Draws a network of species associations. By default edge widths are proportional to association strength and edge color reflects assocation type (blue, positive; red, negative).
}
\usage{
plot_netassoc_network(network, layout = layout.auto(network),
plot_netassoc_network(network, layout = layout.fruchterman.reingold(network),
vertex.label = V(network)$name,
vertex.color = NA,
vertex.shape = "none",
Expand Down

0 comments on commit ffd2145

Please sign in to comment.