2.1 Response model

Let Y be the outcome variable, X be the binary exposure variable, and Z be a vector of covariates. We consider the case where these variables are linked by a generalized linear model. That is, given {X,Z}, outcome variable Y follows an exponential family distribution

where a( · ), b( · ) and c( · ) are known functions, φ is the dispersion parameter, and ξ is the canonical parameter. Here the lower case letters y,x and z stand for realizations of Y, X and Z, respectively, and the dependence on x and z of the right-hand side of (1) is suppressed in the notation as shown below.

This distribution form gives that the conditional mean and variance of Y, given {X, Z}, are, respectively, the first and second derivatives of b(ξ) with

where b′(⋅) and b′′(⋅) represent the first and second derivatives of b(·), respectively.

To explicitly reflect the influence of the covariates on the outcome variable, we further postulate the conditional mean E(Y|X,Z) via the regression model

where g( · ) is assumed to be known and monotone, and β=(β0,βx,βzT)T is the vector of regression parameters. The objective is to estimate the parameter β and to highlight this, in the following development we often explicitly show β when presenting an estimating function or a score function but suppress other parameters in the notation.

If the variables are all correctly measured, then inference about β may be based on the likelihood function:

where the conditional probability density or mass function fY|X,Z(y|X,Z) is determined by (1) and (2).

Let SY|X,Z(β;y,X,Z)=(∂/∂β)logfY|X,Z(y|X,Z) and I(β)=E{−(∂/∂βT)SY|X,Z(β;Y,X,Z)} denote the score function and the expected information matrix, respectively.

2.2 Misclassification probabilities

Suppose that covariates Z are perfectly measured, but X is subject to misclassification with a surrogate measurement X∗. Let p00=P(X∗=0|X=0,Z) be the specificity and let p11=P(X∗=1|X=1,Z) be the sensitivity. Let π˜=P(X=1|Z) be the prevalence of the exposure variable. Here for simplicity, the dependence on Z is suppressed in the symbols p00, p00, and π˜.

Then the conditional probability of X∗ given Z is given by

Applying the Bayesian Theorem gives the conditional probabilities

and

these probabilities are called the positive predictive value (PPV) and negative predictive value (NPV), respectively. For ease of notation, we let p11∗ and p00∗ denote the PPV and NPV, respectively.

In subsequent development, we assume that misclassification is nondifferential with

where fY|X,X∗,Z(y|x,x∗,z) represents the conditional probability density or mass function of Y given {X,X∗,Z}. This assumption basically says that conditional on the true covariates X and Z, the surrogate measurement X∗ is not predictive for the value of the response variable Y.

3.1 Induced likelihood method

The induced probability density or mass function of Y given {X∗,Z} is given by

where the probability density or mass function fY|X,Z(y|X=x,Z) is the same as that appears in (3), and P(X=x|X∗,Z) is determined by (4).

Let SY|X∗,Z∗(β;y,X∗,Z)=(∂/∂β)logfY|X∗,Z(y|X∗,Z) be the score function. Then the expected information matrix is given by

where the expectation is taken with respect to the joint probability density or mass function f(y,x∗,z) of Y and {X∗,Z}.

This function is determined by

where fY|X,Z(y|x,z) is given by (3), fX∗|X,Z(x∗|x,z) and fX|Z(x|z) are respectively determined by (12) and (13) in Section 3.3, and fZ(z) is the marginal probability density or mass function of Z.

The induced likelihood method was explored by various authors for different models. For instance, [32] employed this approach for the Cox proportional hazards models where continuous covariates are subject to measurement error. Spiegelman et al. [11] investigated this strategy for logistic regression with covariate misclassification and measurement error.

3.2 Subtraction correction method

It is often tempting to ignore the difference between X∗ and X, and conducting the naive analysis is unfortunately far too common in practice. That is, the naive analysis uses the score function SY|X,Z(β;y,X,Z) under the true model and replaces X with X∗ to form the estimating function for β. Although this substitution gives a computable estimating function SY|X,Z(β;y,X∗,Z), the substitution destroys the unbiasedness of the original score function SY|X,Z(β;y,X,Z). That is, the expectation of SY|X,Z(β;y,X∗,Z) is no longer zero, yielding that the resulting estimator for β is not necessarily consistent.

To obtain a consistent estimator of β, we modify function SY|X,Z(β;y,X∗,Z) by subtracting its expectation. Define

then the estimating function Usub∗(β;y,X∗,Z) is both unbiased and computable. As a result, estimation of β can be carried out using Usub∗(β;y,X∗,Z).

Given the setup in Section 2, Usub∗(β;y,X∗,Z) is specifically given by

where fY|X∗,Z(y|x∗,z) is determined by (6), and dη(y) is used to indicate an integral or a summation for all the possible values of Y, depending on whether Y is continuous or discrete.

Let Isub∗(β)=E{Usub∗(β;Y,X∗,Z)Usub∗T(β;Y,X∗,Z)} and Jsub∗(β)=E{(∂/∂βT)Usub∗(β;Y,X∗,Z)}. Then the Godambe information matrix (e.g. [3], Section 1.3) for the estimating function is given by

While the idea of the subtraction correction method is conceptually simple, this scheme has not been widely employed. As discussed by [33], the subtraction correction method requires evaluation of the expectation E{SY|X,Z(β;Y,X∗,Z)}, which is often difficult to express in an analytic form for many settings. Under the additive hazards model with continuous covariates subject to measurement error for survival data, [34] implemented this strategy and obtained a tractable estimating function form. Here, we leverage the unique properties of discrete variables and successfully implement the subtraction correction scheme to handle misclassification-prone covariates.

3.3 Expectation correction method

Define

where the expectation is taken with respect to the conditional distribution of X given {Y,X∗,Z}.

Function Uexp∗(β;Y,X∗,Z) can be invoked to estimate β because it is both computable and unbiased; the unbiasedness of Uexp∗(β;Y,X∗,Z) follows immediately from using the law of total expectation and the unbiasedness of the score function SY|X,Z(β;Y,X,Z):

where the expectations are evaluated with respect to the models for the corresponding distributions indicated by the associated random variables.

Given the model setup in Section 2, the conditional probability mass function of X∗ given {X,Z} is

and the conditional probability mass function of X given Z is

Then the conditional probability mass function of X given {Y,X∗,Z} is given by

where fY|X,Z(y|x,z) is given by (3). As a result, Uexp∗(β;y,X∗,Z) is given by

Let Iexp∗(β)=E{Uexp∗(β;Y,X∗,Z)Uexp∗T(β;Y,X∗,Z)} and Jexp∗(β)=E{(∂/∂βT)Uexp∗(β;Y,X∗,Z)}. Then the Godmabe information matrix of such an estimating function is given by

When error-prone covariates are continuous, [35] developed the expectation correction method to accommodate measurement error when replicate measurements are available. Here, we explore this strategy to account for bias due to misclassification-prone discrete covariates.

3.4 Corrected score method

In contrast to the expectation correction method, we describe the corrected score method. The idea is to find a computable estimating function, denoted by Ucor∗(β;Y,X∗,Z), so that a conditional expectation of this function can recover the score function SY|X,Z(β;Y,X,Z). That is, as long as

then using Ucor∗(β;y,X∗,Z) would produce a consistent estimator for β under certain regularity conditions.

Given the model setup of the response and misclassification processes, it is easily seen that setting

would satisfy (16). Here, S(X = k) is the value of SY|X,Z(β;Y,X,Z) evaluated at X = k for k = 0,1.

Let Icor∗(β)=E{Ucor∗(β;Y,X∗,Z)Ucor∗T(β;Y,X∗,Z)} and Jcor∗(β)=E{(∂/∂βT)Ucor∗(β;Y,X∗,Z)}. Then the Godambe information matrix of such an estimating function is given by

The corrected score method was proposed by [36] for generalized linear measurement error models with mis-measured continuous covariates. It has been successfully used for various regression models, such as Normal, Poisson and Gamma regression models ([2], Ch. 7). Extensions of this strategy to a range of settings have been discussed by several authors, including [12, 20, 37, 38]. Here, we use the corrected score scheme to develop unbiased estimating functions for estimation pertinent to regression models with misclassified discrete covariates.

3.5 Efficiency comparison

Standard estimating function theory shows that each of the preceding estimation methods results in a consistent estimator for regression parameter β, provided regularity conditions. The consistency of the estimators is ensured by the unbiasedness of the induced score function in Section 3.1 and the estimating functions in Sections 3.2–3.4. It is noted that the corrected score approach is valid under fewer modeling assumptions than the other three methods – in general, the corrected score approach does not require the distributional specification of the true covariate process but the other methods do need this. So in this sense, the corrected score approach is more robust than the other methods, as discussed by Carroll et al. [2] and Yi [3], among others.

To make a fair comparison among the four methods in Sections 3.1–3.4 in terms of the relative performance, here we assume that the parameters associated with the misclassification process are known. This assumption enables the four methods to share the same model assumptions; their differences are solely reflected from the different formulation perspectives. A general development is provided in section 4 that deals with the case where such an assumption is invalid. In Appendix A, we show the following result.

It should be noted that the equivalence of the induced likelihood and the expectation correction methods holds only when the parameters (π˜,p00,p11) are known. This relationship no longer holds when the parameters (π˜,p00,p11) are estimated from a validation study, as demonstrated in the simulations reported in Section 6.

Asymptotically, the induced likelihood method and the expectation correction method are the most efficient. In terms of relative performance of the subtraction correction and the corrected score methods, compared to the induced likelihood method, we calculate the relative efficiency using the information carried with each type of estimating functions. For each component βj of β, we define

where Isub,j(β), Icor,j(β) and Ij∗(β) represent the jth diagonal element, respectively, corresponding to Isub(β) in (10), Icor(β) in (18) and I∗(β) in (7). Then, the ratios of Rsub,j and Rcor,j can be used to compare the efficiency of the subtraction correction and the corrected score methods, relative to the induced likelihood method, for estimators of the parameter component βj.

To gain a direct understanding of the efficiency differences among the preceding estimation methods, we consider a logistic regression response model, where the outcome variable Y is binary and function g( · ) in (2) is the logit function g(t)=log{t/(1−t)}, given by

for y = 0,1. In Appendix B, we provide detailed expressions for the following efficiency comparisons.

Now, we numerically investigate efficiency differences between the preceding estimation methods under a number of scenarios with various degrees of misclassification probabilities and the prevalence of the outcome as well as different magnitudes of covariate effects. We consider the case where Z is a binary variable with probability 0.2 for assuming value 1, sensitivity p11=0.6 or 0.8, and specificity p00=0.6 or 0.8. In all these settings, we take β0=log1.2, leading to P(Y=1)≈60%, and βz=log1.2. Let βx be −log20 or log2.0.

Regarding the prevalence π˜=P(X=1|Z), we consider two cases. In the first case, we assume that X and Z are independent, and consider two values for π˜: 0.3 and 0.7. The results are displayed in Table 1. In the second case, we assume that X and Z are correlated, and the prevalence π˜ is given by

with (α0,αz)=(1,0.5), leading to the correlation between X and Z being about 15%. The results are summarized in Table 2.

It is seen that the subtraction correction method is more efficient than the corrected score method, and is almost as efficient as the induced likelihood method; the efficiency of this method does not seem to change with the misclassification degree or the values of βx we consider. For a given value of βx, the efficiency of the corrected score method on estimation of both βx and βz decreases as the degree of misclassification increases, regardless of the correlation strength between X and Z; such a misclassification effect on reducing the efficiency of the corrected score method is exacerbated when βx assumes a large absolute value.

As commented earlier, the corrected score method is more robust than other methods in that its results are not affected by the distributional specification of the covariate processes. The trade-off between the robustness and efficiency for these methods would be a key to decide which method is more suitable for a specific application.

4.1 Induced likelihood method with validation data

With a main study/internal validation study design, the subjects in the main study contribute through the conditional distribution of {Yi,Xi∗} given Zi while the subjects in the validation study contribute through the conditional distribution of {Yi,Xi,Xi∗} given Zi. Therefore, under the nondifferential misclassification mechanism, the likelihood function for the entire data is given by

where Li∗=fY,X∗|Z(yi,xi∗|zi) is proportional to (6), Li is given by (3), and Liϑ is calculated from (20).

Let SI(θ)=logLI(θ)/∂θ be the score function of the induced likelihood for the observed data. Then solving

gives an estimator of θ,θˆI, where 0 represents the zero vector of the same dimension of θ. In the following development, the symbol 0 may represent the number zero as well as a zero vector or a zero matrix whose dimension is clear from the context.

Let Siβ=∂logLi/∂β and Siϑϑ=∂logLiϑ/∂ϑ. Define Siθ∗=(Siβ∗T,Siϑ∗T)T, where Siβ∗=∂logLi∗/∂β, and Siϑ∗=∂logLi∗/∂ϑ. Under regularity conditions and when the ratio m/n approaches a positive constant ρ as n→∞, θˆI is a consistent estimator of θ, and n(θˆI−θ) has an asymptotic normal distribution with mean zero and covariance matrix ΣI given by ΣI=AI−1, where

This result can be proved by modifying standard likelihood theory; a sketch is given in Appendix C.

On the other hand, under the main study/external validation study design, the likelihood function is given by

where Li∗ is determined by (6), and Liϑ is given by (20). Then solving

gives an estimator of θ,θˆE, where the parameters are explicitly expressed in the score functions Siθ∗ and Siϑϑ.

Under regularity conditions and when the ratio m/n approaches a positive constant ρ as n→∞, θˆE is a consistent estimator of  θ, and n(θˆE−θ) has an asymptotic normal distribution with mean zero and covariance matrix ΣE given by ΣE=11+ρAE−1, where

This result can be proved by modifying standard likelihood theory; a sketch is given in Appendix D, which also shows why the asymptotic variances of θˆE and θˆI are different.

4.2 Unbiased estimating equations for main study/validation study designs

In this subsection, we explore three unbiased estimating equations methods for correcting for bias due to misclassification when validation data are available. Specifically, we consider the main study/internal validation study design and the main study/external validation study design, and develop inference procedures based on the subtraction correction, expectation correction and corrected score methods, as in Section 3.

First, we describe inference procedures using the subtraction correction method. With the main study/internal validation study design, the estimating function is given by

where Usub,i∗ is given by (9), and Siβ(β) is the score function Siβ with the dependence on the parameter β explicitly spelled out. Estimation of  θ is then obtained by solving

for  θ; let θˆsub,I denote the resulting estimator.

The formulation of (23) reflects different uses of the information from the main study and the validation study. The first term of (23) includes the contributions from the main study; its zero subvector shows that this term contributes no information on estimation of ϑ, the parameter which governs the misclassification process. The second term of (23) involves the measurements of (Yi,Xi∗,Xi,Zi) from the validation sample; this sample contributes information for the estimation of both β and ϑ, as suggested by the dependence on β of Siβ(β) and the dependence on ϑ of Siϑϑ(ϑ).

Since the main study subjects contribute no data for estimating ϑ, this estimation procedure is identical to the two-stage procedure: (1) use the likelihood score function Siϑϑ(ϑ) based on the validation data {(xi∗,xi,zi):i∈V} to obtain an estimator ϑˆsub,I of the misclassification parameters ϑsub,I, and (2) solve

to obtain the estimator βˆsub,I of βsub,I.

Under regularity conditions and that the ratio m/n approaches a positive constant ρ as n→∞, θˆsub,I is a consistent estimator of  θ, and n(θˆsub,I−θ) has an asymptotic normal distribution with mean zero and covariance matrix Σsub,I given by Σsub,I=Asub,I−1Bsub,I(Asub,I−1T), where

On the other hand, under the main study/external validation study design, the estimating function is given by

where Usub,i∗(β,ϑ) is given by (9) with the dependence on the parameters spelled out. Then solving

gives an estimator of θsub,E. Let θˆsub,E denote this estimator.

Under regularity conditions and that the ratio m/n approaches a positive constant ρ as n→∞, θˆsub,E is a consistent estimator of θsub,E, and n(θˆsub,E−θ) has an asymptotic normal distribution with mean zero and covariance matrix Σsub,E given by Σsub,E=11+ρAsub,E−1Bsub,E(Asub,E−1T), where