################################################### ### chunk number 1: ################################################### library(emdbook) library(plotrix) ## for corner.loc ## try(detach("package:nlme")) library(bbmle) library(lattice) library(MASS) library(Hmisc,warn=FALSE) source("chapskel.R") ################################################### ### chunk number 2: ################################################### binomNLL1 = function(p,k,N) { -sum(dbinom(k,prob=p,size=N,log=TRUE)) } ################################################### ### chunk number 3: ################################################### data(ReedfrogPred) x = subset(ReedfrogPred,pred=="pred" & density==10 & size=="small") k = x$surv ################################################### ### chunk number 4: ################################################### op=par(mfrow=c(1,2),lwd=2,bty="l",las=1,cex=1.5) p.mle <- sum(k)/40 mlik <- exp(sum(dbinom(k,size=10,prob=p.mle,log=TRUE))) pvec <- seq(0,1,length=100) liks <- sapply(pvec,function(x)sum(dbinom(k,size=10,prob=x,log=TRUE))) plot(pvec,exp(liks),type="l",xlab="Predation probability\nper capita", ylab="Likelihood",log="y",ylim=c(1e-20,1),axes=FALSE) axis(side=1,at=seq(0,1,by=0.25)) ## if axTicks worked ... loc2 = seq(-20,0,by=5) axis(side=2,at=10^loc2, do.call("c",lapply(loc2, function(z)substitute(expression(10^x),list(x=z))))) abline(v=p.mle,lty=2) abline(h=mlik,lty=2) par(xpd=NA) text(p.mle,10^(par("usr")[4]+0.5), expression(hat(p)==0.75)) text(0,mlik*10,expression(L[max]==5.1 %*% 10^-4),adj=0) box() par(xpd=TRUE) plot(table(factor(k,levels=0:10))/4,xlab="# of successes", ylab="Probability") points(0:10,dbinom(0:10,size=10,prob=0.75),pch=16,col="darkgray") par(op) ################################################### ### chunk number 5: ################################################### ptry <- function() { x <- rbinom(4,prob=0.75,size=10) sum(dbinom(x,prob=0.75,size=10,log=TRUE)) } rprobs <- replicate(1000,ptry()) binom.pval <- sum(rprobs 0] warning("levels truncated to positive values only") } if (is.null(confstr)) { confstr <- paste(format(100 * pchisq(levels^2, 1)), "%", sep = "") } mlev <- max(levels) * 1.05 nm <- names(obj) ## opar <- par(mar = c(5, 4, 1, 1) + 0.1) for (i in seq(along = nm)) { ## This does not need to be monotonic ## cat("**",i,obj[[i]]$par.vals[,i],obj[[i]]$z,"\n") ## browser() yvals <- obj[[i]]$par.vals[,nm[i],drop=FALSE] sp <- splines::interpSpline(yvals, obj[[i]]$z, na.action=na.omit) bsp <-splines:: backSpline(sp) ## x0 <- predict(bsp,0)$y if (missing(xlim)) xlim <- predict(bsp, c(-mlev, mlev))$y if (is.na(xlim[1])) xlim[1] <- min(yvals) if (is.na(xlim[2])) xlim[2] <- max(yvals) if (sqrVal) { if (!add) { if (missing(ylim)) ylim <- c(0,mlev^2) plot(I(z^2) ~ par.vals[, i], data = obj[[i]], xlab = xlabs[i], ylim = ylim, xlim = xlim, ylab = ylab, type = "n", ...) } avals <- rbind(as.data.frame(predict(sp)), data.frame(x = drop(yvals), y = obj[[i]]$z)) avals$y <- avals$y^2 lines(avals[order(avals$x), ], col = col.prof, lty=lty.prof) } else { if (absVal) { if (!add) { if (missing(ylim)) ylim <- c(0,mlev) plot(abs(z) ~ par.vals[, i], data = obj[[i]], xlab = nm[i], ylim = ylim, xlim = xlim, ylab = expression(abs(z)), type = "n", ...) } avals <- rbind(as.data.frame(predict(sp)), data.frame(x = yvals, y = obj[[i]]$z)) avals$y <- abs(avals$y) lines(avals[order(avals$x), ], col = col.prof, lty=lty.prof) } else { if (!add) { if (missing(ylim)) ylim <- c(-mlev,mlev) plot(z ~ par.vals[, i], data = obj[[i]], xlab = nm[i], ylim = ylim, xlim = xlim, ylab = expression(z), type = "n", ...) } lines(predict(sp), col = col.prof, lty=lty.prof) } } abline(v = x0, h=0, col = col.minval, lty = lty.minval) for (i in 1:length(levels)) { lev <- levels[i] confstr.lev <- confstr[i] ## Note: predict may return NA if we didn't profile ## far enough in either direction. That's OK for the ## "h" part of the plot, but the horizontal line made ## with "l" disappears. pred <- predict(bsp, c(-lev, lev))$y ## horizontal if (absVal || sqrVal) levs=rep(lev,2) else levs=c(-lev,lev) if (sqrVal) lines(pred, levs^2, type = "h", col = col.conf, lty = lty.conf) else lines(pred, levs, type = "h", col = col.conf, lty = 2) ## vertical pred <- ifelse(is.na(pred), xlim, pred) if (sqrVal) lines(pred, rep(lev^2,2), type = "l", col = col.conf, lty = lty.conf) else if (absVal) lines(pred, rep(lev, 2), type = "l", col = col.conf, lty = lty.conf) else { lines(c(x0,pred[2]), rep(lev, 2), type = "l", col = col.conf, lty = lty.conf) lines(c(pred[1],x0), rep(-lev, 2), type = "l", col = col.conf, lty = lty.conf) } if (plot.confstr) { if (sqrVal) text(labels=confstr.lev,x=x0,y=lev^2,col=col.conf) else text(labels=confstr.lev,x=x0,y=lev^2,col=col.conf) } } } ## par(opar) } op=par(lwd=2,bty="l",las=1,cex=1.5) ## attach(ReedfrogFuncresp) sqrprofplot(p2,sqrVal=TRUE,axes=FALSE,conf=0.95,col.conf="gray", col.prof="black",col.minval=NA,ylim=c(0,40),xlim=c(0.3,0.75), xlab="Attack rate (a)", ylab="Negative log-likelihood") axis(side=1) ## convert scale (deviance) ## scale = 2*(L-Lmin) ## L = scale/2+Lmin ylocs <- c(47,seq(50,75,by=5)) axis(side=2,at=2*(ylocs+logLik(m2)),labels=ylocs) box() aslice = sapply(bsurf$x, function(a)binomNLL2(c(a,coef(m2)["h"]), N=ReedfrogFuncresp$Initial, k=ReedfrogFuncresp$Killed)) lines(bsurf$x,2*(aslice+logLik(m2)),lty=3,xpd=NA) points(coef(m2)["a"],0,pch=16) points(pt3["a"],2*(pt3[3]+as.numeric(logLik(m2))),pch=2) points(pt4["a"],2*(pt4[3]+as.numeric(logLik(m2))),pch=8) par(xpd=NA) text(c(0.36,0.72),c(28,11),c("slice","profile"),cex=1, adj=0) par(xpd=FALSE) ##detach(ReedfrogFuncresp) par(op) ################################################### ### chunk number 52: ################################################### op=par(lwd=2,bty="l",las=1,cex=1.5) ## frog data par(mar=c(5,4,2,4)+0.1) plot(ReedfrogFuncresp$Initial,ReedfrogFuncresp$Killed, xlab="Initial density",ylab="Number killed") tmpf <- function(pars,...) { with(as.list(pars),curve(a*x/(1+a*h*x),add=TRUE,...)) } tmpf(coef(m2)) tmpf(pt3[1:2],lty=2) tmpf(pt4[1:2],lty=3) avals <- c(coef(m2)["a"],pt3["a"],pt4["a"]) hvals <- c(coef(m2)["h"],pt3["h"],pt4["h"]) hts <- avals*104/(1+avals*hvals*100) points(rep(110,3),hts,pch=c(16,2,8),xpd=NA) text(rep(114,3),hts,c("MLE","profile","slice"),xpd=NA,adj=0) ################################################### ### chunk number 53: eval=FALSE ################################################### ## library(rgl) ## bsurf2 = bsurf ## bsurf2$z = pmin(bsurf2$z,55) ## lev <- -logLik(m2)+qchisq(0.95,c(1,2))/2 ## cL = contourLines(bsurf$x,bsurf$y,bsurf$z, ## levels = lev) ## ll <- -logLik(m2) ## persp3d(bsurf$x,bsurf$y,bsurf2$z,col="gray",alpha=0.8) ## grid3d(side="x-") ## grid3d(side="z-") ## grid3d(side="y-") ## ## add contours? ## lines3d(cL[[1]]$x,cL[[1]]$y,rep(lev[1],length(cL[[1]]$x)),col="red",size=2) ## lines3d(cL[[2]]$x,cL[[2]]$y,rep(lev[2],length(cL[[2]]$x)),col="blue",size=2) ## lines3d(prof@profile$a$par.vals[,"a"],prof@profile$a$par.vals[,"h"], ## ll+prof@profile$a$z^2/2,size=2,col="green") ## lines3d(prof@profile$h$par.vals[,"a"],prof@profile$h$par.vals[,"h"], ## ll+prof@profile$h$z^2/2,size=2,col="yellow") ## ## etc.: various bugs to work out ################################################### ### chunk number 54: ################################################### source("sqrprofplot.R") op=par(mfrow=c(1,2),lwd=2,bty="l",las=1,cex=1.5) par(mar=c(5,4,2,4)+0.1,mgp=c(2,1,0)) sqrprofplot(p2,sqrVal=TRUE,col.minval="gray",col.conf="gray", col.prof="black",conf=c(0.95,0.99), xlab="Attack rate (a)", ylab=expression(paste(Delta,"Negative log-likelihood")),axes=FALSE, ylim=c(0,8)) par(bty="u") box() axis(side=1) axis(side=2,at=c(0,4,8),labels=c(0,2,4)) axis(side=4,at=qchisq(c(0.95,0.99),df=1), labels=c(expression(frac(chi[1]^2*(0.95),2)), expression(frac(chi[1]^{2}*(0.99),2))),cex.axis=0.8) hts = qchisq(c(0.95,0.99),1) arrows(conflim95[1,1],hts[1],conflim95[1,2],hts[1],code=3) text(coef(m2)["a"],hts[1]+0.3,"95%") arrows(conflim99[1,1],hts[2],conflim99[1,2],hts[2],code=3) text(coef(m2)["a"],hts[2]+0.3,"99%") sqrprofplot(p2h,sqrVal=TRUE,col.minval="gray",col.conf="gray", col.prof="black",conf=c(0.95,0.99), xlab="Handling time (h)", ylab=expression(paste(Delta,"Negative log-likelihood")),axes=FALSE, ylim=c(0,8)) box() axis(side=1) axis(side=2,at=c(0,4,8),labels=c(0,2,4)) axis(side=4,at=qchisq(c(0.95,0.99),df=1), labels=c(expression(frac(chi[1]^2*(0.95),2)), expression(frac(chi[1]^{2}*(0.99),2))),cex.axis=0.8) arrows(conflim95[2,1],hts[1],conflim95[2,2],hts[1],code=3) text(coef(m2)["h"],hts[1]+0.3,"95%") arrows(conflim99[2,1],hts[2],conflim99[2,2],hts[2],code=3) text(coef(m2)["h"],hts[2]+0.3,"99%") par(op) ################################################### ### chunk number 55: ################################################### alpha=c(0.95,0.95,0.99,0.99) cval = qchisq(alpha,1)/2 cval2 = -logLik(m2)+cval var = rep(c("$a$","$h$"),2) vals <- rbind(conflim95,conflim99) rownames(vals) <- NULL ctab = data.frame(alpha,cval,cval2,var,vals) ctab[,-4] = signif(ctab[,-4],3) ctab[2,1:3]=ctab[4,1:3]=NA colnames(ctab) <- c("$\\alpha$","$\\frac{\\chi_1^2(\\alpha)}{2}$", "$-\\llik+\\frac{\\chi_1^2(\\alpha)}{2}$", "variable","lower","upper") latex(ctab,file="",table.env=FALSE,title="",rowname=NULL) ################################################### ### chunk number 56: ################################################### attach(ReedfrogFuncresp) ################################################### ### chunk number 57: ################################################### p2 = profile(m2) confint(p2) ################################################### ### chunk number 58: ################################################### detach(ReedfrogFuncresp) ################################################### ### chunk number 59: ################################################### binomNLL2.a = function(p,N,k,a) { h = p[1] p = a/(1+a*h*N) -sum(dbinom(k,prob=p,size=N,log=TRUE)) } ################################################### ### chunk number 60: ################################################### avec = seq(0.3,0.8,length=100) aprof = numeric(100) for (i in 1:100) { aprof[i] = optim(binomNLL2.a, par=0.02,k=ReedfrogFuncresp$Killed,N=ReedfrogFuncresp$Initial, a = avec[i], method="BFGS")$value } ################################################### ### chunk number 61: ################################################### prof.lower = aprof[1:which.min(aprof)] prof.avec = avec[1:which.min(aprof)] ################################################### ### chunk number 62: ################################################### approx(prof.lower,prof.avec,xout=-logLik(m2)+qchisq(0.95,1)/2) ################################################### ### chunk number 63: ################################################### x = subset(ReedfrogPred,pred=="pred" & density==10 & size=="small") k = x$surv op=par(lwd=2,bty="l",las=1,cex=1.5) curve(dbeta(x,shape1=sum(k)+1,shape2=40-sum(k)+1), xlab="Predation probability\nper capita", ylab="Probability density", from=0.4,to=1,yaxs="i") c1 = tcredint("beta",list(shape1=sum(k)+1,shape2=40-sum(k)+1), verbose=TRUE) v = with(as.list(c1),seq(lower,upper,length=100)) w = dbeta(v,shape1=sum(k)+1,shape2=40-sum(k)+1) polygon(c(v,rev(v)),c(w,rep(0,length(w))),col="gray") curve(dbeta(x,shape1=sum(k)+1,shape2=40-sum(k)+1),add=TRUE) abline(h=c1["p"],lty=2) qs = qbeta(c(0.025,0.975),shape1=sum(k)+1,shape2=40-sum(k)+1) v2 = seq(0.4,qs[1],length=100) w2 = dbeta(v2,shape1=sum(k)+1,shape2=40-sum(k)+1) polygon(c(v2,rev(v2)),c(w2,rep(0,length(w2))),density=10) v3 = seq(qs[2],1,length=100) w3 = dbeta(v3,shape1=sum(k)+1,shape2=40-sum(k)+1) polygon(c(v3,rev(v3)),c(w3,rep(0,length(w3))),density=10) text(0.75,2.1,"95%\ncredible\ninterval") text(0.5,1.4,"2.5% tails") arrows(c(0.5,0.5),c(1.2,1.2),c(0.58,0.88), c(0.27,0.30),angle=15) par(op) rm(x) rm(k) ################################################### ### chunk number 64: ################################################### post1 = with(frogpred1.bugs$sims.list,kde2d(a,h,n=40)) dx = diff(post1$x[1:2]) dy = diff(post1$y[1:2]) sz = sort(post1$z) c1 = cumsum(sz)*dx*dy c95 = approx(c1,sz,xout=0.05) dens.h = density(frogpred1.bugs$sims.list$h, from=post1$y[1],to=post1$y[length(post1$y)],n=length(post1$y)) dens.a = density(frogpred1.bugs$sims.list$a, from=post1$x[1],to=post1$x[length(post1$x)],n=length(post1$x)) frog.coda <- lump.mcmc.list(as.mcmc.bugs(frogpred1.bugs)) cred.h = HPDinterval(frog.coda[,"h"]) cred.a = HPDinterval(frog.coda[,"a"]) ## plot(sz,c1,type="l") ## abline(v=qval$y) ## abline(h=qval$x) ## sum(sz[sz>qval$y])*dx*dy ## get contour lines, discard edgy ones cl1 = contourLines(post1$x,post1$y,post1$z,level=c95)[[5]] mean.a <- mean(frogpred1.bugs$sims.list$a) mean.h <- mean(frogpred1.bugs$sims.list$h) wmode <- which(post1$z==max(post1$z),TRUE) mode.a <- post1$x[wmode[1]] mode.h <- post1$y[wmode[2]] ################################################### ### chunk number 65: ################################################### op=par(lwd=2,bty="l",las=1,cex=1.5) nf <- layout(matrix(c(2,1,0,3), 2, 2, byrow=TRUE), widths=c(1,3),heights=c(3,1)) par(cex=1.5) plot(cl1$x,cl1$y,type="l", xlab="Attack rate",ylab="",lwd=1, xlim=range(post1$x),ylim=range(post1$y)) mtext("Handling time",side=2,at=0.018,line=3,cex=1.5,las=0) points(mean.a,mean.h,lwd=1) points(mode.a,mode.h,pch=2,lwd=1) points(coef(m2)[1],coef(m2)[2],pch=3,lwd=1) text(c(0.55,0.53,0.49),c(0.0176,0.0149,0.0174), c("mean","mode","MLE"),adj=0,cex=0.7) contour(bsurf$x,bsurf$y,bsurf$z,levels=-logLik(m2)+qchisq(0.95,2)/2, add=TRUE, drawlabels=FALSE,lty=2) legend("topright", c("bivariate credible region", "bivariate confidence region"), lty=1:2, lwd=2,bty="n") oldmar = par(mar=c(5.1,0.5,4.1,0.5))$mar plot(-dens.h$y,dens.h$x,type="l",axes=FALSE,xlab="",ylab="",xaxs="i", xlim=c(-81,0)) axis(side=4,labels=FALSE) axis(side=1,at=c(-80,0),labels=c(80,0)) dh <- subset(data.frame(x=dens.h$x,y=dens.h$y), dens.h$x>cred.h[,"lower"] & dens.h$xcred.a[,"lower"] & dens.a$x>= op=par(lwd=2,bty="l",las=1,cex=1.5,mar=c(0,0,0,0)) spc = 0.03 ## y-spacing: was 0.025 ht=0.014 ## rectangle height? was 0.012 wid=0.08 cex1 = 0.7 plot(0:1,0:1,axes=FALSE,type="n",ann=FALSE) for (i in 1:nlevel) { ## draw segments p = 0.05 if (i