Authors: Qing Xi Ooi (Pharmetheus AB, Uppsala, Sweden), Elodie Plan (Pharmetheus AB, Uppsala, Sweden), Martin Bergstrand (Pharmetheus AB, Uppsala, Sweden)
Categories: Tutorial
Source: CPT: Pharmacometrics & Systems Pharmacology
Doi: 10.1002/psp4.13278
A tutorial on pharmacometric Markov models
Authors: Qing Xi Ooi, Elodie Plan, Martin Bergstrand
Abstract
The Markov chain is a stochastic process in which the future value of a variable is conditionally independent of the past, given its present value. Data with Markovian features are characterized frequent observations relative to the expected changes in values, many consecutive same‐category or similar‐value observations at the individual level, and a positive correlation observed between the current and previous values for that variable. In drug development and clinical settings, the data available commonly present Markovian features and are increasingly often modeled using Markov elements or dedicated Markov models. This tutorial presents the main characteristics, evaluations, and applications of various Markov modeling approaches including the discrete‐time Markov models (DTMM), continuous‐time Markov models (CTMM), hidden Markov models, and item‐response theory model with Markov sub‐models. The tutorial has a specific emphasis on the use of DTMM and CTMM for modeling ordered‐categorical data with Markovian features. Although the main body of this tutorial is written in a software‐neutral manner, annotated NONMEM code for all key Markov models is included in the Supplementary Information.
OVERVIEW
The Markov process is named after the mathematician Andrey Andreyevich Markov (1856–1922), who was the first that studied this systematically. The Markov process is a stochastic process in which the future value of a variable is conditionally independent of the past, given its present value. If the random variable, Yij, for the ith subject at the jth occasion takes on discrete values, kij, the Markov process is characterized by(1)PYij=kij|Yi0=ki0,Yi1=ki1,…,Yij−1=kij−1=PYij=kij|Yij−1=kij−1and the sequence Yi0,Yi1,Yi2,…,Yij with the Markov property is called a Markov chain. In the Markov process, only knowledge of the current state is needed to predict the future state, which is independent of all previous states. In the example above, the Markov chain is of first order since the predicted value depends only on the value immediately preceding it. Conversely, higher‐order Markov chains are processes in which earlier occurring values are also considered. Markov processes find application in various real‐world settings, ranging from evolutionary biology, disease spread, population dynamics, as well as weather predictions, cruise control systems, mobility patterns, queue systems and even consumer behavior patterns, voting patterns, stock market prices, and currency exchange rates. It is an important part in many machine learning or artificial intelligence applications, natural language processing, and search engine algorithms.
In pharmacometric analyses, the data available commonly present Markov features and are increasingly often modeled using Markov elements or dedicated Markov models (based on a search on PubMed using the keywords “markov pharmacometrics”). In many cases, Markov elements can be added to different types of models, such as continuous models, count models, and ordered‐categorical models; in these cases, considering first‐order Markov elements is often sufficient. Unlike Markov elements, dedicated Markov models have transitions from one state to another, for example, transitions from no response to partial response or from full response to partial response, modeled as transition probabilities or as transition rate constants.
This tutorial aims to introduce the foundations of population analysis of data with Markovian features to pharmacometricians who are already familiar with basic concepts of population analysis. The tutorial describes first the Markov elements, then various Markov models, and concludes with a discussion on key considerations in Markov modeling as well as recommendations on the Markov model of choice to analyze data with Markovian features. The main focus of this tutorial will be on the analysis of ordered‐categorical data with Markovian features, with a specific emphasis on the discrete‐time Markov model (DTMM) and continuous‐time Markov model (CTMM). Markov modeling can be performed using the main nonlinear mixed‐effects software. Although the main body of this tutorial is written in a software‐neutral manner, annotated NONMEM code for all key Markov models is included in Supplementary Information [Link], [Link], [Link], [Link], [Link].
MARKOV ELEMENTS
Markov elements in continuous models
Compared to discrete models, pharmacometrics models for a continuous variable make less assumptions of independence between subsequent model predictions. However, the default assumption is that the random effects attributable to the within‐subject variability, which includes both the residual unexplained variability (RUV) and inter‐occasion variability (IOV), are uncorrelated between subsequent observations. Relaxation of this assumption of independence by characterizing the Markovian nature of these random effects for consecutive timepoints is exemplified below, using RUV.
The RUV is often assumed to be independent and identically distributed by default. Non‐independence between successive RUVs is commonly observed when sampling is frequent, such as in electroencephalogram evaluations, or when the structural and/or inter‐individual (IIV) model is misspecified, for example, due to mechanistic deficiency of the model or response‐adaptive designs. The correlation between consecutive timepoints of a continuous variable is commonly referred to as serial correlation or autocorrelation. The application of Markov modeling in clinical studies with response‐adaptive design is further explored in the discussion section.
Serial correlation in RUVs can be visualized by plotting current versus previous residual values, either as population or individual values (e.g., individual conditional weighted residuals), if time intervals are equidistant; else, if time intervals are not equidistant, one can plot residual values versus time, by subject. Based on these plots, the Markovian features are considered prominent if a systematic, and often positive, trend is observed in the current versus previous residual values plot or when subject‐level residuals of the same sign are found clustering at timepoints in close proximity in the longitudinal residual plots.
One of the consequences of ignoring the serial correlations in RUVs is the inflation of Type I error when conducting the likelihood ratio test, potentially leading to false identification of covariate effects, such as the placebo effect. ^1^ Other consequences are the possible bias in the IIV and/or RUV estimates and less realistic simulations.
Serial correlations in RUVs can be modeled using the first‐order autoregressive (AR(1)) model. ^2^ , ^3^ In the AR(1) model, the correlation between the RUV of two consecutive timepoints (tj−1 and tj) is assumed to decrease exponentially with the increased time elapsed between the observations, as (2)corrRUVtj−1,RUVtj=e−tj−1−tjtcorrwhere tcorr is a parameter determining the rate at which the correlation decreases with time. Based on Equation 2, the correlation between RUVs decreases over time at a rate determined by tcorr and for any other given tj−1 and tj values where tj−1≠tj, a high tcorr corresponds to a strong correlation and vice versa. For better interpretability, tcorr may be reparameterized as a half‐life parameter. Alternative models that have been suggested for describing such autocorrelated within‐subject variability include the dynamic IOV model ^4^ and the stochastic differential equations (SDEs). ^5^ Examples of the use of continuous models with Markov elements include modeling of Crohn's Disease Activity Index (CDAI) using AR(1) model ^6^ and temperature variations using SDEs. ^7^
Markov elements in count models
Count data can be defined as numbers of units, and can take integer values of 0 or higher, theoretically until infinity. In drug development, count data are typically number of events per time unit, for example, number of migraines per week, or number of gastrointestinal symptoms per day, or other types of measurement such as number of words remembered in Alzheimer's Disease Assessment Scale–Cognitive Subscale (ADAS‐Cog) tool, or antidrug antibody titers in the immunogenicity assay. During the pharmacodynamics modeling step, the preferred models that are developed to describe this type of data and make inference are the count models, such as the Poisson model or other models of this family. While intuitively the count obtained during a time unit j would be dependent on the count observed within the time window j−1, count models carry an assumption of independence between events, that is, events are memoryless. However, this type of assumption may lead to a poor description of the data or poor simulation properties. To avoid this, count models can be augmented with the introduction of Markov elements (Equation 3):(3)PYij=n|Yij−1=fn;Mi;Yij−1where Mi is the expectation for the Poisson distribution for the ith individual and n is the count value. Markov elements in count models can include the exact preceding count, treated for example with a maximum effect (Emax) function; alternatively, they can include a derived value, such as a dichotomized value (e.g., presence or absence of an event), which serves as an input into, for example, a mixture function. This has been described in a previously published tutorial. ^8^
Examples of the use of count models with Markov elements include modeling of epilepsy seizure count, ^9^ , ^10^ , ^11^ daily pain scores in diabetic patients with neuropathic pain, ^12^ and multiple sclerosis contrast‐enhancing lesions. ^13^
Markov elements in ordered‐categorical models
An ordered‐categorical variable is a variable with a set order of categories and for which the distances between the categories are unknown; this type of variable may also be created from the categorization of a continuous variable. A binary variable is the simplest form of the ordered‐categorical variable. Examples of an ordered‐categorical variable include non‐responder/responder, no/partial/full responses, and no/mild/moderate/severe adverse events (AEs). For the remainder of this tutorial, the following ordered‐categorical variable, or minor variations of it, will be used as an example of response for 50 subjects receiving placebo or 100 mg of a hypothetical drug; three ordered categories will be (i) no response, (ii) partial response, and (iii) full response, which are denoted by a value of 1, 2, and 3, respectively.
The presence of Markovian features in ordered‐categorical data can often be gauged based on both the mechanistic understanding of the variable of interest and the information on relative sampling frequency. To explore the Markovian feature of ordered‐categorical data, one can consider the visualizations shown in Figure 1. Data with Markovian features are characterized by frequent observations relative to the expected changes in values (Figure 1, Panel b), many consecutive same‐category observations at the individual level (Figure 1, Panel b), and a positive correlation observed for current versus previous ordered‐categorical variable values (Figure 1, Panel d). Figure 1 was generated using observations from a typical subject; however, similar plots as those in panels c and d may be generated using population data to evaluate the Markovian features of such type of data.
Various models have been proposed to model ordered‐categorical data, including proportional odds (PO) models ^13^ , ^14^ , ^15^ , ^16^ and bounded‐integer models. ^17^ , ^18^ , ^19^ Markov elements may be added to these models to account for the Markovian features in the data. A PO model considers the cumulative logits, which are the logarithmic odds of falling into a category when being less than or equal to a certain threshold value k, that is, LogitPYij≤k. All cumulative logits are modeled except LogitPYij≤K, where K is the categorical value of the highest order, since LogitPYij≤K=1. The use of the cumulative logits accounts for the ordered nature of the variable and ensures that the probabilities of all possible values for the variable sum up to one. In the current tutorial, the logit transformation is used as an example to map probability values (ranging from zero to one) into real numbers (ranging from negative infinity to positive infinity), prior to regression of the IIV or covariate effects. Other link functions, such as the probit transformation, may also be considered. In addition, while the description in this tutorial uses LogitPYij≤k, it is worth noting that the model may also be structured to consider the LogitPYij≥k instead.
In the most simple form, the cumulative logits are functions of their respective intercept terms, that is, LogitPYij≤k=γk, where γk=logPYij≤k1−PYij≤k. However, to ensure that PYij≤k≥PYij≤k−1, γk is typically reparameterized as γk=α1+a2+…+αK−1, with α1 being estimated without any constraint and successive α2,…,αK−1 constrained to be zero or positive. The cumulative logits may also include a vector of additive functions, fθx, of a vector x of p predictors, x1ij,…,xpij, such as drug exposure and demographics, for the ith individual at the jth occasion and with the predictor‐logarithmic odds relationships governed by a vector of slope parameters θ. Noteworthy, the PO assumption is made in the PO model whereby the effect of a predictor is the same on the logarithmic odds scale for different levels of the response.(4)LogitPYij≤k=γk+fθ,x
Since the cumulative logit ranges from negative infinity to positive infinity, the corresponding cumulative probability that ranges from zero to one can be obtained using an expit (inverse‐logit) transformation of the relevant cumulative (5)PYij≤k=11+e−LogitPYij≤k
The probability for each category that ranges from zero to one may in turn be derived based on the cumulative probabilities as (6)PYij=k=PYij≤1ifk=1PYij≤k−PYij≤k−1if1<k<K1−PYij≤K−1ifk=K
The Markov element may be introduced by assuming that one or more of the threshold‐specific intercept parameters (αk), or slope parameters that define the covariate‐parameter relationships (θk), may vary according to the previous dependent variable value(s), that is, αk∣kij−1 and/or θk∣kij−1, where k is the threshold value and kij−1 is the previous observation. Using a variable with three ordered categories as an (7)mψMarkov,kij−1=Ikij−1=1·ψMarkov1 α1∣kij−1=α1+mψMarkov,kij−1 LogitPYij≤1=α1∣kij−1+fθ,x LogitPYij≤2=α1∣kij−1+α2+fθ,xwhere Ikij−1=1 is an indicator variable for the previous response being equal to 1 and ψMarkov1 is a horizontal shift parameter on the logit scale if Ikij−1=1=1 and is estimated without constraint. If ψMarkov is positive, the observation of the previous response being equal to 1 increases LogitPYij≤1 and LogitPYij≤2. This in turn increases PYij=1 and decreases both PYij=2 and PYij=3, thereby enforcing the Markovian nature of the observed data that kij−1=1 increases the chance of kij=1. In addition, mψMarkov,kij−1 represents a generic function to account for the Markov elements. More complex Markov elements may also be considered, such as the Markovian influence of more levels of previous response (e.g., mψMarkov,kij−1=Ikij−1=1·ψMarkov1+Ikij−1=2·ψMarkov2 +Ikij−1=3·ψMarkov3); the inclusion of Markov elements on other intercept or slope parameters; accounting for the attenuation of Markov elements with lapsed time if time intervals of observations are not equidistant.
Examples of the use of ordered‐categorical models with Markov elements include modeling of hand‐and‐foot syndrome severity, ^13^ diarrhea/acne occurrence, ^14^ dizziness, ^15^ and sedation score, ^16^ all using the logistic or PO model, as well as modeling of the Likert neuropathic pain scale ^17^ and symptom scales specific to lower urinary tract symptoms, ^19^ using the bounded‐integer model.
One key disadvantage associated with PO models is the PO assumption that all ordered categories are affected by the predictors to the same extent on the logarithmic odds scale. In case of assumption violation, the differential odds model, ^20^ DTMM, or CTMM may be considered.
MARKOV MODELS
There are different types of Markov models, as summarized in Figure 2. In Markov models, the transitions from one state to another state are modeled.
Discrete‐time Markov model (DTMM)
Theory and basic model structure
In a DTMM, the dependent variable transitions from one state to another with a certain probability at each distinct, equidistant point in time (i.e., discrete‐time step); this probability is dependent only on the immediately preceding state and not on all other previous states. The corresponding mathematical representation is shown in Equation 1 with j representing discrete timepoints.
In a DTMM, the probabilities of all possible pairwise transitions between states can be estimated, fixed, or derived. Given that transition probabilities are modeled, the probability of each state is necessarily conditioned on the previous state and governed by the transition probabilities estimated, fixed, or derived. In a DTMM, the impact of time since the last observation on the transition probabilities (and the resulting probability of each state) is assumed to be negligible.
An absorbing state (A), a state that is impossible to leave once it is reached, may be included in a DTMM. Once the absorbing state is reached, then PA|A=1 and other transition probabilities from this state to other states are necessarily zero. The absorbing state is often included in a DTMM (and CTMM) to account for study dropout. ^21^ , ^22^ Other applications of the absorbing state include modeling gastrointestinal transit ^23^ and disease progression. ^22^ Using dropout (D) as an example, the transition probabilities to the dropout state, that is, PD|1, PD|2, and PD|3 in Panel b of Figure 2, are usually estimated unless the probability to dropout is negligible or deterministic, the latter being most often enforced by study protocol based on certain criteria (e.g., when a patient with progressive disease must be taken off the treatment).
For a variable with two categories (1 and 2), a total of four transition probabilities (P(1|1), P(2|1), P(1|2), P(2|2)) are relevant; only two parameters (one of P(1|1) and P(2|1) plus one of P(1|2) and P(2|2)) need to be estimated or fixed, and the remaining two are derived given that P(1|1) + P(2|1) = 1 and P(2|2) + P(1|2) = 1. These transition probabilities may be estimated directly with the parameter space constrained to be between zero and one. However, to explore the IIV or covariate effect, a logit transformation can be applied to each typical probability for the transition from the immediately preceding state, kj−1, to the current state, kj; in addition, covariate effect(s) fθx (or IIV) may be included as an additive on the logit scale. Further details on the inclusion of covariate effects are provided in the covariate analysis section below.(8)LogitPYij=kj|Yij−1=kj−1=logPtypYij=kj|Yij−1=kj−11−PtypYij=kj|Yij−1=kj−1+fθ,xwhere Ptyp is the transition probability for a subject with typical covariate values.
In the next step, the resulting term is back‐transformed using an expit function, generating a transition probability that is constrained between zero and one.(9)PYij=kj|Yij−1=kj−1=11+e−LogitPYij=kj|Yij−1=kj−1
For a variable with three or more categories, the cumulative logits are considered in the same manner as in the PO model. Using a variable with three ordered categories (named 1, 2, and 3) as an example, a total of nine transition probabilities (P(1|1), P(2|1), P(3|1), P(1|2), P(2|2), P(3|2), P(1|3), P(2|3), P(3|3)) are relevant, six intercept terms will need to be estimated or fixed, and the remaining probabilities will be derived.(10)LogitPYij≤kj|Yij−1=kj−1=γkj∣kj−1+fθ,x|kj−1 PYij≤kj|Yij−1=kj−1=11+e−LogitPYij≤kj|Yij−1=kj−1 PYij=kj|Yij−1=kj−1=PYij≤kj|Yij−1=kj−1ifkj=1PYij≤kj|Yij−1=kj−1−PYij≤kj−1|Yij−1=kj−1if1<kj<Kj1−PYij≤kj−1|Yij−1=kj−1ifkj=Kj
Here, γkj∣kj−1 is the sum of the relevant intercept terms, that is, αkj∣kj−1, with the first term estimated without constraint and successive terms constrained to be zero or positive. In addition, fθ|x,kj−1 is a function describing the effect of various predictors on the relevant cumulative logit term.
A DTMM differs from a PO model in that the cumulative logit for the transition probabilities, which are conditional probabilities, are considered rather than the cumulative logit for the probability of being in specific categories, which are marginal probabilities. Accordingly, using a variable with three ordered categories as an example, six cumulative logits for P(1|1), P(2|1), P(1|2), P(2|2), P(1|3), and P(2|3) are modeled in a DTMM in contrast with the two cumulative logits for P(1) and P(2) modeled in a PO model. While the probability modeled may also be conditional on the previous observation in a PO model, with the inclusion of the Markov element, the probabilities modeled are necessarily conditional for DTMM. Depending on how the models are set up and how the covariate effects are considered, the PO model with Markov elements may coincide with the DTMM, for instance when all possible Markov elements are added on all the intercept terms and the PO assumption holds. ^24^ , ^25^ Furthermore, the covariate effects are commonly implemented as conditional on the previous state in a DTMM, while in a PO model the same covariate effects apply regardless of the previous state, unless otherwise implemented using the Markov elements. Another key difference between a PO model with Markov elements and a DTMM is that in the latter the relaxation of the PO assumption is often applied by default. In PO models, the assumption that the effect of a predictor is constant for different levels of response is usually assumed unless otherwise explicitly relaxed in the differential odds models. ^20^ In other words, in a DTMM, the function describing the predictor‐logarithmic odds relationship usually differs based on the threshold values kj, that is, fθ,x|kj,kj−1. For instance, the effect of drug concentrations on treatment response may be different on LogitPYij≤2|Yij−1=kj−1 compared to LogitPYij≤1|Yij−1=kj−1, indicating variating drug effects on the different threshold values, for example, no‐response category vs. partial response category.
For DTMM, the initial estimate for the transition probabilities may be set according to the observed transition probabilities in the absence of active treatment. The intensity of the Markov property for the resulting DTMM is governed by the relative magnitude of the estimate for the transition probabilities. For instance, a P(1|1) estimate larger than the P(2|1) and P(3|1) estimates indicates a strong Markov property with a greater tendency for consecutive same‐state observations of 1 rather than cross‐state observations given a previous state of 1.
Covariate analysis
Covariate effects are usually included as an additive on the logit scale for the transition probability parameters.(11)LogitPcov=LogitθP+fθcov,xwhere θP is the typical value of the parameter P for a subject i with covariate value(s) equivalent to that used for centering, Pcov is the parameter value for a subject given the individual covariate value(s), θcov is the vector of parameter value(s) defining the covariate‐parameter relationship(s), x is the vector of individual covariate value(s) for a certain covariate, and fθcov,x is the function defining the covariate‐parameter relationship.
A linear model (Equation 12) is commonly used and fit‐for‐purpose; in addition, other models, including an Emax model (Equation 13) for drug effects, may be considered on the logit scale. Conversely, a power model or an exponential model is usually not appropriate since the resulting covariate effects are necessarily greater than zero, thereby precluding the possibility of a covariate effect being inhibitory on the base probability.(12)fθcov,x=θslope·x−xcent (13)fθcov,x=θEmax·xθEC50+xwhere θslope is the slope parameter, θEmax is the maximum covariate effect, θEC50 is the covariate value that gives half the maximal covariate effect, and xcent is the covariate value used for centering (typically the median covariate value).
For categorical covariates, a covariate model with the estimation of a shift parameter θshift is often used (Equation 14).(14)fθcov,x=0ifx=xcent0+θshiftifx≠xcent
IIV model
The IIV is usually included as an additive on the logit scale for the probability parameters and ranges from negative infinity to positive infinity. This ensures that after the back‐transformation of the logit, using the expit function, the resulting individual probability values are necessarily constrained between zero and one, as intended.(15)LogitPi=LogitPcov+ηi Pi=11+e−LogitPiwhere Pi is the individual value of the parameter and ηi is a normally distributed random variable with a mean of zero and standard deviation ωP. The variance–covariance matrix of the ηi is denoted as Ω.
For parameters with a lower bound of zero, for example, θslope and θEC50, the exponential IIV model is normally (16)Pi=Pcov·eηi
Covariance between ηi can be considered if supported by the data.
Starting the Markov chain and data requirement
In the analysis, the unobserved preceding ordered‐categorical value before the first observation is usually assumed to be the same as the first observed value.
No additional data formatting is required for the specific population analysis using DTMM in NONMEM. A trimmed and annotated example of a NONMEM data set is shown in Supplementary Information S2 and the full data set is provided in Supplementary Information S7.
A discussion of NONMEM estimation methods, with estimation setting examples, is provided in Supplementary Information S1.
Model evaluation
Various visual predictive checks (VPCs) may be generated to evaluate the fit of a DTMM. The recommended plots for the evaluation of a Markov model are shown in Figure 3. Most often, it is desirable to evaluate the ability of the Markov model to characterize the observed longitudinal proportion of subjects in each ordered category (Figure 3, Panel a) as well as in each unique transition category (Figure 3, Panel b); this can be shown with or without stratification on key covariates of interest, for example, treatment or dose allocation. VPCs may also be generated for selected metric(s) of interest (Figure 3, Panel c) such as mean number of days with full response or AE. If the time to the first event is of interest, a VPC for Kaplan–Meier plots may be considered. As for the evaluation of all pharmacometric models, the model evaluation strategy should be adapted keeping in mind the goal for applying the Markov model.
Applications
Examples of the use of a DTMM include modeling of sleep stages, ^26^ improvement in rheumatoid arthritis, ^21^ muscle spasm grade, ^27^ fatigue level and hand and foot syndrome severity, ^25^ and central‐nervous‐system‐specific side effects. ^24^
Overall, a DTMM is straightforward and easy to communicate (compared to a CTMM). However, a basic DTMM assumes that the effect of time since the last observation on transition probabilities is negligible. When the time intervals between observations are equidistant (e.g., daily observations) this does not really matter but it should be noted that under these circumstances the model should only be used to simulate the same type of equidistant observations. DTMM can be extended to account for the time effects. However, such an approach is often not parsimonious or identifiable. In addition, explicit parameters are needed for all observed transitions between states and the application of DTMM is also limited to the prediction of transitions with observations only. For systems with several possible states, many parameters need to be estimated in the DTMM; in addition, the description of covariate effects is complex and not parsimonious, and multiple possible ways to include covariate effects need to be considered, for instance, on which transition probabilities to include covariate effects. A CTMM should be considered in the following if the time‐equidistant assumption is violated; when one wants to simulate data with non‐equidistant time intervals or with equidistant time intervals that differ from those observed; or when the variable to be modeled has many ordered categories.
Continuous‐time Markov model (CTMM)
Theory and basic model structure
The CTMM is based on a compartmental structure with one compartment allocated for each unique or lumped ordered‐categorical value(s). The probability of these ordered‐categorical values is modeled in CTMM. The amount in each compartment gives the probability of the corresponding ordered‐categorical value and the probability in all possible compartments sums up to one. The Markov property is underpinned by the continuous transfers of probabilities between compartments; the rate of probability distribution is governed by the first‐order transition rate constants between adjacent compartments (λ), which are constrained to be equal to or greater than zero. The probability of each state is hence necessarily conditioned on the previous state and the transition rate constant parameters in the CTMM. In a CTMM, the transition of the dependent variable from one state to another is not restricted to occur only at distinct, equidistant timepoints. Time is regarded as a continuous variable. Thus, the impact of time since the last observation on the probability of each state is considered in a CTMM, where the Markov property is the strongest at zero time since the last observation and the influence of the preceding observation decreases as the time between observations increases.
A CTMM may be implemented as algebraic equations, for instance, for a variable with two categories. However, the CTMM is most often described by a system of ordinary differential equations (ODEs).(17)dPYij=kdt=λk+1,k·PYij=k+1−λk,k+1·PYij=kifk=1λk−1,k·PYij=k−1+λk+1,k·PYij=k+1−λk,k−1+λk,k+1·PYij=kif1<k<Kλk−1,k·PYij=k−1−λk,k−1·PYij=kifk=K
With the default parameterization, all the transition rate constants are estimated or fixed. The CTMM may be reparameterized as equilibrium time (ET) and steady‐state probabilities (PSS) by assuming equilibrium of the system of ODEs (with no absorbing state), where the rate of probability input is equal to that of the output. PSS is the hypothetical probability of observing a certain ordered‐categorical value if the time since the last disequilibrium is infinite. PSS also assumes a closed system without trap states, such as the dropout state. In Figure 4, the probability of each ordered‐categorical value should stabilize if dropout is not possible. These stabilized probabilities should correspond to the PSS estimates and the time needed for the stabilization should correspond to the ET estimates. Hence, PSS estimates can be interpreted as steady‐state probabilities for the respective states, conditioned on no dropout. During parameter estimations, ET is constrained to be equal to or greater than zero; PSS, as with all probabilities, is constrained to be between zero and one. As described before, for parameter estimation, an expit transformation may be considered, and further reparameterization of PSS, for instance, as cumulative probabilities may be required for a variable with three or more ordered categories.(18)λk−1,k=1ETk−1,k·1+PSS,k−1PSS,k λk,k−1=λk−1,k·PSS,k−1PSS,k PSS,1=1−PSS,2+⋯+PSS,k+⋯+PSS,K
This reparameterization offers an alternative not only drug effects can be included on steady‐state probabilities, which is potentially more parsimonious than considering the drug effects on separate transition rate constants, but the model can also be simplified to a minimal CTMM (mCTMM), as detailed in later sub‐sections.
Similar to the case of a DTMM, absorbing state(s) may be included for the CTMM, for instance, to account for dropout, with the transition rate constant parameters estimated or fixed. Theoretically, with the presence of the absorbing state(s), the system will always be in disequilibrium. However, it is worth noting that the inclusion of the absorbing state(s) does not preclude using the ET and PSS parameterization, but these parameters should rather be interpreted as parameters at a steady state without the interference of the absorbing state(s).
For CTMM, the initial estimate for ET and PSS may be set according to the longitudinally observed proportion of subjects for each possible ordered‐categorical value (thick solid lines in panel a of Figure 3) and by assuming the absence of dropout. If the transition rate constant parameterization is used for CTMM, the initial estimates for these transition rate constants may be derived based on the initial estimates for ET and PSS using Equation 18.
For CTMM, the intensity of the Markov property is governed by the relative magnitude of the estimated/derived transition rate constants. For instance, a large λ11 value relative to λ12 and λ13 represents a strong Markov property with a greater tendency for consecutive same‐state observations of 1 rather than cross‐state observations given a previous state of 1.
Covariate analysis
In CTMM, covariate effects may be considered on λ, ET, and PSS, depending on the parameterization chosen. Since individual λ and ET must be greater than or equal to zero, multiplicative (proportional) covariate effects typically implemented for the pharmacokinetic (PK) parameters (in fact, any parameters that take on values greater than or equal to zero) can be applied (Equation 19).
For continuous covariates, a power model (Equation 20) or an exponential model (Equation 21) are commonly considered. An exponential model with logarithmically transformed covariate values and centering values is equivalent to the power model. The exponential model is less error‐prone than the power model since zero and negative covariate values are allowed; however, the exponential model tends to fail when individuals have covariate values that largely differ from the centering value.(19)Pcov=θP·fθcov,x (20)fθcov,x=xxcentθpower (21)fθcov,x=eθexp·x−xcentwhere θpower and θexp are the parameters defining the covariate‐parameter relationships respectively in the power model and the exponential model.
For categorical covariates, a model with a fractional difference to the most common category is often implemented (Equation 22).(22)fθcov,x=1ifx=xcent1+θshiftx≠xcent
For the specific case of θP being fixed to zero, an additive covariate effect model (Equation 23) should be considered instead of the multiplicative (proportional) covariate effect model (Equation 19). Fixing of selected θP to zero is common when the occurrence of a drug‐induced categorical adverse event is modeled for study arms in which a limited number of adverse events is expected to occur, such as in the placebo arm.(23)Pcov=θP+fθcov,x
For continuous covariates, linear model (Equation 12) and Emax model (Equation 13) are commonly evaluated, with the former case requiring special consideration to avoid Pcov resulting in a negative value. Likewise, for categorical covariates, Equation 14 can be implemented.
Covariate analysis for PSS, a probability parameter ranging between zero and one, can be performed in the same manner as the covariate analysis described for the probability parameters in DTMM. With respect to covariate analysis, a CTMM has the advantage of being able to account for covariate values that change during the time between observations, thereby allowing for a more accurate characterization of the effect of such time‐varying covariates. Overall, the CTMM is also more parsimonious than the DTMM in accounting for covariate effects (such as drug effects), which may be added on fewer parameters. For further parsimony, covariate parameters, such as θEC50, may be shared between different parameters with θEmax still being estimated independently.
IIV model
IIV model development strategy proceeds in the same manner as described in the DTMM section.
Starting the Markov chain and data requirement
There are various ways to initiate the Markov chain or the compartments. In some cases, this is straightforward, for instance when a sleep study starts at bedtime or when only patients with severe pruritus are randomized. Overall, the initiations may be based on inclusion and/or exclusion criteria, screening data, first available observation of the ordered‐categorical variable, or simply assumption. During initialization, the probability of the assumed state is set to one and the probability of other states is set to zero. Between initialization and/or observations, the probability in all compartments will re‐distribute as determined by the relevant transition rate constants. When the next observation is made, the probability in each compartment will be reset. However, at all times, the sum of the probabilities of all states is one. An example of probability changing over time, as described by a CTMM with an absorbing state for a typical subject is illustrated in Figure 5.
In terms of data programming, an additional requirement for the implementation of a CTMM in NONMEM is to reset the compartment amounts after each observation. The exact implementation of this step depends on various factors, including whether simulations need to be conducted and if there are compartments that should not be reset (e.g., PK compartments). A trimmed and annotated example of a NONMEM data set is shown in Supplementary Information S3 and the full data set is provided in Supplementary Information S8.
A discussion of NONMEM estimation methods with estimation setting examples, is provided in Supplementary Information S1.
Model evaluation
The recommended model evaluation for a CTMM are the same as those listed for a DTMM.
Applications
Examples of the use of a CTMM include modeling of improvements in rheumatoid arthritis, ^28^ muscle spasm grade, ^27^ gastrointestinal tablet transit, ^23^ tumor response, ^22^ proteinuria grade, ^29^ and fibrosis development in nonalcoholic steatohepatitis. ^30^
In terms of model application, the CTMM provides a natural framework to simulate the probability of each state, either for a typical subject or for a population, with or without resetting the probabilities in compartments when observations are made. Plots of these simulation results, stratified or colored by the covariate of interest, for instance dosing regimens, are often useful to help address questions central to the analysis. One such example is illustrated in Figure 4. If needed, the CTMM may be applied to simulate the values of a longitudinal ordered‐categorical variable. The data may be simulated in the same manner as that shown for the creation of the VPC (see Supplementary Information S3 for the relevant code). However, such simulations are more computationally intensive compared to direct simulations of probabilities.
A CTMM offers a few advantages over a DTMM. Unlike the DTMM, the CTMM naturally includes the effect of time since the last observation on the probability of each state and it can easily be applied to model ordered‐categorical data with non‐equidistant time intervals. Likewise, the CTMM can be applied when the aim is to simulate data with non‐equidistant time intervals or equidistant intervals that are different compared to those observed. As detailed earlier, unlike the DTMM, the CTMM is also able to account for time‐varying covariates with values that change during the time between observations. In addition, the CTMM is also more parsimonious than the DTMM in the following for systems with more than two states, excluding the absorbing state, and when accounting for covariate effects (such as drug effects) that may be added on fewer parameters.
Nevertheless, a few disadvantages are associated with the CTMM compared to the DTMM. The CTMM is often implemented as a system of ODEs, hence necessitating the use of the ODE solver, which usually translates into a longer time for parameter estimation. In addition, the implementation of the CTMM requires additional data programming for re‐initialization of compartments after observations. Finally, the CTMM is not as straightforward and easy to communicate to non‐pharmacometricians as the DTMM.
Minimal continuous‐time Markov model (mCTMM)
The mCTMM is a special case of CTMM. If the CTMM is reparameterized using ET and PSS, the number of ET to be estimated is the total number of states minus one. In mCTMM, the ET and PSS parameterization is typically used with the added assumption that the transition rate between two consecutive states is independent of the state. Consequently, only one ET needs to be estimated and it is somewhat an average of the multiple ETs in the original CTMM; therefore, the estimated ET is heuristically referred to as the mean ET (MET). The mCTMM is fully defined by Equation 12 and by assuming MET=ET12=…=ETK−1K. The mCTMM should be considered for cases in which a model more parsimonious than a CTMM is needed and when the assumption inherent to mCTMM ^31^ , ^32^ is fulfilled. Considering that a mCTMM is a nested version of a CTMM, a likelihood ratio test can be conducted to compare mCTMM to CTMM as a way to verify the MET=ET12=⋯=ETK−1K assumption.
Hidden Markov model
Theory and basic model structure
The Markov models described so far in this tutorial can model processes, or transitions, undergone by the observed data in a direct fashion. However, there may be situations in which these processes are underlying, thus influencing the observed data in an indirect manner. In these cases, hidden Markov models can be useful and they have been found to be powerful in explaining patterns' changes, attributed to switches in hidden states. ^33^ Hidden Markov models are hence a class of models characterizing the relationship between observed and hidden variables, the latter representing for example an underlying and unmeasurable disease status.
The implementation of the hidden Markov models is similar to that of the standard Markov models described above, with the difference that the Markov model will be complemented by a function describing the observed variable, conditioned on the predicted current Markov state. The estimation, adapted to the DTMM, uses the forward algorithm, summing over all the probabilities of each state at each position, as outlined in Equation 24 where the number of states (i.e., number of Markov components) is assumed to be 2.(24)Lj=∑12Pj P10=δ1·fθ,x P2o=δ1·fθ,x P1j=π11·P1j−1Lj−1+π21·P2j−1Lj−1·fθ,x P2j=π12·P1j−1Lj−1+π22·P2j−1Lj−1·fθ,x
Here the likelihood is noted L and the number of states is two. The probability to start at time 0 in either state is determined by the emission probabilities P10 and P2o, and the probability to transition at time j is defined by the probabilities P1j and P2j. fθ,x is a function describing the effect of various predictors on the relevant term.
The most likely hidden states chain can be retrieved in a post hoc fashion, using the Viterbi algorithm. Relevant NONMEM code is provided in Supplementary Information S6.
In the last decade, Markov models have also been brought to population modeling with the addition of stochasticity, that is, mixed hidden Markov models (MHMMs). Subsequently, they have been extended to handle two correlated (continuous) variables, that is, bivariate MHMM. ^34^
Applications
Examples of the use of hidden Markov models include modeling of epilepsy seizures, ^35^ migraine attacks, ^36^ sleep in rats, ^37^ and development of immunogenicity. ^38^
Item‐response theory (IRT) model with Markov sub‐models
Theory and basic model structure
Another type of data presenting a categorical aspect and possible candidate for Markov modeling is composite score data, where a total score is composed of multiple sub‐scores. These data can be instruments to measure the evolution of a disease and are largely represented by the family of patient‐reported outcomes (PRO). In this case, the observations are usually a series of categorical data that are best described by individual sub‐models joined by a common underlying latent variable, an approach called IRT modeling. ^39^ Current published literature on the applications of this approach mostly involve PO models and assume the independence between observations. However, there may be cases in which such as assumption is violated and, by ignoring the Markovian patterns, the risk of model misspecification and an inflated number of score changes are increased. These cases typically occur when frequent observations are collected, which is common nowadays since the use of digital technology and access to electronic devices enables the recording of health status on a frequent basis. Therefore, it is advisable to reflect on the dependence between subsequent scores in the IRT model by utilizing Markov models as sub‐models (DTMM, CTMM, or mCTMM). Using mCTMM as an example, one can consider the estimation of the steady‐state probabilities PSS in IRT modeling, which is the steady‐state probability for subject i to have a score of less than or equal to k for item j whereby 0<k≤K. Instead of estimating the PSS directly as shown in Equation 12, PSS in IRT modeling is formulated as a function of a latent variable for subject i (Di), for instance, a variable reflecting the ability of subject i to perform a certain (25)PSSYij≤k=11+e−aj·bj,k−Di PSSYij≤k−1=11+e−aj·bj,k−1−Di PSSYij=k=PSSYij≤kifk=1PSSYij≤k−PSSYij≤k−1if1<k<K1−PSSYij≤k−1ifk=Kwhere aj is the discrimination parameter (which is related to the slope of the expit function) and bj,k is the difficulty parameter (which corresponds to where the inflection in the expit function is centered on the x axis) for the item j.
Applications
To the best of our knowledge, at the time of writing, the only example of the use of IRT with Markov models as sub‐models in the pharmacometrics literature is in chronic obstructive pulmonary disease where the EXACT® score was modeled with an mCTMM. ^32^
DISCUSSION
Markov modeling or modeling of Markov elements is often necessary for ordered or non‐ordered‐categorical data (including cases where dropout data are present) when observations are frequent (leading to dependence between observations), or when there are many consecutive same‐state observations. These cases are often not well described, predicted, or simulated, especially at the individual level, unless the Markovian nature of the data is considered.
Figure 6 shows a decision tree for choosing the type of Markov model suitable to the analysis of ordered‐categorical data with Markovian features. This decision tree is meant to be used only as a general guidance since there may be other factors emerging during modeling that require pharmacometricians to adapt the Markov modeling strategy, for instance, parameter identifiability, practicalities around run time, and model fit.
In theoretical terms, it is always recommended to model the true nature of the data; however, a low occurrence of data in some categories may cause model uncertainties. Therefore, prior to initiating Markov modeling, the observed ordered‐categorical data should be scrutinized for resorting to lump consecutive categories. For example, should the absolute number of observations in a certain state be low or derive from very few subjects, adjacent categories can be lumped together, for instance, moderate nausea and severe nausea. In addition, the number of transitions to model may be reduced by considering whether all transitions to and from a state are possible. To facilitate this evaluation, modelers are suggested to generate a table summarizing the number and percentage of observations in each state and the unique transitions. It is worth noting that estimations of transition probabilities in a DTMM are directly informed by the exact transitions of interest, whereas for a CTMM, the estimations of the rate constants are informed by multiple transitions. It is worth noting that in CTMM transition rate parameters can be estimated even if the number of the specific transition is limited; for instance, the estimation of λ12 and λ23 can be informed mainly by 1‐to‐3 transitions even with limited 1‐to‐2 and 2‐to‐3 transitions. In addition, Markov model parameters may be fixed, most commonly to zero. This decision can be based on the aforementioned criteria, on the graphical analysis and/or understanding of the mechanism of drug action and underpinning physiology, or on the importance of the relevant transitions. Further cases in which relevant Markov model parameters may be fixed to zero include a negligible number of transitions, leading to a final estimate close to zero and a statistically non‐significant change in the objective function value between the Markov model with the parameter fixed and the model with the parameter estimated. Moreover, simplifications of the Markov model by lumping adjacent categories, reducing number of transitions to model, or fixing parameter values may be explored during model development if the Markov model estimated appears to be unstable and/or if parameter estimates are associated with poor precision; nevertheless, the model should not be over‐simplified to the extent that its application is impaired. Imprecision in some parameter estimates may be tolerated if the relevant transitions are of interest and the uncertainty in the parameter estimates may be accounted for in decision‐making during model application. In some cases, the data set available for Markov modeling is extremely large due to very frequent observations. Despite it is desirable to use all available data for model development, unless estimation time is prohibitive, in some cases the model development may need to be performed on a subset of the data. In such cases, well‐defined and unbiased data selection should be used (e.g., including only daily rather than hourly data, e.g., using the exact observed value at selected fixed timepoints or using summary‐level data such as maximum daily AE score). Nonetheless, the exact approach used requires careful considerations with views from different functions that are accounted for. Furthermore, when possible, the final model should preferably be re‐estimated on the full data.
Higher‐order Markov models, taking into account more than one preceding observation, should theoretically always provide a better description of the data compared to first‐order Markov models, due to the added predictive value from more “history”. However, higher‐order Markov models require the estimation of many more parameters, which may not always be supported by the data at hand, and their implementation is more complex as well as the results more difficult to communicate. Generally, the first‐order Markov models are often sufficient.
The Markov models come in different forms with adaptations considered to reach a fit‐for‐purpose model. The CTMM is typically implemented with transition rate constants estimated only between adjacent states. However, the inclusion of non‐adjacent direct transitions in CTMM has also been reported. ^27^ , ^40^ These transitions between non‐adjacent states can be thought of as non‐Markovian in nature and they represent a relaxation of the traditional CTMM stating that all transitions are necessarily Markov. The inclusion of transitions between non‐adjacent states should only be considered as a secondary option to reach an adequate model fit, after establishing that alternative model structures, and in particular the IIV model, have been sufficiently explored.
Multiple absorbing states may be considered. ^22^ , ^23^ The CTMM may also be adapted to include multiple transit compartments ^41^ to represent a certain ordered‐categorical value when there is a minimum residence time associated with the value. In one example, published by Bergstrand et al., a CTMM was used to describe a gastrointestinal tablet a chain of transit compartments was included to describe a prolonged residence time in the proximal and distal small intestine. ^23^
Bi‐modality for transition probabilities in a DTMM or transition rate constants in a CTMM, which is not explained by covariates, may be modeled using a mixture function. Alternatively, a two‐part Markov model ^42^ , ^43^ may be in the first part, a logistic model is used to describe and predict the dichotomous variable defining whether a subject experienced at least an event or no event; in the second part, a Markov model is developed to characterize the ordered‐categorical data with Markovian features only based on data from subjects who experienced at least one event. In this case, the logistic model and the Markov model are used in sequence to generate realistic simulations.
Recently, a new group of models referred to as multistate models, has been applied to describe the therapeutic response in oncology. ^44^ , ^45^ , ^46^ , ^47^ In the authors' opinion, these models are a form of the CTMM. However, some of these published models feature examples of transitions with a time‐dependent hazard (Gompertz or Weibull distribution) in a way that, to the authors' knowledge, has not been described for other nonlinear mixed‐effects CTMM models. ^44^ , ^45^ , ^46^
Many clinical studies incorporate adaptive interventions, such as dosing adjustments based on observed treatment response (efficacy and/or adverse reactions). This is a common practice in chemotherapy, where doses may be modified in response to life‐threatening AEs like thrombocytopenia, neutropenia, and anemia. ^48^ In other areas, such as diabetes management with insulin, the dose is titrated to achieve a desired target response. ^49^ The evaluation of exposure‐response relationships based on data from studies with adaptive dose adjustments requires careful consideration. Empirical exposure‐response relationships derived from steady‐state exposure and response may not accurately reflect the true causal relationship between exposure and response. ^50^ Longitudinal pharmacokinetic‐pharmacodynamic modeling, under appropriate conditions, can accurately characterize the underlying exposure‐response relationship. ^51^ Markov models, which focus on transitions between different response states, are less susceptible to bias introduced by dose adjustments if the true sequence of events is recorded in the analysis data set (e.g., first transition from No‐AE to AE, followed by dose reduction). Logistic regression models, even when incorporating Markov elements, are not as robust to these biases. Markov models are also valuable for modeling dosing adjustments that do not follow a strict deterministic algorithm. Clinical study protocols often allow dose adjustments “at the discretion of the investigator”, meaning that the investigator (possibly in collaboration with the patient) determines whether certain dose adjustments are necessary. In such cases, a model describing the probability of an adjusted dose at the next occasion is crucial for conducting realistic clinical trial simulations.
In terms of the stochastic aspect of the Markov models, the residual unexplained variability is not relevant since the likelihood of the data is modeled in place of the data themselves. Given different observed profiles between subjects, for instance, in their response to treatment, the estimation of the variance of IIV in one or more parameters may be possible. Depending on the available data, the variance of IIV may be associated with poor precision and its estimation may be associated with very long estimation time. In these cases, a Markov model without IIV or with limited inclusion of IIV, as those published by Xu et al. ^40^ and by Svensson et al., ^22^ may be considered and may already be fit for the intended purpose of the model application.
In the PO model with or without Markov element(s) and DTMM, it is often desirable to derive the probability term (e.g., those that include an absorbing state) that is somewhat distinct from other probability terms. The decision on which probability terms should be estimated and which should be derived can also be informed by which parameters are affected by covariates. Usually, the probability for those categories on which covariate effects are tested is estimated instead of being derived. Furthermore, the probability terms for transitions with more observations are likely to be estimated with better precision than those for transitions with less observations in the analysis data set.
An alternative to Markov modeling is the repeated time‐to‐event (RTTE) analysis, where transitions to different states are considered as events, and specific hazard parameters are estimated for the different transitions. The RTTE models resemble the CTMM ones except that time‐to‐event is considered and that the hazard function is estimated instead of the transition rate constants parameters. One example is the application of RTTE analysis to model the repeated time‐to‐start and repeated time‐to‐end of migraine episode, conditional on the current migraine severity status. ^52^
Ignoring Markov properties can have various consequences. ^53^ In terms of parameter estimation, failing to account for Markov properties often leads to an over‐valuation of the information content in the data; this in turn, makes the application of hypothesis testing inappropriate and leads to an overestimation of the precision of the parameter estimates. Most often, the IIV is over‐estimated and the structural model may be misspecified. In addition, models that failed to consider Markov properties will result in simulations with inflated number of transitions, inflated number of extreme value occurrences, and under‐predicted duration in the same state. This means that individual longitudinal pharmacodynamic endpoint values cannot be realistically simulated. Hence, it is recommended to always conduct a model evaluation according to the intended purpose of the model; for instance, one can generate VPCs evaluating the model's ability to predict the selected metric(s) of interest (see Figure 3). Finally, the optimal design results may also be inaccurate if the Markov properties are not accounted for. ^54^ , ^55^
CONCLUSIONS
In summary, Markov models (or Markov elements) are useful to handle dependencies between observations. Markov models can be applied to different types of data with Markovian features and Markov elements can be added to several types of models. Overall, the CTMM approach appears to be most widely applicable as it is robust to different assumptions, relatively parsimonious, and works generally well.
AUTHOR CONTRIBUTIONS
Q.X.O., E.P., and M.B. wrote the manuscript.
FUNDING INFORMATION
No funding was received for this work.
CONFLICT OF INTEREST STATEMENT
The authors declared no competing interests for this work.