{'fraction_non_alphanumeric': 0.03545938082465791, 'fraction_numerical': 0.018775287539209894, 'mean_word_length': 4.484976935544196, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A new multivariate measurement error model with zero-inflated dietary data, and its application to dietary assessment

In the United States the preferred method of obtaining dietary intake data is the 24-hour dietary recall, yet the measure of most interest is usual or long-term average daily intake, which is impossible to measure. Thus, usual dietary intake is assessed with considerable measurement error. Also, diet represents numerous foods, nutrients and other components, each of which have distinctive attributes. Sometimes, it is useful to examine intake of these components separately, but increasingly nutritionists are interested in exploring them collectively to capture overall dietary patterns. Consumption of these components varies widely: some are consumed daily by almost everyone on every day, while others are episodically consumed so that 24-hour recall data are zero-inflated. In addition, they are often correlated with each other. Finally, it is often preferable to analyze the amount of a dietary component relative to the amount of energy (calories) in a diet because dietary recommendations often vary with energy level. The quest to understand overall dietary patterns of usual intake has to this point reached a standstill. There are no statistical methods or models available to model such complex multivariate data with its measurement error and zero inflation. This paper proposes the first such model, and it proposes the first workable solution to fit such a model. After describing the model, we use survey-weighted MCMC computations to fit the model, with uncertainty estimation coming from balanced repeated replication.

1. Introduction. This paper presents statistical models and methodology to overcome a major stumbling block in the field of dietary assessment. More nutritional background is provided in Section 2: a summary of the key conceptual issues follows: • Nutritional surveys conducted in the United States typically use 24-hour (24 h) dietary recalls to obtain intake data, that is, an assessment of what was consumed in the past 24 hours. • Because dietary recommendations are intended to be met over time, nutritionists are interested in "usual" or long-term average daily intake. • Dietary intake is thus assessed with considerable measurement error.
• Consumption patterns of dietary components vary widely; some are consumed daily by almost everyone, while others are episodically consumed so that 24-hour recall data are zero-inflated. Further, these components are correlated with one another. • Nutritionists are interested in dietary components collectively to capture patterns of usual dietary intake, and thus need multivariate models for usual intake. • These multivariate models for usual intakes, taking into account episodically consumed foods, do not exist, nor do methods exist for fitting them.
One way to capture dietary patterns is by scores, although our work is not limited to scores. The Healthy Eating Index-2005(HEI-2005, described in detail in Section 2, is a scoring system based on a priori knowledge of dietary recommendations, and is on a scale of 0-100. Ideally, it consists of the usual intake of 6 episodically consumed and thus 24 h-zero inflated foods, 6 daily-consumed dietary components, adjusts these for energy (caloric) intake, and gives a score to each component. The total score is the sum of the individual component scores. Higher scores indicate greater compliance with dietary guidelines and, therefore, a healthier diet. Here are a few questions that nutritionists have not been able to answer, and that our approach can address: • What is the distribution of the HEI-2005 total score, and what % of Americans are eating a healthier diet defined, for example, by a total score exceeding 80? • What is the correlation between the individual score on each dietary component and the scores of all other dietary components? • Among those whose total HEI-2005 score is >50 or ≤50, what is the distribution of usual intake of whole grains, whole fruits, dark green and orange vegetables and legumes (DOL) and calories from solid fats, alcoholic beverages and added sugars (SoFAAS)? • What % of Americans exceed the median score on all 12 HEI-2005 components?
In this paper, to answer public health questions such as these that can have policy implications, we build a novel multivariate measurement error model for estimating the distributions of usual intakes, one that accounts for measurement error and zero-inflation, and has a special structure associated with the zero-inflation. Previous attempts to fit even simple versions of this model, using nonlinear mixed effects software, failed because of the complexity and dimensionality of the model. We use survey-weighted Monte Carlo computations to fit the model with uncertainty estimation coming from balanced repeated replication. The methodology is illustrated using the HEI-2005 to assess the diets of children aged 2-8 in the United States. This work represents the first analysis of joint distributions of usual intakes for multiple food groups and nutrients.
The paper is outlined as follows. In Section 2 we give the background for the data we observe. In particular, we provide more information about the HEI-2005. Section 3 describes our model which is a highly nonlinear, zero-inflated, repeated measures model with multiple latent variables. The model also has a patterned covariance matrix with structural zeros and ones. We derive a parameterization that allows estimated covariance matrices to be actual covariance matrices. We also define technically what we mean by usual intake, and illustrate the use of simulation methods used to answer the questions posed above, as well as many others.
Section 4 describes our estimation procedure. Previous attempts using nonlinear mixed effects models to estimate the distribution of episodically consumed food groups [Tooze et al. (2006); Kipnis et al. (2009)] do not work here because of the high dimensionality of the problem. We instead develop a Monte Carlo strategy based on the idea of Gibbs sampling; although because of sampling weights, we treat the method as a frequentist (non-Bayesian) one. This section describes some of the basics of the methodology; the full technical details of implementation are given in the Appendix.
Section 5 describes the analysis of the HEI-2005 components using the 2001-2004 National Health andNutrition Examination Survey (NHANES) for children ages 2-8. Important contextual points arise because of the nature of the data. For example, if whole grains are consumed, then necessarily total grains are consumed with probability one, a restriction that a naive use of our model cannot handle. We develop a simple novel device to uncouple consumption variables that are tightly linked in this way. Finally, in this section we provide the first answers to the four questions we have posed. In Section 6 we discuss various additional aspects of the problem and the data analysis. Concluding remarks and a policy application are given in Section 7.

Data and the HEI-2005 scores.
Here we give more detail about the nutrition context that motivates this work.
In surveys conducted in the United States, the preferred method of obtaining intake data is the 24-hour dietary recall because it limits respondent burden and facilitates accurate reporting; yet the measure of greatest interest is "usual" or long-term average daily intake. Thus, dietary intake is assessed with considerable measurement error. Also, diets are comprised of numerous foods, nutrients and other components, each of which may have distinctive attributes and effects on nutritional health. Sometimes, it is useful to examine intake of these components separately, but increasingly nutritionists are interested in exploring them collectively to capture patterns of dietary intake. Consumption patterns of these components vary widely; some are consumed daily by almost everyone, while others are episodically consumed so that 24-hour recall data are zero-inflated. In addition, these various components are often correlated with one another. Finally, it is often preferable to analyze the amount of a dietary component relative to the amount of energy (calories) in a diet because dietary recommendations often vary with energy level, and this approach provides a way of standardizing dietary assessments.
One of the US Department of Agriculture's (USDA's) strategic objectives is "to promote healthy diets" and it has developed an associated performance measure, the Healthy Eating Index-2005(HEI-2005 Table 1 for a list of these components and the standards for scoring, and see  for details. Intakes of each For saturated fat, density is 9 × 100 usual saturated fat (grams) divided by usual calories, that is, the percentage of usual calories coming from usual saturated fat intake. For SoFAAS, the density is the percentage of usual intake that comes from usual intake of calories, that is, the division of usual intake of SoFAAS by usual intake of calories. Here, "DOL" is dark green and orange vegetables and legumes. Also, "SoFAAS" is calories from solid fats, alcoholic beverages and added sugars. The total HEI-2005 score is the sum of the individual component scores.
food or nutrient, represented by one of the 12 components, are expressed as a ratio to energy intake, assessed, and ascribed a score. The HEI-2005 is used to evaluate the diets of Americans to assess compliance with the 2005 Dietary Guidelines, yet use of the HEI-2005 is limited by the challenges described above. Until recently, there have been no solutions to these challenges, so published evaluations have been limited to analyses of mean scores for the population and various subgroups. Freedman et al. (2010) have described a method of estimating the population distribution of a single component of HEI-2005, and the prevalence of high or low scores on that component; but there has been to date no satisfactory way to determine the prevalence of high or low total HEI-2005 scores, considering all of its interrelated components simultaneously. In addition, answers to the complex questions posed in the Introduction remain unavailable. This paper aims to provide a means to do these crucial evaluations.
The 12 HEI-2005 components represent 6 episodically consumed food groups (total fruit, whole fruit, total vegetables, dark green and orange vegetables and legumes or DOL, whole grains and milk), 3 daily-consumed food groups (total grains, meat and beans and oils) and 3 other daily-consumed dietary components (saturated fat; sodium; and calories from solid fats, alcoholic beverages and added sugars, or SoFAAS). The classification of food groups as "episodically" and "daily" consumed is based on the number of individuals who report them on 24 h recalls. If there are only a few zeros for a component, we treat that as a daily-consumed food, and replace all zeros with 1/2 the minimum value of the nonzeros for that food. However, the crucial statistical aspect of the data is that six of the food groups are zeroinflated. The percentages of reported nonconsumption of total fruit, whole fruit, whole grains, total vegetables, DOL and milk on any single day are 17%, 40%, 42%, 3%, 50% and 12%, respectively.
We are interested in the usual intake of foods for children aged 2-8. The data available to us, described in more detail in Section 5, came from the National Health and Nutrition Examination Survey, 2001-2004. The data used here consisted of n = 2,638 children, each of whom had a survey weight w i for i = 1, . . . , n. In addition, one or two 24 h dietary recalls were available for each individual. Along with the dietary variables, there are covariates such as age, gender, ethnicity, family income and dummy variables that indicate a weekday or a weekend day, and whether the recall was the first or second reported for that individual.
Using the 24 h recall data reported, for each of the episodically consumed food groups, two variables are defined: (a) whether a food from that group was consumed; and (b) the amount of the food that was reported on the 24 h recall. For the 6 daily-consumed food groups and nutrients, only one variable indicating the consumption amount is defined. In addition, the amount of energy that is calculated from the 24 h recall is of interest. The number of 7 dietary variables for each 24 h recall is thus 12 + 6 + 1 = 19. The observed data are Y ijk for the ith person, the jth variable and the kth replicate, j = 1, . . . , 19 and k = 1, . . . , m i . In the data set, at most two 24 h recalls were observed, so that • Y i,2ℓ−1,k = Indicator of whether dietary component #ℓ is consumed, with ℓ = 1, 2, 3, 4, 5, 6. • Y i,2ℓ,k = Amount of food #ℓ consumed. This equals zero, of course, if none of food #ℓ is consumed, with ℓ = 1, 2, 3, 4, 5, 6. • Y i,ℓ+6,k = Amount of nonepisodically consumed food or nutrient #ℓ, with ℓ = 7, 8, 9, 10, 11, 12. • Y i,19,k = Amount of energy consumed as reported by the 24 h recall.

Model and methods.
3.1. Basic model description. Our model is a generalization of work by Tooze et al. (2006) and Kipnis et al. (2009) for a single food and Kipnis et al. (2011) and  for a single food and nutrient. Observed data will be denoted as Y , and covariates in the model will be denoted as X. As is usual in measurement error problems, there will also be latent variables, which will be denoted by W .
We use a probit threshold model. Each of the 6 episodically consumed foods will have 2 sets of latent variables, one for consumption and one for amount, while the 6 daily-consumed foods and nutrients as well as energy will have 1 set of latent variables, for a total of 19. The latent random variables are ε ijk and U ij , where (U i1 , . . . , U i,19 ) = Normal(0, Σ u ) and (ε i1k , . . . , ε i,19,k ) = Normal(0, Σ ε ) are mutually independent. In this model, food ℓ = 1, . . . , 6 being consumed on day k is equivalent to observing the binary Y i,2ℓ−1,k , where where g tr (y, λ) = √ 2{g(y, λ)−µ(λ)}/σ(λ), g(y, λ) is the usual Box-Cox transformation with transformation parameter λ, and {µ(λ), σ(λ)} are the sample mean and standard deviation of g(y, λ), computed from the nonzero food data. This standardization is simply a convenient device to improve the numerical performance of our algorithm without affecting the conclusions of our analysis. The reported consumption of daily consumed foods or nutrients ℓ = 7, . . . , 12 is modeled as Finally, energy is modeled as As seen in (3.2)-(3.4), different transformations (λ 1 , . . . , λ 13 ) are allowed to be used for the different types of dietary components; see Section A.12.
So that we can handle any number of episodically consumed dietary components and any number of daily consumed components, suppose that there are J episodically consumed dietary components, and K daily consumed dietary components, and in addition there is energy. Then the restrictions defined above lead to the covariance matrix The difficulty with parameterizations of (3.6) is that the cells that are not constrained to be 0 or 1 cannot be left unconstrained, otherwise (3.6) need not be a covariance matrix, that is, positive semidefinite.
We have developed an unconstrained parameterization that results in the structure (3.6). Consider an unconstrained lower triangular matrix V and define Σ ε = V V T . This is positive semidefinite and therefore qualifies Σ ε as a proper covariance matrix. The form of V is To achieve the desired pattern (3.6), we derive the following four restrictions: The third restriction can be ensured by the further parameterization where q = 2, 3, . . . , J − 1, |r t | ≤ 1, t = 1, . . . , J − 1, and |θ s | ≤ π, s = 1, . . . , Similarly, the fourth restriction can be further expressed by setting 3.3. The use of sampling weights. As described in the Appendix, we used the survey sample weights from NHANES both in the model fitting procedure and, after having fit the model, in estimating the distributions of usual intake.
While not displayed here, we redid the model fitting calculations without weighting, because the covariates we use are major players in determining the sampling weights, hence, it is reasonable to believe that the model in Section 3 holds both in the sample and in the population. When we did this, the parameter estimates were essentially unchanged.
Thus, we use the sampling weights only for estimation of the population distributions. We actually did this for the purpose of handling the clustering in the sample design. For such a complex statistical procedure as ours, we knew we could not do theoretical standard errors, so we thought about the bootstrap, and realized that putting together a bootstrap for the complex survey would be nearly impossible. However, we already had developed a set of Balanced Repeated Replication (BRR) weights [Wolter (1995)]; see Section 5.7 for details. These BRR weights have the property that, in the frequentist survey sampling sense, they appropriately reflect the clustering in the standard error calculations.
Of course, the use of sampling weights in the modeling provide unbiased estimates of the (super) population parameters of interest. In addition, the use of sampling weights in the distribution estimation provides an estimated distribution that is representative of the US population, not just the sample.
3.4. Distribution of usual intake and the HEI-2005 scores. We assume here that estimates of Σ u , Σ ε and β j for j = 1, . . . , 19 have been constructed; see Section 4. Here we discuss what we mean by usual intake for an individual, how to estimate the distribution of usual intakes, how to convert usual intakes into HEI-2005 scores, and how to assess uncertainty.
Consider the first episodically consumed dietary component, a food group, with reporting being done on a weekend. Set X i1,wkend and X i2,wkend to be the versions of X i1k and X i2k where the dummy variable has the indicator of the weekend and that the recall is the first one. Following Kipnis et al. (2009), we define the usual intake for an individual on the weekend to be the expectation of the reported intake conditional on the person's random effects U i . Let the (q, p) element of Σ ε be denoted as Σ ε,q,p . As in Kipnis et al. (2009) Detailed formulas for this are given in Appendix A.11. Then, following the convention of Kipnis et al. (2009), the person's usual intake of the first episodically consumed dietary component on the weekend is defined as . Similarly, let X i1,wkday and X i2,wkday be as above but the dummy variable is appropriate for a weekday. Then the person's usual intake of the first episodically consumed food group on weekdays is defined as . Finally, the usual intake of the first episodically consumed food for the individual is T i1 = (4T i1,wkday + 3T i1,wkend )/7, since Fridays, Saturdays and Sundays are considered to be weekend days. Usual intake for the other episodically consumed food groups is defined similarly.
A person's usual intake of a daily-consumed food group/nutrient and energy on the original scale is defined similarly. Consider, for example, energy, which is the 13th dietary component and the 19th set of terms in the model. Let X i,19,wkend and X i,19,wkday be the versions of X i,19,k where the dummy variable has the indicator of the weekend or weekday, respectively, and that the recall is the first one. Then Similar formulae are used for the other daily-consumed foods and nutrients.
Finally, the energy-adjusted usual intakes and the HEI-2005 scores are then obtained as in Table 1, using the estimated usual intakes of the dietary components.
To find the joint distribution of usual intakes of the HEI-2005 scores, it is convenient to use Monte Carlo methods. Recall that w i is the sampling weight for individual i. Let B be a large number: we set B = 5000. Generate b = 1, . . . , B observations U bi = Normal(0, Σ u ) and then obtain T bi = (T biℓ ) 13 ℓ=1 by replacing U ij in their formulae by U bij . With appropriate sample weighting, the T bi can be used to estimate joint and marginal distributions. Thus, for example, consider the total HEI-2005 score, which is a deterministic function of the usual intakes, say, G( T i ). Its cumulative distribution function is estimated as Frequentist standard errors of derived quantities such as mean, median and quantiles can be estimated using the Balanced Repeated Replication (BRR) method [Wolter (1995)]; see Section 5.7 for details.
4. Comments on the approach to estimation. Our model (3.2)-(3.4) is a highly nonlinear, mixed effects model with many latent variables and nonlinear restrictions on the covariance matrix Σ ε . As seen in Section 3.4, we can estimate relevant distributions of usual intake in the population if we can estimate Σ u , Σ ε and β j for j = 1, . . . , 19. We have found that working within a pseudo-likelihood Bayesian paradigm is a convenient way to do this computation. We emphasize, however, that we are doing this only to get frequentist parameter estimates based on the well-known asymptotic equivalence of frequentist likelihood estimators and Bayesian posterior means, and especially the consistency of both [Lehmann and Casella (1998)]. We are specifically not doing Bayesian posterior inference, since valid Bayesian inference in a complex survey such as NHANES is an immensely challenging task, and because frequentist estimation and inference are the standard in the nutrition community. Kipnis et al. (2009) were able to get estimates of parameters separately for each food group using the nonlinear mixed effects program NLMIXED in SAS with sampling weights. While this gives estimates of β j for j = 1, . . . , 19, it only gives us parts of the covariance matrices Σ u and Σ ε , and not all the entries. Using the 2001-2004 NHANES data, we have verified that our estimates and the subset of the parameters that can be estimated by one food group at a time using NLMIXED are in close agreement, and that estimates of the distributions of usual intake and HEI-2005 component scores are also in close agreement. We expect this because of the rather large sample size in our data set.  have shown that even considering a single food group plus energy is a challenge for the NLMIXED procedure, both in time and in convergence, and using this method for the entire HEI-2005 constellation of dietary components is impossible.
Full technical details of the model fitting procedure are given in Appendices A.1-A.10.
Of course, our model has assumptions, for example, additivity and homoscedasticity on a transformed scale for observed and latent variables, normality of person-specific random effects and normality of day-to-day variability on the transformed scale. These assumptions are clearly not exactly correct, although our marginal model-checking suggests to us that they are mostly not disastrously wrong. Some reasons for this conclusion include the facts that we reproduce the marginal distributions of the components, that comparison with 24 h recalls shows differences that decrease when moving from one 24 h recall to two 24 h recalls, that q-q plots of the data are fairly satisfactory, etc. Thinking, as we do, of our work as a first step, and not a last step, it would be extremely interesting to make the model more general, for example, skew-normal, skew-t or Dirichlet process distributions after transformation, and possibly directly modeling heteroscedasticity. Such generalizations will require effort to implement, but will speak to the robustness of the results and would be a useful future step.

Empirical work.
5.1. Basic analysis. We analyzed data from the 2001-2004 National Health and Nutrition Examination Survey (NHANES) for children ages 2-8. The study sample consisted of 2638 children, among whom 1103 children have two 24 h recalls and the rest have only one. We used the dietary intake data to calculate the 12 HEI-2005 components plus energy. In addition, besides age, gender, race and interaction terms, two covariates were employed, along with an intercept. The first was a dummy variable indicating whether or not the recall was for a weekend day (Friday, Saturday or Sunday) because food intakes are known to differ systematically on weekends and weekdays. The second was a dummy variable indicating whether the 24 h recall was the first or second such recall, the idea being that there may be systematic differences attributable to the repeated administration of the instrument.

Contextual information.
When we ran our program based on the variables in Table 1, the results were disastrous. Mixing of the MCMC sampler was very poor, with long sojourns in different regions.
The reason for this failure to converge depends on the context of the dietary variables. For example, whole grains are a subset of total grains. Thus, if someone consumes any whole grains, then necessarily, with probability 1.0, that person also consumes total grains. Such a restriction cannot be handled by our model, because it would force one of the random effects U to equal infinity. A similar thing happens for energy. Calories coming from saturated fat are a subset of total calories, as are calories from SoFAAS, so there is a restriction that total calories must be greater than calories from saturated fat and also greater than calories from SoFAAS. Since the latter sum makes up a significant portion of calories, this restriction is not something that our model can handle well.
Luckily, there is an easy and natural context-based solution. Instead of using total grains in the model, we used grains that are not whole grains, that is, refined grains, thus decoupling whole grains and total grains, and removing the restriction mentioned above. Similarly, instead of using total fruit, we use fruit that is not whole fruits, that is, fruit juices. Additionally, instead of using total vegetables, we use total vegetables excluding dark green and orange vegetables and legumes. Finally, instead of total energy, we use total energy minus the sum of energy from saturated fat (11% of mean energy) and from SoFAAS (35% of mean energy). We recognize that there is overlap of energy from saturated fat and energy from solid fat, but this has no impact on our analysis since total energy has sources other than these two. An alternative, of course, would have been to simply use total energy minus energy from SoFAAS, This is sufficient to estimate the distributions of interest. If, for example, in the new data set T i1 represents usual intake of nonwhole fruits, and T i2 is usual intake of whole fruits, then the usual intake of total fruits is T i1 + T i2 . Similar remarks apply for total grains and total vegetables.
With these new variables, our model mixed well and gave reasonable looking answers that, as mentioned in Section 4, give similar results to other methods employed with smaller parts of the data set.

Estimation of the HEI-2005 scores.
In the Introduction we posed 4 questions to which answers had not been possible previously. The first open question concerned the distribution of the HEI total score. Along the way toward this, Table 2 presents the energy-adjusted distributions of the dietary components used in the HEI-2005. Table 3 presents the distributions of the HEI-2005 individual component scores and the total score, with a graphical view given in Figure 1. Table 3 presents the first estimates of the distribution of HEI-2005 scores for a vulnerable subgroup of the population, namely, children aged 2-8 years. A previous analysis of 2003-2004 NHANES data, looking separately at 2-5 year olds and 6-11 year olds, was limited to estimates of mean usual HEI-2005 scores [59.6 and 54.7, respectively; see Fungwe et al. (2009)]. The mean scores noted here are comparable to those and reinforce the notion that children's diets, on average, are far from ideal. However, this analysis provides a more complete picture of the state of US children's diets. By including the scores at various percentiles, we estimate that only 5% of children have a score of 69 or greater and another 10% have scores of 41 or lower. While not in the table, we also estimate that the 99th percentile is 74. This analysis  suggests that virtually all children in the US have suboptimal diets and that a sizeable fraction (10%) have alarmingly low scores (41 or lower.) We have also considered whether our multivariate model fitting procedure gives reasonable marginal answers. To check this, we note that it is possible to use the SAS procedure NLMIXED separately for each component to fit 16 S. ZHANG ET AL. For each component score, the first line = estimate from our model, while the second line is its BRR-estimated standard error. The total score is the sum of the individual scores. Here, "DOL" is dark green and orange vegetables and legumes. Also, "SoFAAS" is calories from solid fats, alcoholic beverages and added sugars. Further discussion of the size of the BRR-estimated standard errors is given in the supplementary material ]. procedure. The only case where there is a mild discrepancy is in the estimated variability of the energy-adjusted usual intake of oils, likely caused by the NLMIXED procedure itself, which has an estimated variance 9 times greater than our estimated variance. Of course, it is the distribution of the HEI-2005 total score that cannot be estimated by analysis of one component at a time.
There are other things that have not been computed previously that are simple by-products of our analysis. For example, the correlations among energy-adjusted usual intakes involving episodically consumed foods have not been estimated previously, but this is easy for us; see Table 4. The estimated correlation of −0.64 between energy-adjusted total fruit and energyadjusted SoFAAS, and the −0.47 correlation between DOL and SoFAAS are surprisingly high. 5.4. Component scores and other scores. As described in the Introduction, an open problem has been to estimate the correlation between the individual score on each dietary component and the scores of all other dietary components. In their Table 3,  consider this problem, but of course they did not have a model for usual energy adjusted intakes, and instead they used a single 24 h recall. In Table 5 we show the resulting correlations using (a) a single 24 h recall; (b) the mean of two 24 h recalls for those who have two 24 h recalls; and (c) our model for usual intake. The numbers for the former differ from that of  because we are considering here a different population than do they. A striking and not unexpected aspect of this table is that for those components with nontrivial correlations, the correlations all increase as one moves from a single 24 h 18 S. ZHANG ET AL.  Here TF = Total Fruits, WF = Whole Fruits, TV = Total Vegetables, WG = Whole Grains, TG = Total Grains, SatFat = Saturated Fat. Here, "DOL" is dark green and orange vegetables and legumes. Also, "SoFAAS" is calories from solid fats, alcoholic beverages and added sugars. The column labeled "Two 24 h" is the naive analysis that uses the mean of the two 24 h recalls, while the column labeled "First 24 h" is the naive analysis that uses the first 24 h recall. The column labeled "Model" is our analysis, and the column labeled "BRR s.e." is the estimated standard error of our estimates. Here, "DOL" is dark green and orange vegetables and legumes. Also, "SoFAAS" is calories from solid fats, alcoholic beverages and added sugars.
recall to the mean of two 24 h recalls and then finally to estimated usual intake. Thus, for example, the correlation between the HEI-2005 score for total fruit and its difference with the total score is 0.38 for a single 24 h recall, 0.44 for the mean of two 24 h recalls and then finally 0.62 for usual intake.

Distributions of intakes for subsets of HEI total scores.
A third open question is as follows: among those whose total HEI-2005 score is >50 or ≤50, what is the distribution of energy-adjusted usual intake of whole grains, whole fruits, dark green and orange vegetables and legumes (DOL) and calories from solid fats, alcoholic beverages and added sugars (SoFAAS)? This follows naturally from our method. Following (3.8), let G 1 ( T bi ) be energy adjusted usual intake and let G 2 ( T bi ) be the HEI total score. Then the distributions in question for when the total HEI-2005 score is >50 can be estimated as The results are provided in Table 6, with a graphical view in Figure 2. The results show that those who have poorer diets with usual HEI-2005 total score ≤ 50 are consistently eating poorer diets, that is, less whole fruits, less whole grains and less DOL, but higher SoFAAS. Here, "DOL" is dark green and orange vegetables and legumes. Also, "SoFAAS" is calories from solid fats, alcoholic beverages and added sugars. Units of measurement are given in Table 2. The solid lines are for those whose usual HEI-2005 total score is ≤50, that is, poorer diets, while the dashed lines are for those whose usual HEI-2005 total score is >50, that is, better diets.
5.6. Dietary consistency. We stated in the Introduction that it is interesting to understand the percentage of children whose usual intake HEI score exceeds the median HEI score on all 12 HEI components. Those median scores, say, (κ 1 , . . . , κ 12 ), are estimated in Table 3. If G j ( T bi ) is the HEI component score for episodically consumed food j, then following (3.8) the quantity in question can be estimated as We estimate that the percentage is 6%, woefully small. The percentage of children whose usual intake HEI score exceeds the median HEI score on all 12 HEI components is 0.24%. Figure 3 gives the estimated probabilities of exceeding the κ percentile on all 12 HEI components simultaneously, for κ = 1, 2, . . . , 99.  Tables 2 and 3. The BRR weights are only used in variance calculations. Once we have estimated some quantity, say, θ, from the sample using sample weight, we will need to compute the same quantity using, in succession, the 32 BRR weights. This will give us 32 estimates θ 1 , θ 2 , . . . , θ 32 . The BRR estimate for the variance of θ is (32 × 0.49) −1 32 p=1 ( θ p − θ) 2 . The 32 in the denominator is for the 32 different estimates from the 32 different sets of weights, and the 0.49 is the square of the perturbation factor used to construct the BRR weight sets [Wolter (1995)].
6. Further discussion of the analysis.
6.1. Never consumers. An aspect of the modeling that we have not discussed is the possibility that some people never, ever consume an episodically consumed dietary component. Our model does not allow for this, for general reasons and for reasons that are specific to our data analysis.
It is in principle possible to add an additional modeling step for nonconsumers, via fixed effects probit regression, but we do not think this is a practical issue in our case, for two reasons: • The first is that the HEI-2005 is based on 6 episodically consumed dietary components, namely, total fruit, whole fruit, whole grains, total vegetables, DOL and milk, the latter of which includes cheese, yogurt and soy beverages. None of these are "lifestyle adverse," unlike, say, alcohol. While 40% of the responses for whole fruits, for example, equal zero, the percentage of children who never eat any whole fruits at all is likely to be minuscule.
• Even if one disputes whether there are very few individuals who never consume one of the dietary components, then it necessarily follows that we have overestimated the HEI-2005 total scores, and, hence, the estimates of the proportion of individuals with alarmingly low HEI scores are deflated, and not inflated. The reason is that our model suggests everyone has a positive usual intake of the 6 episodically consumed dietary components.
Since the HEI-2005 score components are nondecreasing functions of usual intake of the episodically consumed dietary components, this would mean that we overestimate the HEI-2005 total score.

6.2.
Computing and data. Our programs were written in Matlab. The programs, along with the NHANES data we used, are available in the Annals of Applied Statistics online archive. Although a much smaller amount of computing effort yields similar results, using 70,000 MCMC steps with a burn-in of 20,000 takes approximately 10 hours on a Linux server.
We also estimated the Monte Carlo standard error which is defined by Flegal, Haran and Jones (2008) as σ g / √ n, where n is the total of iterations, and n = ab, where a is the number of blocks and b is the block size, and whereȲ g(X i ) for j = 1, . . . , a.
The batch means estimate of σ 2 g is The ratio of the Monte Carlo standard error to the estimated standard deviation of the estimated parameters averages 3.4% for Σ u and 1.7% for β.
Because of the public health importance of the problem, the National Cancer Institute has contracted for the creation of a SAS program that performs our analysis. It will allow any number of episodically and daily consumed dietary components. The first draft of this program, written independently in a different programming language, gives almost identical results to what we have obtained, at least suggesting that our results are not the product of a programming error.
7. Discussion. 7.1. Transformations. In Appendix A.12 we describe how we estimated the transformation parameters as a separate component-wise calculation. We have done some analyses where we simultaneously transform each component, and found very little difference with our results. However, the comp-uting time to implement this is extremely high, because of the fact that different transformations make data on different scales, so we have to compute the usual intakes at each step in the MCMC, and not just at the end. 7.2. What have we learned that is new. There are many important questions in dietary assessment that have not been able to be answered because of a lack of multivariate models for complex, zero-inflated data with measurement errors and a lack of ability to fit such multivariate models. Nutrients and foods are not consumed in isolation, but rather as part of a broader pattern of eating. There is reason to believe that these various dietary components interact with one another in their effect on health, sometimes working synergistically and sometimes in opposition. Nonetheless, simply characterizing various patterns of eating has presented enormous statistical challenge. Until now, descriptive statistics on the HEI-2005 have been limited to examination of either the total scores or only a single energy-adjusted component at a time. This has precluded characterization of various patterns of dietary quality as well as any subsequent analyses of how such patterns might relate to health.
This methodology presented in this paper presents a workable solution to these problems which has already proven valuable. In May 2010, just as we were submitting the paper, a White House Task Force on Childhood Obesity created a report. They had wanted to set a goal of all children having a total HEI score of 80 or more by 2030, but when they learned we estimated only 10% of the children ages 2-8 had a score of 66 or higher, they decided to set a more realistic target. The facility to estimate distributions of the multiple component scores simultaneously will be important in tracking progress toward that goal.

7.3.
In what other arenas will our work have impact? There are many other important problems where multivariate models such as ours will be important. One such problem arises when studying the relationship between multiple dietary components or dietary patterns and health outcomes. Traditionally, for cost reasons, large cohort studies have used a food frequency questionnaire (FFQ) to measure dietary intake, sometimes with a small calibration study including short-term measures such as 24 h recalls. However, there is a new web-based instrument called the Automated Self-administered 24-hour Dietary Recall (ASA24 TM ) (see http://riskfactor.cancer.gov/ tools/instruments/asa24), which has been proposed to replace or at least supplement the FFQ and which is currently undergoing extensive testing. The dietary data we will see then is what we have called Y ijk , that is, 24 h recall data. In order to correct relative risk estimates for the measurement error inherent in the ASA24 TM , regression calibration ] will almost certainly be the method of choice, as it is in most of nutritional 24 S. ZHANG ET AL.
epidemiology. This method attempts to produce an estimate of the regression of usual intake on the observed intakes, and then to use these estimates in Cox and logistic regression for the health outcome. In order to perform this regression, a multivariate measurement error model will be required, since the regression is on all the observed dietary intake components in the regression model measured by the ASA24 TM , and not on each individual component. Our methodology is easily extended to address this problem.

APPENDIX: DETAILS OF THE FITTING PROCEDURE
In this Appendix we give the full details of the model fitting procedure.
A.1. Notational convention. In our example, age was standardized to have mean 0.0 and variance 1.0, to improve numerical stability.
A.2. Prior distributions. Because the data were standardized, we used the following conventions: • The prior for all β j were normal with mean zero and variance 100.
• The prior for Σ u was exchangeable with diagonal entries all equal to 1.0 and correlations all equal to 0.50. There were 21 degrees of freedom in the inverse Wishart prior, that is, m u = 21. Thus, the prior is IW{(m u − 19 − 1)Σ u,prior , m u }. We experimented with this prior by using zero correlation, and the results were essentially unchanged. • The prior for r k is Uniform[−1, 1]. Set the initial value: r k = 0, k = 1, . . . , 5.
• For the rest of the nondiagonal v ij 's which could not be determined by the restrictions, we used Uniform[−3, 3] priors. Set the initial values to be 0.
The constraints on Σ ε are nonlinear, and our parameterization enforces them easily without having to have prior distributions for the original parameterization that satisfy the nonlinear constraints.
The key thing that makes things work well with the other components of the matrix V with Σ ε = V V T is that we have standardized the data as described in Appendix A.1. With this standardization, things become much nicer. For example, the variance of the ε's for energy is 19 j=1 v 2 19,j . However, since the sample variance for energy is standardized to equal 2.0, we simply just need to make priors for v 19,j be uniform on a modest range to have real flexibility.
• Finally, we then updated W ik by a single application of the updates given in Appendix A.9.
A.4. Complete data loglikelihood. Let J = 19. The complete data include the indicators of whether a food was consumed, the W variables and the random effect U variables. The loglikelihood of the complete data is (1 − r 2 q )} We used Gibbs sampling to update this complete data loglikelihood, the details for which are given in subsequent appendices. The weights w i are integers and are used here in a pseudo-likelihood fashion. One can also think of this as expanding each individual into w i individuals, each with the same observed data but different latent variables. For computational convenience, since we are only asking for a frequentist estimator and not doing full Bayesian inference, the latent variables in the process are generated once for each individual. Estimates of Σ u , Σ ε and β j for j = 1, . . . , J were computed as the means from the Gibbs samples. Once again, we emphasize that we are not doing a proper Bayesian analysis, but only using MCMC techniques to obtain a frequentist estimate, with uncertainty assessed using the frequentist BRR method.