library(jagsUI) #Requires jags to be installed (https://sourceforge.net/projects/mcmc-jags/) #Read in any of the data files dat<-read.csv("GADW_final_dat.csv") sink("full_model.jags") cat(" model { # Priors and constraints # mean parameter values (regression coefficients # note: JAGS uses precision (1/var), precision = 0.25, variance = 4, SD = 2 b0.mu ~ dnorm(0,0.25) # intercepts b1.mu ~ dnorm(0,0.25) # POND-change in pond density between YearT and Year T-1 b2.mu ~ dnorm(0,0.25) # PDSI- Palmer Drought Severity Index b3.mu ~ dnorm(0,0.25) # upland pair density (BPUP) b4.mu ~ dnorm(0,0.25) # intraspecific breeding pairs per pond (BPPO) b5.mu ~ dnorm(0,0.25) # SEOW-Short-eared owls per BBS route (vole index) b6.mu ~ dnorm(0, 0.25) # SEOW density in previous year (SOL1) b7.mu ~ dnorm(0, 0.25) # upland pair density in previous year (BPL1) b8.mu ~ dnorm(0, 0.25) # Year (linear trend) b9.mu ~ dnorm(0, 0.25) # Site Latitude b10.mu ~ dnorm(0, 0.25) #Site Longitude # Priors for SD to describe stratum-specific variation in regression coefficients b0.sigma ~ dunif(0,5) # numerically indexed as above b1.sigma ~ dunif(0,5) b2.sigma ~ dunif(0,5) b3.sigma ~ dunif(0,5) b4.sigma ~ dunif(0,5) b5.sigma ~ dunif(0,5) b6.sigma ~ dunif(0,5) b7.sigma ~ dunif(0,5) b8.sigma ~ dunif(0,5) # convert SD to precision (tau) by inverting and squaring b0.tau <- pow(b0.sigma,-2) b1.tau <- pow(b1.sigma,-2) b2.tau <- pow(b2.sigma,-2) b3.tau <- pow(b3.sigma,-2) b4.tau <- pow(b4.sigma,-2) b5.tau <- pow(b5.sigma,-2) b6.tau <- pow(b6.sigma,-2) b7.tau <- pow(b7.sigma,-2) b8.tau <- pow(b8.sigma,-2) # priors for random year and banding location (site) effects eta.yr.sigma ~ dunif(0,5) eta.site.sigma ~ dunif(0,5) eta.yr.tau <- pow(eta.yr.sigma,-2) eta.site.tau <- pow(eta.site.sigma,-2) # Likelihood for (i in 1:strata){ # where i indexes to the maximum number of strata b0[i] ~ dnorm(b0.mu,b0.tau) # generate strata specific regression coefficients b1[i] ~ dnorm(b1.mu,b1.tau) b2[i] ~ dnorm(b2.mu,b2.tau) b3[i] ~ dnorm(b3.mu,b3.tau) b4[i] ~ dnorm(b4.mu,b4.tau) b5[i] ~ dnorm(b5.mu,b5.tau) b6[i] ~ dnorm(b6.mu,b6.tau) b7[i] ~ dnorm(b7.mu,b7.tau) b8[i]~ dnorm(b8.mu,b8.tau) for (j in 1:maxyrs){ # where j indexes to number of years surveyed in each stratum eta.yr[i,j] ~ dnorm(0,eta.yr.tau) # generate random year effects within each stratum for (k in 1:maxsites){ # where k indexes to the maximum number of banding sites eta.site[i,j,k] ~ dnorm(0,eta.site.tau) # generate random banding site effects } # close sites loop } # close years loop } # close strata loop # Observation model for (k in 1:nrows){ # k indexes each line of data (from a banding site within a year within a stratum) # predict proportion of juveniles (Pjuv) based on stratum-specific random covariates, # fixed effects (lat, long), and random effects of year/stratum and site/year/stratum logit(Pjuv[k]) <- b0[stratum[k]] + b1[stratum[k]]*POND[k] + b2[stratum[k]]*PDSI[k] + b3[stratum[k]]*BPUP[k] + b4[stratum[k]]*BPPO[k]+ b5[stratum[k]]*SEOW[k] + b6[stratum[k]]*SOL1[k] + b7[stratum[k]]*BPL1[k] + b8[stratum[k]]*YEAR[k]+ b9.mu*Lat[k] + b9.mu*Long[k] + eta.yr[stratum[k],yr[k]] + eta.site[stratum[k],yr[k],site[k]] # compare observed juveniles to predicted juveniles based on above model HY_bands[k] ~ dbinom(Pjuv[k],bands[k]) # Juv ~ Binomial(P=Pjuv, N=Juv+Ad) #Fit assessments: Residuals predicted[k]<-Pjuv [k] #Age ratio predicted from data based on estimated coefficients resid[k]<-(HY_bands[k]/bands[k])-predicted[k] # Residuals: observed-predicted #Gelman-Rubin Fit statistic derivation chi2[k]<-pow(((HY_bands[k]/bands[k])-Pjuv[k]),2)/sqrt(Pjuv[k]+0.00001) # Chi squared: (Observed-predicted)2/sqrt(predicted + a very small number to make it non-zero) HYBands_new[k]~dbinom(Pjuv[k],bands[k]) # new dataset based on predicted age ratio and observed number of bands new.Age.rat[k]<-HYBands_new[k]/bands[k] # New age ratios chi2.new[k]<-pow(((HYBands_new[k]/bands[k])-Pjuv[k]),2)/sqrt(Pjuv[k]+0.00001) #Chi squared for new data } # Add discrepancy measures for dataset fit<-sum(chi2[]) # cumulative discrepancy for model relative to data. fit.new<-sum(chi2.new[])# cumulative discrepancy for model relative to new data. } # end bugs model ",fill=TRUE) sink() # Data bugs.data <- list(strata=max(dat$Strata_fac), stratum=dat$Strata_fac, yr=dat $Year_fac, maxyrs=max(dat $Year_fac), site=dat$Site_fac, maxsites=max(dat$Site_fac), PDSI=dat$ST_wPDSI, POND=dat$ST_PD1Diff, BPUP=dat$ST_BPDenUP, BPPO=dat$ST_BPDenPond, SEOW=dat$ST_SEOWT, SOL1=dat$ST_SEOWTL1, BPL1=dat$ST_TPop.Den.TM1, Lat=dat$ST.Lat, Long=dat$ST.Long, YEAR=dat$ST_Year, HY_bands=dat$HY_bands, bands=dat$TBands, nrows=nrow(dat)) # Parameters to save. parameters <- c("b0.mu","b0.sigma", "b1.mu", "b1.sigma", "b2.mu", "b2.sigma", "b3.mu", "b3.sigma", "b4.mu", "b4.sigma", "b5.mu", "b5.sigma", "b6.mu", "b6.sigma", "b7.mu", "b7.sigma", "b8.mu", "b8.sigma", "b9.mu","b10.mu","eta.yr.sigma", "eta.site.sigma","eta.yr", "eta.site", "resid", "b0", "b1", "b2", "b3", "b4", "b5", "b6", "b7","b8","fit","fit.new") out.model <- jagsUI(bugs.data, inits=NULL, parameters, "full_model.jags", n.chains=3, n.thin=2,n.iter=20000,n.burnin=100, n.adapt=2000) ##Bayesian p value good if close to 0.5, bad if close to 0 or 1 mean(out.model$sims.list$fit.new>out.model$sims.list$fit) # #==================================================== # Partitioning variance using the mallard model as an example # Marginal and Conditional R2GLMM per Johnson (2014) Methods in Ecology and Evolution 5 (944-946) # and Nakagawa & Schielzeth (2013) Methods in Ecology and Evolution 4 (133-142) # Their approach was designed to estimate a pseudo-R2 for generalized models with random intercepts # and slopes and can be adapted for our bayesian modelling approach. # Marginal and conditional R2 (Table 2 in Nakagawa and Schielzeth (2013) are calculated based # on variance components and distribution-specific adjustment ((π^2)/3 in the case of a logit model). # The components are: # σ2f= variance of fixed effects components, # σ2γ= variance of group-specific effects, # σ2α = variance of individual-specific effects, # σ2e = variance of residual #1. Calculate random effects variance (σ2γ +σ2α) based on the estimated random effect # sigmas saved by the model. This is a simplification of the covariance-matrix-based # method for finding the total random effects covariance provided by Johnson (2014) eqns 10-11. # The simple sum of the random effects variance was cross validated against Johnson’s variance-covariance matrix approach for all species. # List of random effect sigma objects from coda object vars.sigl<-list(out.model$sims.list$b0.sigma, out.model$sims.list$eta.yr.sigma, out.model$sims.list$eta.site.sigma, out.model$sims.list$b1.sigma, out.model$sims.list$b2.sigma, out.model$sims.list$b3.sigma, out.model$sims.list$b4.sigma, out.model$sims.list$b5.sigma, out.model$sims.list$b6.sigma, out.model$sims.list$b7.sigma, out.model$sims.list$b8.sigma ) # create empty matrix for each random effect to hold variance (sigma2) with one row for each # random effect (three intercept random effects (for strata, year within strata and banding site within year and stratum) and 8 slopes) RE.s2mat<-matrix(NA,ncol=1, nrow=11) rownames(RE.s2mat)<-c("x.S", "x.SY","x.SYB", "x.PDD1", "x.PDSI", "x.BPUP", "x.BPPO", "x.SEOW", "x.SOL1", "x.BPL1", "x.YEAR") # create list of variance (sigma2) of each random effect in matrix so that can be used in # matrix multiplication later: for(i in 1:nrow(RE.s2mat)){ RE.s2mat[i,1]<-(mean(vars.sigl[[i]]))^2 } Sl<-colSums(RE.s2mat) # Total Random Effect variance- same result as eqn 11 in Johnson 2014 #can be achieved with var-cov matrix as well # We can also look at subsets of the random effects variance Sl.t<-sum(RE.s2mat[c(2),]) #only year random effect Sl.s<-sum(RE.s2mat[c(1,3:10),]) #leave out year random effect Sl.site<-RE.s2mat[3,] #only site random effect #-------------------------------------------------------- # 2. Mean fixed effects variance(σ2f): # save coefficients out.coef<-c(out.model$mean$b0.mu, out.model$mean$b1.mu, out.model$mean$b2.mu, out.model$mean$b3.mu, out.model$mean$b4.mu, out.model$mean$b5.mu, out.model$mean$b6.mu, out.model$mean$b7.mu, out.model$mean$b8.mu, out.model$mean$b9.mu, out.model$mean$b10.mu ) #Matrix of 1s of length data to provide intercepts for matrix algebra below x<-matrix(1, ncol=1, nrow=nrow(dat)) # covariate data matrix xmat.FE<-cbind(x[,1], dat$ST_PD1Diff, dat$ST_wPDSI,dat$ST_BPDenUP,dat$ST_BPDenPond, dat$ST_SEOWT,dat$ST_SEOWTL1,dat$ST_TPop.Den.TM1, dat$ST_Year, dat$ST.Lat, dat$ST.Long) # Fixed effects variance following code from Johnson is the variance of product of coefficient matrix and data matrix Sf <- var(xmat.FE %*% out.coef) Sf.t<-var(xmat.FE[,1:7]%*%out.coef[1:7]) # Fixed effects var without lat long or year Sf.t2<-var(xmat.FE[,1:8]%*%out.coef[1:8]) # Fixed effects var with year but not lat long # 3. Variance of the residuals (σ2e): Se<- var(colMeans(out.model$sims.list$resid)) # Distribution-specific variance for binomial following Nakagawa & Schielzeth (2013) Table 2 Sd<-(3.14159^2)/3 # Total variance (denominator in Marginal and Conditional R2GLMM equations) Tvar<-Sl+Se+Sf+Sd # Marginal R2 (Just fixed effects in numerator); Johnson(2014) eqn 1 R2M<-Sf/Tvar # Conditional R2 (fixed and random effects in numerator); Johnson(2014) eqn 2 R2C<-(Sf+Sl)/Tvar