R Script illustrating basic SAD simulation, analysis and fitting

sads-test.R

# sads-test.R

package=function(somepackage)
{
  cpackage <- as.character(substitute(somepackage))
  if(!require(cpackage,character.only=TRUE)){
    install.packages(cpackage)
    library(cpackage,character.only=TRUE)
  }
}

package(sads)
package(vegan)

op=par(mfrow=c(4,2),ask=TRUE)

# Fully observed case
comm=rsad(S=1000,frac=1,sad="lnorm",meanlog=5,sdlog=2)
barplot(comm,xlab="Species",ylab="Abundance")
comm = comm[order(-comm)]

tad=as.data.frame(table(comm))
names(tad)=c("Abundance","# taxa")
barplot(tad[,2],names.arg=tad[,1],xlab="Abundance",ylab="# species",main="TAD")
tad$Abundance=as.numeric(as.character(tad$Abundance))

oc=octav(comm)
plot(oc,main="Preston plot")
plot(rad(comm),"Rank abundance")

# Fractionally observed case
comm=rsad(S=1000,frac=0.0002,sad="lnorm",meanlog=5,sdlog=2)
barplot(comm,xlab="Species",ylab="Abundance")
comm = comm[order(-comm)]

tad=as.data.frame(table(comm))
names(tad)=c("Abundance","# taxa")
barplot(tad[,2],names.arg=tad[,1],xlab="Abundance",ylab="# species",main="TAD")
tad$Abundance=as.numeric(as.character(tad$Abundance))

oc=octav(comm)
plot(oc,main="Preston plot")
plot(rad(comm),"Rank abundance")

# Estimate number of unobserved species
print(estimateR(comm))

# Fit TAD/SAD
mod = rad.lognormal(comm) # vegan
mod
#plot(mod)

mod = fitsad(comm) # sads
mod
plot(mod)

# poisson-log model will fit sensibly to a fractional sample...
mod=fitsad(comm,"poilog")
# should be 5+log(0.0001)=-4.2 and 2
print(mod)
plot(mod)
# this seems to work quite well

par(op)

# eof