Accessing MCMC chain values/real estimates
Posted: Thu Oct 30, 2025 8:24 amI ran a multistate model within Program MARK using MCMC estimation (as Multistrata in RMark doesn't appear to support MCMC) and am now trying to get the individual chain values as real estimates via the .BIN file from Program MARK. There used to be a button in Program MARK that would give all the individual beta and real estimates for each iteration of the chains but it seems to be gone (If there is a way to get an old version of Program MARK that will output a CSV of the chain values - please let me know!). As an alternative, I am now trying to get the real estimate via the .BIN file in RMark.
Model Overview: 4 digit capture histories, 5 states
Parameterized as S(state*time)p(state*time)psi(state to state) meaning there should be 50 parameters (output later said 43 parameters estimated).
mlogit used for psi, S and p default (logit)
Ran 3 chains but other MCMC settings were default (uninformed priors, 4000 tuning, 1000 burnin, 10000 stored).
The summary of the chains looks fine (convergence, reasonable estimates, etc), and now I'm trying to access the individual chain values as real parameter estimates.
I initially tried to determine if there were 50 vs. 43 parameters in the BIN file given when I saw in the output file.
Assuming that the columns are parameters and that the chains are stacked vertically, I used this code and divided by 43 and 50
con <- file("MCMC.BIN", "rb")
raw_vals <- readBin(con, numeric(), n = 1e7)
Dividing by 50 resulted in an even number of values per column (31,000 rows per column). I would have expected only 30,000 rows per column (10,000 iterations x 3 chains) but I assume that the 1000 burnin samples were also stored.
I was able to make and export this as a csv. using the following code.
mcmc_matrix <- matrix(raw_vals, ncol = 50, byrow = TRUE)
mcmc_matrix
library(coda)
mcmc_all <- mcmc(mcmc_matrix)
mcmc_df <- as.data.frame(mcmc_all)
write.csv(mcmc_df, "full_mcmc_chains.csv", row.names = FALSE)
The dataframe I created appeared to be beta estimates, so I then attempted to convert to real estimates.
I applied inverse logit to the first 30 columns (my S and p parameters)
inv_logit <- function(x) 1 / (1 + exp(-x))
However, when I generate means of these transformed columns, the means seem off from the real parameter estimates provided in the Program MARK output. For example, the 3 chains in the output provide mean estimates of 0.7-0.65 for the first survival parameter. The mean of the column from my conversion of the .BIN file has an mean of 0.5.
Does what I've done so far seem correct or is there an easier way to be getting the chain values as real estimates?
If it does seem correct, suggestions on the conversion of the psi columns to real estimates?
Any help would be appreciated. Thanks all
NOTE:I did also try create.mark.mcmc but had no luck with it, even when I used a 32x version of R to be sure it matched Program MARK.
Model Overview: 4 digit capture histories, 5 states
Parameterized as S(state*time)p(state*time)psi(state to state) meaning there should be 50 parameters (output later said 43 parameters estimated).
mlogit used for psi, S and p default (logit)
Ran 3 chains but other MCMC settings were default (uninformed priors, 4000 tuning, 1000 burnin, 10000 stored).
The summary of the chains looks fine (convergence, reasonable estimates, etc), and now I'm trying to access the individual chain values as real parameter estimates.
I initially tried to determine if there were 50 vs. 43 parameters in the BIN file given when I saw in the output file.
Assuming that the columns are parameters and that the chains are stacked vertically, I used this code and divided by 43 and 50
con <- file("MCMC.BIN", "rb")
raw_vals <- readBin(con, numeric(), n = 1e7)
Dividing by 50 resulted in an even number of values per column (31,000 rows per column). I would have expected only 30,000 rows per column (10,000 iterations x 3 chains) but I assume that the 1000 burnin samples were also stored.
I was able to make and export this as a csv. using the following code.
mcmc_matrix <- matrix(raw_vals, ncol = 50, byrow = TRUE)
mcmc_matrix
library(coda)
mcmc_all <- mcmc(mcmc_matrix)
mcmc_df <- as.data.frame(mcmc_all)
write.csv(mcmc_df, "full_mcmc_chains.csv", row.names = FALSE)
The dataframe I created appeared to be beta estimates, so I then attempted to convert to real estimates.
I applied inverse logit to the first 30 columns (my S and p parameters)
inv_logit <- function(x) 1 / (1 + exp(-x))
However, when I generate means of these transformed columns, the means seem off from the real parameter estimates provided in the Program MARK output. For example, the 3 chains in the output provide mean estimates of 0.7-0.65 for the first survival parameter. The mean of the column from my conversion of the .BIN file has an mean of 0.5.
Does what I've done so far seem correct or is there an easier way to be getting the chain values as real estimates?
If it does seem correct, suggestions on the conversion of the psi columns to real estimates?
Any help would be appreciated. Thanks all
NOTE:I did also try create.mark.mcmc but had no luck with it, even when I used a 32x version of R to be sure it matched Program MARK.