HW 02 - Data Wrangling

Hanke Zheng

Download and Read in the data

First download and then read in with data.table:fread()

# download individual/regional data from github
download.file("https://raw.githubusercontent.com/USCbiostats/data-science-data/master/01_chs/chs_individual.csv", "chs_individual.csv", method="libcurl", timeout = 60)
individual <- data.table::fread("chs_individual.csv")

download.file("https://raw.githubusercontent.com/USCbiostats/data-science-data/master/01_chs/chs_regional.csv", "chs_regional.csv", method="libcurl", timeout = 60)
regional <- data.table::fread("chs_regional.csv")

library(data.table)

## Warning: package 'data.table' was built under R version 4.0.2

library(dtplyr)

## Warning: package 'dtplyr' was built under R version 4.0.2

library(dplyr)

There are 1200 raws and 23 coloumns in the individual dataset; 12 raws and 27 variables in the regional dataset. ### Merge

## merge the two dataset by location
ind_region <- merge(x = individual, y = regional, by.x = "townname", by.y = "townname", all.x = TRUE, all.y = FALSE) 

count_n <- merge(x = individual, y = regional, by.x = "townname", by.y = "townname", all.x = TRUE, all.y = FALSE) %>%nrow()

The # of the merged dataset is the same as that of the individual dataset - no duplicates.

Manipulate the data

# check the missing values for bmi
message("Missing: ", ind_region[is.na(bmi), .N])

## Missing: 89

# impute data using the average within the variables “male” and “hispanic"
ind_region[, bmi_imp := fcoalesce(bmi, mean(bmi, na.rm = TRUE)),
    by = .(male, hispanic)]
ind_region$bmi[is.na(ind_region$bmi)] <- ind_region$bmi_imp

## Warning in ind_region$bmi[is.na(ind_region$bmi)] <- ind_region$bmi_imp: number
## of items to replace is not a multiple of replacement length

# check the missing values for fev
message("Missing: ", ind_region[is.na(fev), .N])

## Missing: 95

# impute data using the average within the variables “male” and “hispanic"
ind_region[, fev_imp := fcoalesce(fev, mean(fev, na.rm = TRUE)),
    by = .(male, hispanic)]
ind_region$fev[is.na(ind_region$fev)] <- ind_region$fev_imp

## Warning in ind_region$fev[is.na(ind_region$fev)] <- ind_region$fev_imp: number
## of items to replace is not a multiple of replacement length

# Create a new categorical variable named “obesity_level” using the BMI measurement (underweight BMI<14; normal BMI 14-22; overweight BMI 22-24; obese BMI>24). 
ind_region$obesity_level <- ifelse(ind_region$bmi>= 14 & ind_region$bmi<= 22, "normal", 
                         ifelse(ind_region$bmi<14, "underweight",
                                ifelse(ind_region$bmi > 24, "obese", "overweight")))
table(ind_region$obesity_level)

## 
##      normal       obese  overweight underweight 
##         959         109          95          37

range(ind_region$bmi)

## [1] 11.29640 41.26613

# check the missing values for binary variables 'smoke' and 'gasstove'
message("Missing: ", ind_region[is.na(smoke), .N])

## Missing: 40

message("Missing: ", ind_region[is.na(gasstove), .N])

## Missing: 33

# Create another categorical variable named “smoke_gas_exposure” that summarizes “Second Hand Smoke” and “Gas Stove.” The variable should have four categories in total.
ind_region$smoke_gas_exposure <- ifelse(ind_region$smoke== 1 & ind_region$gasstove == 1, "exposed to both", 
                         ifelse(ind_region$smoke==1 & ind_region$gasstove != 1, "exposed to smoke",
                        ifelse(ind_region$smoke !=1 & ind_region$gasstove == 1, "exposed to gas stove", "not exposed to both")))
table(ind_region$smoke_gas_exposure)

## 
##      exposed to both exposed to gas stove     exposed to smoke 
##                  151                  739                   36 
##  not exposed to both 
##                  214

Summary Table

#Create four summary tables showing the average (or proportion, if binary) and sd of “Forced expiratory volume in 1 second (ml)” and asthma indicator by town, sex, obesity level, and “smoke_gas_exposure.”
# convert to data frame
ind_region_df <- as.data.frame(ind_region)
# by town
ind_region_df %>%
  group_by(townname) %>%
  summarise(fev_ave = mean(fev),
            fev_sd = sd(fev),
            prop_asthma = mean(asthma, na.rm = TRUE),
            asthma_sd = sd(asthma, na.rm = TRUE)
  )

## `summarise()` ungrouping output (override with `.groups` argument)

## # A tibble: 12 x 5
##    townname      fev_ave fev_sd prop_asthma asthma_sd
##    <chr>           <dbl>  <dbl>       <dbl>     <dbl>
##  1 Alpine          2096.   299.       0.113     0.319
##  2 Atascadero      2083.   327.       0.255     0.438
##  3 Lake Elsinore   2039.   311.       0.126     0.334
##  4 Lake Gregory    2090.   322.       0.152     0.360
##  5 Lancaster       2012.   342.       0.165     0.373
##  6 Lompoc          2035.   355.       0.113     0.319
##  7 Long Beach      1982.   322.       0.135     0.344
##  8 Mira Loma       2004.   350.       0.158     0.367
##  9 Riverside       1993.   281.       0.11      0.314
## 10 San Dimas       2025.   320.       0.172     0.379
## 11 Santa Maria     2026.   328.       0.134     0.342
## 12 Upland          2034.   363.       0.121     0.328

# by sex
ind_region_df %>%
  group_by(male) %>%
  summarise(fev_ave = mean(fev),
            fev_sd = sd(fev),
            prop_asthma = mean(asthma, na.rm = TRUE),
            asthma_sd = sd(asthma, na.rm = TRUE)
  )

## `summarise()` ungrouping output (override with `.groups` argument)

## # A tibble: 2 x 5
##    male fev_ave fev_sd prop_asthma asthma_sd
##   <int>   <dbl>  <dbl>       <dbl>     <dbl>
## 1     0   1971.   329.       0.121     0.326
## 2     1   2101.   313.       0.173     0.378

# by obesity level
ind_region_df %>%
  group_by(obesity_level) %>%
  summarise(fev_ave = mean(fev),
            fev_sd = sd(fev),
            prop_asthma = mean(asthma, na.rm = TRUE),
            asthma_sd = sd(asthma, na.rm = TRUE)
  )

## `summarise()` ungrouping output (override with `.groups` argument)

## # A tibble: 4 x 5
##   obesity_level fev_ave fev_sd prop_asthma asthma_sd
##   <chr>           <dbl>  <dbl>       <dbl>     <dbl>
## 1 normal          2004.   310.      0.141      0.349
## 2 obese           2249.   328.      0.198      0.400
## 3 overweight      2221.   308.      0.163      0.371
## 4 underweight     1715.   307.      0.0811     0.277

# by smoke_gas_exposure
ind_region_df %>%
  group_by(smoke_gas_exposure) %>%
  summarise(fev_ave = mean(fev),
            fev_sd = sd(fev),
            prop_asthma = mean(asthma, na.rm = TRUE),
            asthma_sd = sd(asthma, na.rm = TRUE)
  )

## `summarise()` ungrouping output (override with `.groups` argument)

## # A tibble: 5 x 5
##   smoke_gas_exposure   fev_ave fev_sd prop_asthma asthma_sd
##   <chr>                  <dbl>  <dbl>       <dbl>     <dbl>
## 1 exposed to both        2023.   308.       0.130     0.338
## 2 exposed to gas stove   2032.   327.       0.148     0.355
## 3 exposed to smoke       2046.   371.       0.171     0.382
## 4 not exposed to both    2056.   335.       0.148     0.356
## 5 <NA>                   2018.   346.       0.149     0.360

4. Exploratory Data Analysis

1. What is the association between BMI and FEV (forced expiratory volume)?
1. What is the association between smoke and gas exposure and FEV?
1. What is the association between PM2.5 exposure and FEV?

Facet plot showing scatterplots with regression lines of BMI vs FEV by “townname”.

library(leaflet)

## Warning: package 'leaflet' was built under R version 4.0.2

library(tidyverse)

## Warning: package 'tidyverse' was built under R version 4.0.2

## ── Attaching packages ───────────────────────────────────────────────────────── tidyverse 1.3.0 ──

## ✓ ggplot2 3.3.1     ✓ purrr   0.3.4
## ✓ tibble  3.0.1     ✓ stringr 1.4.0
## ✓ tidyr   1.1.0     ✓ forcats 0.5.0
## ✓ readr   1.3.1

## ── Conflicts ──────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::between()   masks data.table::between()
## x dplyr::filter()    masks stats::filter()
## x dplyr::first()     masks data.table::first()
## x dplyr::lag()       masks stats::lag()
## x dplyr::last()      masks data.table::last()
## x purrr::transpose() masks data.table::transpose()

library(ggplot2)

ind_region_df %>%
  filter(!(townname %in% NA)) %>%
  ggplot(mapping=aes(x=bmi, y=fev, color=townname)) +
  geom_point()+
  stat_smooth(method=lm)+
  facet_wrap(~townname)

## `geom_smooth()` using formula 'y ~ x'

A positive association was observed between BMI and FEV based on the scatterplots.

#Stacked histograms of FEV by BMI category and FEV by smoke/gas exposure. Use different color schemes than the ggplot default.

ind_region_df %>%
  ggplot()+
  geom_histogram(mapping=aes(x=fev,fill=obesity_level))+
  # change the default palette
  scale_fill_brewer(palette = "PuOr")+
  labs(title="FEV by BMI category and by obesity level", x="FEV", y="count")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

  theme()

##  Named list()
##  - attr(*, "class")= chr [1:2] "theme" "gg"
##  - attr(*, "complete")= logi FALSE
##  - attr(*, "validate")= logi TRUE

ind_region_df %>%
  filter(!(smoke_gas_exposure %in% NA)) %>%
  ggplot()+
  geom_histogram(mapping=aes(x=fev,fill=smoke_gas_exposure))+
  # change the default palette
  scale_fill_brewer(palette = "Zissou")+
  labs(title="FEV by BMI category and by smoke/gas exposure", x="FEV", y="count")

## Warning in pal_name(palette, type): Unknown palette Zissou

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

  theme()

##  Named list()
##  - attr(*, "class")= chr [1:2] "theme" "gg"
##  - attr(*, "complete")= logi FALSE
##  - attr(*, "validate")= logi TRUE

FEV is normally distributed by obesity level and exposure status.

Barchart of BMI by smoke/gas exposure.

ind_region_df %>%
  filter(!(smoke_gas_exposure %in% NA)) %>%
  ggplot()+
  geom_bar(mapping=aes(x=smoke_gas_exposure,fill=obesity_level))+
  # change the default palette
  scale_fill_brewer(palette = "PuOr")+
  labs(title="BMI level by smoke/gas exposure", x="Exposure Status", y="count")

  theme()

##  Named list()
##  - attr(*, "class")= chr [1:2] "theme" "gg"
##  - attr(*, "complete")= logi FALSE
##  - attr(*, "validate")= logi TRUE

The majority of people in different exposure status are within normal BMI level.

#Statistical summary graphs of FEV by BMI and FEV by smoke/gas exposure category.

ind_region_df %>%
ggplot() +
  geom_boxplot(mapping=aes(y=fev, fill=obesity_level))

ind_region_df %>%
    filter(!(smoke_gas_exposure %in% NA)) %>%
    ggplot() +
    geom_boxplot(mapping=aes(y=fev, fill=smoke_gas_exposure))

#A leaflet map showing the concentrations of PM2.5 mass in each of the CHS communities.

library(leaflet)
leaflet(ind_region_df) %>%
  addProviderTiles('OpenStreetMap') %>%
  addCircles(lat=~lat, lng=~lon, color="green",opacity=1,
             fillOpacity=1, radius=~(pm25_mass*300))

The PM2.5 level is relatively high in LA areas in comparative to areas in the northern California.

Choose a visualization to examine whether PM2.5 mass is associated with FEV.

ind_region[,fev_ave :=mean(fev),by=townname]
ind_region %>%
  ggplot(mapping = aes(x=pm25_mass, y=fev_ave))+
  geom_point()+
  stat_smooth(method=lm)+
  labs(titles="FEV VS PM2.5 mass", x="PM2.5 mass", y="FEV")

## `geom_smooth()` using formula 'y ~ x'

Using the average FEV for each town, a negatvie association between PM2.5 and FEV was observed.

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

HW2.md

HW2.md

HW 02 - Data Wrangling

Download and Read in the data

Manipulate the data

Summary Table

4. Exploratory Data Analysis

Facet plot showing scatterplots with regression lines of BMI vs FEV by “townname”.

Barchart of BMI by smoke/gas exposure.

Choose a visualization to examine whether PM2.5 mass is associated with FEV.

Files

HW2.md

Latest commit

History

HW2.md

File metadata and controls

HW 02 - Data Wrangling

Download and Read in the data

Manipulate the data

Summary Table

4. Exploratory Data Analysis

Facet plot showing scatterplots with regression lines of BMI vs FEV by “townname”.

Barchart of BMI by smoke/gas exposure.

Choose a visualization to examine whether PM2.5 mass is associated with FEV.