Unhealthy Days and Exercise

Research Topic

The level of obesity in the United States has a been a major concern in the media. And the level of mental health in the United States is starting to be a topic of debate. The research topic for this analysis is going to try and see what determines the amount of unhealthy days a person has per month. The independent variable is going to be the amount of access a person has to exercise opportunities. The population in this analysis are all the counties in all the states within the United States. The hypothesis for this analysis is that the higher the amount of exercise opportunities a person has the lower the amount of unhealthy days.

The Data

This data was taken from Social Explorer under the Health Data for 2016. The categories that were chosen were “Quality of Life” and “Diet and Exercise.” The Quality of Life category had “physically unhealthy days per month” and “mentally unhealthy days per month” which were combined together for this analysis to become the variable “Unhealthyday.” The exercise variable came from the category “Diet and Exercise.” Only certain variables were selected for this analysis from the “Diet and Exercise” category. The variable for States is going to be “StateNumber.” The exercise variable is measured by the “Percentage of population with adequate access to locations for physical activity.”

library(nlme) 
library(dplyr)
library(magrittr)
library(tidyr)
library(haven)
library(lmerTest)
library(ggplot2)
library(texreg)
library(readr)
health <- read_csv("/Users/paulkim/Downloads/CountyHealth.csv")
health$exercised <- ifelse(health$Exercise < 50, 
c(0), c(1)) 
head(health)

This R chunk is simply installing and calling the necessary packages. Also the data is read into R for analysis. The “Exercise” variable is turned into two different categories, if the variable for exercise is over 50 then the variable is equal to “1” and if the variable for exercise is less than 50 than the variable is equal to “0.” The head function is used to view the data and to see if the data was imported into R correctly.

length(unique(health$StateNumber))

[1] 51

The length function is showing the total amount of States in this analysis. There are going to be 51 states in this analysis because the District of Columbia is included in the “StateNumber” variable.

health %>% 
  group_by(StateNumber) %>% 
  summarise(n_county = n())

This R chunk is displaying how many counties are in each state.

Ecological Analysis

states <- health %>% 
  group_by(StateNumber) %>% 
  summarise(mean_p = mean(UnhealthyDays, na.rm = TRUE), mean_s = mean(exercised, na.rm = TRUE))
head(states)

The above code is demonstrating the ecological analysis of the data by taking the average of the “UnhealthyDays” variable and the average of the “exercised” variable. In this ecological analysis, the counties are ignored and the data is only looking at the state level analysis. So this analysis is looking at the state level amount of unhealthy day by the proportion of people who have access to exercise opportunities above and below 50 percent. This analysis is based off of the 51 different states.

ecoreg <- lm(mean_p ~ mean_s, data = states)
summary(ecoreg)

Call:
lm(formula = mean_p ~ mean_s, data = states)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.5069 -0.6917  0.1328  0.6436  1.9154 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   8.7152     0.5351  16.287   <2e-16 ***
mean_s       -1.9015     0.7133  -2.666   0.0104 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.017 on 49 degrees of freedom
Multiple R-squared:  0.1267,    Adjusted R-squared:  0.1088 
F-statistic: 7.107 on 1 and 49 DF,  p-value: 0.01037

The intercept is describing states with no exercise opportunities. The mean_s value is showing that states with more than 50% exercise opportunities will have -1.9015 unhealthy days per month than compared to those states with less than 50% exercise opportunity. This analysis is flawed because of an ecological fallacy. This analysis took the averages of different states and not the county level so assumptions cannot be made on the county level only on the state level.

Complete Pooling Method

cpooling <- lm(UnhealthyDays ~ exercised, data = health)
summary(cpooling)

Call:
lm(formula = UnhealthyDays ~ exercised, data = health)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.1521 -0.9521 -0.0117  0.8883  4.2883 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  7.85211    0.04068  193.01   <2e-16 ***
exercised   -0.54039    0.04997  -10.81   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.3 on 3027 degrees of freedom
  (74 observations deleted due to missingness)
Multiple R-squared:  0.0372,    Adjusted R-squared:  0.03689 
F-statistic:   117 on 1 and 3027 DF,  p-value: < 2.2e-16

The intercept here is showing that counties that have less than 50% exercise opportunity will have 7.85211 unhealthy days. This analysis is also showing that counties that have more than 50% exercise opportunity have -0.54039 unhealthy days than compared to those with less than 50% exercise opportunity. This means that as the amount of exercise opportunity increases the amount of unhealthy days decreases. The problem with the complete pooling analysis is that it ignores that the different counties come from different states and just assumes that each county comes from the same environment. This method does not look at the individual counties as having any weight.

No-pooling Model

The no-pooling model is giving out 51 different regression models for each state in the data set. Random Intercept Model

dcoef <- health %>% 
    group_by(StateNumber) %>% 
    do(mod = lm(UnhealthyDays ~ exercised, data = .))
coef <- dcoef %>% do(data.frame(intc = coef(.$mod)[1]))
ggplot(coef, aes(x = intc)) + geom_histogram()

The intercepts are showing that there is variation in the amount of unhealthy days per month between the different states. The intercepts are based off of the states with less than 50% exercise opportunity. The amount of unhealthy days ranges from about 5.5 to 10.

Random Slope Model

dcoef <- health %>% 
    group_by(StateNumber) %>% 
    do(mod = lm(UnhealthyDays ~ exercised, data = .))
coef <- dcoef %>% do(data.frame(exercised = coef(.$mod)[2]))
ggplot(coef, aes(x = exercised)) + geom_histogram()

The slope is showing the difference in unhealthy days based off the less than and greater than 50% exercise opportunity. The slope is also showing variation in the amount of unhealthy days per month between the different counties based on the percentage of exercise opportunity. The slopes here are showing the states that have over 50% exercise opportunity and how this will affect the amount of unhealthy days people will have. Many of the slopes here are showing 0 and negative numbers. This means that the amount of unhealthy days a state has will decrease as the amount of exercise opportunities increases. But the effect that the amount of exercise opportunities have on individual states are different due to the different environments of the varying states.

Partial Pooling

Random Intercept Model

m1_lme <- lme(UnhealthyDays ~ exercised, data = health, random = ~1|StateNumber, method = "ML", na.action = na.omit)
summary(m1_lme)

Linear mixed-effects model fit by maximum likelihood
 Data: health 
       AIC      BIC    logLik
  6571.307 6595.371 -3281.654

Random effects:
 Formula: ~1 | StateNumber
        (Intercept)  Residual
StdDev:    1.036961 0.6881355

Fixed effects: UnhealthyDays ~ exercised 
                Value  Std.Error   DF   t-value p-value
(Intercept)  7.618359 0.14831887 2977  51.36473       0
exercised   -0.374522 0.02832672 2977 -13.22152       0
 Correlation: 
          (Intr)
exercised -0.137

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-3.77381940 -0.62136187 -0.06369155  0.51851411  6.78511412 

Number of Observations: 3029
Number of Groups: 51 

The random intercept model here is showing random intercepts for the different states. The first portion of the lme function is estimating fixed effect parameters for this model and then the random effect is being determined by the “StateNumber.” The analysis here is making the “StateNumber” intercepts follow a normal distribution of 0 with a standard deviation of 1.036961. The fixed effect analysis is showing that states with less than 50% exercise opportunity will have 7.618359 unhealthy days. And that states with more than 50% exercise opportunity will have 0.374522 less unhealthy days compared to those states with less than 50% exercise opportunity.

Random Slope Model

m2_lme <- lme(UnhealthyDays ~ exercised, data = health, random = ~ exercised|StateNumber, method = "ML", na.action = na.omit)
summary(m2_lme)

Linear mixed-effects model fit by maximum likelihood
 Data: health 
       AIC      BIC    logLik
  6549.412 6585.508 -3268.706

Random effects:
 Formula: ~exercised | StateNumber
 Structure: General positive-definite, Log-Cholesky parametrization
            StdDev    Corr  
(Intercept) 1.0673804 (Intr)
exercised   0.2180573 -0.285
Residual    0.6813568       

Fixed effects: UnhealthyDays ~ exercised 
                Value  Std.Error   DF  t-value p-value
(Intercept)  7.627097 0.15376284 2977 49.60299       0
exercised   -0.395558 0.04671568 2977 -8.46734       0
 Correlation: 
          (Intr)
exercised -0.32 

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-3.76389208 -0.61906124 -0.06509881  0.49691916  6.65643924 

Number of Observations: 3029
Number of Groups: 51 

The fixed effects give the common trend shared by all of the 51 states in this observation. The 7.627097 intercept is subject to state level variation. So what a typical state looks like is based off of the fixed effect numbers. The average amount of unhealthy days for states that have less than 50% exercise opportunity is 7.627097 and state with more than 50% exercise opportunity tend to have 0.395558 less unhealthy days than those with less than 50% exercise opportunity.

AIC(ecoreg, cpooling, m1_lme, m2_lme)

models are not all fitted to the same number of observations

This R chunk is helping to select the model that best fits this analysis. According to these comparisons, model “ecoreg” is the best fitting model because it has the lower AIC value. But as mentioned earlier, the Ecological Analysis is not a good method to follow because of the ecological fallacy. The ecological fallacy means for this analysis that assumptions cannot be made on the county level only on the state level. Model “m2_lme” is the Random Slope Model. Model “m2_lme” is the next best fitting model because it has the second lowest AIC.

Findings

The Ecological Analysis, the Complete Pooling Method, the Random Intercept Model and the Random Slope Model are all showing that counties and states with less than 50% exercise opportunity will always have some amount of unhealthy days but as the amount of exercise opportunity goes above 50% the amount of unhealthy days will decrease. The increase in exercise opportunity was shown to decrease the amount of unhealthy days in all analyses conducted. The affect that exercise opportunity has on the amount of unhealthy days varies from method to method. The Random Slope Method is showing the intercept to be 7.627097 and states with more than 50% exercise opportunity have 0.395558 less unhealthy days than those with less than 50% exercise opportunity.

MerTools

library(lme4)
library(merTools)
library(readr)
health <- read_csv("/Users/paulkim/Downloads/CountyHealth.csv")

Parsed with column specification:
cols(
  County = col_character(),
  CountyState = col_character(),
  Column1 = col_character(),
  StateNumber = col_integer(),
  UnhealthyDays = col_double(),
  HealthCareCosts = col_double(),
  LimitedHealthyFood = col_double(),
  Exercise = col_double(),
  Obese = col_double(),
  PhysicallyInactive = col_double(),
  ChildrenFreeLunch = col_double()
)

head(health)

This R chunk is reading in the data and checking if the data is read in correctly.

n1 <- lmer(UnhealthyDays ~ Exercise + (1|Exercise), data=health)
fastdisp(n1)

lme4::lmer(formula = UnhealthyDays ~ Exercise + (1 | Exercise), 
    data = health)
            coef.est coef.se
(Intercept)  8.21     0.06  
Exercise    -0.01     0.00  

Error terms:
 Groups   Name        Std.Dev.
 Exercise (Intercept) 0.08    
 Residual             1.29    
---
number of obs: 3029, groups: Exercise, 3003
AIC = 10170.1

reEx <- REsim(n1)

Error: class(merMod) %in% c("lmerMod", "glmerMod", "blmerMod", "bglmerMod") is not TRUE

The linear mixed effect model is showing that counties with no exercise opportunity will have 8.21 unealthy days. But for every percentage increase in exercise opportunity, there is 0.01 less unhealthy days. The REsim function is creating random effect values to be used in a simulation.

p1 <- plotREsim(reEx)
p1

This graph is showing that there is little variance in the effect that exercise has on the amount of unhealthy days. There does not seem to be any major outliers. The effect range is staying close to 0 even as the amount of exercise varies.