CODEWITHSUNDEEP
rplotvegan
Pedro Rocha
Pedro Rocha

Reputation: 79

How to add CI whiskers to a rarefaction curve

I need to make some rarefaction curves and I wish them to display the whiskers at the edges of the confidence interval bars, whereas the default is to simply display the bars without whiskers:

library(vegan) 
data("dune")
result <- specaccum(dune)
plot(result, lwd=2)

Default rarefaction curve Default rarefaction curve

I've tried to add some whiskers using the arrows function, but the results from the specaccum function only include the standard deviation. So I ended up with half the job done:

samples <- result$sites
error <- result$sd
richness <- result$richness
arrows(samples, richness-error, samples, richness+error, angle=90, code=3, length=0.05)

Rarefaction curve with sd whiskers Rarefaction curve with sd whiskers

From what I've searched, the most common approach would be to convert the confidence intervals into a shaded area (by using the argument ci.type="polygon") and then add a boxplot to the plotted curve. However, this leads to a very busy image that I'd rather avoid.

Does anyone have a more elegant solution?

Upvotes: 1

Views: 1586

Answers (3)

Jari Oksanen
Jari Oksanen

Reputation: 3682

The vegan repository in github has now an option for drawing those short horizontal bars to the error bars with argument ci.length. The default is zero (no horizontal bars) to preserve the old behaviour.

Upvotes: 0

Gavin Simpson
Gavin Simpson

Reputation: 174778

You forgot the multiplier (see argument ci in ?plot.specaccum). What you drew were for a ~68% confidence interval. Multiplying by 2 (ci = 2) gives an approximate 95% confidence interval, which is what plot.specaccum draws by default.

Including the (default) multiplier in a modification of the code you used

plot(result)
with(result, arrows(sites, richness - (2 * sd), sites, richness + (2 * sd),
                    angle = 90, code = 3, length = 0.05))

we get:

enter image description here

You can ignore the warning; the standard error of the last data point drawn is zero

> result$sd
 [1] 2.3510636 1.8763851 1.5722711 1.4469584 1.3901594 1.3530349 1.3164796
 [8] 1.2749034 1.2282010 1.1763410 1.1193437 1.0564537 0.9874094 0.9115998
[15] 0.8286890 0.7380921 0.6333903 0.5139710 0.3570714 0.0000000

and arrow() is just warning you that it won't draw a length 0 arrow.

Upvotes: 0

Ulluchu
Ulluchu

Reputation: 41

you could try to add a plotCI-plot:

library(vegan) 
library(plotrix) 
data("dune")
result <- specaccum(dune)

plot(result)
plotCI(result$sites,result$richness,result$sd*2,err="y", lwd=2,add=TRUE, pch=NA)

I admit this is also not the most elegant option as well but it works.

Upvotes: 1

Related Questions