### R code from vignette source 'sads_intro.Rnw' ################################################### ### code chunk number 1: R setup ################################################### options(width=60, continue=" ") ################################################### ### code chunk number 2: installation (eval = FALSE) ################################################### ## install.packages('sads') ################################################### ### code chunk number 3: load-sads ################################################### library(sads) ################################################### ### code chunk number 4: install-dev-version (eval = FALSE) ################################################### ## library(devtools) ## install_github(repo = 'piLaboratory/sads', ref= 'dev', build_vignettes = TRUE) ################################################### ### code chunk number 5: load-sads-package ################################################### library(sads) ################################################### ### code chunk number 6: Loading datasets ################################################### data(moths)# William's moth data used by Fisher et al (1943) data(ARN82.eB.apr77)# Arntz et al. benthos data data(birds)# Bird census used by Preston (1948) ################################################### ### code chunk number 7: Tabulating species in octaves ################################################### (moths.oc <- octav(moths)) (arn.oc <- octav(ARN82.eB.apr77)) ################################################### ### code chunk number 8: Ploting-octaves ################################################### plot(moths.oc) ################################################### ### code chunk number 9: Biomass-octave-plot ################################################### plot(arn.oc) ################################################### ### code chunk number 10: octaves-relative-frequencies ################################################### plot(moths.oc, prop = TRUE, border=NA, col=NA) lines(octav(birds), mid = FALSE, prop = TRUE, col="red") lines(octav(moths), mid = FALSE, prop = TRUE) legend("topright", c("Preston's birds", "Fisher's moths"), col=c("red", "blue"), lty=1, bty="n") ################################################### ### code chunk number 11: Rank-abundance tables ################################################### head(moths.rad <- rad(moths)) head(arn.rad <- rad(ARN82.eB.apr77)) ################################################### ### code chunk number 12: radplot1 ################################################### plot(moths.rad, ylab="Number of individuals") ################################################### ### code chunk number 13: radplots ################################################### plot(arn.rad, ylab="Biomass") ################################################### ### code chunk number 14: rads-relative-abundances ################################################### plot(moths.rad, prop = TRUE, type="n") lines(rad(birds), prop = TRUE, col="red") lines(rad(moths), prop = TRUE) legend("topright", c("Preston's birds", "Fisher's moths"), col=c("red", "blue"), lty=1, bty="n") ################################################### ### code chunk number 15: Fitting a log-series model ################################################### (moths.ls <- fitsad(moths,'ls')) ################################################### ### code chunk number 16: Operations on fitsad object ################################################### summary(moths.ls) coef(moths.ls) logLik(moths.ls) AIC(moths.ls) ################################################### ### code chunk number 17: Profiling and intervals ################################################### moths.ls.prf <- profile(moths.ls) likelregions(moths.ls.prf) #likelihood intervals confint(moths.ls.prf) ################################################### ### code chunk number 18: Ploting-profiles ################################################### par(mfrow=c(1,2)) plotprofmle(moths.ls.prf)# log-likelihood profile plot(moths.ls.prf)# z-transformed profile par(mfrow=c(1,1)) ################################################### ### code chunk number 19: Plot-of-predicted-values ################################################### par(mfrow=c(2,2)) plot(moths.ls) par(mfrow=c(1,1)) ################################################### ### code chunk number 20: Fitting two other models ################################################### (moths.pl <- fitsad(x=moths, sad="poilog"))#default is zero-truncated (moths.ln <- fitsad(x=moths, sad="lnorm", trunc=0.5)) # lognormal truncated at 0.5 ################################################### ### code chunk number 21: Model selection table ################################################### AICtab(moths.ls, moths.pl, moths.ln, base=TRUE) ################################################### ### code chunk number 22: Predicted values for octaves ################################################### head(moths.ls.oc <- octavpred(moths.ls)) head(moths.pl.oc <- octavpred(moths.pl)) head(moths.ln.oc <- octavpred(moths.ln)) ################################################### ### code chunk number 23: Octaves-plot ################################################### plot(moths.oc) lines(moths.ls.oc, col="blue") lines(moths.pl.oc, col="red") lines(moths.ln.oc, col="green") legend("topright", c("Logseries", "Poisson-lognormal", "Truncated lognormal"), lty=1, col=c("blue","red", "green")) ################################################### ### code chunk number 24: Predicted values - radplots ################################################### head(moths.ls.rad <- radpred(moths.ls)) head(moths.pl.rad <- radpred(moths.pl)) head(moths.ln.rad <- radpred(moths.ln)) ################################################### ### code chunk number 25: Rad-plots ################################################### plot(moths.rad) lines(moths.ls.rad, col="blue") lines(moths.pl.rad, col="red") lines(moths.ln.rad, col="green") legend("topright", c("Logseries", "Poisson-lognormal", "Truncated lognormal"), lty=1, col=c("blue","red", "green")) ################################################### ### code chunk number 26: grasslands dataset ################################################### head(grasslands) ################################################### ### code chunk number 27: grasslands breakpoints ################################################### (grass.brk <- c(0,1,3,5,seq(15,100, by=10),100) ) ################################################### ### code chunk number 28: grasslands histogram ################################################### grass.h <- hist(grasslands$mids, breaks = grass.brk, plot = FALSE) ################################################### ### code chunk number 29: grasslands histogram data ################################################### data.frame(midpoint = grass.h$mids, N.spp = grass.h$counts) ################################################### ### code chunk number 30: grasslands model fit ################################################### grass.e <- fitsadC(grass.h, 'exp') # Exponential grass.g <- fitsadC(grass.h, 'gamma') # Pareto grass.l <- fitsadC(grass.h, 'lnorm') # Log-normal grass.p <- fitsadC(grass.h, 'pareto') # Pareto grass.w <- fitsadC(grass.h, 'weibull') # Weibull ################################################### ### code chunk number 31: grasslands model selection ################################################### AICctab(grass.e, grass.g, grass.l, grass.p, grass.w, weights = TRUE, base = TRUE) ################################################### ### code chunk number 32: grasslands_hist_plot ################################################### plot(grass.h, main = "", xlab = "Abundance class") ################################################### ### code chunk number 33: grasslands coverpred ################################################### ## Predicted by each model grass.e.p <- coverpred(grass.e) grass.g.p <- coverpred(grass.g) grass.l.p <- coverpred(grass.l) grass.p.p <- coverpred(grass.p) grass.w.p <- coverpred(grass.w) ################################################### ### code chunk number 34: grasslands_coverpred_plot ################################################### ## Plot plot(grass.h, main = "", xlab = "Abundance class", xlim = c(0,40)) ## Adds predicted points points(grass.e.p, col = 1) points(grass.g.p, col = 2) points(grass.l.p, col = 3) points(grass.p.p, col = 4) points(grass.w.p, col = 5) legend("topright", legend = c("Exponential", "Gamma", "Log-normal", "Pareto", "Weibull"), col = 1:5, bty = "n", lty = 1 , pch = 1) ################################################### ### code chunk number 35: rsad-example1 ################################################### set.seed(42)# fix random seed to make example reproducible (samp1 <- rsad(S = 10, frac = 0.1, sad = "lnorm", coef=list(meanlog = 3, sdlog = 1.5), zeroes=TRUE, ssize = 2)) ################################################### ### code chunk number 36: rsad-example2 ################################################### (samp2 <- rsad(S = 100, frac=0.1, sad="lnorm", list(meanlog=5, sdlog=2))) ################################################### ### code chunk number 37: rsad-poilog-fit ################################################### (samp2.pl <- fitsad(samp2, "poilog")) ## checking correspondence of parameter mu coef(samp2.pl)[1] - log(0.1) ################################################### ### code chunk number 38: rsad-repeated-samples ################################################### results <- matrix(nrow=75,ncol=2) for(i in 1:75){ x <- rsad(S = 100, frac=0.1, sad="lnorm", list(meanlog=5, sdlog=2)) y <- fitsad(x, "poilog") results[i,] <- coef(y) } results[,1] <- results[,1]-log(0.1) ################################################### ### code chunk number 39: rsads_bias ################################################### ##Mean of estimates apply(results,2,mean) ## relative bias (c(5,2)-apply(results,2,mean))/c(5,2) ################################################### ### code chunk number 40: rsads_bias ################################################### ##Mean of estimates apply(results,2,sd) ## relative precision apply(results,2,sd)/apply(results,2,mean) ################################################### ### code chunk number 41: rsads-bias-plots ################################################### par(mfrow=c(1,2)) plot(density(results[,1]), main=expression(paste("Density of ",mu))) abline(v=c(mean(results[,1]),5), col=2:3) plot(density(results[,2]), main=expression(paste("Density of ",sigma))) abline(v=c(mean(results[,2]), 2), col=2:3) par(mfrow=c(1,1))