Behavioral Risk Factor Surveillance System 2013 Exploratory Data Analysis
In this project, we carry out exploratory analysis of the BRFSS-2013 data set by setting out research questions, and then exploring relationship between identified variables to answer those questions. To know more about BRFSS and the dataset, visit this link.
The project was completed as a part of Duke University's 'Introduction to Probability and Data' online course on Coursera, the first of the Statistics with R Specialization.
Setup
Load packages
1library(ggplot2)2library(dplyr)Load data
1load("brfss2013.RData")The Data
The BRFSS-2013 dataset was sampled from the non-institutionalised adult population (i.e. 18 years and older) residing in the US. The data was collected through landline and cellular-telephone based surveys.
Disproportionate stratified sampling, which is more efficient than simple random sampling, was used for the landline sample (source). The cellular sample was generated from randomly selected respondents, with an equal probability of selection.
As random sampling was used for both data collection methods, the data for the sample is generalizable to the population. On the other hand, as this is an observational study, it won't be possible to make causal inferences from the data.
Research questions
Research question 1:
Are non-smoking heavy drinkers, generally healthier than regular smokers, who are not heavy drinkers?
While researching this, we're trying to explore the impact of consuming alcohol vs smoking tobacco on a person's health and see which is worse.
Research question 2:
Do people who sleep fewer hours than average person, also have more than days with poor mental health?
Research has suggested that inadequate sleep has a negative effect on a person's overall health. Here we try to determine if it also has a negative effect on their mental health.
Research question 3:
Are people who have completed higher levels of education, more likely to consume fruits and vegetables once or more in a day?
We might assume that educated people live a healthier lifestyle i.e. exercising or eating nutritious food. We'll try and figure out if that's the case here by comparing education levels with fruit and vegetable consumption.
Exploratory data analysis
Research question 1:
Are non-smoking heavy drinkers, generally healthier than regular smokers, who are not heavy drinkers?
We'll be using the following variables for this question:
- genhlth: Respondent's health, in general
- _rfsmok3: Is the respondent a current smoker?
- _rfdrhv4: Is the respondent a heavy drinker?
Type of the variables we're dealing with:
1str(select(brfss2013,genhlth,X_rfsmok3,X_rfdrhv4))## 'data.frame': 491775 obs. of 3 variables:## $ genhlth : Factor w/ 5 levels "Excellent","Very good",..: 4 3 3 2 3 2 4 3 1 3 ...## $ X_rfsmok3: Factor w/ 2 levels "No","Yes": 1 1 2 1 1 1 1 2 1 1 ...## $ X_rfdrhv4: Factor w/ 2 levels "No","Yes": 1 1 2 1 1 1 1 1 1 1 ...All of the above are categorical variable. General health of a person is defined in 5 levels, while a person is or isn't a heavy drinker or a smoker.
To begin, let's check out our selected variables individually.
genhlth: General Health
1total_obs <- nrow(brfss2013)2
3
4brfss2013 %>%5 group_by(genhlth) %>%6 summarise(count=n(),percentage=n()*100/total_obs)7
8## # A tibble: 6 × 39## genhlth count percentage10## <fctr> <int> <dbl>11## 1 Excellent 85482 17.382339512## 2 Very good 159076 32.347313313## 3 Good 150555 30.614610314## 4 Fair 66726 13.568400215## 5 Poor 27951 5.683696816## 6 NA 1985 0.403639917
18ggplot(brfss2013, aes(x=genhlth)) + geom_bar() + ggtitle('General Health of Respondents') + xlab('General Health') + theme_bw()
Around 80% of the respondents in our dataset are in good health or better, and most of the people have 'Very good' health. There are some missing (NA) values too which we'll deal with later as they don't make much sense with our analysis.
_rfsmok3: Currently a smoker?
According to the codebook, respondents who have replied 'Yes', now smoke every day or some days; while those who replied 'No' have either never smoked in their lifetimes or don't smoke now.
1brfss2013 %>%2 group_by(X_rfsmok3) %>%3 summarise(count=n(),percentage=n()*100/total_obs)## # A tibble: 3 × 3## X_rfsmok3 count percentage## <fctr> <int> <dbl>## 1 No 399786 81.294494## 2 Yes 76654 15.587210## 3 NA 15335 3.118296
ggplot(brfss2013, aes(x=X_rfsmok3)) + geom_bar() + ggtitle('Smoking Status of Respondents') + xlab('Currently a smoker?')+ theme_bw()
More than 81% of the respondents are not current smokers, though they might have smoked earlier in their lifetimes.
_rfdrhv4: Heavy drinker?
The heavy drinker variable is defined as adult men having more than two drinks per day and adult women having more than one drink per day).
1brfss2013 %>%2 group_by(X_rfdrhv4) %>%3 summarise(count=n(),percentage=n()*100/total_obs)## # A tibble: 3 × 3## X_rfdrhv4 count percentage## <fctr> <int> <dbl>## 1 No 442359 89.951502## 2 Yes 25533 5.192009## 3 NA 23883 4.856489
ggplot(brfss2013, aes(x=X_rfdrhv4)) + geom_bar() + ggtitle('Drinking Habits of Respondents') + xlab('Heavy Drinker?') +theme_bw()
Only about 5% of the respondends in our dataset are heavy drinkers.
Now to answer our original question, we can make things a bit easier for ourselves by creating a new categorical variable to categorise a person as: 'Smoker', 'Heavy Drinker', 'Both' or 'None.'
1brfss2013 <- brfss2013 %>%2 mutate(smoke_alc = ifelse(X_rfdrhv4 == 'Yes',3 ifelse(X_rfsmok3 == 'Yes','Both','Heavy Drinker'),4 ifelse(X_rfsmok3 == 'Yes','Current Smoker','None')))Let's check out the distribution of our new variable:
1brfss2013 %>%2 group_by(smoke_alc) %>%3 summarise(count=n(),percentage=n()*100/total_obs)## # A tibble: 5 × 3## smoke_alc count percentage## <chr> <int> <dbl>## 1 Both 8144 1.656042## 2 Current Smoker 66000 13.420772## 3 Heavy Drinker 17269 3.511565## 4 None 374377 76.127701## 5 <NA> 25985 5.283920
ggplot(brfss2013,aes(x=smoke_alc)) + geom_bar() + ggtitle('Drinking and Smoking Habits of Respondents') + xlab('Drinker or Smoker?') +theme_bw()
About 76% of the respondents don't smoke or drink heavily. Around 13.4% are current smokers, and about 3.5% drink heavily. We'll be focusing on the last two.
A good way to represent the counts of two categorical variables is a contingency table.
1rq1_table <- table(brfss2013$smoke_alc,brfss2013$genhlth)2
3 rq1_table## ## Excellent Very good Good Fair Poor## Both 998 2428 2957 1273 446## Current Smoker 6637 17160 22372 12828 6741## Heavy Drinker 4140 6729 4648 1356 349## None 69056 125370 112123 47522 18884It's a little hard to look at the number and quickly understand what proportions of Current Smokers or Heavy Drinkers have the better health. So we'll calculate those proportions, and make sure that the rows sum to 1. So we're calculating the proportions of health across drinker or smokers.
1prop.table(rq1_table,1)## ## Excellent Very good Good Fair Poor## Both 0.12317946 0.29967909 0.36497161 0.15712170 0.05504814## Current Smoker 0.10096139 0.26103623 0.34032067 0.19513828 0.10254343## Heavy Drinker 0.24039020 0.39072117 0.26988735 0.07873650 0.02026478## None 0.18515907 0.33615316 0.30063412 0.12742020 0.05063345Now we have a sense of what's going on. Let's visualize the table through a mosaic plot.
1mosaicplot(prop.table(rq1_table,1),main='Drinking and/or Smoking vs General Health', xlab='Drinking and/or Smoking status', ylab='General Health')
Looking at the summary statistics and the visualization, we can see that, compared to the 'Current Smoker' category, there is a higher proportion of 'Heavy Drinkers' with 'Excellent' or 'Very good' health. Even though there are proportionally more smokers with 'Good' health, heavy drinkers have lower 'Fair' or 'Poor' health, something we can consider below par.
Hence, it looks like smokers have poorer health than heavy drinkers.
Research question 2:
Do people who sleep fewer hours than average person, also have more than days with poor mental health?
For this, we have to look at the relationship between the variables:
- sleptim1: On average, the hours of sleep a person gets in a 24-hour period.
- menthlth: Out of 30, number of days the mental health of a person wasn't good.
Checking out the type of variables that we're dealing with:
1str(select(brfss2013,sleptim1,menthlth))## 'data.frame': 491775 obs. of 2 variables:## $ sleptim1: int NA 6 9 8 6 8 7 6 8 8 ...## $ menthlth: int 29 0 2 0 2 0 15 0 0 0 ...Both of the above variables are continuous integers, but we can also think of time slept as a categorical variable, and then calculate the average number of days with poor mental health for a person who gets that much amount of sleep, to answer our question. We'll do that in a bit.
sleptim1
Taking a look at how the sleptim1 variable is distributed.
1ggplot(brfss2013,aes(x=sleptim1)) + geom_bar()## Warning: Removed 7387 rows containing non-finite values (stat_count).
The plot extends to more than 400 on the x-axis, which is a bit suspicious. We'll have to check for unrealistic values for sleptim1 in our dataset. Let's filter for observations where the time slept is above the 24 hour period.
1brfss2013 %>%2 filter(sleptim1>24) %>%3 select(sleptim1)## sleptim1## 1 103## 2 450As suspected, there are two unrealistic values there, and we'll have to filter for them when doing our analyis or visualisation. One thing we can do to avoid repeated filters for this section is to make a new clean dataframe.
1rq2_brfss2013 <- brfss2013 %>%2 filter(sleptim1 <= 24)Also, since there are only 24 (25 if you count 0) possible values for the amount of time slept, we can consider this variable as a factor. Attempting the plot again:
1ggplot(rq2_brfss2013,aes(x=as.factor(sleptim1))) + geom_bar() + ggtitle('Amount of Sleep of Respondents') + xlab('Hours slept') + theme_bw()
Better. It looks like most people get 6-8 hours of sleep. What is the average hours of sleep for our data?
1rq2_brfss2013 %>%2 summarise(avg_sleep = mean(sleptim1))## avg_sleep## 1 7.050986So we can consider 7 hrs to be the optimum/average amount of sleep for our population.
menthlth
1ggplot(rq2_brfss2013, aes(x=menthlth)) + geom_histogram()## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 7771 rows containing non-finite values (stat_bin).
Looks like we might have the same problem here as with the previous variable. Removing impossible outliers:
1rq2_brfss2013 <- rq2_brfss2013 %>%2 filter(menthlth <= 30)3
4 ggplot(rq2_brfss2013, aes(x=menthlth)) + geom_histogram() + ggtitle('Mental Health of Respondents') + xlab('Number of days with poor mental health (out of 30)') + theme_bw()## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
To answer our question, we can look at the average mental health of groups of people who get different hours of sleep.
1slep_ment <- rq2_brfss2013 %>%2 group_by(hours_slept = as.factor(sleptim1)) %>%3 summarise(avg_poor_mental = mean(menthlth), count=n()) 4
5 slep_ment## # A tibble: 24 × 3## hours_slept avg_poor_mental count## <fctr> <dbl> <int>## 1 1 11.651376 218## 2 2 13.286550 1026## 3 3 12.265064 3369## 4 4 9.772228 13834## 5 5 6.280284 32699## 6 6 3.933364 104493## 7 7 2.206032 140794## 8 8 2.217405 138934## 9 9 2.739771 23414## 10 10 4.686248 11831## # ... with 14 more rowsVisualizing the data:
1ggplot(slep_ment, aes(x=hours_slept,y=avg_poor_mental)) + geom_bar(stat='identity') + ggtitle('Do people who sleep irregularly have poor mental health?') + xlab('Hours Slept') + ylab('No. of days with poor mental health (out of 30)') + theme_bw()
Looking at the summary statistics and the bar graph, we can see that people who get around 6-9 hours of sleep per day have considerably lower number of days with poor mental health, though it's even better if you get 7 or 8 hours of sleep.
People with lesser, or even more than average hours of sleep have more days of poor mental health. The observation with 23 hours of sleep and 30 days of poor mental health is caused by outliers with scarce data, as are others, possibly.
We can check this as well:
1slep_ment %>%2 filter(as.integer(hours_slept) > 12)## # A tibble: 12 × 3## hours_slept avg_poor_mental count## <fctr> <dbl> <int>## 1 13 9.448454 194## 2 14 8.567198 439## 3 15 8.424855 346## 4 16 6.719547 353## 5 17 8.428571 35## 6 18 10.318750 160## 7 19 10.000000 12## 8 20 8.548387 62## 9 21 7.500000 2## 10 22 8.200000 10## 11 23 30.000000 3## 12 24 7.218750 32So it does seem that there is a relationship between inadequate sleep and mental health, though we cannot be absolutely sure if one directly causes the other. But, people who sleep lower or more than average, are also those who suffer from poorer mental health than those who sleep adequately.
Research question 3:
Are people who have completed higher levels of education, more likely to consume fruits and vegetables once or mor in a day?
For this we'll be using the following variables:
- _educag: Computed level of education completed.
- _frtlt1: Consume fruit 1 or times per day.
- _veglt1: Consume vegetables 1 or times per day.
1str(select(brfss2013,X_educag,X_frtlt1,X_veglt1))## 'data.frame': 491775 obs. of 3 variables:## $ X_educag: Factor w/ 4 levels "Did not graduate high school",..: 4 3 4 2 4 4 2 3 4 2 ...## $ X_frtlt1: Factor w/ 2 levels "Consumed fruit one or more times per day",..: 1 2 2 2 2 1 1 2 1 2 ...## $ X_veglt1: Factor w/ 2 levels "Consumed vegetables one or more times per day",..: 2 1 1 1 1 1 1 1 1 1 ...All three of these are categorical variables, so we can answer our question through a contingency table.
Before proceeding, we'll change the name of the factor levels, as they make the outplot and visualizations look better. For this purpose, we;ll assume that college and technical school the same thing in _educag variable.
1levels(brfss2013$X_frtlt1) <- c('Once or more a day','Less than once a day')2 levels(brfss2013$X_veglt1) <- c('Once or more a day','Less than once a day')3
4 levels(brfss2013$X_educag)[3] <- c('Attended college')5 levels(brfss2013$X_educag)[4] <- c('Graduated college')Beginning with an examination of the variables on their own.
_educag
1brfss2013 %>%2 group_by(X_educag) %>%3 summarise(count=n(), percentage=n()*100/total_obs)## # A tibble: 5 × 3## X_educag count percentage## <fctr> <int> <dbl>## 1 Did not graduate high school 42213 8.5838036## 2 Graduated high school 142968 29.0718316## 3 Attended college 134196 27.2880891## 4 Graduated college 170118 34.5926491## 5 NA 2280 0.46362671ggplot(brfss2013, aes(x=X_educag)) + geom_bar() + ggtitle('Education Level of Respondents') + xlab('Completed Education Level') + theme_bw() + theme(axis.text.x = element_text(angle = 30, hjust = 1))
Around 90% of the respondents have graduated high school or higher, while 8.5% did not graduate high school. In the context of our question, we expect that as we go up the completed level of education, we'll see higher proportions of people who consume fruits and vegetables. We'll come back to that later.
_frtlt1
1brfss2013 %>%2 group_by(X_frtlt1) %>%3 summarise(count=n(), percentage=n()*100/total_obs)## # A tibble: 3 × 3## X_frtlt1 count percentage## <fctr> <int> <dbl>## 1 Once or more a day 291729 59.321641## 2 Less than once a day 171343 34.841747## 3 NA 28703 5.8366121ggplot(brfss2013, aes(x=X_frtlt1)) + geom_bar() + ggtitle('Fruit Consumption by Respondents') + xlab('Fruit Consumption') + theme_bw() + theme(axis.text.x = element_text(angle = 30, hjust = 1))
More people in our dataset, around 59% consume fruits one or more times per day than those who don't, around 34.8%.
_veglt1
1brfss2013 %>%2 group_by(X_veglt1) %>%3 summarise(count=n(), percentage=n()*100/total_obs)## # A tibble: 3 × 3## X_veglt1 count percentage## <fctr> <int> <dbl>## 1 Once or more a day 359834 73.170454## 2 Less than once a day 101777 20.695847## 3 NA 30164 6.1336991ggplot(brfss2013, aes(x=X_veglt1)) + geom_bar() + ggtitle('Vegetable Consumption by Respondents') + xlab('Vegetable Consumption') + theme_bw() + theme(axis.text.x = element_text(angle = 30, hjust = 1))
Around 73% people in our dataset consume vegetables one or more times per day, while around 20.7% don't. We can also see that there are more people that consume veggies (one or more times per day), than there are who consume fruits.
Answering our question, we'll build a 2-way contingency table to count frequencies of completed education level with both fruit and vegetable consumption.
1ct_f <- table(brfss2013$X_educag,brfss2013$X_frtlt1)2
3 prop.table(ct_f,1)## ## Once or more a day Less than once a day## Did not graduate high school 0.5352677 0.4647323## Graduated high school 0.5734143 0.4265857## Attended college 0.6220397 0.3779603## Graduated college 0.7052274 0.29477261mosaicplot(prop.table(ct_f,1), main='Completed Education vs Fruit Consumption', xlab='Completed Education Level', ylab='Fruit Consumption')
In both the proportional frequency table, and the mosaic, we can see an increase in the proportion of people who consume fruits, as we move to increased levels of completed education, vs those who don't. Let's see if this also holds true for vegetable consumption.
1ct_v <- table(brfss2013$X_educag,brfss2013$X_veglt1)2
3 prop.table(ct_v,1)## ## Once or more a day Less than once a day ## Did not graduate high school 0.6435618 0.3564382 ## Graduated high school 0.7123442 0.2876558 ## Attended college 0.7909306 0.2090694 ## Graduated college 0.8585728 0.14142721mosaicplot(prop.table(ct_v,1), main='Completed Education vs Vegetable Consumption', xlab='Completed Education Level', ylab='Vegetable Consumption')
There was already a higher proportion of people who consumed vegetables than those who consumed fruits, but the increase in proportions, depending on completed education level is still evident here.
So, we can say that people with higher education levels are more likely to adopt healthy eating habits.