################################################### ### chunk number 1: ################################################### setwd("~/Documents/Projects/Book/draft/Chap10") rm(list=ls()) badfonts = FALSE .detach.stuff = function() { s1 = grep("package|Autoloads",search()) nattach = length(search()) xattach = search()[c(-1,-s1)] for (i in xattach) eval(substitute(detach(i),list(i=i))) } .detach.stuff() ## ps.options(family="NimbusSan") ps.options(colormodel="cmyk") options(width=60, SweaveHooks=list(fig=function() par(mar=c(5,4,.5,1.5))) ) ## palette(gray((0:8)/8)) library(primer) trellis.par.set(canonical.theme(color=FALSE)) ################################################### ### chunk number 2: BSsuc ################################################### bsdiv <- read.csv('C3rank.csv') matplot(bsdiv[,1], bsdiv[,c(2:4, 7,9,11)], type='l', lty=1:3, lwd=rep(1:2, each=3), col=1, log="y", xlab="Abundance Rank", ylab="Relative Abundance") legend('topright', paste("Year ", c(1,2,4,15,25,35), sep=""), ncol=2, lty=1:3, lwd=rep(1:2, each=3), bty='n') ?moths ################################################### ### chunk number 3: ################################################### dens <- data.frame(Salwho=c(1,1,2,3), Fravir=c(21,8,13,5)) row.names(dens) <- LETTERS[1:4] dens ################################################### ### chunk number 4: dens ################################################### plot(dens, type='n'); text(dens, row.names(dens)) ################################################### ### chunk number 5: ################################################### dens[1,]/sum(dens[1,]) ################################################### ### chunk number 6: ################################################### dens[,1]/sum(dens[,1]) ################################################### ### chunk number 7: ################################################### x <- dens[2,1] - dens[3,1] ################################################### ### chunk number 8: ################################################### y <- dens[2,2] - dens[3,2] ################################################### ### chunk number 9: ################################################### sqrt(x^2 + y^2) ################################################### ### chunk number 10: ################################################### (alldists <- dist(dens)) ################################################### ### chunk number 11: dens3 ################################################### dens[["Manoff"]] <- c(11,3,7,5) (spp3 <- dist(dens)) ################################################### ### chunk number 12: ################################################### pairs(dens)# not shown library(scatterplot3d) sc1 <- scatterplot3d(dens, type='h', pch="", xlim=c(0,5), ylim=c(0, 25), zlim=c(0,15)) text(sc1$xyz.convert(dens), labels=rownames(dens)) ################################################### ### chunk number 13: ################################################### dens$Acolyc <- c(16, 0, 9,4) ################################################### ### chunk number 14: mds1 ################################################### library(vegan) mdsE <- metaMDS(dens, distance="euc", autotransform=FALSE, trace=0) plot(mdsE, display="sites", type="text") ################################################### ### chunk number 15: mds2 ################################################### mdsB <- metaMDS(dens, distance="bray", autotransform=FALSE, trace=0) plot(mdsB, display="sites", type="text") ################################################### ### chunk number 16: ################################################### (dens.RA <- t( apply(dens, 1, function(sp.abun) sp.abun/sum(sp.abun) )) ) ################################################### ### chunk number 17: ################################################### (mins <- apply(dens.RA[1:2,], 2, min)) ################################################### ### chunk number 18: ################################################### sum(mins) *100 ################################################### ### chunk number 19: ################################################### (shared <- apply(dens[1:2,], 2, function(abuns) all(abuns != 0) ) ) ################################################### ### chunk number 20: ################################################### (Rs <- apply(dens[1:2,], 1, function(x) sum(x>0)) ) ################################################### ### chunk number 21: ################################################### 2*sum(shared)/sum(Rs) ################################################### ### chunk number 22: ################################################### H <- function(x) {x <- x[x>0] p <- x/sum(x) - sum (p * log(p) ) } Sd <- function(x) {p <- x/sum(x) 1 - sum(p^2)} h1 <- H(x=rep(.25,4)); h2<- H(x=c(50,75,75));h3<- H(c(20,20,20,140)); h4<- H(x=200) sd1 <- Sd(x=rep(.25,4)); sd2<- Sd(x=c(50,75,75));sd3<- Sd(c(20,20,20,140)); ################################################### ### chunk number 23: ################################################### H <- function(x) {x <- x[x>0] p <- x/sum(x) - sum (p * log(p) ) } Sd <- function(x) {p <- x/sum(x) 1 - sum(p^2)} ################################################### ### chunk number 24: ################################################### h1 <- H(x=rep(.25,4)); h2<- H(x=c(50,75,75));h3<- H(c(1,1,1,197)); h4<- H(x=200) sd1 <- Sd(x=rep(.25,4)); sd2<- Sd(x=c(50,75,75));sd3<- Sd(c(1,1,1,197)); ################################################### ### chunk number 25: EqAb ################################################### Rs <- 2:20 ComsEq <- sapply(Rs, function(R) (1:R)/R) Hs <- sapply(ComsEq, H); Sds <- sapply(ComsEq, Sd) plot(Rs, Hs, type='l', ylab="Diversity", ylim=c(0,3)); lines(Rs, Sds, lty=2); legend("right", c(expression(italic("H")), expression(italic("S"["D"]))), lty=1:2, bty='n') ################################################### ### chunk number 26: Ab90 ################################################### Coms90 <- sapply(Rs, function(R) {p <- numeric(R) p[1] <- .9 p[2:R] <- .1/(R-1) p} ) Hs <- sapply(Coms90, H); Sds <- sapply(Coms90, Sd) plot(Rs, Hs, type='l', ylim=c(0,.7)); lines(Rs, Sds, lty=2) ################################################### ### chunk number 27: ################################################### s1 <- matrix(c(1,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,1,1), nr=6) colnames(s1) <- c("Sp.A", "Sp.B", "Sp.C");s1 ################################################### ### chunk number 28: simpson1 ################################################### library(scatterplot3d) s13d <- scatterplot3d(jitter(s1, .3), type='h', angle=60, pch=c(1,1,2,2,3,3), # main="Simpson's Div = 0.67", xlim=c(-.2, 1.4), ylim=c(-.2, 1.4), zlim=c(-.2, 1.4)) s13d$points3d(x=1/3, y=1/3, z=1/3, type='h', pch=19, cex=2) ################################################### ### chunk number 29: ################################################### (centroid1 <- colMeans(s1)) ################################################### ### chunk number 30: ################################################### (SS <- sum( sapply(1:3, function(j) (s1[,j] - centroid1[j]))^2) ) ################################################### ### chunk number 31: ################################################### SS/6 ################################################### ### chunk number 32: ################################################### p <- c(2,2,2)/6 1-sum(p^2) ################################################### ### chunk number 33: simpson2 ################################################### s2 <- matrix(c(1,1,1,1,0,0,0,0,0,0,1,0,0,0,0,0,0,1), nr=6) colnames(s2) <- c("Sp.A", "Sp.B", "Sp.C");s2 centroid2 <- colMeans(s2) s23d <- scatterplot3d(jitter(s2, .2), type='h', angle=60, pch=c(1,1,1,1,2,3), # main="Simpson's Div = 0.50", xlim=c(-.2, 1.4), ylim=c(-.2, 1.4),zlim=c(-.2, 1.4)) s23d$points3d(x=4/6, y=1/6, z=1/6, type='h', pch=19, cex=2) ################################################### ### chunk number 34: simpson3 ################################################### s3 <- matrix(c(1,1,1,1,1,0,0,0,0,0,0,1,0,0,0,0,0,0), nr=6) colnames(s3) <- c("Sp.A", "Sp.B", "Sp.C");s3 (centroid3 <- colMeans(s3)) s33d <- scatterplot3d(jitter(s3, .3), type='h', angle=60, pch=c(1,1,1,1,1,2), xlim=c(-.2, 1.4), ylim=c(-.2, 1.4), zlim=c(-.2, 1.4) # main="Simpson's Div = 0.28" ) s33d$points3d(x=5/6, y=1/6, z=0, type='h', pch=19, cex=2) ################################################### ### chunk number 35: rarefaction eval=FALSE ################################################### ## ### Here I play around with stuff ## library(vegan) ## data(BCI); ## N <- sort(colSums(BCI), decr=TRUE) ## subs1 <- c( seq(50, 400, by=50), sum(BCI[1,])) ## rar1 <- rarefy(BCI[1,], MARGIN=1, sample=subs1, se=T) ## ## subs2 <- c( seq(50, 400, by=50), sum(BCI[2,])) ## rar2 <- rarefy(BCI[2,], MARGIN=1, sample=subs2, se=T) ## ## plot(subs1, rar1[1,], type='b', ylim=c(0, 100), ylab="Individual Rarefied Richness", ## xlab="No. of Individuals") ## lines(subs1, rar1[1,]+2*rar1[2,], lty=3) ## lines(subs1, rar1[1,]-2*rar1[2,], lty=3) ## lines(subs2, rar2[1,], type="b") ## ################################################### ### chunk number 36: ################################################### library(vegan) data(BCI); bci <- BCI[seq(5,50, by=5),] ################################################### ### chunk number 37: est ################################################### N <- colSums(bci) NoOfInds <- c(seq(500, 4500, by=500), sum(N)) rar3 <- rarefy(N, sample=NoOfInds, se=T, MARG=2) plot(NoOfInds, rar3[1,], ylab="Species Richness", axes=FALSE, xlab="No. of Individuals", type='n', ylim=c(0,250), xlim=c(500,7000) ) axis(1, at=1:5*1000); axis(2); box(); text(2500, 200, "Individual-based rarefaction (10 plots)") lines(NoOfInds, rar3[1,], type='b') lines(NoOfInds, rar3[1,]+2*rar3[2,], lty=3) lines(NoOfInds, rar3[1,]-2*rar3[2,], lty=3) ace <- estimateR(N) segments(6000, ace['S.ACE']-2*ace['se.ACE'], 6000, ace['S.ACE']+2*ace['se.ACE'], lwd=3) text(6000, 150, "ACE estimate", srt=90, adj=c(1,.5) ) chaoF <- specpool(bci) segments(6300, chaoF[1,'Chao'] - 2*chaoF[1,'Chao.SE'], 6300, chaoF[1,'Chao'] + 2*chaoF[1,'Chao.SE'], lwd=3, col='grey') text(6300, 150, "Chao2 estimate", srt=90, adj=c(1,.5) ) points(6700, dim(BCI)[2], pch=19, cex=1.5) text(6700, 150, "Richness observed in 50 ha", srt=90, adj=c(1,.5) ) ################################################### ### chunk number 38: ################################################### N <- colSums(bci) subs3 <- c(seq(500, 4500, by=500), sum(N)) rar3 <- rarefy(N, sample=subs3, se=T, MARG=2) ################################################### ### chunk number 39: ################################################### plot(subs3, rar3[1,], ylab="Species Richness", axes=FALSE, xlab="No. of Individuals", type='n', ylim=c(0,250), xlim=c(500,7000) ) axis(1, at=1:5*1000); axis(2); box(); text(2500, 200, "Individual-based rarefaction (10 plots)") ################################################### ### chunk number 40: ################################################### lines(subs3, rar3[1,], type='b') lines(subs3, rar3[1,]+2*rar3[2,], lty=3) lines(subs3, rar3[1,]-2*rar3[2,], lty=3) ################################################### ### chunk number 41: ################################################### ace <- estimateR(N) segments(6000, ace['S.ACE']-2*ace['se.ACE'], 6000, ace['S.ACE']+2*ace['se.ACE'], lwd=3) text(6000, 150, "ACE estimate", srt=90, adj=c(1,.5) ) ################################################### ### chunk number 42: ################################################### chaoF <- specpool(bci) segments(6300, chaoF[1,'Chao'] - 2*chaoF[1,'Chao.SE'], 6300, chaoF[1,'Chao'] + 2*chaoF[1,'Chao.SE'], lwd=3, col='grey') text(6300, 150, "Chao2 estimate", srt=90, adj=c(1,.5) ) ################################################### ### chunk number 43: ################################################### points(6700, dim(BCI)[2], pch=19, cex=1.5) text(6700, 150, "Richness observed in 50 ha", srt=90, adj=c(1,.5) ) ################################################### ### chunk number 44: bci1 ################################################### data(BCI); N <- sort(colSums(BCI), decr=TRUE) hist(N, main=NULL) ################################################### ### chunk number 45: bci2 ################################################### hist(log(N), xlab="Log Density Classes", main=NULL) m.spp <- mean(log(N)); sd.spp <- sd(log(N)); R <- length(N) curve(dnorm(x,m.spp,sd.spp)*R, 0, 8, add=T) ################################################### ### chunk number 46: bci3 ################################################### plot(log(N), type='b', ylim=c(0,8), main=NULL, xlab="Species Rank", ylab="Log Density") ranks.lognormal <- R * (1-pnorm( log(N), m.spp, sd.spp)) lines(ranks.lognormal, log(N)) ################################################### ### chunk number 47: ################################################### N <- sort(colSums(BCI), decr=TRUE) f1 <- sort(cumsum(1/(R:1)), decr=TRUE) Nt <- sum(N) NMac <- Nt*f1/sum(f1) ################################################### ### chunk number 48: ################################################### d <- N[1]/Nt Ngeo.f <- d*(1-d)^(0:(R-1)) Ngeo <- Nt * Ngeo.f ################################################### ### chunk number 49: ################################################### library(untb) alpha <- optimal.theta(N) ################################################### ### chunk number 50: ################################################### sei <- function( t=1 ) exp(-t)/t ### standard expo integral alpha <- optimal.theta(N) ranks.logseries <- sapply(N, function(x) { n <- x * log(1 + alpha/Nt ) f <- integrate(sei, n, Inf) fv <- f[["value"]] alpha * fv } ) ################################################### ### chunk number 51: other ################################################### plot(1:R, N, ylim=c(.1,2000), xlim=c(1,301), axes=FALSE, log='y') axis(2); axis(1, 1 + seq(0, 300, by=50));box() lines(1:R, NMac, lty=4) lines(1:R, Ngeo, lty=3) lines(ranks.logseries, N, lty=2) lines(ranks.lognormal, N, lty=1) legend("topright", c("Broken Stick", "Log-series", "Log-normal", "Geometric"), lty=c(4,2,1,3), bty='n' ) ################################################### ### chunk number 52: metacomm eval=FALSE ################################################### ## library(untb) ## par(mar=c(1,1,1,1)) ## display.untb(start=rep(10, 900), prob=0.1, gens=10000) ## #box() ## polygon(c(6.5,13.5,13.5,6.5), y=c(6.5,6.5,13.5,13.5), lwd=6) ################################################### ### chunk number 53: driftfig ################################################### out <- untb(rep(1:10, 90), prob=0, gens=1000, D=9, keep=TRUE) matplot(1:1000,species.table(out), type='l', lwd=1.5, ylab="Population Density", xlab="Generation") ################################################### ### chunk number 54: ################################################### library(untb) a <- untb(start=rep(1:20, 25), prob=0, gens=2500, D=9, keep=TRUE) ################################################### ### chunk number 55: ################################################### ( a2 <- a[901:903, 1:10 ] ) ################################################### ### chunk number 56: ################################################### times <- c(1,50,2500) sppcolors <- sapply(times, function(i) grey((a[i,]-1)/20) ) ################################################### ### chunk number 57: communities ################################################### layout(matrix(1:3, nr=1)) par(mar=c(1,1,3,1)) for(j in 1:3){ plot(c(1,20), c(1,25), type="n", axes=FALSE) points(rep(1:20,25), rep(1:25,each=20), pch=19, cex=2, col=sppcolors[,j]) title(paste("Time = ", times[j], sep="")) } ################################################### ### chunk number 58: comms ################################################### layout(matrix(1:9, nr=3, byrow=TRUE)) par(mar=c(1,1,3,1)) for(j in 1:3){ plot(c(1,20), c(1,25), type="n", axes=FALSE) points(rep(1:20,25), rep(1:25,each=20), pch=19, cex=1.5, col=sppcolors[,j]) title(paste("Time = ", times[j], sep="")) } par(mar=c(5,4,1,1)) for(i in times){ plot(as.count(a[i,]), ylim=c(1,max(as.count(a[2000,]))), xlim=c(0,20)) } par(mar=c(6,4,1,3)) out <- lapply(times, function(i) preston(a[i,]) ) bins <- matrix(0, nrow=3, ncol=length(out[[3]]) ) for(i in 1:3) bins[i,1:length(out[[i]])] <- out[[i]] bins colnames(bins) <- names(preston(a[times[3],])) for(i in 1:3){ par(las=2) barplot(bins[i,], ylim=c(0,20), xlim=c(0,8), ylab="No. of Species", xlab="Abundance Category" ) } ################################################### ### chunk number 59: ################################################### layout(matrix(1:3, nr=1)) for(i in times){ plot(as.count(a[i,]), ylim=c(1,max(as.count(a[times[3],]))), xlim=c(0,20)) title(paste("Time = ", i, sep="")) } ################################################### ### chunk number 60: ################################################### out <- lapply(times, function(i) preston(a[i,]) ) bins <- matrix(0, nrow=3, ncol=length(out[[3]]) ) for(i in 1:3) bins[i,1:length(out[[i]])] <- out[[i]] bins colnames(bins) <- names(preston(a[times[3],])) ################################################### ### chunk number 61: ################################################### layout(matrix(1:3, nr=1)) for(i in 1:3){ par(las=2) barplot(bins[i,], ylim=c(0,20), xlim=c(0,8), ylab="No. of Species", xlab="Abundance Category" ) } ################################################### ### chunk number 62: spptable ################################################### sppmat <- species.table(a) matplot(1:times[3], sppmat[1:times[3],], type='l', ylab="Population Density") ################################################### ### chunk number 63: ################################################### cvtimes <- seq(2, 2500, by=10) CV.within <- sapply(cvtimes, function(i) {# Use sapply with i as a row counter ## For each row i, use apply to caluculate the variance for each ## population from 1:i cvs <- apply(sppmat[1:i,],2, function(x) sd(x)/mean(x) ) mean(cvs) # get the mean among populations (for each i). } ) #popvar.among <- apply(sppmat, 1, var) ################################################### ### chunk number 64: vartime ################################################### plot( cvtimes, CV.within, type='l' ) ################################################### ### chunk number 65: ################################################### library(vegan) data(BCI) n <- colSums(BCI) par(las=2, mar=c(7,4,1,1)) xs <- plot(preston(rep(1:225,n), n=12, original=TRUE)) ################################################### ### chunk number 66: ################################################### v1 <- volkov(sum(n), c(48, 0.1), bins=TRUE) points(xs, v1[1:12], type='b') ################################################### ### chunk number 67: ################################################### v2 <- volkov(sum(n), c(48, 0.14), bins=TRUE) xs <- plot(preston(rep(1:225,n), n=12, original=TRUE), col=0, border=0) axis(1, at=xs, labels=FALSE) points(xs, v2[1:12], type='b', lty=2, col=2) ################################################### ### chunk number 68: ################################################### v4 <- volkov(sum(n), c(571, 0.002), bins=TRUE) points(xs, v4[1:12], type='b', lty=4, col=2) ################################################### ### chunk number 69: bcipreston ################################################### par(las=2, mar=c(5,4,1,1)) xs <- plot(preston(rep(1:225,n), n=12, original=TRUE)) points(xs, v1[1:12], type='b') ################################################### ### chunk number 70: bciJabot ################################################### par(las=2, mar=c(5,4,1,1)) xs <- plot(preston(rep(1:225,n), n=12, original=TRUE), col=0, border=0) axis(1, at=xs, labels=FALSE) points(xs, v2[1:12], type='b') points(xs, v4[1:12], type='b', lty=2) ################################################### ### chunk number 71: ################################################### library(maps) library(ellipse) N <- 1000 crd <- map(, xlim=c(12,42), ylim=c(-5,5)) coords <- matrix(NA, nr=N, ncol=2) for(i in 1:N) coords[i,] <- { xr <- range(crd$x[!is.na(crd$x)]); yr <- range(crd$y[!is.na(crd$y)]) x <- runif(1, xr[1], xr[2]); y <- runif(1,yr[1], yr[2]) c(x,y) } ### Pick letters relative to latitude a <- (coords[,1] + .1 - min(coords[,1] ) )*2 b <- max(a) - a +.1 inds <- sapply(1:N, function(i) { sample(1:26, 1, replace=T, prob=dbeta(seq(0.001,.999,length=26), a[i], b[i])) }) inds <- sapply(1:N, function(i) { sample(LETTERS, 1, replace=T, prob=dbeta(seq(0.001,.999,length=26), a[i], b[i])) }) rws <- rep(c(-2,2), each=6) cls <- rep(seq(15,40,5), 2) x1s <- coords[1,] - rws[1] pts <- ellipse(0, centre=c(0,0), level=.5) sqrt(apply(pts^2,1,sum)) #quartz(height=7, width=7) ################################################### ### chunk number 72: landscape ################################################### library(ellipse) #install.packages("maps", dep=TRUE) library(maps) map.text('world', xlim=c(12,42), ylim=c(-5,5)) map.axes() for(i in 1:12) { pts <- ellipse(0, centre=c(cls[i], rws[i]), level=.5, type="l" ) lines(pts) } points(coords, pch=inds, cex = .5) ################################################### ### chunk number 73: beta ################################################### a <- seq(1,100, by=.1) g <- max(a) bA <- g-a bM <- g/a matplot(a, cbind(bA, bM), type='l', xlab=quote(bar(italic(alpha))), ylab=quote(italic(beta)), col=1) ################################################### ### chunk number 74: ################################################### (gamma <- dim(BCI)[2]) (alpha.bar <- mean( specnumber(BCI) )) ################################################### ### chunk number 75: beta ################################################### (beta.A <- gamma - alpha.bar) (beta.M <- gamma/alpha.bar) ################################################### ### chunk number 76: moths ################################################### moths <- read.table("Summervillle_CristMothsALL.txt", header=TRUE) rowsnct <- c(4,7,9,11,6,3,5,1,10,12,8,2) rowswap<- c(16, 20,21,14,17,19,15,18,13) #identify mothsnct <- moths[rowsnct,] mothswap <- moths[rowswap,] library(maps) #wap <- map.text("state",xlim=c(-83, -82.3), ylim=c(39.4,39.7)) #map.text("state",xlim=c(-85.1, -84.5), ylim=c(39.4,39.6)) #text(moths[['long']], moths[['lat']], moths[['site']], cex=.5) #identify(wap) #locator() map.text("state",xlim=c(-85.5, -82), ylim=c(38.9,39.95)) text(-85.5, 39.9, labels="Indiana", adj=c(0,1)) text(-83, 39.9, labels="Ohio", adj=c(1,1)) text(-84.5, 38.5, labels="Kentucky", adj=c(0,0)) map.axes() points(moths[['long']], moths[['lat']], pch=1) n <- dim(mothsnct)[1] x <- cos( seq(pi, 1.25*2*pi, length=n) ) y <- sin(seq(pi, 1.25*2*pi, length=n) ) x2 <- .3*x + mean(mothsnct[['long']]) y2 <- .25*y + mean(mothsnct[['lat']]) segments(mothsnct[['long']], mothsnct[['lat']], x2, y2) x3 <- .35*x + mean(mothsnct[['long']]) y3 <- .3*y + mean(mothsnct[['lat']]) text(x3,y3, mothsnct[['spp']]) n <- dim(mothswap)[1] x <- cos( seq(pi, 1.25*2*pi, length=n) ) y <- sin(seq(pi, 1.25*2*pi, length=n) ) x2 <- .3*x + mean(mothswap[['long']]) y2 <- .25*y + mean(mothswap[['lat']]) x2 <- .3*x + -82.531 y2 <- .25*y + 39.584 segments(mothswap[['long']], mothswap[['lat']], x2, y2) x3 <- .35*x + -82.531 y3 <- .3*y + 39.584 text(x3,y3, mothswap[['spp']], adj=c(.75,0.5)) ################################################### ### chunk number 77: eval=FALSE ################################################### ## data(moths) ################################################### ### chunk number 78: ################################################### a1 <- mean(moths[['spp']]) ################################################### ### chunk number 79: ################################################### a2 <- sum(c(NCT=179, WAP=173) * c(12,9)/21) g <- 230 ################################################### ### chunk number 80: ################################################### b1 <- a2-a1 b2 <- g-a2 abg <- c(a1=a1, b1=b1, a2=a2, b2=b2, g=g) abg a1 + b1 + b2 == g ################################################### ### chunk number 81: partplot ################################################### par(las=2) bp <- barplot(abg[c(1,2,4,5)], hor=T, xlim=c(0, 250), space=1, names=expression(alpha[1],beta[1],beta[2],gamma)) text(abg[1]/2, bp[1], c("Mean site\nrichness"), cex=.8) text(abg[2]/2, bp[2], c("Species not in sites,\nbut within ecoregion"), cex=.8) text(abg[4]/2, bp[3], c("Species not\nin ecoregions"), cex=.8) text(abg[5]/2, bp[4], c("Regional richness (all sampled species)"), cex=.8) ################################################### ### chunk number 82: SARs1 ################################################### A <- 10^1:10; c <- 1.5; z <- 0.25 curve(c*x^z, 10, 10^10, n=500, ylab="No. of Species", xlab="Area (ha)") ################################################### ### chunk number 83: SARs2 ################################################### A <- 10^1:10; c <- 1.5; z <- 0.25 curve(log(c,10) + z*x, 1, 10, ylab=quote(log[10]("No. of Species")), xlab=quote(log[10]("Area (ha)")) ) ################################################### ### chunk number 84: SARmoths ################################################### moths <- read.table("Summervillle_CristMothsALL.txt", header=TRUE) names(moths) plot(log(spp,10) ~ log(area,10), moths) mod <- lm(log(spp,10) ~log(area,10), data=moths) abline(mod) try(rm("a")); try(rm("z")) mod.nonlin <- nls(spp ~ a*area^z, start=list(a=1, z=.2),data=moths) curve(log(coef(mod.nonlin)[1],10) + x*coef(mod.nonlin)[2], 0, 3, add=TRUE, lty=2) ################################################### ### chunk number 85: ################################################### plot(log(spp,10) ~ log(area,10), moths) mod <- lm(log(spp,10) ~log(area,10), data=moths) abline(mod) ################################################### ### chunk number 86: ################################################### mod.nonlin <- nls(spp ~ a*area^z, start=list(a=1, z=.2),data=moths) curve(log(coef(mod.nonlin)[1],10) + x*coef(mod.nonlin)[2], 0, 3, add=TRUE, lty=2) ################################################### ### chunk number 87: ################################################### confint(mod) confint(mod.nonlin) ################################################### ### chunk number 88: IE ################################################### I0 <- log(1) b <- .1 curve(exp(I0 - b*x),0, 50, xlab='No. of Species (R)', ylab="Rate (I or E)") d <- .01 curve(exp(d*x)-1, 0, 50, add=TRUE) deltaR <- function(R) { (exp(I0-b*R) - (exp(d*R)-1))^2 } S1 <- optimize(f=deltaR, interval=c(1,50))$minimum I1 <- exp(I0-b*S1) segments(S1, -.1, S1, I1) I0 <- log(1/2) curve(exp(I0-b*x), 0,50, add=TRUE, lty=2) d <- .014 curve(exp(d*x)-1, 0, 50, add=TRUE, lty=2) S2 <- optimize(f=deltaR, interval=c(1,50))$minimum I2 <- exp(I0-b*S2) segments(S2, -.1, S2, I2, lty=2) ################################################### ### chunk number 89: ################################################### I0 <- log(1) b <- .1 curve(exp(I0 - b*x),0, 50, xlab='No. of Species (R)', ylab="Rate (I or E)") ################################################### ### chunk number 90: ################################################### d <- .01 curve(exp(d*x)-1, 0, 50, add=TRUE) ################################################### ### chunk number 91: ################################################### deltaR <- function(R) { (exp(I0-b*R) - (exp(d*R)-1))^2 } ################################################### ### chunk number 92: ################################################### optimize(f=deltaR, interval=c(1,50))[["minimum"]] ################################################### ### chunk number 93: ################################################### I0 <- log(1/2) curve(exp(I0 - b*x), 0,50, add=TRUE, lty=2) ################################################### ### chunk number 94: ################################################### d <- .014 curve(exp(d*x)-1, 0, 50, add=TRUE, lty=2) ################################################### ### chunk number 95: ################################################### optimize(f=deltaR, interval=c(1,50))[["minimum"]] ################################################### ### chunk number 96: ################################################### plot(spp ~ area, data=moths, ylim=c(30,230), #log='y', xlab="Area (ha)", ylab="No. of Species (R)") curve(coef(mod.nonlin)[1] * x^coef(mod.nonlin)[2], 0, max(moths[['area']]), add=TRUE, lty=2, lwd=2) #text(150, 200, # bquote(italic("R") == .(round(coef(mod.nonlin)[1],1)) * # italic("A")^.(round(coef(mod.nonlin)[2],2)) ) # ) abline(h=g, lty=3) text(275, g, quote(gamma), adj=c(.5,1.5), cex=1.5) ################################################### ### chunk number 97: predR ################################################### (MaxR <- predict(mod.nonlin, list(area=max(moths[['area']]) )) ) ################################################### ### chunk number 98: ################################################### b.area <- MaxR - a1 b.eco <- a2-(b.area+a1) b.geo <- g - a2 ################################################### ### chunk number 99: ################################################### abline(h=g, lty=3) abline(h=b.eco+b.area+a1, lty=3) abline(h=b.area+a1, lty=3) abline(h=a1, lty=3) ################################################### ### chunk number 100: SARpart ################################################### par(mar=c(5,4,1,3), las=1) plot(spp ~ area, data=moths, ylim=c(0,250), #log='y', xlab="Area (ha)", ylab="No. of Species (R)") cols <- hcl(c(150), alpha=c(1,.5,.3,.1)) polygon(c(-20,320,320,-200), c(0,0,a1,a1), border=FALSE, col=cols[1]) polygon(c(-20,320,320,-200), c(0,0,a1+b.area,a1+b.area), border=FALSE, col=cols[2]) polygon(c(-20,320,320,-200), c(0,0,a1+b.area+b.eco,a1+b.area+b.eco), border=FALSE, col=cols[3]) polygon(c(-20,320,320,-200), c(0,0,g,g), border=FALSE, col=cols[4]) box() points(moths[['area']], moths[['spp']]) curve(coef(mod.nonlin)[1] * x^coef(mod.nonlin)[2], min(moths[['area']])/2, max(moths[['area']]), add=TRUE) abline(h=g, lty=3) abline(h=b.eco+b.area+a1, lty=3) abline(h=b.area+a1, lty=3) abline(h=a1, lty=3) text(225, g, quote(gamma), adj=c(-.5,0), cex=1.5) arrows(225, 0, 225, g, code=3, length=.05) #title(sub=quote(beta["replace"]==beta['eco']+beta['geo'])) #mtext(quote(beta["replace"]==beta['eco']+beta['geo'])) #text(175, 300, quote(beta["replace"]==beta['eco']+beta['geo']), cex=1.2) #arrows(200, MaxR, 200, g, code=3, length=.05) mtext(quote(alpha["1"]), side=4, at=a1, line=.5, cex=1.2) mtext(quote(alpha["2"]), side=4, at=a2, line=.5, cex=1.2) mtext(quote(italic("R"["max"])), side=4, at=a1+b.area, line=.5) #text(280, a1, quote(alpha["1"]), adj=c(0,0), cex=1.2) #text(280, a2, quote(alpha["2"]), adj=c(0,0), cex=1.2) text(150, a2+b2/2, quote(beta["geo"]), adj=c(1.1,0.5), cex=1.2) arrows(150, a1+b.area+b.eco, 150, g, code=3, length=.05) text(100, a1+b1/2, quote(beta["eco"]), adj=c(1.1,0.5), cex=1.2) arrows(100, a1+b.area, 100, a1+b.area+b.eco, code=3, length=.05) text(25, a1+b.area/2, quote(beta["area"]), adj=c(1.1,.5), cex=1.2) arrows(32, a1, 32, a1 + b.area, code=3, length=.05) text(32, a1/2, quote(bar(alpha)), adj=c(1.2,.5), cex=1.2) arrows(32, 0, 32, a1, code=3, length=.05)