index.Rmd

---
title: "Open Case Studies: Exploring global patterns of obesity across rural and urban regions"
css: style.css
output:
  html_document:
    includes:
       in_header: GA_Script.Rhtml
    self_contained: yes
    code_download: yes
    highlight: tango
    number_sections: no
    theme: cosmo
    toc: yes
    toc_float: yes
  pdf_document:
    toc: yes
  word_document:
    toc: yes
---

<style>
#TOC {
  background: url("https://opencasestudies.github.io/img/icon-bahi.png");
  background-size: contain;
  padding-top: 240px !important;
  background-repeat: no-repeat;
}
</style>

<!-- Open all links in new tab-->  
<base target="_blank"/>  
<div id="google_translate_element"></div>

<script type="text/javascript" src='//translate.google.com/translate_a/element.js?cb=googleTranslateElementInit'></script>

<script type="text/javascript">
function googleTranslateElementInit() {
  new google.translate.TranslateElement({pageLanguage: 'en'}, 'google_translate_element');
}
</script>



```{r setup, include=FALSE}
knitr::opts_chunk$set(include = TRUE, comment = NA, echo = TRUE,
                      message = FALSE, warning = FALSE, cache = FALSE,
                      fig.align = "center", out.width = '90%')
library(here)
library(knitr)
library(magrittr)
remotes::install_github("benmarwick/wordcountaddin", type = "source", dependencies = TRUE)
remotes::install_github("alistaire47/read.so")
library(wordcountaddin)
library(read.so)

rmarkdown:::perf_timer_reset_all()
rmarkdown:::perf_timer_start("render")
```

 
#### {.outline }

```{r, echo = FALSE}
knitr::include_graphics(here::here("img", "main_plot.png"))
```

####

#### {.disclaimer_block}

**Disclaimer**: The purpose of the [Open Case Studies](https://opencasestudies.github.io){target="_blank"} project is **to demonstrate the use of various data science methods, tools, and software in the context of messy, real-world data**. A given case study does not cover all aspects of the research process, is not claiming to be the most appropriate way to analyze a given dataset, and should not be used in the context of making policy decisions without external consultation from scientific experts. 

####

#### {.license_block}

This work is licensed under the Creative Commons Attribution-NonCommercial 3.0 [(CC BY-NC 3.0)](https://creativecommons.org/licenses/by-nc/3.0/us/){target="_blank"}  United States License.

####

#### {.reference_block}

To cite this case study please use:

Wright, Carrie and Meng, Qier and Jager, Leah and Taub, Margaret and Hicks, Stephanie. (2020). https://github.com/opencasestudies/ocs-bp-rural-and-urban-obesity. Exploring global patterns of obesity across rural and urban regions (Version v1.0.0).

####

To access the GitHub Repository for this case study see here: https://github.com/opencasestudies/ocs-bp-rural-and-urban-obesity.

You may also access and download the data using our `OCSdata` package. To learn more about this package including examples, see this [link](https://github.com/opencasestudies/OCSdata). Here is how you would install this package:

```{r, eval=FALSE}
install.packages("OCSdata")
```

This case study is part of a series of public health case studies for the [Bloomberg American Health  Initiative](https://americanhealth.jhu.edu/open-case-studies).

***

The total reading time for this case study is calculated via [koRpus](https://github.com/unDocUMeantIt/koRpus) and shown below: 

```{r, echo=FALSE}
readtable = text_stats("index.Rmd") # producing reading time markdown table
readtime = read.so::read.md(readtable) %>% dplyr::select(Method, koRpus) %>% # reading table into dataframe, selecting relevant factors
  dplyr::filter(Method == "Reading time") %>% # dropping unnecessary rows
  dplyr::mutate(koRpus = paste(round(as.numeric(stringr::str_split(koRpus, " ")[[1]][1])), "minutes")) %>% # rounding reading time estimate
  dplyr::mutate(Method = "koRpus") %>% dplyr::relocate(koRpus, .before = Method) %>% dplyr::rename(`Reading Time` = koRpus) # reorganizing table
knitr::kable(readtime, format="markdown")
```

***

**Readability Score: **

A readability index estimates the reading difficulty level of a particular text. Flesch-Kincaid, FORCAST, and SMOG are three common readability indices that were calculated for this case study via [koRpus](https://github.com/unDocUMeantIt/koRpus). These indices provide an estimation of the minimum reading level required to comprehend this case study by grade and age. 

```{r, echo=FALSE}
rt = wordcountaddin::readability("index.Rmd", quiet=TRUE) # producing readability markdown table
df = read.so::read.md(rt) %>% dplyr::select(index, grade, age) %>%  # reading table into dataframe, selecting relevant factors
  tidyr::drop_na() %>% dplyr::mutate(grade = round(as.numeric(grade)), # dropping rows with missing values, rounding age and grade columns
                                     age = round(as.numeric(age))
                                     )
knitr::kable(df, format="markdown")
```

***

Please help us by filling out our survey.


<div style="display: flex; justify-content: center;"><iframe src="https://docs.google.com/forms/d/e/1FAIpQLSfpN4FN3KELqBNEgf2Atpi7Wy7Nqy2beSkFQINL7Y5sAMV5_w/viewform?embedded=true" width="1200" height="700" frameborder="0" marginheight="0" marginwidth="0">Loading…</iframe></div>


## **Motivation**

Body mass index ([BMI](https://www.cdc.gov/healthyweight/assessing/bmi/adult_bmi/index.html){target="_blank"}) is often used as a proxy for **adiposity** (the condition of having excess body fat) and is measured as an individual's weight in kilograms (kg) or pounds (lbs) divided by the individual's height in meters ($m^2$) squared.

There are different categories of weight defined using different BMI thresholds as defined by the [World Health Organization](http://www.euro.who.int/en/health-topics/disease-prevention/nutrition/a-healthy-lifestyle/body-mass-index-bmi){target="_blank"} (WHO), the [Centers for Disease Control and Prevention](https://www.cdc.gov/obesity/adult/defining.html){target="_blank"}, and the [National Institutes of Health](https://www.nhlbi.nih.gov/sites/default/files/media/docs/obesity-evidence-review.pdf){target="_blank"}. 
For example, these are the categories defined by the WHO: 

```{r,  echo=FALSE, out.width='100%'}
knitr::include_graphics(here::here("img","obesity_classes.jpeg"))
```

###### [[source](https://www.ncbi.nlm.nih.gov/books/NBK535456/bin/bmi_WHO.jpg){target="_blank"}]

The following chart can help you identify your estimate of BMI:
```{r,  echo=FALSE, out.width='100%'}
knitr::include_graphics(here::here("img","BMIchart.jpeg"))
```

###### [[source](http://diet.mayoclinic.org/diet/eat/what-is-your-bmi?xid=nl_MayoClinicDiet_20160426){target="_blank"}]


Recently, an [article](https://www.nature.com/articles/s41586-019-1171-x.pdf){target="_blank"} published in [*Nature*](https://www.nature.com){target="_blank"} evaluated and compared the BMI of populations in rural and urban communities around the world: 

```{r,  echo=FALSE, out.width='100%'}
knitr::include_graphics(here::here("img", "paper.png"))
```

#### {.reference_block}
NCD Risk Factor Collaboration (NCD-RisC). Rising rural body-mass index is the main driver of the global obesity epidemic in adults. *Nature* **569**, 260–264 (2019). 

####

The article challenged the widely-held view that increased urbanization was one of the major reasons for increased global obesity rates. 
This view came about because many countries around the world have shown increased urbanization levels in parallel with increased obesity rates. 
In this article, however, the NCD-RisC argued that this might not be the case and that in fact for most regions around the world, BMI measurements are increasing in rural populations just as much, if not more so, than urban populations. 
Furthermore, this study suggested that obesity has particularly increased in female populations in rural regions: 

> "We noted a persistently higher rural BMI, especially for women."

In this case study, we will evaluate the data reported in this article to explore regional and gender specific differences in the obesity rates around the world in 1985 and 2017. 
Most importantly we will test if there is a difference in obesity rates between rural and urban communities and if there has been a change in obesity over time for these regions, particularly for women. 
To do this we will test if there is a difference between the average BMI for each group.


## **Main Questions**
***

#### {.main_question_block}

<b><u> Our main questions are: </u></b>

1) Is there a difference between rural and urban BMI estimates around the world? In particular, what does this difference look like for women?
2) How have BMI estimates changed from 1985 to 2017? In particular, what does this change over time look like for women?
3) How do different countries vary in terms of estimates of BMI? In particular, how does the United States compare to the rest of the world?

####

## **Learning Objectives** 
***

The skills, methods, and concepts that students will be familiar with by the end of this case study are:

<u>**Data science Learning Objectives:**</u>

1. Importing data from a PDF (`pdftools`)  
2. Subsetting and filtering data (`dplyr`)  
3. Working with character strings (`stringr`)  
4. Reshaping data into different formats (`tidyr`)  
5. Applying functions to all columns of a tibble (`purrr`)  
6. Creating data visualizations (`ggplot2`) with labels (`ggrepel`)  
7. Combining multiple plots (`cowplot` and `patchwork`)  

<u>**Statistical Learning Objectives:**</u> 

1. Familiarity with the use of Quantile-Quantile plots to assess normality   
2. Define and understand the utility of alpha and the p value  
3. Describe the difference between nonparametric and parametric tests  
4. Be able to identify paired data  
5. Implementation of a [paired $t$-test](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5579465/){target="_blank"}  
6. Interpretation of a [paired $t$-test](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5579465/){target="_blank"}  
7. Implementation of a [Wilcoxon signed-rank test](http://www.biostathandbook.com/wilcoxonsignedrank.html){target="_blank"}
8.Interpretation of a [Wilcoxon signed-rank test](http://www.biostathandbook.com/wilcoxonsignedrank.html){target="_blank"}
9. Understanding of the need for multiple testing correction  


In this case study, we will walk you through importing data from a pdf, cleaning, wrangling and visualizing the data, and <b> comparing two groups </b>.  

We will use well-established and commonly used packages, including `stringr`, `tidyr`, `dplyr`, `purrr`, and `ggplot2` from the [`tidyverse`](https://www.tidyverse.org/){target="_blank"}. 
Specifically, the tidyverse is

> "an opinionated collection of R packages designed for data science. All packages share an underlying philosophy and common APIs".

Another way of putting it is that it's a set of packages that are useful specifically for data manipulation, exploration, and visualization with a common philosophy.
While some students may be familiar with previous packages from the R programming language, these packages make data science in R especially efficient.

```{r, echo = FALSE, out.width = "20%"}
knitr::include_graphics("https://tidyverse.tidyverse.org/logo.png")
```

We will begin by loading the packages that we will need with a brief explanation:

```{r}
library(here)
library(pdftools)
library(stringr)
library(readr)
library(dplyr)
library(tibble)
library(magrittr)
library(glue)
library(purrr)
library(tidyr)
library(ggplot2)
library(ggrepel)
library(cowplot)
library(patchwork)
library(OCSdata)
```


<u>**Packages used in this case study:** </u>

 Package   | Use in this case study                                                                         
---------- |-------------
[here](https://github.com/jennybc/here_here){target="_blank"}       | to easily load and save data with relative paths
[pdftools](https://cran.r-project.org/web/packages/pdftools/pdftools.pdf){target="_blank"}   | to read a text from pdf into R   
[stringr](https://stringr.tidyverse.org/articles/stringr.html){target="_blank"}    | to manipulate the text data
[readr](https://readr.tidyverse.org/){target="_blank"}      | to manipulate the text data within the pdf into individual lines  
[dplyr](https://dplyr.tidyverse.org/){target="_blank"}      | to arrange/filter/select subsets of the data 
[tibble](https://tibble.tidyverse.org/){target="_blank"}     | to create data objects that we can manipulate with `dplyr`/`stringr`/`tidyr`/`purrr`
[magrittr](https://magrittr.tidyverse.org/articles/magrittr.html){target="_blank"}   | to use the `%<>%` piping operator
[glue](https://www.tidyverse.org/blog/2017/10/glue-1.2.0/){target="_blank"}  | to paste or combine character strings and data together
[purrr](https://purrr.tidyverse.org/){target="_blank"}      | to perform functions on all columns of a tibble
[tidyr](https://tidyr.tidyverse.org/){target="_blank"}      | to convert data from 'wide' to 'long' format
[ggplot2](https://ggplot2.tidyverse.org/)    | to make visualizations with multiple layers
[ggrepel](https://cran.r-project.org/web/packages/ggrepel/vignettes/ggrepel.html){target="_blank"}    | to allow labels in figures not to overlap
[cowplot](https://cran.r-project.org/web/packages/cowplot/vignettes/introduction.html){target="_blank"} and [patchwork](https://github.com/thomasp85/patchwork){target="_blank"} | to allow plots to be combined 
[OCSdata](https://github.com/opencasestudies/OCSdata){target="_blank"} | to access and download OCS data files


___

The first time we use a function, we will use the `::` to indicate which package we are using (e.g. `stringr::str_detect()`). 
Unless we have overlapping function names, this is not necessary, but we will include it here to be informative about where the functions we will use come from.

## **Context**
***

Previous work has shown higher BMI (>30) is associated with [higher rates of all-causes of mortality](https://www.nhlbi.nih.gov/sites/default/files/media/docs/obesity-evidence-review.pdf){target="_blank"}, as well as increased rates of type 2 diabetes, cancer, heart disease, and stroke. 

However, there appears to be caveats associated with the use of a threshold like this to evaluate public health risk in various populations.

According to this [report](https://www.who.int/nutrition/publications/bmi_asia_strategies.pdf){target="_blank"}:

>...the associations between BMI, percentage of body fat, and body fat distribution differ across populations. 

> In particular, in some Asian populations a specific BMI reflects a higher percentage of body fat.

> Some Pacific populations also have a lower percentage of body fat at a given BMI.

> The associations of BMI and comorbidities are probably not stable within populations over time.

##### [[source]](https://www.who.int/nutrition/publications/bmi_asia_strategies.pdf)

## **Limitations**
***

The measurement of BMI has some [limitations](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4890841/pdf/nt-50-117.pdf){target="_blank"} that are well recognized, as it does not account for the composition of body mass, the location of body fat, or the contribution of body frame size. 
However, [BMI has been a useful health indicator](https://journals.lww.com/acsm-healthfitness/Fulltext/2016/07000/THE_BENEFITS_OF_BODY_MASS_INDEX_AND_WAIST.8.aspx#pdf-link){target="_blank"} for risk for many diseases and conditions particularly when combined with other risk factor information.

While [gender](https://www.genderspectrum.org/quick-links/understanding-gender/){target="_blank"} and [sex](https://www.who.int/genomics/gender/en/index1.html){target="_blank"} are not actually binary, the data presented that is used in this analysis only contain data for groups of individuals described as men or women. 


## **What are the data?**
***
We will be using data located within a table of the [supplementary material](https://static-content.springer.com/esm/art%3A10.1038%2Fs41586-019-1171-x/MediaObjects/41586_2019_1171_MOESM1_ESM.pdf){target="_blank"} for the NCD-RisC paper referenced above. 

This is a pdf that can be found freely available online. Here is a screenshot of the first few rows of this table: 


```{r, echo = FALSE}
knitr::include_graphics(here::here("img","first_page.png"))
```

You can see that the data contain average BMI estimates for both men and women in countries at the national level, as well as the average BMI estimates for the rural and urban areas of these countries in both 1985 and 2017.

The data within the parentheses are the 95% [credible intervals](https://en.wikipedia.org/wiki/Credible_interval){target="_blank"} (CIs) for the average BMI estimates. The authors provide these CIs as a guide to understand how likely the estimate is for the true population mean BMI. A wider range suggests that the estimate is less accurate, as there are more possible values for the true mean with credible evidence. Notice that the mean BMI values are "age-standardised," which means that the mean values were adjusted for the different age distributions of the various countries so that the countries can be more fairly compared. For example, if one country has a population that is considerably younger, the mean BMI might be quite low (as younger individuals tend to have lower BMI values due to faster metabolic rates). One might falsely conclude that people in that country generally have lower BMI values than people in most other countries, however the lower overall BMI mean might simply be due to the fact that the country has a younger population than those other countries.


## **Data Import**
***
First let's download the data from the PDF file. 
We use `file.exists()` from base R to check if the file we want to download already exists. 
If not, then we download it using the function `utils::download.file()`. We also use the `here()` function of the `here` package to look specifically in the `raw` subdirectory of the `data` directory of the directory where our open RStudio project file is located. This helps us, as we do not need to write the full path to the document. This makes our code more reproducible for someone else, who may have a different overall file structure on their computer but downloads our case study files from GitHub and therefore has the same file structure specifically for the files included in our GitHub repository. Thus our commands about file locations still work on their local machine.

```{r}
if(!file.exists(here("data", "raw", "paper_supplement.pdf"))){
  url <- "https://static-content.springer.com/esm/art%3A10.1038%2Fs41586-019-1171-x/MediaObjects/41586_2019_1171_MOESM1_ESM.pdf"
  utils::download.file(url, here("data", "raw", "paper_supplement.pdf"))
}
```


We will save our data to a directory within our working directory called data. We will create subdirectories within this directory to organize our data. We can use the `here` function from the `here` package to make this process easier. The `here()` function allows us to specify the path or location of the document that we want to import, starting from the directory where an `.Rproj` file is located after creating a project in RStudio. In this case, we will import our files within subdirectories of  a directory called `raw` of the `data` directory. (*Note the next chunk of code will only work for you if you pull the repository from GitHub and set up your file structure in the same way.*) If you had trouble downloading the data from the orginal sources you can download them from our [GitHub repository for this case study](https://github.com/opencasestudies/ocs-bp-rural-and-urban-obesity/tree/master/data) or use the `OCSdata` package:

```{r, eval=FALSE}
# install.packages("OCSdata")
library(OCSdata)
raw_data("ocs-bp-rural-and-urban-obesity", outpath = getwd())
# This will save the pdf file in a "OCSdata/data/raw/" subfolder 
# in your current working directory
```

***

<details> <summary> Click here to see more about creating new projects in RStudio. </summary>

You can create a project by going to the File menu of RStudio like so:


```{r, echo = FALSE, out.width="60%"}
knitr::include_graphics(here::here("img", "New_project.png"))
```

You can also do so by clicking the project button:

```{r, echo = FALSE, out.width="60%"}
knitr::include_graphics(here::here("img", "project_button.png"))
```

See [here](https://support.rstudio.com/hc/en-us/articles/200526207-Using-Projects) to learn more about using RStudio projects and [here](https://github.com/jennybc/here_here) to learn more about the `here` package.

</details>
***
</details>

***

Now that we have downloaded the PDF file, we will read it in to R using the `pdftools` package:

```{r}
pdf_obesity <- pdftools::pdf_text(here("data", "raw", "paper_supplement.pdf"))
```


Let's take a look at the data -- the `summary()` function from base R helps us to look at the structure of R objects.

```{r}
summary(pdf_obesity)
```

We can see that we have 63 elements that are character strings. 
You may also notice that the original PDF has 63 pages. Let's take a look at some of these elements and compare them to the pages in the original document.

```{r}
pdf_obesity[1] # this looks like the first page
pdf_obesity[2] # this looks like the second page
```

```{r, eval = TRUE}
str(pdf_obesity[52], nchar.max = 1600)
```
Yes, this looks like data in the table we want. 
We are using the `str()` function to print just a portion of the data. The `str()` function compactly displays the structure of arbitrary R objects. 
The `nchar.max` argument lets us define the maximal number of characters to show. 
Note `print(pdf_obesity[52])` would show the entire contents of this page of the table. 

Here is the PDF version again:
```{r, echo = FALSE}
knitr::include_graphics(here::here("img", "first_page.png"))
```

We can see that the output looks pretty similar to the pages of the pdf, but the spacing is a bit awkward. 
Note that the way the data is displayed is partially influenced by the width setting of the RStudio window.  

Now, we are interested in a supplementary Table 3, which has multiple pages and includes the same header on each page. 
We can use that to determine what elements of our `pdf_obesity` character strings include our table. 
We will use the `str_detect()` function from the `stringr` package to search for the elements that contain text that is consistently in the header. 
The output of of this function will show which elements of the object (in this case pages of the pdf) include this pattern indicated as a `TRUE` or `FALSE`.

```{r}
stringr::str_detect(pattern = "Age-standardised mean BMI", 
                    string = pdf_obesity)
```
We note, the "Age-standardised mean BMI" is part of the header in the table on every page.
The code above shows the pages that contain our table of interest (as `TRUE` values).

Now we will extract just the data for the table now and call it `rural_urban`. 
To do this we will use the `str_subset()` function of the `stringr` package.

```{r}
rural_urban <- stringr::str_subset(pattern = "Age-standardised mean BMI", 
                                   string = pdf_obesity)

summary(rural_urban) 
```

Here, we see there are 10 pages worth of elements in our `rural_urban` object. 

Let's check the first and last page:

```{r, eval=TRUE}
str(rural_urban[1], nchar.max = 1600)
```
Using `[1]` shows just the first page of the 10 pages that contain the table. 
Also, notice the new line `"\n"` values that indicate new lines. 

```{r, echo = FALSE, eval = FALSE}
knitr::include_graphics(here::here("img","first_page_in_R.png"))
```


```{r, echo = FALSE}
knitr::include_graphics(here::here("img","first_page.png"))
```

This looks the same as the beginning... how about the end?

```{r}
str(rural_urban[10], nchar.max = 1800)
```

```{r, echo = FALSE }
knitr::include_graphics(here::here("img","last_page.png"))
```


Great! Our `rural_urban` object looks like it contains the entire Supplementary 3 table, as both the beginning and the end include the data we expected.

Now we will save this imported data for later use like so in a subdirectory of the `data` directory called `imported`:

```{r}
save(rural_urban, file = here::here("data", "imported", "imported_data.rda"))
```

## **Data Wrangling**
***
If you have been following along but stopped, we could load our imported data from the "data" directory like so:

```{r}
load(file = here::here("data", "imported", "imported_data.rda"))
```

***
<details> <summary> If you skipped the data import section click here. </summary>

First you need to install and load the `OCSdata` package:

```{r, eval=FALSE}
install.packages("OCSdata")
library(OCSdata)
```

Then, you may load the imported data using the following code:

```{r, eval=FALSE}
imported_data("ocs-bp-rural-and-urban-obesity", outpath = getwd())
load(here::here("OCSdata", "data", "imported", "imported_data.rda"))
```

If the package does not work for you, alternatively, an RDA file (stands for R data) of the data can be found [here](https://github.com//opencasestudies/ocs-bp-rural-and-urban-obesity/master/data/imported) or slightly more directly [here](https://raw.githubusercontent.com/opencasestudies/ocs-bp-rural-and-urban-obesity/master/data/imported/impoted_data.rda). Download this file and then place it in your current working directory within a subdirectory called `imported` within a subdirectory called `data` to copy and paste our code. We used an RStudio project and the [`here` package](https://github.com/jennybc/here_here) to navigate to the file more easily. 

```{r}
load(here::here("data", "imported", "imported_data.rda"))
```

***
<details> <summary> Click here to see more about creating new projects in RStudio. </summary>

You can create a project by going to the File menu of RStudio like so:


```{r, echo = FALSE, out.width="60%"}
knitr::include_graphics(here::here("img", "New_project.png"))
```

You can also do so by clicking the project button:

```{r, echo = FALSE, out.width="60%"}
knitr::include_graphics(here::here("img", "project_button.png"))
```

See [here](https://support.rstudio.com/hc/en-us/articles/200526207-Using-Projects) to learn more about using RStudio projects and [here](https://github.com/jennybc/here_here) to learn more about the `here` package.

</details>
***
</details>
***

At this point we have large strings now for each page of the table, but this is not very convenient to work with. 
Now, we will wrangle the data into a more usable form. 
Ideally, we would like our data to be in some sort of tabular form.



### **Separate the data into lines**
***
First, it would be useful to separate the data that is currently in 10 giant string chunks into individual lines or rows of the table. To do this we can use the `read_lines()` function of the `readr` package which will separate the lines based on the `"\n"` new line characters.
```{r}
summary(rural_urban)
rural_urban <- readr::read_lines(rural_urban)
summary(rural_urban) 
# now we have 461 lines
```
We can see that each string in the `rural_urban` object is now a single line from the table.

For example, line 6 shows the data for women in Afghanistan.
```{r}
rural_urban[6] 
```


### **Removing excess white-space**
***

We also have a lot of white-space. Let's get rid of the excess white spaces using `str_squish()` also from the `stringr` package.

```{r}
rural_urban <- stringr::str_squish(rural_urban) 
head(rural_urban)
```

Now it is much easier to see the data.

If we look at the part of our data that contains the end of the first page of the table and the start of the second we can see that the header information is repeated, as well as a line with the page number and an empty line, and a line that says "2 2". The "2 2" comes from importing the PDF file, these are the superscript "2"s in the $m^2$ strings in the table header - no need to worry about this. You may also notice some values that are `na*` and we will deal with this later.

#### {.scrollable }
```{r}
rural_urban[44:56] # scroll through!
```
####

Although the header was necessary on all of the pages of the pdf version of the table, we only need that information once in our data.

### **Removing unnecessary repeated header information**
***
So, let's remove all the header information and the page number lines from the `rural_urban` object, then we will make a single line header for the beginning.
One way to do this is to find all lines that include either "Women" or "Men" and only keep this data.

First, let's see how many lines include "Women" or "Men". 
We can do this by first identifying the lines that contain these patterns and then check the length of the resulting vector for the line numbers that contain these patterns. 
Note, the `"|"` indicates that we are looking for either "Women" or "Men".

#### {.scrollable }
```{r}
#scroll through the output!
stringr::str_which(string = rural_urban, 
                   pattern = "Women|Men")
```
####

```{r}
length(stringr::str_which(string = rural_urban, 
                          pattern = "Women|Men"))

```

Note that `str_which()` is case sensitive, so it would not work to use "women" as the pattern, and using "men" would return the lines that contain "Wo**men**" or "Ye**men**" etc.  Try this out with `pattern = "women"` and `pattern = "men"` to see the result!

OK, so this looks correct. 
This includes most lines but there are gaps where the header is located. It looks like out of the original 461 lines, there are 400 lines in our table that aren't headers.

```{r}
rural_urban <- str_subset(string = rural_urban, 
                          pattern = "Women|Men")
```

We can check our data now using either `head()` from base R or `glimpse()` from the `tibble`  and `dplyr` package. (Yup, the same function is available from both packages!)  The `head()` function shows us the first rows or lines of the data, while the `glimpse()` function provides us information about the total size of the object and shows us the first line or row.

```{r}
dplyr::glimpse(rural_urban)
```

Now our `rural_urban` object 400 lines instead of 461, like it did before when it contained the redundant header information.

```{r}
head(rural_urban)
```

Great! So now our data looks much better but we need to add back our header and we would like this to only be a single line to make it easy to transform our data into a table or table-like object.

As a reminder, here is what our header used to look like:

```{r, echo = FALSE,  out.width = '100%'}
knitr::include_graphics(here::here("img","last_page.png"))
```

We will add a new header after first dealing with some additional spacing and missing data issues.

### **Dealing with spacing**
***
First let's try splitting our header-less data into columns based on spaces using the `str_split()` function, where we specify that we are splitting the data based on the pattern of a space (the space is included in quotes). Referring to [str_split() documentation](https://www.rdocumentation.org/packages/stringr/versions/1.4.0/topics/str_split), the `simplify` argument is set to `FALSE` by default, where a list of character vectors is returned. If it is set to `TRUE` as we will be doing here, a character matrix is returned.

Here, we will take a look at just the first 10 lines using the `[1:10,]` syntax: 


#### {.scrollable }
```{r}
#scroll through the output!
str_split(string = rural_urban, 
          pattern = " ", 
          simplify = TRUE)[1:10,] 
```
####

This almost worked, but unfortunately country names that have spaces will be a problem. 
For example, we can see that "American Samoa" has been divided into two columns and all subsequent columns are shifted.

Let's try this another way.  

As we can see by looking at our `rural_urban` data, the country information is only present when the sex is Women. Let's try to extract the country information from the rows that contain data about Women.

```{r}
head(rural_urban)
```
To do this, we start by noticing that the sex always starts with either a capital "W" if the sex is female. Note, we need to use a space before the "W" otherwise we will split some of the country names if the names starts with "W". 

Here, we will also introduce the concept of piping, which uses the `%>%` from the `magrittr` package, which allows us to pipe the output from one function to the input for another function.
This is really useful when we have multiple steps, which we will show soon.

First, we will select just the data for Women and Men separately:

```{r}

Women <- str_subset(string = rural_urban, 
                    pattern = "Women") 
Men <- str_subset(string = rural_urban, 
                    pattern = "Men")


head(Women)
head(Men)
```


And then split the `Women` data based on the pattern `pattern = " Women"` (Note that the space before the word "Women" is included - so that it is not within the country data).
Here we take our `Women` dataset and pipe it into the first argument of the `str_split()` function, where then we can split on the pattern `" Women"`.  We then pipe through the `unlist()` function so that the results are a vector instead of a list.  We assign this new result to the `country_split` object.

```{r}
country_split <- 
  Women %>%
  stringr::str_split(pattern = " Women") %>%
  unlist() 

head(country_split)
```

However, we can also do the same thing in a single step, where we replace the `Women` dataset with the code we used to create that dataset.  So first, we select the "Women" and then we split the data by " Women" to create a new dataset all within the same chunk of code.
```{r}
country_split <- 
  str_subset(string = rural_urban,
             pattern = "Women") %>%
  stringr::str_split(pattern= " Women") %>%
  unlist()

head(country_split)
```

Now, we can see that `Country` is always the odd rows and the  BMI data for Women is the even rows. 
We can select the odd rows using the <b> Modulus (Remainder from division) operator</b>. 
This operator looks like this `%%`. 
All odd values such as 3, 7, or 15 when divided by 2 would leave a remainder of 1. 
We will select just these rows using the `filter()` function of `dplyr`.

```{r}
country <- 
  tibble(country_split) %>% 
  dplyr::filter(row_number() %% 2 == 1) 
```

We can take a look to make sure that all the country names look as expected. You can scroll through the country names if you hover over the data.

```{r, eval= FALSE}
country
```

#### {.scrollable }
```{r, echo = FALSE}
as.data.frame(country)
```
####

Looks good!  We have 200 countries represented in our dataset.

Great! Now we have a list of the countries that can be used for both the male and female data. 


Remember the even rows are the rows with the BMI data for women.
These row values have a remainder of 0 when divided by 2.
Let's select that data too.

```{r}
Women_BMI <- 
  tibble(country_split) %>% 
  dplyr::filter(row_number() %% 2 == 0)

head(Women_BMI)
```
This looks similar to our `Men` data (although slightly different):

```{r}
head(Men)
```

Notice how the data in each row is in quotes. This is because these are character strings. We can check this by using the base `class()` function.

```{r}
class(Men)
```

Whereas the `Women_BMI` data is a tibble.
```{r}
class(Women_BMI)
```

Let's make the `Men` data into a tibble as well:
```{r}
Men <- tibble(Men)
```

It's always a good idea to check that your data objects are the size you expect when wrangling. 
We can do so with the `dim()` function, which shows us the dimensions of data objects.

```{r}
dim(Women_BMI)
dim(Men)
head(Women_BMI)
head(Men)
```
Great! There are 200 rows of character strings with the BMI data like we expected to match the 200 countries in our data.

We do however need to add the word `" Women"` back to the `Women_BMI` data.
To do this we will use the `glue()` function of the `glue` package. 

Here is a simple example of how this function glues two things together.

```{r}
where <-"the moon"

glue::glue("the cow jumped over {where}")

```

Now let's add  the word "Women" to each row of the `Women_BMI` data using `glue`.

First let's pull the `country_split` data out as a vector to make this simpler when we `glue`. To do this we will use the `pull()` function of the `dplyr` package. This function will allow us to extract or pull out the single variable called `country_split` which contains the BMI values for women from within the tibble called `Women_BMI` which we just created by grabbing only the even rows of our original `country_split` tibble.

Then, let's use the `first()` function of the `dplyr` package to get a reminder of what just the first row of the `country_split` column of the `Women_BMI` data looks like.


```{r}
women_data <- Women_BMI %>%
  pull(country_split)

dplyr::first(women_data)
```

Now we will paste the word `"Women"` in the beginning of each value in the vector using the `glue()` function.

Note that there will be a space in between.
This is because the space is already in front of the first BMI value.

```{r}
Women <- glue("Women{women_data}")
head(Women)
```

You could also do all of this in one step but its less obvious what is happening.

Here we can pull the `country_split` data out of the `Women_BMI` object and then paste `"Women"` in front of each string in one command.
```{r, eval = FALSE}
Women <- glue("Women{pull(Women_BMI,country_split)}")
```



If we try splitting our data by space again, will it have the expected number of columns? What about the rows that contain `na*` values?

Let's just take the `Men` data that contain `na*` values. 
This column is called `rural_urban`. 

```{r}
unknown <- 
  Men%>% 
  filter(str_detect(pattern ="na\\*", 
                    string = Men)) 

head(unknown)
```

Now we can try splitting by a space
```{r}
tibble::as_tibble(str_split(pull(unknown, Men), " ",
                            simplify = TRUE))

```

So close! Notice that the `"na*"` values have shifted the subsequent values within the columns because typically there is a space between the BMI and the credible intervals.
Here we can see this data in our original pdf:

```{r, echo = FALSE, out.width="100%"}
knitr::include_graphics(here::here("img","missing_pdf.png"))
```


### **Dealing with NA values**
***

We need to replace our `na*` values with something that includes a space so that when we separate our data by space, we will have two values instead of one when we have an `na*`. 
Therefore, `na* na*` should work.

So, in order to use the functions within the `stringr` package, our `Men` data needs to be of `character` class, so let's check that now.

```{r}
head(Men)
```
The `<chr>` indicates that our `Men` data is already `character`. But if this is not the case, we could use the ` as.character()` function that allows us to change factor data within to the character class.

We are also changing the structure so that the data is not within a column called `Men`. 
```{r}
Men <- pull(Men, Men)
head(Men)
```

```{r}
Men <- 
  stringr::str_replace_all(string = Men, 
                           pattern = "na\\*", 
                           replacement = "na\\* na\\*") 

# The \\ are necessary because the * is a special character 
# The * would typically indicate any possible value, 
# but here we actually want a "*" instead
# Thus the double backslash does that for us
# Here we are replacing all occurrences of the na* values 
#(thus str_replace_all instead of str_replace) with na* na*. 
```

Let's check that it worked...
```{r}
Men[20:30]
```

Great!

Now for the Women data object
```{r}
class(Women)
# the female data is already character type
Women <- 
  stringr::str_replace_all(string = Women, 
                           pattern = "na\\*", 
                           replacement = "na\\* na\\*") 
Women[20:30]
```

Great, now we can split our data by spaces.

### **Splitting the data**
***
We will split our data finally by spaces using the `str_split()` function of the `stringr` package. We can check the rows with NA values by looking at rows 20-30.
```{r}
Men <- tibble::as_tibble(str_split(Men, " ", 
                                   simplify = TRUE))
Men[20:30,]
```
Great! It looks good!

Now for the `Women` data.

```{r}
Women <- as_tibble(str_split(Women, " ", 
                             simplify = TRUE))
head(Women)
Women[20:30,] 
```

We can see that our `na` values look correct. 

Looks good!
 
 
### **Adding new header**
***
We can see from our pdf and our object called `header` what the header was like in the original pdf document. We need to add column names to our data based on the original data.

We'll wait to add `Country` to the header until we add the country information to our data.

```{r}
new_header <- c("Sex","National_BMI_1985", 
                "National_BMI_1985_CI", "Rural_BMI_1985", 
                "Rural_BMI_1985_CI", "Urban_BMI_1985",
                "Urban_BMI_1985_CI", "National_BMI_2017",
                "National_BMI_2017_CI","Rural_BMI_2017",
                "Rural_BMI_2017_CI", "Urban_BMI_2017", 
                "Urban_BMI_2017_CI")
```

Let's change the names of our columns of our tibbles to this new header for our `Men` and `Women` data

```{r}
names(Women) <- new_header
names(Men) <- new_header

head(Women)
head(Men)
```

Now we will add our country data to both our `Men` and `Women` tibbles using the `bind_cols()` function. This will add the `country` as a new column to the `Women` data object on the left. 

```{r}
Women <- dplyr::bind_cols(country, Women)
head(Women)
```


Similarly, for the data about Men, this will add the `country` as a new column to the `Men` data object on the left. 
```{r}
Men <- dplyr::bind_cols(country, Men) 
head(Men)
```


### **Changing variable names**
***
Here, we will change the variable name in the `country` dataset to `country` (currently it is called `country_split`). 
We will also introduce the concept of the assignment pipe. 
In this case our pipe operator looks like this `%<>%`. 
Using this additional pipe requires another package, called `magrittr`. 
The other simpler pipe options from this package are loaded with `tidyverse` (if you used `library(tidyverse)` which loads most tidyverse packages), but not this version. 

The `>` portion of the pipe still behaves like a normal pipe, while the `<` portion of the pipe makes an assignment to whatever the `<`is pointing to, just like when we use the typical assignment operator `<-`. 

```{r}
# library(magrittr) 
# We can't use the `%<>%` unless we load the magrittr package
# We have already done this but we include this for illustrative purposes.
# Here we will just use the traditional assignment strategy

Women <- dplyr::rename(Women, Country = country_split) 
```
Here, we have renamed the `country_split` variable to be called `Country`. 

Here, we reassign `Men` using the pipe strategy.
```{r}
Men %<>% rename(Country = country_split) 
```

This is equivalent to
```{r, eval=FALSE}
Men <- rename(Men, Country = country_split)
```
We have renamed the `country_split` variable to to be called `Country`. 
We have also reassigned `Men` to the same data object. which has the `Country` variable renamed. 

Let's take a look at our data for `Men` and `Women`:

```{r}
head(Women)
head(Men)
```
### **Joining the data**
***
Now that both of these tibbles have the same structure, we can combine our `Men` and `Women` data using the `bind_rows()` function of `dplyr`. This will combine all of the rows based on their similar column structure.
All of the male data  will be first and all the female data will be second.

```{r}
BMI <- dplyr::bind_rows(Men, Women)
```

Let's check the size of our BMI data... it should have 400 rows (obs). We can see the end of the `BMI` data using the  `tail()` function:
```{r}
str(BMI)
tail(BMI)
```

### **Sorting the data**
***

Now, let's sort the data alphabetically by `Country` using the `arrange()` function from the `dplyr` package:

```{r}
BMI <- dplyr::arrange(BMI, Country)
head(BMI)
```

Our data is looking great! 

Now, we might want to make sure that our observations for each variable look the way we want. 
In other words, if we want to make plots about National BMI in 1985 then we would need our values to be numeric. 
Looking at our BMI data using `str()`, we can see the type of data for each of our variables listed just after variable name and the `":"` colon. 

```{r, warning=FALSE, eval}
str(BMI)
```

Looks like none of our BMI data is actually `numeric` (or `num`), but of the class `character` (or `chr`). 
Let's change that now.

To start, we could change these values to be numeric with 6 lines of code like this, where the double brackets`[[]]` and quotes are used to select the specific variable within the `BMI` tibble. In this case we need to use the brackets instead of `pull()` because we can't start a command with `pull()` to reassign values in a variable.

```{r, eval = TRUE}
#pull(BMI, National_BMI_1985)%<>% as.numeric #This will not work
BMI[["National_BMI_1985"]] %<>% as.numeric
BMI[["Rural_BMI_1985"]] %<>% as.numeric
BMI[["Urban_BMI_1985"]] %<>% as.numeric

BMI[["National_BMI_2017"]] %<>% as.numeric
BMI[["Rural_BMI_2017"]] %<>% as.numeric
BMI[["Urban_BMI_2017"]] %<>% as.numeric
```

And if we did this we would see these variables show as `"num"` which stands for a numeric type of variable.
```{r}
str(BMI)
```

Or we can use a more automated way with the `map()` function of the `purrr` package:
```{r}
BMI_numeric <- 
  purrr::map((BMI %>% 
                select(-matches("CI|Sex|Country"))),
             as.numeric) %>% 
  as_tibble()
```

The `map()` function of the `purrr` package allows us to apply the `as.numeric()` function to all the selected columns of `BMI` dataset.
Above, we have selected those that don't contain `CI`, `Sex`, or `Country` in the column name. 
If you are familiar with `apply()`, this function is very similar.

Now we need to add the `Country` and `Sex` variables back to our new numeric dataset. To do this we will use the `mutate()` function from the `dplyr` package. This function allows us to create and modify variables.

```{r}
BMI_numeric %<>% 
  mutate(Country = pull(BMI, Country), 
         Sex = pull(BMI, Sex))
head(BMI_numeric)
```

In this case our numeric data is of a class `dbl`. 
This means double-precision floating point numbers (with decimals).
The alternative numeric option is the integer class. 

### **Transforming to tidy data**
***

It is generally useful to get the data in what is called <b>long</b> format for other analyses, and particularly for plotting. 

For a more detailed description about this, please see the Data Visualization section of this [case study](https://www.opencasestudies.org/ocs-bp-opioid-rural-urban/#Rural_and_Urban_areas/#:~:text=%20the%20data%20was%20presented%20in%20a%20format%20that%20is%20called%20long%20format.){target="_blank"}.

Basically, the long format has more rows and fewer columns.
Thus, we can combine some columns together.
In our case, we can put all the different BMI data together in one long column.
We will create a new column that tells us what each row of the new BMI column represents.
In other words, it will tell us what the original column was.

To do this we will use the `pivot_longer()` function of the `tidyr` package. The `data` argument specifies what data we want to work with, the `cols` argument specifies the names of the columns or variables that we want to combine and change to one long format variable, and the `names_to` argument specifies the name for the new variable that we are creating. The `:` is used to indicate that we want to combine the columns or variables between `National_BMI_1985` and `Urban_BMI_2017`.

```{r}
BMI_long <-
  tidyr::pivot_longer(data = BMI_numeric, 
                      cols = National_BMI_1985:Urban_BMI_2017, 
                      names_to = "class_BMI")

head(BMI_long)
```

We really want to promote the idea of working with [tidy data](https://r4ds.had.co.nz/tidy-data.html){target="_blank"}. Tidy data is data that is structured in a way that makes it especially usable. One of the rules of tidy data, is that each column or variable should have only one type of data. Currently we have two types of data in the `class_BMI` variable: the type of BMI measurement ( Rural, Urban, or National), and the year (1985 and 2017).  We can separate these two data types of the `class_BMI` column using the `separate()` function of the `tidyr` package to replace the `class_BMI` column with two new columns.

```{r}
BMI_long %<>% 
  tidyr::separate(class_BMI, 
                  c("Region", NA, "Year")) 
```

Here, we separate the `class_BMI` column of the `BMI_long` tibble. 
This is based on the underscores (“_”) and creates two new columns. 
The first column will be called `Region`. 
It will contain the 1st part of the `class_BMI` data before the 1st underscore. 
The 2nd column will be called `Year`.
It will contain the 3rd part of the data based on underscores.
The `NA` in our code is for removing the middle part of the column (“BMI”). Thus, this part will not be used to create any new columns.

```{r}
head(BMI_long)
```

Let's see how the dimensions of the `BMI` data have changed in long format:
```{r}
dim(BMI_numeric)
dim(BMI_long)
```
Let's rename the `value` column to something more informative. We can use the `rename()` function of the `dplyr` package to do this. Notice that the value we want to replace is listed second after the `=` sign.

```{r}

BMI_long %<>% 
  rename(BMI = value)
head(BMI_long)

```
#### {.question_block}
<u>Question opportunity:</u> Why exactly are there 2,400 rows now?

####

Great! Our data is very usable now in this format! Now let's save our wrangled data for future use. We will save it as an r compatible file (a .rda file), as well as csv files, as these are useful to give to collaborators. We will use the `write_csv()` function of the `readr` package to do this.This requires that each tibble be saved one at a time.

```{r, eval = FALSE}
save(BMI_numeric, 
     BMI_long, 
     file = here::here("data",  "wrangled", "wrangled_data.rda"))
readr::write_csv(BMI_numeric, 
                 path = here::here("data", "wrangled", "BMI_numeric.csv"))
readr::write_csv(BMI_long, 
                 path = here::here("data", "wrangled", "BMI_long.csv"))

```


## **Data Exploration**
***

If you have been following along but stopped, we could load our wrangled data like so:

```{r}
load(here::here("data", "wrangled", "wrangled_data.rda"))
```

***
<details> <summary> If you skipped the data wrangling section click here. </summary>

First you need to install and load the `OCSdata` package:

```{r, eval=FALSE}
install.packages("OCSdata")
library(OCSdata)
```

Then, you may load the wrangled data using the following code:

```{r, eval=FALSE}
wrangled_rda("ocs-bp-rural-and-urban-obesity", outpath = getwd())
load(here::here("OCSdata", "data", "wrangled", "wrangled_data.rda"))
```

If the package does not work for you, alternatively, an RDA file (stands for R data) of the data can be found [here](https://github.com//opencasestudies/ocs-bp-rural-and-urban-obesity/tree/master/data/wrangled) or slightly more directly [here](https://raw.githubusercontent.com/opencasestudies/ocs-bp-rural-and-urban-obesity/master/data/wrangled/wrangled_data.rda). Download this file and then place it in your current working directory within a subdirectory called "wrangled" within a subdirectory called "data" to copy and paste our code. We used an RStudio project and the [`here` package](https://github.com/jennybc/here_here) to navigate to the file more easily.

```{r}
load(here::here("data", "wrangled", "wrangled_data.rda"))
```


***
<details> <summary> Click here to see more about creating new projects in RStudio. </summary>

You can create a project by going to the File menu of RStudio like so:


```{r, echo = FALSE, out.width="60%"}
knitr::include_graphics(here::here("img", "New_project.png"))
```

You can also do so by clicking the project button:

```{r, echo = FALSE, out.width="60%"}
knitr::include_graphics(here::here("img", "project_button.png"))
```

See [here](https://support.rstudio.com/hc/en-us/articles/200526207-Using-Projects) to learn more about using RStudio projects and [here](https://github.com/jennybc/here_here) to learn more about the `here` package.

</details>
***
</details>
***

Now that our data is clean and in a format that we can work with, we can start to take a look at the data and how different groups might compare.

Going back to one of our original questions, 

> Is there a difference between rural and urban BMI estimates around the world?

We want to use a statistical test to identify if there are differences in the average BMI estimates between e.g. rural and urban groups. 
First, let's explore our data a bit. 

### **General summary**
***


```{r}
summary(BMI_numeric)
```
We can see that the mean BMI estimates are larger in 2017.
However, is this a meaningful increase?
This is difficult to determine without a statistical test.

Lets look at the data stratified by `Sex`:

```{r}
BMI_numeric %>% 
  filter(Sex == "Women") %>% 
  summary()
```

```{r}
BMI_numeric %>% 
  filter(Sex == "Men") %>% 
  summary()
```

It looks like mean BMIs have increased in all regions for both men and women, but let's look deeper with some statistical tests.

### **Distribution**
***

In order to apply a statistical test to compare the means, one of the first things to do is to explore the frequency of the different observed values. 
One way to summarize the frequency of different observed values is the <b>frequency distribution</b>, which can be shown in a table or a plot. 
See [here](http://onlinestatbook.com/2/introduction/distributions.html){target="_blank"} for more information about distributions. 

There are other types of distributions too, such as <b>probability distributions</b>, which describe the probability the occurrence of different possible outcomes. 
Probability distributions can be defined for discrete data (e.g. flipping a coin as heads or tails) and continuous data (e.g. BMI estimates). 
If we are talking about discrete data, the relative frequency of the different categories (heads or tails) can be used to define a <b>discrete probability distribution</b>. 
We can simply assign a probability to each category, e.g. $P(heads) = 0.5$ and $P(tails) = 0.5$. 

<b>Continuous probability distributions</b> work in a similar way, except instead of assigning a probability to each category or type of observation (e.g. heads or tails), we assign a probability to an interval of continuous values (e.g. BMI between 22 and 26). 
This is because it's not useful to construct a distribution that assigns a probability to each possible outcome (e.g. a BMI of 22.11, 22.12, etc) with such high precision. In this case, we would have to assign a very small probability to every possible BMI outcome. 

The <b>normal (or Gaussian) distribution</b> is a type of <b>theoretical continuous probability distribution</b> that is frequently used to approximate the distribution of many naturally occurring measurements, such as that of BMI. 
We say that a random quantity $X$ (e.g. BMI) is normally distributed with mean $m$ and standard deviation $s$ if the proportion of values in a given interval $(a,b)$ (e.g. BMI between 22 and 26) can be computed using this mathematical formula: 

$$P(a < x < b) = \int_a^b \frac{1}{\sqrt{2\pi} s} e^{-\frac{1}{2}\big(\frac{x-m}{s}\big)^2}  $$

Translating this to our case study, if our BMI data follow a normal distribution with mean of $m = 24$ and standard deviation of $s = 1$, then the distribution of different BMI outcomes should be centered around that mean value of 24 and would look something like this:

```{r}
norm_BMI_ex_data <- 
  tibble(norm_data = rnorm(n = 200, mean = 24, sd = 1))
hist(pull(norm_BMI_ex_data, norm_data))
```

Here we are simulating 200 observations from a normal distribution with a mean of $m = 24$ and standard deviation of $s = 1$ using the `rnorm()` function. 

Alternatively to using base R, we can plot the frequency of simulated observations using the `geom_hist()` function of the `ggplot2` R package. 
The [ggplot2](https://ggplot2.tidyverse.org/){target="_blank"} package creates plots by using layers.
Notice in the following code how there is a plus sign between the `ggplot()` function and the `geom_histogram()` function. 
With `ggplot2` we select what data we would like to plot using the first function (`ggplot()`) and then we add on additional layers of complexity (these layers can even involve different data). 

***
<details> <summary> Click here for additional information on using `ggplot2`.</summary>

The reasons [`ggplot2`](http://ggplot2.tidyverse.org) is generally intuitive for beginners is the use of [grammar of graphics](http://vita.had.co.nz/papers/layered-grammar.html) or the `gg` in `ggplot2`. 
The idea is that you can construct many sentences by learning just a few nouns, adjectives, and verbs. 
There are specific "words" that we will introduce, and once you are comfortable with them, you will be able to create (or "write") hundreds of different plots.

**_What is the `ggplot()` function?_**

As explained by Hadley Wickham: 

> the grammar tells us that a statistical graphic is a mapping from data to aesthetic attributes (colour, shape, size) of geometric objects (points, lines, bars). The plot may also contain statistical transformations of the data and is drawn on a specific coordinates system.

**_`ggplot2` Terminology_**

* **ggplot** - the main function where you specify the dataset and variables to plot (this is where we define the `x` and
`y` variable names)
* **geoms** - geometric objects
    * e.g. `geom_point()`, `geom_bar()`, `geom_line()`, `geom_histogram()`
* **aes** - aesthetics
    * shape, transparency, color, fill, line types
* **scales** - define how your data will be plotted
    * continuous, discrete, log, etc

The function `aes()` is an aesthetic mapping function inside the `ggplot()` object. 
We use this function to specify plot attributes (e.g. `x` and `y` variable names) that will not change as we add more layers.  

Anything that goes in the `ggplot()` object becomes a global setting. 
From there, we use the `geom` objects to add more layers to the base `ggplot()` object. 
These will define what we are interested in illustrating using the data.

</details> 
***

```{r}
  ggplot2::ggplot(norm_BMI_ex_data, aes(x = norm_data)) + 
  ggplot2::geom_histogram() 
```

Here, we are using the column called `norm_data` within `norm_BMI_ex_data`.
As we noted, with `ggplot2` we must first specify the data to plot with the `ggplot()` function. 
Then, the `geom_*` function specifies what type of plot to create (e.g. `geom_histogram()` create a histogram). 

We could also make this same plot by piping the `norm_BMI_ex_data` data into the `ggplot()` function:
```{r}
norm_BMI_ex_data %>% 
  ggplot2::ggplot(aes(x = norm_data)) + 
  ggplot2::geom_histogram() 
```

We will see later how we can add many layers to plots with `ggplot2`.


Now instead of using simulated data from a normal distribution, let's make a histogram of our BMI data:
```{r}
BMI_long %>% 
  ggplot(aes(x = BMI)) + 
  geom_histogram() 
```

Interesting! The data looks like it has one mode, but it does not look completely symmetric. In fact it looks like it might be what is called right skewed. See [here](http://onlinestatbook.com/2/glossary/skew.html){target="_blank"} for more information about skewed distributions.

To check we can calculate the mean and median of the BMI data using the base `summarize()` function of the `dplyr` package. If the mean and median are very similar, then we will know that the data isn't skewed. Note that it is necessary to filter out the NA values first.

```{r}
BMI_long %>%
  filter(!is.na(BMI))%>%
  summarize(mean = mean(BMI), median = median(BMI))
```
OK, so they are essentially the same, mean= 24.87 and median = 24.9, so this does not fit the hallmark of right-skewed data where the the mean is greater than the median. 


However so far we have made a histogram of the BMI data in the long format which combines BMI estimates for different subgroups of people, so rural women, rural men, urban women, urban men, etc. are all combined. If we want to make group comparisons, we need to look at the distributions for the specific groups that we are evaluating. So, we need to pick a particular population of individuals and graph their BMI, for example just the national estimates for women and for men in 1985.  

The easiest way to do this is to use some other functions of the `ggplot2` package, particularly the `facet_wrap()` function, which will allow us to look at differences in the distribution of the BMI estimates stratified by `Year`, `Sex`, and `Region`.
We can sequentially divide our plots by deeper levels using multiple variables and the `+` plus sign within `facet_wrap()`.

```{r}
BMI_long %>%
  ggplot(aes(x=BMI)) +
  geom_histogram() +
  facet_wrap(~ Sex + Year + Region) 
```

Some of these plots look pretty similar to a normal distribution, like the `Women 2017 Urban` plot (although it might be right-skewed, again we could check the mean and median to determine that).
However some of the other plots, like the `Men 2017 Urban` plot, look very different. 
These plots show what is called a [bimodal distribution](http://onlinestatbook.com/2/introduction/distributions.html){target="_blank"}, meaning that there appear to be two peaks (or modes), or in other words two different values (and values near each of these) that are the most frequent.

### **Quantile-Quantile Plots**
***

Another useful plot to create when exploring distributions is something called a "quantile-quantile" plot (or <b>Q-Q plot</b> for short). When we talk about a <b>quantile</b>, we are talking about dividing up the distribution of the data into roughly equal portions where roughly the same number of observations fall into each portion. For example, if you divide your data into 100 quantiles, you can think about this as percentiles, but you could also divide your data into 10 quantiles and these would be called deciles. 

<b>Why Q-Q plots?</b> This plot allows us to compare the quantiles of two distributions together: (1) quantiles of a known theoretical distribution (like the normal distribution) compared to (2) quantiles of the distribution of our data.  If the quantiles from these two distributions line up in the plot, then that is a visual piece of evidence that our data follow that theoretical distribution (like the normal distribution).  

<b>How does this work?</b> To do this we will plot the quantiles of our data on the y-axis and the quantiles of the theoretical normal distribution on the x-axis. If the quantiles line up then we can say that our data is fairly normal. See [here](http://onlinestatbook.com/2/advanced_graphs/q-q_plots.html){target="_blank"} for more information about Q-Q plots.

Using the `stat_qq()` function of the `ggplot2` package, we can easily create a Q-Q plot for our data randomly sampled from a normal distribution ("sample") and compare it to the quantiles from a normal distribution ("theoretical").  The default comparison distribution for these functions is the normal distribution, so we don't need to specify it in our code.
```{r}
norm_BMI_ex_data %>%
  ggplot(aes(sample = norm_data)) +
  stat_qq() + 
  stat_qq_line()
```

The `stat_qq_line()` function is used to add a line (computes the slope an intercept) on the plot. If the points lie on this straight line, this is evidence that the data have a normal distribution, not that the data have a particular normal distribution.
Not surprisingly, we see that the points mostly fall on the line, indicating that the quantiles are fairly similar between the observed and theoretical data.

There are some points that are a bit further from the line as we get to the extreme quantiles.

Note that this plot is comparing our `norm_BMI_ex_data` which was randomly sampled from a normal distribution with mean 24 and standard deviation 1 to a *standard normal distribution* with mean 0 and standard deviation 1, which is why the values on the x-axis are -3 to 3 instead of 21 to 27. 

As expected we see that about half the points are below our mean of 24 on the y-axis.  

We can specify a specific mean and standard deviation for our theoretical distribution in our Q-Q plot; here we specify a normal distribution with mean 24 and standard deviation 1 (so that the theoretical distribution looks more like our data):
```{r}
norm_BMI_ex_data %>%
  ggplot(aes(sample = norm_data)) +
  stat_qq(distribution = qnorm, dparams = c(mean = 24, sd = 1)) + 
  stat_qq_line(distribution = qnorm, dparams = c(mean = 24, sd = 1))
```
Notice that the plot looks the same except for the values on the x-axis have changed, meaning our interpretation will be the same -- it looks like the values are consistent with a normal distribution.


Now, let's take a look at our actual BMI data:
```{r}
BMI_long %>%
  ggplot(aes(sample =BMI)) +
  stat_qq() + 
  stat_qq_line() +
  facet_wrap(~ Sex + Year + Region)
```

This is a little hard to see, some plots suggest those subsets of the data are similar to a normal distribution, but others seem to be skewed. Let's take a closer look.

For the sake of simplicity, we are going to focus on the subset of data for the women. 
This is because women were identified in the paper to have fastest increases in BMI over time.
If we want to compare two subsets of women (e.g. 1985 versus 2017 or rural versus urban), then we need to look at the distribution of each of these subsets.

```{r}
BMI_long %>%
  filter(Sex == "Women") %>%
  ggplot(aes(sample = BMI)) +
  stat_qq() + 
  stat_qq_line() +
  facet_wrap(~ Year + Region)
```

We can see that at the extremes of our quantiles, for most of our data, our tails are not very similar to the theoretical distribution. 
We can can also see that our data is right- or positive-skewed, because the extreme values are above the line.

The rural data looks more normal than the urban data. 
You might be asking, why do we care if the data look more normal or not? 
It turns out if you want to ask questions like "Is there a difference in the average BMI estimates from women in 1985 compared to women in 2017?", we will need to use some of type of statistical test. 
The most commonly used statistical tests that can answer that question rely on an assumption about the data you have, namely that your data have a similar shape as a theoretical normal distribution. 
So, if we want to use those powerful statistical tests that rely on this assumption of "normality", then we need to check that assumption before using those tests. 

We can conclude from these plots that the assumption of normality appears to be violated in our data.


## **Data Analysis**
***

If you have been following along but stopped, we could load our wrangled data like so:

```{r}
load(here::here("data", "wrangled", "wrangled_data.rda"))
```

***
<details> <summary> If you skipped the data wrangling section click here. </summary>

First you need to install and load the `OCSdata` package:

```{r, eval=FALSE}
install.packages("OCSdata")
library(OCSdata)
```

Then, you may load the wrangled data using the following code:

```{r, eval=FALSE}
wrangled_rda("ocs-bp-rural-and-urban-obesity", outpath = getwd())
load(here::here("OCSdata", "data", "wrangled", "wrangled_data.rda"))
```

If the package does not work for you, alternatively, an RDA file (stands for R data) of the data can be found [here](https://github.com//opencasestudies/ocs-bp-rural-and-urban-obesity/tree/master/data/wrangled) or slightly more directly [here](https://raw.githubusercontent.com/opencasestudies/ocs-bp-rural-and-urban-obesity/master/data/wrangled/wrangled_data.rda). Download this file and then place it in your current working directory within a subdirectory called "wrangled" within a subdirectory called "data" to copy and paste our code. We used an RStudio project and the [`here` package](https://github.com/jennybc/here_here) to navigate to the file more easily.

```{r}
load(here::here("data", "wrangled", "wrangled_data.rda"))
```


***
<details> <summary> Click here to see more about creating new projects in RStudio. </summary>

You can create a project by going to the File menu of RStudio like so:


```{r, echo = FALSE, out.width="60%"}
knitr::include_graphics(here::here("img", "New_project.png"))
```

You can also do so by clicking the project button:

```{r, echo = FALSE, out.width="60%"}
knitr::include_graphics(here::here("img", "project_button.png"))
```

See [here](https://support.rstudio.com/hc/en-us/articles/200526207-Using-Projects) to learn more about using RStudio projects and [here](https://github.com/jennybc/here_here) to learn more about the `here` package.

</details>
***
</details>
***


In the previous section, we hinted at the existence of statistical tests to be able to assess if there is a difference in the average BMI estimates between two groups (e.g. national BMI estimates in rural areas and national BMI estimates in urban areas). 
One such test is the [Student's $t$-test](https://stattrek.com/statistics/dictionary.aspx?definition=two-sample%20t-test){target="_blank"}, which can be used to determine if [two group means are different](http://onlinestatbook.com/2/estimation/difference_means.html){target="_blank"}.
However, as mentioned before, this test depends on certain assumption of our data. 

In this section, we are going to talk about this test and other tests that do not require the same assumptions about our data. 

### **Hypothesis testing**
***

Let's remind ourselves of one of our original questions, 

> Is there a difference between rural and urban BMI estimates around the world?

In hypothesis testing, we are interested in comparing two different hypotheses: a "null" hypothesis (can be thought of like a baseline e.g. the means between two groups are the same) compared to an "alternative" hypothesis (e.g. the means between two groups are different). We are going to ask if there is enough evidence in our data to reject the null hypothesis.  

Let's try to formalize this a bit. 

Suppose we want to test whether there is a difference between rural and urban BMI estimates using the data from 1985.  In our case, then, we want to test whether the mean BMI estimate among the rural areas is the same as the mean BMI estimate in the urban areas. If we call the true unknown means of the two groups $\mu_U$ and $\mu_R$, for the urban and rural areas, respectively, then we can define the <b>null hypothesis</b> that there is no difference in the two means:

$$ H_0: \mu_U = \mu_R $$  

In contrast, we also define an <b>alternative hypothesis</b> that there is a difference between the mean BMI estimates: 

$$ H_a: \mu_U \neq \mu_R $$

The idea behind a hypothesis test is that we assume the null hypothesis is true and we use our data to help us identify if there is enough evidence to _reject the null hypothesis_. 

This is similar to the idea of assuming that individuals are not guilty until proven otherwise.
If there is not enough evidence in the data, then we say we "fail to reject the null hypothesis".

You might be asking, "But we do not know what the true group means are?" 
That is correct. Instead, we must estimate these means based on the data that we have.
For example, in our data, we can estimate the average BMI estimates for women in rural and urban regions in 1985. Here we will use the `summarize()` function from the `dplyr` package again, this time using the argument `na.rm = TRUE` to remove NA values.

```{r}
BMI_long %>% 
  filter(Year == 1985, Sex == "Women", Region != "National") %>% 
  group_by(Region) %>%
  summarize(avg_bmi = mean(BMI, na.rm = TRUE)) 
```

Here we see that that in 1985, the average BMI estimate for women across urban areas is 24.6 compared to 23.6 across rural areas.  It looks like BMI estimates are, on average, higher in urban areas compared to rural areas.  However, now we want to use a statistical test to assess if this difference is statistically meaningful or a difference of this size could be found due to random chance. 

### **Two-sample $t$-test**
***

The two-sample $t$-test is a common way to determine if the means of two groups are different. It does this by calculating a observed test-statistic, $t$, that looks at the difference in the sample means of the two groups divided by an estimate of the variability in the difference in sample means.  Here $s_{pooled}$ is standard deviation of the data when pooling the samples together.
$$
t = \frac{\bar{X}_U-\bar{X}_R}{s_{pooled}\sqrt{\frac{1}{n_U} + \frac{1}{n_R}}}
$$
If the null hypothesis is true, and the means for urban and rural areas are the same, then the sample means in the numerator will be similar and this test statistic will be close to 0.  However, if the test statistic is not close to 0, this provides evidence the means for the urban and rural areas are difference, and so we may choose to reject the null hypothesis.  

To formally decide whether to reject the null hypothesis, we compare the observed test statistic to the distribution of possible test statistics that could occur if the null hypothesis were true using what we call a **p-value**.  The p-value tells us the probability of observing a test statistic as or more extreme than the one we observed if the null hypothesis is true.  If the p-value is small, say less than 0.05, then our observed test statistic is pretty rare if the null hypothesis is true.  This would lead us to reject the null hypothesis and conclude that the means of the two groups are different. We will talk more about p-values, and in particular our choice of 0.05 as a comparison for the p-value, a little later.

Here we perform a two-sample $t$-test comparing estimate for rural women in 1985 to estimates for urban women in 1985.  We use the `var.equal = TRUE` option here to use the pooled standard deviation $s_{pooled}$ from the formula above.
```{r}
t.test(pull(filter(BMI_long, Sex == "Women", 
                   Year == "1985", 
                   Region == "Rural"), BMI), 
       pull(filter(BMI_long, Sex == "Women", 
                   Year == "1985", 
                   Region == "Urban"), BMI), 
       var.equal = TRUE)
# means are different - p value <.05 reject the null: no difference in the means
```

In this case the p-value for the test is 0.00012, which is quite small.  Here would would reject the null hypothesis that the means in the two groups are equal and conclude there are differences in mean BMI estimates for women in urban and rural areas.  

But wait!  We need to check the assumptions of the statistical test we used to be sure it is valid.  The two-sample $t$-test relies on several <u><b>assumptions</b></u>:

1) The data for each group is [normally distributed](http://onlinestatbook.com/2/introduction/distributions.html){target="_blank"}.
2) The [variance](https://stattrek.com/statistics/dictionary.aspx?definition=variance){target="_blank"} of both groups is similar.
3) The observations from the two groups are [independent](https://www.stat.cmu.edu/~cshalizi/36-220/lecture-5.pdf){target="_blank"} (meaning that observations do not influence each other).
4) The observations within each group are [independent](https://www.stat.cmu.edu/~cshalizi/36-220/lecture-5.pdf){target="_blank"} (meaning that observations do not influence each other).


If these assumptions are violated, this doesn't necessarily mean we can't perform a statistical test. It just means we may need to consider the following options:

1) If the data isn't normally distributed, we may need to consider either a [transformation](https://www.statisticshowto.datasciencecentral.com/transformation-statistics/){target="_blank"} of the data to make it more normally distributed or a [nonparametric](https://www.mayo.edu/research/documents/parametric-and-nonparametric-demystifying-the-terms/doc-20408960){target="_blank"} test which doesn't rely on the normality assumption.
2) If the variances of the two groups aren't similar, we can perform a [Welch's $t$-test](https://www.statisticshowto.datasciencecentral.com/welchs-test-for-unequal-variances/){target="_blank"} (also called the unequal variance $t$-test) to account for difference variances.  This modifies the denominator of the test statistic to no longer use the standard deviation from pooling the samples.  We can do this in `R` by using the `var.equal = FALSE` option in the `t.test()` function.
3) If the data in the two samples are not independent and are instead what we call paired, we should use the [Paired $t$-test](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5579465/){target="_blank"}.
4) If the data within a group is not independent, we would need to consider methods beyond the scope of this case study.

Are these assumptions met in our case?  We've already seen, using Q-Q plots, that these BMI estimates may not be normally distributed.  We haven't checked the assumption of equal variance yet.  But it seems clear that the last assumption (of independent samples) is violated in our case!  The BMI data for rural areas comes from the same countries and the BMI data for urban areas, meaning the samples are not independent and are instead what we call [paired](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5579465/){target="_blank"}. We need to account for this because paired values may be more similar to one another because they come from the same country. Note that this would also be true for the male and female measurements from the same country or the values in the same countries from 1985 and later in 2017. 

### **Paired $t$-test**
***

When we perform a paired $t$-test, instead of looking at the means of the two samples separately, we combine observations in a pair (from the same country in our case) to calculate a difference in values between the two samples.  For example, we are really testing the `bmi_diff` column in the data below:

```{r}
BMI_long %>% 
  filter(Year == 1985, Sex == "Women", Region != "National") %>% 
  group_by(Country) %>%
  summarize(bmi_rural = BMI[Region == "Rural"],
            bmi_urban = BMI[Region == "Urban"],
            bmi_diff = bmi_rural - bmi_urban) 
```

This means our hypotheses are slightly different from the two-sample $t$-test. In this case we are testing the differences among the pairs of observations and how close these differences are to zero. (A mean difference of zero would mean the BMI estimates are the same in the rural and urban areas.) If we let $\mu_d$ represent the true mean difference between the paired observations of the two groups, then our null hypothesis is that the mean of the differences is equal to zero:
$$
                                        H_0 : \mu_d = 0
$$                                        

And the alternative hypothesis is that the mean of the differences is not equal to zero:
$$
                                        H_a : \mu_d \neq 0
$$                                      

We can perform a paired $t$-test by adding the option `paired = TRUE` to the `t.test()` command.  Note that we don't need to specify anything about equality of variances now, because there is only one variance, the variance of the paired differences!
```{r}
t.test(pull(filter(BMI_long, Sex == "Women", 
                   Year == "1985", 
                   Region == "Rural"), BMI), 
       pull(filter(BMI_long, Sex == "Women", 
                   Year == "1985", 
                   Region == "Urban"), BMI), 
       paired = TRUE)
# means are different - p value <.05 reject the null: no difference in the means
```

#### {.question_block}
<u>Question opportunity:</u> 
Looking at the t value, was global BMI lower in Rural or Urban areas in 1985?

####

Let's come back to the other $t$-test assumptions now.  With a paired $t$-test, we no longer need to assess the assumption about equal variance, since we are really testing the mean of one sample (of paired differences) rather than comparing two samples.  However, we still have to check normality and independence of observations.  We can rephrase these now to be about the paired differences rather than the two separate samples:

1) The paired differences are [normally distributed](http://onlinestatbook.com/2/introduction/distributions.html){target="_blank"}.
2) The paired differences are [independent](https://www.stat.cmu.edu/~cshalizi/36-220/lecture-5.pdf){target="_blank"} (meaning that observations do not influence each other).

If we assume that measurements between different countries are independent, this second assumption is not violated.  Therefore, we only now need to check normality of these paired differences.

We can start by making a Q-Q plot of the paired differences.
```{r}
BMI_long %>% 
  filter(Year == 1985, Sex == "Women", Region != "National") %>% 
  group_by(Country) %>%
  summarize(bmi_rural = BMI[Region == "Rural"],
            bmi_urban = BMI[Region == "Urban"],
            bmi_diff = bmi_rural - bmi_urban)  %>%
  ggplot(aes(sample = bmi_diff)) +
  stat_qq() +
  stat_qq_line()
```

We see that these paired differences do not appear to be normally distributed because the points do not lie on the straight line.  We can also see this by making a histogram of these paired differences:
```{r}
BMI_long %>% 
  filter(Year == 1985, Sex == "Women", Region != "National") %>% 
  group_by(Country) %>%
  summarize(bmi_rural = BMI[Region == "Rural"],
            bmi_urban = BMI[Region == "Urban"],
            bmi_diff = bmi_rural - bmi_urban)  %>%
  ggplot(aes(x = bmi_diff)) +
  geom_histogram()
```

This histogram appears to show show a bimodal distribution, meaning two groups of countries -- most countries have a rural BMI estimates lower than urban BMI estimates (`bmi_diff` < 0), while a smaller group of countries have urban BMI estimates higher than urban BMI estimates (`bmi_diff` >0).

First we will try transform our data to make it more normally distributed. One way to do this is to take the logarithm of the data values. Logarithms make extreme values less extreme, so can be helpful when data has a right skew.

```{r}
# create a new dataset with the paired differences
# create a log version of the difference variable using mutate()
BMI_diff_log <- BMI_long %>%
  filter(Year == 1985, Sex == "Women", Region != "National") %>% 
  group_by(Country) %>%
  summarize(bmi_rural = BMI[Region == "Rural"],
            bmi_urban = BMI[Region == "Urban"],
            bmi_diff = bmi_rural - bmi_urban) %>%  
  mutate(log_diff=log(bmi_diff))
```

Notice the warning message that `In log(bmi_diff) :  NaNs produced`.  This is because some of our differences were less than zero, and the logarithm of a negative number is not defined.  We can fix this issue by adding the same positive number, say 5, to each difference before taking the logarithm:
```{r}
BMI_diff_log <- BMI_long %>%
  filter(Year == 1985, Sex == "Women", Region != "National") %>% 
  group_by(Country) %>%
  summarize(bmi_rural = BMI[Region == "Rural"],
            bmi_urban = BMI[Region == "Urban"],
            bmi_diff = bmi_rural - bmi_urban) %>%  
  mutate(log_diff=log(bmi_diff + 5))
```

Let's look at a Q-Q plot and histogram of this transformed data:
```{r}
BMI_diff_log %>% 
  ggplot(aes(sample = log_diff)) +
  stat_qq() +
  stat_qq_line()

BMI_diff_log %>% 
  ggplot(aes(x = log_diff)) +
  geom_histogram()
```

The transformed data still does not appear to be normally distributed and the bimodal nature of the distribution is still apparent in the histogram.  What should we do?  The $t$-test is fairly **robust** to violations of the normality assumption if the sample size is relatively large, due to what is called the [central limit theorem](https://www.analyticsvidhya.com/blog/2019/05/statistics-101-introduction-central-limit-theorem/){target="_blank"}, which states that as samples get larger, the sample mean has an approximate normal distribution.  

With our sample size of 200, we could consider using the paired $t$-test even with our non-normal data.  However, if we had a smaller sample, we might still need to consider another alternative since the log-transformed data still doesn't look normal.  We might want to consider a non-parametric test in that case.

### **Nonparametric tests**
***

When the normality assumption of a paired $t$-test is violated, we can consider a nonparametric alternative instead.  What's a non-parametric test?  <b>Parametric tests</b> are based on assumptions about the distribution of the underlying data, while <b>nonparametric tests</b> do not rely on these distributional assumptions. 

See [here](https://www.mayo.edu/research/documents/parametric-and-nonparametric-demystifying-the-terms/doc-20408960){target="_blank"} for more information about the difference between these two classes of tests.

The nonparametric alternative to the paired $t$-test is called the [Wilcoxon signed-rank test](http://www.biostathandbook.com/wilcoxonsignedrank.html){target="_blank"}.  This test doesn't require the same normality assumption as the $t$-test, so can be considered with the data does not appear to be normally distributed and particularly when the number of samples is low. The null hypothesis for this test is that the **median** of the paired differences is 0.

To do a Wilcoxon test, we use the `wilcox.test()` function:
```{r}
wilcox.test(pull(filter(BMI_long, Sex == "Women", 
                        Year == "1985", 
                        Region == "Rural"), BMI),
           pull(filter(BMI_long, Sex == "Women", 
                       Year == "1985", 
                       Region == "Urban"), BMI),
           paired = TRUE)
```

With this nonparametric test, we still see a small p-value, indicating that BMI estimates between urban and rural areas are different. We would conclude that there's a difference in the median BMI estimates between these types of regions.

### **Multiple Testing Correction**
***

Let's take a look back at the questions we wanted to answer:
1) Is there a difference between rural and urban BMI estimates around the world? In particular, how does this look for females?
2) How have BMI estimates changed from 1985 to 2017? Again, In particular, how does this look for females?
3) How do different countries compare for BMI estimates? In particular, how does the United States compare to the rest of the world?

It's important to note that ultimately we wanted to test 4 different tests to answer question 1 and 2:

1) Is there a difference in mean BMI for women between Rural and Urban areas in 1985?
2) Is there a difference in mean BMI for women between Rural and Urban areas in 2017?
3) Is there a difference in mean BMI for women between 1985 and 2017 in Rural areas?
4) Is there a difference in mean BMI for women between 1985 and 2017 in Urban areas?

So far we have only looked at the first test.  Let's perform all 4 of these tests using the non-parametric Wilcoxon signed-rank test:
```{r}
wilcox.test(pull(filter(BMI_long, Sex == "Women", 
                        Year == "1985", 
                        Region == "Rural"), BMI),
           pull(filter(BMI_long, Sex == "Women", 
                       Year == "1985", 
                       Region == "Urban"), BMI),
           paired = TRUE)

wilcox.test(pull(filter(BMI_long, Sex == "Women", 
                        Year == "2017", 
                        Region == "Rural"), BMI),
           pull(filter(BMI_long, Sex == "Women", 
                       Year == "2017", 
                       Region == "Urban"), BMI),
           paired = TRUE)

wilcox.test(pull(filter(BMI_long, Sex == "Women", 
                        Year == "1985", 
                        Region == "Rural"), BMI),
           pull(filter(BMI_long, Sex == "Women", 
                       Year == "2017", 
                       Region == "Rural"), BMI),
           paired = TRUE)

wilcox.test(pull(filter(BMI_long, Sex == "Women", 
                        Year == "1985", 
                        Region == "Urban"), BMI),
           pull(filter(BMI_long, Sex == "Women", 
                       Year == "2017", 
                       Region == "Urban"), BMI),
           paired = TRUE)
```
In each of these tests, the very small p-value that would lead us to conclude that there is a difference between the two groups we are comparing.  However, when we test multiple questions like this, we need to correct for the multiple tests that we are performing.

The more we test data, the more likely we are to see a significant finding just by [random chance](https://www.stat.berkeley.edu/~mgoldman/Section0402.pdf){target="_blank"} even when the finding is in fact not real. It all comes down to probability.

#### Probability

You may recall that probability is the chance that an event will happen and it ranges from 0 to 1. 
You can also think of it as a percentage of 0% (probability =0) to 100% (probability = 1) of chance that an event will happen.

```{r, echo = FALSE, out.width='100%'}
include_graphics(here::here("img", "probability-line.svg"))

```



[[source](https://www.mathsisfun.com/data/images/probability-line.svg){target="_blank"}]

If we consider flipping a coin, the probability of an event happening (heads) is always equal to  1 ( or 100% chance) minus the probability of the event not happening (tails instead of heads) because the total probability can not exceed 100%.

        probability of 1 (100% chance) = the probability of heads (50% chance)
                                         + the probability of tails (50% chance)
                                       
                                     1 = P(heads) + P(tails)
            probability of one outcome = 1 - the probability of the other outcome
                              P(heads) = 1 - P(tails) 
                                        or
                              P(tails) = 1 - P(heads)

If we think about the probability of getting significant or nonsignificant results we can think of it like this:

                                     1 = P(significant results) 
                                          + P(no significant results)
                P(significant results) = 1 - P(no significant results)

### **Alpha**
***

Alpha is the amount of risk that you are willing to accept that a statistical test result is actually a false positive, in other words we reject the null hypothesis when the null hypothesis is true. This is also called [type 1 error](https://web.ma.utexas.edu/users/mks/statmistakes/errortypes.html){target="_blank"}. A significance threshold of .05 is customary and it means that we are accepting a 5% chance that our result is actually a false positive.

                                 alpha = probability of false significant results 
                                           or
                                     α = P(Making a type 1 error)
                                            
Remember that in our statistical tests if the p-value < alpha it is considered significant and if p-value > alpha it is considered not significant


But what really is the [p-value](https://towardsdatascience.com/p-values-explained-by-data-scientist-f40a746cfc8){target="_blank"}?

### **p-values**
***

The p in p-value stands for <b>probability</b> - the probability of observing a test statistic (for example the t in our student $t$-tests based on the means of our groups of comparison) as or more extreme than what we observe if the null hypothesis is true (continuing with the t test example, the distribution of possible t statistics if there is no difference between groups). Therefore a p-value of 0.02 means that there is a 2 percent chance that our data may also look the way it does if the null hypothesis is true. This does not mean there is a 2% chance that the null hypothesis is true. It only means that there is a 2% chance we would observe data as extreme as ours if there really is no difference between the groups. 

So then if alpha is the threshold for the p-value for obtaining false significant results then the probability of not making incorrect conclusions is 1 - alpha:

                 P(Not making a type 1 error) = 1 - α

```{r}
#P(Not making a type 1 error) = 1 - α
1-.05 

```
  
OK so if we use an alpha of .05 we are accepting that 95% of the time we will not make a Type 1 error and 5 % of the time we will just by random chance.              
Taking this one step further:

                     P(Making an error) = 1 - P(Not making an error)
                     P(Making an error) = 1 - (1 - α)
                                  
Here we can see that this checks out:                             
```{r}
#P(Making an error) = 1 - (1 - α)
1-(1-.05)
```

So what happens if we perform multiple tests?

The probability of not making a type 1 error would remain the same for each test. Therefore we need to multiply the probabilities together each time to determine the over all probability of making an error across multiple tests. See [here](http://mathforum.org/library/drmath/view/74065.html){target="_blank"} about why we multiply probabilities together.

        P(Not making an error in m tests) = (1 - α)^m
    P(Making at least 1 error in m tests) = 1 - (1 - α)^m

Let's consider if we performed 10 different statistical tests and if we as usual considered the significance threshold alpha of .05:

So the probability of getting 1 significant result with 10 tests is:                                   
```{R}
#P(Making at least 1 error in 10 tests) = 1 - ((1 - α)(1 - α)(1 - α)(1 - α)(1 - α)(1 - α)(1 - α)(1 - α)(1 - α)(1 - α))
#this is the same as:
#P(Making at least 1 error in 10 tests) = 1 - (1 - α)^10
 1- (1-.05)^10
```

So there is a 40% chance that that there will be a significant finding simply due to random chance alone.

What about 100 tests?

```{R}

#P(Making at least 1 error in 100 tests) = 1 - (1 - α)^100
 1- (1-.05)^100
```
                        
Yikes!! That is almost a 100% chance that there will be a significant finding simply due to chance alone!

Much of this explanation is described in this [lecture](https://www.gs.washington.edu/academics/courses/akey/56008/lecture/lecture10.pdf){target="_blank"}.


One way to correct for this is multiple testing issue is using the [Bonferroni method](http://mathworld.wolfram.com/BonferroniCorrection.html){target="_blank"}. 

In this method we would divide our significance threshold (generally 0.05) by the number of tests.

```{r}
.05/4 # 4 tests
```

Our new significance threshold is now 0.0125. 
Thus our p-values should be less than this value for us to reject the null that there is no difference in means.
In our 4 Wilcoxon tests, our p-values were less than 0.0125.
So we see a significant difference in the means of the groups after multiple testing correction for our four different tests.

## **Data Visualization**
***

If you have been following along but stopped, we could load our wrangled data like so:

```{r}
load(here::here("data", "wrangled", "wrangled_data.rda"))
```

***
<details> <summary> If you skipped the data wrangling section click here. </summary>

First you need to install and load the `OCSdata` package:

```{r, eval=FALSE}
install.packages("OCSdata")
library(OCSdata)
```

Then, you may load the wrangled data using the following code:

```{r, eval=FALSE}
wrangled_rda("ocs-bp-rural-and-urban-obesity", outpath = getwd())
load(here::here("OCSdata", "data", "wrangled", "wrangled_data.rda"))
```

If the package does not work for you, alternatively, an RDA file (stands for R data) of the data can be found [here](https://github.com//opencasestudies/ocs-bp-rural-and-urban-obesity/tree/master/data/wrangled) or slightly more directly [here](https://raw.githubusercontent.com/opencasestudies/ocs-bp-rural-and-urban-obesity/master/data/wrangled/wrangled_data.rda). Download this file and then place it in your current working directory within a subdirectory called "wrangled" within a subdirectory called "data" to copy and paste our code. We used an RStudio project and the [`here` package](https://github.com/jennybc/here_here) to navigate to the file more easily.

```{r}
load(here::here("data", "wrangled", "wrangled_data.rda"))
```


***
<details> <summary> Click here to see more about creating new projects in RStudio. </summary>

You can create a project by going to the File menu of RStudio like so:


```{r, echo = FALSE, out.width="60%"}
knitr::include_graphics(here::here("img", "New_project.png"))
```

You can also do so by clicking the project button:

```{r, echo = FALSE, out.width="60%"}
knitr::include_graphics(here::here("img", "project_button.png"))
```

See [here](https://support.rstudio.com/hc/en-us/articles/200526207-Using-Projects) to learn more about using RStudio projects and [here](https://github.com/jennybc/here_here) to learn more about the `here` package.

</details>
***
</details>
***


Again we will utilize `ggplot2` to create plots to look at the directional and magnitude of the differences in BMI we are interested in. If you need additional information please see the [Data Exploration section](https://www.opencasestudies.org/ocs-bp-rural-and-urban-obesity/#Distribution/#:~:text=Click%20here%20for%20additional%20information%20on%20using%20ggplot2.){target="_blank"} earlier in this case study. The top two lines of the code for the following plots, filter the data to only specific values of interest. Then we layer what is called a jitter on top of a box plot. A jitter is essentially a dot plot but with some variation on the location of the points so that they do not line up vertically which can make the individual points difficult to see.

Our first main question was: Is there a difference between rural and urban BMI estimates around the world?
Let's look at the national mean BMI estimates for each of the years:
```{r}
BMI_long %>% 
  filter(Sex == "Women", 
         Year == "1985", 
         Region != "National") %>%
  ggplot(aes(x = Region, y = BMI)) +
  geom_boxplot() +
  geom_jitter(width = .3) 
# the width determines how wide the points are plotted

BMI_long %>% 
  filter(Sex == "Women", 
         Year == "2017", 
         Region != "National") %>%
  ggplot(aes(x = Region, y = BMI)) +
  geom_boxplot() +
  geom_jitter(width = .3)

```

Our 2nd main question was: How have BMI estimates changed from 1985 to 2017?
Let's look at the change in rural and urban mean BMI estimates over time:
```{r}
#this time we will add titles for Rural and Urban using ggtitle()
BMI_long %>% 
  filter(Sex == "Women", 
         Year %in% c("1985", "2017"), 
         Region == "Rural") %>%
  # filtering this way allows us to keep data from both years
  ggplot(aes(x = Year, y = BMI)) +
  geom_boxplot() +
  geom_jitter(width = .3) +
  ggtitle("Rural")

BMI_long %>% 
  filter(Sex == "Women", 
         Year %in% c("1985", "2017"), 
         Region == "Urban") %>%
  ggplot(aes(x = Year, y = BMI)) +
  geom_boxplot() +
  geom_jitter(width = .3) +
  ggtitle("Urban")

```

Let's put the plots together to see how the change over the years differs between the regions.
We will use again use `facet_wrap()` to do this:

```{r}

BMI_long %>% 
  filter(Sex == "Women", 
         Year %in% c("1985", "2017"), 
         Region %in% c("Rural", "Urban")) %>%
  # filtering this way allows us to keep data from both years
  ggplot(aes(x = Year, y = BMI)) +
  geom_boxplot() +
  geom_jitter(width = .3) +
  facet_wrap(~ Region)

# we would also get national data without `Region %in% c("Rural", "Urban"))`

```
Indeed the change in BMI over time looks bigger in the rural areas!


One third main questions was: How do the different countries compare?  Or in other words what do the individual dots represent in our box plots?
We will take a look using `geom_label()`:

```{r}

BMI_long %>% 
  filter(Sex == "Women", 
         Year %in% c("1985", "2017"), 
         Region == "Rural") %>%
  ggplot(aes(x = Year, y = BMI)) +
  geom_boxplot() +
  geom_jitter(width = .3) +
  geom_label(aes(label = Country))

```

If we include all country names this is a bit too much... so perhaps we should focus on just the extreme BMI values using `filter()` function of the`dplyr` package. We will also use the `ggrepel` package to have our labels not overlap each other.

Part of the third question was: How does the United States of America compare? So let's label the data points for the United states.

```{r}

BMI_long %>% 
  filter(Sex == "Women", 
         Year %in% c("1985", "2017"), 
         Region == "Rural") %>%
  ggplot(aes(x = Year, y = BMI)) +
  geom_boxplot() +
  geom_jitter(data=BMI_long %>% 
              filter(Sex == "Women", 
                     Year %in% c("1985", "2017"), 
                     Region == "Rural") %>%
              filter(BMI>31 | BMI<19 | Country == "United States of America"), 
              aes(x =Year, y =BMI),
              width = .02) +
  ggrepel::geom_text_repel(data=BMI_long %>% 
              filter(Sex == "Women", 
                     Year %in% c("1985", "2017"), 
                     Region == "Rural") %>%
              filter(BMI>31 | BMI<19 | Country == "United States of America"), 
              aes(x =Year, y =BMI, label = Country))
              
              
```

And let's fill the box plots with color and outline in black:
```{r}

BMI_long %>% 
  filter(Sex == "Women", 
         Year %in% c("1985", "2017"), 
         Region == "Rural") %>%
  ggplot(aes(x = Year, y = BMI)) +
  geom_boxplot(outlier.shape = NA, color = "black", aes(fill = Year)) +
  geom_jitter(data=BMI_long %>% 
              filter(Sex == "Women", 
                     Year %in% c("1985", "2017"), 
                     Region == "Rural") %>%
              subset(BMI>31 | BMI<19 | Country == "United States of America"), 
              aes(x =Year, y =BMI),
              width = .02) +
  ggrepel::geom_text_repel(data=BMI_long %>% 
              filter(Sex == "Women", 
                     Year %in% c("1985", "2017"), 
                     Region == "Rural") %>%
              subset(BMI>31 | BMI<19 | Country == "United States of America"), 
              aes(x =Year, y =BMI, label = Country))

```


### **Overall differences**
***

Let's take a look at all the data together, including the data for men:
```{r}
ggplot(BMI_long, aes(x = Year, y = BMI, col = Region)) +
  geom_boxplot(outlier.shape = NA, color = "black" , aes(fill = Region)) + 
  facet_grid(~ Sex)   
```

That's useful, but let's look at the individual points and include our country labels, to do so lets change our United States of America label to USA:

```{r}
BMI_long[["Country"]] <-BMI_long[["Country"]] %>%
  str_replace( pattern = "United States of America", replacement = "USA")


ggplot(BMI_long, aes(x = Year, y = BMI, col = Year)) + 
  geom_boxplot(outlier.shape = NA, color = "black" , aes(fill = Sex)) + 
  facet_grid(~ Sex + Region) + 
  geom_hline(yintercept=30, linetype="dashed", color = "red", size =1) + 
  # This will add a horizontal dashed line to indicate the obesity BMI threshold
  geom_jitter(data=subset(BMI_long), 
              aes(x =Year, y =BMI), 
              width = .2, size =2, shape =21, color = "black", fill = "gray") + 
  # This will add the individual country data points
  # The shape 21 allows for a different fill and outline color
  # The width determines how wide the jitter points are plotted
  geom_jitter(data=subset(BMI_long, Country == "USA"), 
              aes(x =Year, y =BMI), 
              width = .02, size =12, shape =21, color = "black", fill = "gray") + 
  # This will add points that are larger for the USA data
  geom_text(data=subset(BMI_long,  Country == "USA"), 
            aes(x =Year, y =BMI,label=Country), 
            color = "black") + 
  # This will add USA labels to the USA points
  theme(legend.position = "none", 
  # This is useful for removing the legend
        axis.text.x = element_text(size = 15,angle = 30), 
        # this changes the size and angle of the x axis point labels 
        axis.text.y = element_text(size = 20), 
        axis.title.y = element_text(size =15), 
        axis.title.x = element_text(size =15), 
        strip.text.x = element_text(size = 15)) +
  # This changes the size of x axis labels for the facet 
        ggtitle( "Differences in BMI Over Time and Across Region Type and Gender") +
  # Add a  plot title
        scale_fill_manual(values=c("dodgerblue", "orchid2"))

```

Here we can see that overall BMI appears to be increasing globally over time. Additionally we can see that this is occurring not just in urban areas, but also in rural areas. The US is consistently above the median in all strata of the data. In general, the female data shows higher BMI values than the male data. The rural USA BMI appears to be higher than the urban BMI for both men and women. Many countries have average BMI estimates above the obesity threshold of 30. Thus it appears that education and outreach programs for weight management should focus on both rural and urban areas and both genders. Education and assistance for women may be especially helpful. 

### **Differences in rate of change**
***

How does the rate of change in BMI differ between groups? Which group might especially need attention?

First let's calculate the differences in BMI from 2017 and 1985 and add this to our BMI_long data object. This will create new variables called `Rural_difference,` `Urban_difference`, and `National_difference`.
```{r}
BMI_numeric  %<>% 
  mutate(Rural_difference = Rural_BMI_2017 - Rural_BMI_1985) %>%
  mutate(Urban_difference = Urban_BMI_2017 - Urban_BMI_1985) %>%
  mutate(National_difference = National_BMI_2017 - National_BMI_1985)

BMI_diff_long <- BMI_numeric %>% 
  select(Country:National_difference) %>% 
  gather(key = Type, value = Difference , Rural_difference:National_difference)


BMI_diff_long <- BMI_numeric %>% 
  select(Country: National_difference) %>% 
  pivot_longer( names_to = "Difference" , cols = Rural_difference:National_difference)


head(BMI_diff_long)
```

Let's replace "United states of America" with "USA" and make a plot with this data to compare the change in BMI:

```{r}

BMI_diff_long[["Country"]] <-BMI_diff_long[["Country"]] %>%
  str_replace( pattern = "United States of America", replacement = "USA")


ggplot(BMI_diff_long, aes(x = Difference, y = value, col = Difference)) + 
  geom_boxplot(outlier.shape = NA, color = "black" , aes(fill = Difference)) + 
  facet_grid(~ Sex) +
  geom_jitter(data=subset(BMI_diff_long), 
              aes(x =Difference, y =value), 
              width = .2, size = 2, shape = 21, color = "black", fill = "gray") + 
  # This will add the individual country data points
  # The shape 21 allows for a different fill and outline color
  # The width determines how wide the jitter points are plotted
  geom_jitter(data=subset(BMI_diff_long, Country == "USA"), 
              aes(x = Difference, y = value), 
              width = .02, size = 12, shape = 21, color = "black", fill = "gray") + 
  # This will add points that are larger for the USA data
  geom_text(data=subset(BMI_diff_long,  Country == "USA"), 
            aes(x =Difference, y = value, label = Country), color = "black") + 
  # This will add USA labels to the USA points
  theme(legend.position = "none", 
  # This is useful for removing the legend
        axis.text.x = element_text(size = 15,angle = 30), 
        # this changes the size and angle of the x axis point labels 
        axis.text.y = element_text(size = 20), 
        axis.title.y = element_text(size = 15), 
        axis.title.x = element_text(size = 15), 
        strip.text.x = element_text(size = 15)) +
        # this changes the size of x axis labels for the facet
        ggtitle( "Change in BMI Over Time and Across Region Type and Gender")

```

We can now see that the rate of change from 1985 to 2017 appears to be larger in the women compared to the men in all regions. The group with the largest increase in the USA is the women living in rural areas.



Let's check the difference with some statistical tests.  Remember that these rural differences are measured on the same countries as the urban differences, so we will need to use paired tests to make these comparisons.  We can either consider the paired $t$-test or the Wilcoxon signed-rank test.  To decide between them, we'll look at the normality of the paired differences in change over time.

Let's compare the change over time in rural communities to the change in urban communities among women and separately among men.

```{r}
# lets look at the normality of these paired differences
# first we need to calculate the paired differences of change over time
BMI_diff_long %>%
  filter(Difference != "National_difference") %>%
  group_by(Country, Sex) %>%
  summarize(change_rural = value[Difference == "Rural_difference"],
            change_urban = value[Difference == "Urban_difference"],
            change_diff = change_rural - change_urban) %>%
  ggplot(aes(x=change_diff)) +
  geom_histogram() +
  facet_wrap(~ Sex)
# interesting...the spread of the difference in change over time is much 
# broader for women than for men
# for women, the differences look reasonably normal, spread between -1 and +1
# for men there appear to be two groups of countries; one with difference in change over time
# centered around 0 and the other with difference in change over time between 1 and 2.


# Let's look at Q-Q plots:
BMI_diff_long %>%
  filter(Difference != "National_difference") %>%
  group_by(Country, Sex) %>%
  summarize(change_rural = value[Difference == "Rural_difference"],
            change_urban = value[Difference == "Urban_difference"],
            change_diff = change_rural - change_urban) %>%
  ggplot(aes(sample = change_diff)) +
  stat_qq() + # to add the line we need to also include stat_qq_line()
  stat_qq_line() +
  facet_wrap(~ Sex)

# the distribution of differences in change over time does not appear normal
# for either men or women
```

Based on the non-normality of these paired differences, we will use a Wilcoxon signed-rank test to compare changes over time between urban and rural areas separately for men and women:
```{r}
wilcox.test(pull(filter(BMI_diff_long, 
                    Sex == "Women",  
                    Difference == "Rural_difference"), value), 
        pull(filter(BMI_diff_long,
                    Sex == "Women",  
                    Difference == "Urban_difference"), value),
        paired = TRUE)

wilcox.test(pull(filter(BMI_diff_long,
                    Sex == "Men",
                    Difference == "Rural_difference"), value), 
        pull(filter(BMI_diff_long,
                    Sex == "Men", 
                    Difference == "Urban_difference"), value),
        paired = TRUE)
```
Using a p-value threshold of 0.05, we would conclude there's a difference in change over time between urban and rural areas for women (p = 0.00327) but not for men (p = 0.08256).  

However, now we have performed six comparisons (4 earlier) so we should apply our multiple testing correction and use a new p-value threshold of:
```{r}
.05/6
```

With the threshold adjusted for multiple comparisons, we would still conclude there's a difference in the change over time between urban and rural areas for women but not for men.

Interestingly, although we don't see a difference in change over time for men, our histogram of the differences between urban and rural areas for men highlighted two groups of countries -- one centered around zero (no difference) and one above 1 (showing a difference):  
```{r}
BMI_diff_long %>%
  filter(Difference != "National_difference") %>%
  group_by(Country, Sex) %>%
  summarize(change_rural = value[Difference == "Rural_difference"],
            change_urban = value[Difference == "Urban_difference"],
            change_diff = change_rural - change_urban) %>%
  ggplot(aes(x=change_diff)) +
  geom_histogram() +
  facet_wrap(~ Sex)
```

In fact for these countries, the difference between rural and urban areas for men is higher than the difference for women for almost all of the countries!

We can look to see which countries have such a large difference in change over time between rural and urban areas for men:
```{r}
BMI_diff_long %>%
  filter(Sex == "Men", Difference != "National_difference") %>%
  group_by(Country) %>%
  summarize(change_rural = value[Difference == "Rural_difference"],
            change_urban = value[Difference == "Urban_difference"],
            change_diff = change_rural - change_urban) %>%
  filter(change_diff > 1)
```


For women, we did see a difference in the change over time between urban and rural areas.  There appear to be specific countries where BMI shows a particular increase especially for women. Which countries are those? How does that compare with the US? Clearly the US is among the countries with the highest differences. 

```{r}
BMI_diff_long %>% 
  filter(Country == "USA")

```

In the US, focus should be placed on <b> both urban and rural </b> women to improve this public health issue. 

Here we can see the countries that have the largest differences in BMI from 1985-2017:

```{r}
BMI_diff_long %>%
  filter(value > 4.9) %>%
  arrange(- value)
#here we sorted by the highest to lowest values of BMI Difference
```

However it is important to see what the mean BMI values were for these countries in 2017. 
It could be that the average was underweight in 1985... let's take a look.

```{r}
BMI_long %>% 
  filter(Country %in% pull(filter(BMI_diff_long, value > 4.9), Country), 
                    Sex %in% pull(filter(BMI_diff_long, value > 4.9), Sex),
                    Year =="2017") %>%
  arrange(- BMI)

#Here is the data just for the changes in BMI in Egypt
filter(BMI_diff_long, Country == "Egypt") %>% 
  arrange(- value)

```

Thus rural women in Egypt had the greatest increase in BMI from 1985 to 2017 in this data (5.9) and the average BMI is now over the obesity threshold of 30.

The data suggests that rural women in Egypt and other countries may especially benefit from dietary resources and our interventions and programs to assist with weight management.



## **Summary**
***

### **Summary Plot**
***
Now let's make a plot that summarizes our findings. To do this we will simplify the other plots we made and then combine them together.

```{r}

# simplified national means plot
Means_plot<-BMI_long %>% 
  filter(Sex %in% c("Men", "Women"), 
         Year %in% c("1985", "2017"), 
         Region == "National") %>%
ggplot(aes(x = Sex, y = BMI)) + 
  geom_boxplot(outlier.shape = NA, color = "black" , aes(fill = Sex)) + 
 scale_fill_manual(values=c("dodgerblue", "orchid2")) +
  facet_grid(~ Year) + 
  geom_hline(yintercept=30, linetype="dashed", color = "red", size =1) + 
  geom_jitter(data=BMI_long %>%
                filter(Sex %in% c("Men", "Women"), 
                       Year %in% c("1985", "2017"), 
                       Region == "National"),
                aes(x =Sex, y =BMI), 
                width = .2, size =2, shape =21, 
                color = "black", fill = "gray") +
  geom_jitter(data=subset(BMI_long %>%
                            filter(Sex %in% c("Men", "Women"), 
                                   Year %in% c("1985", "2017"), 
                                   Region == "National", 
                                   Country == "USA")), 
                             aes(x =Sex, y =BMI), 
                             width = .02, size =12, shape =21, 
                             color = "black", fill = "gray") + 
  # This will add points that are larger for the USA data
  geom_text(data=subset(BMI_long%>%
                          filter(Sex %in% c("Men", "Women"), 
                                 Year %in% c("1985", "2017"), 
                                 Region == "National", Country == "USA")),
                          aes(x =Sex, y =BMI,label=Country), 
                          color = "black") + 
  # This will add USA labels to the USA points
  theme_linedraw() +
  # This will make the background of the plot white and the facet labels black
  theme(legend.position = "none", 
  # This is useful for removing the legend
        axis.text.x = element_text(size = 15,angle = 30, vjust = 0.5), 
        # this changes the size and angle of the x axis point labels 
        axis.text.y = element_text(size = 20), 
        axis.title.y = element_text(size =15), 
        axis.title.x = element_text(size =15), 
        strip.text.x = element_text(size = 15))
   
#Simplified difference plot 
BMI_diff_long[["Difference"]] <-BMI_diff_long[["Difference"]]%>%
  str_replace( pattern = "_difference", replacement = "")

Diff_plot<-BMI_diff_long %>% 
  filter(Sex %in% c("Men", "Women"), 
         Difference != "National") %>%
ggplot( aes(x = Difference, y = value, col = Difference)) + 
  geom_boxplot(outlier.shape = NA, color = "black" , aes(fill = Sex)) + 
  scale_fill_manual(values=c("dodgerblue", "orchid2")) +
  facet_grid(~ Sex) +
  geom_jitter(data=BMI_diff_long %>% 
    filter(Sex %in% c("Men", "Women"), 
           Difference != "National"), 
    aes(x =Difference, y = value), 
    width = .2, size = 2, shape = 21, 
    color = "black", fill = "gray") + 
  geom_jitter(data=subset(BMI_diff_long %>% 
    filter(Sex %in% c("Men", "Women"), 
           Difference != "National", 
           Country == "USA")), 
    aes(x = Difference, y = value), 
    width = .02, size =12, shape =21, 
    color = "black", fill = "gray") + 
  # This will add points that are larger for the USA data
  geom_text(data=subset(BMI_diff_long %>% 
    filter(Sex %in% c("Men", "Women"),
           Difference != "National"), 
           Country == "USA"),
    aes(x = Difference, y = value,label=Country), 
    color = "black") + 
  # This will add USA labels to the USA points
  theme_linedraw() +
    # This will make the background of the plot white and the facet labels black
  theme(legend.position = "none", 
  # This is useful for removing the legend
        axis.text.x = element_text(size = 15,angle = 30, vjust = 0.5), 
        # this changes the size and angle of the x axis point labels 
        axis.text.y = element_text(size = 20), 
        axis.title.y = element_text(size =15), 
        axis.title.x = element_text(size =15), 
        strip.text.x = element_text(size = 15))
  # this changes the size of x axis labels for the facet

#add labels
Diff_plot<-Diff_plot + 
        labs(title = "Change in BMI by region", 
             x = "", 
             y = "Change in BMI \n (1985 to 2017)") +
        theme(title = element_text (size = 12, face = "bold"))

Means_plot <-Means_plot + 
        labs(title = "Mean National BMI over time", 
             x = "", 
             y = "Mean BMI") +
        theme(title = element_text (size = 12, face = "bold"))



obesity_text<-tibble(Year=c(1985),BMI=c(31),Sex=c("Men"),label=c("Obesity"))

Means_plot <-Means_plot + 
             geom_text(data = obesity_text,
                       label=pull(obesity_text,label), 
                       color = "red", 
                       aes( fontface ="bold.italic", size = 13))

cowplot::plot_grid(Means_plot, Diff_plot, labels = c("A", "B"))
```

We could make a similar plot with the `patchwork` package:
```{r}
Means_plot +labs(tag = 'A') + Diff_plot+ labs(tag = 'B') +
  plot_layout(guides = 'collect', ncol = 2)
```

Great, now we have put two plots together using the `plot_grid()` function of the `cowplot` package. This way we can clearly communicate two messages. The first being that BMI has increased over time globally and that many countries including the United States of America are approaching a mean BMI that is above the obesity threshold of 30. We can also see that women on average have larger BMI values than males, but that both genders show increased levels over time. In the second plot we can see that the increase in BMI is not just happening in  urban communities, but in both rural and urban communities and particularly in women. 

Our plot visually explores all of our main questions:
1) Is there a difference between rural and urban BMI estimates around the world? In particular, how does this look for females?
2) How have BMI estimates changed from 1985 to 2017? Again, In particular, how does this look for females?
3) How do different countries compare for BMI estimates? In particular, how does the United States compare to the rest of the world?

Now let's save our plot as a png file using the `png()` function:

```{r}
png(here::here("img", "main_plot.png"), height = 1500, width = 2300, res = 300)
cowplot::plot_grid(Means_plot, Diff_plot, labels = c("A", "B"))
dev.off()
```

You might think that these two plots look a bit odd together, that if you look at the change in the boxplots in part A, that the data for Men shows a bigger change, however, we can't tell from this plot how the individual countries are changing. Part B of the figure looks at the specific change for each country.

We could check our plot using the following code:
```{r}
BMI_diff_long %>%
  group_by(Sex, Difference) %>%
  filter(!is.na(value))%>%
  filter(Country == "USA") %>%
  summarize(mean = mean(value), median = median(value))

BMI_diff_long %>%
  group_by(Sex, Difference) %>%
  filter(!is.na(value))%>%
  summarize(mean = mean(value), median = median(value))

BMI_long %>%
  group_by(Sex, Year, Region) %>%
  filter(!is.na(BMI))%>%
  filter(Country == "USA") %>%
  summarize(mean = mean(BMI), median = median(BMI))

BMI_long %>%
  group_by(Sex, Year, Region) %>%
  filter(!is.na(BMI))%>%
  summarize(mean = mean(BMI), median = median(BMI))

BMI_numeric %>%
  group_by(Sex)%>%
  summarize(Nat_mean1985 = mean(National_BMI_1985), 
            Nat_mean2017 = mean(National_BMI_2017))

BMI_numeric %>%
  group_by(Sex)%>%
  summarize(Nat_med1985 = median(National_BMI_1985), 
            Nat_med2017 = median(National_BMI_2017))
```


We can see that the largest changes occur for Women.

Scroll through the output!

#### {.scrollable }
```{r}
BMI_diff_long %>%
  filter(value>3.5) %>%
  arrange(-value) %>%
  # this allows us to show the full output
  print(n = 1e3)
```
####

### **Synopsis
***

We have evaluated BMI average estimates from 200 different countries around the world. To do so we imported data from a pdf using the `pdftools` package. We used `tidyverse` packages such as `dplyr`, `stringr`, and `tidy` to clean the data and get it in a workable format. Our statistical analysis focused on evaluating differences in BMI in females around the world across time and between rural and urban areas. We found a significant difference both between years and between the type of community among women globally using [paired $t$-tests](https://stattrek.com/hypothesis-test/paired-means.aspx#example1){target="_blank"} and [nonparametric](https://www.mayo.edu/research/documents/parametric-and-nonparametric-demystifying-the-terms/doc-20408960){target="_blank"} tests. Thus BMI has increased in women since 1985. Although BMI estimates are significantly higher in Urban areas compared to rural areas, BMI estimates have increased in both regions. We learned that there are assumptions and considerations to keep in mind when performing tests that compare two groups. 

<u> We learned that the two sample $t$-test relies on:</u>

1) normality of the data for both groups (this is not as much of an issue if the number of observations is relatively large total n>30)
2) equal variance between the two groups (make sure you do the correct test if the data is not normal)
3) independent observations

We learned that we can evaluate if our data is [normally distributed](https://www.physiology.org/doi/full/10.1152/advan.00064.2017){target="_blank"} by plotting the distribution and by creating Q-Q plots.

<u>We learned that if our data is not normally distributed, we can consider these options:</u>

1) We can still perform a $t$-test if our n is large
2) We can transform the data before performing a $t$-test

We learned that if our data has unequal variance, we can consider using the [Welch's $t$-test](https://www.statisticshowto.datasciencecentral.com/welchs-test-for-unequal-variances/){target="_blank"}.


Importantly, we learned that when we have [paired data](https://www.statisticshowto.datasciencecentral.com/paired-data/){target="_blank"} (are observations are not independent) and  we intend to compare means, we should perform a paired test ([paired $t$-test](https://stattrek.com/hypothesis-test/paired-means.aspx#example1){target="_blank"} or the nonparametric Wilcoxon signed rank test. Our data used in this case study was paired because observations were taken for the same countries for different categories of populations - thus, we wanted to compare these populations within the same country (male vs female, rural vs urban etc.). Other examples would be in a case-control study (if samples are matched for various demographics) or a study with repeated measures (for example, symptom measures from the same individual before and after a treatment). We learned that [paired $t$-test](https://stattrek.com/hypothesis-test/paired-means.aspx#example1){target="_blank"} tests a different hypothesis the two-sample $t$-test. In this case we are testing the differences among the pairs of observations and how close these differences are to zero, thus there is only one variance, the variance of the paired differences and therefore, we do not need to worry about the equality of variance assumption like we do with the two-sample $t$test.
We learned that we can perform a paired $t$-test by adding the option `paired = TRUE` to the `t.test()` command.

The paired $t$-test assumes that:

1) The paired differences are [normally distributed](http://onlinestatbook.com/2/introduction/distributions.html){target="_blank"}, which we can evaluate using Q-Q plots or plotting the distribution of the differences.
2) The paired differences are [independent](https://www.stat.cmu.edu/~cshalizi/36-220/lecture-5.pdf){target="_blank"} (meaning that observations do not influence each other).

We learned that if the normality assumption is violated, we can use the Wilcoxon signed rank test, or if our data is large enough, due to the central limit theorem, the paired $t$-test should still be robust to deviations from normality.

Using the `ggplot2` package we were able to visualize trends in the data. Importantly, the largest increase appears to be in the rural areas particularly for women. We were also able to identify how the US compared to other countries and which countries showed the largest increase in BMI over time. Analyses like this are important for defining which groups could benefit the most from interventions, education, and policy changes when attempting to mitigate public health challenges. You can see in the [article](https://www.nature.com/articles/s41586-019-1171-x.pdf){target="_blank"} that this data came from that many additional considerations would be involved to perform such an analysis.


## **Suggested Homework**
***

Students can evaluate the change in BMI over time using the global data available for each year between 2015 and 2017. The data can be found [here](http://www.ncdrisc.org/downloads/bmi/NCD_RisC_Lancet_2017_BMI_age_standardised_world.csv){target="_blank"}.

## **Additional Information**
***

### **Helpful Links**
***

[BMI](https://www.cdc.gov/healthyweight/assessing/bmi/adult_bmi/index.html){target="_blank"}  
[Long and Wide Data Formats](https://opencasestudies.github.io/ocs-healthexpenditure/ocs-healthexpenditure.html){target="_blank"}  
[Distributions](http://onlinestatbook.com/2/introduction/distributions.html){target="_blank"}
[Normal Distribution](http://onlinestatbook.com/2/introduction/distributions.html){target="_blank"}
[Skewed Distributions](http://onlinestatbook.com/2/glossary/skew.html){target="_blank"} 
[Bimodal Distribution](http://onlinestatbook.com/2/introduction/distributions.html){target="_blank"} 
[ggplot2](https://opencasestudies.github.io/ocs-healthexpenditure/ocs-healthexpenditure.html){target="_blank"}    
[Q-Q Plots](http://onlinestatbook.com/2/advanced_graphs/q-q_plots.html){target="_blank"}  
[Student $t$-test](https://stattrek.com/statistics/dictionary.aspx?definition=two-sample%20t-test)   
[Paired Data](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5579465/){target="_blank"}  
[Welch's $t$-test](https://www.statisticshowto.datasciencecentral.com/welchs-test-for-unequal-variances/){target="_blank"}    
[Parametric and Nonparametric Methods](https://www.mayo.edu/research/documents/parametric-and-nonparametric-demystifying-the-terms/doc-20408960)  
[Variance](https://stattrek.com/statistics/dictionary.aspx?definition=variance){target="_blank"}  
[Balanced Study Design](https://www.statisticshowto.datasciencecentral.com/balanced-and-unbalanced-designs/){target="_blank"}  
[Independent Observations](https://www.stat.cmu.edu/~cshalizi/36-220/lecture-5.pdf){target="_blank"}  
[Transformation](https://www.statisticshowto.datasciencecentral.com/transformation-statistics/){target="_blank"}  
[Permutation/Resampling Methods](https://jhu-advdatasci.github.io/2019/lectures/21-resampling-techniques.html){target="_blank"}   
[Central Limit Theorem](https://www.analyticsvidhya.com/blog/2019/05/statistics-101-introduction-central-limit-theorem/){target="_blank"}  
[Mood's Two-Sample Scale Test](https://files.eric.ed.gov/fulltext/ED065559.pdf){target="_blank"}  
[Wilcoxon Signed Rank Test](http://www.biostathandbook.com/wilcoxonsignedrank.html)   
[Wilcoxon Rank Sum Test](http://sphweb.bumc.bu.edu/otlt/mph-modules/bs/bs704_nonparametric/BS704_Nonparametric4.html){target="_blank"}  
[Two-sample Kolmogorov-Smirnov Test](https://www.itl.nist.gov/div898/software/dataplot/refman1/auxillar/ks2samp.htm){target="_blank"}  
[Type 1 Error](https://web.ma.utexas.edu/users/mks/statmistakes/errortypes.html){target="_blank"}  
[p-value](https://towardsdatascience.com/p-values-explained-by-data-scientist-f40a746cfc8){target="_blank"}  
[Multiple Testing](https://www.gs.washington.edu/academics/courses/akey/56008/lecture/lecture10.pdf){target="_blank"}    
[Bonferroni Method of Multiple Testing Correction](http://mathworld.wolfram.com/BonferroniCorrection.html){target="_blank"}

<u>**Packages used in this case study:** </u>

 Package   | Use in this case study                                                                         
---------- |-------------
[here](https://github.com/jennybc/here_here){target="_blank"}       | to easily load and save data with relative paths
[pdftools](https://cran.r-project.org/web/packages/pdftools/pdftools.pdf){target="_blank"}   | to read a text from pdf into R   
[stringr](https://stringr.tidyverse.org/articles/stringr.html){target="_blank"}    | to manipulate the text data
[readr](https://readr.tidyverse.org/){target="_blank"}      | to manipulate the text data within the pdf into individual lines  
[dplyr](https://dplyr.tidyverse.org/){target="_blank"}      | to arrange/filter/select subsets of the data 
[tibble](https://tibble.tidyverse.org/){target="_blank"}     | to create data objects that we can manipulate with `dplyr`/`stringr`/`tidyr`/`purrr`
[magrittr](https://magrittr.tidyverse.org/articles/magrittr.html){target="_blank"}   | to use the `%<>%` piping operator
[glue](https://www.tidyverse.org/blog/2017/10/glue-1.2.0/){target="_blank"}  | to paste or combine character strings and data together
[purrr](https://purrr.tidyverse.org/){target="_blank"}      | to perform functions on all columns of a tibble
[tidyr](https://tidyr.tidyverse.org/){target="_blank"}      | to convert data from 'wide' to 'long' format
[ggplot2](https://ggplot2.tidyverse.org/)    | to make visualizations with multiple layers
[ggrepel](https://cran.r-project.org/web/packages/ggrepel/vignettes/ggrepel.html){target="_blank"}    | to allow labels in figures not to overlap
[cowplot](https://cran.r-project.org/web/packages/cowplot/vignettes/introduction.html){target="_blank"} and [patchwork](https://github.com/thomasp85/patchwork){target="_blank"} | to allow plots to be combined 


### **Session info**
***

```{r}
sessionInfo()
```

**Estimate of RMarkdown Compilation Time: **

```{r, echo=FALSE}
rmarkdown:::perf_timer_stop("render")
pts = rmarkdown:::perf_timer_summary()
cat("About", round(pts$time[1]/1000 + 5), "-", round(pts$time[1]/1000 + 15),"seconds")
```

This compilation time was measured on a PC machine operating on Windows 10. This range should only be used as an estimate as compilation time will vary with different machines and operating systems. 

### **Acknowledgments**
***

We would like to acknowledge [Jessica Fanzo](https://www.jhsph.edu/faculty/directory/profile/3380/jessica-fanzo) for assisting in framing the major direction of the case study.

We would like to acknowledge [Michael Breshock](https://mbreshock.github.io/) for his contributions to this case study and developing the `OCSdata` package.

We would also like to acknowledge the [Bloomberg American Health Initiative](https://americanhealth.jhu.edu/) for funding this work. 

<script type='text/javascript' id='clustrmaps' src='//cdn.clustrmaps.com/map_v2.js?cl=080808&w=a&t=tt&d=a1Ews2eHEWlY8v8IFJHWTERegYZfMmtyCxPiYu5y9DM&co=ffffff&cmo=3acc3a&cmn=ff5353&ct=808080'></script>