################################################### ### chunk number 1: ################################################### setwd("~/Documents/Projects/Book/draft/Chap06") 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)) #source("DMBpplane.R") library(deSolve); library(lattice) ################################################### ### chunk number 2: ################################################### LH.df <- read.table("LynxHare.txt", header=FALSE) names(LH.df) <- c("Year", "Hare", "Lynx") summary(LH.df) quartz(width=11/2.54, height=5/2.54); par(cex=.6) matplot(LH.df[["Year"]], LH.df[,2:3], type="l", ylab="Hare and Lynx in Thousands", xlab="Year", main="The Hudson Bay Company Pelt Cycle") legend(1900, 150, c("Snowshoe Hare", "Lynx"), lty=1:2, bty="n") dev.print(pdf, "LH.pdf") ################################################### ### chunk number 3: ################################################### predpreyLV <- function(t, y, params) { H <- y[1]; P <- y[2] with(as.list(params), { dH.dt <- b*H - a*P*H dP.dt <- e*a*P*H - s*P return( list(c(dH.dt, dP.dt)) ) }) } ################################################### ### chunk number 4: FR ################################################### a <- .1; w <- .1; D <- w/a; curve(a*x, 0,2, xlab="Prey Density", ylab="Prey Killed Per Predator") curve(w*x/(D+x), 0,2, add=TRUE, lty=2) curve(w*x^2/(D^2+x^2), 0,2, add=TRUE, lty=3) ################################################### ### chunk number 5: FRH ################################################### curve(w*x^2/(D^2+x^2)/x, 0, 2, ylim=c(0,a), lty=3, xlab="Prey Density", ylab="Prey Killed Per Predator Per Prey") curve(w*x/(D+x)/x, 0,2, lty=2, add=TRUE) curve(a*x/x, 0,2, add=TRUE, lty=1) legend('right', c("Type I", "Type II", "Type III"), lty=1:3, bty='n') ################################################### ### chunk number 6: ################################################### b <- 0.5; a <- 0.01; (Hi <- b/a) e <- 0.1; s <- 0.2 ; (Pi <- s/(e*a) ) ################################################### ### chunk number 7: LVPPiso ################################################### plot(c(0,2*Pi), c(0, 2 * Hi), type = "n", xlab = "H", ylab = "P") abline(h = Hi) text(Pi, Hi, "Prey isocline", adj=c(1.3,-.3)) arrows(x0=c(.85*Pi,1.15*Pi), y0=c(.3*Hi, 1.7*Hi), x1=c(1.15*Pi,.85*Pi), y1=c(.3*Hi, 1.7*Hi), len=.2) abline(v = Pi, lty=2) text(Pi, Hi, "Predator isocline", adj=c(1.1,-.2), srt=90 ) arrows(x0=c(.3*Pi, 1.7*Pi), y0=c(1.15*Hi, .85*Hi), x1=c(.3*Pi, 1.7*Pi), y1=c(.85*Hi, 1.15*Hi), lty=2, len=.2) text(x=c(2*Pi, 0,0, 2*Pi), y=c(2*Hi,2*Hi,0,0), 1:4, cex=2) ################################################### ### chunk number 8: ################################################### Jac <- matrix( c(0, -s/e, e*b, 0), byrow=TRUE, nr=2) eigen(Jac)[["values"]] ################################################### ### chunk number 9: ################################################### params1 <- c(b=b, a=a, s=s, e=e) Time <- seq(0,100, by=.1) # Set time here LV.out <- ode(c(H0=25,P0=5), Time, predpreyLV, params1) ################################################### ### chunk number 10: LVTimeDyn ################################################### matplot(Time, (LV.out[,2:3]), type="l", ylab="Population Size") ################################################### ### chunk number 11: ################################################### LV.out2 <- ode(c(H0=500,P0=15), Time, predpreyLV, params1) LV.out3 <- ode(c(H0=300,P0=50), Time, predpreyLV, params1) ################################################### ### chunk number 12: LVPhase ################################################### plot(LV.out[,2], LV.out[,3], type='l', ylab="P", xlab="H") points(25,5, cex=1.5, pch=19) arrows(x0=c(1300, -20, 500), y0=c(125, 175, 0), x1=c(650, -20, 950), y1=c(200, 100, 2), length=.1) lines(LV.out2[,2], LV.out2[,3]) points(500,15, cex=1.5) lines(LV.out3[,2], LV.out3[,3]) points(300,50, cex=1.5, pch=2) abline(h=b/a, lty=2); abline(v=s/(e*a), lty=2) ################################################### ### chunk number 13: RM model ################################################### predpreyRM <- function(t, y, p) { H <- y[1]; P <- y[2] with(as.list(p), { dH.dt <- b*H*(1-alpha*H) - w*P*H / (D + H) dP.dt <- e*w*P*H / (D + H) - s*P return( list(c(dH.dt, dP.dt)) ) }) } ################################################### ### chunk number 14: ################################################### b<-.8; e<-0.07; s<-.2; w <- 5; D <- 400; alpha <- .001 H <- 0:(1/alpha) Hiso <- expression(b/w*(D + (1-alpha * D)*H - alpha* H^2)) HisoStable <- eval(Hiso) ################################################### ### chunk number 15: rmtraj1 ################################################### p.RM <- c(b=b, alpha=alpha, e=e, s=s, w=w, D=D) Time <- 150 RM1 <- ode(c(900,120), 1:Time, predpreyRM, p.RM) ################################################### ### chunk number 16: RMisoclines ################################################### plot(H, HisoStable, type="l", ylab="P", xlab="H", ylim=c(0, 150)) abline(v=s*D/(e*w-s), lty=2) arrows(RM1[-Time,2],RM1[-Time,3], RM1[-1, 2], RM1[-1,3], length=.1) ################################################### ### chunk number 17: ################################################### dhdt <- expression(b*H*(1-alpha*H) - w*P*H/(D+H)) dpdt <- expression(e*w*P*H/(D+H) -s*P) ################################################### ### chunk number 18: ################################################### RMjac1 <- list( D(dhdt, 'H'), D(dpdt, 'H'), D(dhdt, 'P'), D(dpdt, 'P') ) ################################################### ### chunk number 19: ################################################### H <- s*D/(e*w-s) ################################################### ### chunk number 20: ################################################### P <- eval(Hiso) ################################################### ### chunk number 21: ################################################### ( RM.jac2 <- matrix(sapply(RMjac1, function(pd) eval(pd)), nrow=2) ) eigen(RM.jac2)[["values"]] ################################################### ### chunk number 22: RM.jac ################################################### RM.jac2 ################################################### ### chunk number 23: RM.eig ################################################### eigen(RM.jac2)[["values"]] ################################################### ### chunk number 24: ################################################### p.RM["alpha"] <- alpha <- .0005; Time <- 100 H <- 0:(1/alpha) RM2 <- ode(c(500,110), 1:Time, predpreyRM, p.RM) ################################################### ### chunk number 25: ParadoxPP ################################################### plot(H, eval(Hiso), type="l", ylab="P", xlab="H", ylim=c(0, max(RM2[,3]) ) ) lines(0:1000, HisoStable, lty=3) abline(v=s*D/(e*w-s), lty=2) arrows(RM2[-Time,2],RM2[-Time, 3], RM2[-1, 2], RM2[-1, 3], length=.1) ################################################### ### chunk number 26: ################################################### H <- with(as.list(p.RM), s*D/(e*w-s)) alphas <- 1/seq(H, 3000, by=10) ################################################### ### chunk number 27: ################################################### RM.jacList <- lapply(1:length(alphas), function(i) { alpha <- alphas[i]; P <- eval(Hiso) matrix(sapply(RMjac1, function(pd) eval(pd)), nrow=2) } ) ################################################### ### chunk number 28: ################################################### L1 <- sapply(RM.jacList, function(J) max(Re(eigen(J)[["values"]])) ) ################################################### ### chunk number 29: RMLambda ################################################### plot(1/alphas, L1, type='l', xlab=quote(italic(K)==1/alpha), ylab=quote(lambda[1])) abline(h=0, lty=3) abline(v=H, lty=2) ################################################### ### chunk number 30: ################################################### time <- 20; R <- 3; a <- .005 HPs <- data.frame(Hs <- numeric(time), Ps <- numeric(time)) ################################################### ### chunk number 31: ################################################### Pst <- log(R)/a; Hst <- Pst*R/(R-1) HPs[1, ] <- c(Hst+1, Pst) ################################################### ### chunk number 32: ################################################### for(t in 1:(time-1)) HPs[t+1,] <- {H <- R * HPs[t,1] * exp( -a * HPs[t,2] ) P <- HPs[t,1]*(1-exp(-a*HPs[t,2])) c(H,P) } ################################################### ### chunk number 33: NBdyn ################################################### plot(HPs[,1], HPs[,2], type="n", xlab="H", ylab="P"); points(Hst, Pst); arrows(HPs[-time,1], HPs[-time,2], HPs[-1,1], HPs[-1,2], length=.05) ################################################### ### chunk number 34: ################################################### totals <- c(1066, 176, 48, 8, 5) ################################################### ### chunk number 35: ################################################### obs <- rep(0:4, totals) ################################################### ### chunk number 36: ################################################### H <- sum( totals ) No.attacks <- sum(obs) mean.attacks <- mean(obs) ################################################### ### chunk number 37: ################################################### attack.events <- sample(x=1:H, size=No.attacks, replace=TRUE) ################################################### ### chunk number 38: ################################################### N.attacks.per.host <- table(attack.events) ################################################### ### chunk number 39: ################################################### (tmp <- table( N.attacks.per.host ) ) ################################################### ### chunk number 40: ################################################### not.att <- H - sum(tmp) ################################################### ### chunk number 41: ################################################### obs.sim <- rep(0:max(N.attacks.per.host), c(not.att, tmp)) table(obs) table(obs.sim) ################################################### ### chunk number 42: ################################################### n <- 1000 unatt.sim<- sapply(1:n, function(j) { ## Repeat n times... ## simulate the same number of attacks. host.sim.att <- sample(x=1:H, size=No.attacks, replace=TRUE) ## subtract the number of attacked hosts from the total in the host population ## The number of attacks on each attacked host attacks.per <- table(host.sim.att) ## the number of attacked hosts n.attd <- length( attacks.per) ## the number not attacked H - n.attd } ) ################################################### ### chunk number 43: PoisSim ################################################### hist(unatt.sim, xlab="Simulated # of Unattacked Hosts", prob=TRUE, xlim=c(1000, 1070) ) abline(v=1066, lty=3) abline(v=exp(-mean.attacks), lty=2) ################################################### ### chunk number 44: ################################################### (I <- var(obs)/mean(obs)) ################################################### ### chunk number 45: ################################################### 1-pchisq(I*H, df=H-1) ################################################### ### chunk number 46: ################################################### getk <- function(CV2, mu){mu/(mu*CV2-1)} getk(7.33, -.24) getk(7.77, -.46) getk(2, -.025) getk(.1,.395) ################################################### ### chunk number 47: nbinom ################################################### nb.dat <- cbind( Random = dnbinom(0:5, mu=.5, size=10^10), "k=10"=dnbinom(0:5, mu=.5, size=10), "k=1"=dnbinom(0:5, mu=.5, size=1), "k=0.01"=dnbinom(0:5, mu=.5, size=.01) ) matplot(0:5, nb.dat , type="b", pch=1:4, col=1, ylim=c(0,1), xlab="Attacks per Hosts", ylab="Probability") legend("topright", rev(colnames(nb.dat)), pch=4:1, lty=4:1, bty='n') mtext(quote(mu== 0.5), padj=2) ################################################### ### chunk number 48: discstab ################################################### F.H <- expression(R*H*(1+a*P/k)^-k - H ) F.P <- expression(H - H*(1+a*P/k)^-k - P) F.H.H <- D(F.H, "H"); F.H.P <- D(F.H, "P") F.P.H <- D(F.P, "H"); F.P.P <- D(F.P, "P") ################################################### ### chunk number 49: ################################################### k <- 10^seq(-1, 1, by=.01) R <- 3 a <- .005 HPeigs <- sapply(k, function(ki) { k <- ki P <- k*(R^(1/k) - 1)/a H <- P * R/(R-1) jac <- matrix( c(eval(F.H.H), eval(F.H.P), eval(F.P.H), eval(F.P.P)), nrow=2, byrow=TRUE ) eigen(jac)[["values"]] } ) ################################################### ### chunk number 50: LK2 ################################################### modmaxs <- apply( HPeigs, 2, function(lambdas) {i <- which.max(Mod(lambdas)) sign(Re(lambdas[i])) * Mod(lambdas[i]) }) plot(k, modmaxs, type='l', ylab=quote("Stability "*(lambda[1])) ) abline(h=-1, lty=3) ################################################### ### chunk number 51: ################################################### th <- seq(-pi, pi, len=100) z <- exp(1i*th) ################################################### ### chunk number 52: mod ################################################### par(pty="s") plot(z, type="l") points(0,0, pch=3); points(HPeigs[,100]) arrows(x0=c(0,0), y0=c(0,0), x1=Re(HPeigs[,100]), y1=Im(HPeigs[,100])) ################################################### ### chunk number 53: ################################################### time <- 20; R <- 3; a <- .005; k <- .6 HP2s <- data.frame(Hs <- numeric(time), Ps <- numeric(time)) ################################################### ### chunk number 54: ################################################### P2st <- k*(R^(1/k) - 1)/a; H2st <- P2st * R/(R-1) P2st;H2st HP2s[1, ] <- c(1000,500) ################################################### ### chunk number 55: ################################################### for(t in 1:(time-1)) HP2s[t+1,] <- { H <- R * HP2s[t,1] * (1+a*HP2s[t,2]/k)^(-k) P <- HP2s[t,1] - HP2s[t,1] * (1+a*HP2s[t,2]/k)^(-k) c(H,P) } ################################################### ### chunk number 56: NBdynStable ################################################### plot(HP2s[,1], HP2s[,2], type="l", xlab="H", ylab="P"); points(H2st, P2st); arrows(HP2s[-time,1], HP2s[-time,2], HP2s[-1,1], HP2s[-1,2], length=.05) ################################################### ### chunk number 57: SIRmodel ################################################### SIR <- function(t, y, p) { {S <- y[1]; I <- y[2]; R <- y[3]} with( as.list(p), { dS.dt <- -B*I*S dI.dt <- B*I*S - g*I dR.dt <- g*I return( list(c(dS.dt, dI.dt, dR.dt)) ) } ) } ################################################### ### chunk number 58: ################################################### N <- 10^4; I <- R <- 1; S <- N - I - R parms <- c(B=.01, g=4) ################################################### ### chunk number 59: ################################################### months <- seq(0,3, by=0.01) require(deSolve) SIR.out <- data.frame( ode(c(S,I,R), months, SIR, parms) ) ################################################### ### chunk number 60: SIR1 ################################################### matplot(months, SIR.out[,-1], type='l', lty=1:3, col=1) legend('right', c('R', 'I', 'S'), lty=3:1, col=3:1, bty='n') ################################################### ### chunk number 61: SIRmodel ################################################### SIRf <- function(t, y, p) { {S <- y[1]; I <- y[2]; R <- y[3]; N<- S+I+R} with( as.list(p), { dS.dt <- -B*I*S/N dI.dt <- B*I*S/N - g*I dR.dt <- g*I return( list(c(dS.dt, dI.dt, dR.dt)) ) } ) } ################################################### ### chunk number 62: ################################################### R <- 0; S <- I <- 1000; Ss <- Is <- seq(1, S, length=11); N <- S+I+R betaD <- 0.1; betaF <- betaD*N ################################################### ### chunk number 63: ################################################### mat1 <- sapply(Is, function(i) betaD * i * Ss) mat2 <- sapply(Is, function(i) betaF * i * Ss / (i + Ss + R) ) ################################################### ### chunk number 64: TransFuncs ################################################### layout(matrix(1:2, nr=1)) persp(mat1, theta=20, phi=15, r=10, zlim=c(0,betaD*S*I), main="Density Dependent", xlab="I", ylab="S", zlab="Transmission Rate") persp(mat2, theta=20, phi=15, r=10, zlim=c(0,betaF*S*I/N), main="Frequency Dependent", xlab="I", ylab="S", zlab="Transmission Rate") ################################################### ### chunk number 65: ################################################### S <- 4^(0:4); I <- 1; parmsf <- c(B=1, g=0); parmsd <- c(B=1/16, g=0) ################################################### ### chunk number 66: ################################################### Months <- seq(0, 8, by=0.1) outd <- sapply(S, function(s) {out <- ode(c(s,I,R), Months, SIR, parmsd) out[,3]/apply(out[,2:4], 1, sum) } ) outf <- sapply(S, function(s) {out <- ode(c(s,I,R), Months, SIRf, parmsf) out[,3]/apply(out[,2:4], 1, sum) } ) #TR <- sapply(S, function(s) {R <- s/2; parmsf["B"]*s*I/(s+I+R)}) ################################################### ### chunk number 67: SIRdd ################################################### matplot(Months, outd, type='l', col=1, ylab="Prevalence (I/N)") ################################################### ### chunk number 68: SIRfd ################################################### matplot(Months, outf, type='l', col=1, ylab="Prevalence (I/N)") legend('bottomright', legend=S, lty=1:length(S), bty='n') #plot(S, TR) ################################################### ### chunk number 69: ################################################### SIRbd <- function(t, y, p) { S <- y[1]; I <- y[2]; R <- y[3] with( as.list(p), { dS.dt <- b*(S+I+R) - B*I*S - m*S dI.dt <- B*I*S - g*I - m*I dR.dt <- g*I - m*R return( list(c(dS.dt, dI.dt, dR.dt)) ) } ) } ################################################### ### chunk number 70: ################################################### N <- 10^6 R <- 0; I <- 1; S <- N - I - R ################################################### ### chunk number 71: ################################################### g <- 1/(13/365) ################################################### ### chunk number 72: ################################################### b <- 1/50 ################################################### ### chunk number 73: ################################################### age <- 5 R0 <- 1 + 1/(b*age) ################################################### ### chunk number 74: ################################################### B <- R0 * (g + b) / N ################################################### ### chunk number 75: ################################################### parms <- c(B = B, g = g, b = b, m=b) years <- seq(0,30, by=.1) ################################################### ### chunk number 76: ################################################### SIRbd.out <- data.frame(ode(c(S=S,I=I,R=R), years, SIRbd, parms, hmax=.01)) ################################################### ### chunk number 77: SIR2 ################################################### matplot(SIRbd.out[,1], sqrt(SIRbd.out[,-1]), type='l', col=1, lty=1:3, ylab="sqrt(No. of Individuals)", xlab='Years') legend('right', c('S','I','R'), lty=1:3, bty='n') ################################################### ### chunk number 78: Ross ################################################### data(ross) plot(CumulativeDeaths ~ Week, data=ross) ################################################### ### chunk number 79: ################################################### sirLL=function(logit.B, logit.g, log.N, log.I0) { parms <- c(B=plogis(logit.B), g=plogis(logit.g)); x0 <- c(S=exp(log.N), I=exp(log.I0), R=0) Rs <- ode(y=x0, ross$Week, SIR, parms, hmax=.01)[,4] SD <- sqrt(sum( (ross$CumulativeDeaths -Rs)^2)/length(ross$Week) ) -sum(dnorm(ross$CumulativeDeaths, mean=Rs, sd=SD, log=TRUE)) } ################################################### ### chunk number 80: ################################################### load("rossfits.Rdata") ################################################### ### chunk number 81: eval=FALSE ################################################### ## require(bbmle) ## fit <- mle2(sirLL, start=list(logit.B=qlogis(1e-5), logit.g=qlogis(.2), ## log.N=log(1e6), log.I0=log(1) ), method="Nelder-Mead") ################################################### ### chunk number 82: sumfit ################################################### summary(fit) ################################################### ### chunk number 83: eval=FALSE ################################################### ## fit2 <- mle2(sirLL, start=as.list(coef(fit)), ## fixed=list(log.N=coef(fit)[3], log.I0=coef(fit)[4]), method="Nelder-Mead") ################################################### ### chunk number 84: sumfit2 ################################################### summary(fit2) ################################################### ### chunk number 85: eval=FALSE ################################################### ## pr2 <- profile(fit2) ################################################### ### chunk number 86: profileSIR ################################################### par(mar=c(5,4,4,1)) plot(pr2) ################################################### ### chunk number 87: eval=FALSE ################################################### ## save(fit,fit2,pr2, file="rossfits.Rdata") ################################################### ### chunk number 88: ################################################### p <- as.numeric(c(plogis(coef(fit2)[1:2]), exp(coef(fit2)[3:4])) ); p ################################################### ### chunk number 89: ################################################### inits <- c(S = p[3], I=p[4], R=0) params <- c(B=p[1], g=p[2]) SIR.Bombay <- data.frame( ode(inits, ross$Week, SIR, params) ) ################################################### ### chunk number 90: FittedBombay ################################################### matplot(SIR.Bombay[,1], SIR.Bombay[,3:4], type='l', col=1) points(ross$Week, ross$CumulativeDeaths) legend('topleft', c("I", "R"), lty=1:2, bty='n') ################################################### ### chunk number 91: ################################################### (CIs <- confint(pr2)) ################################################### ### chunk number 92: ################################################### (gs <- as.numeric(plogis(CIs[2,]) )) ################################################### ### chunk number 93: ################################################### 7*1/gs