################################################### ### chunk number 1: ################################################### ## badfonts <- .Platform$OS.type=="unix" badfonts <- FALSE library(Hmisc,warn.conflicts=FALSE) source("chapskel.R") ################################################### ### chunk number 2: ################################################### ## read in and attach seed predation data library(emdbook) data(SeedPred) seedsize <- c(0.11,0.03,0.21,0.16,0.35,0.01,0.52,0.27) ss.levels <- levels(SeedPred$species)[order(seedsize)] SeedPred$species <- factor(SeedPred$species,levels=ss.levels) seedsub = subset(SeedPred,seeds>0 & !is.na(taken)) attach(seedsub,warn.conflicts=FALSE) ## grab useful? libraries invisible(require(lattice,quietly=TRUE)) invisible(require(plotrix,quietly=TRUE,warn.conflicts=FALSE)) ## compute none/any values for number of seeds taken by species t1 <- table(taken>0,species) ## total number of obs per species tot <- table(species) ## number of "successes" and proportion for all species succ <- t1[2,] props <- succ/tot ## restrict to just the first two species spp <- c("pol","psd") succ <- succ[spp] tot <- tot[spp] vals <- t1[2:1,spp] ## reverse order ## MLE estimates of binomial p p.mle <- succ/tot ptot.mle <- sum(succ)/sum(tot) obs.pratio <- p.mle["pol"]/p.mle["psd"] ## let's be frequentist: Fisher exact test F1 <- fisher.test(vals) F1.1 <- fisher.test(vals,alternative="greater") ## likelihood: restricted and full likelihood pr <- sum(succ)/sum(tot) L0 <- dbinom(sum(succ),size=sum(tot),prob=pr,log=TRUE) L0 <- sum(dbinom(succ,size=tot,prob=pr,log=TRUE)) ## ??? sum(log-likelihoods) neq log(pooled data)??? L0x <- dbinom(succ,size=tot,prob=pr,log=TRUE) L0bx <- lchoose(tot,succ)+succ*log(pr)+log(1-pr)*(tot-succ) L0b <- lchoose(sum(tot),sum(succ))+sum(succ)*log(pr)+log(1-pr)*(sum(tot)-sum(succ)) ## correct! L1 <- sum(dbinom(succ,size=tot,prob=succ/tot,log=TRUE)) g1 <- glm(t(vals)~factor(1:2),family="binomial") g2 <- glm(t(vals)~1,family="binomial") ## anova(g1,g2,test="Chisq") ## ??? logLik(g1) matches L1 ## logLik(g2) = -22.23, L0 = -2.84 ## plogis(coef(g2)) ## plogis(coef(g1)[1]) ## plogis(sum(coef(g1))) ## likelihood ratio test: LRT.p <- pchisq(2*(L1-L0),df=1,lower.tail=FALSE) OR.to.PR <- function(OR,p1) { O1 <- p1/(1-p1) O2 <- O1*OR p2 <- O2/(1+O2) p2/p1 } pr.conf <- OR.to.PR(F1$conf,p.mle[2]) ## calculate num taken from odds ratio, ## given margins v1,v2,v3 OR.to.nums <- function(OR,vals) { v1 <- colSums(vals)[1] ## total pol v2 <- colSums(vals)[2] ## total psd v3 <- rowSums(vals)[1] ## total num taken } ## calculate odds ratio given number in [1,1] x.to.OR <- function(x,vals) { v1 <- colSums(vals)[1] ## total pol v2 <- colSums(vals)[2] ## total psd v3 <- rowSums(vals)[1] ## total num taken m <- matrix(nrow=2,ncol=2) m[1,1] <- x m[2,1] <- v1-x m[1,2] <- v3-x m[2,2] <- v2-m[1,2] (m[1,1]/m[2,1])/(m[1,2]/m[2,2]) } ## PR = P1/P2 = (m[1,1]/v1)/(m[1,2]/v2) ## m[1,2] = v3-m[1,1] ## PR = (m[1,1]/v1)/((v3-m[1,1])/v2) ## ((v3-m[1,1])/v2)*PR = m[1,1]/v1 ## v3/v2*PR = m[1,1]*(1/v1+PR/v2) ## m[1,1] = v3/v2*PR/(1/v1+PR/v2) ## calculate number in [1,1] given prob. ratio, ## marginals PR.to.x <- function(PR,vals) { v1 <- colSums(vals)[1] ## total pol v2 <- colSums(vals)[2] ## total psd v3 <- rowSums(vals)[1] ## total num taken v3/v2*PR/(1/v1+PR/v2) } ors <- sapply(0:40,x.to.OR,vals=vals) prs <- sapply(ors,OR.to.PR,p1=vals[1,1]/colSums(vals)[1]) or2 <- pmax(ors,1/ors) smvals <- (0:40)[ors<1 & or2>F1$estimate] lgvals <- vals[1,1]:40 ################################################### ### chunk number 3: ################################################### vals2 = rbind(vals,colSums(vals)) dimnames(vals2)[[1]] <- c("any taken ($t$)","none taken","total ($N$)") latex(vals2,file="",title="",table.env=FALSE) ################################################### ### chunk number 4: ################################################### op <- par(cex=1.5,lwd=2,bty="l",yaxs="i") par(mar=c(5,4,4,2)+0.1) ## stuff for Fisher's test: one-sided test ## vals maxv <- 40 xvec <- 0:maxv v1 <- colSums(vals)[1] ## total pol v2 <- colSums(vals)[2] ## total psd v3 <- rowSums(vals)[1] ## total num taken hvec <- dhyper(xvec,v1,v2,v3,log=TRUE)/log(10) ## convert to log-10 units hvec2 <- dhyper(xvec,v1,v2,v3) ##phyper(29,125,439,52,lower.tail=FALSE) pval <- fisher.test(vals,alternative="greater")$p.val pval2 <- fisher.test(vals,alternative="two.sided")$p.val pval3 <- pval2-pval op <- par(cex=1.5,lwd=2,bty="l",mgp=c(2.5,1,0),yaxs="i") ## barplot(hvec) plot(xvec,hvec,type="n",xlab="Number of pol stations with any seeds taken", ylab="",ylim=c(min(hvec),min(hvec)+diff(range(hvec))*1.1), axes=FALSE,xlim=c(-5,40)) rect(xvec-0.5,rep(min(hvec),length(xvec)),xvec+0.5,hvec, border="black",lwd=1, col=rep(c("darkgray","white","darkgray"), c(length(smvals), length(xvec)-length(smvals)-length(lgvals), length(lgvals)))) ## plot(xvec,hvec2,type="l",log="y") if (!badfonts) mtext("Probability",side=2,at=-10,line=3,las=0,cex=1.5) axis(side=1,at=seq(0,40,by=10)) axis(side=2,at=c(-5,-10,-15),labels=NA) prvals <- seq(0,10,by=2) axis(side=3,at=PR.to.x(prvals,vals=vals),labels=prvals, line=0.5,cex.axis=0.8) v0 <- vals[1,1] obs=v0 obsPR=3.62 par(xpd=NA) arrows(26,4.34,26,1.8,length=0.15,angle=15) ## 1.585 text(26,5.0, ## 4.82, paste("obs=",obsPR,sep=""),cex=0.5) par(xpd=TRUE) ## kluge mtext(side=3,at=PR.to.x(10,vals=vals),text="10", line=1.5,cex=0.8*1.5) mtext(side=3,at=20,line=3,"Probability ratio") par(las=1) invisible(lapply(seq(-5,-15,by=-5), function(x) { mtext(text=substitute(10^x,list(x=x)), side=2,at=x,line=par("mgp")[2],cex=1.5) })) v0 <- vals[1,1] x1 <- lgvals ##polygon(c(x1,rev(x1)),c(hvec[xvec>=v0],rep(par("usr")[3],length(x1))), ## col="gray",border=NA) x2 <- smvals v5 <- max(smvals) ## polygon(c(x2,rev(x2)),c(hvec[xvec<=v5],rep(par("usr")[3],length(x2))), ## col="gray",border=NA) ## lines(xvec,hvec) box(lwd=2) u1=par("usr")[3] u2=par("usr")[4] segments(v0,u1,v0,-2.7,lty=2) segments(v5,u1,v5,u2,lty=2) par(xpd=NA) text(v0,1,paste("observed\n=",v0,sep=""),pos=1) do.call("text",list(40,-8,scinot(pval,"expression",digits=2,pref="p="), cex=0.5)) do.call("text",list(-2,-2,scinot(pval3,"expression",digits=2,pref="p="), cex=0.5)) arrows(40,-9,36,-15,angle=15) arrows(-2,-3,1,-6,angle=15) par(op) ################################################### ### chunk number 5: ################################################### or = vals[1,1]*vals[2,2]/(vals[1,2]*vals[2,1]) ################################################### ### chunk number 6: ################################################### tprobs <- c(0.05,0.04) bprobs <- dbinom(sum(succ),prob=tprobs,size=sum(tot)) ml <- dbinom(sum(succ),prob=sum(succ)/sum(tot),size=sum(tot)) ################################################### ### chunk number 7: ################################################### op <- par(mfrow=c(2,1),las=1,lwd=2,cex=1.5,bty="l") par(mar=c(0,4,2,2)+0.2) curve(dbinom(sum(succ),size=sum(tot),prob=x),from=0.01,to=0.2,ylim=c(0,0.06), axes=FALSE, ylab="Likelihood",xlab="") ## xlab="Prob. any seeds\ntaken in one trial", axis(side=2,at=seq(0,0.06,by=0.03)) box() abline(v=sum(succ)/sum(tot),lty=2) par(mar=c(4,4,1,2)+0.2) corner.label2("a","topleft",inset=0.025) curve(dbinom(sum(succ),size=sum(tot),prob=x,log=TRUE),from=0.01,to=0.2, xlab="", axes=FALSE, ylab="Log-likelihood",ylim=c(-66,0)) mtext(side=1,expression(paste("P(seeds taken), ",p[s])), at=0.1,line=2.5,cex=1.5) axis(side=2,at=seq(0,-60,by=-20)) box() par(xpd=NA) corner.label2("b","topleft",inset=0.025) par(xpd=FALSE) axis(side=1) par(las=1) ## invisible(lapply(seq(-5,-25,by=-10), ## function(x) { mtext(text=substitute(10^x,list(x=x)), ## side=2,at=10^(x),line=par("mgp")[2],cex=1.5) })) #mtext(side=2,at=1e-15,"Likelihood",line=3.2,las=0,cex=1.5) #box() abline(v=sum(succ)/sum(tot),lty=2) par(op) ################################################### ### chunk number 8: ################################################### clip <- function(x) { min(1,max(0,x)) } plik1 <- function(r) { optimize(f=function(x) {-sum(dbinom(x=succ, size=tot,prob=c(x, clip(x*r)),log=TRUE))}, interval=c(0.01,0.4))$objective } ## use odds ratio instead of probability ratio plik1.OR <- function(r) { optimize(f=function(x) { or2 <- (x/(1-x))*r; p2 <- or2/(1+or2); -sum(dbinom(x=succ,size=tot,prob=c(x,p2),log=TRUE))}, interval=c(0.05,0.4))$objective } ## range of probability ratios to try ## changed range -- species 1 now has higher number taken rvec <- seq(0.01,1,length=300) rvec2 <- seq(1.5,7,length=100) ## negative log-likelihood lprof <- sapply(rvec,plik1) lprof2 <- sapply(1/rvec2,plik1) ## likelihood profile on odds ratio lprof.OR <- sapply(rvec,plik1.OR) fit.mle <- optim(fn=function(p) { -sum(dbinom(x=succ,size=tot,prob=c(p[1],clip(p[1]*p[2])), log=TRUE)) },par=c(0.2,2)) ##library(stats4) ##invisible(require(nlme,quietly=TRUE)) ##library(methods) invisible(require(bbmle,quietly=TRUE)) fit2.mle <- mle2(minuslogl=function(p,r) { -sum(dbinom(x=succ,size=tot,prob=c(p,clip(p*r)), log=TRUE)) },start=list(p=0.2,r=2)) ci.mle <- confint(fit2.mle,quietly=TRUE) ################################################### ### chunk number 9: ################################################### lprior <- function(prob,ratio) { 1 } ## improper prior in both directions plik2 <- function(prob,ratio) { exp(sum(dbinom(x=succ,size=tot,prob=c(prob,clip(prob*ratio)),log=TRUE))) } ## odds ratio: ## prob1 + odds ratio -> new odds ratio = prob1/(1-prob1)*oddsratio ## new probability = OR/(1+OR) plik2.OR <- function(prob,ratio) { or2 <- (prob/(1-prob))*ratio ## changed from *ratio to /ratio p2 <- or2/(1+or2) exp(sum(dbinom(x=succ,size=tot,prob=c(prob,p2),log=TRUE))) } pnum <- function(P) lprior(P[1],P[2])*plik2(P[1],P[2]) library(adapt) denom <- adapt(ndim=2,lower=c(0,0),upper=c(1,10),functn=pnum,eps=0.005)$value plik2B <- function(prob,ratio) { sapply(prob,plik2,ratio=ratio) } plik2B.OR <- function(prob,ratio) { sapply(prob,plik2.OR,ratio=ratio) } plik3 <- function(r) { integrate(plik2B,lower=0,upper=1,ratio=r)$value } plik3B <- function(p) { sapply(p,plik3) } plik3.OR <- function(r) { integrate(plik2B.OR,lower=0,upper=1,ratio=r)$value } ## what probability ratio corresponds to a given odds ratio??? ratiopost <- sapply(rvec,plik3) ratiopost.OR <- sapply(rvec,plik3.OR) ################################################### ### chunk number 10: eval=FALSE ################################################### ## plot(rvec,ratiopost) ## par(new=TRUE) ## plot(rvec,exp(-lprof),col=2,type="l") ## abline(v=p.mle[2]/p.mle[1],lty=2) ## ## plot(rvec,ratiopost.OR) ## par(new=TRUE) ## plot(rvec,exp(-lprofC),col=2,type="l") ## abline(v=1.9535) ################################################### ### chunk number 11: ################################################### densfac <- 1/sum(ratiopost)/diff(rvec)[1] ## find credible interval ... ## find lower and upper values for which prob dens = target value lims <- function(pdens) { lower <- uniroot(function(p) {plik3(p)*densfac-pdens}, interval=c(0.0001,rvec[which.max(ratiopost)]))$root upper <- uniroot(function(p) {plik3(p)*densfac-pdens}, interval=c(rvec[which.max(ratiopost)],3))$root c(lower,upper) } ## find area between target values limarea <- function(pdens) { intlim <- lims(pdens) (integrate(plik3B,intlim[1],intlim[2])$value)*densfac } ## find credible interval credint <- function(level) { u <- uniroot(function(x) {limarea(x)-level}, interval=c(1,1e-5)) intlim <- lims(u$root) c(intlim,plik3(intlim[1])*densfac,limarea(u$root)) } cred1 <- credint(0.95) bayesmean <- integrate(function(x)plik3B(x)*x,0.001,0.6)$value*densfac ## range covers plausible range bayesmode <- optimize(plik3B,interval=c(0.001,6),maximum=TRUE)$maximum ## bayeshyp <- integrate(plik3B,1,Inf)$value*densfac ################################################### ### chunk number 12: ################################################### op <- par(las=1,lwd=2,cex=1.5,bty="l",mar=c(5,4,4,4)+0.1) plot(rvec2,-lprof2,type="l",xlab="",ylab="Log-likelihood") mtext(side=1,"Ratio of (pol prob)/(psd prob)",line=2,cex=1.5) abline(h=-fit.mle$value,lty=2) ## plot(rvec2,exp(-lprof),type="l",col=2) abline(v=p.mle["pol"]/p.mle["psd"],lty=2) abline(v=1/ci.mle[2,],lty=3) abline(h=logLik(fit2.mle)+c(0,-qchisq(0.95,1)/2),lty=3) ##abline(v=1,lty=4) par(xpd=NA) text(1/coef(fit2.mle)[2],-4.5,"MLE") text(12,logLik(fit2.mle),"maximum\nlikelihood",adj=1) text(12,logLik(fit2.mle)-qchisq(0.95,1)/2-0.2,"95% cutoff",adj=1) text(1/ci.mle[2,],-4.5,c("upper","lower")) ##text(1,-17,"null",cex=1.5) ##arrows(c(1.3,2.0),-12,ci.mle[2,],-12,lwd=2) par(xpd=FALSE) ################################################### ### chunk number 13: ################################################### ## shade interior of credible interval? load("ch1bayes.RData") op <- par(cex=1.5,lwd=2,bty="l",mgp=c(2.5,1,0),las=1,yaxs="i") plot(rvec,ratiopost.norm,type="l", xlab="Ratio of (pol prob)/(psd prob)", ylab="Probability density",xlim=c(0.9,10),ylim=c(0,0.45)) abline(h=cred1["p"]) ##abline(v=cred1[1:2],lty=2) ##abline(v=1,lty=3) abline(v=bayesmode,lty=1) abline(v=bayesmean,lty=3) labht <- 0.45 ## text(bayesmean,1.0,"mean",cex=1.5) ## text(bayesmode,0.9,"mode",cex=1.5) rcred <- cred1["p"] x1 <- c(rvec[rveccred1["upr"]]) polygon(c(x2,rev(x2)),c(cred1["p"],ratiopost.norm[rvec>cred1["upr"]],rep(0,length(x2))), col="gray",border=NA) lines(rvec,ratiopost.norm) ## redraw ##text(1,1,"null",cex=1.5) ##text(1.7,0.1,"credible\ninterval",cex=1.5) ##arrows(c(1.4,2.0),0.08,cred1[1:2],0.08,lwd=2) arrht <- 0.03 arrows(cred1[1],arrht,cred1[2],arrht,angle=15,length=0.18,code=3) text(mean(cred1[1:2]),0.015,"credible interval",cex=0.65) par(xpd=NA) labgap <- 0.5 text(bayesmode-labgap,labht,"mode",adj=1) arrows(bayesmode-labgap*0.8,labht,bayesmode,labht,angle=15,length=0.15) text(bayesmean+labgap,labht,"mean",adj=0) arrows(bayesmean+labgap*0.8,labht,bayesmean,labht,angle=15,length=0.15) box() par(op) ################################################### ### chunk number 14: ################################################### detach(seedsub) ################################################### ### chunk number 15: ################################################### 2*8 sqrt(25) ################################################### ### chunk number 16: ################################################### x = sqrt(36) x ################################################### ### chunk number 17: eval=FALSE ################################################### ## install.packages("plotrix") ################################################### ### chunk number 18: eval=FALSE ################################################### ## install.packages(c("ellipse","plotrix")) ################################################### ### chunk number 19: eval=FALSE ################################################### ## install.packages("plotrix",repos=NULL) ################################################### ### chunk number 20: eval=FALSE ################################################### ## mypkgdir="c:/Documents and Settings/Bolker/Desktop/Rpkgs" ## # options(repos="http://cran.us.r-project.org") ## install.packages("plotrix",destdir=mypkgdir,lib=mypkgdir) ################################################### ### chunk number 21: eval=FALSE ################################################### ## library(plotrix,lib=mypkgdir) ################################################### ### chunk number 22: ################################################### ## TO DO: ignore packages? (weaver) r1 <- system("grep \"require(\" chap[0-9].Rnw chap[0-9][0-9].Rnw | sed -e 's/^.*require(\\([^),]*\\).*/\\1/'",intern=TRUE) r2 <- system("grep \"library(\" chap[0-9].Rnw chap[0-9][0-9].Rnw | sed -e 's/^.*library(\\([^),]*\\).*/\\1/'",intern=TRUE) r3 <- c(r1,r2) r3 <- unique(sort(r3[-grep("\\\\",r3)])) i1 <- installed.packages() omitpkgs <- rownames(i1)[!is.na(i1[,"Priority"])] omitpkgs <- c(omitpkgs,"Hmisc","weaver","pkg") r4 <- r3[!(r3 %in% omitpkgs)] ## OVERRIDE: add chron and gtools/gdata r4 <- c("adapt","bbmle","chron","coda","ellipse","emdbook","gplots","gtools","gdata", "MCMCpack","odesolve","plotrix","R2WinBUGS","reshape","rgl","scatterplot3d") cat(r4,"\n",fill=50) ################################################### ### chunk number 23: eval=FALSE ################################################### ## help.start() ################################################### ### chunk number 24: ################################################### set.seed(101) ################################################### ### chunk number 25: ################################################### frogs = c(1.1,1.3,1.7,1.8,1.9,2.1,2.3,2.4,2.5,2.8, 3.1,3.3,3.6,3.7,3.9,4.1,4.5,4.8,5.1,5.3) tadpoles = rnorm(n=20,mean=2*frogs,sd=0.5) ################################################### ### chunk number 26: ################################################### tadpoles ################################################### ### chunk number 27: ################################################### plot(frogs,tadpoles) abline(a=0,b=2) ################################################### ### chunk number 28: ################################################### plot(frogs,tadpoles) abline(a=0,b=2) ################################################### ### chunk number 29: ################################################### log_tadpoles = log(tadpoles) plot(frogs,log_tadpoles) ################################################### ### chunk number 30: ################################################### n=1:length(frogs) plot(n,frogs) ################################################### ### chunk number 31: ################################################### h = rev((1:6)/6) xr = seq(0,0.9,length=20) xpos = 1.2 op = par(xpd=TRUE,mar=c(0,0,0,5)) plot(0:1,range(h),type="n",axes=FALSE,ann=FALSE) points(xr,rep(h[1],20),pch=1:20) text(xpos,h[1],"pch: point type",adj=1) palette(gray((0:19)/19)) points(xr,rep(h[2],20),col=1:20,pch=16) palette("default") text(xpos,h[2],"col: point color",adj=1) points(xr,rep(h[3],20),cex=seq(0.01,2,length=20)) text(xpos,h[3],"cex: point size",adj=1) text(xr,rep(h[4],20),c(letters[1:5],LETTERS[1:5],0:9)) text(xpos,h[4],"text",adj=1) ##text(seq(0.1,0.9,length=20),rep(0.6,20),letters[1:20],cex=seq(0.01,2,length=20)) nl = 6 nlsp = 0.1 lend = 0.95 ##segments(seq(0,0.9,length=nl),rep(h[5],2*nl), ## seq(1/nl,1,length=nl),rep(h[5],2*nl), ## want to go from 0 to lend, in nl segments, ## using a fraction fr of the space for interspaces ## each segment + space takes up lend/nl ## fr = 0.7 seg1 = seq(0,by=lend/nl,length=nl) seg2 = seq(lend/nl,to=lend,by=lend/nl)-(lend*(1-fr)/nl) segments(seg1,rep(h[5],nl), seg2, rep(h[5],nl), lty=1:nl,lwd=3) text(xpos,h[5],"lty: line type",adj=1) nv = 10 segments(seq(0,1/(1+nlsp),length=nv),rep(h[6],nv), seq(0.05,0.95,length=nv),rep(h[6],nv), lwd=1:10) text(xpos,h[6],"lwd: line width",adj=1) par(op) ################################################### ### chunk number 32: ################################################### mean(tadpoles) sd(tadpoles) summary(tadpoles) ################################################### ### chunk number 33: ################################################### cor(frogs,tadpoles) ################################################### ### chunk number 34: ################################################### cor.test(frogs,tadpoles) ################################################### ### chunk number 35: eval=FALSE ################################################### ## help.search("correlation")