Optimal designs for epidemiologic longitudinal studies with binary outcomes

Abstract

Alternating presence and absence of a medical condition in human subjects is often modelled as an outcome of underlying process dynamics. Longitudinal studies provide important insights into research questions involving such dynamics. This article concerns optimal designs for studies in which the dynamics are modelled as a binary continuous-time Markov process. Either one or both the transition rate parameters in the model are to be estimated with maximum precision from a sequence of observations made at discrete times on a number of subjects. The design questions concern the choice of time interval between observations, the initial state of each subject and the choice between number of subjects versus repeated observations per subject. Sequential designs are considered due to dependence of the designs on the model parameters. The optimal time spacing can be approximated by the reciprocal of the sum of the two rates. The initial distribution of the study subjects should be taken into account when relatively few repeated samples per subject are to be collected. A study with a reasonably large size should be designed in more than one phase because there are then enough observations to be spent in the first phase to revise the time spacing for the subsequent phases.

Keywords

binary outcome discrete observation Fisher information longitudinal data optimal design sequential designs

1 Introduction

Continuous-time stochastic processes with a binary state space are often used to model the alternating presence and absence of medical conditions in human subjects. For example, the within-host dynamics of a recurrent asymptomatic infection can be described in terms of repeated transitions between states of being colonised and non-colonised.^1,2 Likewise, some parasites are present as repeated presence and absence of the infection.³ When the transition rates between the two states are constant, the process is Markovian and also known as the alternating Poisson process.⁴ An obvious inference problem concerns the estimation of the rates from longitudinal observations of the process.

In practice, it is not always possible to observe the actual transitions between the two states continuously in time. A study designed to estimate the transition rates usually gives discrete-time observations of the underlying process for each study subject, i.e. the states of the subject i at a number of (pre)fixed times τ_ij, j = 1,…, K. The estimation of the two rates in a binary model using such longitudinal data has been addressed previously by several authors.^3–7

The precision in estimation of the rates is of particular interest. Using the maximum likelihood method, the precision is determined by the Fisher information which in turn depends on the number of subjects in the study (N), the number of observations per subject (K), the observation times ({#x003C4;_ij}) and the initial distribution (ρ₁).^5,8 While some of these quantities may have to be fixed in advance by practical constraints, others remain subject to optimisation when designing the study. In the context of discrete observation of a continuous-time process, there is particular interest in the optimal choice of the times of observation.^8,9

Albert and Brown⁸ considered the choice of optimal sampling times for the binary Markov process in connection with estimation of the two rates as well as their ratio. For an arbitrary initial distribution, they studied optimal observation times when K = 2, and when K > 2, their analysis was confined to the stationary initial distribution and a special choice with ρ₁ = 0. The efficiency of discrete observation with respect to continuous observation was considered by Albert and Brown⁸ and Cook.¹⁰

There are several related models with similar design questions, based on discrete observation of an underlying continuous-time process. These include the simple birth process,¹¹ the simple death process¹² and the simple birth-death process.¹³ Processes with multiple states have been considered, e.g. for the analysis of epidemic data.¹⁴ In all these studies, optimisation is based on the Fisher information about the unknown rates. Bayesian approaches to optimal design have also been considered in the context of binary sequence data⁹ and the susceptible-infected epidemic processes.¹⁵

Designs based on the Fisher information depend on unknown values of the model parameters. Approaches to overcome the parameter dependency include the use of an initial best guess for the parameter values,¹⁶ sequential procedures to update parameter values using previous observations,¹⁶ use of a Bayesian prior distribution for the parameters^9,15 and application of the maximin approach, where the information is maximised over a parametric region.¹⁷

In this article, we discuss designing studies with discrete-time observations of a binary Markov process based on the Fisher information. Our main focus is in determining the sampling times ({#x003C4;_ij}) for a given number of observations per subject (K). The analysis generalises the corresponding results in Albert and Brown⁸ to the case K > 2. In addition, we determine optimal initial distributions (ρ₁) and consider designs for the estimation of one of the two rates in the model. We also consider optimal designs in sequential studies in which sampling times are revised in the course of the study, and address the optimal split of observations between the two phases in a sequential study.

The structure of this article is as follows. The next section defines the binary Markov model and describes the maximum likelihood estimation of the parameters. Section 3 reviews design questions and methodological approaches relevant for discretely observed processes. Based on these, Section 4 presents different scenarios for studies with binary outcomes and guidelines for their design and analysis. Illustrative examples on estimation of acquisition and clearance rates of a bacterial pathogen, Streptococcus pneumoniae, are given in Section 5. Section 6 provides the concluding remarks.

2 Model and maximum likelihood estimation

Let {#x003B6;(τ), τ ≥ 0} be a continuous-time process with state-space {0, 1}, where 1 and 0 denote the presence and absence of the condition in question. Let the process ζ(τ) be governed by two constant rates (per time unit):

λ_{rs} = lim_{Δ τ \to 0} P (ζ (τ + Δ τ) = s | ζ (τ) = r) Δ τ, r, s \in {0, 1}, r \neq s .

Let X_ij be a discrete-time observation ζ(τ_ij) of subject i at time τ_ij. The transition probabilities for two consecutive observations with time spacing t_ij = τ_ij − τ_i(j−1) are

P_{rs} (t_{ij}) = P (X_{ij} = s | X_{i (j - 1)} = r; λ) = λ_{rs} λ_{01} + λ_{10} (1 - \exp - (λ_{01} + λ_{10}) t_{ij}),

where r, s ∈ {0, 1}, r ≠ s, λ = (λ₀₁, λ₁₀) and P_rr(t_ij) = 1 − P_rs(t_ij). For N independent realisations of {#x003B6;(τ); τ ≥ 0} at K equally spaced times {X_ij, t_ij = t, i = 1,…, N, j = 1,…, K}, the likelihood function for λ is

Π_{i = 1}^{N} (P (X_{i 1}) Π_{j = 2}^{K} P (X_{ij} | X_{i (j - 1)}; λ)) = Π_{i = 1}^{N} P (X_{i 1}) L_{i} (λ; X_{i}) .

Note that X_ij depends only on X_i(j−1), since with constant rates {#x003B6;(τ), τ ≥ 0} is a Markov process. The rationale of restricting the analysis to equally spaced times is explained later in this article. The model of X_i1 is left unspecified here and the inference about λ is based on the likelihood

L (λ; X) = Π_{i = 1}^{N} L_{i} (λ; X_{i}) .

(1)

Denoting by $1_{rs}^{ij}$ the indicator of the event (X_i(j−1) = r, X_ij = s), r, s ∈ {0, 1}, i = 1,…, N, j = 1,…, K, the logarithmic likelihood function is

l (λ; X) = \sum_{i = 1}^{N} \log (L_{i} (λ; X_{i})) = \sum_{i = 1}^{N} \sum_{j = 2}^{K} \sum_{r, s = 0}^{1} 1_{rs}^{ij} \log (P_{rs} (t)) .

The maximum likelihood estimators of the transition probabilities are

{\overset{\land}{P}}_{01} (t) = N_{01} N_{01} + N_{00}, {\overset{\land}{P}}_{10} (t) = N_{10} N_{10} + N_{11},

where

N_{rs} = \sum_{i = 1}^{N} \sum_{j = 2}^{K} 1_{rs}^{ij}

. The maximum likelihood estimators of the corresponding rates are

{\overset{\land}{λ}}_{rs} = - {\overset{\land}{P}}_{rs} (t) t ({\overset{\land}{P}}_{01} (t) + {\overset{\land}{P}}_{10} (t)) \log (1 - {\overset{\land}{P}}_{01} (t) - {\overset{\land}{P}}_{10} (t)), r, s \in {0, 1}, r \neq s .

Asymptotically, the maximum likelihood estimates exist, and are normally distributed and unbiased. In finite sample estimation, some complications may occur as estimates and do not exist if ${\overset{\land}{P}}_{01} (t) + {\overset{\land}{P}}_{10} (t) \geq 1 .$ A probabilistic condition for this states that non-existence is likely if t(λ₀₁ + λ₁₀) is large.^4,5,7 Specifically, t(λ₀₁ + λ₁₀) is large if t is such that probability P₀₁(t) is close to the corresponding stationary probability, λ₀₁/(λ₀₁ + λ₁₀).

Finite samples can also lead to biased estimates, especially if the time spacing is far from optimal. The smallest positive values for ${\overset{\land}{P}}_{01} (t)$ and ${\overset{\land}{P}}_{10} (t)$ are 1/(1 + N₀₀) and 1/(1 + N₁₁), respectively. If the time spacing is so short that the true probabilities P₀₁(t) and P₁₀(t) are even smaller, the smallest non-zero estimates of the rates are higher than the true values. Similarly, if the time spacing is too long, the true transition probabilities P₀₁(t) and P₁₀(t) are close to the stationary probabilities but ${\overset{\land}{P}}_{01} (t)$ and ${\overset{\land}{P}}_{10} (t)$ are again limited by the sample size and the existence condition. If the estimates exist, they are too low. We refer to Brown et al.⁴ for a detailed discussion on the problem of bias in estimating the rates using finite samples.

The Fisher information matrix for λ in the binary Markov model is derived in Kalbfleisch and Lawless,⁵ and is given as follows

I (λ; X) = E (- H (λ; X)) = \sum_{i = 1}^{N} E (- \nabla_{λ} \nabla_{λ}^{T} \log (L_{i} (λ; X_{i}))) = \sum_{j = 2}^{K} \nabla_{λ} P T (t) diag (E (N_{0} (τ_{j - 1})) P_{01} (t) P_{00} (t), E (N_{1} (τ_{j - 1})) P_{10} (t) P_{11} (t)) \nabla_{λ} P (t),

(2)

where ▿_λ = (∂/∂λ₀₁, ∂/∂λ₁₀)^T is the partial differential operator,

\nabla_{λ} P (t) = (\partial P_{01} (t) \partial λ_{01} \partial P_{01} (t) \partial λ_{10} \partial P_{10} (t) \partial λ_{01} \partial P_{10} (t) \partial λ_{10}),

▿_λ P^T(t) is the transpose of the above matrix, diag(a, b) is a 2 × 2 diagonal matrix with diagonal elements a and b, and H(λ; X) is the Hessian matrix of the logarithmic likelihood function. The expected number of subjects in state r at time τ_j−1, E(N_r(τ_(j−1))) is given by

E (N_{r} (τ_{j - 1})) = N [P_{0 r} ((j - 1) t) (1 - ρ_{1}) + P_{1 r} ((j - 1) t) ρ_{1}],

(3)

where the initial distribution ρ₁ = P(X_i1 = 1), i = 1,…, N.

Recall that the asymptotic covariance matrix V(λ; X) of $({\overset{\land}{λ}}_{01}, {\overset{\land}{λ}}_{10})$ is the inverse of the information matrix I(λ; X) and that both matrices depend also on quantities N, K, t and ρ₁. In principle, the Fisher information could be derived for any continuously differentiable function of the model parameters. For instance, log-parametrisation could enhance the asymptotic properties of the estimators. In the case of a reasonably large sample size, the asymptotic performance of the maximum likelihood estimators is known to be satisfactory for absolute rates as well.

3 Design issues and approaches

Optimisation of designs is often based on some statistical criterion derived from the Fisher information. In this article, designs are determined under the A- and A_i-optimality criteria. The A-optimality criterion minimises the trace of the covariance matrix V(λ; X), i.e. the sum of the asymptotic variances of both parameter estimators. When the criterion is applied to only one of the model parameters (A_i-optimality), this reduces to the minimisation of the asymptotic variance of its estimator.

Consider a binary Markov model with two constant rates, λ₀₁ and λ₁₀, as the parameters of interest. A study is to be planned with a design d ∈ , where is a space of equidistant designs, i.e. each d ∈ is characterised by a constant time spacing between consecutive observations. This restriction on the design space is not too severe. Specifically, when the first observation follows the stationary distribution, the optimal equidistant design is more efficient than any non-equidistant one.⁸ When the first observation does not follow the stationary distribution, there are non-equidistant designs which are more efficient until stationarity is attained.¹² However, it is argued in this article that for the optimal equidistant design the process attains stationarity for small values of K.

Two designs d₁ and d₀ can be compared using their asymptotic relative efficiency (ARE). Under the A-optimality criterion, we define the ARE of d₁ with respect to d₀ as

ARE (d_{1}, d_{0}; λ, X) = 1 / tr (V_{d_{1}} (λ; X)) 1 / tr (V_{d_{0}} (λ; X)) = tr (V_{d_{0}} (λ; X)) tr (V_{d_{1}} (λ; X)),

where tr(V_d(λ; X)) is the trace of the covariance matrix under design d. If design d₁ is more efficient than design d₀, ARE(d₁, d₀; λ, X) is greater than 1, and vice versa.

In a two-phase study, N subjects have been followed and K₀ (≥ 2) observations per subject, $X_{K_{0}}^{obs} = (X_{ij}; j = 1, \dots, K_{0}, i = 1, \dots, N)$ , have already been made with time spacing t₀. Design questions for K future observations in the second phase include when and how the time spacing t should be revised and the selection of n subjects out of N. Denote the second-phase observations as $X_{K}^{new} = (X_{ij}; j = K_{0} + 1, \dots, K_{0} + K, i \in S)$ , where S = {i₁,…, i_n} ⊆ {1,…, N}.

The search for optimal t and S in the second phase can be based on the total information contained in the data $X = {X_{K_{0}}^{obs}, X_{K}^{new}}$ ,

I (λ; X) = E [E [- \partial 2 \log (L (λ; X_{K_{0}}^{obs})) \partial λ 2 + - \partial 2 \log (L (λ; X_{K}^{new} | X_{K_{0}}^{obs})) \partial λ 2 | X_{K_{0}}^{obs}]] = E (- H (λ; X_{K_{0}}^{obs})) + E [I (λ; X_{K}^{new} | X_{K_{0}}^{obs})],

where the outer expectation is with respect to the marginal distribution of

X_{K_{0}}^{obs} .

The first term in the equation is the information contained in the first-phase data, and the second term is the expected information contained in the second-phase data conditional on the first-phase data. Because the second phase is designed after the first-phase data are available, we take

X_{K_{0}}^{obs}

to be fixed and hence base the design of the second phase on the observed total information

I (λ; X) = - H (λ; X_{K_{0}}^{obs}) + I (λ; X_{K}^{new} | X_{K_{0}}^{obs}) .

The first term carries the influence of the observed information. The method including the observed information has previously been employed in Karvanen et al.¹⁸

In the above expressions, λ is the true parameter value which, however, is unknown. The design has to be based on a given value of λ which, after the first phase, we take to be the maximum likelihood estimate $\overset{\land}{λ}$ . The information used in the design of the second phase is thus

I (\overset{\land}{λ}; X) = - H (\overset{\land}{λ}; X_{K_{0}}^{obs}) + I (\overset{\land}{λ}; X_{K}^{new} | X_{K_{0}}^{obs}) .

(4)

Because of the Markov property, the conditional distribution of

X_{K}^{new}

depends only on X_{iK
₀}, i ∈ S, and hence the expectation (3) now reads as

E (N_{r} (τ_{(j - 1)})) = n [P_{0 r} ((j - 1) t) (1 - ρ_{1}^{new}) + P_{1 r} ((j - 1) t) ρ_{1}^{new}],

where

ρ_{1}^{new}

is the proportion of subjects in S with initial state 1.

Conducting a study in two phases gives the opportunity to improve the design by use of estimates $(\overset{\land}{λ})$ from the first phase but poses the practical question of when to revise the design. The optimal split between the two phases can be searched by maximising the ARE of the two-phase study with respect to the corresponding one-phase study. The evaluation of ARE then needs to take into account that the second-phase time spacing is a random variable (see the Appendix for computational details). The optimal split of observations between the two phases in a two-phase study is considered in Section 4.2. The above ideas can be generalised to an m-phase study in which data $X_{(m - 1)}^{obs}$ have been observed from m − 1 phases with time spacings t₀,…, t_m−2.

4 Optimal designs: scenarios and guidelines

4.1 One-phase studies

We consider scenarios in which the design needs to be based on hypothesised values of the two rates. Designs based on such hypothesised values are referred to as locally optimal.^17,19 We determine A-optimal and A₁-optimal (minimisation of variance of ${\overset{\land}{λ}}_{01}$ ) designs, corresponding to situations in which the interest is to estimate either both or one of the rates, respectively. Results for A₂-optimal designs (minimisation of variance of ${\overset{\land}{λ}}_{10}$ ) are obtained similarly by symmetry. We consider the dependence of the optimal time spacing t_opt on the initial distribution (ρ₁) and the number of samples per subject (K). In addition, we consider the joint optimality of t_opt and ρ₁ and discuss the influence of K and the number of subjects N on the precision of the rate estimators.

4.1.1 Optimal time spacing (t_opt) with the stationary initial distribution

When the binary process evolves in time, the proportion of subjects in state 1 approaches the corresponding stationary probability. Thus, the stationary distribution is a natural choice for an initial distribution ρ₁ for study designs in which no known factors have had influence on the evolution of the process.

In the case of stationary ρ₁, the optimal time spacing is not affected by K or N. Figure 1 shows that the optimal spacings under both A- and A₁-optimality criteria are almost equal and only depend on the ratio of the rates λ₁₀/λ₀₁ and the time scale. If the ratio is not far from 1, t_opt in the time units of 1/λ₀₁ is well approximated by 1/(1 + λ₁₀/λ₀₁). When the ratio approaches 0, t_opt becomes almost independent of the ratio and tends to 1.59 (in the units of 1/λ₀₁). This is the optimal sampling time for an exponential random variable with rate λ₀₁ and, curiously, also the optimal observation period for a simple death process with the same rate when state 1 is the absorbing state.^12,16 When the ratio λ₁₀/λ₀₁ approaches ∞, the results are symmetric for λ₁₀ and can be read from the figure by interchanging the roles of λ₀₁ and λ₁₀.

4.1.2 Optimal time spacing (t_opt) with an arbitrary initial distribution

In some cases the initial distribution may be known or approximately conjectured in advance. For instance, in studies of bacterial colonisation among newborns,¹ the study subjects have, by biological knowledge, a minimal probability for the presence of the bacterium at the initiation of the study. In such and similar situations, the initial distribution favours one of the states and is clearly not the stationary one.

In the case when ρ₁ is not stationary, t_opt depends on ρ₁, K and the optimality criterion. Under A₁-optimality (Figure 2, upper panels), t_opt increases monotonically with ρ₁ irrespective of the ratio λ₁₀/λ₀₁. This reflects the fact that learning about λ₀₁ requires the transition 1 → 0 to occur first if the initial state is 1. Under A-optimality (Figure 2, lower panels), t_opt is a convex function of ρ₁ because now both rates are of interest. Specifically, when the ratio is 1, t_opt becomes longer when ρ₁ approaches either 0 or 1 for similar reasons as under A₁-optimality. When the ratio is greater than 1, t_opt tends to be shorter with increasing ρ₁, and vice versa. Under both optimality criteria, the t_opt curves converge towards the line corresponding to t_opt under stationarity as K increases. Concrete examples follow in Section 5.

4.1.3 Joint optimal time spacing (t_opt) and initial distribution $(ρ_{1}^{opt})$

Under some scenarios, the initial distribution ρ₁ can be taken as a parameter subject to optimisation. This obviously requires some type of prior knowledge or control of ρ₁. For instance, in Hill et al.,¹ newborns starting pathogen-free were followed to study the acquisition rate (λ₀₁), in which case, the initial state is 0 (ρ₁ = 0). One could have alternatively initiated the study later after birth, in which case, the initial distribution would be most likely close to the stationary one. In another instance, if the states of the subjects are measured once, a selection of whom to follow can be made. The approach that the rate of clearance (λ₁₀) of bacterial colonisation would be learned well by following subjects that initiate as colonised (i.e. in state 1) was employed in O'Brien et al.²⁰ Finally, examples in which the choice of initial distribution can be made directly include animal studies,²¹ in which the pathogen is inoculated to the subjects, giving complete control over the initial state. Similar challenge studies with humans can be conducted with special consideration on ethics.²²

Panel A in Table 1 presents the optimal initial distribution,

ρ_{1}^{opt}

, when t and ρ₁ are optimised jointly. Under A-optimality,

ρ_{1}^{opt}

is determined by the ratio of the two rates. The optimal initial probability of the rarer state is close to 1 whenever K > 2. By contrast, under A₁-optimality (for λ₀₁)

ρ_{1}^{opt}

tends always towards 0, because starting from state 0 obviously gives the most information about rate λ₀₁. The impact of the initial distribution on ARE is significant only when K is small. For example, when λ₁₀/λ₀₁ = 1, the ARE of a design with any ρ₁ ∈ [0, 1] with respect to the optimal design with

ρ_{1}^{opt}

is greater than 0.5, 0.7, 0.8 and 0.9, for K = 3, 4, 5, 6, respectively.

Table 1.

Optimal initial distribution under A- and A₁-optimality. (Panel A) One-phase studies: The optimal initial probability of state 1 ( $ρ_{1}^{opt}$ ) for various values of samples per subject (K) and the ratio of the rates. (Panel B) Two-phase studies: The optimal proportion of subjects in state 1 ( $ρ_{1}^{new}$ ) for the second phase under six different estimates $({\overset{\land}{λ}}_{01}, {\overset{\land}{λ}}_{10})$ of the rates, their variances $Var (\overset{\land}{λ}) = (Var ({\overset{\land}{λ}}_{01}), Var ({\overset{\land}{λ}}_{10}))$ and covariance $Cov (\overset{\land}{λ})$

			K
	λ₁₀/λ₀₁	Var $(\overset{\land}{λ})$ ,Cov $(\overset{\land}{λ})$	2	3	4	5	6
Panel A: One-phase design
A-optimal	0.2		0.33	0.08	0.00	0.00	0.00
	1		0.50	0.50	0.50	0.50	0.50
	5		0.67	0.92	1.00	1.00	1.00
A₁-optimal	0.2		0.16	0.00	0.00	0.00	0.00
	1		0.17	0.05	0.01	0.00	0.00
	5		0.17	0.17	0.18	0.18	0.18
Panel B: Two-phase design ( ${\overset{\land}{λ}}_{01} = 1.0)$
A-optimal	0.2	(0.004,0.005),0.002	0.50	0.21	0.00	0.00	0.00
	1	(0.008,0.043),0.014	0.92	0.81	0.79	0.79	0.79
	5	(0.023,0.710),0.103	1.00	1.00	1.00	1.00	1.00
A₁-optimal	0.2	(0.004,0.005),0.002	0.18	0.00	0.00	0.00	0.00
	1	(0.008,0.043),0.014	0.32	0.15	0.00	0.00	0.00
	5	(0.023,0.713),0.103	0.33	0.27	0.25	0.24	0.24
A-optimal	0.2	(0.030,0.001),0.002	0.00	0.00	0.00	0.00	0.00
	1	(0.045,0.009),0.015	0.08	0.19	0.21	0.21	0.21
	5	(0.110,0.140),0.100	0.42	0.72	1.00	1.00	1.00
A₁-optimal	0.2	(0.030,0.001),0.002	0.00	0.00	0.00	0.00	0.00
	1	(0.045,0.009),0.015	0.00	0.00	0.00	0.00	0.00
	5	(0.110,0.140),0.100	0.00	0.00	0.00	0.00	0.00

According to these results, the approach in O'Brien et al.²⁰ to restrict the follow-up to those in state 1 to learn about λ₁₀ is appropriate. As a practical illustration, we note that in Trotter and Gay,² two types of studies were employed in which the initial distribution was either based on a random sample of subjects, or selection of subjects was made towards initial state 1. The upper limit of a 95% confidence interval for λ₀₁ from a study in which the subjects started from state 1 with three repeated observations was 110 times the point estimate. By contrast, for λ₁₀, the corresponding ratio was only 1.7. Such a difference did not occur in those studies in which the initial proportion of subjects in state 1 was ∼0.4.

4.1.4 Optimal number of repeated samples per subject (K)

In this section, we study the choice of the number of repeated samples (K) under the assumption that the total number of observations (NK) is fixed. Such a question might arise, for instance, when the initial distribution departs from the stationary one and the cost of data collection depends heavily on the cost per sample rather than the cost of introducing new subjects to the study.

With a fixed total number of observations (NK = 105), we compare two studies d₁ and d₂ with K = 5 and K = 7, respectively, to a reference study d₀ with K = 3. Figure 3 presents ARE(d₁, d₀; λ, X) and ARE(d₂, d₀; λ, X) as functions of ρ₁ for different values of the ratio λ₁₀/λ₀₁. Note that ARE(d₂, d₁; λ, X) = ARE(d₂, d₀; λ, X)/ARE(d₁, d₀; λ, X). Figure 3(A) shows that ARE(d₁, d₀; λ, X) > 1 and ARE(d₂, d₁; λ, X) > 1 for any ρ₁ and so it is more efficient to take more samples per subject. The same applies when the ratio is below 1 and ρ₁ tends towards the corresponding stationary value and above (Figure 3B and C). Only when ρ₁ tends to 0, the optimal initial distribution according to Panel A in Table 1, it is more efficient to take more subjects than repeated samples per subject. The figures corresponding to ratios greater than 1 are mirror images of those for ratios less than 1.

The recommendation about favouring repeated samples over increasing the number of subjects obviously relies on the setting in which the rates can be taken to be constant over time and across study subjects. In practice, the assumptions might not hold. Most notably, there could be heterogeneity across subjects in the rates. In such a case, if the simple model with constant rates is used, a reasonable number of subjects, depending on the setting, would be needed to average over such variation. Other practical issues potentially affecting the choice between N and K are the constancy of the rates over time and the feasibility of a study with an extended follow-up per subject.

4.2 Two-phase studies and sequential designs

Given the possibility that the rates can be estimated during the study, sequential designs can be applied to improve the efficiency of estimation. In a two-phase study, the design of the first phase is similar to a one-phase study. Here, we consider optimal designs for the second phase using the observed total information (4).

4.2.1 Optimal time spacing (t_opt) for the second phase

We simulated first-phase observations

X_{K_{0}}^{obs}

, and for each realisation, searched t_opt for the second phase. In most situations, t_opt was well approximated by the reciprocal of the sum of the estimated rates, similar to one-phase studies (Section 4.1). Only when one rate was estimated with a clearly higher precision than the other one in the first phase, and only a small number of repeated samples per subject was to be collected in the second phase (K < 3), there was sometimes, a tendency for t_opt to be large, i.e. towards independent observations in the second phase (Table 2). This reflects the fact that independent observations give the best information about the ratio of the rates.⁸

Table 2.

Influence of the observed information from the first phase on the time spacing in the second phase when one of the rates is estimated with higher precision than the other. Scenario: Ten subjects are to be followed in two phases. The estimates from the first phase are ${\overset{\land}{λ}}_{01} = {\overset{\land}{λ}}_{10} = 1$ with precisions 1 and 14. The table presents the optimal time spacing for K = 3 and different values of $ρ_{1}^{new}$ (the initial distribution for the second phase). For comparison, the optimal time spacing without including the observed information is also shown ( $t'_{opt}$ )

$ρ_{1}^{new}$	0.0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1.0
t _opt	0.62	0.64	0.66	0.69	0.74	0.81	∞	∞	∞	∞	∞
$t'_{opt}$	0.62	0.56	0.53	0.51	0.50	0.49	0.50	0.51	0.53	0.56	0.62

4.2.2 Joint optimal time spacing (t_opt) and initial distribution $(ρ_{1}^{opt})$ for the second phase

In two-phase designs, the choice of $ρ_{1}^{new}$ for the second phase can be a part of the design question as the subjects have been under observation. In practice, the optimal value of $ρ_{1}^{new}$ cannot be always employed but only pursued, as the subset of subjects cannot be selected arbitrarily.

In a one-phase study, the optimal initial distribution is determined by the ratio of the two rates and the optimality criterion (Panel A in Table 1). In a two-phase study, the optimal initial distribution $ρ_{1}^{new}$ depends also on the observed information (Panel B in Table 1), i.e. the precision with which the rates have been estimated from the first phase. Under A-optimality, when the observed data are such that one of the rates (e.g. λ₀₁) is estimated with a lower precision in the first phase, the state from which that rate is learned better (0 for λ₀₁) should be favoured as the initial state for the subjects in the second phase. Under A₁-optimality, the optimal $ρ_{1}^{new}$ is again close to 0, regardless of the observed data or ratio of the two rates.

4.2.3 When the time spacing should be revised?

The robustness of a first-phase design to discrepancies between the initial guess of the rates and their true values can be assessed with the ARE. Figure 4 shows ARE(d_guess, d_true; λ^true, X), evaluated at the true values of the rates λ^true, where d_a denotes the optimal design with time spacing $t_{opt}^{a}$ obtained using the parameter values λ^a. As a function of $t_{opt}^{guess} / t_{opt}^{true} \approx (λ_{01}^{true} + λ_{10}^{true}) / (λ_{01}^{guess} + λ_{10}^{guess})$ , the ARE does not depend on the actual values of the true rates and is almost independent of their ratio and ρ₁. For example, in the stationary case, doubling the spacing from the truly optimal one reduces ARE from 1 to 0.6 whereas halving reduces it to 0.78. These scenarios correspond to guessing the sum of the rates approximately 0.5 or 2 times the true sum, respectively. This observation gives rise to the question about how efficiency of estimation could be improved if there is no information to limit the rates within a reasonable range a priori.

For simplicity, assume that both ρ₁ for the first phase and $ρ_{1}^{new}$ for the second one are stationary. The precision of the estimators in each phase then depends only on the number of transitions in that phase. To study when the time spacing should be revised, two-phase studies with the same size, total number of transitions N(K₀ − 1) + nK, but with different sizes of the first phase (N(K₀ − 1)) are compared. A study in which all (N(K₀ − 1) + nK) transitions are observed in a single phase is taken as a reference study. For detailed methods, see the Appendix.

Figure 5 shows the ARE of a two-phase study with respect to the one-phase study as a function of the size of the first phase and the ratio of the time spacing $t_{opt}^{guess}$ in the first phase to the truly optimal spacing $t_{opt}^{true}$ . The study size was assumed to be 500 and the size of the first phase was at least 4% of it. When $t_{opt}^{guess} / t_{opt}^{true} \leq 0.5$ or $t_{opt}^{guess} / t_{opt}^{true} \geq 2$ , the ARE≥ 1 for nearly all sizes of the first phase and so conducting the study in two phases is more efficient, irrespective of the ratio of the two rates. Then, one should spend approximately 30–40% of the observations in the first phase. For larger study sizes, a smaller proportion in the first phase is optimal, e.g. 20% for a study size of 1000 transitions. By contrast, for smaller study sizes or ratios $t_{opt}^{guess} / t_{opt}^{true}$ further from unity, a larger proportion in the first phase is optimal.

5 Examples on colonisation of S. pneumoniae

S. pneumoniae (pneumococcus) is a bacterium that colonises the nasopharynx of healthy individuals. Colonisation is a pre-requisite for pneumococcal disease and thus of significant clinical importance. It is common in young children and especially newborns have been studied to determine the rate of acquisition (λ₀₁) during the first year of life.^23,24,1 The time spacings between repeated measurements for this purpose have been roughly 1, 0.5 and 0.5 months in Syrjänen et al.,²³ Granat et al.²⁴ and Hill et al.,¹ respectively.

Consider the optimal design for a study to estimate a constant rate of pneumococcal acquisition (λ₀₁) in newborns under three different scenarios. The scenarios are characterised by the prevalence of colonisation, π₁ = 0.2, 0.5, 0.8, reflecting epidemiologically different conditions.^23,24,1 The clearance rate (λ₁₀) for pneumococcus has been estimated around 0.5–1 per month in several different settings. To design the study, we thus adopt value λ₁₀ = 1 for the clearance rate and determine the acquisition rate through the stationary probability λ₀₁/(λ₀₁ + λ₁₀) = π₁, which yields λ₀₁ = 0.25, λ₀₁ = 1 and λ₀₁ = 4 (per month) for the three scenarios, respectively.

The corresponding optimal time spacings, would the initial distribution (ρ₁) be stationary, are obtained from Figure 1, being 0.93, 0.50 and 0.23 (months) in the three scenarios, respectively. Because newborns can be taken to be in state 0 at birth (ρ₁ = 0), and the rate of interest is that of acquisition, the optimal time spacings are slightly shorter. For example, if four samples per individual are to be taken (K = 4), the optimal time spacings are 0.7, 0.4 and 0.15 (months), respectively (Figure 2). Efficiency gain from adjusting the time spacing with respect to the initial distribution is less than 20% (Figure 4). These time spacings are shorter compared to those in Syrjänen et al.,²³ Granat et al.²⁴ and Hill et al.,¹ but the efficiency loss is significant only for the third scenario with π₁ = 0.8 (Figure 4).

As 0 is the optimal initial state to estimate rate λ₀₁ (Panel A in Table 1), a study conducted among newborns improves the precision in estimation of the acquisition rate. By contrast, the clearance rate (λ₁₀) might be estimated with low precision. Suppose the study was to be continued in order to improve the precision of the clearance rate. For this purpose, a good initial state is 1 ( $ρ_{1}^{new}$ close to 1). Especially in scenario 1, the prevalence of state 1 at the Kth measurement (with K = 4, say) might be low and the choice of ρ₁ close to 1 would thus not be possible. If several repeated samples are still to be collected in the second phase (Section 4.2), the individuals should be followed up with a slightly longer time spacing (Figure 2). Alternatively, if only one or two samples in the second phase are available, they should be taken as far as possible to ensure stationarity. If a large proportion of individuals are in state 1 at the Kth measurement, which might be the case especially in scenario 3, the best option would be to follow up only those individuals who are in state 1 (B in Table 1).

It is known that different serotypes of S. pneumoniae have varying acquisition rates and possibly different clearance rates. If the interest is in estimation of serotype-specific acquisition and clearance rates, a multi-state model would be needed. These types of models are not considered in this article.

6 Concluding remarks

We have studied factors that influence the optimal designs for discretely observed binary Markov processes. The time spacing in a one-phase design with stationary initial distribution ρ₁ can be based on the approximation $t_{opt}^{appr} = (λ_{01} + λ_{10}) - 1$ . When the ratio of the rates tends to 0 or ∞, the optimal time spacing actually tends to 1.59/max(λ₁₀, λ₀₁) (Figure 1). However, the ratio of $t_{opt}^{appr}$ to this is always between 1/1.59 and 1, and it then follows that the loss in efficiency in using the approximate time spacing is marginal (Figure 4).

The time spacing $t_{opt}^{appr}$ seems to be a reasonable approximation even for a non-stationary initial distributions because the efficiency gain using optimal time spacing with a non-stationary initial distribution is marginal. The process also approaches stationarity with a small number of observations per subject (K). This follows from the fact that the distance $| π_{1} - ρ_{1} ({Kt}_{opt}^{appr}) |$ between the stationary probability of state 1, π₁, and the marginal probability of state 1 at the Kth observation using time spacing $t_{opt}^{appr}$ is at most $\exp - K (λ_{01} + λ_{10}) t_{opt}^{appr} = \exp - K .$

For the stationary initial distribution, equidistant designs are efficient among all discrete designs.⁸ As a consequence of the fast convergence to stationarity, this also holds for most non-stationary situations. The exceptions can be viewed as special cases of selecting the initial distribution. For example, if all subjects would be known to be in state 0, and λ₁₀/λ₀₁ was notably less than 1, it would be efficient to let the process evolve towards stationarity before collecting any samples. In most cases, however, restriction to the equidistant designs appears justifiable.

Often some knowledge about the studied phenomenon is available a priori, so that a range for the rates can be conjectured. Otherwise two-phase studies should be considered. Enough samples should then be collected in the first phase to ensure high probability for improved design for the second phase. The choice of optimal time spacing and initial distribution after the first phase of a sequential design is based on the estimates from the previous phases and their precision. In most cases, the optimal time spacing can be approximated by the reciprocal of the sum of the rates. The precision with which the rates are estimated from previous phases needs to be taken into account in the choice of $ρ_{1}^{new}$ (Panel B in Table 1).

We have considered only one of the simplest models, in which transitions between two states are governed by constant rates. However, transition rates may also depend on the history of the process and more general models, such as Weibull, allowing time-dependent rates could be more appropriate. Study subjects may also be prone to heterogeneity which can sometimes be adjusted for through the use of covariates. With unaccounted heterogeneity, frailty models would be more appropriate. Other modifications to the model, such as taking into account the sensitivity of detection of the underlying process, could also break down the simple Markovian dependence. Another extension could be to (Markov) models with more than two states. Such models obviously bring in additional complications since the transition probabilities do not allow simple expressions of Section 2.

Footnotes

Acknowledgements

This study was supported as a part of the research of the PneumoCarr Consortium funded by a grant (37875) from the Bill and Melinda Gates Foundation through the Grand Challenges in Global Health Initiative. The authors would like to thank the reviewers for their comments that helped to improve the manuscript.

Appendix

References

Hill

Cheung

Akisanya

Sankareh

Lahai

Greenwood

. Nasopharyngeal carriage of Streptococcus pneumoniae in Gambian infants: a longitudinal study. Clin Infect Dis 2008; 46: 807–814.

Trotter

Gay

. Analysis of longitudinal bacterial carriage studies accounting for sensitivity of swabbing: an application to Neisseria meningitidis. Epidemiol Infect 2003; 130: 201–205.

Nagelkerke

NJD

Chunge

Kinoti

. Estimation of parasitic infection dynamics when detectability is imperfect. Stat Med 1990; 9: 1211–1219.

Brown

Solomon

Stephens

. Estimation of parameters of zero-one processes by interval sampling. Operat Res 1977; 25: 493–505.

Kalbfleisch

Lawless

. The analysis of panel data under a Markov assumption. J Am Stat Assoc 1985; 80: 863–871.

Auranen

Arjas

Leino

Takala

AKT

. Transmission of pneumococcal carriage in families: a latent Markov process model for binary data. J Am Stat Assoc 2000; 95: 1044–1053.

Bladt

Sørensen

. Statistical inference for discretely observed Markov jump processes. J R Stat Soc B Stat Math 2005; 67: 395–410.

Albert

Brown

. The design of a panel study under an alternating poisson process assumption. Biometrics 1991; 47: 921–932.

Quintana

Müller

. Optimal sampling for repeated binary measurements. Can J Stat 2004; 32: 73–84.

10.

Cook

. Information and efficiency consideration in planning studies based on two-state Markov processes. J Stat Res 2000; 34: 161–178.

11.

Becker

Kersting

. Design problems for the pure birth process. Adv Appl Prob 1983 15: 255–273.

12.

Pagendam

Pollett

. Locally optimal designs for the simple death process. J Stat Plan Infer 2010; 140: 3096–3105.

13.

Pagendam

Pollett

. Optimal sampling and problematic likelihood in a simple population model. Environ Model Assess 2009; 14: 759–767.

14.

Hwang

Brookmeyer

. Design of panel studies for disease progression with multiple stages. Lifetime Data Anal 2003; 9: 261–274.

15.

Cook

Gavin

Gilligan

. Optimal observation times in experimental epidemic processes. Biometrics 2008; 64: 860–868.

16.

Abdelbasit

Plackett

. Experimental design for binary data. J Am Stat Assoc 1983; 78: 90–98.

17.

Tekle

Tan

FES

Berge

MPF

. Maximin D-optimal designs for binary longitudinal responses. Comput Stat Data Anal 2008; 52: 5253–5262.

18.

Karvanen

Kulathinal

Gasbarra

. Optimal designs to select individuals for genotyping conditional on observed binary or survival outcomes and non-genetic covariates. Comput Stat Data Anal 2009; 53: 1782–1793.

19.

Chernoff

. Locally optimal designs for estimating parameters. Ann Math Stat 1953; 24: 586–602.

20.

O'Brien

Millar

Zell

Bronsdon

Weatherholtz

Reid

. Effect of pneumococcal conjugate vaccine on nasopharyngeal colonization among immunized and unimmunized children in a community-randomized trial. J Infect Dis 2007; 196: 1211–1220.

21.

Smith

Adler

McAuley

Gutenkunst

Ribeiro

McCullers

. Effect of 1918 PB1-F2 expression on influenza a virus infection kinetics. PLOS Comput Biol 2011; 7: e1001081.

22.

Hope

McMillan

. Challenge studies of human volunteers: ethical issues. J Med Ethics 2004; 30: 110–116.

23.

Syrjänen

Kilpi

Kaijalainen

Herva