Skip to content
GitLab
Commit b49d4181 authored by Georges Kunstler's avatar Georges Kunstler
Browse files

change in mean species and genus sd

parent db0de601
master BRIEFCOM ctfs package test.dplyr
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 103 additions and 43 deletions
+103 -43
R/FUN.TRY.R
+ 8
− 0
View file @ b49d4181
......@@ -351,3 +351,11 @@ if(genus){
return(data.t)
}
}
## function to compute mean sd over species or genus with more than Nobs.lim observation
fun.sd.mean.sp.or.genus <- function(traits,species,genus=FALSE,Nobs.lim)
{
data.temp <- fun.sd.sp.or.genus((traits),species,genus=genus)
return(mean(data.temp[data.temp[,"nobs"]>Nobs.lim & !is.na(data.temp[,"sd"]),"sd"]))
}
TRY.R
+ 95
− 43
View file @ b49d4181
......@@ -133,12 +133,6 @@ species.tab2$Latin_name_syn<- trim.trailing(species.tab2$Latin_name_syn)
## tnrs(c("Quercus robur"), getpost="POST", source_ = "NCBI")
## tnrs(c("Quercus robur"), getpost="POST", source_ = "MSW3")
## http://taxosaurus.org/submit?query=Fagus+sylvatica
## {"status":"OK","names":[{"matchCount":2,"matches":[{"acceptedName":"","sourceId":"iPlant_TNRS","score":"1","matchedName":"Fagus sylvatica","annotations":{"Authority":""},"uri":""},{"acceptedName":"Fagus sylvatica","sourceId":"NCBI","score":"1","matchedName":"Fagus sylvatica","annotations":{},"uri":"http://www.ncbi.nlm.nih.gov/taxonomy/28930"}],"submittedName":"Fagus sylvatica"}],"metadata":{"spellcheckers":[{"name":"NCBI","description":"NCBI Spell Checker","annotations":{},"uri":"http://www.ncbi.nlm.nih.gov/","sourceId":1,"publication":"http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2578899/","call":"python2.6 tnrs_spellchecker/ncbi_spell.py","rank":1}],"sources":[{"status":"200: OK","name":"NCBI","description":"NCBI Taxonomy","uri":"http://www.ncbi.nlm.nih.gov/taxonomy","annotations":{},"sourceId":"NCBI","publication":"Federhen S. The Taxonomy Project.2002 Oct 9 [Updated 2003 Aug 13]. In: McEntyre J., Ostell J., editors. The NCBI Handbook [Internet]. Bethesda (MD): National Center for Biotechnology Information (US);2002.","rank":3,"code":"ICZN,ICN,ICNB"},{"status":"200: OK","name":"iPlant Collaborative TNRS v3.1","description":"The iPlant Collaborative TNRS provides parsing and fuzzy matching for plant taxa.","uri":"http://tnrs.iplantcollaborative.org/","annotations":{"Authority":"Author attributed to the accepted name (where applicable)."},"sourceId":"iPlant_TNRS","publication":"Boyle, B. et.al. The taxonomic name resolution service: an online tool for automated standardization of plant names. BMC Bioinformatics. 2013, 14:16. doi:10.1186/1471-2105-14-16. If you use TNRS results in a publication, please also cite The Taxonomic Name Resolution Service; http://tnrs.iplantcollaborative.org; version 3.1.","rank":2,"code":"ICN"},{"status":"200: OK","name":"Mammal Species of the World v3.0","description":"Mammal Species of the World, 3rd edition (MSW3) is a database of mammalian taxonomy. Our adaptor searches the indexed database for both exact and loose mathces","uri":"http://www.bucknell.edu/msw3/","annotations":{"Authority":"Don E. Wilson & DeeAnn M. Reeder (editors). 2005. Mammal Species of the World. A Taxonomic and Geographic Reference (3rd ed)"},"sourceId":"MSW3","publication":"Don E. Wilson & DeeAnn M. Reeder (editors). 2005. Mammal Species of the World. A Taxonomic and Geographic Reference (3rd ed)","rank":4,"code":"ICZN"}],"sub_date":"Thu Jul 18 16:36:23 2013","resolver_version":"1.2.0","jobId":"f7e6e1feba55b8f5d41a208630385630"}}
## http://taxosaurus.org/submit?query=Quercus+robur
## {"status":"OK","names":[{"matchCount":2,"matches":[{"acceptedName":"Quercus robur","sourceId":"iPlant_TNRS","score":"1","matchedName":"Quercus robur","annotations":{"Authority":"(Ten.) A. DC."},"uri":"http://www.tropicos.org/Name/50280607"},{"acceptedName":"Quercus robur","sourceId":"NCBI","score":"1","matchedName":"Quercus robur","annotations":{},"uri":"http://www.ncbi.nlm.nih.gov/taxonomy/38942"}],"submittedName":"Quercus robur"}],"metadata":{"spellcheckers":[{"name":"NCBI","description":"NCBI Spell Checker","annotations":{},"uri":"http://www.ncbi.nlm.nih.gov/","sourceId":1,"publication":"http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2578899/","call":"python2.6 tnrs_spellchecker/ncbi_spell.py","rank":1}],"sources":[{"status":"200: OK","name":"NCBI","description":"NCBI Taxonomy","uri":"http://www.ncbi.nlm.nih.gov/taxonomy","annotations":{},"sourceId":"NCBI","publication":"Federhen S. The Taxonomy Project.2002 Oct 9 [Updated 2003 Aug 13]. In: McEntyre J., Ostell J., editors. The NCBI Handbook [Internet]. Bethesda (MD): National Center for Biotechnology Information (US);2002.","rank":3,"code":"ICZN,ICN,ICNB"},{"status":"200: OK","name":"iPlant Collaborative TNRS v3.1","description":"The iPlant Collaborative TNRS provides parsing and fuzzy matching for plant taxa.","uri":"http://tnrs.iplantcollaborative.org/","annotations":{"Authority":"Author attributed to the accepted name (where applicable)."},"sourceId":"iPlant_TNRS","publication":"Boyle, B. et.al. The taxonomic name resolution service: an online tool for automated standardization of plant names. BMC Bioinformatics. 2013, 14:16. doi:10.1186/1471-2105-14-16. If you use TNRS results in a publication, please also cite The Taxonomic Name Resolution Service; http://tnrs.iplantcollaborative.org; version 3.1.","rank":2,"code":"ICN"},{"status":"200: OK","name":"Mammal Species of the World v3.0","description":"Mammal Species of the World, 3rd edition (MSW3) is a database of mammalian taxonomy. Our adaptor searches the indexed database for both exact and loose mathces","uri":"http://www.bucknell.edu/msw3/","annotations":{"Authority":"Don E. Wilson & DeeAnn M. Reeder (editors). 2005. Mammal Species of the World. A Taxonomic and Geographic Reference (3rd ed)"},"sourceId":"MSW3","publication":"Don E. Wilson & DeeAnn M. Reeder (editors). 2005. Mammal Species of the World. A Taxonomic and Geographic Reference (3rd ed)","rank":4,"code":"ICZN"}],"sub_date":"Thu Jul 18 16:44:45 2013","resolver_version":"1.2.0","jobId":"a00b4210bdfeb8c7e6d2f6bf7f10df04"}}
## lapply(unique(species.tab2$Latin_name_syn)[1:10],FUN=tnrs)
......@@ -236,46 +230,97 @@ saveRDS(data.ifn.species.try.noout ,file="./data/process/data.ifn.species.try.no
data.ifn.species.try.noout <- readRDS("./data/process/data.ifn.species.try.noout.rds")
#####################################################################
#### assume that the SD is equal mean species if less than 10 obs same for genus
#### assume that the SD is equal mean species
########
## USE TABLE 5 in Kattge et al. 2011
### LMA species sd log 0.09
## The table 5 in Kattge et al. 2011 GCB provides estimation of mean species sd
### SLA species sd log 0.09
### Nmass species sd log 0.08
### Seed Mass sd log 0.13
#### SPECIES mean sd
# SEE
## sd.log.SLA <- 0.09 ### based on Kattge et al. 2011
## sd.log.Nmass <- 0.08 ### based on Kattge et al. 2011
## sd.log.Seed.Mass <- 0.13 ### based on Kattge et al. 2011
## sd.log.LL <- 0.03 ### based on Kattge et al. 2011
######
# for wood density no value reportedin Kattge et al so need to compute mean sd for species withe more than 5 obs
data.sd.WD.log <- fun.sd.sp.or.genus(log10(TRY.DATA.FORMATED$StdValue.Stem.specific.density..SSD.),TRY.DATA.FORMATED$AccSpeciesName)
sd.log.WD <- mean(data.sd.WD.log[data.sd.WD.log[,2]>4 & !is.na(data.sd.WD.log[,2]),1])
### HOWEVER INCLUDING SPECIES WITH FEW observation ( less than 5 or 10) seems to lead to an underestimation
## see code below
sd.log.SLA <- 0.09 ### based on Kattge et al. 2011
sd.log.Nmass <- 0.08 ### based on Kattge et al. 2011
sd.log.Seed.Mass <- 0.13 ### based on Kattge et al. 2011
sd.log.LL <- 0.03 ### based on Kattge et al. 2011
### so I rpopose to use rather sd of species with nore than 5 obs
lm.obj <- (lm(log10(TRY.DATA.FORMATED$StdValue.Stem.specific.density..SSD.)~ TRY.DATA.FORMATED$AccSpeciesName))
sd.log.WD <-sd(residuals(lm.obj))
lm.obj <- (lm(log10(TRY.DATA.FORMATED$StdValue.Leaf.specific.area..SLA.)~ TRY.DATA.FORMATED$AccSpeciesName))
sd.log.SLA <-sd(residuals(lm.obj))
lm.obj <- (lm(log10(TRY.DATA.FORMATED$StdValue.Leaf.lifespan)~ TRY.DATA.FORMATED$AccSpeciesName))
sd.log.LL <-sd(residuals(lm.obj))
#######
##### COMPUTE GENUS MEAN SD IN THIS TRY DATA EXTRACTION AS NOT REPORTED IN Kattge et al. 2011
lm.obj <- (lm(log10(TRY.DATA.FORMATED$StdValue.Seed.mass)~ TRY.DATA.FORMATED$AccSpeciesName))
sd.log.Seed.Mass <-sd(residuals(lm.obj))
data.sd.WD.log.genus <- fun.sd.sp.or.genus(log10(TRY.DATA.FORMATED$StdValue.Stem.specific.density..SSD.),TRY.DATA.FORMATED$AccSpeciesName,genus=TRUE)
sd.log.WD.genus <- mean(data.sd.WD.log.genus[data.sd.WD.log.genus[,"nobs"]>10 & !is.na(data.sd.WD.log.genus[,"nobs"]),"sd"])
lm.obj <- (lm(log10(TRY.DATA.FORMATED$StdValue.Leaf.nitrogen..N..content.per.dry.mass)~ TRY.DATA.FORMATED$AccSpeciesName))
sd.log.Nmass <-sd(residuals(lm.obj))
data.sd.SLA.log.genus <- fun.sd.sp.or.genus(log10(TRY.DATA.FORMATED$StdValue.Leaf.specific.area..SLA.),TRY.DATA.FORMATED$AccSpeciesName,genus=TRUE)
sd.log.SLA.genus <- mean(data.sd.SLA.log.genus[data.sd.SLA.log.genus[,"nobs"]>10 & !is.na(data.sd.SLA.log.genus[,"nobs"]),"sd"])
##
Nobs.lim <- 5
sd.log.WD <- fun.sd.mean.sp.or.genus(log10(TRY.DATA.FORMATED$StdValue.Stem.specific.density..SSD.),
TRY.DATA.FORMATED$AccSpeciesName,genus=FALSE,Nobs.lim=Nobs.lim)
data.sd.LL.log.genus <- fun.sd.sp.or.genus(log10(TRY.DATA.FORMATED$StdValue.Leaf.lifespan),TRY.DATA.FORMATED$AccSpeciesName,genus=TRUE)
sd.log.LL.genus <- mean(data.sd.LL.log.genus[data.sd.LL.log.genus[,"nobs"]>10 & !is.na(data.sd.LL.log.genus[,"nobs"]),"sd"])
sd.log.SLA <- fun.sd.mean.sp.or.genus(log10(TRY.DATA.FORMATED$StdValue.Leaf.specific.area..SLA.),
TRY.DATA.FORMATED$AccSpeciesName,genus=FALSE,Nobs.lim=Nobs.lim)
sd.log.LL <- fun.sd.mean.sp.or.genus(log10(TRY.DATA.FORMATED$StdValue.Leaf.lifespan),
TRY.DATA.FORMATED$AccSpeciesName,genus=FALSE,Nobs.lim=Nobs.lim)
data.sd.Seed.Mass.log.genus <- fun.sd.sp.or.genus(log10(TRY.DATA.FORMATED$StdValue.Seed.mass),TRY.DATA.FORMATED$AccSpeciesName,genus=TRUE)
sd.log.Seed.Mass.genus <- mean(data.sd.Seed.Mass.log.genus[data.sd.Seed.Mass.log.genus[,"nobs"]>10 &
!is.na(data.sd.Seed.Mass.log.genus[,"nobs"]),"sd"])
sd.log.Seed.Mass <- fun.sd.mean.sp.or.genus(log10(TRY.DATA.FORMATED$StdValue.Seed.mass),
TRY.DATA.FORMATED$AccSpeciesName,genus=FALSE,Nobs.lim=Nobs.lim)
data.sd.Nmass.log.genus <- fun.sd.sp.or.genus(log10(TRY.DATA.FORMATED$StdValue.Leaf.nitrogen..N..content.per.dry.mass),
TRY.DATA.FORMATED$AccSpeciesName,genus=TRUE)
sd.log.Nmass.genus <- mean(data.sd.Nmass.log.genus[data.sd.Nmass.log.genus[,"nobs"]>10 & !is.na(data.sd.Nmass.log.genus[,"nobs"]),"sd"])
sd.log.Nmass <- fun.sd.mean.sp.or.genus(log10(TRY.DATA.FORMATED$StdValue.Leaf.nitrogen..N..content.per.dry.mass)
,TRY.DATA.FORMATED$AccSpeciesName,genus=FALSE,Nobs.lim=Nobs.lim)
## check how sd is related to minimum number of observation per species against simul
vec.sd.theo <- c()
vec.sd.theo2 <- c()
vec.sd
nsp <- table(data.sd.SLA.log[,"nobs"])[-(1:2)]
nobs <- as.numeric(names(nsp))
i <- 1
for ( nlow in 1:(length(nsp)-15))
{
vec.sd.theo[i] <- mean(unlist(sapply(nlow:length(nsp),FUN=function(i,nsp,nobs) replicate(nsp[i],expr=sd(rnorm(nobs[i],mean=0,sd=1))),nsp,nobs)))
vec.sd.theo2[i] <-sd(unlist(sapply(nlow:length(nsp),FUN=function(i,nsp,nobs) replicate(nsp[i],expr=(rnorm(nobs[i],mean=0,sd=1))),nsp,nobs)))
i <- i+1
}
plot(nobs[1:(length(nsp)-15)],vec.sd.theo,type="b",xlim=c(0,50))
lines(nobs[1:(length(nsp)-15)],vec.sd.theo2,type="b",xlim=c(0,50))
abline(h=1,col="red")
#######
##### COMPUTE GENUS MEAN SD IN THIS TRY DATA EXTRACTION AS NOT REPORTED IN Kattge et al. 2011 for all species with more than 10 OBSERVARION
Nobs.lim <- 5
sd.log.WD.genus <- fun.sd.mean.sp.or.genus(log10(TRY.DATA.FORMATED$StdValue.Stem.specific.density..SSD.),
TRY.DATA.FORMATED$AccSpeciesName,genus=TRUE,Nobs.lim=Nobs.lim)
sd.log.SLA.genus <- fun.sd.mean.sp.or.genus(log10(TRY.DATA.FORMATED$StdValue.Leaf.specific.area..SLA.),
TRY.DATA.FORMATED$AccSpeciesName,genus=TRUE,Nobs.lim=Nobs.lim)
sd.log.LL.genus <- fun.sd.mean.sp.or.genus(log10(TRY.DATA.FORMATED$StdValue.Leaf.lifespan),
TRY.DATA.FORMATED$AccSpeciesName,genus=TRUE,Nobs.lim=Nobs.lim)
sd.log.Seed.Mass.genus <- fun.sd.mean.sp.or.genus(log10(TRY.DATA.FORMATED$StdValue.Seed.mass),
TRY.DATA.FORMATED$AccSpeciesName,genus=TRUE,Nobs.lim=Nobs.lim)
sd.log.Nmass.genus <- fun.sd.mean.sp.or.genus(log10(TRY.DATA.FORMATED$StdValue.Leaf.nitrogen..N..content.per.dry.mass)
,TRY.DATA.FORMATED$AccSpeciesName,genus=TRUE,Nobs.lim=Nobs.lim)
#############
......@@ -288,14 +333,21 @@ genus.names <- paste(c("Leaf.N","Seed.mass","SLA","Wood.Density","Leaf.Lifespan"
,"genus",sep=".")
sd.vec.sp <- c(sd.log.Nmass,sd.log.Seed.Mass,sd.log.SLA,sd.log.WD,sd.log.LL)
names(sd.vec.sp) <- c("sdlog10.sp.Nmass","sdlog10.sp.Seed.Mass","sdlog10.sp.SLA","sdlog10.sp.WD","sdlog10.sp.LL")
sd.vec.genus <- c(sd.log.Nmass.genus,sd.log.Seed.Mass.genus,sd.log.SLA.genus,sd.log.WD.genus,sd.log.LL.genus)
names(sd.vec.genus) <- c("sdlog10.gs.Nmass","sdlog10.gs.Seed.Mass","sdlog10.gs.SLA","sdlog10.gs.WD","sdlog10.gs.LL")
## save mean species and genus sd
saveRDS(sd.vec.sp,file="./data/process/sd.vec.sp.rds")
saveRDS(sd.vec.genus,file="./data/process/sd.vec.genus.rds")
### function to select obs with less than Nthresh obs
fun.select.sd.with.to.few.obs.sp <- function(data,sd.names,nobs.names,genus.names,Nthreshold=2){
fun.select.sd.with.to.few.obs.sp <- function(data,sd.names,nobs.names,genus.names,Nthreshold=10){
(data[[nobs.names[i]]]<Nthreshold & !is.na(data[[nobs.names[i]]])
& !data[[genus.names[i]]]) }
fun.select.sd.with.to.few.obs.genus <- function(data,sd.names,nobs.names,genus.names,Nthreshold=2){
fun.select.sd.with.to.few.obs.genus <- function(data,sd.names,nobs.names,genus.names,Nthreshold=10){
(data[[nobs.names[i]]]<Nthreshold & !is.na(data[[nobs.names[i]]])
& data[[genus.names[i]]]) }
####
......@@ -303,13 +355,13 @@ data.TRY.sd.update <- data.ifn.species.try.noout
for (i in 1:length(nobs.names)){
print( sd.names[i])
print("species")
print(fun.select.sd.with.to.few.obs.sp(data=data.TRY.sd.update,sd.names,nobs.names,genus.names,Nthreshold=10))
data.TRY.sd.update[[sd.names[i]]][fun.select.sd.with.to.few.obs.sp(data=data.TRY.sd.update,sd.names,nobs.names,genus.names,Nthreshold=10)] <-
print(fun.select.sd.with.to.few.obs.sp(data=data.TRY.sd.update,sd.names,nobs.names,genus.names,Nthreshold=Nobs.lim))
data.TRY.sd.update[[sd.names[i]]][fun.select.sd.with.to.few.obs.sp(data=data.TRY.sd.update,sd.names,nobs.names,genus.names,Nthreshold=Nobs.lim)] <-
sd.vec.sp[i]
print("genus")
print(fun.select.sd.with.to.few.obs.genus(data=data.TRY.sd.update,sd.names,nobs.names,genus.names,Nthreshold=10))
data.TRY.sd.update[[sd.names[i]]][fun.select.sd.with.to.few.obs.genus(data=data.TRY.sd.update,sd.names,nobs.names,genus.names,Nthreshold=10)] <-
print(fun.select.sd.with.to.few.obs.genus(data=data.TRY.sd.update,sd.names,nobs.names,genus.names,Nthreshold=Nobs.lim))
data.TRY.sd.update[[sd.names[i]]][fun.select.sd.with.to.few.obs.genus(data=data.TRY.sd.update,sd.names,nobs.names,genus.names,Nthreshold=Nobs.lim)] <-
sd.vec.genus[i]
}
......@@ -318,13 +370,13 @@ saveRDS(data.TRY.sd.update,file="./data/process/data.TRY.sd.update.rds")
###
# plot sd to show mark
pdf("./figs/sd.traits.pdf")
r <- barplot(sd.vec.sp ,names.arg=c("Leaf.N","SM","SLA","WD","LL"),las=2,ylim=c(0,0.6),ylab="sd log10")
r <- barplot(sd.vec.sp ,names.arg=c("Leaf.N","SM","SLA","WD","LL"),las=2,ylim=c(0,0.9),ylab="sd log10")
points(r[,1],sd.vec.genus,col="red",pch=16,cex=2)
for (i in 1:length(nobs.names)){
sd.obs <- data.TRY.sd.update[[sd.names[i]]][!data.TRY.sd.update[[genus.names[i]]]]
points(rep(r[i,1],length(sd.obs)),sd.obs)
points(r[i,1],sd.vec.genus[i],col="red",pch=16,cex=2)
sd.obs <- data.TRY.sd.update[[sd.names[i]]][data.TRY.sd.update[[genus.names[i]]]]
points(rep(r[i,1],length(sd.obs)),sd.obs,col="red",pch=4)
## sd.obs <- data.TRY.sd.update[[sd.names[i]]][!data.TRY.sd.update[[genus.names[i]]]]
## points(rep(r[i,1],length(sd.obs)),sd.obs)
## sd.obs <- data.TRY.sd.update[[sd.names[i]]][data.TRY.sd.update[[genus.names[i]]]]
## points(rep(r[i,1],length(sd.obs)),sd.obs,col="red",pch=4)
print(sd.obs)
}
dev.off()
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment