################################################### ### chunk number 1: ################################################### setwd("~/Documents/Projects/Book/draft/Chap01") options(width=75, SweaveHooks=list(fig=function() par(mar=c(5,4,.5,1.5))) ) 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") ## palette(gray((0:8)/8)) library(lattice) library(primer) trellis.par.set(canonical.theme(color=FALSE)) #sparrows <- read.csv("SSparrow.csv") ################################################### ### chunk number 2: Melospiza1 ################################################### #myQ() #par(xpd=TRUE) data(sparrows) plot(Count~Year, data=sparrows[1:6,],type="b") #dev.print(pdf, "Melospiza1.pdf") ################################################### ### chunk number 3: DI ################################################### ssgr <- with(sparrows[1:6,], Count[2:6]/Count[1:5]) plot(ssgr ~ Count, data=sparrows[1:5,], ylab="Growth Rate") #dev.print(pdf, "DI.pdf") ################################################### ### chunk number 4: lily ################################################### N <- c(1, 3,9,27,81) year <- 2001:2005 plot(year, N) ################################################### ### chunk number 5: ################################################### rates = N[2:5] / N[1:4] rates ################################################### ### chunk number 6: ################################################### N0 <- 1 lambda <- 2 time <- 0:10 ################################################### ### chunk number 7: ################################################### Nt <- N0 * lambda^time Nt ################################################### ### chunk number 8: ################################################### N0 <- c(10, 20, 30) lambda <- 2 time <- 0:4 ################################################### ### chunk number 9: ################################################### Nt.s <- sapply(N0, function(n) n*lambda^time) Nt.s ################################################### ### chunk number 10: NtInitN ################################################### matplot(time, Nt.s, pch=1:3) ################################################### ### chunk number 11: NtInitNlog ################################################### matplot(time, Nt.s, log="y", pch=1:3) ################################################### ### chunk number 12: ################################################### N0 <- 100 time <- 0:3 lambdas <- c(0.5, 1, 1.5) ################################################### ### chunk number 13: ################################################### N.all <- sapply(lambdas, function(x) N0*x^time) ################################################### ### chunk number 14: NtLambda ################################################### matplot(time, N.all, xlab="Years", ylab="N", pch=1:3) abline(h = N0, lty=3) text(.5, 250, expression(lambda > 1.0), cex=1.2) text(.5, 20, expression(lambda < 1.0), cex=1.2) ################################################### ### chunk number 15: ################################################### t <- 5 SS6 <- sparrows[1:(t+1),] ################################################### ### chunk number 16: ################################################### SSgr <- SS6$Count[2:(t+1)] / SS6$Count[1:t] lam.A <- sum(SSgr)/t lam.G <- prod(SSgr)^(1/t) ################################################### ### chunk number 17: SSaves ################################################### N0 <- SS6$Count[1] plot(0:t, SS6$Count, ylab="Projected Population Size") lines(0:t, N0*lam.A^(0:t), lty=2 ) lines(0:t, N0*lam.G^(0:t), lty=1 ) legend(0,70, c("Arithmetic Ave.", "Geometric Ave."), title="Projections Based On:", lty=2:1, bty='n', xjust=0) ################################################### ### chunk number 18: e1 ################################################### n <- 0:100 N0 <- 1 rd <- 1 ################################################### ### chunk number 19: e2 ################################################### N1 <- N0 * (1+rd/n)^n ################################################### ### chunk number 20: e ################################################### plot(n, N1/N0, type="l") text(50, 2, "For n = 100,") text(50, 1.6, bquote((1 + frac("r"["d"],"n"))^"n" == .(round(N1[101]/N0,3)) ) ) ################################################### ### chunk number 21: ################################################### r <- c(-.03, -.02, 0, .02, .03) N0 <- 2 t <- 1:100 cont.mat <- sapply(r, function(ri) N0*exp(ri*t) ) ################################################### ### chunk number 22: cont ################################################### layout(matrix(1:2, nrow=1)) matplot(t, cont.mat, type="l", ylab="N", col=1) legend("topleft", paste(rev(r)), lty=5:1, col=1, bty="n", title="r") matplot(t, cont.mat, type="l", ylab="N", log="y", col=1) ################################################### ### chunk number 23: ################################################### m.time <- function(r, m=2) { log(m)/r } ################################################### ### chunk number 24: ################################################### rs <- c(0,1,2) m.time(rs) ################################################### ### chunk number 25: SongSparrowSim ################################################### names(sparrows) ################################################### ### chunk number 26: ################################################### attach(sparrows) ################################################### ### chunk number 27: ################################################### quartz(,7,3.5, dpi=100) layout(matrix(1:2, nr=1)) ################################################### ### chunk number 28: ################################################### plot(Count ~ Year, type="b") ################################################### ### chunk number 29: ################################################### obs.R <- Count[-1]/Count[-length(Count)] ################################################### ### chunk number 30: ################################################### plot(obs.R ~ Year[-length(Count)]) abline(h=1, lty=3) ################################################### ### chunk number 31: ################################################### dev.print(pdf, "AllSS.pdf") ################################################### ### chunk number 32: ################################################### years <- 50 ################################################### ### chunk number 33: ################################################### set.seed(3) sim.Rs <- sample(x=obs.R, size=years, replace = TRUE) ################################################### ### chunk number 34: ################################################### output <- numeric(years+1) ################################################### ### chunk number 35: ################################################### output[1] <- Count[Year==max(Year)] ################################################### ### chunk number 36: ################################################### for( t in 1:years ) output[t+1] <- {output[t] * sim.Rs[t]} ################################################### ### chunk number 37: ################################################### myQ() ################################################### ### chunk number 38: ################################################### plot(0:years, output, type="l") ################################################### ### chunk number 39: ################################################### dev.print(pdf, "OneSim.pdf") ################################################### ### chunk number 40: ################################################### sims = 10 sim.RM <- matrix( sample(obs.R, sims * years, replace=TRUE), nrow=years, ncol=sims) ################################################### ### chunk number 41: ################################################### output[1] <- Count[Year==max(Year)] outmat <- sapply(1:sims, function(i) { for( t in 1:years ) output[t+1] <- output[t] * sim.RM[t,i] output} ) ################################################### ### chunk number 42: ################################################### matplot(0:years, outmat, type="l", log="y") ################################################### ### chunk number 43: ################################################### dev.print(pdf, "TenSim.pdf") ################################################### ### chunk number 44: ################################################### PopSim <- function(Rs, N0, years=50, sims=10) { sim.RM = matrix( sample(Rs, size=sims * years, replace =TRUE), nrow=years, ncol=sims) output <- numeric(years+1) output[1] <- N0 outmat <- sapply(1:sims, function(i) { for( t in 1:years ) output[t+1] <- round( output[t] * sim.RM[t,i], 0 ) output} ) return(outmat) } ################################################### ### chunk number 45: ################################################### system.time( output <- PopSim(Rs=obs.R, N0=43, sims=1000) ) ################################################### ### chunk number 46: ################################################### dim(output) ################################################### ### chunk number 47: ################################################### N.2053 <- output[51,] summary(N.2053, digits=6) ################################################### ### chunk number 48: ################################################### quantile(N.2053, prob=c(0.0275, .975) ) ################################################### ### chunk number 49: ################################################### quartz(,7,4) par(mfrow=c(1,2)) ################################################### ### chunk number 50: ################################################### hist(N.2053, main="N") hist(log10(N.2053+1), main="log(N+1)" ) abline(v=log10(quantile(N.2053, prob=c(0.0275, .975) )+1), lty=3) ################################################### ### chunk number 51: ################################################### dev.print(pdf, "histsim.pdf") ################################################### ### chunk number 52: ################################################### logOR <- log(obs.R) n <- length(logOR) t.quantiles <- qt( c(0.025, 0.975), df=n-1) ################################################### ### chunk number 53: ################################################### se <- sqrt( var(logOR) / n ) CLs95 <- mean( logOR ) + t.quantiles * se ################################################### ### chunk number 54: ################################################### R.limits <- exp(CLs95) R.limits ################################################### ### chunk number 55: ################################################### N.Final.95 <- Count[Year==max(Year)]*R.limits^50 round(N.Final.95) ################################################### ### chunk number 56: ################################################### myQ() ################################################### ### chunk number 57: ################################################### qqplot(qt(ppoints(n), df=n-1), scale(logOR)) qqline(scale(logOR)) ################################################### ### chunk number 58: ################################################### dev.print(pdf, "tCheck.pdf")