Table Of Contents

Load packages

library(dplyr)
library(ggplot2)
library(statsr)

Load data

Make sure your data and R Markdown files are in the same directory. When loaded your data file will be called brfss2013. Delete this note when before you submit your work.

load("brfss2013.RData")

Part 1: Data

For this survey, the data is collected from a randomly selected adult in a household. But for the cellular telephone version of the survey, the data is collected from the people who are living in private residence or college housing.

This process clearly involves random sampling. Hence the result obtained from the data will be more generalizable

There is no random assignment of treatment to the random sampling, to take the survey. So the overall survey is not causal

Part 2: Research questions

Research quesion 1:

As the entire survey is about Behavioural risk factor, the first question is does this survey gives any data about one being healthy based on various input like food (fruits, vegetables, greens vegetables, orange vegetables) and physical exercises of any kind (based on how much time involved for it).

Healthy can be described based on depression level, blood pressure, diabetes and kidney disease.

So this involves variables
  1. food
  2. exercise
  3. health conditions

Research quesion 2:

Does income play a major role in having one or more personal health professionals and taking decision on whether to see doctor or not?

Variables Involved
  1. income
  2. has personal health professional
  3. medcost

Research quesion 3:

And with regards to mental health, does income play a role in it?

Variables involved
  1. menthlth
  2. income2

Part 3: Exploratory data analysis

Research quesion 1:

For Food, we’ll consider the variables fruit1 (cumulative frequency : 93.13%), fvbeans (cumulative frequency : 92.38%), fvgreen (cumulative frequency : 92.85%), fvorang (cumulative frequency : 92.57%) and vegetab1 (cumulative frequency : 92.09%) . With respect to all the cumulative frequencies (with proper values), our conclusion will be a generalized one.

Since the values are different as per day consumption, per week consumption and per month consumption, we will change it to a common base as no of fruits consumed in month as per below R code. So, We’ll create 5 new columns fruit_consumption, fvbeans_consumption, fvgreen_consumption, fvorang_consumption and vegetab1_consumption, which will contain the per month consumption.

brfss2013 <- brfss2013 %>% mutate(fruit_consumption = 
                                    ifelse(fruit1>=101 & fruit1<=199, (fruit1-100)*30,
                                           ifelse(fruit1>=201 & fruit1<=299, (fruit1-200)*4.5, 
                                                  ifelse(fruit1>=300 & fruit1<=399, fruit1-300, 
                                                         ifelse(is.na(fruit1),0,0)))
))
brfss2013 <- brfss2013 %>% mutate(fvbeans_consumption = 
                                    ifelse(fvbeans>=101 & fvbeans<=199, (fvbeans-100)*30,
                                           ifelse(fvbeans>=201 & fvbeans<=299, (fvbeans-200)*4.5, 
                                                  ifelse(fvbeans>=300 & fvbeans<=399, fvbeans-300, 
                                                         ifelse(is.na(fvbeans),0,0)))
))
brfss2013 <- brfss2013 %>% mutate(fvgreen_consumption = 
                                    ifelse(fvgreen>=101 & fvgreen<=199, (fvgreen-100)*30,
                                           ifelse(fvgreen>=201 & fvgreen<=299, (fvgreen-200)*4.5, 
                                                  ifelse(fvgreen>=300 & fvgreen<=399, fvgreen-300, 
                                                         ifelse(is.na(fvgreen),0,0)))
))
brfss2013 <- brfss2013 %>% mutate(fvorang_consumption = 
                                    ifelse(fvorang>=101 & fvorang<=199, (fvorang-100)*30,
                                           ifelse(fvorang>=201 & fvorang<=299, (fvorang-200)*4.5, 
                                                  ifelse(fvorang>=300 & fvorang<=399, fvorang-300, 
                                                         ifelse(is.na(fvorang),0,0)))
))
brfss2013 <- brfss2013 %>% mutate(vegetab1_consumption = 
                                    ifelse(vegetab1>=101 & vegetab1<=199, (vegetab1-100)*30,
                                           ifelse(vegetab1>=201 & vegetab1<=299, (vegetab1-200)*4.5, 
                                                  ifelse(vegetab1>=300 & vegetab1<=399, vegetab1-300, 
                                                         ifelse(is.na(vegetab1),0,0)))
))

As we are going to concentrate overall, we’ll have a new column overall_fvo_consumption to know the overall consumption.

brfss2013 <- brfss2013 %>% 
  mutate(overall_fvo_consumption = fruit_consumption+fvbeans_consumption+fvgreen_consumption+fvorang_consumption+vegetab1_consumption)

Now, we’ll create a new column exercise_time_in_mins to store the overall time spent in physical work out. We’ll make use of the variable exerhmm1 (cumulative frequency : 43.25%, so we’ll take it as they didn’t do exercise) for this.

brfss2013 <- brfss2013 %>% 
  mutate(exercise_time_in_mins = ifelse(!is.na(exerhmm1),(floor(exerhmm1/100)*60)+((exerhmm1/100)-floor(exerhmm1/100))*100,0))

Now that we have both food and exercise variable, we’ll work on health variables bphigh4 (cumulative frequency : 99.71%), addepev2 (cumulative frequency : 99.89%), chckidny (cumulative frequency : 99.95%) and diabete3 (cumulative frequency : 99.83%). We’ll create 4 new variables for these 4 categorical variables with values as 1/0 (1 - have disease, 2 - no disease).

brfss2013 <- brfss2013 %>% mutate(bp_level = ifelse(bphigh4 == "Yes" | bphigh4 == "Told borderline or pre-hypertensive", 1, 0))
brfss2013 <- brfss2013 %>% mutate(dep_level = ifelse(addepev2 == "Yes", 1, 0))
brfss2013 <- brfss2013 %>% mutate(kid_level = ifelse(chckidny == "Yes", 1, 0))
brfss2013 <- brfss2013 %>% mutate(dia_level = ifelse(diabete3 == "Yes" | diabete3 == "No, pre-diabetes or borderline diabetes", 1, 0))

We’ll create a common variable health_level to denote the overall health as the sum of all 4 above newly created variables. Min value denotes good health and bigger value denotes some health problem.

brfss2013 <- brfss2013 %>% mutate(health_level = bp_level+dep_level+kid_level+dia_level)

Let’s take a summary stats of the variables and save it in health_data

health_data <- brfss2013 %>% 
  filter(!is.na(overall_fvo_consumption), !is.na(health_level), !is.na(exercise_time_in_mins)) %>% 
  group_by(health_level) %>% 
  summarise(avg_food=mean(overall_fvo_consumption), avg_ex=mean(exercise_time_in_mins))

The content of health_data is

health_data
## # A tibble: 5 x 3
##   health_level avg_food avg_ex
##          <dbl>    <dbl>  <dbl>
## 1            0     93.9   50.8
## 2            1     87.9   43.8
## 3            2     83.2   36.1
## 4            3     80.2   28.7
## 5            4     78.3   24.5

2 observations made are

1. Let’s check the plot of Health vs Excise

ggplot(data = health_data, aes(x=health_level, y=avg_ex)) + geom_line()

Conclusion : With more time spent in physical exercise, health will be good

2. Let’s check the plot of Health vs Excise

ggplot(data = health_data, aes(x=health_level, y=avg_food)) + geom_line()

Conclusion : With more and more of fruits/veggies, health will be good

Research quesion 2:

Does income play a major role in having one or more personal health professionals and taking decision on whether to see doctor or not?

Also to conclude if the people of money conscious or health conscious or medical expense is too high

Let’s take a look at the comparison of income and if people has 1 or more personal medical professionals. From below, it is obvious that income plays a level for people who want to have 1 or more health care professionals. Also due to income, more people with less income has no professional health professional compare to the people with more income.

personal_doc_table = table(brfss2013$income2, brfss2013$persdoc2)
mosaicplot(personal_doc_table, "Income vs # Personal Heath Professionals", xlab="Income", ylab="Personal Health Professionals")

Let’s take a look at the comparison of income and medical cost which made people to not see doctor. From below, we can infer that even the people with more money hesitates to see doctors. So it means that the medical expenses are high

mosaic_table = table(brfss2013$income2, brfss2013$medcost)
mosaicplot(mosaic_table, "Income vs Seeing Doctor", xlab="Income", ylab="Seeing Doctor")

Let’s look at the together comparison. There is a small conclusion that even with people with more income and has more that 1 medical professional, still they opt not to see the doctor

mosaic_table = table(brfss2013$income2, brfss2013$persdoc2, brfss2013$medcost)
mosaicplot(mosaic_table, "Income Category vs (persdoc2 vs medcost) ", xlab="Income category", ylab="persdoc2")

We need to look deeper. So we’ll create a new data health_income_persdoc2_medcost with only 3 deciding variables income2 (cumulative frequency : 85.48%), persdoc2 (cumulative frequency : 99.87%) and medcost (cumulative frequency : 99.94%). New data will be grouped by income, persdoc2 and then medcost. We’ll have one more column probability which holds the probability on whether a particular group (with income2 and persdoc2) will choose to go to see doctor or not.

health_income_persdoc2_medcost <- brfss2013 %>% 
  filter(!is.na(medcost), !is.na(income2), !is.na(persdoc2)) %>% 
  group_by(income2, persdoc2, medcost) %>% 
  summarize(n=n())

health_income_persdoc2_medcost <- health_income_persdoc2_medcost %>% mutate(probability=0)
health_income_persdoc2_medcost
## # A tibble: 48 x 5
## # Groups:   income2, persdoc2 [24]
##    income2           persdoc2      medcost     n probability
##    <fct>             <fct>         <fct>   <int>       <dbl>
##  1 Less than $10,000 Yes, only one Yes      3855           0
##  2 Less than $10,000 Yes, only one No      12639           0
##  3 Less than $10,000 More than one Yes       576           0
##  4 Less than $10,000 More than one No       1717           0
##  5 Less than $10,000 No            Yes      2820           0
##  6 Less than $10,000 No            No       3566           0
##  7 Less than $15,000 Yes, only one Yes      3883           0
##  8 Less than $15,000 Yes, only one No      14615           0
##  9 Less than $15,000 More than one Yes       576           0
## 10 Less than $15,000 More than one No       2060           0
## # ... with 38 more rows

We have to create one more data health_income_persdoc2 to know the no of people in the income2 and persdoc2 specific group

health_income_persdoc2 <- brfss2013 %>% 
  filter(!is.na(medcost), !is.na(income2), !is.na(persdoc2)) %>% 
  group_by(income2, persdoc2) %>% 
  summarize(n=n())
health_income_persdoc2
## # A tibble: 24 x 3
## # Groups:   income2 [?]
##    income2           persdoc2          n
##    <fct>             <fct>         <int>
##  1 Less than $10,000 Yes, only one 16494
##  2 Less than $10,000 More than one  2293
##  3 Less than $10,000 No             6386
##  4 Less than $15,000 Yes, only one 18498
##  5 Less than $15,000 More than one  2636
##  6 Less than $15,000 No             5446
##  7 Less than $20,000 Yes, only one 23800
##  8 Less than $20,000 More than one  3011
##  9 Less than $20,000 No             7767
## 10 Less than $25,000 Yes, only one 29611
## # ... with 14 more rows

Now we’ll calculate the probability and save it in health_income_persdoc2_medcost based on the group count calculated & saved in health_income_persdoc2

for(i in 1:dim(health_income_persdoc2_medcost)[1]){
    
  inc     <- data.frame(lapply(health_income_persdoc2_medcost[i,1], as.character), stringsAsFactors=FALSE)
  persdoc <- data.frame(lapply(health_income_persdoc2_medcost[i,2], as.character), stringsAsFactors=FALSE)
    
  health_income_persdoc2_medcost <- health_income_persdoc2_medcost %>% 
          mutate(probability = ifelse((income2 %in% inc & persdoc2 %in% persdoc), n/(health_income_persdoc2 %>% filter(income2 %in% inc, persdoc2 %in% persdoc))$n,probability))
}
health_income_persdoc2_medcost
## # A tibble: 48 x 5
## # Groups:   income2, persdoc2 [24]
##    income2           persdoc2      medcost     n probability
##    <fct>             <fct>         <fct>   <int>       <dbl>
##  1 Less than $10,000 Yes, only one Yes      3855       0.234
##  2 Less than $10,000 Yes, only one No      12639       0.766
##  3 Less than $10,000 More than one Yes       576       0.251
##  4 Less than $10,000 More than one No       1717       0.749
##  5 Less than $10,000 No            Yes      2820       0.442
##  6 Less than $10,000 No            No       3566       0.558
##  7 Less than $15,000 Yes, only one Yes      3883       0.210
##  8 Less than $15,000 Yes, only one No      14615       0.790
##  9 Less than $15,000 More than one Yes       576       0.219
## 10 Less than $15,000 More than one No       2060       0.781
## # ... with 38 more rows

Now that we have all required data, let see some stats

health_income_persdoc2_medcost %>%filter(medcost == "Yes", persdoc2 == "More than one") %>% arrange(desc(probability))
## # A tibble: 8 x 5
## # Groups:   income2, persdoc2 [8]
##   income2           persdoc2      medcost     n probability
##   <fct>             <fct>         <fct>   <int>       <dbl>
## 1 Less than $10,000 More than one Yes       576      0.251 
## 2 Less than $15,000 More than one Yes       576      0.219 
## 3 Less than $20,000 More than one Yes       538      0.179 
## 4 Less than $25,000 More than one Yes       614      0.175 
## 5 Less than $35,000 More than one Yes       524      0.132 
## 6 Less than $50,000 More than one Yes       510      0.102 
## 7 Less than $75,000 More than one Yes       350      0.0694
## 8 $75,000 or more   More than one Yes       394      0.0433
health_income_persdoc2_medcost %>%filter(medcost == "Yes", persdoc2 == "Yes, only one") %>% arrange(desc(probability))
## # A tibble: 8 x 5
## # Groups:   income2, persdoc2 [8]
##   income2           persdoc2      medcost     n probability
##   <fct>             <fct>         <fct>   <int>       <dbl>
## 1 Less than $10,000 Yes, only one Yes      3855      0.234 
## 2 Less than $15,000 Yes, only one Yes      3883      0.210 
## 3 Less than $20,000 Yes, only one Yes      4502      0.189 
## 4 Less than $25,000 Yes, only one Yes      4788      0.162 
## 5 Less than $35,000 Yes, only one Yes      4361      0.122 
## 6 Less than $50,000 Yes, only one Yes      4171      0.0893
## 7 Less than $75,000 Yes, only one Yes      3184      0.0621
## 8 $75,000 or more   Yes, only one Yes      2898      0.0310
health_income_persdoc2_medcost %>%filter(medcost == "Yes", persdoc2 == "No") %>% arrange(desc(probability))
## # A tibble: 8 x 5
## # Groups:   income2, persdoc2 [8]
##   income2           persdoc2 medcost     n probability
##   <fct>             <fct>    <fct>   <int>       <dbl>
## 1 Less than $10,000 No       Yes      2820      0.442 
## 2 Less than $15,000 No       Yes      2344      0.430 
## 3 Less than $20,000 No       Yes      2995      0.386 
## 4 Less than $25,000 No       Yes      2851      0.341 
## 5 Less than $35,000 No       Yes      2225      0.252 
## 6 Less than $50,000 No       Yes      1788      0.187 
## 7 Less than $75,000 No       Yes      1057      0.122 
## 8 $75,000 or more   No       Yes       758      0.0577

In all the above cases, even people with less income preferred to have more than 1 personal health professional. But still they don’t want to see the doctor due to the cost of the medical. This trend becomes less evident with people with more income.

Conclusion : Even people with low income, they are all health conscious. But the medical cost is high, which causes people to not see the doctor. It is evident from 44.2% of low level income people didn’t go to see doctor even though they don’t have any personal health professional.

Research quesion 3:

As a final question, with regards to mentally healthy, does income play a role in it?

Let’s take a guess, people with moderate level of income should be mentally healthy compared to low level and high level incomes. Because low level income people will be stressed with the money and high level income people worrying about making more money.

For this, we’ll consider IQR to have a robust factor to measure the happiness. We’ll create a new data mental_health_data with the variables income2 (cumulative frequency : 85.48%) and menthlth (cumulative frequency : 98.24%).

mental_health_data <- brfss2013 %>% 
  filter(!is.na(income2), !is.na(menthlth)) %>% 
  group_by(income2) %>% 
  summarize(iqr = IQR(menthlth)) %>% 
  arrange(desc(iqr))

The data we obtained is

mental_health_data
## # A tibble: 8 x 2
##   income2             iqr
##   <fct>             <dbl>
## 1 Less than $10,000    15
## 2 Less than $15,000    10
## 3 Less than $20,000     5
## 4 Less than $25,000     3
## 5 Less than $35,000     2
## 6 Less than $50,000     2
## 7 Less than $75,000     1
## 8 $75,000 or more       1

Plot :

ggplot(data=mental_health_data, aes(x=income2, y=iqr)) + geom_bar(stat="identity")

Conclusion :So to be mentally healthy you need not have much income