Hi Nick, Yes, you are right. There's one small bug on my code. The 'if' within the for-loop is wrong. Try it now with the code below.
rrarefy.custom <- function (x, sample, rep.param=F) { if (!identical(all.equal(x, round(x)), TRUE)) stop("function is meaningful only for integers (counts)") x <- as.matrix(x) if (ncol(x) == 1) x <- t(x) if (length(sample) > 1 && length(sample) != nrow(x)) stop(gettextf("length of 'sample' and number of rows of 'x' do not match")) sample <- rep(sample, length = nrow(x)) colnames(x) <- colnames(x, do.NULL = FALSE) nm <- colnames(x) if (!rep.param && any(rowSums(x) < sample)) warning("Some row sums < 'sample' and are not rarefied") for (i in 1:nrow(x)) { if (!rep.param && sum(x[i, ]) <= sample[i]) next row <- sample(rep(nm, times = x[i, ]), sample[i], replace = rep.param) row <- table(row) ind <- names(row) x[i, ] <- 0 x[i, ind] <- row } x } I have test it now before pasting it. Best, *Luisfo Chiroque* /PhD Student | PhD Candidate IMDEA Networks Institute/ http://fourier.networks.imdea.org/people/~luis_nunez/ <http://fourier.networks.imdea.org/people/%7Eluis_nunez/> On 09/07/2016 10:27 PM, Nick Pardikes wrote: > Hey Luisfo, > > This looks great, however I still get the same plot as before (seen > below). The output looks the same. Here is the figure that was > generated from this code: > > rrarefy.custom <- function (x, sample, rep.param=F) > { > if (!identical(all.equal(x, round(x)), TRUE)) > stop("function is meaningful only for integers (counts)") > x <- as.matrix(x) > if (ncol(x) == 1) > x <- t(x) > if (length(sample) > 1 && length(sample) != nrow(x)) > stop(gettextf("length of 'sample' and number of rows of 'x' do not > match")) > sample <- rep(sample, length = nrow(x)) > colnames(x) <- colnames(x, do.NULL = FALSE) > nm <- colnames(x) > if (!rep.param && any(rowSums(x) < sample)) > warning("Some row sums < 'sample' and are not rarefied") > for (i in 1:nrow(x)) { > if (rep.param && sum(x[i, ]) <= sample[i]) > next > row <- sample(rep(nm, times = x[i, ]), sample[i], replace = rep.param) > row <- table(row) > ind <- names(row) > x[i, ] <- 0 > x[i, ind] <- row > } > x > } > > raredata <- rarecurve(rrarefy.custom(netdata, sample=100,rep.param=T), > label=F, col=rgb(0, 0, 1, 0.1)) > > However, I like what you did to the rrarefy function to add the sample > with replacement option. > > On Wed, Sep 7, 2016 at 8:17 AM, Luisfo <luisf...@yahoo.es> wrote: >> Hi Nick, >> >> If you use the following >> raredata <- rarecurve(rrarefy(netdata, sample=100), label=F, col=rgb(0, >> 0, 1, 0.1)) >> should work for any sample size, e.g. sample=100. >> However, you will have a 'warning' if you don't have samples enough, because >> it has not replacement. >> >> If you type 'rrarefy' on the R console (without brackets), or any other >> function name, you will see the R code of the function. >> rrarefy uses the function 'sample()' for sampling, but has no option for >> replacement. >> I did the following. I created my custom rrarefy function from the original. >> rrarefy.custom <- function (x, sample, rep.param=F) >> { >> if (!identical(all.equal(x, round(x)), TRUE)) >> stop("function is meaningful only for integers (counts)") >> x <- as.matrix(x) >> if (ncol(x) == 1) >> x <- t(x) >> if (length(sample) > 1 && length(sample) != nrow(x)) >> stop(gettextf("length of 'sample' and number of rows of 'x' do not >> match")) >> sample <- rep(sample, length = nrow(x)) >> colnames(x) <- colnames(x, do.NULL = FALSE) >> nm <- colnames(x) >> if (!rep.param && any(rowSums(x) < sample)) >> warning("Some row sums < 'sample' and are not rarefied") >> for (i in 1:nrow(x)) { >> if (rep.param && sum(x[i, ]) <= sample[i]) >> next >> row <- sample(rep(nm, times = x[i, ]), sample[i], replace = rep.param) >> row <- table(row) >> ind <- names(row) >> x[i, ] <- 0 >> x[i, ind] <- row >> } >> x >> } >> You can check the differences with the original code if you type 'rrarefy' >> on the R console. >> >> So now, if you type the following >> raredata <- rarecurve(rrarefy.custom(netdata, sample=100,rep.param=T), >> label=F, col=rgb(0, 0, 1, 0.1)) >> you will have the desired behaviour. >> >> WARNING: I do not understand about rarefunction curves or communities in >> your context. So, be careful when resampling. It might not be statistically >> correct. >> >> Regards, >> Luisfo Chiroque >> PhD Student | PhD Candidate >> IMDEA Networks Institute >> http://fourier.networks.imdea.org/people/~luis_nunez/ >> >> On 09/07/2016 12:07 AM, Nick Pardikes wrote: >> >> I am currently having difficulty producing a graph using rarecurve in the >> vegan package. I have produced rarefaction curves (seen below) using the >> following code. >> >> >> library(vegan) >> >> myMat <- round(matrix(rlnorm(2000), 50)) #creates distribution of >> communities >> >> netdata <- as.data.frame(myMat) #generates a matrix of communities (rows), >> species (columns) >> >> raredata <- rarecurve(netdata, label=F, col=rgb(0, 0, 1, 0.1)) #uses >> rarecurve to plot a rarefaction for each individual community (n=50) >> >> >> However I would like to produce a graph in which all rarefaction curves end >> at the same sample size. For example, in this graph it would be great to >> extend the x-axis (sample size) to 100 and have all curves end at this >> point. Is there any way to use rarecurve to resample a community (row) with >> replacement the same number of times for all 50 communities? With >> replacement is important because the communities differ greatly in their >> size (number of species). >> >> >> I understand that rarefaction is useful to compare communities with >> different sample efforts, but I would still like to generate the figure. My >> actual data has 5000 simulated communities that differ greatly in matrix >> size and number of samples. >> >> >> Thank you in advance for your help and suggestions. >> >> >> Cheers, >> >> Nick >> >> >> >> ______________________________________________ >> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see >> https://stat.ethz.ch/mailman/listinfo/r-help >> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html >> and provide commented, minimal, self-contained, reproducible code. >> >> > > [[alternative HTML version deleted]] ______________________________________________ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.