R语言跑程序到处结果 崩了
R语言崩了
NTPD.nocons.fun=function(seedling,pdmatrix,sample,nsim=999,phylo)
{
sp.seedling=as.character(seedling$sp2)
focalsp=sp.seedling[match(colnames(sample),seedling$tag)]
MINPD= MINPD.nocons.fun (focalsp=focalsp,pdmatrix=pdmatrix, sample=sample)
MINPD.null=MINPD.nocons.null.fun(focalsp=focalsp,sample=sample,nsim=nsim,phylo=phylo)
NTPD =numeric()
for(i in 1:length(focalsp))
{
MINPD.null.mean=mean(MINPD.null[i,])
MINPD.null.sd=sd(MINPD.null[i,])
NTPD[i]=(MINPD[i]-MINPD.null.mean)/MINPD.null.sd
}
return(NTPD)
}
MINPD.nocons.fun =function(focalsp,pdmatrix,sample)
{
MINPD.nocons=numeric()
for(i in 1:length(focalsp))
{
spinsample=rownames(sample)[sample[,i]!=0]
spinsample.abund=sample[,i][sample[,i]!=0]
focalsp.pos=which(rownames(pdmatrix)==focalsp[i])
spinsample.pos=match(spinsample,colnames(pdmatrix))
pdtofocal=numeric()
for(j in 1:length(spinsample.pos))
{
pdtofocal[j]=pdmatrix[focalsp.pos,spinsample.pos[j]]
}
h<-which(spinsample==focalsp[i])
if (length(h)>0) {
pdtofocal.nocons=pdtofocal[-h]
MINPD.nocons[i]=min(pdtofocal.nocons)}
else
MINPD.nocons[i]=min(pdtofocal)#包括同种邻体的话，有同种邻体出现时，这里肯定是0
}
return(MINPD.nocons)
}
#MINPD.nocons.null.fn函数：计算在999次零模型条件下的MINPD
MINPD.nocons.null.fun=function(focalsp,sample,nsim,phylo)#phylo为系统发育树
{
require(picante)
MINPD.null=matrix(nrow=length(focalsp),ncol=nsim)
for(j in 1:nsim)
{
phylo.tem <- tipShuffle(phylo)
pdmatrix.tem=cophenetic(phylo.tem)
#pdmatrix.tem=pdmatrix.tem/(-max(pdmatrix.tem))+1
MINPD.null[,j]=MINPD.nocons.fun(focalsp=focalsp,pdmatrix=pdmatrix.tem,sample=sample)
}
return(MINPD.null)
}
###############################################
##THESE TOP TWO FUNCTIONS ARE NEEDED TO MAKE ##
######THE CIRCLE SAMPLE FILES FOR APd calculation ####
####COPY AND PASTE THESE TWO FUNCTIONS INTO R##
###############################################
###############################################
###############################################
###############################################
###############################################
#####################
#seedling: just including tree and shrub seedlings，得到每个幼苗个体同一个样方内的所有幼苗邻体的sample，得到一个矩阵，rownames:出现在所有sample里的幼苗邻体种名sp,colnames: focaldata的个体的tag
######################
sample.seneigh.allfocal.fun=function(seedling, radius,plotdim=c(200,200),minheight,include.cons,phylo)
#seedling data.frame (year 1), including tag (unique for each seedling),sp, gx,gy,height#
{
seedling.focal=seedling
colname.output=seedling.focal$tag
seedling.neighbour=seedling.focal #tree.neighbour is data.frame
focaldata=as.matrix(seedling.focal)#class(focaldata[1,1]):"character"
#apply sample.seneigh.onefocal.fun to each row in focaldata (each row is an seedling individual)
output=apply(focaldata,1,sample.seneigh.onefocal.fun,neighbour=seedling.neighbour,include.cons)
abund=table(seedling.neighbour$sp2)
rownames(output)=names(abund)
colnames(output)=colname.output
return(output)#output is matrix
}
#####################################################################
# sample.seneigh.onefocal.fun: function to produce a sample matrix with 1 row, called in apply above
#####################################################################
sample.seneigh.onefocal.fun=function(fdata,neighbour,include.cons)#fdata为一个目标幼苗个体,neighbour is all focaldata
{
quadrat=as.numeric(fdata["qua"])
tag.focal=as.character(fdata["tag"])
sp.focal=as.character(fdata["sp"])
#subset neighbour around focal individual in the same seedling plot
neighbour$sp<-as.factor(neighbour$sp)
square=subset(neighbour,qua==quadrat)
#remove focal individual
square=square[-which(square$tag==tag.focal),]
if (include.cons=="TRUE")
{
species.in.circle=table(square$sp)
sam=matrix("NA",nrow=1,ncol=length(species.in.circle))
if (nrow(square)>0) #if there are no neighbors, it will return NA's#
{
sam <- t(as.matrix(species.in.circle))
}
}
if (include.cons=="FALSE")
{
square=subset(square,sp!=sp.focal)
species.in.circle=table(square$sp)
sam=matrix("NA",nrow=1,ncol=length(species.in.circle))
if (nrow(square)>0) #if there are no neighbors, it will return NA's#
{
sam <- t(as.matrix(species.in.circle))
}
}
return(sam)#sam is matrix
}
#APdobs_NTPD.fn函数：计算NTPD观察值,长度为focalsp的长度
APdobs_NTPD.fn =function(focalsp,pdmatrix,sample)
#focalsp为目标个体的种名，pdmatrix为系统发育距离矩阵(标准化成0-1,参考Liza的课件)，spinsample为周围邻体的其他种类
{
pdtofocal.min =numeric()
for(i in 1:length(focalsp))
{
spinsample=rownames(sample)[sample[,i]!=0]
spinsample.abund=sample[,i][sample[,i]!=0]
focalsp.pos=which(rownames(pdmatrix)==focalsp[i])
spinsample.pos=match(spinsample,colnames(pdmatrix))
pdtofocal=numeric()
for(j in 1:length(spinsample.pos))
{
pdtofocal[j]=pdmatrix[focalsp.pos,spinsample.pos[j]]
}
pdtofocal.min[i]=min(pdtofocal)
}
return(pdtofocal.min)
}
##APdobs.fn函数：计算null model 平均值
APdobs.fun=function(focalsp,pdmatrix,sample)#focalsp为目标个体的种名，pdmatrix为系统发育距离矩阵(标准化成0-1,参考Liza的课件)，spinsample为周围邻体的其他种类
{
pdtofocal.mean.all=numeric()
for(i in 1:length(focalsp))
{
spinsample=rownames(sample)[sample[,i]!=0]
spinsample.abund=sample[,i][sample[,i]!=0]
focalsp.pos=which(rownames(pdmatrix)==focalsp[i])
spinsample.pos=match(spinsample,colnames(pdmatrix))
pdtofocal=numeric()
for(j in 1:length(spinsample.pos))
{
pdtofocal[j]=pdmatrix[focalsp.pos,spinsample.pos[j]]*spinsample.abund[j]
}
h<-which(names(spinsample.abund)==focalsp[i])
if (length(h)>0){
focal.abund=spinsample.abund[h]
pdtofocal.mean=sum(pdtofocal)/(sum(spinsample.abund)-focal.abund) }
if(length(h)==0){
pdtofocal.mean=sum(pdtofocal)/sum(spinsample.abund) }
pdtofocal.mean.all[i]=pdtofocal.mean
}
return(pdtofocal.mean.all)
}
#APdnull.fn函数：计算NTPD在999次零模型条件下的NTPD
APdnull.fun=function(focalsp,sample,nsim,phylo)#phylo为系统发育树
{
require(picante)
pdtofocal.mean.all.null=matrix(nrow=length(focalsp),ncol=nsim)
for(j in 1:nsim)
{
phylo.tem <- tipShuffle(phylo)
pdmatrix.tem=cophenetic(phylo.tem)
pdmatrix.tem=pdmatrix.tem/(-max(pdmatrix.tem))+1
pdtofocal.mean.all=APdobs.fun(focalsp=focalsp,pdmatrix=pdmatrix.tem,sample=sample)
pdtofocal.mean.all.null[,j]=pdtofocal.mean.all
}
return(pdtofocal.mean.all.null)
}
APd.allfocalsp.fun=function(focalsp=focalsp,pdmatrix=pdmatrix,sample,nsim=999,phylo=phylo.xtbg)
{
pdtofocal.min.all= APdobs_NTPD.fn (focalsp=focalsp,pdmatrix=pdmatrix, sample=sample)
pdtofocal.mean.null.all=APdnull.fun(focalsp=focalsp,sample=sample,nsim=nsim,phylo=phylo)
NTPD =numeric()
for(i in 1:length(focalsp))
{
APdnull.mean=mean(pdtofocal.mean.null.all[i,])
APdnull.sd=sd(pdtofocal.mean.null.all[i,])
NTPD[i]=(pdtofocal.min.all[i]-APdnull.mean)/APdnull.sd
}
return(NTPD)
}
############steps#########################
library(picante)
setwd("C:/Users/Administrator/Desktop/幼苗_1")
phylo<-read.tree("ym_mytree.txt")
d1<-read.csv("1.csv",header=T,fileEncoding="GBK")##
d2=read.csv('NCI_1.csv',header=T,fileEncoding="GBK")
d2$sp2=NA
d2$sp2=d1[match(d2$tag,d1$tag),]$sp
seedling.focal=d2[is.element(d2$sp2,phylo$tip.label),]#只保留系统发育树上有的种类
phylo$tip.label=d2[match(d2$sp2,phylo$tip.label),]$sp
seedling.focal$sp=NA
x=read.csv('x.csv',header=T,fileEncoding="GBK")
seedling.focal$sp=x[match(seedling.focal$tag,x$tag),]$sp
nrow(seedling.focal)
seedling.focal$gx=NA
seedling.focal$gy=NA
qua=read.csv('qua_location.csv',header=T,fileEncoding="GBK")
seedling.focal$gx=qua[match(seedling.focal$qua,qua$qua),]$gx
seedling.focal$gy=qua[match(seedling.focal$qua,qua$qua),]$gy
pdmatrix<-read.csv("phy-dis.csv",header=T,fileEncoding="GBK")#系统发育距离矩阵
#sp.seedling=seedling.focal$sp
#length(sp.seedling)
#focalsp=sp.seedling
rownames(pdmatrix)=pdmatrix$X
pdmatrix=pdmatrix[,-1]
#colnames(pdmatrix)=rownames(pdmatrix)
sample=sample.seneigh.allfocal.fun(seedling.focal, radius=3,plotdim=c(200,200),minheight=10,include.cons=T,phylo) ###用于计算幼苗邻体指数
results <-NTPD.nocons.fun(seedling=seedling.focal,pdmatrix=pdmatrix,sample,nsim=999,phylo=phylo)
seedling.focal$s_ntpd=NA
seedling.focal$s_ntpd=results
which(is.na(seedling.focal$s_ntpd))
write.table(seedling.focal, "ntpd.seedling.csv")
在结果导出时崩啦 怎么解决呀
- 这篇文章：R语言捕获错误 也许有你想要的答案，你可以看看
如果你已经解决了该问题, 非常希望你能够分享一下解决方案, 写成博客, 将相关链接放在评论区, 以帮助更多的人 ^-^