Dynamic Occupancy Model

MCMC algorithms generate Markov Chains which simulate the posterior distribution. The idea behind the Geweke diagnostic is to mimic a two-sample test of means and determine whether the first 10% and last 50% are of the same posterior distribution. If the chain has converged to a stationary distribution, the means of the two sections should be equal and the Geweke diagnostic approaches a normal distribution with a Z-score test statistic. The Geweke diagnostic plots are generated by dividing the chain into segments. The Z-scores are then calculated for the chain using all the iterations at first, then again for the iterations minus the first segment and then, minus the second segment and so on. If a large amount of the Z-scores fall outside the range \(\pm 1.96\), i.e. are extreme and fall within the tails of a standard normal distribution, the chain has not converged.

There are eleven parameters that required estimation in the dynamic occupancy model, including three regression coefficients for both the colonization, survival and occupancy (in the first year) probabilities. As well as, two regression coefficients for the detection probability.

For each the parameter within the chain, the \(\hat{R}\) value is inspected. The \(\hat{R}\) is the potential scale reduction factor and indicates whether within- and between-chain estimates agree. The notion behind \(\hat{R}\) is that if the chains have not reached the same target distribution, the within chain variance will be much lower than the variance of the chain combined. If chains have mixed well, the \(\hat{R}\) value should be below 1.1. Additionally, the trace plots are generated and Geweke diagnostic plots inspected. Regression coefficients have been summarised in main body of project.

R-hat values

All parameters have R-hat values below 1.05 and thus the chains have mixed well. Note the small effective sample sizes (ESS) for ecah parameter and slightly larger ESS for colonisation coefficients (beta.lgam).

Model3_no_res$summary
##                        mean           sd          2.5%           25%
## alpha[1]         -0.2769345  0.009738728    -0.2961775 -2.834579e-01
## alpha[2]          0.3300550  0.010910343     0.3087419  3.227497e-01
## beta.lpsi1[1]     0.0217937  0.131426967    -0.2274415 -6.801882e-02
## beta.lpsi1[2]    -4.7700144  0.404077545    -5.6082366 -5.035755e+00
## beta.lpsi1[3]    -4.1794014  0.343980517    -4.8920578 -4.406586e+00
## beta.lphi[1]      2.6533775  0.124487390     2.4234445  2.565762e+00
## beta.lphi[2]     -1.1352697  0.203230320    -1.5250771 -1.274949e+00
## beta.lphi[3]     -1.3624004  0.153605423    -1.6672004 -1.463932e+00
## beta.lgam[1]     -3.0042292  0.094642508    -3.1929013 -3.066366e+00
## beta.lgam[2]     -1.7221745  0.102815574    -1.9275287 -1.791120e+00
## beta.lgam[3]     -1.1584731  0.101319288    -1.3596852 -1.226734e+00
## deviance      62323.8827198 66.946578828 62194.6745025  6.227834e+04
##                         50%           75%         97.5%     Rhat n.eff overlap0
## alpha[1]      -2.768993e-01    -0.2703857    -0.2581184 1.000359  6267        0
## alpha[2]       3.299953e-01     0.3373572     0.3514584 1.000345  5572        0
## beta.lpsi1[1]  1.823052e-02     0.1071812     0.2906577 1.000555 25500        1
## beta.lpsi1[2] -4.761974e+00    -4.4912589    -4.0038617 1.004727   939        0
## beta.lpsi1[3] -4.175700e+00    -3.9365311    -3.5421334 1.004443   943        0
## beta.lphi[1]   2.648897e+00     2.7381319     2.9052417 1.002356   908        0
## beta.lphi[2]  -1.137203e+00    -1.0023862    -0.7185801 1.005661   385        0
## beta.lphi[3]  -1.363282e+00    -1.2585681    -1.0636318 1.005371   386        0
## beta.lgam[1]  -3.001940e+00    -2.9404624    -2.8190568 1.003417   616        0
## beta.lgam[2]  -1.721595e+00    -1.6526205    -1.5221614 1.001053  4461        0
## beta.lgam[3]  -1.157740e+00    -1.0899276    -0.9619152 1.001781  3163        0
## deviance       6.232364e+04 62368.6463199 62455.5221203 1.001473  1392        0
##                       f
## alpha[1]      1.0000000
## alpha[2]      1.0000000
## beta.lpsi1[1] 0.5556471
## beta.lpsi1[2] 1.0000000
## beta.lpsi1[3] 1.0000000
## beta.lphi[1]  1.0000000
## beta.lphi[2]  1.0000000
## beta.lphi[3]  1.0000000
## beta.lgam[1]  1.0000000
## beta.lgam[2]  1.0000000
## beta.lgam[3]  1.0000000
## deviance      1.0000000

Occupancy in first year

Trace plots of the parameter estimates for occupancy probability in the first year show no apparent anomalies, but there is evidence for mild serial correlation between successive draws, more so for \(\beta_1\) and \(\beta_2\).

psi1<-data.frame(Model3_no_res$sims.list$beta.lpsi1)
psi1$CHAIN<-rep(seq(1,3,by=1),each=8500)
psi1$ITER<-rep(seq(1,8500,by=1),3)
colnames(psi1)[c(1,2,3)]<-c("beta[0]","beta[1]","beta[2]")
psi1<-as.matrix(psi1)
psi1_mcmc = post_convert(psi1)

# traceplots
traplot(psi1_mcmc,"beta[0]",style="plain",col=c("deepskyblue2","firebrick1","slategray"),
        greek=T,main="",pty="s")

traplot(psi1_mcmc,"beta[1]",style="plain",col=c("deepskyblue2","firebrick1","slategray"),greek=T,main="")

traplot(psi1_mcmc,"beta[2]",style="plain",col=c("deepskyblue2","firebrick1","slategray"),greek=T,main="")

Some of the Geweke Diagnostic plots for the slope regression coefficients of occupancy probability (in the first year) are suggestive of convergence, given that only a few z-scores are found to be significant and are larger than the absolute value of 1.96. However, this does not hold for all chains and parameters as some plots have many z-scores outside the range of \(\pm 1.96\) which suggest that additional iterations are necessary.

geweke.plot(psi1_mcmc)

Survival

The trace plots for the regression coefficients of survival probability are not severely problematic but are suggestive of marginal mixing and autocorrelation among samples. These plots indicate that the chains are not traversing the distribution as rapidly as possible and in practice, one should increase the number of iterations specified.

surv<-data.frame(Model3_no_res$sims.list$beta.lphi)
surv$CHAIN<-rep(seq(1,3,by=1),each=8500)
surv$ITER<-rep(seq(1,8500,by=1),3)
colnames(surv)[c(1,2,3)]<-c("beta[0]","beta[1]","beta[2]")
surv<-as.matrix(surv)
surv_mcmc = post_convert(surv)

# traceplots
traplot(surv_mcmc,"beta[0]",style="plain",col=c("deepskyblue2","firebrick1","slategray"),
        greek=T,main="",pty="s")

traplot(surv_mcmc,"beta[1]",style="plain",col=c("deepskyblue2","firebrick1","slategray"),greek=T,main="")

traplot(surv_mcmc,"beta[2]",style="plain",col=c("deepskyblue2","firebrick1","slategray"),greek=T,main="")

Some of the Geweke Diagnostic plots for the slope regression coefficients of survival probability are suggestive of convergence, given that only a few z-scores are found to be significant and are larger than the absolute value of 1.96. However, this does not hold for all chains and parameters as some plots have many z-scores outside the range of \(\pm 1.96\) which suggest that additional iterations are necessary.

geweke.plot(surv_mcmc)

Colonisation

Given the larger effective sample sizes observed for the coefficients of the colonisation submodel, there is an overall improvement in the trace plots which show reduced serial correlation between subsequent samples and increased mixing of chains.

col<-data.frame(Model3_no_res$sims.list$beta.lgam)
col$CHAIN<-rep(seq(1,3,by=1),each=8500)
col$ITER<-rep(seq(1,8500,by=1),3)
colnames(col)[c(1,2,3)]<-c("beta[0]","beta[1]","beta[2]")
col<-as.matrix(col)
col_mcmc = post_convert(col)

# traceplots
traplot(col_mcmc,"beta[0]",style="plain",col=c("deepskyblue2","firebrick1","slategray"),
        greek=T,main="",pty="s")

traplot(col_mcmc,"beta[1]",style="plain",col=c("deepskyblue2","firebrick1","slategray"),greek=T,main="")

traplot(col_mcmc,"beta[2]",style="plain",col=c("deepskyblue2","firebrick1","slategray"),greek=T,main="")

All Geweke Diagnostic plots are indicative of convergence.

geweke.plot(col_mcmc)