## ---- echo=FALSE, warning=FALSE------------------------------------------ myecho <- TRUE myeval <- TRUE library(knitr) ## http://stackoverflow.com/questions/11062497/how-to-avoid-using-round-in-every-sexpr # knit_hooks$get("inline") ## Then customize it to your own liking inline_hook <- function(x) { if(is.numeric(x)) x <- round(x, 0) paste(as.character(x), collapse=", ") } knit_hooks$set(inline = inline_hook) opts_chunk$set( tidy=TRUE, size = 'small' ) ## ---- echo=myecho, eval=myeval------------------------------------------- xbill <- read.table("lure trials.csv",header=T,sep=",") attach(xbill) ## ---- echo=myecho, eval=myeval------------------------------------------- w <- 850 ## ---- echo=myecho, eval=myeval, comment=NA------------------------------- habitat <- factor(habitat) behavcode <- factor(behavcode) model1 <- glm(response~dist+numbirds+day+time+habitat+behavcode,family=binomial) summary(model1) anova(model1) ## ---- echo=myecho, eval=myeval, comment=NA------------------------------- model2 <- glm(response~dist,family=binomial) summary(model2) ## ---- echo=myecho, eval=myeval, fig=TRUE--------------------------------- fits <- fitted(model2) plot(dist,response,xlim=c(0,w),xlab="Distance from lure",ylab="Probability of response") dist2 <- 0:w phat <- predict.glm(model2,newdata=data.frame(dist=dist2),type="response") lines(dist2,phat) ## ---- echo=myecho, eval=myeval, fig=TRUE--------------------------------- pi_r <- dist2/sum(dist2) ## ---- echo=myecho, eval=myeval, comment=NA------------------------------- # Pa is prob of response, integrating r from 0 to w, assuming that phat is # a function of distance alone (no z's in model) (Pa <- sum(phat*pi_r)) ## ---- echo=myecho, eval=myeval------------------------------------------- detections <- read.table("mainsurveydetections.csv",sep=",",header=T) attach(detections) ## ---- echo=myecho, eval=myeval------------------------------------------- kmain <- length(point) a <- kmain*pi*(w/1000)^2 A <- 3505.8 # size in km2 of study area ## ---- echo=myecho, eval=myeval, comment=NA------------------------------- (Nscot <- sum(nscottish)*A/(Pa*a)) ## ---- echo=myecho, eval=myeval, comment=NA------------------------------- index <- 1:kmain nboot <- 999 bNscot <- vector(length=nboot) m <- length(dist) tindex <- 1:m bdist <- vector(length=m) bresponse <- vector(length=m) for (i in 1:nboot) { # bootstrap trials btindex <- sample(tindex,m,replace=TRUE) for (l in tindex) { bdist[l] <- dist[btindex[l]] bresponse[l] <- response[btindex[l]] } bmodel <- glm(bresponse~bdist,family=binomial) bphat <- predict.glm(bmodel,newdata=data.frame(bdist=dist2),type="response") bPa <- sum(bphat*pi_r) #bootstrap points rindex <- sample(index,kmain,replace=TRUE) bNscot[i] <- sum(nscottish[rindex])*A/(bPa*a) } #calculate bootstrap standard error seNscot <- sd(bNscot) #calculate bootstrap percentile confidence limits bNscot <- sort(bNscot) alpha.over.2 <- 0.025 loNscot <- bNscot[round(alpha.over.2*(nboot+1))] hiNscot <- bNscot[round((1-alpha.over.2)*(nboot+1))] (c(Nscot,seNscot,loNscot,hiNscot))