Explanatory Variable

• The explanatory variable is GENHLTH, which is self-reported level of general health response to the following question.

Question: Would you say that in general your health is…

• The General Health variable is categorical variable that has been stored as an integer. In the study, nearly all categorical variables were stored in a numeric format.

• This categorical variable has an ordinal scale because its values of Excellent through Poor have an order to them. However, it is not an interval scale, so we cannot assume that the intervals between adjacent values are of equal length.

Response Variable

• The response variable is the Body Mass Index originally labeled as X_BMI5. The World Health Organization (WHO) describes the BMI as a “simple index of weight-for-height that is commonly used to classify underweight, overweight and obesity in adults.”

• It is a numeric and continuous variable (although the study only records it to 2 decimal places) that was multiplied by 100 and stored as an integer.

• The BRFSS calculated the respondents’ BMI by converting the reported heights and weights to meters and kilograms and then dividing the weight by the square of the height.

• The WHO lists the primary classification of the BMI as follows:

WHO: BMI

Generalizability:

The population of interest is adult U.S. residents. Since random sampling was used, the findings can be generalized to the population.

Weighting was used to account for residents without telephone access.

Causality:

Because this is an observation study and not an experimental one, no random assignment was used, and therefore, these data cannot be used to establish causal links between the variables of interest.

Entire Sample

options(scipen = 999)
kable(t(round(describe(bmi_gen_health$BMI), 4)))

By General Health Response

stats_by_grp <- describeBy(bmi_gen_health$BMI, bmi_gen_health$reported_general_health)

kable(sapply(stats_by_grp, function(x) rbind(x)))

Distribution of the General Health Responses

ggplot(bmi_gen_health, aes(x = reported_general_health)) + geom_bar()

Boxplots by Reported General Health

boxplot(bmi_gen_health$BMI ~ bmi_gen_health$reported_general_health)

BMI Histogram & Q-Q plot

The BMI distribution is right-skewed.

ggplot(bmi_gen_health, aes(x = BMI)) + geom_histogram()

qqnorm(bmi_gen_health$BMI)
qqline(bmi_gen_health$BMI)

Check Conditions

• Independence: Yes, the data are from a stratified random sample.

• Sample Size: Yes, the sample size is very large.

• Normal Distribution: The distribution is right-skewed, but this condition can be relaxed since we have a massive sample size.

95% Confidence interval for the entire sample

I used the T-distribution since the BMI distribution is right-skewed. However, the sample size and degrees of freedom are so large that it will approach the normal distribution.

n <- nrow(bmi_gen_health)
df <- n - 1
xbar <- mean(bmi_gen_health$BMI)
sign_level <- 0.05
s <- sd(bmi_gen_health$BMI)
tstar <- qt(sign_level/2, df, lower.tail = F)
se <- s/sqrt(n)
me <- se * tstar
c(xbar - me, xbar + me)

## [1] 28.02216 28.06319

We are 95% confident that the BMI population mean is between 28.02216 and 28.06319.

Check Conditions for ANOVA

• Independence: Yes, the data are from a random sample from less than 10% of the population. Also, the groups are independent of each other.

• Approximately normal: The Q-Q plots do not look very normal at all. However, since the sample size is so large, perhaps the normality requirement can be set aside.

op <- par(mfrow = c(2, 3))

for (i in levels(bmi_gen_health$reported_general_health)) {
    tmp <- subset(bmi_gen_health$BMI, bmi_gen_health$reported_general_health == 
        i)
    hist(tmp, xlab = "Reported General Health", main = i)
}
par(op)

op2 <- par(mfrow = c(2, 3))

for (i in levels(bmi_gen_health$reported_general_health)) {
    tmp <- subset(bmi_gen_health$BMI, bmi_gen_health$reported_general_health == 
        i)
    qqnorm(tmp, xlab = "Reported General Health", main = i)
    qqline(tmp)
}
par(op2)

# code citation:
# http://stackoverflow.com/questions/20781663/qqnorm-plotting-for-multiple-subsets

• Constant variance: From the 5 box plots above and the standard deviation statistics, we can see that the variability increases from the Excellent to Poor health groups. This may undermine the results of ANOVA.

t(t(sapply(stats_by_grp, function(x) x$sd)))

##                 [,1]
## 1 Excellent 5.247704
## 2 Very Good 5.797408
## 3 Good      6.761064
## 4 Fair      7.790429
## 5 Poor      8.620629

ANOVA Hypotheses

\(H_0:\) The average BMI is equal across all levels of self-reported general health.

\(H_A:\) The average BMI varies by the level of self-reported general health. At least one pair of group averages is more different than we would expect to see by chance.

ANOVA Analysis

The F-statistics is quite large at 5312, which yields a p-value near zero. Consequently, we would reject the null hypothesis.

aov_output <- aov(BMI ~ reported_general_health, bmi_gen_health)

summary(aov_output)

##                             Df   Sum Sq Mean Sq F value
## reported_general_health      4   893717  223429    5312
## Residuals               404025 16992204      42        
##                                      Pr(>F)    
## reported_general_health <0.0000000000000002 ***
## Residuals                                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Check conditions

• Independence: Yes, the data are from a random sample of less than 10% of the population. There is independence between groups.

• Sample Size: Yes, the sample size is very large.

• Normal Distribution: The distribution is right-skewed, but this condition can be relaxed since we have a massive sample size.

Hypotheses

For each comparison, the null hypotheis is that there is no difference between the means, and the alternative hypotheis is that there is a difference.

\(H_0: \mu_diff = 0\)

\(H_A: \mu_diff \neq 0\)

Analysis

For 9 out of the 10 comparisons, the p-value is less than the Bonferroni-corrected significance level of 0.005. Thus, we would reject the null hypothesis that there is no difference in means for these instances.

H0 <- 0

health_levels <- levels(bmi_gen_health$reported_general_health)
# number of groups
k <- length(health_levels)

K <- k * (k - 1)/2
alpha <- 0.05

# Bonferroni correction
alpha_star <- alpha/K

combos <- data.frame(t(combn(health_levels, 2)), stringsAsFactors = F)

results <- data.frame()

for (i in 1:nrow(combos)) {
    v1 <- subset(bmi_gen_health$BMI, bmi_gen_health$reported_general_health == 
        combos[i, 1])
    v2 <- subset(bmi_gen_health$BMI, bmi_gen_health$reported_general_health == 
        combos[i, 2])
    n1 <- length(v1)
    n2 <- length(v2)
    xbar1 <- mean(v1)
    xbar2 <- mean(v2)
    s1 <- sd(v1)
    s2 <- sd(v2)
    xbar_diff <- xbar1 - xbar2
    
    se_diff <- sqrt(s1^2/n1 + s2^2/n2)
    df_diff <- min(n1, n2) - 1
    
    t_score <- (xbar_diff - H0)/se_diff
    pvalue <- pt(-abs(t_score), df_diff) * 2
    
    results <- rbind(results, cbind(combos[i, ], xbar_diff, t_score, pvalue, 
        pvalue < alpha_star))
}

kable(results)