Effective Sample Size and Monte Carlo Error in Joint Models

Why Draw Count Is Not the Whole Story

A joint model may produce thousands of post-warmup MCMC draws for an association parameter, a baseline-hazard parameter, or a random-effects variance. Those draws are not generally independent: nearby iterations can carry similar information. Consequently, the number of recorded draws alone does not tell us how precisely the simulation estimates a posterior mean or quantile.

Effective sample size (ESS) expresses the information in a correlated MCMC sequence as an approximate number of independent draws. Monte Carlo standard error (MCSE) estimates the simulation uncertainty in a reported posterior summary. These quantities answer a different question from posterior standard deviation or credible-interval width. Posterior uncertainty describes what remains uncertain about the model parameter; MCSE describes how much the finite simulation may have affected the numerical summary.

Definition: For a scalar posterior quantity $\theta$, let $S=MN$ be the total number of post-warmup draws from $M$ chains, each with $N$ retained iterations. If $\rho_k$ is the lag-$k$ autocorrelation of the stationary sequence for $\theta$, its integrated autocorrelation time is $$ \tau_\theta = 1 + 2\sum_{k=1}^{\infty}\rho_k. $$ The corresponding effective sample size is $$ \operatorname{ESS}_\theta = \frac{S}{\tau_\theta}. $$

In applications, the infinite sum is estimated from the observed chains using a suitable truncation or spectral-variance procedure. The precise estimator varies across software, especially in how it combines multiple chains and handles noisy estimates at larger lags. ESS is therefore an estimate, not a literal count of independent parameter values. Under positive serial dependence it is often smaller than $S$; negative dependence can make it larger.

The formula applies to a particular scalar quantity. A joint model has many parameters and derived summaries, and each may have a different dependence pattern. For example, the ESS for a longitudinal slope need not be the ESS for the association parameter, and neither necessarily equals the ESS for a hazard ratio calculated from that parameter. Compute or request ESS for the quantity being summarized.

From Effective Sample Size to Monte Carlo Error

Suppose the chains are exploring the target posterior distribution adequately and the sequence for $\theta$ is stationary, with finite variance and dependence that decays sufficiently for a central limit approximation. Then the variance of the simulated sample mean is approximately the posterior variance divided by ESS. If $s_\theta$ is the posterior standard deviation estimated from the draws, the estimated MCSE of the posterior mean is

\widehat{\operatorname{MCSE}}(\bar{\theta}) \approx \frac{s_\theta}{\sqrt{\operatorname{ESS}_\theta}}.

This relation can be justified directly. Write the retained sequence as $\theta_1,\ldots,\theta_S$, with posterior variance $\sigma_\theta^2$ and lag-$k$ autocorrelation $\rho_k$. The variance of its sample mean is

\operatorname{Var}(\bar{\theta}) = \frac{\sigma_\theta^2}{S} \left[1+2\sum_{k=1}^{S-1}\left(1-\frac{k}{S}\right)\rho_k\right].

When the chain is long enough that the weighted sum is well approximated by the integrated autocorrelation time, this becomes $\sigma_\theta^2\tau_\theta/S$. Since $\operatorname{ESS}_\theta=S/\tau_\theta$, the variance is approximately $\sigma_\theta^2/\operatorname{ESS}_\theta$. Taking the square root and replacing the unknown posterior standard deviation with $s_\theta$ gives the MCSE formula.

Under the same approximation, a rough 95% Monte Carlo error range for an estimated posterior mean is the estimate plus or minus $1.96$ MCSE. This is an interval for simulation error in the numerical estimate, not a 95% posterior credible interval for $\theta$. In practice, judge MCSE relative to the precision needed for the scientific report—for example, relative to the posterior standard deviation, the displayed decimal places, or a meaningful effect-size threshold.

Worked Example: MCSE for a Joint-Model Association Mean

Suppose four chains each retain 1,000 draws for an association parameter $\alpha$, which links the longitudinal marker trajectory to the event hazard. Thus $S=4(1{,}000)=4{,}000$. For illustration, suppose the estimated average autocorrelations at lags 1 through 4 are $0.50, 0.25, 0.10,$ and $0.05$, and the estimator's positive-sequence truncation retains these four lags and no later ones. The estimated integrated autocorrelation time is $$ \hat{\tau}_{\alpha}=1+2(0.50+0.25+0.10+0.05)=1+2(0.90)=2.80. $$ Therefore, $$ \widehat{\operatorname{ESS}}_{\alpha}=\frac{4{,}000}{2.80}=1{,}428.6. $$

Suppose the posterior mean of $\alpha$ is $0.84$ per one-unit increase in the marker, and its posterior standard deviation is $0.30$ per marker unit. The estimated MCSE of the mean is $$ \widehat{\operatorname{MCSE}}(\bar{\alpha}) \approx \frac{0.30}{\sqrt{1{,}428.6}} = \frac{0.30}{37.80} =0.00794 $$ per marker unit. The approximate 95% Monte Carlo error range is $0.84\pm1.96(0.00794)$, or $0.824$ to $0.856$ per marker unit. This suggests the simulation error in the mean is small relative to the posterior standard deviation of $0.30$. It does not say that the posterior uncertainty in $\alpha$ is only about $0.008$.

Derived Effects and Posterior Quantiles

Joint-model estimates are often reported on a transformed scale. For instance, if $\alpha$ is a log-hazard association coefficient, $\exp(\alpha)$ is a hazard ratio for a one-unit marker increase, under the model's specified association structure. For a nonlinear transformation $g(\theta)$, calculate $g$ for each posterior draw and summarize those transformed draws. The transformed sequence has its own dependence and therefore its own ESS; it should not automatically inherit the ESS of $\theta$.

Worked Example: MCSE for a Posterior Mean Hazard Ratio

Suppose a hypothetical joint-model analysis retains 4,000 draws of $\alpha$. The posterior mean and standard deviation of the draw-by-draw hazard ratios $\exp(\alpha)$ are 1.42 and 0.38, respectively. Assume the ESS calculated for the transformed hazard-ratio draws is 1,250. The MCSE for the posterior mean hazard ratio is $$ \widehat{\operatorname{MCSE}}\{\overline{\exp(\alpha)}\} \approx \frac{0.38}{\sqrt{1{,}250}} = \frac{0.38}{35.36} =0.0107. $$

The approximate 95% Monte Carlo error range for the estimated mean hazard ratio is $1.42\pm1.96(0.0107)$, or $1.399$ to $1.441$. The estimate is dimensionless. This range describes numerical simulation error around the posterior mean hazard ratio, not a credible interval for the hazard ratio itself. The draw-by-draw transformation matters: in general, the mean of $\exp(\alpha)$ is not equal to the exponential of the mean of $\alpha$.

For a posterior quantile, the mean-based MCSE formula is not the appropriate direct calculation. If $q_p$ is the posterior $p$-quantile and the posterior density $f$ is positive and smooth near it, an independent sample of size $n$ has approximate quantile variance $p(1-p)/\{n f(q_p)^2\}$. For MCMC, replace $n$ with an effective sample size appropriate to estimating that quantile. In particular, the relevant dependence is that of the indicator sequence $I(\theta_t\leq q_p)$, not necessarily the dependence summarized by the ordinary mean ESS for $\theta$. Quantile-specific ESS or a quantile MCSE reported by reliable software is preferable to substituting a generic ESS without checking its definition.

Worked Example: MCSE for an Upper Posterior Quantile

Suppose the reported 95th posterior percentile of a log-hazard association coefficient is $q_{0.95}=0.62$ per marker unit. For illustration, assume the estimated posterior density at this quantile is $0.20$ per coefficient unit and the quantile-specific ESS is 2,500. The approximate MCSE is $$ \widehat{\operatorname{MCSE}}(q_{0.95}) \approx \frac{\sqrt{0.95(1-0.95)/2{,}500}}{0.20} =\frac{\sqrt{0.000019}}{0.20} =\frac{0.004359}{0.20} =0.0218 $$ per marker unit. Thus the reported percentile has appreciable simulation uncertainty on the scale of a few hundredths of a coefficient unit. The calculation depends on the density estimate and quantile-specific ESS, so it is an approximation rather than a guarantee of accuracy.

Assessing Precision and Recognizing Failure Modes

ESS and MCSE should be interpreted alongside chain diagnostics. Between-chain agreement, assessed with $\hat{R}$, addresses whether chains appear to have reached a compatible distribution; it does not by itself show that a reported mean or tail quantile has small simulation error. Conversely, a large pooled ESS is not reassuring if chains disagree or have not explored the target distribution. First check that the sampling behavior supports a stationary posterior summary, then assess ESS and MCSE for the particular quantities to report.

A useful reporting check is to ask whether additional simulation could plausibly change the conclusion or the stated precision. If an MCSE is small compared with the posterior standard deviation and with the number of decimal places reported, the estimate of a mean may be adequately precise for that purpose. If it is not, more effective draws are needed. For tail probabilities or quantiles, assess tail-relevant ESS and MCSE rather than relying only on a good mean ESS. In joint models, review important components separately: an association estimate, a random-effects scale, and a derived survival probability can have substantially different Monte Carlo precision.

Several alternatives are useful but answer different questions. Running longer chains can increase ESS when sampling is otherwise functioning well. Running additional chains can help assess between-chain agreement and may add information, but merely adding chains does not fix a chain that fails to explore the posterior. Thinning by retaining every $k$th draw can reduce storage, but usually discards information and does not create more effective draws per unit of computation. Reparameterization or sampler changes may improve efficiency, but choosing among them requires diagnosing the sampling problem rather than treating ESS as a substitute for model or sampler assessment.

The main limitations are that ESS is estimated, depends on the summary being considered, and relies on the chains representing the target posterior. MCSE formulas use large-sample approximations; quantile MCSE is especially sensitive to how the density and tail dependence are estimated. A numerical MCSE cannot account for model misspecification, inappropriate priors, or uncertainty in the data-generating assumptions. It quantifies simulation error conditional on the fitted model and the sampling procedure.

Common pitfall: “There are 4,000 stored draws, so the posterior mean has the precision of 4,000 independent observations.” This ignores dependence. Use an ESS estimate for the quantity being summarized, then calculate MCSE from that ESS. Also do not confuse MCSE with posterior standard deviation: the former is uncertainty in the computed summary, while the latter describes posterior spread in the parameter.

A Practical Summary Workflow

Choose the reported quantity.
Identify whether it is a model parameter, a transformation such as a hazard ratio, a posterior mean, or a posterior quantile.

Check the chains.
Use trace plots and between-chain diagnostics to assess whether posterior summaries are interpretable before treating pooled ESS as evidence of precision.

Obtain the relevant ESS.
Use the quantity-specific ESS, and use tail- or quantile-focused diagnostics for quantiles and tail probabilities.

Calculate and report MCSE.
For a mean, use posterior standard deviation divided by the square root of ESS. Compare the resulting MCSE with the precision required for the scientific conclusion.

Check Your Understanding

Try each question, then compare your work with the solutions.

A scalar summary has $S=6{,}000$ retained draws and estimated integrated autocorrelation time $\hat{\tau}=3$. What is its ESS?
A joint-model coefficient has posterior standard deviation $0.24$ and ESS $900$. Calculate the approximate MCSE of its posterior mean.
Why should a hazard-ratio mean generally be calculated by transforming each posterior draw, rather than by exponentiating the mean log-hazard coefficient?
What ESS is most relevant to estimating a posterior 95th percentile: the ESS for the parameter mean, or a quantile-specific ESS? Explain briefly.
State one distinction between posterior standard deviation and MCSE.

Solutions. (1) $\operatorname{ESS}=S/\hat{\tau}=6{,}000/3=2{,}000$. (2) $\widehat{\operatorname{MCSE}}=0.24/\sqrt{900}=0.24/30=0.00800$, in the coefficient's units. (3) Exponentiation is nonlinear, so generally $\overline{\exp(\alpha)}\ne\exp(\bar{\alpha})$; transforming every draw estimates the posterior distribution of the hazard ratio directly. (4) A quantile-specific ESS is more relevant because quantile precision depends on the dependence of tail-indicator draws, which can differ from dependence relevant to the mean. (5) Posterior standard deviation describes posterior spread in the parameter; MCSE describes simulation uncertainty in a computed summary, such as the posterior mean.