Worked Joint Model for a Bounded Biomarker

Why Use a Beta Longitudinal Submodel?

Some biomarkers are proportions or normalized measurements that vary continuously between 0 and 1. Examples might include a hypothetical assay fraction or a normalized tissue measurement. Such values are bounded, and their variability may depend on their mean: observations near 0 or 1 cannot vary in the same way as observations near 0.5. A Gaussian longitudinal model can predict outside the allowable range and may represent this mean-dependent variability poorly.

Beta regression provides a model for measurements strictly inside the unit interval. In a joint model, its subject-specific mean trajectory can also enter an event hazard. This gives a coherent way to describe how the latent longitudinal process and event time are associated, while keeping the biomarker model on its natural scale. The association is statistical, not by itself causal.

Model and Assumptions

Definition: Let $Y_{ij}$ be subject $i$'s biomarker measurement at time $t_{ij}$, with $0<Y_{ij}<1$. A beta-regression longitudinal submodel is $$ Y_{ij}\mid b_i \sim \operatorname{Beta}\{\mu_i(t_{ij})\phi,\,[1-\mu_i(t_{ij})]\phi\}, \qquad \operatorname{logit}\{\mu_i(t)\}=\mathbf{x}_i(t)^\mathsf{T}\boldsymbol{\beta} +\mathbf{z}_i(t)^\mathsf{T}b_i. $$ Here $\mu_i(t)=E\{Y_i(t)\mid b_i\}$, $\phi>0$ is the precision parameter, and $b_i$ is a subject-specific random effect, commonly modeled as $b_i\sim N(0,D)$. The corresponding conditional variance is $\mu_i(t)[1-\mu_i(t)]/(\phi+1)$.

The mean is constrained to the open unit interval by the inverse-logit link. Precision controls how tightly observations are distributed around that mean: larger $\phi$ gives smaller variance at a fixed mean. In the version used here, $\phi$ is constant across subjects and times. Extensions can allow precision to vary, but that changes the model and should be supported by the data.

For the event component, let $T_i$ be the event time, $\delta_i$ the event indicator, and $\mathbf{w}_i$ baseline covariates. One common current-value association is

h_i(t\mid b_i)=h_0(t)\exp\{\boldsymbol{\gamma}^\mathsf{T}\mathbf{w}_i+\alpha\mu_i(t)\}, \qquad S_i(t\mid b_i)=\exp\left\{-\int_0^t h_i(u\mid b_i)\,du\right\}.

The parameter $\alpha$ links the fitted conditional biomarker mean to the event hazard. For example, a positive $\alpha$ means that, at a given time and covariate values, a larger $\mu_i(t)$ corresponds to a larger conditional hazard. The hazard is an instantaneous rate among subjects still event-free, not the probability of an event at that instant.

This specification rests on several substantive conditions. The beta response distribution requires measurements strictly between 0 and 1; conditional on the random effects and covariates, repeated biomarker values are taken to be independent beta observations. The random-effects distribution and both submodels must be adequate descriptions of the data. Given the shared random effects and included covariates, the event and longitudinal processes are assumed conditionally independent; their modeled association is represented through the common latent effects and the chosen trajectory link. Censoring must be noninformative conditional on the modeled information, and the observation process must be ignorable or modeled separately. These assumptions are not guaranteed merely because the model converges.

How the Two Data Streams Form One Model

For a subject observed until event or censoring time $T_i$, the event contribution conditional on $b_i$ is $h_i(T_i\mid b_i)^{\delta_i}S_i(T_i\mid b_i)$. Multiplying it by the conditional beta densities gives the subject's conditional joint contribution. Integrating over the random-effects distribution accounts for the fact that $b_i$ is unobserved:

L_i=\int \left\{\prod_j f_{\operatorname{Beta}}\bigl(y_{ij}\mid\mu_i(t_{ij}),\phi\bigr)\right\} h_i(T_i\mid b_i)^{\delta_i}S_i(T_i\mid b_i)\,p(b_i)\,db_i.

This expression explains why separately fitting a beta regression and a survival model, then treating their estimates as unrelated, is not generally equivalent to fitting the joint model. The event contribution informs the distribution of plausible subject-specific effects, and those effects also determine the subject's biomarker trajectory. The joint likelihood uses both sources of information while accounting for that shared latent structure.

Worked Example: Specifying and Reading a Mean Trajectory

Suppose, hypothetically, that a normalized biomarker is recorded over time in years. Use a random-intercept beta model with $\operatorname{logit}\{\mu_i(t)\}=-0.40+0.25t+b_i$ and precision $\phi=10$. For one hypothetical subject, take $b_i=0.20$. Find the subject's fitted mean at years 0 and 2.

At year 0, the linear predictor is $-0.40+0.25(0)+0.20=-0.20$. Therefore, $\mu_i(0)=1/(1+\exp(0.20))\approx0.4502$.

At year 2, the predictor is $-0.40+0.25(2)+0.20=0.30$, so $\mu_i(2)=1/(1+\exp(-0.30))\approx0.5744$.

The fitted conditional mean increases by about $0.1242$ of the normalized biomarker scale over these two years. These are means conditional on this subject's random intercept, not raw measurements and not population-average means after integrating over the random-effects distribution. At year 2, the model-implied conditional variance is $0.5744(1-0.5744)/(10+1)\approx0.02222$.

What a Beta Density Contribution Means

If $Y\mid b_i$ has mean $\mu$ and precision $\phi$, its beta shape parameters are $a=\mu\phi$ and $c=(1-\mu)\phi$. Its density at an interior value $y$ is

f_{\operatorname{Beta}}(y\mid\mu,\phi) =\frac{y^{\mu\phi-1}(1-y)^{(1-\mu)\phi-1}} {B\{\mu\phi,(1-\mu)\phi\}}, \qquad 0<y<1.

A density value is not the probability of observing exactly $y$. For a continuous response, probabilities are calculated over intervals by integrating the density. This distinction matters when interpreting the longitudinal factor in a joint likelihood.

Worked Example: Calculating a Beta Measurement Density

Suppose, hypothetically, a measurement has conditional mean $\mu=0.60$, precision $\phi=10$, and observed value $y=0.50$. The shape parameters are $a=0.60(10)=6$ and $c=0.40(10)=4$. Thus $B(6,4)=\Gamma(6)\Gamma(4)/\Gamma(10)=120(6)/362880=1/504$.

The density is $f(0.50)=0.50^{5}0.50^{3}/B(6,4)=0.00390625(504)=1.96875$, or about 1.969 per unit of the biomarker scale. A density greater than 1 is valid; it is not a probability. For example, the probability of a measurement in an interval would require integrating this beta density over that interval, not reading the density at one point as that probability.

Relating the Fitted Trajectory to Event Risk

Under the current-value hazard model, $\exp(\alpha\Delta)$ is the conditional hazard ratio for two subjects whose fitted means at the same time differ by $\Delta$, with baseline covariates and other modeled terms held fixed. This is a comparison of instantaneous hazards. It is not generally a risk ratio, because cumulative risk depends on the hazard over an interval and therefore on the evolving trajectory and baseline hazard.

Worked Example: Conditional Hazard and Survival Along a Trajectory

Continue with the hypothetical subject from the first example. Suppose there are no additional covariates, the baseline hazard is constant at $h_0(t)=0.08$ per year, and $\alpha=1.20$. The conditional hazard is $0.08\exp\{1.20\mu_i(t)\}$ per year. At years 0, 1, and 2, the fitted means are approximately 0.4502, 0.5125, and 0.5744, respectively. The corresponding hazards are approximately $0.08\exp(1.20(0.4502))=0.1373$, $0.08\exp(1.20(0.5125))=0.1480$, and $0.08\exp(1.20(0.5744))=0.1594$ per year.

To illustrate conditional survival through year 2, approximate the cumulative hazard by Simpson's rule using these three values: $\int_0^2 h_i(u\mid b_i)\,du\approx(2/6)[0.1373+4(0.1480)+0.1594]=0.2962$. Thus $S_i(2\mid b_i)\approx\exp(-0.2962)=0.7437$, giving conditional event risk by year 2 of about $1-0.7437=0.2563$. This is a numerical approximation conditional on the specified random effect and model, not a population-average risk.

At year 2, comparing this subject with a hypothetical subject whose fitted mean is 0.10 lower, the conditional hazard ratio is $\exp\{1.20(0.10)\}=\exp(0.12)\approx1.1275$. That is about a 12.75% higher instantaneous hazard for the subject with the larger mean, at that time and under the stated conditions. It does not say that the subject's two-year event risk is 12.75% higher.

Diagnostics, Alternatives, and Limitations

Check whether the data and fitted model support the assumptions before interpreting the association. The beta model is not directly defined for exact zeros or ones. First determine whether such values are genuine boundary outcomes, rounding artifacts, or values produced by the measurement process. Do not silently replace them with arbitrary nearby values; a model for boundary mass or a different response formulation may be needed. For interior values, compare observed and fitted distributions across time and relevant covariate groups, paying attention to mean-dependent spread and skewness. Residual checks should be appropriate to the fitted beta distribution rather than assuming Gaussian residuals.

Also inspect the fitted trajectory for implausible behavior, assess whether the assumed precision is adequate, and check the event model's fit, including its baseline hazard and the form of the biomarker association. Convergence problems, sensitivity to starting values, unstable random-effects estimates, or inadequate numerical integration can make estimates unreliable. A fitted mean trajectory can look reasonable while the event component is misspecified, so both parts require assessment.

Approach	Useful when	Main caution
Beta regression in a joint model	The biomarker is continuous and strictly between 0 and 1, with potentially mean-dependent variance.	It does not directly accommodate exact boundary values; distribution and link assumptions need checking.
Gaussian longitudinal submodel	The response is approximately unbounded on its analysis scale and conditional residual behavior is reasonable.	For a bounded response it may predict outside the support or misrepresent skewness and variance.
Fractional-response or other bounded-outcome model	Values at 0 or 1 occur and the chosen method can accommodate them.	Its likelihood or variance assumptions differ from beta regression, so it is not an interchangeable fit.

A further limitation is interpretive: the fitted association parameter describes the specified conditional model. It does not establish that changing the biomarker would change event risk. The result can also depend on whether the event hazard is linked to the current mean, a slope, or another feature of the latent trajectory. These alternatives answer different scientific questions and should not be conflated.

Common Mistake

Common mistake: Interpreting $\exp(\alpha\Delta)$ as a ratio of event probabilities, or treating a beta density at an observed value as the probability of that value. The correction is to identify the quantity first: $\exp(\alpha\Delta)$ is a conditional hazard ratio under the stated model, while a beta density is a density value. Event probabilities require survival calculations over an interval, and measurement probabilities require integration over a set of response values.

The practical value of this joint specification is that it respects the biomarker's bounded support while connecting its subject-specific mean trajectory to event timing. Its credibility still depends on support, distributional, dependence, censoring, and hazard assumptions; those assumptions determine what the fitted trajectory and event association can reasonably mean.

Check Your Understanding

Answer each question, then compare with the solutions below.

For a beta model with $\mu=0.30$ and $\phi=20$, what are the two beta shape parameters?
If $\alpha=0.8$, what is the conditional hazard ratio for a 0.05 difference in fitted means?
Why is a density of 1.969 at one biomarker value not a probability of 1.969?
Name one reason an exact zero biomarker value cannot be handled by the stated beta density.

Solutions. (1) The shape parameters are $0.30(20)=6$ and $0.70(20)=14$. (2) The ratio is $\exp\{0.8(0.05)\}=\exp(0.04)\approx1.041$, a conditional hazard ratio, not a risk ratio. (3) A density is probability per unit of response scale; probabilities come from integrating density over an interval. (4) The beta density specified here has support only for values strictly between 0 and 1, so a model allowing boundary mass or another suitable approach is required.