PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptNIH Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
J Stat Plan Inference. Author manuscript; available in PMC Feb 1, 2010.
Published in final edited form as:
J Stat Plan Inference. Feb 1, 2009; 139(2): 454.
doi:  10.1016/j.jspi.2008.04.030
PMCID: PMC2719909
NIHMSID: NIHMS80163

Partial covariate adjusted regression

Abstract

Covariate adjusted regression (CAR) is a recently proposed adjustment method for regression analysis where both the response and predictors are not directly observed (Şentürk and Müller, 2005). The available data has been distorted by unknown functions of an observable confounding covariate. CAR provides consistent estimators for the coefficients of the regression between the variables of interest, adjusted for the confounder. We develop a broader class of partial covariate adjusted regression (PCAR) models to accommodate both distorted and undistorted (adjusted/unadjusted) predictors. The PCAR model allows for unadjusted predictors, such as age, gender and demographic variables, which are common in the analysis of biomedical and epidemiological data. The available estimation and inference procedures for CAR are shown to be invalid for the proposed PCAR model. We propose new estimators and develop new inference tools for the more general PCAR setting. In particular, we establish the asymptotic normality of the proposed estimators and propose consistent estimators of their asymptotic variances. Finite sample properties of the proposed estimators are investigated using simulation studies and the method is also illustrated with a Pima Indians diabetes data set.

Keywords: Asymptotic normality, Binning, Confidence intervals, Multiple regression, Varying-coefficient models

1 Introduction

Covariate adjusted regression has been recently proposed to adjust for the distorting effects of a confounder in a regression setting. It was motivated by a common adjustment method in medical and health related studies. The adjustment entails normalization by anthropometric measurements, such as body mass index (BMI) and/or other measures of body configuration, as confounding variables that affect the primary variables of interest. For example, in a study involving haemodialysis patients, it is of interest to examine the relationship between elevated plasma fibrinogen level (a risk factor for cardiovascular disease in haemodialysis patients) and other predictors, such as serum transferrin protein level (Kaysen et al., 2003; Şentürk and Müller, 2005). However, both primary variables, fibrinogen and transferrin protein levels, are known to depend on body mass index, which exerts a confounding effect on the protein measurements. A common approach to adjust for the confounders, like BMI, is to normalize the primary variables of interest by simply dividing (by BMI). Note that this adjustment by division implies that the assumed contamination is of a multiplicative form. Let , X, and U denote the observed fibrinogen concentration, serum transferrin level, and confounder BMI, respectively. Using these notations, the adjusted primary variables that are thought to be free from the confounding effect of BMI are,

equation M1

The basic motivation for the above adjustment is to obtain normalized versions of the observed primary variables by removing the confounder effects, so that the measurements are comparable across patients. Other examples include normalizations by BMI in studies on diabetes, and division of brain volumetric structures by total brain volume in neurological studies (Pinter et al., 2001).

Şentürk and Müller (2005, 2006) proposed a more flexible adjustment, by modeling the confounding through unknown functions of the confounder instead of the confounder itself. This reflects the uncertainty encountered in many applications about the precise nature of the commonly assumed multiplicative relation between the confounder and the variables. For the case of p predictors, Şentürk and Müller model the underlying variables as

equation M2

where they are defined to be the parts of the observed variables, , X1,…, Xp, that are independent of the observable confounder U. In the haemodialysis data example, the latent variables would be defined to be serum protein levels adjusted for body mass index. Here, ϕ1(·),…, ϕp(·) and ψ(·) denote unknown smooth contaminating functions of U. CAR gives consistent estimators of the coefficients in the unobserved regression model which can be expressed as

equation M3

where e is the error term, assumed to be independent of equation M4 and U. The estimation procedure is based on the observed data: the distorted response, , distorted predictors, equation M5, and confounder U.

The main goal of this paper is to construct estimation and inference procedures needed to allow some of the variables as unadjusted/undistorted predictors, denoted by Z1,…, Zs. The proposed underlying regression model is then of the form

equation M6
(1)

where Xr = ϕ(U)Xr and = ψ(U)Y denote the adjusted/distorted predictors and response, respectively, and Zs denotes the unadjusted predictors. The observed data is , Xr, Zs and U. Furthermore, the confounding covariate, U, is allowed to depend on the unadjusted predictors {Zs}. The flexibility of allowing both distorted and undistorted predictors is needed, particularly in the regression analysis of biomedical data. In many of these applications, which motivated our current work, the researcher is directly interested in the effects of Zs = age, gender, obesity measures and/or ethnicity on Y. Therefore, the predictors Zs are unadjusted/undistorted. A specific example is the study of Kaysen et al. (2002) where albumin turnover and protein catabolic rate were adjusted for body surface area (U) via division, while age and gender were among the unadjusted variables. This is an example of model (1) above, where an adjustment method is needed that allows for unadjusted predictors. In this example, model (1) is used to reflect (a) the known dependence of albumin turnover and protein catabolic rate on U and (b) the specific interest in the direct effects of age and gender in the regression relationship. This issue is further discussed in Section 3 in the context of estimation as well as in Section 8.

Under the more general partial covariate adjusted regression (PCAR) setting, formally presented in Section 2, the original CAR estimators (Şentürk and Müller, 2005) for equation M7 are inconsistent and may have an arbitrarily large asymptotic bias as shown in Section 3. We propose alternative estimators that are consistent under this extended CAR setting where the issues of estimation are discussed in Section 3. The proposed PCAR methodology, like CAR, provides consistent estimators not only under multiplicative but also under additive distortion as discussed also in Section 3. The inference procedures developed for the CAR modeling are not valid for the PCAR setting, mainly due to the different dependence structure needed for PCAR. This new structure is explained in detail in Section 2 and Section 3. Thus, we develop new theoretical tools for valid inference in the PCAR model. We derive the asymptotic distributions of the proposed estimators, and present them in Section 4. Consistent estimators of the asymptotic variance are also derived in Section 4. Simulation studies to characterize the finite sample properties of the proposed estimators are summarized in Section 6. The method is further illustrated with a Pima Indians diabetes data set given in Section 5. The proofs of the main results are assembled in Section 7, where some technical conditions and auxiliary results are deferred to the Appendix. We conclude with a brief discussion in Section 8.

We note here that an advantage of CAR and PCAR is that, under the identifiability conditions introduced in Section 2, it yields consistent estimates whether the distortion is multiplicative or additive, i.e. Y = − ψ(U) and Xr = Xr − ϕr(U). (A more detailed discussion of the additive distortion case is given in Section 3.) Additive distortions can be handled by the method of nonparametric partial regression; however, there existed no consistent estimation procedure targeting the γ’s under multiplicative distortion of both the response and the predictors. Also, the proposed distortion setting has similarities with measurement error modeling if one views ψ(U) and ϕr(U) as unobserved errors affecting the response and predictors. There is an extensive literature on additive measurement error modeling, which dates back to Berkson (1950). See Carroll et al. (2006) and references therein for a comprehensive overview. A key difference between the proposed framework and traditional measurement error is that, in the proposed distortion setting, the error is a function of an observable covariate U. Also, since U is observable this information is incorporated into the estimation. While much work is devoted to additive measurement error, work on multiplicative measurement errors is limited. Relevant work includes Hwang (1986) and Iturria, Carroll and Firth (1999) who proposed estimation procedures targeting the regression coefficients under multiplicative measurement error in the predictors. However, the case of multiplicative measurement errors that affect both the predictors and the response has not been considered previously to our knowledge. Therefore, the method and theory presented here may potentially be of interest to the area of measurement error modeling, despite the differences in the methodological motivation.

2 Partial covariate adjusted regression models

We consider the underlying (unobserved) regression model

equation M8
(2)

where Yni, eni, χni = (1, Xn1i,…, Xnpi, Zn1i,…, Znqi)T and α = (γ0,…, γp, δ1,…, δq)T are the response, error, p + q predictors and unknown regression coefficients, respectively. The error variable e has mean zero, variance σ2 and a finite moment that is higher than the 4th moment. The goal is estimation and inference for the parameter vector α of the unobserved regression model (2). Estimation is based on available distorted predictor and response data, namely equation M9, where

equation M10
(3)

The unknown distorting functions equation M11 are assumed to be smooth functions of the confounder, U.

Some constraints on the unknown smooth distortion functions are needed for the identifiability of the estimation problem. A set of reasonable constraints for ψ(·) and {ϕr(·)} is implied by the natural assumption that the mean distorting effect should correspond to no distortion (Şentürk and Müller, 2005), i.e. the means of adjusted variables are the same as the means of the observed variables, E(Xr) = E(Xr) and E() = E(Y). These conditions directly imply that

equation M12
(4)

We consider the following dependence structure. The underlying predictors Xr and the undistorted predictors Zs are allowed to be dependent. The error, e, is assumed to be mutually independent of Xr, Zs, and U. We depart from the original CAR model (Şentürk and Müller, 2005), where U is independent of all the latent predictors, by allowing U to depend on Zs, while still being independent of Xr. We will elaborate further on this important difference of the proposed setting from CAR at the end of Section 3. This is an important flexibility of the proposed method, since the common confounder correlates with all the observed variables in these distortion settings. This is consistent with the assumption that the observed predictors Xr and Zs are dependent on the confounder U, and that the latent variable Xr is defined to be the part of Xr that is independent of U.

The assumption that the underlying predictors, equation M13, and response, Y, are independent of the contaminating variable U is a fundamental assumption for the estimation procedure. It defines the proposed contamination setting through defining the unobserved, underlying variables. This independence assumption cannot be checked in practice since Xr and Y are unobservable. Instead, the question of more relevance in practice is whether the independence conditions help define interpretable latent variables of interest from their observable counterparts. In the haemodialysis data example, the latent variables are defined to be serum protein levels adjusted for body mass index, which are commonly used in medical studies.

We refer to the model described by (2)(4) as the partial covariate adjusted regression (PCAR) model, since only a partial set of the predictors are adjusted for the confounder. Note also that the CAR model is a special case of the PCAR model.

For the estimation, note that it follows from (2) and the mutual independence of {e and U}, {e and (Xr,Zs)}, and {U and Xr}, for r = 1,…, p, s = 1,…, q, that the regression of on [chi] = (1, X1,…, Xp,Z1,…, Zq)T leads to a fully observable varying coefficient model (Cleveland, Grosse and Shyu, 1991; Hastie and Tibshirani, 1993),

equation M14
(5)

where

equation M15
(6)

and ϵ(Uni) = ψ(Uni)eni. For an extensive overview of the estimation procedures proposed for varying coefficient models, see Wu and Yu (2002). Note that in (6), the varying coefficient functions, equation M16, are proportional to the quotient of the original distorting functions, {ψ(·)/ϕr(·)}; both the intercept function, β0(·), and the functions equation M17 are proportional to ψ(·). The constants of proportionality are precisely the underlying regression parameters, {γr, δs}, of interest. These connections allow estimation of the underlying model through the varying coefficient functions. In Section 3 below, we describe an estimation procedure which targets {γr, δs} and mitigates the effects of the distorting functions {ψ(·), ϕr(·)}.

3 Estimation procedure

The estimation of the regression coefficients, γ0, equation M18, in the underlying regression model equation M19 is a two-step procedure. The first step involves estimation of the varying coefficient functions in model (5), namely β0(·), equation M20 using a binning approach. These varying functions are estimable because , Xr, Zs, and U are all observable. The underlying regression coefficients are targeted in the second step, with weighted averages of the estimated β0(·), βr(·) and ηs(·) for γ0, γr and δs, respectively. The estimation makes use of the relations between the varying coefficient functions and the regression coefficients given by (6) and the identifiability conditions (4), as will be described next.

The binning approach for the estimation of the varying coefficient functions involves dividing the support of U into m equidistant bins and then fitting linear regressions of on [chi] using the data falling within each bin. The observed data is the collection of n samples: equation M21. It is assumed that the confounding covariate, U, is bounded below and above, aUb, where a < b are real numbers. In practice a and b would be taken to be mini Uni and maxi Uni, respectively. The estimation procedure initially divides the interval [a, b] into m equidistant intervals, denoted Bn1,…, Bnm and referred to as bins. Let Lnj be the number of Uni’s falling into bin j. Furthermore, let equation M22 be the kth data element in the jth bin, Bnj, where k = 1,…, Lnj. Data elements in any given bin are marked by a prime.

After the initial binning of the data, a linear regression is fitted to the data observed within each bin Bnj, j = 1,…, m. The least squares estimator of the multiple regression of the data in the jth bin is

equation M23
(7)

where the response vector is equation M24 is the Lnj × (p+q+1) data matrix in bin j, with the kth observation equation M25. The estimated regression coefficients in each bin (7) are the estimators of the varying coefficient functions.

In the second step of the estimation procedure, the estimators of the targeted regression parameters, γ0, equation M26, are obtained as weighted averages of the estimators equation M27 from the m bins. The proposed PCAR estimators for γ0, equation M28 are

equation M29
(8)

where equation M30. The weights in (8) depend on the number of data points in each bin, namely Lnj for j = 1,…, m. Note that the estimators proposed in (8) are method of moments estimators targeting E0(U)}, E{βr(U)Xr}/E(Xr) and E{ηs(U)}, respectively. It follows then that they are consistent for the underlying parameters, formally stated in Section 4, since E0(U)} = γ0, Er(U)Xr}/E(Xr) = γr and E{ηs(U)} = δs, by the relations in (6) and the identifiability conditions.

We note that the estimators [gamma with circumflex]n0 and [gamma with circumflex]nr have the same form as the CAR estimators (Şentürk and Müller, 2005), whereas [delta with circumflex]ns are different. Furthermore, a straightforward application of the CAR algorithm yield inconsistent estimators for δs under the more general PCAR model. To see this, denote the original CAR estimators for δs by equation M31. It follows from Şentürk and Müller (2005), that

equation M32

where equation M33. The estimators equation M34 do not target δs, instead they target E{ηs(U)Zs}/E(Zs) = δsE{ψ(U)Zs}/E(Zs) = δsCs, where Cs [equivalent] E{ψ(U)Zs}/E(Zs) = [cov{ψ(U),Zs}/E(Zs)]+1 can get arbitrary large as E(Zs) approaches zero. Note here that since Zs is correlated with U, one can argue that a latent variable can also be created for Zs that would be independent of U, similar to the construction of Xr and Y. An argument can be made that the dependence structure suggests equation M35, for some latent variable equation M36 that is independent of U, which can be ed as equation M37. In this case, Cs further simplifies to equal equation M38 (This insight was provided by a reviewer.) From this last form of Cs, it follows immediately that if (1) the response is not distorted by U, i.e. ψ(U) = 1, or if (2) Zs is independent of U, i.e. equation M39, in which case the CAR estimator for δs would be consistent. For the latter mentioned case (2) of assuming Zs is independent of U, this is not realistic in most data applications, since U is the common confounder that correlates with all the observed variables in these distortion settings. In the data examples given in the Introduction, the undistorted predictor age is not independent of the confounders like body mass index, body weight or height.

For the former mentioned case (1) when the response is undistorted, simpler adjustments given by Hwang (1986) and Iturria et al. (1999), for multiplicative measurement error only in the predictors, would also be applicable. Hwang (1986) proposes a consistent estimator for the regression coefficients by estimating and adjusting for the bias of the regular least squares estimator. The estimation assumes that consistent estimates of the moments of the measurement error are available. Iturria (1999) proposes two estimation methods, where the first considers specific distributional forms for the measurement error and the second also models the distribution of the unobserved predictors. The two approaches of Hwang and Iturria are similar to PCAR in assuming that the error is independent of the unobserved predictors and that it has a mean of one. The difference of PCAR from the previous two approaches is that no knowledge of the distributional forms or the specific moments are assumed. Instead, information from the observed covariate U, of which the measurement error ϕr(U) is a function of, is utilized in the proposed estimation procedure.

The CAR estimators are biased for only the case of undistorted predictors. In other words, if all variables are considered as distorted then CAR provides consistent estimators. This is the difference between CAR and PCAR. While the PCAR set-up allows researchers the choice to consider a subset of predictors as distorted, the CAR set-up requires that all predictors are considered as distorted in order to yield consistent estimates.

As the PCAR set-up allows researchers to consider a subset of predictors as distorted, an important issue is how to determine whether a predictor should be considered as distorted or undistorted. Note that this distinction cannot be made using a statistical/analytical approach via studying the dependence or the relations between the observed variables and U. This is because all observed variables are dependent on U in the model (correlate with U) whether they are distorted or undistorted. Hence, this decision instead should be made by considering the duality between (a) the decision/assumption on undistorted (distorted) variables and (b) the specification of the underlying model; more specifically, the consideration of the predictor choice in the underlying model.

More precisely, determining to include a predictor as distorted or undistorted corresponds to two completely different underlying models of interest, one involving Zs and the other involving equation M40, respectively. Specification of the underlying model will depend on the specific interest of the researcher and the specific context of the application. For example, in the data analysis presented in Section 5, the goal is to uncover the relation between a diabetes marker and diastolic blood pressure adjusted for body mass index. Hence, the diabetes marker and diastolic blood pressure are considered as adjusted/distorted variables. On the other hand, age and triceps skin fold thickness are considered unadjusted/undistorted, as we are interested in the direct effects of age and triceps skin fold thickness on the BMI-adjusted diabetes markers.

We also note here that the PCAR estimators given in (8) are consistent for the parameters of the underlying model (2) also under additive distortion. More precisely, consider the simple case of one distorted and one undistorted predictor. The regression model in (2) simplifies to Y = γ01X1Z+e, where the multiplicative error structure is replaced by the additive error structure, given by = Y + ψa(U) and X = X + ϕa(U). The proposed PCAR estimators given in (8) are consistent even when the multiplicative error is replaced by additive error as described above. The above additive error model leads to the following specific varying coefficient model: E(Ỹ|X,Z,U) = β0(U) + β1(U)X + η1(U)Z, where β0(U) = γ0 − γ1ϕa(U) + ψa(U), β1(U) = γ1 and η1(U) = δ1. The PCAR estimators given in (8), namely [gamma with circumflex]0, [gamma with circumflex]1 and [delta with circumflex]1, target

equation M41

respectively. This holds regardless of the specific error structure, whether it be additive or multiplicative. Furthermore, under the additive distortion model, we have that

equation M42

This follows since Ea(U)} = Ea(U)} = 0 in the additive distortion model, under the identifiability condition of no average distortion, i.e. E() = E(Y) and E(X) = E(X). Thus, the PCAR estimators proposed in (8) are consistent for parameters of the underlying model also under additive distortion structure.

4 Asymptotic properties

We present the asymptotic distribution of the estimators [gamma with circumflex]n0, [gamma with circumflex]nr and [delta with circumflex]ns in (8) when the number of subjects n tends to infinity. As in typical smoothing applications, the number of bins m = m(n) is required to satisfy m → ∞, n/(m log n) → ∞ and equation M43 as n → ∞. We denote convergence in distribution by equation M44 and convergence in probability by equation M45. Also let

equation M46
(9)

denote the least squares estimators of the multiple regression of the unobserved data falling into Bnj, where the vectors equation M47 are defined the same way as equation M48, with equation M49 replacing equation M50, respectively. This quantity is not estimable, but will be used in the proof of the main results.

For the PCAR estimators given in (8) to be well defined, the least squares estimators given in (7) must exist for each bin Bnj, i.e. equation M51. Correspondingly, the estimators in (9) will exist under the condition that equation M52. The following theorems are given under event En, explicitly defined in the Appendix, which summarizes the above outlined conditions for the PCAR estimators to be well defined.

For the following theorems, we define the following notations: λψ = E2(U)}, λϕ = E2(U)}, λψϕr = E{ψ(Ur(U)}, equation M53, χT = (1, X1,…, Xp,Z1,…, Zq), Γ = E(χχT|U), x2133 = Γ−1χχTΓ−1, equation M54 for l = 1,…, p + s + 1 and k = 1,…, Lnj.

Theorem 1

Under the technical conditions (C1)–(C7) in Section 6, on event En with pr(En) → 1 as n → ∞,

equation M55
equation M56

where

equation M57

Theorem 1 establishes the asymptotic normality of the proposed PCAR estimators. The following theorem provides consistent estimators of the asymptotic variances given in Theorem 1.

Theorem 2

Under the technical conditions (C1)–(C7) in Section 6, on event En with pr(En) → 1 as n → ∞,

equation M58

where

equation M59
equation M60

Normalizing by the above consistent variance estimators, it holds that

equation M61

Therefore, the approximate (1−α)100% asymptotic confidence intervals for γr and δs have the endpoints

equation M62
(10)

where zα/2 is the (1 − α/2)th quantile of the standard Gaussian distribution.

Remark

These proposed variance estimators are motivated by the identifiability conditions, the definition of the smooth varying coefficients functions given in (6), Lemma 3 and Lemma 4 (a.). Using the consistency of [beta]nrj and [eta w/ hat]nsj for the values of the functions βr and ηs at the midpoint of the jth bin and the definitions of equation M63, we target the quantities equation M64 with the estimators equation M65, respectively. Furthermore, relying mainly on Lemma 3 and Lemma 4 (a.), we target equation M66 and equation M67, respectively.

5 Application to the Pima Indians diabetes data

We illustrate the proposed partial covariate adjusted regression methodology with an application to the Pima Indians diabetes data set, available at http://www.ics.uci.edu/~mlearn. Obesity is an important contributing factor to diabetes and has been widely studied in the Pima Indians population (Smith et al., 1988; Knowler et al., 1991; Hansen et al., 1998). One-half of adult Pima Indians have diabetes and 95% of those with diabetes are overweight (National Institute of Diabetes and Digestive and Kidney Diseases, http://diabetes.niddk.nih.gov). The available data comes from a larger database, where the subgroup used consists of n = 524 females at least 21 years old and of Pima Indian heritage. (The population lives near Phoenix, Arizona, U.S.A.) An oral glucose tolerance test is one of the diagnostic tests for type II diabetes. The goal is to uncover the underlying, BMI adjusted, regression relation, PGC = γ0 + γ1DBP + δ1Age + δ2TSFT + e, based on the observed plasma glucose concentration (equation M68; from a oral glucose tolerance test), diastolic blood pressure (equation M69), triceps skin fold thickness (TSFT), age and body mass index. We chose to adjust only the main relation of interest, namely the one between plasma glucose concentration (the response) and diastolic blood pressure for body mass index, and included age and triceps skin fold thickness as unadjusted predictors as they are commonly accounted factors in studies on diabetes.

Table 1 gives the regression coefficient estimates for (γ0, γ1, δ1, δ2) using the proposed PCAR method, CAR method, the ordinary least squares (OLS) estimates from regressing the observed equation M70 on (equation M71, Age, TSFT) without adjusting for the confounder BMI, and adjustment via division, i.e. regressing equation M72 on (equation M73, Age, TSFT). The approximate 95% asymptotic confidence intervals for the regression parameters obtained through all three methods are also displayed. The approximate confidence intervals for PCAR estimates were obtained as proposed in (10).

Table 1
Parameter estimates for the regression model PGC = γ0 + γ1DBP + δ1Age + δ2TSFT + e, obtained by least squares regression of = equation M169 (plasma glucose concentration) on = equation M170 (diastolic blood pressure), Z1 = Age and Z2 = ...

The implementation of the binning algorithm allows for merging of sparsely populated bins. Bin widths were chosen such that there are at least (p + q + 1) points, enough to fit the linear regression with (p + q) predictors in each bin. If there were bins with less than (p + q + 1) elements, such bins were randomly merged with neighboring bins. The merging algorithm is randomized to avoid the introduction of any additional bias. It starts by merging bins with no points. If there are more than one such bin, it randomly picks one and merges it with its neighbor of smallest number of points. After merging all the bins with no points, the bins with one point and eventually bins with p + q points are merged. For this example with n = 524 (after the removal of outliers), the average number of points per bin was 15, yielding a total of 34 bins after merging. Note that CAR estimates have been shown to be sufficiently robust regarding different choices of m, under the rate conditions given in Section 4 (Şentürk and Müller, 2006). We have found this property to hold for the proposed PCAR estimates as well, where the range of m values yielding robust estimates for different sample sizes is given explicitly in the next section.

Note that coefficients obtained by adjustment via division are quite different from the other three methods applied. In this adjustment the coefficient of DBP becomes quite pronounced compared to the other two predictors. This is most likely due to the pseudo dependence created between equation M74 and equation M75 via division by the common variable BMI. This is an example of the misleading conclusions that adjustment by division may suggest. In other words, if the original contamination is not exactly multiplication by the confounder (BMI in this example), then normalization by division may create further confounding, or “coupling” (as defined in Archie, 1981), creating a pseudo dependence that does not exist in the original data.

Even though OLS estimates for blood pressure and age are different from the PCAR and CAR estimates, all are found statistically significant at the usual 5% level. Thus, diastolic blood pressure and age are still important predictors of PGC even after adjusted for body mass index. However, using OLS, TSFT is a significant predictor of PGC, but it is not significant using PCAR and CAR at the 5% significance level. This result is not too surprising, since both TSFT and body mass index are indicators of obesity. They are positively correlated (Pearson correlation 0.67). Thus, adjusting for one, the other becomes an insignificant factor for predicting plasma glucose concentration. We note that even though estimation via CAR leads to the same conclusion as PCAR on the significance of the predictors for this analysis, the estimates from these two methods are different for TSFT. This is again to be expected, since CAR estimates are shown to be biased for the undistorted predictors.

6 Numerical studies

To examine the numerical properties of the estimators, we implemented the following simulation studies. The underlying multiple regression model is

equation M76
(11)

where the parameters of interest are (γ0, γ1, γ2, δ)T = (4,−1, 0.3, 3). The error variable is e ~ N(0, .5), and the confounder variable U is generated from a uniform distribution on [2, 6]. We considered the joint distribution of the predictors to be multivariate normal: (X1,X2,Z)T ~ N3(µ,Σ), with a general covariance structure

equation M77

The mean vector is µ = (0.7, 1.2, |U| − 3.5)T, so that the undistorted predictor Z is dependent on U. To simulate the distorted (observed) data, we consider the following distorting functions, ψ(U) = (U + 3)/7, ϕ1(U) = (U + 1)2/26.3333, and ϕ2(U) = (U + 10)/14, satisfying the identifiability constraints that E{ψ(U)} = 1 and Er(U)} = 1. The distorted response and predictors are = ψ(U)Y, X1 = ϕ1(U)X1, and X2 = ϕ2(U)X2. Under this simulation setting, we examine (1) the confidence interval coverage levels based on the asymptotic results and (2) the finite sample bias of the estimators, as well as comparing CAR and PCAR estimators in terms of variance and MSE.

We conducted 1000 Monte Carlo simulation runs for sample sizes n = 100, 150, 350, 800, and 1400 to study the approximate asymptotic confidence intervals given in (10). For the sample sizes n = 100, 150, 350, 800 and 1400, the total number of bins formed were m = 16, 27, 32, 50 and 70. Table 2 summarizes the coverage and interval lengths, averaged over the 1000 simulation runs, for the approximate 95% asymptotic confidence intervals for the parameter vector (γ0, γ1, γ2, δ)T = (4,−1,.3, 3). The numerical study indicates that the estimated non-coverage percentages are close to the target value of 0.05, as the sample size n increases. The estimated interval lengths are decreasing as n increases, as expected.

Table 2
Estimated coverage (in percent) and average interval lengths for the approximate 95% confidence intervals formed for the parameters of the regression model (11), corresponding to sample sizes n = 100, 150, 350, 800, and 1400. Results are based on 1000 ...

We also examined the bias, variance and mean squared error (MSE) of the proposed estimators in comparison to the CAR estimators. For example, the estimated (absolute bias, variance, MSE) values for PCAR estimators at the smallest sample size n = 100 are (0.0112, 0.2223, 0.2224), (0.0120, 0.1392, 0.1393), (0.0028, 0.0421, 0.0421) and (0.0167, 0.0651, 0.0654) for [gamma with circumflex]0, [gamma with circumflex]1, [gamma with circumflex]2 and [delta with circumflex], respectively. These values are averages over 1000 Monte Carlo runs. The results are similar for other sample sizes, where the variance seems to be the dominating factor contributing to the MSE. The estimated (bias, variance, MSE) values for the CAR estimator for δ (even though their asymptotic distributions are different, the three other point estimates for γ0, γ1 and γ2 are the same for the two methods), [delta with circumflex]*, are (1.294, 1.694, 3.368) at the same sample size of n = 100. The multiplicative bias factor of the CAR estimate for δ, shown to be Cs = E{ψ(U)Zs}/E(Zs) in Section 3, is equal to 1.416 for this simulation set-up. As expected, the CAR estimate [delta with circumflex]* is off target for δ = 3 (with a mean of 4.294 at n = 100, and 4.146 at n = 1400). In addition to being biased, note that the CAR estimator [delta with circumflex]* has substantially larger variance relative to the PCAR estimator.

As stated above, for sample sizes n = 100, 150, 350, 800 and 1400, the total number of bins formed were m = 16, 27, 32, 50 and 70, respectively. We carried out additional simulation studies to examine the affect of the number of bins m. The results suggest that the estimators are robust, based on estimated MSE, when m is chosen in the intervals [13, 18], [18, 27], [25, 45], [35, 65] and [60, 90] corresponding to sample sizes n = 100, 150, 350, 800 and 1400. In a given application, the above intervals can give rough guidelines on how the choice of m may change with sample size, although a sensitivity analysis for the choice of m specific to the data would also be informative.

Finally, we note that even though the smallest sample size at which the proposed asymptotic confidence intervals attain reasonable coverage in our simulation study is n = 100 (the coverage is around 5% off the targeted level for n = 100), the proposed PCAR point estimators still yield reasonable (bias, variance and MSE) values: (0.0515, 0.5748, 0.5775), (0.0026, 0.3552, 0.3553), (0.0256, 0.0916, 0.0922) and (0.0205, 0.3328, 0.3332) for γ0 = 4, γ1 = −1, γ2 = 0.3, δ = 3, respectively, at n = 50. These values are estimated under the current simulation set-up with 3 predictor variables for n = 50 and m = 8. For the case of simple linear regression, roughly the same number of bins can be attained with n = 30, and hence n = 30 would be the smallest sample size where CAR and PCAR give reasonable point estimates. For sample sizes smaller than 30, systematic localization via binning may not be fully feasible so the PCAR estimates should be taken with caution. In this case a very rough localization by stratification (by U) into 2–3 groups for crude comparisons is possible.

7 Proofs of the main results

We provide the major steps of the proofs of the main results (Theorem 1 and 2) here and defer the auxiliary results for these proofs to the Appendix, where they are listed as lemma 1 to lemma 4. We introduce the following technical conditions:

  • (C1) The covariate U is bounded below and above, −∞ < aUb < ∞ for real numbers a < b. The density f(u) of U satisfies infa≤u≤b f(u) > c1 > 0, supa≤u≤b f(u) < c2 < ∞ for real c1, c2, and is uniformly Lipschitz continuous, i.e., there exists a real number M such that supa≤u≤b|f(u + c) − f(u)| ≤ M|c| for any real number c.
  • (C2) The variables (e, U, Xr) are mutually independent for r = 1,…, p. In addition, (e, Zs) are assumed to be independent.
  • (C3) For the predictors, sup1≤i≤n,1≤r≤p,1≤s≤q{|Xnri|, |Znsi|} ≤ B for some bound B [set membership] R. In addition, the predictors Xr satisfy the condition that E(Xr) ≠ 0.
  • (C4) Contamination functions ψ(·) and ϕr(·), 1 ≤ r ≤ p, are twice continuously different-ttiable, satisfying Eψ(U) = 1, Eϕr(U) = 1, and ϕr(·) > 0, 1 ≤ r ≤ p.
  • (C5) The matrices Γnj, j = 1,…, m are nonsingular, i.e. ρ = |infj det(Γnj)| > 0, where equation M78 is a Lnj × (p+q+1) undistorted data matrix in bin j, and equation M79 denotes the kth observation.

The technical conditions above are similar to those introduced in Şentürk and Müller (2006), except for the new independence structure outlined in (C2), the boundedness of the undistorted predictors Zs in (C3), and the bin dependent limiting matrices Γnj in (C5), resulting from the dependence structure between Zs and U. Bounded covariates are standard in asymptotic theory for least squares regression, as are conditions (C2) and (C5) (see Lai, Robbins and Wei, 1979). The identifiability conditions stated in (C4) are equivalent to E(|X) = E(Y|X) and E(Xr|Xr) = Xr.

In the proofs of the main results, the following notations will be utilized.

  1. A [small dot in box] B: The Hadamard product of two matrices, A and B, of the same dimension. The matrix A [small dot in box] B is also of the same dimension with (i, j)th element equal to the product of the (i, j)th elements of matrices A and B.
  2. 1a×b: A matrix of size a × b with all entries equal to one.
  3. [theta w/ hat]nj = ([beta]n0j [beta]n1j,…, [beta]npj, [eta w/ hat]n1j,…, [eta w/ hat]npj)T.
  4. equation M80
  5. equation M81
  6. We use equation M82 to denote the matrix equation M83 and Lnj(i) to denote the number of points in the jth bin such that Uni [set membership] Bnj, and equation M84 is the (r, k)th element of the matrix equation M85 for 1 ≤ r ≤ p + q + 1, where equation M86 is the kth element in the ordered sample equation M87.
  7. equation M88.

Proof of Theorem 1

From Lemma 4 (b.), we have that

equation M89
(12)

where equation M90

Lemma 3 together with (12) implies that, on event En,

equation M91
(13)

We first consider the case of r = 0 and show that equation M92 is asymptotically normal. Using Lemma 4, (13), and some algebra, equation M93 can be expressed as

equation M94

Since the above sum is over all bins indexed by j, and over all points within the bins indexed by k, it is equal to the sum over all data points indexed by i, summed up in a random order.

Thus, the above expression for equation M95 can be further simplified to

equation M96
(14)

Therefore, equation M97 is asymptotically equivalent to equation M98 because the second term equation M99 is negligible when equation M100. Next, let Fn0t be the σ-field generated by equation M101. Then {Sn0t, Fn0t, 1 ≤ tn} is a mean zero martingale for n ≥ 1, since E(Sn0t) = 0, E(Sn0,t+1|Fn0t) = Sn0t, and Sn0t is adapted to Fn0t. Furthermore, note that the σ-fields are nested, that is, Fn0t [subset, dbl equals] Fn0,t+1 for all tn. Hence, it follows from Lemma 1 that equation M102 in distribution (McLeish, 1974, Theorem 2.3 and subsequent discussion). This establishes the asymptotic normality of equation M103.

We proceed next to establish the asymptotic normality of equation M104 for r = 1,…, p. Let equation M105 and note that [gamma with circumflex]nr = [nu with circumflex]nr/[nu with circumflex]nr. We first show that

equation M106
(15)

For (15) to hold, by the Cramér-Wold device, it is enough to show the asymptotic normality of

equation M107
(16)

for real a, b. The asymptotic normality of equation M108 will follow from (15) by applying the δ-method with [gamma with circumflex]nr = [nu with circumflex]nr/[nu with circumflex]nr. Again applying Lemma 4 together with (13) and some simple algebra, we can express [nu with circumflex]nr and [nu with circumflex]nr as

equation M109

Thus, using similar simplifications as was done for the case of r = 0 in (14), the linear combination (16), namely equation M110 can be expressed as

equation M111

The second term equation M112 is asymptotically negligible and it is straightforward to verify that equation M113 is a mean zero martingale for n ≥ 1. Analogous to the case of r = 0, described in more details earlier, it follows from Lemma 2 that equation M114. Finally, a direct application of the δ-method gives equation M115 for 1 ≤ rp, where equation M116 is explicitly given in Theorem 1.

The asymptotic normality of equation M117 follows similarly to the case of equation M118, since they have similar forms in (13). (See also definition/notation 4.) The asymptotic variance equation M119 which has a similar form as equation M120, is given explicitly in Theorem 1. This completes the proof of Theorem 1.

Proof of Theorem 2

The following relation holds on event An

equation M121
(17)

It follows from Lemma 4 (a.) and (b.). Utilizing (17) together with (13) gives

equation M122
(18)

By the Law of Large Numbers, (18), and boundedness considerations

equation M123

It has been shown in Şentürk and Müller (2006) that

equation M124

and it can be shown similarly that equation M125 Also, using Lemma 3, Lemma 4 (a.), (b.) and the Law of Large Numbers, we have

equation M126

The estimators of the asymptotic variances given in Theorem 2, in terms of the above quantities, are: (1) equation M127, (2) equation M128, and (3) equation M129. Thus, the first part of Theorem 2 follows by noting that equation M130 and equation M131. Asymptotic confidence intervals given in the second part of the Theorem follow immediately from Theorem 1 and Slutsky’s theorem using the consistent variance estimators.

8 Discussion

In this work we extend covariate adjusted regression (CAR) models to partial covariate adjusted regression (PCAR) models that allow for the specification of the effects of un-adjusted predictors. Asymptotic normality of the proposed estimators are derived. The PCAR (and CAR) estimation approach was designed to estimate the underlying regression relationship directly, bypassing the estimation of the exact distortion forms. Although the current approach leads to estimators that are simple to implement with known asymptotic properties and good finite sample performance, it does not provide direct estimates of the distorting functions. If the primary interest is in the distorting functions then alternative approaches are needed. A potential alternative approach would involve considering refined estimators for the varying coefficient functions to be used in targeting the distorting functions. However, nontrivial work is required and this remains largely an open problem.

Another interesting and relevant issue, brought to light by a reviewer, is an alternative analysis of the diabetes data given in Section 5. The analysis proceeds by considering the observed variable Zs = Age as equation M132, where U = BMI (body mass index) and equation M133 is a latent variable that can be interpreted as “core biological wear”. This implies that the underlying model of interest would involve the latent construct of core biological wear instead of age (defined to be the length of time that a person has lived, which is observable). Such an application is an interesting extension of CAR in the presence of demographic variables, such as age or gender, and may be of key interest to the area of latent variable modeling in general. However, if one is interested in the direct effect of age on the response Y, then the specific model of interest would include age as a predictor. In this case, the PCAR methodology developed here would be appropriate, as the application of CAR estimation would result in inconsistent estimates. In either case, the specific area of application and relevant body of scientific literature can provide guidance to the researcher in choosing the relevant model. For example, if a researcher is unsure of whether to enter age into the model as actually age or as ‘core biological wear’, then the statistician can provide guidance on focusing/clarifying the specific aims of the research hypotheses with the researcher and the available relevant literature. After the specific research hypothesis is defined, which may involve illiciting the relationship (between the predictor and response) of interest to the re-searcher, then a suitable/reasonable PCAR model can be chosen/entertained. Hence, once the hypothesis (relationship of interest) is identified, then the (underlying) model can be specified (containing predictor Zs or equation M134) and the appropriate estimators (PCAR or CAR) can be applied accordingly.

Finally, we note the following regarding the implementation of the binning procedure, in the context of the data analysis. For the theory, we assumed that the support of U is the interval [a, b]. To be able to bin the data with respect to U = BMI and compute the estimators, we take a and b to be the min/max of the data. For any given population under study, a reasonable range can be inferred to define the limits a and b. For adults, we can reasonably set the limits to 14 and 65 BMI, for instance. In our data, the observed min and max are 18.2 and 49.7 and the intervals/bins are between these observed limits. However, if one is applying the binning using the limits 14 and 65, for instance, our estimator weighted by Lj/n will assign zero weight to bins/intervals with no data (as it should), e.g. bins with U between 14 to 18.2 and 49.7 to 65 BMI. Thus, starting the binning at the min/max of the data is approximately equivalent to starting at [a, b].

Acknowledgment

We are grateful to the reviewers for many detailed suggestions which substantially improved the paper as well as to the Editor for careful review. This work is supported by Grant Number UL1 RR024146 from the National Center for Research Resources (NCRR), a component of the National Institutes of Health (NIH) and NIH Roadmap for Medical Research, NIEHS grant P01-ES011269-06, NIH grants UL1RR024922 and RL1AG032115, and National Institute of Child Health and Human Development grant HD036071.

Appendix

In this section we provide the additional lemmas and their proofs utilized earlier for the main results. We begin by formally defining the events under which the two main theorems of Section 4 are given. Summarizing the existence conditions for the PCAR estimators, define the events

equation M135
(19)

where equation M136, ρ is as defined in (C5), equation M137 is the average of the U’s in Bnj, and (Ω,F, P) is the underlying probability space. The estimators in (8) and (9) are well defined on events Ãn and An, respectively. The event En in Theorem 1 and Theorem 2 is defined to be the intersection of An and Ãn, En = AÃn. It is shown following Lemma 4 that pr(En) → 1 as n → ∞.

We next introduce some additional technical conditions that are needed for the proof of Lemma 4 given below:

  • (C6) The functions h1(u) = ∫ xg1(x, u)dx and h2(u) = ∫ xg2(x, u)dx are uniformly Lips-chitz, where g1(·,·) and g2(·,·) are the joint density functions of (χ,U) and (χe,U), respectively.
  • (C7) The error term satisfies E|eτ| < ∞ for τ > 4.

Lemma 1

Under the technical conditions (C1)–(C7), on event An (19), the martingale differences Wn0t satisfy the conditions

equation M138
(a.)
equation M139
(b.)

.

Proof

Let Wn0t = wn0tυn0t, where equation M140, υn0t = γ0ψ(Unt) + ψ(Unt)entκ1k(t)−γ0 [equivalent] α1nt+α2ntent, α1nt. Using (C1), (C3) and (C4), it holds on event An that sup1≤tn1nt| < c1 and sup1≤tn2nt| < c2 for some c1, c2 > 0. Thus, it holds for ϵ > 0 that

equation M141

Now, equation M142 is bounded uniformly in n and t, since ent has finite fourth moment by (C7). Also note that equation M143. Lemma 1 (a.) follows, since equation M144 uniformly in n and t, equation M145 and |ent| being i.i.d. with finite fourth moments.

Next, consider the term equation M146 given in Lemma 1 (b.). It is equal to

equation M147

Using Law of Large Numbers, it holds that equation M148. Since T4 and T5 have expected values zero and variances O(n−1), they are both Op(n−1/2). By Lemma 4 (a.) and the Law of Large Numbers, term T6 is equal to

equation M149

where Γ and x2133 are as defined prior to Theorem 1. Thus, equation M150 and Lemma 1 (b.) follows.

Lemma 2

Under the technical conditions (C1)–(C7), on event An(19), the martingale differences Wnrt satisfy the conditions

equation M151
(a.)
equation M152
(b.)

Proof

Part (a.) of Lemma 2 follows in a similar fashion as part (a.) of Lemma 1. Therefore, we focus on the proof of part (b.). The term equation M153 in Lemma 2 (b.) is equal to

equation M154

Using the Law of Large Numbers, it holds that equation M155. Since T11, T12, T13 and T14 have expected values zero and variances O(n−1), they are all Op(n−1/2). By Lemma 4 (a.) and the Law of Large Numbers, term T15 is equal to

equation M156

Thus

equation M157

where equation M158, and Σr22 = var(Xr). Hence Lemma 2 (b.) follows.

Lemma 3

Under the technical conditions (C1)–(C6), it holds on event En that,

equation M159

where

equation M160
equation M161

The proof follows from Lemma 3 of Şentürk and Müller (2006) by substituting 1 in place of ϕp+1,…, ϕp+q.

Lemma 4

Under the technical conditions (C1)–(C7), for a sequence rn such that equation M162, on event An

equation M163
(a.)
equation M164
(b.)

, where Γnj is assumed to be nonsingular by (C5), and equation M165.

The proof is similar to the proof of Lemma 4 given in Şentürk and Müller (2006). However a key difference is that the limiting term in part (a.), equation M166, contains expectations taken conditional on U. The conditioning on U does not disappear because of the dependence between Zs and U in the case of PCAR.

Proof that pr(En) → 1. The formula given in (32) in Şentürk and Müller (2006) can be extended to supj |det(Θnj) − det(Γnj)| = Op(rn), where rn is as defined in Lemma 4. This implies, on event An, that pr(infj det(Θnj) > ζ) → 1 as n → ∞, where equation M167 and ρ is as defined in (C5). Similarly, it can be shown that pr(minj Lnjp + q) → 0 as n → ∞, where p + q denotes the number of predictors. Thus, pr(A) → 1 as n → ∞. Furthermore, Lemma 3 implies that

equation M168

This shows that pr(infj det(Θ̃nj) > ζ) → 1 as n → ∞, which implies pr(Ãn) → 1 as n → ∞. Thus pr(En) → 1 as n → ∞.

Footnotes

Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

References

  • Archie JP. Mathematical coupling of data: a common source of error. Annals of Surgery. 1981;193:296–303. [PubMed]
  • Berkson J. Are there two regressions? J. Am. Statist. Assoc. 1950;45:164–180.
  • Carroll RJ, Ruppert D, Stefanski LA, Crainiceanu C. Measurement Error in Nonlinear Models: A Modern Perspective. Second Edition. Boca Raton, FL: Chapman & Hall; 2006.
  • Cleveland WS, Grosse E, Shyu WM. Local regression models. In: Chambers JM, Hastie TJ, editors. Statistical Models in S. Pacific Grove: Wadsworth & Brooks; 1991. pp. 309–376.
  • Hanson RL, Ehm MG, Pettitt DJ, Prochazka M, Thompson DB, Timberlake D, Foroud T, Kobes S, Baier L, Burns DL, Almasy L, Blangero J, Garvey WT, Bennett PH, Knowler WC. An autosomal genomic scan for loci linked to type II diabetes mellitus and body-mass index in Pima Indians. Am. J.Hum. Genet. 1998;63:1130–1138. [PubMed]
  • Hastie T, Tibshirani R. Varying coefficient models. J. R. Statist. Soc. B. 1993;55:757–796.
  • Hwang JT. Multiplicative errors-in-variables models with applications to recent data released by the U.S. department of energy. J. Am. Statist. Assoc. 1986;81:680–688.
  • Iturria S, Carroll RJ, Firth D. Polynomial regression and estimating functions in the presence of multiplicative measurement error. J. R. Statist. Soc. B. 1999;61:547–561.
  • Kaysen GA, Dubin JA, Müller HG, Mitch WE, Rosales LM, Levin NW. the Hemo Study Group. Relationship among inflammation nutrition and physiologic mechanisms establishing albumin levels in hemodialysis patients. Kidney Int. 2003;61:2240–2249. [PubMed]
  • Knowler WC, Pettitt DJ, Saad MF, Charles MA, Nelson RG, Howard BV, Bogardus C, Bennett PH. Obesity in the Pima Indians: its magnitude and relationship with diabetes. Am. J. Clin. Nutr. 1991;53:1543S–1551S. [PubMed]
  • Lai TL, Robbins H, Wei CZ. Strong consistency of least-squares estimates in multiple regression 2. J. Mult. Anal. 1979;9:343–361.
  • McLeish DL. Dependent central limit theorems and invariance principles. Ann. Statist. 1974;2:620–628.
  • Pinter JD, Brown WE, Eliez S, Schmitt JE, Capone GT, Reiss AL. Amygdala and hippocampal volumes in children with Down syndrome: A high-resolution MRI study. Neurology. 2001;56:972–974. [PubMed]
  • Şentürk D, Müller HG. Covariate adjusted regression. Biometrika. 2005;92:75–89.
  • Şentürk D, Müller HG. Inference for covariate adjusted reg ression via varying coefficient models. Ann. Statist. 2006;34:654–679.
  • Smith JW, Everhart JE, Dickson WC, Knowler WC, Johannes RS. Using the ADAP learning algorithm to forecast the onset of diabetes mellitus; Proceedings of the Symposium on Computer Applications and Medical Care; 1988. pp. 261–265.
  • Wu CO, Yu KF. Nonparametric varying coefficient models for the analysis of longitudinal data. Int. Statist. Rev. 2002;70:373–393.