搜索
人大经济论坛 附件下载

附件下载

所在主题:
文件名:  community.txt
资料下载链接地址: https://bbs.pinggu.org/a-1224433.html
附件大小:
264 Bytes   举报本内容
各位大虾,本人在用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、本站为非盈利性质的学术交流网站,鼓励和保护原创作品,拒绝未经版权人许可的上传行为。本站如接到版权人发出的合格侵权通知,将积极的采取必要措施;同时,本站也将在技术手段和能力范围内,履行版权保护的注意义务。
(如有侵权,欢迎举报)
二维码

扫码加我 拉你入群

请注明:姓名-公司-职位

以便审核进群资格,未注明则拒绝

GMT+8, 2025-12-29 04:38