[R] how a latent state matrix is updated using package R2WinBUGS

[Uwe Ligges]

On 23.05.2012 17:20, Jean V Adams wrote:
> 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.

R2WinBUGS is just an interface package for R that exports data, calls 
WinBUGS and imports the results.
For BUGS related questions, i.e. in particular the correct specification 
of the model (file), I'd recommend to ask on the BUGS related mailing lists.

Note that the new version of WinBUGS is called OpenBUGS and the is a 
more "direct" R interface called BRugs.

Best,
Uwe Ligges


> 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 at 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.