In this work, we use a subset data from a home-based HIV counseling and testing (HBCT) study that was conducted in rural and peri-urban communities in KwaZulu-Natal Province, South Africa, between November 2011 and June 2012 ²¹ . The data were obtained from the Human Sciences Research Council (HSRC) of South Africa ²¹ . This study was conducted to provide a better understanding of the complexity, severity, and prevalence of non-communicable disease (NCDs) in a community known to have one of the highest rates of HIV incidence and prevalence in the world ²¹ .

Home-based HIV counseling and testing is a cross-sectional, single-site study in South Africa that aims to increase engagement in HIV care by integrating NCDs screening with community-based HIV testing ²² . A random sampling approach was used, where 587 participants over the age of 18 were selected from 50,000 people living in the Mpumuza suburb ²¹ . Anthropometric and biological measures were collected in the survey with the purpose of establishing the prevalence of a range of NCDs and associated risk factors. Eligible individuals participated in a face-to-face interview, physical, psychological and clinical examinations. Persons younger than 18 years living in Mpumuza and all household members not previously enrolled, and members unable to give written consent were excluded from the study. Mobile phones were used for data collection to increase efficiency in data capture and analysis ²¹ .

In our study, we used a subset data consisting of 76 individuals who self-reported the number of cigarettes smoked, fruit and vegetable consumption. We use the dataset to illustrate the multivariate method in modeling the amount of association between body mass index (BMI) and three exposures (smoking, fruit, and vegetable intakes). BMI was measured in kg/ m ², while smoking was measured as the average number of cigarettes smoked per day. Initially, fruit and vegetable intakes were measured in terms of the number of servings consumed per day. It is often assumed that a standard portion of fruit/vegetable weighs about 80g ⁵ . Therefore, for this study, we converted the number of servings to grams per day (g/day) by multiplying the reported number of servings by 80g. The subset data has the following three properties that make it suitable for use in this work: (1) measurement error in the recorded number of cigarettes smoked due to possible misreporting, (2) measurement error in fruits and vegetable consumption due to recall bias, and conversion of the number of servings of fruits and vegetables into grams, and (3) the measurement error in the three exposures is often correlated, for instance, smokers are likely to over-report fruit and vegetable intakes due to their beneficial effects, and to under-report the number of cigarettes they smoke due to the associated harmful effects. Epidemiologically, BMI is used as a risk factor of a health outcome. However, in this study, we model BMI as an outcome as in other several studies, for instance, ^{23– 26} . The subset data is only used to illustrate the method and not to draw inference.

Ethics approval was granted by both HSRC Research Ethics Committee (REC: 1/26/05/11) and the University of Washington Institutional Review Board (48733). Informed written consent was obtained from each participant in the study. Participants were provided with written information on the study (including the study’s background and objectives) and their rights regarding participation and withdrawal at any time.

An interest in epidemiological study could be to investigate the association between BMI and three exposures namely: fruit, vegetable and smoking using the multiple linear regression

Y = β 0 + β X 1 X 1 + β X 2 X 2 + β X 3 X 3 + ϵ , ( 1 )

where Y denotes the BMI, β ₀ is the intercept, β _{X 1} , β _{X 2} and β _{X 3} are the coefficient parameters for the true long-term fruit ( X ₁), vegetable ( X ₂) and cigarette ( X ₃) intake respectively and ϵ is the random error term. In this study, we use vegetable intake and cigarette smoking as confounders and assume that the main interest is in estimating β _{X 1} . In practice, the true intakes are unobservable and, therefore, the intakes recorded in self-reported questionnaires are used. Let W ₁, W ₂ and W ₃ denote the measured versions of X ₁, X ₂ and X ₃, respectively. The use of W _p ’s in place of X _p ’s, ( p = 1, 2, 3), in Equation (1) yields biased estimates β ^ W 1 , β ^ W 2 and β ^ W 3 of β _{X 1} , β _{X 2} and β _{X 3} respectively. Let β ^ W = ( β ^ W 1 , β ^ W 2 , β ^ W 3 ) T .

We assumed that the observed exposures are related to the true exposures with additive measurement error as

W i = α 0 i + α 1 i X i + ϵ W i , i = 1 , 2 , 3 ( 2 )

where ϵ _W = ( ϵ _{W 1} , ϵ _{W 2} , ϵ _{W 3} ) ^⊤, ϵ _W ∼ N(0, Σ _{ϵ W} ); W = ( W ₁, W ₂, W ₃) ^⊤; α ₀ = ( α ₀₁, α ₀₂, α ₀₃) ^⊤, α ₁ = ( α ₁₁, α ₁₂, α ₁₃) ^⊤; with the terms in α ₀ and α ₁ quantifying the constant bias and the proportional scaling bias respectively; ϵ _W is a random error term, ϵ _{W i} is assumed to be independent of the true exposure X _i and the systematic bias components, α _{0 i} and α _{1 i} .

A univariate method. In a univariate case, the bias in the association between an outcome and an exposure is adjusted by dividing the unadjusted association estimate by the attenuation factor ¹⁹ . Attenuation factor ( λ _i ) is defined as λ _i = var( X _i )/ var( W _i ), i.e., the ratio of the variance of the true exposure to the variance of the observed exposure, also referred to as reliability ratio. This method ignores correlations between the errors and also the correlation between the true exposures.

Multivariate method. We propose and describe a general approach for handling p-exposures (p≥2) measured with correlated errors. For simplicity and without loss of generality, we assume that W _i is measured without systematic bias (i.e., α _{0 i} = 0, α _{1 i} = 1 in Equation 2). For multiple exposures measured with correlated errors, the adjusted association estimates can be obtained by pre-multiplying the unadjusted association estimates by the inverse of the transpose of attenuation-contamination matrix as

β ^ X ∗ = ( Λ ^ p T ) − 1 β ^ W ∗ , ( 3 )

where β ^ X ∗ and β ^ W ∗ denotes vectors of true and biased coefficients for the p-exposures respectively and Λ _p denotes a p × p attenuation-contamination matrix ^{19, 27} . The off-diagonal elements in Λ are known as contamination factors while the diagonal elements are called attenuation factors ¹⁴ . Noteworthy, the attenuation factor quantifies the bias in the association between an outcome and an exposure. In contrast, the contamination factor quantifies the effect of measurement error in one exposure variable on the other exposure variable’s estimate. β ^ W ∗ in Equation (3) can be obtained from the observed questionnaire data.

In the multiple exposures case, the estimate of attenuation-contamination matrix Λ ^ p is defined as

Λ ^ p = [ σ ^ X 1 2 σ ^ X 2 X 1 ⋮ σ ^ X p X 1 σ ^ X 1 X 2 σ ^ X 2 2 ⋮ σ ^ X p X 2 ⋯ ⋯ ⋰ ⋯ σ ^ X 1 X p σ ^ X 2 X p ⋮ σ ^ X p 2 ] ︸ Σ ^ X * [ σ ^ W 1 2 σ ^ W 2 W 1 ⋮ σ ^ W p W 1 σ ^ W 1 W 2 σ ^ W 2 2 ⋮ σ ^ W p W 2 ⋯ ⋯ ⋰ ⋯ σ ^ W 1 W p σ ^ W 2 W p ⋮ σ ^ W p 2 ] ︸ Σ ^ W ∗ − 1 − 1 , ( 4 )

where Σ ^ X * is the estimate of covariance matrix of the true exposures, Σ ^ W *− 1 is the inverse of the estimate of covariance matrix of the measured exposures, σ ^ X i 2 is the variance estimate of X _i ( i = 1, 2, ..., p) ; σ ^ X i X j ( j = 1, 2, ..., p; i ≠ j) denotes the covariance estimate between the true exposures; σ W i 2 is the variance estimate of W _i ; σ ^ W i W j ( i ≠ j) is the covariance estimate between the observed exposures.

The elements of the variance-covariance matrix of the observed exposures, Σ W ∗ , are estimated from the observed data. The variances of the true exposures, σ X i 2 ’s, can be estimated using validity coefficients for the questionnaire. According to Kipnis et al. ⁶ , the validity coefficient is given by:

ρ W i X i = c o v ( W i ,     X i ) v a r ( W i ) v a r ( X i ) , i = 1 , 2 , … , p . = σ X i σ W i , ( 5 )

where W _i is assumed to be the measured with error term only and ϵ _{W i} is assumed to be independent of X _i . From Equation (5), we estimate the variance of the true exposures as

σ ^ X i 2 = ( ρ ^ W i X i σ ^ W i ) 2 , ( 6 )

by incorporating external validation information on ρ _{W i X i} . To obtain covariances between the true exposures, one of the following two approaches is used: (i) if external information about the correlation between true exposures (i.e. ρ ^ X i X j ) is available, we obtain covariances between true exposures as follows:

σ ^ X i X j = ρ ^ X i X j σ ^ X i σ ^ X j , i ≠ j , ( 7 )

where σ ^ X i are obtained as shown in Equation (6); (ii) if we can obtain prior information about the correlation between the errors in the observed exposures, ρ ^ ϵ W i ϵ W j , we can solve for σ ^ X i X j by decomposing the covariance of observed exposures into unknown covariance between true exposures and unknown covariance between errors as follows:

σ ^ W i W j = σ ^ X i X j + σ ^ ϵ W i ϵ W j + σ ^ X i ϵ W j ︸ 0 + σ X j ϵ W i ︸ 0 = σ ^ X i X j + ρ ^ ϵ W i ϵ W j σ ^ ϵ W i σ ^ ϵ W j , ( 8 )

where X _i and ϵ _{W j} , X _j and ϵ _{W i} are assumed to be uncorrelated.

From Equation (2) and Equation (6), the estimate of the error variance σ ^ ϵ W i 2 is

σ ^ ϵ W i 2 = σ ^ W i 2 − σ ^ W i 2 ρ ^ W i X i 2 ︸ σ ^ X i 2 , = σ ^ W i 2 ( 1 − ρ ^ W i X i 2 ) , ( 9 )

See Appendix B of the extended data ²⁸ for the proof.

From Equation (8)– Equation (9), the covariances between the true exposures are given by

σ ^ X i X j = σ ^ W i W j − ρ ^ ϵ W i ϵ W j σ ^ W i σ ^ W j ( 1 − ρ ^ W i X i 2 ) ( 1 − ρ ^ W j X j 2 ) , ( i ≠ j ) , ( 10 )

Using the observed data and external information, we can determine all the terms required to estimate the attenuation-contamination matrix, Λ, as shown in Equation (4) and adjust for the bias in the association between the exposures measured with error and the outcome using Equation (3).

We illustrate a method that accounts for uncertainty in the validity measures attributable to heterogeneity in the study populations and in parameter estimation. The proposed Bayesian method applies Markov Chain Monte Carlo (MCMC) estimation approach to combine observed self-reported data and external validation data in adjusting for measurement error in three exposures measured with correlated errors. MCMC is a class of algorithms that samples from the posterior distributions by traversing the parameter space ²⁹ . The posterior distribution is obtained by updating the prior distribution with observed data. The steps for implementing the trivariate method are described below.

We first obtained external information on validity coefficients and generated validity coefficients for use by interpreting the lower and upper limits obtained from the literature as the 95% credible intervals (CIs) of the distribution of possible values respectively. Due to the skewed distribution of validity coefficients, Fisher’s transformation was used to generate the validity coefficients as explained in the next section.

Second, for the observed exposures, we estimated the posterior distribution of the covariance matrix (Σ _W ). The exposures were assumed to follow a multivariate normal distribution with mean and covariance, i.e., W ∼ N ₃( µ _W , Σ _W ). We assumed a weakly informative multivariate normal prior for µ _W as µ _{W prior} ∼ N ₃(0,10 ⁶ I ₃), where I ₃ is a 3 × 3 identity matrix. In a multivariate normal distribution, Σ _W must satisfy two conditions: (1) be positive definite (i.e. W ^T Σ _W W > 0, for all W ) and (2) be a symmetric matrix. The semi-conjugate prior distribution for Σ _W , which has these two properties, is the inverse-Wishart distribution ²⁹ . To minimize the influence of the prior information on the estimate of Σ _W , we considered weakly informative inverse-wishart prior as Σ _{W prior} ∼ IW( I ₃, v), where v = 3 is the degrees of freedom.

Third, using the validity coefficients generated from the external data and the posterior distribution of covariance matrix for observed exposures, we estimated the variance of true intakes, σ ^ X i 2 ( i = 1, 2, 3), using the relationship given in Equation (6) so that

σ ^ X 1 2 = ( ρ ^ W 1 X 1 σ ^ W 1 ) 2 σ ^ X 2 2 = ( ρ ^ W 2 X 2 σ ^ W 2 ) 2 σ ^ X 3 2 = ( ρ ^ W 3 X 3 σ ^ W 3 ) 2 . ( 11 )

The covariances between true intakes ( σ ^ X i X j ; j = 1, 2, 3) were estimated as,

σ ^ X 1 X 2 = σ ^ W 1 W 2 − ρ ^ ϵ W 1 ϵ W 2 σ ^ W 1 σ ^ W 2 ( 1 − ρ ^ W 1 X 1 2 ) ( 1 − ρ ^ W 2 X 2 2 ) σ ^ X 1 X 3 = σ ^ W 1 W 3 − ρ ^ ϵ W 1 ϵ W 3 σ ^ W 1 σ ^ W 3 ( 1 − ρ ^ W 1 X 1 2 ) ( 1 − ρ ^ W 3 X 3 2 ) σ ^ X 2 X 3 = σ ^ W 2 W 3 − ρ ^ ϵ W 2 ϵ W 3 σ ^ W 2 σ ^ W 3 ( 1 − ρ ^ W 2 X 2 2 ) ( 1 − ρ ^ W 3 X 3 2 ) . ( 12 )

by incorporating external validation information on correlation between the errors ( ρ _{ϵ W i ϵ W j} ). We generated the correlation between errors from a plausible range guided by correlation in the observed data and prior expert information on the most likely sign of the correlation between the exposures, as described in the next section.

Having obtained the covariance matrices of the true and observed exposures, we estimated the attenuation-contamination matrix (Λ ₃) from their joint distribution as

Λ ^ 3 = [ σ ^ X 1 2 σ ^ X 1 X 2 σ ^ X 1 X 3 σ ^ X 1 X 2 σ ^ X 2 2 σ ^ X 2 X 3 σ ^ X 1 X 3 σ ^ X 2 X 3 σ ^ X 3 2 ] ︸ Σ ^ X [ σ ^ W 1 2 σ ^ W 1 W 2 σ ^ W 1 W 3 σ ^ W 1 W 2 σ ^ W 2 2 σ ^ W 2 W 3 σ ^ W 1 W 3 σ ^ W 2 W 3 σ ^ W 3 2 ] ︸ Σ ^ W − 1 − 1 , ( 13 )

where Σ ^ X is the estimate of covariance matrix of the three true exposures, Σ ^ W − 1 is the inverse of the estimate of covariance matrix of the three measured-with-error exposures, σ ^ X i 2 is the variance estimate of X _i ( i = 1, 2, 3); σ ^ X i X j ( i ≠ j) denotes the covariance estimate between the true exposures; σ ^ W i 2 is the variance estimate of W _i ; σ ^ W i W j ( i ≠ j) is the covariance estimate between the observed exposures.

Lastly, we fitted a Bayesian multiple linear regression model (hereafter, naive method) to obtain the posterior distributions of the unadjusted coefficient estimates β ^ W = ( β ^ W 1 , β ^ W 2 , β ^ W 3 ) T . In the naive model, we assumed weakly informative normal independent priors by choosing a very small precision (large variance) for the unadjusted coefficient estimates as β _{W i prior} ∼ N(0, 10 ⁶). The adjusted coefficient estimates β ^ X were then obtained from the joint posterior distribution of Λ ^ 3 and β ^ W as

β ^ X = ( Λ ^ 3 T ) − 1 β ^ W . ( 14 )

We implemented the trivariate method in R version 3.6.3 using rjags (version 4-10), coda (version 0.19-3), MCMCpack (version 1.4-9), and mvtnorm (version1.1-1) packages. To facilitate Bayesian estimation of the covariance matrix of the observed exposures (Σ _W ), rjags package was used to provide an interface from R to the JAGS library ³⁰ . JAGS is a gibbs sampler that uses MCMC to draw dependent samples from the posterior distribution of the parameters ³¹ . The Bayesian estimation of Σ _W proceeded in the following steps: (1) defining a model for Σ _W under Bayesian inference using gibbs sampling (BUGS) algorithm in a stand alone file, (2) reading the model file using the jags.model function, (3) updating the model using the update method for jags objects and (4) extracting the posterior samples of the model using the coda.samples function from the coda package.

MCMCregress function from the MCMCpack package was used to generate a posterior density sample from the naive linear regression model ³² . MCMC convergence diagnostics of all the model parameters was done using trace plots and autocorrelation (ACF) plots from the coda package ³³ . See extended data: Appendix C ²⁸ for convergence diagnostics results. For each model, the burn-in iterations were set to 2,000 and 10,000 MCMC iterations were run after the burn-in iterations. Every first sample value was kept in the MCMC simulations by using a thinning interval of 1. When compiling a JAGS model, an initial sampling step may be needed during which the samplers learn their behaviour to maximize their performance ³⁴ . Therefore, the number of iterations for adaptation in the the jags model was set to 500. The results were presented in terms of density plots, posterior mean and median. We compared the results obtained under naive, univariate, and trivariate methods. The R code used for analysis is presented in the extended data ²⁸ .

External information on the validity coefficient and error correlations for fruit, vegetable, and cigarette information was obtained from the literature. According to Kaaks et al. ¹ , the validity coefficient of self-reported fruit intake ranges from 0.33 to 0.79, while that of vegetable intake ranged from 0.30 to 0.60. A meta-analysis study on the validity of questionnaires assessing fruit and vegetable consumption by Collese et al. ² reported validity coefficients of 0.26 for vegetables and 0.49 for fruits. Other similar validation studies reported validity coefficients in the aforementioned ranges for fruits and vegetables ^{3, 4, 35} . Therefore, based on these information we considered a range of 0.3 to 0.8 for fruits and a range of 0.25 and 0.7 for vegetables.

In the Scottish Heart Health Study of 2,849 men and 2,900 women ³⁶ , the correlation between the self-reported number of cigarettes and biochemical measures was reported between 0.67 and 0.72. In a study on the validation of self-reported smoking by analysis of hair for nicotine and cotinine ³⁷ , the validity coefficient between the number of cigarettes smoked per day and nicotine/cotinine levels in hair and plasma was found to be between 0.48 and 0.63, while the correlation between the average number of cigarettes smoked and carboxyhemoglobin was 0.70. In a follow-up study to examine the relationships among self-reported cigarette consumption, exhaled carbon monoxide, and urinary cotinine/creatinine ratio in pregnant women ³⁸ , a validity coefficient in the range of 0.61 to 0.70 was reported. A study by Stram et al. ³⁹ found the correlation between the self-reported number of cigarettes smoked and the true lung dose to be between 0.40 and 0.70, and this range was consistent with the findings from the previously discussed related validation studies. Based on this information, we considered a validity coefficient range of 0.40 and 0.70.

We generated the correlation between errors from plausible ranges that were determined based on the correlation in the observed data and the most probable sign of the correlation among fruits, vegetables, and cigarettes as explained below:

a. Since the correlation coefficient between fruit and vegetable intake in the observed data was positive, we also assumed the error correlation between fruit and vegetables to be mostly positive;

b. An investigation on the correlation coefficient between cigarette smoking and fruits/vegetable intake in the observed data showed a negative correlation coefficient. Based on this and the fact that persons who tend to overstate fruit and vegetable consumption are likely to understate the number of cigarettes smoked, we assumed the error correlation to be mostly negative.

We obtained the upper limits of error correlations by assuming that the error covariance equals the covariance in the observed data and set the lower limit of the error correlation to zero, based on the assumption that the covariance in the observed data equals the covariance between the true intakes ¹⁴ .

Using the range of plausible values obtained from external validation information, we generated the validity coefficients using the Fisher-Z transformation method by assuming that the reported lower and upper limits are 0.05 and 0.95 quantiles of the uncertainty distribution, respectively. Fisher Z-transformation is a commonly used method to transform the sampling distribution of correlation coefficients to become approximately normally distributed ^{40, 41} . The procedure is as outlined below:

We investigated how varying the level of uncertainty assumed for the limits of the validity coefficients reported from literature affected the estimates for fruit, vegetable, and the average number of cigarettes smoked. We also investigated how the estimates varied with the magnitude of the correlation between errors in fruit and vegetable intake, fruit and cigarette smoking, and vegetable and cigarette smoking. This helps determine the estimates’ sensitivity to various magnitudes of CI and the correlation between errors when using the multivariate method.

Data used in this study are made available to the researcher upon registration and agreeing to the terms and conditions of use in the HSRC web site at http://curation.hsrc.ac.za/ Dataset-565-datafiles.phtml.

Figshare: A Method to Adjust for Measurement Error in Multiple Exposures Measured with Correlated Error in the Absence of Internal Validation Study-Supplementary materials. https://doi.org/10.6084/m9.figshare.13147970.v2 ²⁸

The file shows the validity coefficient derivation, Proof for the estimate of error variance, R code for implementing the methods and convergence diagnostics results (i.e. Trace plots and ACF plots for the standard deviation and naive regression coefficient estimates of the fruits, vegetables and average number of cigarettes smoked, with explanation) and the sensitivity analysis results (supporting Tables) for varying the magnitude of error correlation between the exposures.

The extended data are available under the terms of the Creative Commons Zero (CC0) license.