PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Biom J. Author manuscript; available in PMC 2017 March 7.
Published in final edited form as:
Published online 2016 December 16. doi:  10.1002/bimj.201600098
PMCID: PMC5340598
NIHMSID: NIHMS831000

Estimating the Prevalence of Atrial Fibrillation from A Three-Class Mixture Model for Repeated Diagnoses

Abstract

Atrial fibrillation (AF) is an abnormal heart rhythm characterized by rapid and irregular heart beat, with or without perceivable symptoms. In clinical practice, the electrocardiogram (ECG) is often used for diagnosis of AF. Since the AF often arrives as recurrent episodes of varying frequency and duration and only the episodes that occur at the time of ECG can be detected, the AF is often underdiagnosed when a limited number of repeated ECGs are used. In studies evaluating the efficacy of AF ablation surgery, each patient undergo multiple ECGs and the AF status at the time of ECG is recorded. The objective of this paper is to estimate the marginal proportions of patients with or without AF in a population, which are important measures of the efficacy of the treatment. The underdiagnosis problem is addressed by a three-class mixture regression model in which a patient’s probability of having no AF, paroxysmal AF, and permanent AF is modeled by auxiliary baseline covariates in a nested logistic regression. A binomial regression model is specified conditional on a subject being in the paroxysmal AF group. The model parameters are estimated by the EM algorithm. These parameters are themselves nuisance parameters for the purpose of this research, but the estimators of the marginal proportions of interest can be expressed as functions of the data and these nuisance parameters and their variances can be estimated by the sandwich method. We examine the performance of the proposed methodology in simulations and two real data applications.

Keywords: Atrial fibrillation, Latent class model, Mixture model, Two-part model, Zero-inflated binomial

1 Modeling the Atrial Fibrillation

Atrial fibrillation (AF) affects approximately 2.2 million individuals in the United States (Fuster, 2006), particularly those with structural heart disease and the elderly. It is characterized by rapid and irregular heart beat. Each AF episode may last between a few seconds to a few days, and the frequency also varies with subjects. If the AF episodes occur intermittently, they are called paroxysmal AF; if the patient is constantly in AF, this condition is called permanent AF. Most AF episodes occur without any symptoms, while some may be accompanied by perceived heart palpitations, weakness, shortness of breath, or chest pain. AF patients may suffer from tachycardia, low cardiac output from loss of atrial function, atrial and ventricular remodeling, and significantly elevated risk of stroke. Accurately detecting the start and end of the AF episodes is difficult due to the possible intermittent recurrence of the episodes. While devices such as the insertable cardiac monitors can record the AF episodes continuously over a long period of time, these devices require a surgery to be implanted on the patient, and another surgery to be extracted out at the end of the monitoring. Hence, they are not suitable for use in the general at-risk population. In the current clinical practice, physicians often rely on electrocardiogram (ECG) or Holter monitoring to capture the AF episodes. The ECG is usually performed at inpatient visits and lasts for less than 15 minutes; the Holter monitoring device may be carried at home by the patient and provides 24 to 72 hours of continuous monitoring. Since the gap times between consecutive AF episodes may vary in length between a few minutes to many days, in many situations such monitoring methods can be viewed as only providing a “snapshot” of the underlying continuous on-and-off process of AF. Consequently, the AF episodes are easily missed in these snapshots, leading to underdiagnosis or false negative diagnosis. It helps if the monitoring frequency and duration are increased (Arya et al, 2007), e.g., with the 7-day ECG or daily telephonic ECG, but these options are not always feasible in the current clinical practice.

In this paper, we investigate the diagnosis of AF under the following study design. Suppose n patients, indexed by i = 1, 2, …, n, satisfy the inclusion criteria and are enrolled in a study on the risk of AF after a surgical procedure such as the Cox-Maze, with the goal of evaluating whether the surgery reduces or eliminates the AF (Gillinov et al, 2006). After the surgery, each patient is scheduled for ni recurrent clinical visits; the presence of absence of AF at the j-th visit (j = 1, 2, …, ni) is denoted by Yij. The patients may be free of AF, have paroxysmal AF with recurrent episodes of arrythmia, or have permanent AF so that the heart is always in abnormal arrythmia. Within the time frame of the study, a patient may be in one of three states, denoted by Ci, with Ci = 0 indicating no AF, Ci = 1/2 indicating paroxysmal AF, and Ci = 1 indicating permanent AF. Ci is not directly observable, but may be inferred from the observed data Yi = (Yi1, …, yini)T. The objective of this paper is to estimate the marginal prevalence probabilities:

p0=Pr(Ci=0),p1/2=Pr(Ci=1/2),p1=Pr(Ci=1)
(1)

For a given patient population receiving AF treatments, these estimands quantifies the effectiveness of the surgery in reducing or eliminating AF in that population. An effective procedure is expected to increase the probability of no AF (p0) and decrease the probabilities of paroxysmal AF (p1/2) and permanent AF (p1). If we can estimate p0, p1/2 and p1 for a single population, extension to randomized clinical trial or propensity score matched observational studies with multiple treatment groups is straightforward because these marginal prevalence probabilities can be estimated for each group under comparison and the estimators and their estimated variances can be compared or used to form a statistical test for group differences.

Let Mi=j=1niYij denote the number of positive ECG results. A simple method to estimate the marginal prevalence is to claim that a patient is free of AF if Mi = 0, in permanent AF if Mi = ni, and in paroxysmal AF if 0 < Mi < ni. The idea of this approach conforms to some clinical practice where a patient is diagnosed with AF if the AF episodes are captured by the ECGs, and non-AF otherwise. According to this approach, the estimated marginal prevalence probabilities are:

p0=1ni=1n1{Mi=0},p1/2=1ni=1n1{0<Mi<ni},p1=1ni=1n1{Mi=ni}
(2)

These prevalence estimates may be biased. By definition, if a patient is free of AF, Yij [equivalent] 0 and Mi = 0; if a patient is in permanent AF, Yij [equivalent] 1 and Mi = ni; if a patient is in paroxysmal AF, Yij may be 0 or 1 with non-zero probabilities. By chance, some paroxysmal AF patients may have all the Yij’s being 0, which leads them to be incorrectly classified into the non-AF group. Similarly, by chance some paroxysmal AF patients may have all the Yij’s being 1, which leads them to be incorrectly classified into the permanent AF group. Mathematically, p0[p with tilde]0, p1/2[p with tilde]1/2, p1[p with tilde]1. Unless ni is very large for every subject, the equalities in these expressions are unlikely to hold even with large sample size n. This analysis explains a widely observed phenomenon in AF research that AF detection depends on the monitoring intensity: the more frequent the AF is monitored, the higher the estimated AF prevalence becomes (Arya et al, 2007; Gaillard et al 2010; Senatore et al 2006).

The estimands of interest (1) cannot be identified without additional auxiliary information. There are well established prognosis risk factors of AF. In this paper, we consider exploiting their relationship with Ci and Yi (i.e., Mi) to improve the estimation of these marginal prevalence probabilities. Let Xi be a vector of baseline prognostic factors. Suppose we can specify a hierarchical two-stage regression model for the data,

Pr(Ci=cXi;θ)andPr(MiCi=c,Xi;θ)

where c = 0, 1/2, 1 and θ is a generic notation for all the unknown parameters in this model. For the purpose of this research, θ are nuisance parameters, because the estimands of interest are p0, p1/2, p1. We first estimate θ by maximizing the log likelihood of the data {Mi, Xi; i = 1,2, …, n}:

L(θ)=i=1nlogPr(MiXi)=i=1nlog{c{0,1/2,1}Pr(MiCi=c,Xi)Pr(Ci=cXi)}
(3)

Denote the maximum likelihood estimator of θ by [theta w/ hat]. Then we propose the following estimator for the marginal prevalence probabilities p0, p1/2, p1:

p^0=Pr^(Ci=0)=1ni=1nPr(Ci=0Mi,Xi;θ^)p^1/2=Pr^(Ci=1/2)=1ni=1nPr(Ci=1/2Mi,Xi;θ^)p^1=Pr^(Ci=1)=1ni=1nPr(Ci=1Mi,Xi;θ^)
(4)

These are the sample averages of each subject’s conditional probability of Ci given the observed data, evaluated at θ = [theta w/ hat]. We can prove that these estimators are unbiased for p0, p1/2, p1, using E{Pr(Ci = c|Mi, Xi; θ)} = E{1{Ci = c}} = Pr(Ci = c) and the consistency of [theta w/ hat] for the nuisance parameter θ.

In Section 2, we provide details on the model and estimation procedure. Section 3 presents the simulation results. Section 4 includes two real data examples to illustrate the proposed method. Discussion is in Section 5.

2 Estimation

We first describe the estimation of the nuisance parameter θ through an EM algorithm, and then discuss the point and variance estimation of the estimands of interest p0, p1/2, p1. Throughout this paper, we use the following model for Pr(Ci|Xi):

Pr(Ci=0Xi)=11+exp(XiTα)Pr(Ci=1/2Xi)=Pr(Ci=1/2or1Xi)×Pr(Ci=1/2Xi,Ci=1/2or1)=exp(XiTα)1+exp(XiTα)11+exp(XiTβ)Pr(Ci=1Xi)=Pr(Ci=1/2or1Xi)×Pr(Ci=1Xi,Ci=1/2or1)=exp(XiTα)1+exp(XiTα)exp(XiTβ)1+exp(XiTβ)

This is a nested logistic regression model, in which we first model the probability of Ci = 0 vs. Ci = 1/2 or 1, and then model the probability of Ci = 1 given Ci = 1/2 or 1. For notational convenience, we assume that the first element in Xi is 1, corresponding to the intercept. With one more set of parameters, the nested logistic regression model is in general more flexible than some widely used models such as the proportional odds model, but is also less parsimonious and less transparent to interpret. Since we are interested in the marginal probabilities, and all the regression parameters θ are nuisance, we sacrifice some parsimony and interpretability to gain flexibility when making the model choice.

The conditional distribution of Mi given Ci and Xi may be specified as:

Pr(MiCi=0,Xi)=1{Mi=0}Pr(MiCi=1/2,Xi)=(niMi)piMi(1pi)niMi=(niMi)exp{MiXiTγnilog(1+eXiTγ)}pi=exp(XiTγ)1+exp(XiTγ)Pr(MiCi=1,Xi)=1{Mi=ni}
(5)

When the patient does not have AF (Ci = 0), Mi can only be 0; when the patient is in permanent AF (Ci = 1), Mi can only be ni; when the patient is in paroxysmal AF (Ci = 1/2), Mi follows a binomial distribution with probability pi, which depends on Xi. An alternative way to model (5) is to specify the conditional distribution of Yi given Ci = 1/2 and Xi through a regression model for longitudinal binary outcomes, incorporating the time effect and intrasubject correlation. Under that model specification, the estimation of marginal prevalence probabilities via (4) is still valid. However, the purpose of this article is to study the diagnosis of AF during a period of time when the patient’s AF status is stable. For that purpose, the binomial model above is more parsimonious and interpretable.

The nuisance parameters θ = (αT, βT, γT)T can be estimated by maximizing the log likelihood (3) in an EM algorithm. We define the complete data set as {Ci, Mi, Xi; i = 1, 2, …, n}, where Ci is unobserved. There are only three possibilities in the complete data set: (1) Ci = 0, Mi = 0; (2) Ci = 1/2, 0 ≤ Mini; (3) Ci = 1, Mi = ni; the case with Ci = 0, Mi > 0, or Ci = 1, Mi < ni does not exist. Hence, the complete data log likelihood is:

i=1n{1{Ci=0}logPr(Ci=0Xi)+1{Ci=1/2}logPr(Ci=1/2Xi)Pr(MiXi,Ci=1/2)}+1{Ci=1}logPr(Ci=1Xi)}=i=1n{1{Ci=0}(1)log(1+eXiTα)+1{Ci=1/2}[log(eXiTα1+eXiTα11+eXiTβ)+log(niMi)+MiXiTγnilog(1+eXiTγ)]+1{Ci=1}log(eXiTα1+eXiTαeXiTβ1+eXiTβ)}

The following is the EM algorithm:

  • (Initial Step) Obtain initial values for θ. The initial value for α and β can be obtained by fitting a nested logistic regression of 1{Mi = 0}, 1{0 < Mi < ni} and 1{Mi = ni} on Xi. The initial values for γ can be obtained by fitting a binomial regression of Mi on Xi for the subset of data with 0 < Mi < ni.
  • (E-step) Given the current value of θ(m), calculate the conditional expectation of the indicators in the complete data log-likelihood given the observed data Mi and conditional on Xi.
    Pr(Ci=cMi,Xi)=Pr(MiCi=c,Xi)Pr(Ci=cXi)c{0,1/2,1}Pr(MiCi=c',Xi)Pr(Ci=c'Xi),(c=0,1/2,1)
    Denote
    p0i(m)=Pr(Ci=0Mi,Xi;θ(m))p1/2,i(m)=Pr(Ci=1/2Mi,Xi;θ(m))p1i(m)=Pr(Ci=2Mi,Xi;θ(m))
    (6)
  • (M-step) Given p0i(m), p1/2,i(m), and p1i(m), we maximize the following expected log complete data likelihood with respect to θ:
    i=1n{p0i(m)(1)log(1+eXiTα)+p1/2,i(m)[log(eXiTα1+eXiTα11+eXiTβ)+log(niMi)+MiXiTγnilog(1+eXiTγ)]+p1i(m)log(eXiTα1+eXiTαeXiTβ1+eXiTβ)}
    Since the terms involving α, β and γ are linearly additive, we can maximize with respect to α, β and γ separately in three M-steps. Each step involves maximization of a concave function with closed form first and second derivatives, and can be completed using the Newton-Raphson method.
  • Iterate between E-Step and M-Step until the following convergence criteria is satisfied:
    arg max|θ(m)θ(m1)|<ε0
    That is, we iterate till all the elements in the vector |θ(m)θ(m−1)| be less than a pre-specified small positive number ε0. We set ε0 = 10−6 for the numerical study in this paper. If it converges at step K, We denote the final point estimator as [theta w/ hat] = θ(K).

Once the nuisance parameters [theta w/ hat] is obtained, we can estimate the marginal prevalence probabilities (1) according to (4):

p^0=1ni=1np0i(K),p^1/2=1ni=1np1/2,i(K),p^1=1ni=1np1i(K)

where the individual conditional probabilities p0i(K), p1i(K) and p2i(K) can be obtained directly from the last (the K-th) iteration of the EM algorithm.

The asymptotic variance matrix of [theta w/ hat] can be estimated by the inverse of the negative Hessian matrix, obtained by numerically differentiating the log likelihood (3). However, for the purpose of this research, θ is nuisance parameter, and the estimands of interest are p0, p1/2, and p1. We use the following partial M-estimation procedure (Stefanski and Boos, 2002) to estimate the variances of [p with hat]0, [p with hat]1/2, and [p with hat]1. The estimator [p with hat]0, and [theta w/ hat] can be viewed as solutions to the following estimating equation for κ = (p0, θT)T:

0=i=1nϕi(k)=i=1n[Pr(Ci=0|Mi,Xi;θ)p0θlog(c{0,1/2,1}Pr(Mi|Ci=c,Xi)Pr(Ci=c|Xi))]

By the sandwich method, var^(p^0) is the first diagonal entry of the matrix var^(k^)=n1A^n1B^nA^nT, with A^n=n1i=1nϕi(k)/kT|k=k^ and B^n=n1i=1nϕi(k)ϕi(k)T|k=k^. The derivatives can be calculated numerically. Variances of [p with hat]1/2 and [p with hat]1 are estimated in similar ways.

3 Simulation

We conducted simulations to study the performance of the proposed model and method. At each simulation, we generated two standard Gaussian baseline covariates with a correlation coefficient of 0.2. We consider a setting in which the three marginal prevalence probabilities, p0 = Pr(Ci = 0), p1/2 = Pr(Ci = 1/2), and p1 = Pr(Ci = 1), are close to each other and a setting in which p0 and p1 are much smaller than p1/2. Under the first setting, α = (1.2, −1.0, 1.0)T, β = (−0.2, 0.6, −0.8)T, γ = (0.3, 1.2, 1)T; under the second setting, α = (2.5, −1.0, 1.0)T, β = (−2.0, 0.6, −0.8)T, γ = (0.3, 1.2, 1)T. We varied the sample size between n = 250 and 1000, and varied the median number of repeated measures (ni) between 5 (range 1 – 9 in a uniform distribution) and 10 (range 5 – 15). Therefore, we had a total of 8 simulation scenarios. One thousand simulations were run in each setting. We compared the naive method, given by (2) and the proposed estimator (4).

The simulation results are summarized in Table 1. The naive method generally leads to biased estimator of marginal prevalence probabilities and sometimes the bias is quite large, while the proposed method is unbiased in all settings. The bias of the naive method decreases when there are more repeated measures per subject, and increasing the sample size does not help reducing bias. When there is a large number of diagnoses per subject, the chance that all these diagnosis fail to capture any AF episodes on a paroxysmal AF patient is greatly reduced, and so is the chance that these diagnoses capture at least one AF episode on the paroxysmal AF patient. The 95% confidence intervals are correct coverage probabilities, though slight under coverage can be seen in small data sets with small sample size and number of repeated measures per subject. The simulation also shows that the proposed estimators of the nuisance parameters (α, β, γ) are nearly unbiased under all the simulation scenarios (result omitted).

Table 1
Simulation results comparing the naive and the proposed methods in estimating marginal prevalence probabilities p0, p1/2, and p1. The results are reported as the average (Est) and empirical standard deviation (empSD) of the point estimators, the average ...

The proposed mixture regression model is a working model to help identify the estimands in (1), which are otherwise unidentified from the data. Hence, we conducted additional simulation to study the proposed estimators when the working model is misspecified. The simulation was designed in a similar setup as in Setting 1 of Table 1. The sample size is fixed at 500. The median number of repeated measures is 5 (range 1-9), 10 (range 5-15), or 20 (range 10-30). We fit the correctly specified three-part mixture model as well as an incorrectly specified model omitting X2. The result is visualized in Figure 1. The result reconfirms the finding in Table 1 that the correctly specified mixture model produces unbiased estimators. When the model is misspecified, there is bias in the proposed estimator, but the bias is much smaller than the naive estimator in all scenarios. When the number of repeated measures per subject increases, the bias of the naive method decreases, and so does the bias of the proposed method under misspecified model. A heuristical explanation is that when the monitoring frequency increases, it is difficult for a paroxysmal AF patient to have Yi [equivalent] 0. In other words, there is less evidence in the data for a paroxysmal AF patient to have a high probability of being in the other two groups, even when there is some model misspecification. Hence the proposed estimator converges with the naive estimator, whose bias decreases with increasing monitoring frequency. Therefore, based on this result, we recommend the use of the proposed methodology in all scenarios in replacement of the naive method, regardless of the monitoring frequency.

Figure 1
The average estimated marginal prevalence probabilities from correctly specified three-part mixture model (green, solid line), misspecified three-part mixture model (green, dashed line), and the naive method (red), when the median number of repeated measures ...

The simulations in this section were performed using functions written in R (R Core Team 2014). The R code is available for download at the publisher’s website. The initial values of the EM algorithm were chosen using the simple method described in Section 2, and convergence was achieved for all the simulation runs. Analyzing a typical simulated data set with n = 1000 took about 80 seconds on a PC with 3.30 GHz CPU and 16GB memory. The good convergence performance may be partly attributed to the concavity of the target functions in the M-step of the EM algorithm. Previous literature suggests that it is desirable to run the EM algorithm using different initial values to protect against convergence to local, instead of global maximum (Biernachi et al 2003). We intentionally let the initial values deviate from the recommended values in Section 2. The algorithm always converged to the same estimators, unless excessively large deviations were used that led to unreasonable initial values, in which case the algorithm did not converge.

4 Applications

We illustrate the proposed methodology with two real data applications. The first application is from a retrospective observational study on the effectiveness of three different surgical ablation procedures in reducing atrial fibrillation: the Cox-Maze ablation procedure, which is a cut-and-sew procedure, considered as the gold standard, but is more complicated and time consuming; the pulmonary vein isolation, which is the simplest; pulmonary vein isolation alone with added lesions, which is between the other two procedures in terms of complexity. Data on heart rhythm were collected from postoperative ECGs and extracted from the Electronic Health Records (Gillinov et al 2006). The data set includes 413 patients with at least one ECG diagnosis between 6 and 24 months. The ECG data within the first 6 months after surgery have been excluded as patients’ heart rhythms may be affected by many other peri-operative risk factors and are not stable shortly after the open heart surgery. Of these patients, 169 had one ECG, 105 has two, 49 had 3, and the maximum number of ECGs is 14. The naive method estimated that [p with hat]0 = 0.564, [p with hat]1/2 = 0.196, [p with hat]1 = 0.240. The proposed method estimated that [p with hat]0 = 0.379(0.0171), [p with hat]1/2 = 0.481(0.0208), [p with hat]1 = 0.140(0.00682) (the numbers in brackets are standard errors). The proportion of permanent AF patients appears to be overestimated under the naive method. Table 2 shows the estimated regression model parameters. These are nuisance for the purpose of estimating the marginal prevalence probabilities, but can be informative for exploring the association between risk factors and outcomes.

Table 2
The estimated model parameters from the three-part mixture model applied to the first data application. Results are presented as the point estimator, standard error estimator, and p-value from a Wald-type test against the null hypothesis of zero parameter. ...

The second data application is from a prospective randomized study to determine if adding surgical ablation to mitral valve surgery is more effective than surgery alone in reducing AF at 6 and 12 months. In this study, patients underwent weekly trans-telephonic monitoring of their heart rhythm after cardiac surgery. The trans-telephonic monitoring (TTM) system functions in a similar way as the ECG in the physician’s office, but it is smaller and more portable so that patients can use it to monitor their heart rhythm at home. The device can be connected to the home telephone line to transmit the heart rhythm data to the data center of the study. The patients were instructed to make at least one transmission each week, though more frequent transmissions were permitted. There are n = 150 patients with weekly TTM data between 6 and 12 months. The median number of repeated measures is 23, minimum 1, and maximum 46. The naive method estimated that [p with hat]0 = 0.38, [p with hat]1/2 = 0.28, [p with hat]1 = 0.34. The proposed method estimated that [p with hat]0 = 0.375(0.0168), [p with hat]1/2 = 0.287(0.0128), [p with hat]1 = 0.339(0.0152). Unlike the first example, the results from the proposed method and the naive method are very similar. As shown in the simulations, when the number of repeated measures per subject is large, the discrepancy between the proposed and the naive methods diminishes. In other words, the naive method produces nearly unbiased result with large number of repeated measures. From a clinical perspective, this finding highlights the importance and usefulness of the modern technology such as the TTM system in the diagnosis and monitoring of patients at risk of AF. Table 3 shows the estimated model parameters. As an additional sensitivity analysis, we artificially reduced the number of diagnoses per subject in this data set by randomly sampling a fraction of the repeated measures with a probability of retention [var phi], and plotted the change in estimated marginal prevalence probabilities from the naive and the proposed methods in Figure 2. The estimated probabilities from the proposed method are less sensitive to the monitoring frequency than those from the naive method.

Figure 2
The estimated prevalence probabilities by the naive (dashed lines) and proposed (solid lines) methods, as the probability of retention ω decreases, causing a reduction in the number of repeated measures per subject from the original TTM data set ...
Table 3
The estimated model parameters from the three-part mixture model applied to the second data application. Results are presented as the point estimator, standard error estimator, and p-value from a Wald-type test against the null hypothesis of zero parameter. ...

5 Discussion

In this paper, we proposed a new statistical methodology to address one aspect of the widely recognized diagnosis problem in atrial fibrillation research, that the estimated prevalence of AF increases with the monitoring intensity (Arya et al 2007; Senatore et al 2006). This phenomenon arises because AF episodes occur intermittently and may be easily missed when the monitoring frequency is not adequate. Estimating the prevalence of people with no AF, paroxysmal AF and permanent AF in a given population is of scientific interest in many contexts including the evaluation of the efficacy of AF ablation surgeries. If these probabilities are estimated based on the observed proportions such as in the naive method, the results will be biased unless a large number of repeated diagnoses are made, such as in the TTM data example. Hence, the results of this article highlight the importance of using modern, high frequency monitoring technology such as the TTM system. However, high frequent monitoring is not always available in clinical practice. When the monitoring frequency is not large, we propose to use a model-based approach. Our simulation shows that the proposed approach outperforms the naive approach and produces unbiased result regardless of the monitoring frequency. Even when the model is moderately misspecified, the proposed method can still achieve substantial bias reduction and outperforms the naive method. Our real data applications show that when the monitoring frequency is high as in the second example, the proposed method and the naive method give similar results; when the monitoring frequency is low as in the first data example, there are considerable differences between the two methods but the proposed method is better justified than the naive method. In addition, the result of the proposed method is less sensitive to the monitoring frequency than the naive method, as demonstrated in the real data experiment with the second data example (Figure 2). Based on all the results above, we recommend the proposed methodology in replacement of the naive method in all situations, regardless of the monitoring frequency.

The marginal prevalence probabilities are not identifiable with only the AF outcome Y. They are identifiable in our approach because of the use of the three part mixture model. This is analogous in a perspective to the cure model in survival analysis (Farewell, 1982), where the cured patients are not expected to have the clinical event of interest and are indistinguishable from the censored patients without additional modeling assumptions. The covariates X can be thought as auxiliary variables, and we make use of them to recover the lost information in the missing data, i.e., the group membership C of each patients in the no AF, paroxysmal AF, and permanent AF groups. Among all the patients with Y [equivalent] 0, some are expected to have higher risks than others according to their baseline variables X. This information is exploited by the three-part mixture model and is reflected in its calculation of the conditional probability of Ci given the observed data via equation (4), which leads to a correction of the bias from the naive method. There are other examples in the statistical literature that use auxiliary variables to recover information lost due to missing data. For example, Faucett et al (2002) used time-dependent covariates as auxiliary variables to impute the time to event data missing due to censoring. There is some literature on zero-inflated binomial regression, which applies to counts data with a binomial distribution with excessive zero counts (Hall 2000; Hall and Zhang 2004). While the proposed methodology in this paper uses a similar EM algorithm to find the estimators, it differs from the literature in the novel motivating application, a focus on the estimation of the marginal prevalence instead of regression coefficients, and the use of three-class mixture model instead of two-part mixture model.

Supplementary Material

R code

readme

Click here to view.(1016 bytes, txt)

simulation

Acknowledgments

This research is funded by NIH grants R01HL103552 and P30CA016672. We want to thank the two referees and the Associate Editor for constructive comments that improved this paper.

Footnotes

Conflict of interest None.

Supporting Information

“code.zip” in the Supporting Information contains the code required to generate the simulation results in this paper together with instructions for their use. Available at the publisher’s website.

References

  • Arya A, Piorkowski C, Sommer P, Kottkamp H, Hindricks G. Clinical implications of various follow up strategies after catheter ablation of atrial fibrillation. Pacing Clin Electrophysiol. 2007;30(4):458–62. [PubMed]
  • Biernachi C, Celeux G, Govaert G. Choosing starting values for the EM algorithm for getting the highest likelihood in multivariate Gaussian mixture models. Computational Statistics and Data Analysis. 2003;41:561–575.
  • Farewell VT. The use of mixture models for the analysis of survival data with long-term survivors. Biometrics. 1982;38:1041–1046. [PubMed]
  • Faucett CL, Schenker N, Taylor JM. Survival analysis using auxiliary variables via multiple imputation, with application to AIDS clinical trial data. Biometrics. 2002;58(1):37–47. [PubMed]
  • Fuster V. ACC/AHA/ESC 2006 guidelines for the management of patients with atrial fibrillation. Circulation. 2006;114(7):e257–354. [PubMed]
  • Gaillard N, Deltour S, Vilotijevic B, Hornych A, Crozier S, Leger A, Frank R, Samson Y. Detection of paroxysmal atrial fibrillation with transtelephonic EKG in TIA or stroke patients. Neurology. 2010;74(21):1666–70. [PubMed]
  • Gillinov AM, Bhavani S, Blackstone EH, Rajeswaran J, Svensson LG, Navia JL, Pettersson BG, Sabik JF, 3rd, Smedira NG, Mihaljevic T, McCarthy PM, Shewchik J, Natale A. Surgery for permanent atrial fibrillation: impact of patient factors and lesion set. Annals of Thoracic Surgery. 2006;82(2):502–13. [PubMed]
  • Hall DB. Zero-inflated Poisson and binomial regression with random effects: a case study. Biometrics. 2000;56(4):1030–1039. [PubMed]
  • Hall DB, Zhang Z. Marginal models for zero inflated clustered data. Statistical Modeling. 2004;4:161–180.
  • R Core Team. R: A language and environment for statistical computing. R Foundation for Statistical Computing; Vienna, Austria: 2014. URL http://www.R-project.org/
  • Senatore G, Stabile G, Bertaglia E, Donnici G, De Simone A, Zoppo F, Turco P, Pascotto P, Fazzari M. Role of transtelephonic electrocardiographic monitoring in detecting short-term arrhythmia recurrences after radiofrequency ablation in patients with atrial fibrillation. J Am Coll Cardiol. 2006;45(6):873–6. [PubMed]
  • Stefanski LA, Boos DD. The calculus of M-estimation. The American Statistician. 2002;56:29–38.