**Win-Vector Blog » R**, and kindly contributed to R-bloggers)

The following article is getting quite a *lot* of press right now: David Just and Brian Wansink (2015). Fast Food, Soft Drink, and Candy Intake is Unrelated to Body Mass Index for 95% of American Adults. Obesity Science & Practice, forthcoming (upcoming in a new pay for placement journal). Obviously it is a popular contrary position (some coverage: here, here, and here).

I thought I would take a peek to learn about the statistical methodology (see here for some commentary). I would say the kindest thing you can say about the paper is: its problems are not statistical.

At this time the authors don’t seem to have supplied their data preparation or analysis scripts and the paper “isn’t published yet” (though they have had time for a press release), so we have to rely on their pre-print. Read on for excerpts from the work itself (with commentary).

After excluding the clinically underweight and morbidly obese, consumption of fast food, soft drinks or candy was not positively correlated with measures of BMI.

(Eliminate enough outcome variation and there is no variation to measure/explain.)

We restrict our sample to adults, defined as age 18 or older, who completed two 24-hour dietary recall surveys.

(It plausibly takes more than two days of measurements to get a good image of long term eating habits. Also most “food regulation”, a topic these authors have written on, is targeted at children. So for a useful public policy analysis it would have been nice to leave them in.)

We focus on eating episode rather than amount eaten because it is less subject to recall bias.

(Breaking the actual relation between eating and health, by leaving out amount. Also some effective diets advise more sittings of much smaller portions.)

We compare average eating episodes within food and across BMI categories.

(I am guessing this means they are modeling BMI category code instead of the BMI number. There are only about 3 BMI category codes left after “excluding the clinically underweight and morbidly obese.” Again eliminate variation in the measured outcome, and nothing will correlate to it.)

Missing data were omitted from the analysis …

(Just dropping missing data is not likely to work with interview data, unless you truly believe censoring is completely independent of health, diet, and health/diet interactions.)

Likewise, those with normal BMIs consume an average of 1.1 salty snacks over two days, while overweight, obese, and morbidly obese consume an average 0.9, 1.0, and 0.9 salty snacks, respectively.

(Uh, I thought we were “excluding the clinically underweight and morbidly obese.” I guess this is a different analysis. But here is a statistical issue: it really doesn’t look like the independent variable (“salty snacks”) is varying. So you are not going to be able to see if it drives an outcome. And since there isn’t a complete methods section I really wonder if the analysis is really looking at the claimed underlying data, or just looking at aggregate values.)

From: Table 1. Average Instances of Consumption in 48 Hours of Various Food Items, Sorted by BMI

(I’m not a statistician, but a negative p-value? Maybe that is some variation of z? But the weird values are not just in one column. Is all this just off one ANOVA table? Also, Why not try a linear regression on BMI score using non-grouped data, or a logistic regression on BMI category?)

Also when the input (or “independent”) variables are not known to be independent of each other ANOVA is variable order dependent! Usually this is handled by experiment design- but in this case we are observing eating patterns, not assigning them.

Some R code showing the effect is given below. Notice all of the x’s have the same relation to y, but the ANOVA analysis assigns effect in variable order. It does not make any sense to say “x1 is significant, but x10 is not” as the F-scores are not about each variable in isolation.

```
set.seed(6326)
d <- data.frame(y=rnorm(100))
for(i in 1:10) {
d[[paste('x',i,sep='')]] <- d$y + rnorm(nrow(d))
}
anova(lm(y~x1+x2+x3+x4+x5+x6+x7+x9+x10,data=d))
```

```
## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## x1 1 70.643 70.643 640.0173 < 2.2e-16 ***
## x2 1 22.647 22.647 205.1824 < 2.2e-16 ***
## x3 1 5.285 5.285 47.8821 6.425e-10 ***
## x4 1 6.588 6.588 59.6906 1.491e-11 ***
## x5 1 2.382 2.382 21.5771 1.155e-05 ***
## x6 1 3.027 3.027 27.4269 1.063e-06 ***
## x7 1 0.494 0.494 4.4757 0.03714 *
## x9 1 1.914 1.914 17.3441 7.137e-05 ***
## x10 1 0.376 0.376 3.4048 0.06830 .
## Residuals 90 9.934 0.110
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```