| 所在主题: | |
| 文件名: community.txt | |
| 资料下载链接地址: https://bbs.pinggu.org/a-1224433.html | |
| 附件大小: | |
|
各位大虾,本人在用r语言做论文,想求通过ses命令求得sd等值,但运行下列命令后,不知道如何看到具体结果,求助啊????????{:soso_e127:}
命令如下 `ses.pd` <- function (samp, tree, null.model = c("taxa.labels", "richness", "frequency", "sample.pool", "phylogeny.pool", "independentswap", "trialswap"), runs = 999, iterations = 1000, ...) { pd.obs <- as.vector(pd(samp, tree, ...)$PD) null.model <- match.arg(null.model) pd.rand <- switch(null.model, taxa.labels = t(replicate(runs, as.vector(pd(samp, tipShuffle(tree), ...)$PD))), richness = t(replicate(runs, as.vector(pd(randomizeMatrix(samp,null.model="richness"), tree, ...)$PD))), frequency = t(replicate(runs, as.vector(pd(randomizeMatrix(samp,null.model="frequency"), tree, ...)$PD))), sample.pool = t(replicate(runs, as.vector(pd(randomizeMatrix(samp,null.model="richness"), tree, ...)$PD))), phylogeny.pool = t(replicate(runs, as.vector(pd(randomizeMatrix(samp,null.model="richness"), tipShuffle(tree), ...)$PD))), independentswap = t(replicate(runs, as.vector(pd(randomizeMatrix(samp,null.model="independentswap", iterations), tree, ...)$PD))), trialswap = t(replicate(runs, as.vector(pd(randomizeMatrix(samp,null.model="trialswap", iterations), tree, ...)$PD))) ) pd.rand.mean <- apply(X = pd.rand, MARGIN = 2, FUN = mean, na.rm=TRUE) pd.rand.sd <- apply(X = pd.rand, MARGIN = 2, FUN = sd, na.rm=TRUE) pd.obs.z <- (pd.obs - pd.rand.mean)/pd.rand.sd pd.obs.rank <- apply(X = rbind(pd.obs, pd.rand), MARGIN = 2, FUN = rank)[1, ] pd.obs.rank <- ifelse(is.na(pd.rand.mean),NA,pd.obs.rank) data.frame(ntaxa=specnumber(samp),pd.obs, pd.rand.mean, pd.rand.sd, pd.obs.rank, pd.obs.z, pd.obs.p=pd.obs.rank/(runs+1),runs=runs, row.names = row.names(samp)) } 前期还有一部分命令,主要是想得到pd的sd,大虾也可以参考一下 samp<-read.csv("d:\\127comm110.csv",row.names=1) phy<-read.tree("d:\\127tree.new") traits<-read.table("d:\\127traits.csv",header=TRUE,row.names=1,sep=",") plot(phy) plot(x=phy,type="phylogram",show.node.label=TRUE, cex=.75) prunedphy<-prune.sample(samp,phy) samp<-samp[,prunedphy$tip.label] traits<-traits[prunedphy$tip.label,] par(mfrow=c(2,3),cex=.75) #side by side plots for(i in row.names(samp)) { plot(prunedphy,show.tip.label=TRUE,main=i) tiplabels(tip=which(samp[i, ]>0), pch=19,col="blue",cex=.7) } par(mfrow=(c(1,1))) par(mfrow=c(1,1),cex=1) #just one plot plot(prunedphy, show.tip.label=TRUE, show.node.label=TRUE, cex=.7) tiplabels(tip=which(traits$habit=="tree"),pch=19,cex=1,col="blue") tiplabels(tip=which(traits$habit=="shrub"),pch=19,cex=1,col="yellow") tiplabels(tip=which(traits$habit=="vine"),pch=19,cex=1,col="red") plot(prunedphy, show.tip.label=TRUE, tip.color=as.numeric(traits$habit),edge.width=3) phydist<-as.data.frame(cophenetic(prunedphy)) species<-row.names(phydist) p<-cbind(traits,species,phydist) pd.result<-pd(samp,prunedphy,include.root=TRUE) #ses.mpd(samp,phydist,null.model="taxa.labels",abundance.weighted=FALSE,runs=999) #ses.mntd(samp,phydist,null.model="taxa.labels",abundance.weighted=FALSE,runs=999) psv.result<-psv(samp,prunedphy,compute.var=TRUE) psv.result pd.result cdnt<-comdistnt(samp,phydist) cdnt library(cluster) cclust<-hclust(cdnt) plot(cclust) pd.result<-pd(samp,prunedphy,include.root=TRUE,compute.var=TRUE) psv.result<-psv(samp,prunedphy,compute.var=TRUE) |
|
熟悉论坛请点击新手指南
|
|
| 下载说明 | |
|
1、论坛支持迅雷和网际快车等p2p多线程软件下载,请在上面选择下载通道单击右健下载即可。 2、论坛会定期自动批量更新下载地址,所以请不要浪费时间盗链论坛资源,盗链地址会很快失效。 3、本站为非盈利性质的学术交流网站,鼓励和保护原创作品,拒绝未经版权人许可的上传行为。本站如接到版权人发出的合格侵权通知,将积极的采取必要措施;同时,本站也将在技术手段和能力范围内,履行版权保护的注意义务。 (如有侵权,欢迎举报) |
|
京ICP备16021002号-2 京B2-20170662号
京公网安备 11010802022788号
论坛法律顾问:王进律师
知识产权保护声明
免责及隐私声明