I'm trying to understand how a latent state matrix is updated by the MCMC
iterations in a WinBUGS model, using the package R2WinBUGS and an example
from Kery and Schaub's (2012) book, "Bayesian Population Analysis Using
WinBUGS". The example I'm using is 7.3.1. from a chapter on the
Cormack-Jolly-Seber model. Some excerpted code is included at the end of
this message; the full code is available at
http://www.vogelwarte.ch/downloads/files/publications/BPA/bpa-code.txt
The latent state of individual i on occasion t is stored in the z matrix
where rows index individuals (owls that are marked and released) and
columns index capture occasions. Each value in the matrix represents the
latent state for individual i at occasion t: z[i,t]=0 if individual i is
dead at time t, and =1 if individual i is alive at time t.
In the example, a matrix of known values for z is created from the capture
histories and provided as data; I will call it known.z. And a matrix of
initial values (where z is unknown) is also created and provided; I will
call it init.z. The dimensions of these two matrices are 250 individuals
by 6 capture occasions. However, if I save z as a parameter of interest,
the dimensions of the last saved z matrix from the last chain,
last.z <- cjs.c.c$last.values[[cjs.c.c$n.chains]]$z,
are 217 by 5. Why are the dimensions different? What happened to the
other 33 rows (individuals) and 1 column (occasion)? I thought perhaps
that the missing rows corresponded to those rows where all the latent
states are known, but that does not appear to be the case. There were no
individuals with 6 known latent states, and only 4 (not 33) with 5:
table(apply(!is.na(known.z), 1, sum))
0 1 2 3 4 5
162 39 27 14 4 4
Also, how can I verify that the known values of z are maintained in the
iterated z matrices? Even with the loss of 33 rows, those individuals
with 5 known 1=alive latent states don't seem to line up.
seq(known.z[,1])[apply(known.z[,-1], 1, paste, collapse="")=="11111"]
[1] 2 6 17 46
seq(last.z[,1])[apply(last.z, 1, paste, collapse="")=="11111"]
[1] 1 26 110 112 115 116
I would appreciate any insight you might have to offer. I am experienced
with R, but relatively inexperienced with WinBUGS and the R2WinBUGS
package. I am using R 2.15.0, R2WinBUGS 2.1-18, and WinBUGS 1.4.3 on
Windows 7.
Thanks.
Jean
> init.z[1:10, ]
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] NA 0 0 0 0 0
[2,] NA NA NA NA NA NA
[3,] NA NA 0 0 0 0
[4,] NA 0 0 0 0 0
[5,] NA 0 0 0 0 0
[6,] NA NA NA NA NA NA
[7,] NA 0 0 0 0 0
[8,] NA NA 0 0 0 0
[9,] NA NA NA 0 0 0
[10,] NA NA NA 0 0 0
> known.z[1:10, ]
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] NA NA NA NA NA NA
[2,] NA 1 1 1 1 1
[3,] NA 1 NA NA NA NA
[4,] NA NA NA NA NA NA
[5,] NA NA NA NA NA NA
[6,] NA 1 1 1 1 1
[7,] NA NA NA NA NA NA
[8,] NA 1 NA NA NA NA
[9,] NA 1 1 NA NA NA
[10,] NA 1 1 NA NA NA
> last.z[1:10, ]
[,1] [,2] [,3] [,4] [,5]
[1,] 1 1 1 1 1
[2,] 0 0 0 0 1
[3,] 0 0 0 0 1
[4,] 0 0 0 0 0
[5,] 0 0 0 0 1
[6,] 1 1 1 0 0
[7,] 0 0 0 0 0
[8,] 0 0 0 0 1
[9,] 0 0 0 0 0
[10,] 0 0 0 0 0
model {
# Priors and constraints
for (i in 1:nind){
for (t in f[i]:(n.occasions-1)){
phi[i,t] <- mean.phi
p[i,t] <- mean.p
} #t
} #i
mean.phi ~ dunif(0, 1) # Prior for mean survival
mean.p ~ dunif(0, 1) # Prior for mean recapture
# Likelihood
for (i in 1:nind){
# Define latent state at first capture
z[i,f[i]] <- 1
for (t in (f[i]+1):n.occasions){
# State process
z[i,t] ~ dbern(mu1[i,t])
mu1[i,t] <- phi[i,t-1] * z[i,t-1]
# Observation process
y[i,t] ~ dbern(mu2[i,t])
mu2[i,t] <- p[i,t-1] * z[i,t]
} #t
} #i
}
# Call WinBUGS from R
cjs.c.c <- bugs(
data = list(z = known.z, <snipped>),
inits = list(z = init.z, <snipped>),
parameters.to.save = c("z", <snipped>),
<snipped>
)
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
