# Read data pop.data=(read.csv("Logistic-Data-II.csv")) names(pop.data)=c("Year","Population Size", "Growth Rate") # Generate initial values for 3 chains inits=list( list(K=1200, r.adj=.5, r=0.21, sigma=.03), #chain 1 list(K=1300, r.adj=.1,r=0.12, sigma=.06), #chain 2 list(K=1250, r.adj=.7, r=0.2, sigma=.01) #1/0.01 #chain 3 ) n.xy = nrow(pop.data) #create vector of N's for calculating dNdt (See part d of Ex 1) N=seq(0,1500,10) #substitute 1 for 0 in the N vector N[N==0]=1 # Data must be a list data=list(n=n.xy, x=as.numeric(pop.data$"Population Size"), y=as.numeric(pop.data$"Growth Rate"), N=N, switch=1)# switch controls centering library(rjags) ##call to JAGS. # Run centered regression. jm=jags.model("Logistic example centered JAGS.R",data=data, inits,n.chains=length(inits), n.adapt = 5000) update(jm, n.iter=5000) zm=coda.samples(jm,variable.names=c("r.adj","r", "K" ,"sigma"),n.iter=10000) summary(zm) #Create dataframe for results df=as.data.frame(rbind(zm[[1]],zm[[2]],zm[[3]])) # Run uncentered regression. Note switch =0. data=list(n=n.xy, x=as.numeric(pop.data$"Population Size"),y=as.numeric(pop.data$"Growth Rate"),N=N, switch=0) jm=jags.model("Logistic example centered JAGS.R",data=data, inits,n.chains=length(inits), n.adapt = 5000) update(jm, n.iter=5000) zm1=coda.samples(jm,variable.names=c("K", "r.adj" ,"r","sigma","maxgrowth","dNdt"), n.iter=10000) zj=jags.samples(jm,variable.names=c("K", "r.adj" ,"r","sigma","maxgrowth","dNdt"), n.iter=10000) # Put results in dataframe df1=as.data.frame(rbind(zm1[[1]],zm1[[2]],zm1[[3]])) # Part C. Compare autocorrelation of centered and uncentered data for a.adj. ############################################################################ par(mfrow=c(2,1)) plot(acf(df$r.adj),main="Autocorrelation centered data", xlim=c(0,20)) plot(acf(df1$r.adj),main="Autocorrelation uncentered data", xlim=c(0,20)) # part D. ###################################################################### # Estimate and plot max.growth. max(df1$maxgrowth) mean(df1$maxgrowth) densityplot(df1$maxgrowth) # Plot median and quantiles of population growth rates, dN/dt b=summary(zj$dNdt,mean)$stat b=summary(zj$dNdt,quantile,c(.025,.5,.975))$stat par(mfrow=c(1,1)) plot(N,b[2,], xlab="N", ylab="Population growth rate",typ="l", ylim=c(-40,80)) lines(N,b[1,],lty="dashed") lines(N,b[3,],lty="dashed")