Dear R-guru's,
I've been working on a script that performs a rarefaction of variance on
a 41 eigenvectors and their respective principal coordinates. The principal
coordinates are derived from a PCoA performed on a discrete morphological
character matrix. The rarefaction is meant to analyze changing variance
as taxa (n=60) are added, to determine if variance estimates between
subdivisions of the 60 taxon data set are the result of differences in
included taxa, or a product of different sized samples. I've been pounding
away at this script since Spring, and got a different version assessing
range
to work relatively easily.
The problem I have run into with the "Rarefaction Variance" script
is that I can get it to work correctly analyzing values by axes (columns;
eigenvectors),
but not rows (taxa), which is what I need. I have tried transposing the
data set, placing
the taxa in the columns and eigenvectors in the rows, to no avail (so far).
Here is the script:
##Rarefaction Variance ##
multi_vectors<-read.table("multi_vectors.txt",row.names=1, header=T)
mv<-multi_vectors[,c(1:41)]
attach(mv)
###Importing the eigenvector output from the PCoA analysis###
rarevar=data.frame("ntaxa"=1:60)
z=data.frame("sumvariances"=matrix(nrow=60, ncol=41))
z1=data.frame("sumvariances"=matrix(nrow=41, ncol=1))
z2=data.frame("sumvariances"=matrix(nrow=60, ncol=41))
for (l in 1:41)
{
for (j in 1:60)
{
z[j,l]=0
z1[j,1]=0
z2[j,l]=0
}
}
###Creating empty matrices to deposit outputs of the loops###
mvt=t(mv)
###Transposing the "mv" matrix so that variance of the taxa can be
assessed, ###
###instead of the eigenvectors.###
for (j in 2:60)
{
s1<-tmv[,sample(60,j,replace=FALSE)]
###this portion randomly samples the columns of taxa, starting with 2 and
continuing up to 60##
for (k in 1:41)
{
clm=s1[k,j]
v=var(clm)
z[j,k]=sum(v)
}
}
###This portion is where everything falls apart. As taxa are added,
variance should ###
###increase to a point, then level off. I have yet to achieve this, even
with this latest###
###of many iterations.###
###After variance values are determined within each iteration of "j",###
###they are to be stored in one of the "z" matrices I built at the
beginning of this script.###
###Additionally, the script would iteratively sum the variance values
across 100-1000 replications ###
###of the rarefaction code. After doing so, the mean of this sum would be
used to build a rarefaction curve.###
Below I have included the functioning Range Rarefaction script, to give
some idea of what I am going for.
Thank you for any input or guidance you can provide.
- David Levering
##Rarefaction: Range##
multi_vectors<-read.table("multi_vectors.txt",row.names=1, header=T)
mv<-multi_vectors[,c(1:41)]
attach(mv)
rarerange=data.frame("ntaxa"=1:60)
z=data.frame("ranges"=matrix(nrow=60, ncol=41))
z1=data.frame("sumranges"=matrix(nrow=60, ncol=1))
z2=data.frame("sumranges"=matrix(nrow=60, ncol=41))
##plot(30,110)##
for (l in 1:41)
{
for (j in 1:60)
{
z[j,l]=0
z1[j,1]=0
z2[j,l]=0
}
}
for (l in 1:100)
{
for (j in 2:60)
{
s1<-mv[sample(60,j,replace=FALSE),]
for (k in 1:41)
{
clm=s1[,k]
m=max(clm)
n=min(clm)
range=m-n
z[j,k]=range
print(z)
}
z1[j,]=(sum(z[j,]))/41
print(z1)
}
}
z=cbind(z,z1)
plot(z[,42])
write.csv(Z1[l],file="range_rarefaction2.csv",row.names=T)
--
David Levering
Department of Zoology
Life Science West
Oklahoma State University
Stillwater, OK 74078
[[alternative HTML version deleted]]
_______________________________________________
R-sig-phylo mailing list
[email protected]
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
