################################################### ### chunk number 1: ################################################### setwd("~/Documents/Projects/Book/draft/Chap03") 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=75, 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: ################################################### data(sparrows) attach(sparrows) obs.Rs <- Count[-1]/Count[-length(Count)] ################################################### ### chunk number 3: MelospizaDDG1 ################################################### plot(Year,Count, type="b") ################################################### ### chunk number 4: MelospizaDDG2 ################################################### plot(obs.Rs ~ Count[-length(Count)], ylab=expression("N"[t+1]/"N"[t]), xlab="Sparrow Count", type="n") text(Count[-length(Count)], obs.Rs, substr(as.character(Year), 3,4 ) ) abline(mod <- lm(obs.Rs ~ Count[-length(Count)] ), lty=2 ) #newd <- 20:120 #cts <- Count[-length(Count)] #mod <- loess(obs.Rs ~ cts ) #lines(newd, predict( mod, newdata=data.frame(cts=newd) )) ################################################### ### chunk number 5: ################################################### dlogistic <- function(alpha=0.01, rd=1, N0=2, t=15) { ## make a holder for initial N and all subsequent N. N <- c(N0, numeric(t)) ## Use a for loop to create the population dynamic. We ## CANNOT use lapply or sapply, because each N depends on ## the previous N. for(i in 1:t) N[i+1] <- {N[i] + rd * N[i] * (1 - alpha * N[i])} return(N) } ################################################### ### chunk number 6: ################################################### Nts <- dlogistic() ################################################### ### chunk number 7: DLG1 ################################################### t <- 15 a <- 0.01 plot(0:t, Nts) abline(h=1/a, lty = 3) ################################################### ### chunk number 8: ################################################### total.incr <- Nts[1:t+1] - Nts[1:t] per.capita.incr <- total.incr/Nts[1:t] ################################################### ### chunk number 9: PGI ################################################### plot(Nts[1:t], total.incr) ################################################### ### chunk number 10: PPGI ################################################### plot(Nts[1:t], per.capita.incr) ################################################### ### chunk number 11: ################################################### N0s <- c(0, runif(30) * 1.1 * 1/a) N <- sapply(N0s, function(n) dlogistic(N0=n) ) ################################################### ### chunk number 12: DLGinit ################################################### matplot(0:t,N, type="l", lty=1, lwd=.75, col=1) text(t, 1/a, expression(italic("K") == 1/alpha), adj=c(1,0)) ################################################### ### chunk number 13: ################################################### a.s <- 1/runif(30, min=50, max=1000) N <- sapply(a.s, function(a) dlogistic(alpha=a, t=15) ) ################################################### ### chunk number 14: DLGalpha ################################################### matplot(0:t, N, type="l", ylim=c(0,1000), lty=1, lwd=.75, col=1) text(8, 1/min(a.s), bquote(italic(alpha) == .(round(min(a.s),3))), adj=c(1,0.5)) text(10, 1/max(a.s), bquote(italic(alpha) == .(round(max(a.s),3))), adj=c(0,1.2)) ################################################### ### chunk number 15: ################################################### rd.v <- seq(1.3, 2.8, by=.3) t <- 15 Ns <- data.frame(sapply(rd.v, function(r) dlogistic(rd=r, t=t) )) ################################################### ### chunk number 16: DLGrd ################################################### matplot(0:t, Ns, type="l", col=1) ################################################### ### chunk number 17: ################################################### tmp <- data.frame(rd=as.factor(rd.v), t(Ns)) ################################################### ### chunk number 18: ################################################### Ns2 <- reshape(tmp, varying=list(2:ncol(tmp)), idvar="rd", v.names="N", direction="long") str(Ns2) ################################################### ### chunk number 19: DDGr2 ################################################### library(lattice) print( xyplot(N ~ time|rd, data=Ns2, type="l", layout=c(3,2,1), col=1 ) ) ################################################### ### chunk number 20: ################################################### num.rd <- 201 rd.s <- seq(1,3, length=num.rd) t <- 400 ################################################### ### chunk number 21: eval=FALSE ################################################### ## num.rd <- 1001 ## rd.s <- seq(1,3, length=num.rd) ## t <- 1000 ################################################### ### chunk number 22: ################################################### tmp <- sapply(rd.s, function(r) dlogistic(rd=r, N0=99, t=t) ) ################################################### ### chunk number 23: ################################################### tmp.s <- stack( as.data.frame(tmp) ) names(tmp.s) <- c("N", "Old.Column.ID") tmp.s$rd <- rep(rd.s, each=t+1) tmp.s$time <- rep(0:t, num.rd) ################################################### ### chunk number 24: Bifurcation ################################################### N.bif <- subset(tmp.s, time > 0.5 * t ) plot(N ~ rd, data = N.bif, pch=".", xlab=quote("r"["d"])) ################################################### ### chunk number 25: ################################################### N.init <- c(97,98,99) t <- 30 Ns <- sapply(N.init, function(n0) dlogistic(rd=2.7, N0=n0, t=t) ) ################################################### ### chunk number 26: DDGChaosInitN ################################################### matplot(0:t,Ns, type="l", col=1) ################################################### ### chunk number 27: ################################################### B.N <- expression(-a*N^2 + e*N - f) D.N <- expression(g*N^2) ################################################### ### chunk number 28: ################################################### a <-1/1600; e <- 80*a; f <- .2; g <- 0.00004 N <- 0:100 ################################################### ### chunk number 29: DD1 ################################################### plot(N, eval(B.N), type="l", ylab="B(N), D(N)") lines(N, eval(D.N), lty=2) abline(h=0, lty=3) legend("bottom", c("Effect on Birth Rate", "Effect on Death Rate"), lty=1:2, bty='n') ################################################### ### chunk number 30: ################################################### myQ() ################################################### ### chunk number 31: ################################################### plot(N, eval(B.N) + eval(D.N), type='l', ylim=c(-0.4,1), ylab="Net Density Dependence, F(N)") abline(h=0, lty=3) ################################################### ### chunk number 32: ################################################### curve(1-0.01*x, 0,100, lty=2, add=T) legend('topright', c("Generalized", "Logistic"), lty=1:2, bty='n') ################################################### ### chunk number 33: ################################################### dev.print(pdf, "DD2.pdf") ################################################### ### chunk number 34: ################################################### myQ(w=6) ################################################### ### chunk number 35: ################################################### pop.growth.rate <- expression( r * N * (1 - alpha * N) ) r <- 1; alpha <- 0.01; N <- 0:120 ################################################### ### chunk number 36: ################################################### plot(N, eval(pop.growth.rate), type="l", ylab="Population Growth Rate (dN/dt)", xlab="N") abline(h=0) legend('topright', "r=1", lty=1) ################################################### ### chunk number 37: ################################################### N <- c(0, 10, 50, 100, 115) points(N, eval(pop.growth.rate), cex=1.5) text(N, eval(pop.growth.rate), letters[1:5],adj=c(.5,2)) ################################################### ### chunk number 38: ################################################### arrows(20,2, 80,2, length=.1, lwd=3) arrows(122,-2,109,-2, length=.1, lwd=3) ################################################### ### chunk number 39: ################################################### dev.print(pdf, "Stability1.pdf") ################################################### ### chunk number 40: Stability2 ################################################### mi <- 99.4; ma <- 100.6 N <- seq(mi,ma, length=30) plot(N, eval(pop.growth.rate), type="l", ylab="Population Growth Rate (dN/dt)", xlab="N") abline(h=0) r <- .5 lines(N, eval(pop.growth.rate), lty=2); r <- 1 segments(c(100-.3, 100+.3), c(-.05,-.05), c(100-.3,100+.3), c(.05,.05)) text(99.7,-.15, quote("N*"-"x")) text(100.3,.15, quote("N*"+"x")) points(100, 0) text(100, -.1, "N*") legend('topright', c("r=1", "r=0.5"), lty=1:2) ################################################### ### chunk number 41: ################################################### myQ( mr = c( 1, 3, .25, 1.25) ) par( mgp = c( 2, .75, 0 ) ) x <- 1 plot(c(-1,8), 1.2*c(-x,x), type='n', ylab="N", xlab='', axes=FALSE) #box() axis(2, at=c(-x,0,x), labels=c(expression("N*"-"x"),"N*", expression("N*"+"x")) ) abline(h=0) t <- 1 points(c(0, t), c(1, x*exp(-1*t) ), pch=c(19,1)) curve(1*exp(-1*x), 0, .98, add=T, lty=3) arrows(0,-x-.2, 7,-x-.2) mtext("Time", side=1) text(0, 0, "Perturbation", srt=90, adj=c(-.1,-.1) ) arrows(0,0,0,.8, length=.1) text(3, .7, expression("Recovery at rate"=="e"^"-r")) points(c(0+2, t+2), c(-1, -x*exp(-1*t) ), pch=c(19,1)) x2 <- seq(0,.98, length=30) y <- 0 - 1 * exp(-1*x2) lines(x2+2, y, lty=3) text(2, 0, "Perturbation", srt=90, adj=c(1.02,-.1) ) arrows(2, 0, 2,-.8, length=.1) text(5, -.7, expression("Recovery at rate"=="e"^"-r")) ################################################### ### chunk number 42: ################################################### dev.print(pdf, "Stability3.pdf") ################################################### ### chunk number 43: ################################################### dF.dN <- deriv(pop.growth.rate,"N") N <- c(0, 1/alpha) eval(dF.dN) ################################################### ### chunk number 44: ################################################### clogistic <- function(times, y, parms) { n <- y[1]; r <- parms[1]; alpha <- parms[2] dN.dt <- r * n * (1 - alpha * n) return( list( c(dN.dt) ) ) } ################################################### ### chunk number 45: ################################################### prms <- c(r=1, alpha=.01); init.N <- c(1) t.s <- seq(0.1, 10, by=.1) ################################################### ### chunk number 46: ################################################### library(deSolve) out <- ode(y=init.N, times=t.s, clogistic, parms=prms) ################################################### ### chunk number 47: ContDyn1 ################################################### plot(out[,1], out[,2], type='l', xlab="Time", ylab="N") ################################################### ### chunk number 48: ################################################### outmat <- matrix(NA, nrow=length(t.s), ncol=20) for(j in 1:20) outmat[,j] <- { y <- runif(n=1, min=0, max=120) prms <- c(r=runif(1, .01,2), alpha=0.01 ) ode(y, times=t.s, clogistic, prms)[,2] } ################################################### ### chunk number 49: ContDyn2 ################################################### matplot(t.s, outmat, type='l', col=1, ylab="All Populations") ################################################### ### chunk number 50: ################################################### thetalogistic <- function(times, y, parms) { n <- y[1] with(as.list(parms), { dN.dt <- r * n * (1-(alpha*n)^theta) return( list( c(dN.dt) ) ) } ) } ################################################### ### chunk number 51: ThetaDD ################################################### r <- .75; alpha <- 0.01; theta <- c(.5, 1, 2); N <- 0:110 theta.out <- sapply(theta, function(th) {1-(alpha*N)^th}) matplot(N, theta.out, type='l', col=1) abline(h=0) legend('topright', legend=paste("theta =", c(2, 1, 0.5) ), lty=3:1) ################################################### ### chunk number 52: ################################################### myQ() ################################################### ### chunk number 53: ################################################### thetaGR.out <- sapply(theta, function(th) {r*N*(1-(alpha*N)^th)}) matplot(N, thetaGR.out, type='l', col=1) abline(h=0) legend('bottomleft', legend=paste("theta =", c(2, 1, 0.5) ), lty=3:1) ################################################### ### chunk number 54: eval=FALSE ################################################### ## r <- 3.4 ## lines(N, r*N*(1-(alpha*N)^theta[1])) ## text(12, 24, "r=3.4") ################################################### ### chunk number 55: ################################################### dev.print(pdf, "ThetaGR.pdf") ################################################### ### chunk number 56: ThetaN ################################################### prms <- c(r <- 0.75, alpha <- 0.01, theta=1) thetaN <- sapply(theta, function(th) {prms["theta"] <- th ode(y=1, t.s, thetalogistic, prms)[,2] } ) matplot(t.s, thetaN, type='l') legend('topleft', legend=paste("theta =", c(2, 1, 0.5) ), lty=3:1, bty='n') ################################################### ### chunk number 57: ################################################### myQ() ################################################### ### chunk number 58: ################################################### r <- .5; alpha <- 0.01 N <- 0:105 plot(N, eval(pop.growth.rate), type="l", ylim=c(-3, 15), ylab="dN/dt and FN") abline(h=0) ################################################### ### chunk number 59: ################################################### F <- r/2 abline(c(0,F), lty=2) ################################################### ### chunk number 60: ################################################### dev.print(device=pdf, file="MSY.pdf") ################################################### ### chunk number 61: MSY2 ################################################### pgr.H <- expression(r * N * (1 - alpha * N)-F*N) N <- 0:55 plot(N, eval(pgr.H), type="l", ylab="dN/dt (growth - harvesting)") abline(h=0) ################################################### ### chunk number 62: eval=FALSE ################################################### ## library(nlme) ## library(lattice) ## ClostExp <- groupedData(No.per.ml ~ Day|ID, tmp) ################################################### ### chunk number 63: ################################################### library(nlme) library(lattice) data(ClostExp) summary(ClostExp) ################################################### ### chunk number 64: ################################################### trellis.device(width=8, height=6, color=FALSE) ################################################### ### chunk number 65: ################################################### xyplot(No.per.ml ~ Day|Nutrients, ClostExp, groups=rep, type='b', scales=list(relation="free"), auto.key=list(columns=4, lines=T) ) ################################################### ### chunk number 66: ################################################### dev.print(pdf, "ClostRaw.pdf") ################################################### ### chunk number 67: ################################################### subset(ClostExp, Nutrients=="high" & No.per.ml > 1000) subset(ClostExp, Nutrients=="low" & No.per.ml > 100) ################################################### ### chunk number 68: ################################################### Hi.c <- subset(ClostExp, Nutrients=="high" & rep=="c") ################################################### ### chunk number 69: ################################################### n <- nrow(Hi.c) N.change <- Hi.c$No.per.ml[-1]/Hi.c$No.per.ml[-n] interval <- diff(Hi.c$Day) pgr <- log(N.change)/interval ################################################### ### chunk number 70: ################################################### myQ() ################################################### ### chunk number 71: ################################################### Nt <- Hi.c$No.per.ml[-n] plot(pgr ~ Nt) mod1 <- lm(pgr ~ Nt) abline(mod1) ################################################### ### chunk number 72: ################################################### dev.print(pdf, "PGR.pdf") ################################################### ### chunk number 73: ################################################### summary(mod1) ################################################### ### chunk number 74: ################################################### EachPop <- lapply( split(ClostExp, list(ClostExp$Nutrients, ClostExp$rep)), function(X) { n <- nrow(X) N.change <- (X$No.per.ml[-1]/X$No.per.ml[-n]) interval <- diff(X$Day) data.frame(Nutrients=as.factor(X$Nutrients[-n]), rep = as.factor(X$rep[-n]), pgr = log(N.change)/interval, Nt = X$No.per.ml[-n]) }) ################################################### ### chunk number 75: ################################################### AllPops <- NULL for(i in 1:length(EachPop) ) AllPops <- rbind(AllPops, EachPop[[i]]) ################################################### ### chunk number 76: ################################################### trellis.device(width=8, height=6, color=FALSE) ################################################### ### chunk number 77: ################################################### xyplot(pgr ~ Nt|rep*Nutrients, AllPops, layout=c(4,2,1), #scales=list(x=list(rot=90, alternating=F, relation="free")), scales=list(x=list(rot=90)), panel=function(x,y){ panel.grid(h=-1, v=-4) panel.xyplot(x,y, type=c('p','r') ) }) ################################################### ### chunk number 78: ################################################### dev.print(pdf, "pgrall.pdf") ################################################### ### chunk number 79: ################################################### AllPops$ID <- with(AllPops, Nutrients:rep) modSlope <- lme(pgr ~ Nt + Nutrients + Nt:Nutrients, data=AllPops, random=~1|ID ) ################################################### ### chunk number 80: ################################################### trellis.device(width=4, height=4, color=FALSE) ################################################### ### chunk number 81: ################################################### xyplot(resid(modSlope) ~ fitted(modSlope)) ################################################### ### chunk number 82: ################################################### dev.print(pdf, "pgrres1.pdf") ################################################### ### chunk number 83: ################################################### trellis.device(width=8, height=3, color=FALSE) ################################################### ### chunk number 84: ################################################### qqmath( ~resid(modSlope)|ID, data=AllPops,layout=c(8,1,1) ) ################################################### ### chunk number 85: ################################################### dev.print(pdf, "pgrres2.pdf") ################################################### ### chunk number 86: ################################################### trellis.device(width=4, height=4, color=FALSE) ################################################### ### chunk number 87: ################################################### qqmath(~ranef(modSlope)) ################################################### ### chunk number 88: ################################################### dev.print(pdf, "pgrranqq.pdf") ################################################### ### chunk number 89: ################################################### modSlope2 <- update(modSlope, weights=varExp()) ################################################### ### chunk number 90: ################################################### anova(modSlope, modSlope2) ################################################### ### chunk number 91: ################################################### anova(modSlope2) ################################################### ### chunk number 92: ################################################### summary(modSlope2)$tTable ################################################### ### chunk number 93: ################################################### cfs <- fixef(modSlope2); cfs - cfs[2]/cfs[1] - (cfs[2]+cfs[4])/(cfs[1]+ cfs[3]) ################################################### ### chunk number 94: ################################################### ilogistic <- function(t, alpha, N0, r) { N0*exp(r*t) / (1 + alpha*N0*(exp(r*t) - 1)) } ################################################### ### chunk number 95: ################################################### myQ() ################################################### ### chunk number 96: ################################################### plot(No.per.ml ~ Day, ClostExp, subset=Nutrients=="high") curve(ilogistic(x, alpha=-cfs[2]/cfs[1], N0=6, r=cfs[1] ), 1, 60, add=T, lty=2) ################################################### ### chunk number 97: ################################################### dev.print(pdf, "hi.pdf") ################################################### ### chunk number 98: fit.list ################################################### Cmod.list <- nlsList(No.per.ml ~ ilogistic(Day, alpha, N0, r)|ID, data=ClostExp, start= c(alpha=1/100, N0=5, r=.15), control=list(maxiter=1000)) ################################################### ### chunk number 99: "q" ################################################### quartz(, 6,3) ################################################### ### chunk number 100: ################################################### plot(coef(Cmod.list), layout=c(3,1,1)) ################################################### ### chunk number 101: ################################################### dev.print(pdf, 'coefs.pdf') ################################################### ### chunk number 102: "resids.Comd.list" ################################################### plot(Cmod.list) ################################################### ### chunk number 103: ################################################### Cmod.all <- nlme( Cmod.list , random=pdDiag(form=alpha + N0 + r~1), weights=varPower() ) ################################################### ### chunk number 104: "cfs2" ################################################### cfs2 <- fixef(Cmod.all);cfs2 ################################################### ### chunk number 105: ################################################### Cmod.all2 <- update(Cmod.all, fixed = list( alpha ~ Nutrients, N0 + r ~1), start=c(cfs2[1], 0, cfs2[2], cfs2[3]) ) ################################################### ### chunk number 106: ################################################### summary(Cmod.all2) ################################################### ### chunk number 107: ################################################### Cmod.all3 <- update(Cmod.all2, random=N0~1) anova(Cmod.all2, Cmod.all3) ################################################### ### chunk number 108: ################################################### intervals(Cmod.all3)$fixed ################################################### ### chunk number 109: ################################################### 1/intervals(Cmod.all3)$fixed[1:2, c(1,3)] ################################################### ### chunk number 110: ################################################### trellis.device(width=8, height=6, color=FALSE) ################################################### ### chunk number 111: ################################################### plot(augPred(Cmod.all3, level=0:1), #scales=list(y=list(log=10)), layout=c(4,2,1), scales=list(y=list(relation="free"))) ################################################### ### chunk number 112: ################################################### dev.print(pdf, "dyns.pdf") ################################################### ### chunk number 113: eval=FALSE ################################################### ## prms <- c(r=log(2)/.5, alpha=10^(-7)); init.N <- c(1000) ## t.s <- seq(0.1, 40, by=.1) ## library(deSolve) ## out <- ode(y=init.N, times=t.s, clogistic, parms=prms) ## plot(out[,1], out[,2], type='l', xlab="Time", ylab="N")