#### #### #### Part I: Preparing the workspace and exploring the data #### #### #### ### Part I: Preparing the data rm(list=ls()) getwd() library(rethinking) library(rstan) library(bbmle) library (MASS) dataf <- read.table("dataforstats_final_BB121613.txt", header=T) names (dataf) dataf$depth2 <- dataf$depth^2 dataf$r2 <- 0 dataf$r3 <- 0 dataf$r4 <- 0 dataf$r2[dataf$ranchno==2] <- 1 dataf$r3[dataf$ranchno==3] <- 1 dataf$r4[dataf$ranchno==4] <- 1 dataf$point <- as.factor(dataf$point) colSums(is.na(dataf)) dataf <- subset(dataf, select = -c(upland_elev_m,Point_elev,Point_elev,Maxdepth,maxdepth2,Distanceupland)) par(mfrow = c(2,4)) m0 <- mean(dataf$fishct) m1 <- mean(dataf$fishct[dataf$fishct >0]) Poiss_s0 <- rpois(length(dataf$fishct),m0) Poiss_s1 <- rpois(length(dataf$fishct),m1) Neg_bin_s0 <- rnbinom(length(dataf$fishct),mu=m0,size=0.58) Neg_bin_s1 <- rnbinom(length(dataf$fishct),mu=m1,size=0.58) Neg_bin_s0_100 <- rnbinom(length(dataf$fishct),mu=m0,size=100) Neg_bin_s1_100 <- rnbinom(length(dataf$fishct),mu=m0,size=0.19) b=seq(0,max(dataf$fishct)+30,1) hist(dataf$fishct,breaks = b,ylim=c(0,700),main="all data", xlim =c(0,20)) abline(v=m0,col="red") abline(v=m1,col="blue") hist(dataf$fishct[dataf$fishct>0],breaks = b,ylim=c(0,700),main="all data >0", xlim =c(0,20)) abline(v=m0,col="red") abline(v=m1,col="blue") hist(Poiss_s0,breaks = b,ylim=c(0,700),main="Poisson", xlim =c(0,20)) abline(v=m0,col="red") abline(v=m1,col="blue") hist(Poiss_s1,breaks = b,ylim=c(0,700),main="Poisson no zeros", xlim =c(0,20)) abline(v=m0,col="red") abline(v=m1,col="blue") hist(Neg_bin_s0,breaks = b,ylim=c(0,700),main="Negative Binomial 0.58", xlim =c(0,20)) abline(v=m0,col="red") abline(v=m1,col="blue") hist(Neg_bin_s1,breaks = b,ylim=c(0,700),main="Neg Bin no zeros 0.58", xlim =c(0,20)) abline(v=m0,col="red") abline(v=m1,col="blue") hist(Neg_bin_s0_100,breaks = b,ylim=c(0,700),main="Negative Binomial 100", xlim =c(0,20)) abline(v=m0,col="red") abline(v=m1,col="blue") hist(Neg_bin_s1_100,breaks = b,ylim=c(0,700),main="Negative Binomial 0.19", xlim =c(0,20)) abline(v=m0,col="red") abline(v=m1,col="blue") #### #### #### Part II: Model building and selection #### #### #### ## Model 1 Poisson m1 <- ulam( alist( fishct ~ dpois(lambda), log(lambda) <- a + b*depth + c*depth2 + d*r2 + e*r3 + f*r4, c(a,b,c,d,e,f) ~ dnorm(0,10) ), data = dataf,chains =3,iter=10000,warmup=1000,log_lik = TRUE ) precis(m1,digits=3,depth=2) m2 <- ulam( alist( fishct ~ dgampois(pbar,scale), log(pbar) <- a + b*depth + c*depth2 + d*r2 + e*r3 + f*r4, c(a,b,c,d,e,f) ~ dnorm(0,10), scale ~ dcauchy(0,2) ), data = dataf, constraints =list(theta ="lower=0"), iter=4000,warmup=1000,chains =3,log_lik = TRUE ) precis(m2,digits=3,depth=2) m3 <- ulam( alist( fishct ~ dzipois(pbar,lambda), logit(pbar) <- z1 + z2*depth, log(lambda) <- a + b*depth + c*depth2 + d*r2 + e*r3 + f*r4, c(a,b,c,d,e,f) ~ dnorm(0,10), c(z1,z2) ~ dnorm(0,10) ), data = dataf, iter=4000,warmup=1000,chains =3,log_lik = TRUE ) precis(m3,digits=3,depth=2) par(mfrow=c(1,1)) plot(m3) pairs(m3, horInd=c(1:8), verInd=c(1:8)) post3 <- extract.samples(m3) par(mfrow=c(2,5)) for(i in names(post3)){ dfpost3 <- data.frame(post3[i]) plot(as.numeric(rownames(dfpost3)), dfpost3[,1], type = "l", xlab = "Index", ylab = i) } compare(m1,m2,m3) ## Fig.5 ## x<-seq(0,max(dataf$depth),1) name <-c(" Ald","Bir","Pal","Wil") dataf$ranchn <- as.factor(dataf$ranchn) lr <- levels(dataf$ranchn) par(mfrow=c (1,4)) for (k in 1:4) { plot(dataf$depth,(dataf$fishct+1),type="n",log="y", ylim=c(1,100),xlim=c(1,60), main=name[k],xlab= "depth (cm)",ylab=("number of fish+1"),cex.lab=1.5,cex.main=1.7) depth_dat <- dataf$depth[dataf$ranchn==lr[k]] fish <- dataf$fishct[dataf$ranchn==lr[k]] t <- table(fish,depth_dat) dep <- unique(depth_dat) pez <- unique(fish) ord <- order(dep) dep <- dep[ord] orf <- order(pez) pez <- pez[orf] vv <- (lr==k)*1 for (i in 1: length(dep)) { for (j in 1: length(pez)) { points(dep[i],(pez[j]+1),pch=1,cex=log(t[j,i])+1,col="blue")}} mu <- link(m3,data=data.frame(depth=x,depth2=x^2,r2=vv[2],r3=vv[3],r4=vv[4])) mu.pbar.mean <- apply(mu$pbar,2,mean) mu.lambda.mean <- apply(mu$lambda,2,mean) mu.PI <- apply(mu$lambda,2,PI,prob=0.95) lines(x,mu.lambda.mean, ylim=c(0,1),col="red") shade(mu.PI,x)} ################################################ ## Fig.6 ## par(mfrow=c (1,4)) for (k in 1:4) { plot(seq(0,60,length=10),seq(0,1,length=10),type="n", ylim=c(0,1), xlim=c(0,60), main=name[k],xlab= "depth (cm)",ylab="P(zero)", cex.lab=1.5,cex.main=1.7) depth_dat <- dataf$depth[dataf$ranchn==lr[k]] fish <- dataf$fishct[dataf$ranchn==lr[k]] t <- table(fish,depth_dat) dep <- unique(depth_dat) pez <- unique(fish) ord <- order(dep) dep <- dep[ord] orf <- order(pez) pez <- pez[orf] vv <- (lr==k)*1 prebyden <- t[1,] tab <- colSums(t) probs <-prebyden/tab probs <- as.vector(probs) for (i in 2: length(fish)) { if(!is.na(probs[i])){ points(dep[i],probs[i],pch=16,col="blue")}} mu <- link(m3,data=data.frame(depth=x,depth2=x^2,r2=vv[2],r3=vv[3],r4=vv[4])) mu.pbar.mean <- apply(mu$pbar,2,mean) mu.lambda.mean <- apply(mu$lambda,2,mean) mu.PI <- apply(mu$pbar,2,PI,prob=0.95) lines(x,mu.pbar.mean,col="red") shade(mu.PI,x)}