library(coda) # data = number of recruits per species (rows) and grid cell (columns) # light = log(light) or NCI for all grid cells (vector) # abun = log(abundance) per species (vector) # apply function to data d1=date() res = metrop.RecLight.Gibbs(data,light,scale.burnin=1000,steps=7000,showstep=100,abun) d2=date() ######################### # function for two chains metrop.RecLight.Gibbs=function(data,light,scale.burnin=10000,steps=20000,showstep=100,abun) { whichuse = scale.burnin:steps # initial values and initial stepwidths for hyperparameters IntA1=SlopeA1=SdA1=IntB1=SlopeB1=SdB1=matrix(nrow=steps,ncol=3) IntA1[1,1]=-4 SlopeA1[1,1]=0 IntB1[1,1]=1 SlopeB1[1,1]=0 SdA1[1,1]=1 SdB1[1,1]=1 IntA1[1,3]=0.3 SlopeA1[1,3]=0.05 IntB1[1,3]=0.1 SlopeB1[1,3]=0.05 SdA1[1,3]=0.1 SdB1[1,3]=0.1 IntA2=SlopeA2=SdA2=IntB2=SlopeB2=SdB2=matrix(nrow=steps,ncol=3) IntA2[1,1]=-2 SlopeA2[1,1]=0.5 IntB2[1,1]=0 SlopeB2[1,1]=0.5 SdA2[1,1]=1 SdB2[1,1]=1 IntA2[1,3]=0.3 SlopeA2[1,3]=0.05 IntB2[1,3]=0.1 SlopeB2[1,3]=0.05 SdA2[1,3]=0.1 SdB2[1,3]=0.1 # initial values and initial stepwidths for parameters nospp=dim(data)[1] spa1=spb1=spk1=scalea1=scaleb1=scalek1=matrix(ncol=steps,nrow=nospp) spa2=spb2=spk2=scalea2=scaleb2=scalek2=matrix(ncol=steps,nrow=nospp) spa1[,1]=runif(nospp,-4,0) spb1[,1]=runif(nospp,-1,3) spk1[,1]=1 spa2[,1]=runif(nospp,-4,0) spb2[,1]=runif(nospp,-1,3) spk2[,1]=1 scalea1[,1]=scaleb1[,1]=scalea2[,1]=scaleb2[,1]=0.1 scalek1[,1]=scalek2[,1]=0.1 # prepare tables for covergence tests convergence = matrix(data=NA, ncol=3, nrow=nospp) converged = matrix(data=F, ncol=3, nrow=nospp) convergenceh = matrix(data=NA, ncol=1, nrow=6) convergedh = matrix(data=F, ncol=1, nrow=6) # loop over all steps for(i in 2:steps) { adjust=ifelse((i <= scale.burnin), T, F) # perform one Metropolis step for all hyperparamters IntA1[i,]=metrop1step(func=int.Gibbs,start.param=IntA1[i-1,1],scale.param=IntA1[i-1,3],Slope=SlopeA1[i-1,1], Sd=SdA1[i-1,1], sppar=spa1[,i-1], abun=abun, adjust.scale=adjust) SlopeA1[i,]=metrop1step(func=slope.Gibbs,start.param=SlopeA1[i-1,1],scale.param=SlopeA1[i-1,3],Int=IntA1[i,1], Sd=SdA1[i-1,1], sppar=spa1[,i-1], abun=abun, adjust.scale=adjust) SdA1[i,]=metrop1step(func=sd.Gibbs,start.param=SdA1[i-1,1],scale.param=SdA1[i-1,3],Int=IntA1[i,1],Slope=SlopeA1[i,1], sppar=spa1[,i-1], abun=abun, adjust.scale=adjust) IntB1[i,]=metrop1step(func=int.Gibbs,start.param=IntB1[i-1,1],scale.param=IntB1[i-1,3],Slope=SlopeB1[i-1,1], Sd=SdB1[i-1,1], sppar=spb1[,i-1], abun=abun, adjust.scale=adjust) SlopeB1[i,]=metrop1step(func=slope.Gibbs,start.param=SlopeB1[i-1,1],scale.param=SlopeB1[i-1,3],Int=IntB1[i,1], Sd=SdB1[i-1,1], sppar=spb1[,i-1], abun=abun, adjust.scale=adjust) SdB1[i,]=metrop1step(func=sd.Gibbs,start.param=SdB1[i-1,1],scale.param=SdB1[i-1,3],Int=IntB1[i,1],Slope=SlopeB1[i,1], sppar=spb1[,i-1], abun=abun, adjust.scale=adjust) # 2nd chain IntA2[i,]=metrop1step(func=int.Gibbs,start.param=IntA2[i-1,1],scale.param=IntA2[i-1,3],Slope=SlopeA2[i-1,1], Sd=SdA2[i-1,1], sppar=spa2[,i-1], abun=abun, adjust.scale=adjust) SlopeA2[i,]=metrop1step(func=slope.Gibbs,start.param=SlopeA2[i-1,1],scale.param=SlopeA2[i-1,3],Int=IntA2[i,1], Sd=SdA2[i-1,1], sppar=spa2[,i-1], abun=abun, adjust.scale=adjust) SdA2[i,]=metrop1step(func=sd.Gibbs,start.param=SdA2[i-1,1],scale.param=SdA2[i-1,3],Int=IntA2[i,1],Slope=SlopeA2[i,1], sppar=spa2[,i-1], abun=abun, adjust.scale=adjust) IntB2[i,]=metrop1step(func=int.Gibbs,start.param=IntB2[i-1,1],scale.param=IntB2[i-1,3],Slope=SlopeB2[i-1,1], Sd=SdB2[i-1,1], sppar=spb2[,i-1], abun=abun, adjust.scale=adjust) SlopeB2[i,]=metrop1step(func=slope.Gibbs,start.param=SlopeB2[i-1,1],scale.param=SlopeB2[i-1,3],Int=IntB2[i,1], Sd=SdB2[i-1,1], sppar=spb2[,i-1], abun=abun, adjust.scale=adjust) SdB2[i,]=metrop1step(func=sd.Gibbs,start.param=SdB2[i-1,1],scale.param=SdB2[i-1,3],Int=IntB2[i,1],Slope=SlopeB2[i,1], sppar=spb2[,i-1], abun=abun, adjust.scale=adjust) # convergence test for hyperparameters if(i%%showstep==0 & i>=500) { g <- gelman.diag(mcmc.list(mcmc(IntA1[((i-500)+1):i,1]), mcmc(IntA2[((i-500)+1):i,1])), confidence = 0.95, transform=T, autoburnin=F) if (g[[1]][1] < 1.1 & !convergedh[1]) { convergenceh[1] <- i convergedh[1] <- T } g <- gelman.diag(mcmc.list(mcmc(SlopeA1[((i-500)+1):i,1]), mcmc(SlopeA2[((i-500)+1):i,1])), confidence = 0.95, transform=T, autoburnin=F) if (g[[1]][1] < 1.1 & !convergedh[2]) { convergenceh[2] <- i convergedh[2] <- T } g <- gelman.diag(mcmc.list(mcmc(SdA1[((i-500)+1):i,1]), mcmc(SdA2[((i-500)+1):i,1])), confidence = 0.95, transform=T, autoburnin=F) if (g[[1]][1] < 1.1 & !convergedh[3]) { convergenceh[3] <- i convergedh[3] <- T } g <- gelman.diag(mcmc.list(mcmc(IntB1[((i-500)+1):i,1]), mcmc(IntB2[((i-500)+1):i,1])), confidence = 0.95, transform=T, autoburnin=F) if (g[[1]][1] < 1.1 & !convergedh[4]) { convergenceh[4] <- i convergedh[4] <- T } g <- gelman.diag(mcmc.list(mcmc(SlopeB1[((i-500)+1):i,1]), mcmc(SlopeB2[((i-500)+1):i,1])), confidence = 0.95, transform=T, autoburnin=F) if (g[[1]][1] < 1.1 & !convergedh[5]) { convergenceh[5] <- i convergedh[5] <- T } g <- gelman.diag(mcmc.list(mcmc(SdB1[((i-500)+1):i,1]), mcmc(SdB2[((i-500)+1):i,1])), confidence = 0.95, transform=T, autoburnin=F) if (g[[1]][1] < 1.1 & !convergedh[6]) { convergenceh[6] <- i convergedh[6] <- T } write.table(convergenceh, file="convergenceh.txt", sep="\t") write.table(convergedh, file="convergedh.txt", sep="\t") } # loop over species for(j in 1:nospp) { # perform one Metropolis step for all parameters nextv=metrop1step(func=spa.Gibbs,start.param=spa1[j,i-1],scale.param=scalea1[j,i-1],spb=spb1[j,i-1],spk=spk1[j,i-1], MuA=IntA1[i,1]+SlopeA1[i,1]*abun[j], SdA=SdA1[i,1], light=light, obs=data[j,], adjust.scale=adjust) spa1[j,i]=nextv[1] scalea1[j,i]=nextv[3] nextv=metrop1step(func=spb.Gibbs,start.param=spb1[j,i-1],scale.param=scaleb1[j,i-1],spa=spa1[j,i],spk=spk1[j,i-1], MuB=IntB1[i,1]+SlopeB1[i,1]*abun[j],SdB=SdB1[i,1], light=light, obs=data[j,], adjust.scale=adjust) spb1[j,i]=nextv[1] scaleb1[j,i]=nextv[3] nextv=metrop1step(func=spk.Gibbs,start.param=spk1[j,i-1],scale.param=scalek1[j,i-1],spa=spa1[j,i],spb=spb1[j,i], light=light, obs=data[j,], adjust.scale=adjust) spk1[j,i]=nextv[1] scalek1[j,i]=nextv[3] # 2nd chain nextv=metrop1step(func=spa.Gibbs,start.param=spa2[j,i-1],scale.param=scalea2[j,i-1],spb=spb2[j,i-1],spk=spk2[j,i-1], MuA=IntA2[i,1]+SlopeA2[i,1]*abun[j], SdA=SdA2[i,1], light=light, obs=data[j,], adjust.scale=adjust) spa2[j,i]=nextv[1] scalea2[j,i]=nextv[3] nextv=metrop1step(func=spb.Gibbs,start.param=spb2[j,i-1],scale.param=scaleb2[j,i-1],spa=spa2[j,i],spk=spk2[j,i-1], MuB=IntB2[i,1]+SlopeB2[i,1]*abun[j],SdB=SdB2[i,1], light=light, obs=data[j,], adjust.scale=adjust) spb2[j,i]=nextv[1] scaleb2[j,i]=nextv[3] nextv=metrop1step(func=spk.Gibbs,start.param=spk2[j,i-1],scale.param=scalek2[j,i-1],spa=spa2[j,i],spb=spb2[j,i], light=light, obs=data[j,], adjust.scale=adjust) spk2[j,i]=nextv[1] scalek2[j,i]=nextv[3] # convergence test parameters if(i%%showstep==0 & i>=500) { g <- gelman.diag(mcmc.list(mcmc(spa1[j,((i-500)+1):i]), mcmc(spa2[j,((i-500)+1):i])), confidence = 0.95, transform=T, autoburnin=F) if (g[[1]][1] < 1.1 & !converged[j,1]) { convergence[j,1] <- i converged[j,1] <- T } g <- gelman.diag(mcmc.list(mcmc(spb1[j,((i-500)+1):i]), mcmc(spb2[j,((i-500)+1):i])), confidence = 0.95, transform=T, autoburnin=F) if (g[[1]][1] < 1.1 & !converged[j,2]) { convergence[j,2] <- i converged[j,2] <- T } g <- gelman.diag(mcmc.list(mcmc(spk1[j,((i-500)+1):i]), mcmc(spk2[j,((i-500)+1):i])), confidence = 0.95, transform=T, autoburnin=F) if (g[[1]][1] < 1.1 & !converged[j,3]) { convergence[j,3] <- i converged[j,3] <- T } write.table(convergence, file="convergence.txt", sep="\t") write.table(converged, file="converged.txt", sep="\t") } } cat(i, "\n") } # structure of output thin = (whichuse%%4==0) hyperparams = data.frame(IntA1[thin,1],SlopeA1[thin,1],SdA1[thin,1],IntB1[thin,1],SlopeB1[thin,1],SdB1[thin,1],IntA2[thin,1],SlopeA2[thin,1],SdA2[thin,1],IntB2[thin,1],SlopeB2[thin,1],SdB2[thin,1]) colnames(hyperparams) = c("IntA1","SlopeA1","SdA1","IntB1","SlopeB1","SdB1","IntA2","SlopeA2","SdA2","IntB2","SlopeB2","SdB2") out = list(spa1[,thin], spa2[,thin], spb1[,thin], spb2[,thin], spk1[,thin], spk2[,thin], hyperparams) return(out) } ################################################## # Gibbs probability functions for hyperparameters int.Gibbs=function(Int,Slope,sppar,Sd,abun) { if(Sd<=0) return(-Inf) llike=dnorm(sppar,mean=Int+Slope*abun,sd=Sd,log=T) return(sum(llike)+dnorm(Int, mean=0, sd=100, log=T)) } slope.Gibbs=function(Slope,Int,sppar,Sd,abun) { if(Sd<=0) return(-Inf) llike=dnorm(sppar,mean=Int+Slope*abun,sd=Sd,log=T) return(sum(llike)+dnorm(Slope, mean=0, sd=100, log=T)) } sd.Gibbs=function(Sd,Int,Slope,sppar,abun) { if(Sd<=0) return(-Inf) llike=dnorm(sppar,mean=Int+Slope*abun,sd=Sd,log=T) return(sum(llike)+dlnorm(Sd, meanlog=0, sdlog=1, log=T)) } # Gibbs probability functions for parameters spa.Gibbs=function(spa,spb,spk,MuA,SdA,obs,light) { pred.nr = exp(spa+spb*light) llike <- sum(dnbinom(obs, mu=pred.nr, size=spk, log=T))+dnorm(spa,mean=MuA,sd=SdA,log=T) return(llike) } spb.Gibbs=function(spb,spa,spk,MuB,SdB,obs,light) { pred.nr = exp(spa+spb*light) llike <- sum(dnbinom(obs, mu=pred.nr, size=spk, log=T))+dnorm(spb,mean=MuB,sd=SdB,log=T) return(llike) } spk.Gibbs=function(spk,spa,spb,obs,light) { if(spk<0) return(-Inf) pred.nr = exp(spa+spb*light) llike <- sum(dnbinom(obs, mu=pred.nr, size=spk, log=T))+dunif(spk,min=0,max=100,log=T) return(llike) } ########################################################################################################## # FUNCTION metrop1step # # A generic routine for taking a single Metropolis step given any probability function, func, # an initial parameter, start.param, and a scaling parameter for the step size, scale.param. It returns # the next parameter value as well as an acceptance indicator (0 if the value was accepted, 1 if # rejected), the new scale parameter. If adjust.scale=T, the scale parameter # is adjusted to keep the acceptance rate at around 0.234 # adjust.scale should be T only for steps=origlike) accept = T #return(c(newval,0)) else { likeratio=exp(newlike-origlike) if(runif(1)