################################################### ### chunk number 1: ################################################### library(bbmle) library(emdbook) library(plotrix) library(lattice) library(Hmisc) library(MASS) source("chapskel.R") ################################################### ### chunk number 2: ################################################### sizefun = function(S,phi=20,eps=1,beta=1) { exp(eps*(phi-S))/(1+exp(beta*eps*(phi-S))) } maxxval = function(phi=20,eps=1,beta=1) { phi+log(beta-1)/(beta*eps) } maxyval = function(beta) { (beta-1)^(-1/beta)/(1+1/(beta-1)) } ################################################### ### chunk number 3: ################################################### op = par(cex=1.2,mar=c(4,4,0,1)+0.1, mgp=c(2.5,0.75,0),las=1,bty="l") curve(sizefun(x),from=5,to=40,xlab="Prey size",ylab="Predation risk") curve(sizefun(x,beta=1.1),add=TRUE,lty=2) curve(sizefun(x,beta=3),add=TRUE,lty=3) curve(sizefun(x,beta=1.1,eps=5),add=TRUE,lty=4) text(x=c(19.5,5,6.5,28),y=c(0.9,0.5,0.08,0.5), c(expression(list(beta==1,epsilon==1)),expression(list(beta==1.1,epsilon==1)), expression(list(beta==1.1,epsilon==5)), expression(list(beta==3,epsilon==1))), adj=0) arrows(27.5,0.49,21.5,0.4) par(op) ################################################### ### chunk number 4: ################################################### data(ReedfrogSizepred) attach(ReedfrogSizepred) ################################################### ### chunk number 5: ################################################### modlogist = function(x,eps,beta,phi) { exp(eps*(phi-x))/(1+exp(beta*eps*(phi-x))) } ##ricker = function(x,a,b) { ## a*x*exp(-b*x) ##} powricker = function(x,a,b,alpha) { b*(x/a*exp(1-x/a))^alpha } tricker = function(x,a,b,t,min=0.0001) { ifelse(x=1] = 1-minprop prop } ################################################### ### chunk number 26: ################################################### NLL.rogerscomb = function(a,h,eps,beta,phi,T=combT,P=combP) { if (h<0) return(NA) prob = prop.eaten(combInit,combSize,h,P,T,eps,beta,phi) dprob = dbinom(combKilled,prob=prob,size=combInit,log=TRUE) -sum(dprob) } ################################################### ### chunk number 27: ################################################### startvals = c(list(h=coef(FFR.rogers)["h"]),as.list(coef(FSP.modlogist))) ################################################### ### chunk number 28: ################################################### FPcomb = mle2(NLL.rogerscomb,start=startvals, method="Nelder-Mead") confint(FPcomb,method="quad") ################################################### ### chunk number 29: eval=FALSE ################################################### ## FPci2 = try(confint(FPcomb),silent=TRUE) ## startvals2 = list(h=0.88254539,eps=-0.01068388, ## beta=0.93631205,phi=383.67804375) ## FPcomb2 = mle2(NLL.rogerscomb,start=startvals2, ## method="Nelder-Mead",control=list(parscale=c(1,1,1,500))) ################################################### ### chunk number 30: eval=FALSE ################################################### ## cat(h,eps,beta,phi,"\n") ################################################### ### chunk number 31: eval=FALSE ################################################### ## if (any(!is.finite(prob))) cat("NAs:",h,eps,beta,phi,"\n") ################################################### ### chunk number 32: eval=FALSE ################################################### ## if (any(!is.finite(dprob))) { ## browser() ## } ################################################### ### chunk number 33: ################################################### FPcomb = mle2(NLL.rogerscomb,start=startvals, method="L-BFGS-B", lower=c(0.7,0.5,1,14), upper=c(1.8,2.25,2,20), control=list(parscale=c(1,1,1,10))) FPcomb.ci = confint(FPcomb) ################################################### ### chunk number 34: ################################################### NLL.rogerscomb.bb = function(a,h,eps,beta,phi, theta,T=combT,P=combP) { if (h<0) return(NA) prob = prop.eaten(combInit,combSize,h,P,T,eps,beta,phi) dprob = dbetabinom(combKilled,prob=prob,size=combInit,theta=theta,log=TRUE) -sum(dprob) } FPcomb.bb = mle2(NLL.rogerscomb.bb,start=c(startvals,list(theta=100)), method="L-BFGS-B", lower=c(0.7,0.5,1,14,0.1), upper=c(1.8,2.25,2,20,Inf), control=list(parscale=c(1,1,1,10,100))) theta.ci = confint(FPcomb.bb,parm="theta") ################################################### ### chunk number 35: eval=FALSE ################################################### ## c1 = coef(FSP.modlogist) ## FSP.expprop.mean = modlogist(12.8,c1["eps"],c1["beta"],c1["phi"]) ## FSP.exppropvar = deltavar(exp(eps*(phi-12.8))/(1+exp(beta*eps*(phi-12.8))), ## meanval=coef(FSP.modlogist),Sigma=vcov(FSP.modlogist)) ## FSP.expprop.deltaci = FSP.expprop.mean+c(-1.96,1.96)*sqrt(FSP.exppropvar) ################################################### ### chunk number 36: ################################################### c1 = coef(FSP.modlogist) FSP.expprop.mean = modlogist(12.8,c1["eps"],c1["beta"],c1["phi"]) ################################################### ### chunk number 37: ################################################### c2 = coef(FPcomb) FP.expprop.mean = prop.eaten(N0=10,S=12.8,c2["h"], P=2/0.08,T=3,eps=c2["eps"],beta=c2["beta"],phi=c2["phi"]) ################################################### ### chunk number 38: ################################################### set.seed(1001) FSP.expprop.pars = mvrnorm(5000,mu=c1,Sigma=vcov(FSP.modlogist)) FSP.expprop.val = numeric(5000) for (i in 1:5000) { FSP.expprop.val[i] = modlogist(12.8,FSP.expprop.pars[i,1], FSP.expprop.pars[i,2],FSP.expprop.pars[i,3]) } FSP.expprop.ppi = quantile(FSP.expprop.val,c(0.025,0.975)) ################################################### ### chunk number 39: eval=FALSE ################################################### ## FP.exppropvar = deltavar(prop.eaten(N0=10,S=12.8,h,P=2/0.08,T=3,eps,beta,phi), ## meanval=coef(FPcomb),Sigma=vcov(FPcomb)) ## FP.expprop.deltaci = FP.expprop.mean+c(-1.96,1.96)*sqrt(FP.exppropvar) ################################################### ### chunk number 40: ################################################### FP.expprop.pars = mvrnorm(5000,mu=c2,Sigma=vcov(FPcomb)) FP.expprop.val = numeric(5000) for (i in 1:5000) { FP.expprop.val[i] = prop.eaten(N0=10,S=12.8,P=2/0.08,T=3, h=FP.expprop.pars[i,"h"], eps=FP.expprop.pars[i,"eps"], beta=FP.expprop.pars[i,"beta"], phi=FP.expprop.pars[i,"phi"]) } FP.expprop.ppi = quantile(FP.expprop.val,c(0.025,0.975)) ################################################### ### chunk number 41: ################################################### ##esttab = rbind(c(FSP.expprop.mean,NA,NA), ## c(NA,FSP.expprop.deltaci), ## c(NA,FSP.expprop.ppi), ## c(FP.expprop.mean,NA,NA), ## c(NA,FP.expprop.deltaci), ## c(NA,FP.expprop.ppi)) ## dimnames(esttab) = list(c("size-pred", ## "delta method", ## "PPI", ## "combined", ## "delta method", ## "PPI"), ## c("mean","low","high")) esttab = rbind(c(FSP.expprop.mean,FSP.expprop.ppi), c(FP.expprop.mean,FP.expprop.ppi)) dimnames(esttab) = list(c("size-pred","combined"), c("mean","low","high")) latex(round(esttab,3),file="",title="",table.env=FALSE) ################################################### ### chunk number 42: ################################################### spredarea = with(xpars.Sizepred,P*T/area) frpredarea = with(xpars.Funcresp,P*T/area) spredvol = with(xpars.Sizepred,P*T/vol) frpredvol = with(xpars.Funcresp,P*T/vol) ################################################### ### chunk number 43: ################################################### svec = seq(0,40,by=0.25) FP.expprop = prop.eaten(N0=10,S=svec,c2["h"], P=2/0.08,T=3,eps=c2["eps"],beta=c2["beta"],phi=c2["phi"]) FSP.expprop = with(as.list(coef(FSP.modlogist)), modlogist(svec,eps,beta,phi)) FP.expprop.pars = mvrnorm(5000,mu=c2,Sigma=vcov(FPcomb)) FP.expprop.vals = t(apply(FP.expprop.pars,1, function(x) { with(as.list(x),prop.eaten(N0=10,S=svec,P=2/0.08,T=3, h=h, eps=eps, beta=beta, phi=phi))})) FP.expprop.pars = mvrnorm(5000,mu=c2,Sigma=vcov(FPcomb)) FSP.expprop.vals = t(apply(FSP.expprop.pars,1, function(x) { with(as.list(x),modlogist(svec,eps,beta,phi))})) FSP.env = t(apply(FSP.expprop.vals,2,quantile,c(0.025,0.975))) FP.env = t(apply(FP.expprop.vals,2,quantile,c(0.025,0.975))) ################################################### ### chunk number 44: ################################################### op = par(lwd=2,bty="l",las=1,mgp=c(3,1,0),cex=1.5) svec = seq(0,40,by=0.25) with(ReedfrogSizepred,sizeplot(TBL,Kill/10,ylim=c(0,1), xlab="Total body length", ylab="Proportion eaten")) lines(svec,FP.expprop) lines(svec,FSP.expprop,lty=2) abline(v=12.8,col="gray") ## matlines(svec,FSP.env,lty=2,col="gray") ## matlines(svec,FP.env,lty=1,col="gray") ## c3 = coef(FPcomb2) ## FP.expprop2 = prop.eaten(N0=10,S=svec,c3["h"], ## P=2/0.08,T=3,eps=c3["eps"],beta=c3["beta"],phi=c3["phi"]) ## lines(svec,FP.expprop2,lty=3) ## bit = 0.3 ## plotCI(12.8-bit,FP.expprop.mean,li=FP.expprop.ppi[1],ui=FP.expprop.ppi[2], ## add=TRUE,pch=NA) ## plotCI(12.8+bit,FSP.expprop.mean,li=FSP.expprop.ppi[1],ui=FSP.expprop.ppi[2], ## add=TRUE,pch=NA,slty=2) legend("topright",c("size-pred.","combined"),lty=2:1,bty="n") par(op) ################################################### ### chunk number 45: ################################################### detach(ReedfrogSizepred) detach(ReedfrogFuncresp) ################################################### ### chunk number 46: ################################################### library(emdbookx) data(GobySurvival) attach(GobySurvival) ################################################### ### chunk number 47: ################################################### day1 = d1-1 day2 = ifelse(d2==70,Inf,d2-1) ################################################### ### chunk number 48: eval=FALSE ################################################### ## day1 = pmax(day1,.Machine$double.eps) ################################################### ### chunk number 49: ################################################### xmax = 1 op = par(cex=1.5,lwd=2, mgp=c(2.5,0.75,0),las=1,bty="l") par(mfrow=c(1,2)) shapes = c(1,0.5,0.2) curve(pweibull(x,shape=shapes[1],scale=0.5,lower.tail=FALSE),from=0.001,to=xmax,log="y", ylab="Fraction surviving",xlab="Time",ylim=c(0.1,1),axes=FALSE) axis(side=1,at=c(0,0.5,1)) axis(side=2) box() curve(pweibull(x,shape=shapes[2],scale=0.5,lower.tail=FALSE),add=TRUE,lty=2) curve(pweibull(x,shape=shapes[3],scale=0.5,lower.tail=FALSE),add=TRUE,lty=3) ## abline(h=exp(-1),col=2) ## abline(v=0.5,col=2) text(rep(xmax+0.1,3),pweibull(xmax,shape=shapes,scale=0.5,lower.tail=FALSE), paste("shape=",shapes,sep=""), adj=0,xpd=NA) text(corner.loc(x=1,y=1),adj=1,"scale=0.5",cex=1.5) par(mar=c(5,3,4,4)+0.1) scales = c(0.5,1,2) curve(pweibull(x,shape=0.5,scale=scales[1],lower.tail=FALSE),from=0.001,to=xmax, log="y", axes=FALSE,ylim=c(0.1,1),xlab="Time",ylab="") box() axis(side=1,at=c(0,0.5,1)) curve(pweibull(x,shape=0.5,scale=scales[2],lower.tail=FALSE),add=TRUE,lty=2) curve(pweibull(x,shape=0.5,scale=scales[3],lower.tail=FALSE),add=TRUE,lty=3) text(rep(xmax+0.1,3),pweibull(xmax,shape=0.5,scale=scales,lower.tail=FALSE), paste("scale=",scales,sep=""), adj=0,xpd=NA) ## abline(h=exp(-1)) text(corner.loc(x=1,y=1),adj=1,"shape=0.5",cex=1.5) par(op) ################################################### ### chunk number 50: ################################################### op = par(lwd=2,bty="l",las=1,mgp=c(3,1,0),cex=1.5) dens.cat = factor(ifelse(density>median(density),"high","low"), levels=c("low","high")) qual.cat = factor(ifelse(qual>median(qual),"high","low"), levels=c("low","high")) intcat = interaction(qual.cat,dens.cat) cattab = table(intcat) ## when (average) did they die? meansurv = (d1+d2)/2 ## tabulate morttab = table(meansurv,intcat) ## reverse order morttab = morttab[nrow(morttab):1,] ## get cumulative sum: this is number surviving until day x csurvtab = apply(morttab,2,cumsum) ## divide by total number in category cnsurvtab = sweep(csurvtab,2,cattab,"/") days = as.numeric(rownames(csurvtab)) matplot(days,cnsurvtab,type="s",xlab="Time (days)", ylab="Proportion of cohort surviving", log="y",col=1,xlim=c(0,40)) legend("topright", c("high quality/high density", "low quality/low density", "high quality/low density", "low quality/high density"),bty="n", lty=c(4,1,2,3),cex=0.9) par(op) ################################################### ### chunk number 51: eval=FALSE ################################################### ## numsurv = diag(table(d1,d2)) ## days = c(1,4,8,11,70) ## csurv = nrow(X)-c(0,cumsum(numsurv)) ## 1-exp(log(csurv[2:5]/csurv[1:4])/diff(days)) ################################################### ### chunk number 52: ################################################### NLL.GS.xqdi = function(lscale0,lscale.q,lscale.d,lscale.i, lscale.x2,lscale.x3,lscale.x4,lscale.x5, lshape) { lscalediff = c(0,lscale.x2,lscale.x3,lscale.x4,lscale.x5) scale=exp(lscale0+lscalediff[exper]+ lscale.q*qual+(lscale.d+lscale.i*qual)*density) shape = exp(lshape) -sum(log(pweibull(day2,shape,scale)- pweibull(day1,shape,scale))) } ################################################### ### chunk number 53: ################################################### totmeansurv = mean((d1+d2)/2) startvals.GS = list(lscale0=log(totmeansurv), lscale.x2=0,lscale.x3=0,lscale.x4=0,lscale.x5=0, lscale.q=0,lscale.d=0,lscale.i=0, lshape=0) GS.xqdi = mle2(NLL.GS.xqdi,startvals.GS) ################################################### ### chunk number 54: eval=FALSE ################################################### ## ## this also works! ## dicweib = function(x,shape,scale,log=FALSE) { ## if (is.matrix(x)) { ## day1 = x[,1] ## day2 = x[,2] ## } else { ## day1 = x[1] ## day2 = x[2] ## } ## v = log(pweibull(day2,shape,scale)-pweibull(day1,shape,scale)) ## if (log) v else exp(v) ## } ## ## fexper = factor(exper) ## mle2(cbind(day1,day2)~dicweib(exp(shape),exp(scale)), ## parameters=list(scale~fexper+qual*density), ## start=list(scale=log(totmeansurv),shape=0)) ################################################### ### chunk number 55: ################################################### summary(GS.xqdi) ################################################### ### chunk number 56: ################################################### GS.xqd = mle2(NLL.GS.xqdi,startvals.GS,fixed=list(lscale.i=0)) ################################################### ### chunk number 57: ################################################### GS.xd = mle2(NLL.GS.xqdi,startvals.GS, fixed=list(lscale.i=0,lscale.q=0)) GS.x = mle2(NLL.GS.xqdi,startvals.GS, fixed=list(lscale.i=0,lscale.q=0,lscale.d=0)) ################################################### ### chunk number 58: ################################################### anova(GS.xqdi,GS.xqd,GS.xd,GS.x) ################################################### ### chunk number 59: ################################################### GS.qdi = mle2(NLL.GS.xqdi,startvals.GS,fixed=list(lscale.x2=0,lscale.x3=0,lscale.x4=0,lscale.x5=0)) GS.qd = mle2(NLL.GS.xqdi,startvals.GS, fixed=list(lscale.i=0, lscale.x2=0,lscale.x3=0,lscale.x4=0,lscale.x5=0)) GS.xq = mle2(NLL.GS.xqdi,startvals.GS, fixed=list(lscale.i=0,lscale.d=0)) GS.d = mle2(NLL.GS.xqdi,startvals.GS, fixed=list(lscale.i=0,lscale.x2=0,lscale.x3=0,lscale.x4=0, lscale.x5=0,lscale.q=0)) GS.q = mle2(NLL.GS.xqdi,startvals.GS, fixed=list(lscale.i=0, lscale.x2=0,lscale.x3=0,lscale.x4=0,lscale.x5=0,lscale.d=0)) GS.0 = mle2(NLL.GS.xqdi,startvals.GS, fixed=list(lscale.i=0,lscale.x2=0,lscale.x3=0,lscale.x4=0, lscale.x5=0,lscale.d=0,lscale.q=0)) ################################################### ### chunk number 60: ################################################### a1 = AICtab(GS.xqdi,GS.xqd,GS.qdi,GS.xq,GS.xd,GS.qd,GS.x,GS.q,GS.d,GS.0,delta=TRUE,weights=TRUE) a2 = AICctab(GS.xqdi,GS.xqd,GS.qdi,GS.xq,GS.xd,GS.qd,GS.x,GS.q,GS.d,GS.0, nobs=nrow(GobySurvival),delta=TRUE) b1 = BICtab(GS.xqdi,GS.xqd,GS.qdi,GS.xq,GS.xd,GS.qd,GS.x,GS.q,GS.d,GS.0, nobs=nrow(GobySurvival),delta=TRUE) atab = cbind(a1$df,a1$dAIC,a1$weight,a2$dAICc,b1$dBIC) rownames(atab) = gsub("GS\\.","",attr(a1,"row.names")) ## need to double backslashes colnames(atab) = c("params","$\\Delta \\mbox{AIC}$","AIC weights","$\\Delta \\mbox{AIC}_c$","$\\Delta \\mbox{BIC}$") atab = round(atab[order(atab[,2]),],3) ################################################### ### chunk number 61: ################################################### latex(round(atab,2),file="",table.env=FALSE,title="model") ################################################### ### chunk number 62: ################################################### data(SeedPred) ################################################### ### chunk number 63: ################################################### SeedPred = na.omit(subset(SeedPred,available>0)) attach(SeedPred) ################################################### ### chunk number 64: ################################################### nz = subset(SeedPred,taken>0) ################################################### ### chunk number 65: eval=FALSE ################################################### ## barchart(table(nz$taken,nz$available,nz$dist,nz$species),stack=FALSE) ## barchart(table(nz$taken,nz$species,nz$dist,nz$available),stack=FALSE) ## barchart(table(nz$species,nz$available,nz$dist,nz$taken),stack=FALSE) ## barchart(table(nz$available,nz$dist,nz$taken),stack=FALSE) ## barchart(table(nz$available,nz$species,nz$taken),stack=FALSE) ################################################### ### chunk number 66: ################################################### barchart(table(available,dist,taken),stack=FALSE) ################################################### ### chunk number 67: eval=FALSE ################################################### ## tcumfac = cut(nz$tcum,breaks=c(0,20,40,60,180)) ## barchart(table(nz$available,tcumfac,nz$taken),stack=FALSE) ## barchart(table(available,tcumfac,taken),stack=FALSE) ################################################### ### chunk number 68: ################################################### pcomb = table(nz$available,nz$taken) pcomb = sweep(pcomb,1,rowSums(pcomb),"/") trellis.par.set(canonical.theme(color=FALSE)) print(barchart(pcomb[-1,],stack=FALSE,auto.key=list(x=0.8,y=0.8,corner=c(1,1), title="Taken",reverse.rows=TRUE), xlab="Frequency",ylab="Seeds available")) ################################################### ### chunk number 69: ################################################### dzibinom = function(x,prob,size,zprob,log=FALSE) { logv = log(1 - zprob) + dbinom(x, prob = prob, size = size, log = TRUE) logv = ifelse(x == 0, log(zprob + exp(logv)), logv) if (log) logv else exp(logv) } dzibb = function(x,size,prob,theta,zprob,log=FALSE) { logv = ifelse(x>size, NA,log(1 - zprob) + dbetabinom(x, prob = prob, size = size, theta=theta, log = TRUE)) logv = ifelse(x == 0, log(zprob + exp(logv)), logv) if (log) logv else exp(logv) } ################################################### ### chunk number 70: ################################################### SP.zibb = mle2(taken~dzibb(size=available,prob,theta,plogis(logitzprob)), start=list(prob=0.5,theta=1,logitzprob=0)) print(SP.zibb) ################################################### ### chunk number 71: ################################################### cov2cor(vcov(SP.zibb)) ################################################### ### chunk number 72: ################################################### SP.bb = mle2(taken~dbetabinom(prob,theta,size=available), start=list(prob=0.5,theta=1)) ################################################### ### chunk number 73: ################################################### logLik(SP.bb)-logLik(SP.zibb) ################################################### ### chunk number 74: ################################################### SP.zib = mle2(taken~dzibinom(size=available,prob=p, zprob=plogis(logitzprob)), start=list(p=0.2,logitzprob=0)) ################################################### ### chunk number 75: ################################################### AICtab(SP.zib,SP.zibb,SP.bb,sort=TRUE,weights=TRUE) ################################################### ### chunk number 76: ################################################### comb = table(taken,available) pcomb = sweep(comb,2,colSums(comb),"/") ################################################### ### chunk number 77: ################################################### mtab = matrix(0,nrow=6,ncol=5) for (N in 1:5) { cvals = coef(SP.zibb) mtab[1:(N+1),N] = dzibb(0:N,size=N,prob=cvals["prob"], theta=cvals["theta"], zprob=plogis(cvals["logitzprob"])) } ################################################### ### chunk number 78: ################################################### mtab2 = matrix(0,nrow=6,ncol=5) for (av in 1:5) { mtab2[1:(av+1),av] = dbetabinom(0:av,prob=coef(SP.bb)["prob"], theta=coef(SP.bb)["theta"],size=av) } mtab3 = matrix(0,nrow=6,ncol=5) for (av in 1:5) { mtab3[1:(av+1),av] = dzibinom(0:av,prob=coef(SP.zib)["p"], zprob=plogis(coef(SP.zib)["logitzprob"]), size=av) } ################################################### ### chunk number 79: ################################################### tmpf = function(m) { m2 = m[-1,] sweep(m2,2,colSums(m2),"/") } op = par(mfrow=c(2,2),mar=c(2,2,2,2)) ##mosaicplot(tmpf(pcomb),main="data") barplot(tmpf(pcomb),main="data",axes=FALSE) barplot(tmpf(mtab),main="Z-I beta-binomial",axes=FALSE) barplot(tmpf(mtab2),main="beta-binomial",axes=FALSE) barplot(tmpf(mtab3),main="Z-I binomial",axes=FALSE) par(op) ################################################### ### chunk number 80: ################################################### pval = numeric(5) for (N in 1:5) { obs = comb[1:(N+1),N] prob = mtab[1:(N+1),N] pval[N] = chisq.test(obs,p=prob)$p.value } ################################################### ### chunk number 81: ################################################### latex(t(matrix(format.pval(pval,digits=2,eps=0.001), dimnames=list(1:5,NULL))),file="",title="",table.env=FALSE) ################################################### ### chunk number 82: ################################################### SP.bb.dist = mle2(taken~dbetabinom(prob,size=available,theta), parameters=list(prob~dist-1,theta~dist-1),start=as.list(coef(SP.bb))) ################################################### ### chunk number 83: ################################################### anova(SP.bb,SP.bb.dist) ################################################### ### chunk number 84: ################################################### startvals =list(lprob=qlogis(coef(SP.bb.dist)["prob.dist10"]), ltheta=log(coef(SP.bb.dist)["theta.dist10"])) ################################################### ### chunk number 85: ################################################### SP.bb.dist2 = mle2(taken~dbetabinom(plogis(lprob),exp(ltheta),size=available), parameters=list(lprob~dist,ltheta~dist),start=startvals) ################################################### ### chunk number 86: ################################################### summary(SP.bb.dist2) ################################################### ### chunk number 87: ################################################### SP.bb.probdist = mle2(taken~dbetabinom(plogis(lprob),exp(ltheta),size=available), parameters=list(lprob~dist),start=startvals) ################################################### ### chunk number 88: ################################################### anova(SP.bb,SP.bb.probdist,SP.bb.dist) AICtab(SP.bb,SP.bb.probdist,SP.bb.dist,sort=TRUE,weights=TRUE) ################################################### ### chunk number 89: ################################################### c1 = coef(SP.bb.probdist) plogis(c(c1[1],c1[1]+c1[2])) ################################################### ### chunk number 90: ################################################### SP.bb.sp = mle2(taken~dbetabinom(plogis(lprob),exp(ltheta),size=available), parameters=list(lprob~species,ltheta~species),start=startvals) ################################################### ### chunk number 91: ################################################### s1 = summary(SP.bb.sp) s1@coef = s1@coef[,-3] ## drop z value printCoefmat(s1@coef,has.Pvalue=TRUE) ################################################### ### chunk number 92: ################################################### SP.bb.probsp = mle2(taken~dbetabinom(plogis(lprob),exp(ltheta),size=available), parameters=list(lprob~species),start=startvals) ################################################### ### chunk number 93: ################################################### anova(SP.bb.sp,SP.bb.probsp,SP.bb) AICtab(SP.bb.sp,SP.bb.probsp,SP.bb,sort=TRUE,weights=TRUE) ################################################### ### chunk number 94: eval=FALSE ################################################### ## tmpfun = calc_mle2_function(taken~dbetabinom(plogis(lprob),exp(ltheta),size=available), ## parameters=list(lprob~species-1),start=startvals)$fn ## pars = coef(SP.bb.probsp) ## pars[2:8] = pars[2:8]+pars[1] ## names(pars) = NULL ## do.call("tmpfun",as.list(pars)) ################################################### ### chunk number 95: ################################################### SP.bb.probsp0 = mle2(taken~dbetabinom(plogis(lprob),exp(ltheta),size=available), parameters=list(lprob~species-1),start=startvals, method="L-BFGS-B",lower=rep(-10,9),upper=rep(10,9)) ################################################### ### chunk number 96: eval=FALSE ################################################### ## lprob0 = coef(SP.bb.probsp)["lprob0"] ## lprobdiff = lprob0+coef(SP.bb.probsp)[2:nsp] ## probvals = c(lprob0,lprob0+lprobdiff) ## predprob = plogis(probvals) ################################################### ### chunk number 97: eval=FALSE ################################################### ## sdvals = coef(summary(SP.bb.probsp))[,"Std. Error"] ## sdvals0 = sdvals["lprob0"] ## sdvalsdiff = sqrt(sdvals[1]^2+sdvals[2:nsp]^2) ## sdvals = c(sdvals0,sdvalsdiff) ################################################### ### chunk number 98: eval=FALSE ################################################### ## ci.lo = plogis(probvals-1.96*sdvals) ## ci.hi = plogis(probvals+1.96*sdvals) ################################################### ### chunk number 99: ################################################### predprob = plogis(coef(SP.bb.probsp0))[1:8] SP.bb.ci = plogis(confint(SP.bb.probsp0,method="quad"))[1:8,] ################################################### ### chunk number 100: ################################################### op = par(cex=1.5,mar=c(4,4,0,1)+0.1,bty="l",mgp=c(2.5,0.75,0),las=1) data(SeedPred_mass) plotCI(SeedPred_mass,predprob,li=SP.bb.ci[,1],ui=SP.bb.ci[,2],pch=NA,gap=0.012, xlab="Seed mass",ylab="Removal probability (p)") textrec(SeedPred_mass,predprob,levels(SeedPred$species),cex=0.75, bg="white",expand=1.2) text(SeedPred_mass,predprob,levels(SeedPred$species),cex=0.75) box() par(op) ################################################### ### chunk number 101: eval=FALSE ################################################### ## NLL.bb.probsptrans = function(p) { ## nsp = 8 ## lprob0 = p[1] ## lprobdiff = p[2:nsp] ## lprobtr = p[nsp+1] ## lprobtrdiff = p[(nsp+2):(2*nsp)] ## lprob = c(lprob0,lprob0+lprobdiff, ## lprob0+lprobtr,lprob0+lprobtr+lprobdiff+lprobtrdiff) ## probvec = plogis(lprob) ## ltheta0 = p[2*nsp+1] ## prob = probvec[interaction(SeedPred$species,SeedPred$dist)] ## theta = exp(ltheta0) ## -sum(dbetabinom(taken,size=available,prob=prob,theta=theta,log=TRUE)) ## } ## parnames(NLL.bb.probsptrans) = c(parnames(NLL.bb.sp)[1:nsp], ## paste("tr",parnames(NLL.bb.sp)[1:nsp],sep=""), ## "ltheta") ## svec.probtr = c(svec[1:nsp],rep(0,nsp),svec[nsp+1]) ## names(svec.probtr) = parnames(NLL.bb.probsptrans) ## SP.bb.probsptrans = mle2(NLL.bb.probsptrans,start=svec.probtr,vecpar=TRUE) ################################################### ### chunk number 102: ################################################### SP.bb.probspdist = mle2(taken~dbetabinom(plogis(lprob),exp(ltheta),size=available), parameters=list(lprob~species*dist),start=startvals, method="L-BFGS-B",lower=rep(-10,9),upper=rep(5,9)) ################################################### ### chunk number 103: ################################################### AICtab(SP.bb,SP.bb.probsp,SP.bb.probspdist,SP.bb.sp,SP.bb.probdist,SP.bb.dist, weights=TRUE,sort=TRUE) ################################################### ### chunk number 104: eval=FALSE ################################################### ## mean.prop.taken = tapply(taken/available,tint,mean,na.rm=TRUE) ## sd.prop.taken = tapply(taken/available,tint,sd,na.rm=TRUE) ## n.tint = table(tint) ## se.prop.taken = sd.prop.taken/sqrt(n.tint) ################################################### ### chunk number 105: ################################################### par(mfrow=c(1,2)) op = par(cex=1.5, mar=c(4,4,0,1)+0.1,bty="l",mgp=c(2.5,0.75,0), bty="l",las=1) nint = table(tint) mtaken = tapply(taken/available,tint,mean,na.rm=TRUE) setaken = tapply(taken/available,tint,sd,na.rm=TRUE)/sqrt(nint) par(mfrow=c(1,2)) plotCI(3:10,mtaken,setaken,xlab="Interval (days)",ylab="Proportion taken") corner.label2("a","topleft",inset=0.025) plot(date,tint,xlab="Date",ylab="Interval (days)") corner.label2("b","topleft",inset=0.025) par(op) ################################################### ### chunk number 106: eval=FALSE ################################################### ## NLL.bb.probspdate = function(p) { ## nsp = 8 ## lprob0 = p[1] ## lprobdiff = p[2:nsp] ## lprob = c(lprob0,lprob0+lprobdiff) ## probvec = plogis(lprob) ## ltheta0 = p[nsp+1] ## date = p[nsp+2] ## prob = probvec[SeedPred$species]*exp(-SeedPred$tcum*date) ## theta = exp(ltheta0) ## -sum(dbetabinom(taken,size=available,prob=prob,theta=theta,log=TRUE)) ## } ## parnames(NLL.bb.probspdate) = c(parnames(NLL.bb.sp)[1:nsp],"ltheta","date") ## svec.probspdate = c(svec[1:(nsp+1)],0.001) ## names(svec.probspdate) = parnames(NLL.bb.probspdate) ## SP.bb.probspdate = mle2(NLL.bb.probspdate,start=svec.probspdate,vecpar=TRUE) ################################################### ### chunk number 107: ################################################### SP.bb.probspdate = mle2(taken~dbetabinom(plogis(lprob)*exp(-tcum*date), exp(ltheta),size=available), parameters=list(lprob~species), start=c(startvals,date=0), method="L-BFGS-B",lower=c(rep(-10,9),0),upper=c(rep(5,9),2)) ################################################### ### chunk number 108: ################################################### SP.bb.probdate = mle2(taken~dbetabinom(plogis(lprob)*exp(-tcum*date), size=available,exp(ltheta)), start=c(startvals,date=0), method="L-BFGS-B",lower=c(rep(-10,2),0),upper=c(rep(5,2),2)) ################################################### ### chunk number 109: ################################################### op = par(cex=1.5,mar=c(4,4,0,1)+0.1,bty="l",mgp=c(2,0.75,0),las=1) dates = unique(date) ndate = table(date) mtaken = tapply(taken/available,date,mean,na.rm=TRUE) sdtaken = tapply(taken/available,date,sd,na.rm=TRUE)/sqrt(ndate) plotCI(dates,mtaken,sdtaken,xlab="Date",ylab="",ylim=c(0.001,0.4)) mtext(side=2,"Proportion taken",line=2.5,cex=1.5,las=0) t0 = as.numeric(SeedPred$date[1]) p0 = plogis(coef(SP.bb.probdate)["lprob"]) dpar = coef(SP.bb.probdate)["date"] curve(p0*exp(-dpar*(x-t0)),add=TRUE) par(op) ################################################### ### chunk number 110: ################################################### detach(SeedPred)