[Esm-commits] r39 - in pkg/ESM: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Sep 28 16:00:08 CEST 2010
Author: timotheepoisot
Date: 2010-09-28 16:00:08 +0200 (Tue, 28 Sep 2010)
New Revision: 39
Modified:
pkg/ESM/DESCRIPTION
pkg/ESM/R/.Rhistory
pkg/ESM/R/functions.r
pkg/ESM/man/getspe.Rd
Log:
update
Modified: pkg/ESM/DESCRIPTION
===================================================================
--- pkg/ESM/DESCRIPTION 2010-09-21 10:01:01 UTC (rev 38)
+++ pkg/ESM/DESCRIPTION 2010-09-28 14:00:08 UTC (rev 39)
@@ -1,5 +1,5 @@
Package: ESM
-Version: 1.3.0-01
+Version: 1.4.0-01
Date: 2010-08-29
Title: Ecological Specificity Measures
Author: Timothee Poisot <tpoisot at um2.fr>
Modified: pkg/ESM/R/.Rhistory
===================================================================
--- pkg/ESM/R/.Rhistory 2010-09-21 10:01:01 UTC (rev 38)
+++ pkg/ESM/R/.Rhistory 2010-09-28 14:00:08 UTC (rev 39)
@@ -1,133 +1,27 @@
+load('DynEvolNet.Rdata')#
+NET<-t(empty(NETS[['7C5']]))
+setwd('/Users/Tim/Desktop/Data et tests/DynEvolNet')
+library(ESM)#
+load('DynEvolNet.Rdata')#
+NET<-t(empty(NETS[['7C5']]))
+NET
+image(Overlap(NET))
+getOLsp(Overlap(NET),10)
+getOLvec
+Overlap(NET)
+nrow(Overlap(NET))
+getOLvec(NET)
+length(getOLvec(NET))
+dim(NET)
+OLmat <- Overlap(net)
+OLmat <- Overlap(NET)
+l <- c(1:nrow(OLmat))
+l
+OL <- lapply(l,function(x)getOLsp(OLmat,id=x))
+OL
+setwd('/Users/Tim/esm/pkg/ESM/R')
empty = function(mat) mat[rowSums(mat)>0,colSums(mat)>0]#
#
-get.PDI <- function(fit)#
-{#
- if(length(fit)==1){return(1)}#
- fit <- sort(as.vector(fit),decreasing=TRUE)#
- test <- fit[2:length(fit)]#
- out <- sum(fit[1]-test)/length(test)#
- return(out)#
-}#
-#
-get.RR <- function(pop)#
-{#
- if(length(pop)==1){return(1)}#
- Ni <- sum(pop>0)#
- N <- length(pop)#
- return(1-(Ni-1)/(N-1))#
-}#
-#
-get.HS <- function(pop,ifzero=1e-12)#
-{#
- if(length(pop)==1){return(1)}#
- partiel <- NULL#
- pop <- as.vector(pop)#
- pop[pop==0] <- ifzero#
- for(i in 1:length(pop))#
- {#
- p <- pop[i]/sum(pop)#
- partiel[i] <- p * log(p)#
- }#
- shannon <- ifelse(length(pop)<=1,0,-sum(partiel) / log(length(pop)))#
- return(1-shannon)#
-}#
-#
-get.SSI <- function(occup)#
-{#
- if(length(occup)==1){return(1)}#
- score <- NULL#
- fit <- as.vector(occup)#
- h <- sum(fit>0)#
- H <- length(occup)#
- SSI <- sqrt((H/h-1)/(H-1))#
- return(SSI)#
-}#
-#
-pdi <- function(m)#
-{#
- spe <- vector('numeric',length=nrow(m))#
- for(i in 1:nrow(m))#
- {#
- spe[i] <- get.PDI(m[i,])#
- }#
- return(spe)#
-}#
-#
-hs <- function(m)#
-{#
- spe <- vector('numeric',length=nrow(m))#
- for(i in 1:nrow(m))#
- {#
- spe[i] <- get.HS(as.vector(as.numeric(m[i,])))#
- }#
- return(spe)#
-}#
-#
-ssi <- function(m)#
-{#
- spe <- vector('numeric',length=nrow(m))#
- for(i in 1:nrow(m))#
- {#
- spe[i] <- get.SSI(as.vector(as.numeric(m[i,])))#
- }#
- return(spe)#
-}#
-#
-rr <- function(m)#
-{#
- spe <- vector('numeric',length=nrow(m))#
- for(i in 1:nrow(m))#
- {#
- spe[i] <- get.RR(as.vector(as.numeric(m[i,])))#
- }#
- return(spe)#
-}#
-#
-getspe <- function(mat,measure=pdi,...)#
-{#
- if(max(mat)!=1){mat<-mat/max(mat)}#
- out <- measure(as.matrix(mat),...)#
- names(out) <- rownames(mat)#
- return(out)#
-}#
-#
-cv = function(d) sd(d)/mean(d)#
-last = function(d) d[length(d)]#
-#
-scale = function(v,m=0,M=1)#
-{#
- v <- v-min(v)#
- v <- v/max(v)#
- v <- v*(M-m)#
- v <- v+m#
- return(v)#
-}#
-#
-dmat = function(m,n=2)#
-{#
- e <- (max(m)-min(m))/n#
- cl <- seq(from=0,to=1,by=e)#
- nm <- matrix(0,ncol=ncol(m),nrow=nrow(m))#
- for(i in 1:(length(cl)-1))#
- {#
- nm[(cl[i]<m)&(m<=cl[(i+1)])] <- cl[i]#
- }#
- return(scale(nm,0,1))#
-}
-emptt()
-empty()
-M <- matrix(round(runif(1e6),0),1e3)
-empty(mat)
-M <- matrix(round(rbinom(1e6,1,0.35),0),1e3)
-dim(M)
-M
-M <- matrix(round(rbinom(1e6,1,0.05),0),1e3)
-dim(empty(M))
-empty(M)
-colSums(M)
-M <- matrix(round(rbinom(1e6,1,0.25),0),1e3)
-empty = function(mat) mat[rowSums(mat)>0,colSums(mat)>0]#
-#
pdi <- function(fit)#
{#
if(length(fit)==1){return(1)}#
@@ -137,287 +31,14 @@
return(out)#
}#
#
-get.RR <- function(pop)#
+rr <- function(fit)#
{#
- if(length(pop)==1){return(1)}#
- Ni <- sum(pop>0)#
- N <- length(pop)#
+ N <- length(fit)#
+ if(N==1){return(1)}#
+ Ni <- sum(fit>0)#
return(1-(Ni-1)/(N-1))#
}#
#
-get.HS <- function(pop,ifzero=1e-12)#
-{#
- if(length(pop)==1){return(1)}#
- partiel <- NULL#
- pop <- as.vector(pop)#
- pop[pop==0] <- ifzero#
- for(i in 1:length(pop))#
- {#
- p <- pop[i]/sum(pop)#
- partiel[i] <- p * log(p)#
- }#
- shannon <- ifelse(length(pop)<=1,0,-sum(partiel) / log(length(pop)))#
- return(1-shannon)#
-}#
-#
-get.SSI <- function(occup)#
-{#
- if(length(occup)==1){return(1)}#
- score <- NULL#
- fit <- as.vector(occup)#
- h <- sum(fit>0)#
- H <- length(occup)#
- SSI <- sqrt((H/h-1)/(H-1))#
- return(SSI)#
-}#
-#
-hs <- function(m)#
-{#
- spe <- vector('numeric',length=nrow(m))#
- for(i in 1:nrow(m))#
- {#
- spe[i] <- get.HS(as.vector(as.numeric(m[i,])))#
- }#
- return(spe)#
-}#
-#
-ssi <- function(m)#
-{#
- spe <- vector('numeric',length=nrow(m))#
- for(i in 1:nrow(m))#
- {#
- spe[i] <- get.SSI(as.vector(as.numeric(m[i,])))#
- }#
- return(spe)#
-}#
-#
-rr <- function(m)#
-{#
- spe <- vector('numeric',length=nrow(m))#
- for(i in 1:nrow(m))#
- {#
- spe[i] <- get.RR(as.vector(as.numeric(m[i,])))#
- }#
- return(spe)#
-}#
-#
-getspe <- function(mat,measure=pdi,...)#
-{#
- mat <- empty(mat)#
- mat <- mat/max(mat)#
- out <- unlist(apply(mat,1,measure,...))#
- names(out) <- rownames(mat)#
- return(out)#
-}#
-#
-cv = function(d) sd(d)/mean(d)#
-last = function(d) d[length(d)]#
-#
-scale = function(v,m=0,M=1)#
-{#
- v <- v-min(v)#
- v <- v/max(v)#
- v <- v*(M-m)#
- v <- v+m#
- return(v)#
-}#
-#
-dmat = function(m,n=2)#
-{#
- e <- (max(m)-min(m))/n#
- cl <- seq(from=0,to=1,by=e)#
- nm <- matrix(0,ncol=ncol(m),nrow=nrow(m))#
- for(i in 1:(length(cl)-1))#
- {#
- nm[(cl[i]<m)&(m<=cl[(i+1)])] <- cl[i]#
- }#
- return(scale(nm,0,1))#
-}
-empty = function(mat) mat[rowSums(mat)>0,colSums(mat)>0]#
-#
-pdi <- function(fit)#
-{#
- if(length(fit)==1){return(1)}#
- fit <- sort(as.vector(fit),decreasing=TRUE)#
- test <- fit[2:length(fit)]#
- out <- sum(fit[1]-test)/length(test)#
- return(out)#
-}#
-#
-rr <- function(pop)#
-{#
- if(length(pop)==1){return(1)}#
- Ni <- sum(pop>0)#
- N <- length(pop)#
- return(1-(Ni-1)/(N-1))#
-}#
-#
-hs <- function(pop,ifzero=1e-12)#
-{#
- if(length(pop)==1){return(1)}#
- partiel <- NULL#
- pop <- as.vector(pop)#
- pop[pop==0] <- ifzero#
- for(i in 1:length(pop))#
- {#
- p <- pop[i]/sum(pop)#
- partiel[i] <- p * log(p)#
- }#
- shannon <- ifelse(length(pop)<=1,0,-sum(partiel) / log(length(pop)))#
- return(1-shannon)#
-}#
-#
-ssi <- function(occup)#
-{#
- if(length(occup)==1){return(1)}#
- score <- NULL#
- fit <- as.vector(occup)#
- h <- sum(fit>0)#
- H <- length(occup)#
- SSI <- sqrt((H/h-1)/(H-1))#
- return(SSI)#
-}#
-#
-getspe <- function(mat,measure=pdi,...)#
-{#
- mat <- empty(mat)#
- mat <- mat/max(mat)#
- out <- unlist(apply(mat,1,measure,...))#
- names(out) <- rownames(mat)#
- return(out)#
-}#
-#
-cv = function(d) sd(d)/mean(d)#
-last = function(d) d[length(d)]#
-#
-scale = function(v,m=0,M=1)#
-{#
- v <- v-min(v)#
- v <- v/max(v)#
- v <- v*(M-m)#
- v <- v+m#
- return(v)#
-}#
-#
-dmat = function(m,n=2)#
-{#
- e <- (max(m)-min(m))/n#
- cl <- seq(from=0,to=1,by=e)#
- nm <- matrix(0,ncol=ncol(m),nrow=nrow(m))#
- for(i in 1:(length(cl)-1))#
- {#
- nm[(cl[i]<m)&(m<=cl[(i+1)])] <- cl[i]#
- }#
- return(scale(nm,0,1))#
-}
-getspe(M,ssi)
-system.time(getspe(M))
-system.time(getspe(M,rr))
-system.time(getspe(M,ssi))
-plot(getspe(M,pdi),getspe(M,rr))
-plot(getspe(M,pdi),getspe(M,ssi))
-cv(getspe(M))
-cv(getspe(M,rr))
-cv(getspe(M,ssi))
-cv(getspe(M,hs))
-empty = function(mat) mat[rowSums(mat)>0,colSums(mat)>0]#
-#
-pdi <- function(fit)#
-{#
- if(length(fit)==1){return(1)}#
- fit <- sort(as.vector(fit),decreasing=TRUE)#
- test <- fit[2:length(fit)]#
- out <- sum(fit[1]-test)/length(test)#
- return(out)#
-}#
-#
-rr <- function(pop)#
-{#
- if(length(pop)==1){return(1)}#
- Ni <- sum(pop>0)#
- N <- length(pop)#
- return(1-(Ni-1)/(N-1))#
-}#
-#
-hs <- function(pop,ifzero=1e-12)#
-{#
- if(length(pop)==1){return(1)}#
-# partiel <- NULL#
- pop <- as.vector(pop)#
- pop[pop==0] <- ifzero#
- pp <- pop/sum(pop)#
- partiel <- pp * log(pp)#
-# for(i in 1:length(pop))#
-# {#
-# p <- pop[i]/sum(pop)#
-# partiel[i] <- p * log(p)#
-# }#
- shannon <- ifelse(length(pop)<=1,0,-sum(partiel) / log(length(pop)))#
- return(1-shannon)#
-}#
-#
-ssi <- function(occup)#
-{#
- if(length(occup)==1){return(1)}#
- score <- NULL#
- fit <- as.vector(occup)#
- h <- sum(fit>0)#
- H <- length(occup)#
- SSI <- sqrt((H/h-1)/(H-1))#
- return(SSI)#
-}#
-#
-getspe <- function(mat,measure=pdi,...)#
-{#
- mat <- empty(mat)#
- mat <- mat/max(mat)#
- out <- unlist(apply(mat,1,measure,...))#
- names(out) <- rownames(mat)#
- return(out)#
-}#
-#
-cv = function(d) sd(d)/mean(d)#
-last = function(d) d[length(d)]#
-#
-scale = function(v,m=0,M=1)#
-{#
- v <- v-min(v)#
- v <- v/max(v)#
- v <- v*(M-m)#
- v <- v+m#
- return(v)#
-}#
-#
-dmat = function(m,n=2)#
-{#
- e <- (max(m)-min(m))/n#
- cl <- seq(from=0,to=1,by=e)#
- nm <- matrix(0,ncol=ncol(m),nrow=nrow(m))#
- for(i in 1:(length(cl)-1))#
- {#
- nm[(cl[i]<m)&(m<=cl[(i+1)])] <- cl[i]#
- }#
- return(scale(nm,0,1))#
-}
-empty = function(mat) mat[rowSums(mat)>0,colSums(mat)>0]#
-#
-pdi <- function(fit)#
-{#
- if(length(fit)==1){return(1)}#
- fit <- sort(as.vector(fit),decreasing=TRUE)#
- test <- fit[2:length(fit)]#
- out <- sum(fit[1]-test)/length(test)#
- return(out)#
-}#
-#
-rr <- function(pop)#
-{#
- if(length(pop)==1){return(1)}#
- Ni <- sum(pop>0)#
- N <- length(pop)#
- return(1-(Ni-1)/(N-1))#
-}#
-#
hs <- function(fit,ifzero=1e-12)#
{#
if(length(fit)==1){return(1)}#
@@ -440,10 +61,11 @@
return(SSI)#
}#
#
-getspe <- function(mat,measure=pdi,...)#
+getspe <- function(mat,measure=pdi,normal='whole',...)#
{#
mat <- empty(mat)#
- mat <- mat/max(mat)#
+ if(normal=='species'){mat <- t(apply(mat,1,function(x)x/max(x)))}#
+ if(normal=='whole'){mat <- mat/max(mat)}#
out <- unlist(apply(mat,1,measure,...))#
names(out) <- rownames(mat)#
return(out)#
@@ -471,37 +93,43 @@
nm[(cl[i]<m)&(m<=cl[(i+1)])] <- cl[i]#
}#
return(scale(nm,0,1))#
-}
-system.time(getspe(M,hs))
-rownames(M) <- c(1:1000)
-getspe(M,pdi)
-boxplot(getspe(M,pdi))
-hist(getspe(M,pdi))
-hist(getspe(M,rr))
-hist(getspe(M,ssi))
-hist(getspe(M,hs))
-setwd('/Users/Tim/esm/pkg/ESM/R')
-nullweb = function(ref)#
+}#
+#
+Overlap = function(mat)#
{#
- ref <- empty(ref)#
- ref[ref>0] <- 1#
- margin <- ifelse(ncol(ref)<nrow(ref),2,1)#
- currentweb <- matrix(0,ncol=ncol(ref),nrow=nrow(ref))#
- pc <- colMeans(ref)#
- pr <- rowMeans(ref)#
- if(margin==2)#
+ for(i in 1:nrow(mat))#
{#
- for(i in 1:ncol(ref))#
+ mat[i,] <- mat[i,]/sum(mat[i,])#
+ }#
+ out <- matrix(0,nrow=nrow(mat),ncol=nrow(mat))#
+ colnames(out) <- rownames(mat)#
+ rownames(out) <- rownames(mat)#
+ for(i in 1:(nrow(mat)-1))#
+ {#
+ for(j in (i+1):nrow(mat))#
{#
- currentweb[,i] <- (pc[i]+pr)/2#
+ out[i,j] <- 1-1/2*sum(abs(mat[i,]-mat[j,]))#
}#
- } else {#
- for(i in 1:nrow(ref))#
- {#
- currentweb[i,] <- (pr[i]+pc)/2#
- }#
}#
- return(apply(currentweb,margin,function(x)rbinom(length(x),1,x)))#
+ diag(out) <- rep(1,length(diag(out)))#
+ return(out)#
}#
#
-nullreps = function(ref,rep) replicate(rep,nullweb(ref),FALSE)
+getOLsp = function(OLmat,id)#
+{#
+ byR <- OLmat[id,]#
+ byC <- OLmat[,id]#
+ OlSP <- c(byR[c((id+1):length(byR))],byC[c(1:(id-1))])#
+ return(OlSP)#
+}#
+#
+getOLvec = function(net)#
+{#
+ net <- empty(net)#
+ net <- net/max(net)#
+ OLmat <- Overlap(net)#
+ l <- c(1:nrow(OLmat))#
+ OL <- lapply(l,function(x)getOLsp(OLmat,id=x))#
+ if(!is.null(rownames(OLmat))){names(OL)<-rownames(OLmat)}#
+ return(OL)#
+}
Modified: pkg/ESM/R/functions.r
===================================================================
--- pkg/ESM/R/functions.r 2010-09-21 10:01:01 UTC (rev 38)
+++ pkg/ESM/R/functions.r 2010-09-28 14:00:08 UTC (rev 39)
@@ -11,9 +11,9 @@
rr <- function(fit)
{
- if(length(fit)==1){return(1)}
+ N <- length(fit)
+ if(N==1){return(1)}
Ni <- sum(fit>0)
- N <- length(fit)
return(1-(Ni-1)/(N-1))
}
@@ -39,10 +39,11 @@
return(SSI)
}
-getspe <- function(mat,measure=pdi,...)
+getspe <- function(mat,measure=pdi,normal='whole',...)
{
mat <- empty(mat)
- mat <- mat/max(mat)
+ if(normal=='species'){mat <- t(apply(mat,1,function(x)x/max(x)))}
+ if(normal=='whole'){mat <- mat/max(mat)}
out <- unlist(apply(mat,1,measure,...))
names(out) <- rownames(mat)
return(out)
@@ -70,4 +71,43 @@
nm[(cl[i]<m)&(m<=cl[(i+1)])] <- cl[i]
}
return(scale(nm,0,1))
+}
+
+Overlap = function(mat)
+{
+ for(i in 1:nrow(mat))
+ {
+ mat[i,] <- mat[i,]/sum(mat[i,])
+ }
+ out <- matrix(0,nrow=nrow(mat),ncol=nrow(mat))
+ colnames(out) <- rownames(mat)
+ rownames(out) <- rownames(mat)
+ for(i in 1:(nrow(mat)-1))
+ {
+ for(j in (i+1):nrow(mat))
+ {
+ out[i,j] <- 1-1/2*sum(abs(mat[i,]-mat[j,]))
+ }
+ }
+ diag(out) <- rep(1,length(diag(out)))
+ return(out)
+}
+
+getOLsp = function(OLmat,id)
+{
+ byR <- OLmat[id,]
+ byC <- OLmat[,id]
+ OlSP <- c(byR[c((id+1):length(byR))],byC[c(1:(id-1))])
+ return(OlSP)
+}
+
+getOLvec = function(net)
+{
+ net <- empty(net)
+ net <- net/max(net)
+ OLmat <- Overlap(net)
+ l <- c(1:nrow(OLmat))
+ OL <- lapply(l,function(x)getOLsp(OLmat,id=x))
+ if(!is.null(rownames(OLmat))){names(OL)<-rownames(OLmat)}
+ return(OL)
}
\ No newline at end of file
Modified: pkg/ESM/man/getspe.Rd
===================================================================
--- pkg/ESM/man/getspe.Rd 2010-09-21 10:01:01 UTC (rev 38)
+++ pkg/ESM/man/getspe.Rd 2010-09-28 14:00:08 UTC (rev 39)
@@ -4,21 +4,25 @@
\title{Specificity}
\description{Returns a vector containing the specificity of each organism}
\usage{
-getspe(mat,measure=pdi,...)
+getspe(mat,measure=pdi,normal='whole',...)
}
\arguments{
\item{mat}{A matrix giving the performance of each organism on each resource, with organisms as rows}
\item{measure}{The index of specificity to use. See details, defaults to \kbd{pdi}.}
+\item{normal}{The method to normalize the matrix. Defaults to \kbd{whole}, can be set to \kbd{species}.}
\item{...}{Additional arguments to be passed to \kbd{measure}}
}
\value{
- \item{}{A vector with the specificity of each organism. Uses names if the matrix as such property.}
+ \item{}{A vector with the specificity of each organism. Uses names if the matrix has such property.}
}
-\details{Values of \kbd{measure} can be : \kbd{rr}, the classical resource range; \kbd{ssi}, the species specialization index; \kbd{pdi}, the Paired Differences Index; \kbd{hs}, Shannon's entropy.}
+\details{Values of \kbd{measure} can be : \kbd{rr}, the classical resource range; \kbd{ssi}, the species specialization index; \kbd{pdi} (default), the Paired Differences Index; \kbd{hs}, Shannon's entropy.
+\kbd{normal} is the way to normalize the performance value within the matrix. Defaults to \kbd{whole}, which means that the matrix is divided by its maximal value. Can be set to \kbd{species}, so that the maximal value of each species is set to 1. The performances of the measures were not evaluated with this method of normalization, and it will result in a overestimation of the specificity.
+}
+
\references{For PDI, and this package, cite : Poisot, T, Canard E, Mouquet N & Hochberg, M E "..." XXX.
For HS, cite : Schug J, et al. (2005) Promoter features related to tissue specificity as measured by Shannon entropy. Genome biology 6(4):R33
More information about the Esm-commits
mailing list