MATH2349 Semester 1, 2020. 2016 Suicide Rate by Country and Gender.

Assignment 2

Required packages

library(dplyr)
library(tidyr)
library(ggplot2)
library(MVN)
library(forecast)

Executive Summary

This report describes the pre-processing of two data sets: population projection (United Nations 2019) and suicide rate estimates (World Health Organization 2018).

Initially, unused columns and rows are filtered on both data sets. This report only uses the data for the year 2016 and data for other years is removed. The population data set is found to be untidy, so the original wide data format is converted into a long format.

There are a few countries that have slightly different names in the two data sets (e.g. “North Macedonia” vs. “Republic of North Macedonia”. So the name of these countries are manually updated to match so that the two data sets could be merged by country and gender. Then, the gender variable is converted from character into an unordered factor with two possible values: male and female.

A new variable for the estimate number of suicides per country per gender is calculated from the suicide rate and number of populations. The data frame is then checked for missing values where it is determined that no missing values are found.

The data frame has three numeric variables: suicide rate, population, and number of suicide estimates. The distribution of these three variables is visualised through box plots, where it is determined that all three have a right skewed distribution. Hence, the Tukey’s method is used to look for outliers. However, the determined outliers do not seem to come from data entry errors. They seem to come from countries having extreme values of suicide rate and population (e.g. China and India), so these outliers are not imputed and are kept.

Finally, a Box-Cox transformation is performed on the suicide rate variable to try to convert the distribution into a normal distribution.

Data

The first data set is the United Nation population projection by location and gender from the year 1950 to 2100 (United Nations 2019). The data set has several projection variants where each uses different assumptions on projecting future levels of fertility, mortality, and migration.

setwd("~/Dropbox/data-wrangling/assignment2")

population <- read.csv("WPP2019_TotalPopulationBySex.csv", stringsAsFactors = FALSE)
head(population)

These variables are:

• LocID - The location numeric identifier. This variable will be dropped because the identifier is not useful for our purpose.
• Location - The location name. The location could be a country, region, or subregion (e.g. Central Asia), but this report will only be using countries due to the nature of the second data set.
• VarID - The variant numeric identifier. This variable will be dropped because the identifier is not useful for our purpose.
• Variant - Projection variant name. This report will only be using the medium projection variant because it is the most used (United Nations 2019).
• Time - The year from 1950 to 2100. This report will only be using 2016 because it is the most recent year available in the second data set.
• MidPeriod - The mid period of the data (e.g. 1950.5 indicates July 1950). This variable will be dropped because it is not useful for our purpose.
• PopMale - The male population in thousands.
• PopFemale - The female population in thousands.
• PopTotal - The male and female population in thousands. This variable will be dropped since it can be easily calculated from the male and female if needed.
• PopDensity - Population per square kilometre in thousands. This variable will be dropped because this report does not discuss the density.

# Subset the population data set.
populationFiltered <- population %>%
                      filter(Variant == "Medium") %>%
                      filter(Time == 2016) %>%
                      select(Location, PopMale, PopFemale)
head(populationFiltered)

The second data set is the World Health Organization suicide rate estimate (World Health Organization 2018). The data set contains information on the estimate of suicide rate on each country and gender. When reading the csv file, the first line is skipped because it contains the column header.

suicide <- read.csv("data.csv", stringsAsFactors = FALSE, skip=1)
head(suicide)

These variables are:

• Country - The country name.
• Sex - The gender.
• X2016, X2015, X2010, X2000 - The suicide rate estimate for that year per 100,000 population. Since the year available is not continuous, analysing the trend would not make so much sense. Hence, this report will only use the most recent year available, that is, 2016.

# Since we are only interested in the year 2016, remove columns of other years.
suicideFiltered <- suicide[ , !(names(suicide) %in% c("X2015", "X2010", "X2000"))]

# Change column X2016 to SuicideRate
names(suicideFiltered)[names(suicideFiltered) == "X2016"] <- "SuicideRate"

head(suicideFiltered)

Tidy & Manipulate Data I

The population data set is untidy because the variable gender does not have its own column and the variable population is spread over two columns.

So the wide data frame is converted into a long format with the gather() function by converting the “PopMale” and “PopFemale” columns into “Sex” and “Population”. Then the string “Pop” from “PopMale” and “PopFemale” is removed.

populationTidied <- populationFiltered %>%
                    gather(PopMale, PopFemale, key = "Sex", value = "Population") %>%
                    mutate(Sex = sub("Pop", "", Sex))
head(populationTidied)

Merge

As mentioned in the first section, the population data set contains locations other than countries while the suicide data set contains only countries. So the two data sets will be joined by countries and gender. But first, we will check which countries in the suicide data set that are not in the population data set.

anti_join(suicideFiltered, populationTidied, by = c("Country" = "Location")) %>% distinct(Country)

I looked at the population data set and it turns out these countries actually exist in the population data set under a different name. So let us manually change their names so that they match. Then we merge the two data sets by country and gender. Since all the countries in the suicide data set exists in the population data set, we can either use left join on suicide data set or inner join. We will use inner join.

# Rename countries that have different name on the two data sets.
populationTidied$Location[populationTidied$Location == "Dem. People's Republic of Korea"] <- "Democratic People's Republic of Korea"
populationTidied$Location[populationTidied$Location == "North Macedonia"] <- "Republic of North Macedonia"
populationTidied$Location[populationTidied$Location == "Micronesia (Fed. States of)"] <- "Micronesia (Federated States of)"
populationTidied$Location[populationTidied$Location == "United Kingdom"] <- "United Kingdom of Great Britain and Northern Ireland"

df <- suicideFiltered %>%
      inner_join(populationTidied, c("Country" = "Location", "Sex" = "Sex"))

Understand

The merged data set:

str(df)

'data.frame':   366 obs. of  4 variables:
 $ Country    : chr  "Afghanistan" "Afghanistan" "Albania" "Albania" ...
 $ Sex        : chr  "Male" "Female" "Male" "Female" ...
 $ SuicideRate: num  10.6 2.1 7 4.3 4.9 1.8 14 4.6 0 0.9 ...
 $ Population : num  18187 17196 1471 1416 20482 ...

As seen above, the gender variable has the character type, so we need to convert the gender variable into an unordered factor.

df$Sex <- factor(df$Sex, levels=c("Male","Female"), ordered = FALSE)

Now all the variables are in the correct type.

• Country - The country name.
• Sex - The gender. An unordered factor variable with 2 possible values: male or female.
• SuicideRate - The 2016 estimate suicide rate per 100,000 population.
• Population -The 2016 estimate population in thousands.

str(df)

'data.frame':   366 obs. of  4 variables:
 $ Country    : chr  "Afghanistan" "Afghanistan" "Albania" "Albania" ...
 $ Sex        : Factor w/ 2 levels "Male","Female": 1 2 1 2 1 2 1 2 1 2 ...
 $ SuicideRate: num  10.6 2.1 7 4.3 4.9 1.8 14 4.6 0 0.9 ...
 $ Population : num  18187 17196 1471 1416 20482 ...

Tidy & Manipulate Data II

A new variable for the estimate number of suicide in thousands for the year 2016 will be added. To do this, the population variable is multiplied with suicide rate and divided by 100,000 because the suicide rate is per 100,000 population. Then the result is rounded to 3 decimal places because the variable is in thousands.

df <- df %>% mutate(NumOfSuicide = round(Population * SuicideRate / 100000, 3))
head(df)

The new variable is called “NumOfSuicide”. For example, the estimate number of suicide for Afghanistan male in 2016 is 1,928 men. Please note that the number of suicide is not the actual number of suicide for the year 2016. The variables suicide rate and population are estimates from WHO and UN respectively. Hence, the resulting number of suicide is also an estimate.

Scan I

The is.na() function will be used to determine missing values. It returns TRUE if the value is NA. So we will use it with colSums() to count the number of NA on each column. Determine the number of missing values on each of the columns:

colSums(is.na(df))

     Country          Sex  SuicideRate   Population NumOfSuicide 
           0            0            0            0            0 

It seems that there is no missing values in the data frame.

Next, we will check whether the data is complete in terms of countries and gender. For example, there are no duplicates (e.g. multiple rows for Australia male) or missing rows. That is, the number of rows in the data frame should be equal to the number of unique countries multiplied by 2 (male and female). Which should also be equal to the number of unique country-gender combination.

# Number of unique countries.
df$Country %>% unique %>% length

[1] 183

# Number of unique gender.
df$Sex %>% unique %>% length

[1] 2

# Number of unique country-gender combination.
df %>% unite(CountrySex, Country, Sex) %>% .$CountrySex %>% unique %>% length

[1] 366

# Number of rows in the data frame.
nrow(df)

[1] 366

There are 183 different countries and 2 different levels for gender. 183*2 == 366. Looks like we have a complete data.

Scan II

We have three numeric variables: suicide rate, population, and number of suicide estimates. First, we will determine whether the distribution of these three variables are normal or not with box plots. We will create a long data frame by combining the three numeric variables into a single column so that we could easily plot a faceted box plot.

longDf <- df %>%
          gather(SuicideRate, Population, NumOfSuicide, key = "Var", value = "Value")

# Convert the Var column into a factor so that we can order the box plot.
longDf$Var <- factor(longDf$Var, levels=c("SuicideRate", "Population", "NumOfSuicide"))

# A better description of the variables for the box plots.
levels(longDf$Var) <- c("Suicide per 100,000", "Population in thousands", "Num of Suicide in thousands")

bp <- longDf %>%
      ggplot(aes(x=Var, y=Value)) + 
      geom_boxplot(aes(fill=Sex)) +
      facet_wrap( ~ Var, scales="free")
bp

Suicide rate has a right skewed distribution.

However, the box plots for population and number of suicide are not quite visible because there are few values that are substantially larger than the rest. So we will limit the y axis of the box plot by using the function coord_cartesian(). The function does not exclude data, it only limits the axis drawn. Thus, the median and quartiles will not be affected.

bpPop <- longDf %>%
         filter(Var=="Population in thousands") %>%
         ggplot(aes(x=Var, y=Value)) + 
         geom_boxplot(aes(fill=Sex)) +
         facet_wrap( ~ Var, scales="free") +
         coord_cartesian(ylim = c(0, 100000)) # Limit the y axis to 100,000
bpPop

bpNumOfSuicide <- longDf %>%
                  filter(Var=="Num of Suicide in thousands") %>%
                  ggplot(aes(x=Var, y=Value)) + 
                  geom_boxplot(aes(fill=Sex)) +
                  facet_wrap( ~ Var, scales="free") +
                  coord_cartesian(ylim = c(0, 10)) # Limit the y axis to 10
bpNumOfSuicide

It seems that population and number of suicide also have right skewed distribution. Using the z-score method to determine outliers will not be appropriate because it requires the values to have a normal distribution. So we will use the Tukey’s method.

According to the box plots above, all the outliers are located above the box plot, so we only need to calculate the maximum fence. Moreover, the box plots show that there are substantial difference between male and female suicide rate and number of suicide, so we will look for outliers by gender.

We need to calculate the maximum fence on each of the three numeric variables on each of the gender. The maximum fence is calculated by Q3 + 1.5 * IQR where Q3 is the 75^th percentile and IQR is the inter-quartile range.

summary <- longDf %>%
           group_by(Var, Sex) %>%
           summarise(Max = quantile(Value, 0.75) + 1.5 * IQR(Value))
summary

The “Max” column represents the maximum fence where values greater than this will be considered outliers.

To find the suicide rate outliers, we have to filter the data frame by gender and look for suicide rate that is higher than 33.45 (for male) and 11.60 (for female).

Suicide rate outliers (the outliers are sorted by gender and suicide rate):

df %>%
  filter((Sex == "Male" & SuicideRate > summary$Max[1]) | (Sex == "Female" & SuicideRate > summary$Max[2])) %>%
  arrange(Sex, desc(SuicideRate))

The same step is also applicable in finding the population and number of suicide outliers. But please note that different max fence is used (e.g. summary$Max[3])).

Population outliers (the outliers are sorted by gender and population):

df %>%
  filter((Sex == "Male" & Population > summary$Max[3]) | (Sex == "Female" & Population > summary$Max[4])) %>%
  arrange(Sex, desc(Population))

Number of suicide outliers (the outliers are sorted by gender and number of suicide):

df %>%
  filter((Sex == "Male" & NumOfSuicide > summary$Max[5]) | (Sex == "Female" & NumOfSuicide > summary$Max[6])) %>%
  arrange(Sex, desc(NumOfSuicide))

I think these outliers are not a result of data entry errors, so it would not be good idea to remove them. That is, I think the numbers on these outliers seem to make sense. For example, China and India are the two most populated countries in the world. I also did some quick googling on the number of suicides on several of these countries that are considered outliers and the numbers seem to be accurate. Therefore, I decided not to impute these outliers and leave them as they are.

Transform

In this section, we will apply a transformation on the suicide rate variable to convert the distribution into a more symmetric / normal distribution. As seen on the box plot in one of the earlier section, the suicide rate variable have a right skewed distribution. There are several ways to reduce the skewness of a right skewed distribution, such as log transformation, square root transformation, reciprocal transformation, and Box-Cox transformation. I tried each of these possible methods with different parameters and determined that the Box-Cox transformation seems to work slightly better than the rest for the suicide rate variable. So in this section, we will do a Box-Cox transformation on the suicide rate variable.

We have determined in the previous section that the suicide rate distribution is affected by gender. So we will apply the Box-Cox transformation by gender. First, we will filter the suicide rate variable by gender. Then apply the Box-Cox transformation. Then plot the histogram to visualise the distribution.

maleDf <- df %>% filter(Sex=="Male")
boxCoxMaleSuicideRate <- BoxCox(maleDf$SuicideRate, lambda="auto")
hist(boxCoxMaleSuicideRate)

It seems that there is a value (around -25000) that is very far from the rest of the group. Let us see find the value of this “outlier” by sorting the transformed values and look at the first few values.

boxCoxMaleSuicideRate %>% sort %>% head

[1] -2.437681e+04 -2.231425e-01  7.419486e-01  9.163080e-01  9.163080e-01  1.029641e+00

Let us find the lambda value used:

boxCoxMaleSuicideRate

  [1]  2.360968e+00  1.945988e+00  1.589287e+00  2.639200e+00 -2.437681e+04  2.708201e+00  2.312645e+00  2.856638e+00
  [9]  2.862369e+00  1.458659e+00  1.029641e+00  2.066950e+00  1.704808e+00 -2.231425e-01  3.671501e+00  3.100289e+00
 [17]  2.292643e+00  3.118149e+00  2.624810e+00  2.827478e+00  2.360968e+00  2.907074e+00  2.272232e+00  1.824618e+00
 [25]  2.572748e+00  3.109259e+00  3.140035e+00  3.182420e+00  2.197324e+00  3.292349e+00  2.714846e+00  2.890543e+00
 [33]  2.839244e+00  2.772746e+00  2.066950e+00  2.442469e+00  2.868068e+00  2.632031e+00  2.549578e+00  3.465982e+00
 [41]  2.934033e+00  2.797442e+00  1.974161e+00  2.845075e+00  2.694776e+00  2.708201e+00  2.580353e+00  2.476664e+00
 [49]  2.884971e+00  2.370359e+00  1.974161e+00  3.211055e+00  3.443861e+00  3.109259e+00  3.242808e+00  3.234964e+00
 [57]  2.928699e+00  2.174849e+00  3.035142e+00  2.884971e+00  2.708201e+00  2.549578e+00  2.509728e+00  2.610210e+00
 [65]  2.760166e+00  1.808356e+00  7.419486e-01  1.481650e+00  2.541734e+00  2.186149e+00  3.841903e+00  2.907074e+00
 [73]  1.667764e+00  3.100289e+00  3.077507e+00  2.917945e+00  1.648714e+00  1.589287e+00  1.547612e+00  2.868068e+00
 [81]  2.104225e+00  2.128325e+00  1.163179e+00  3.020612e+00  1.547612e+00  3.691656e+00  2.272232e+00  3.254460e+00
 [89]  9.163080e-01  2.694776e+00  2.557361e+00  3.434229e+00  1.435127e+00  3.122565e+00  2.624810e+00  2.163419e+00
 [97]  3.861035e+00  2.708201e+00  2.351489e+00  2.617536e+00  2.163419e+00  1.280968e+00  2.602829e+00  2.332255e+00
[105]  2.493333e+00  2.525859e+00  2.104225e+00  2.785170e+00  3.148657e+00  2.533828e+00  9.163080e-01  2.639200e+00
[113]  1.840619e+00  2.965453e+00  2.433735e+00  2.557361e+00  2.850873e+00  2.955089e+00  2.442469e+00  2.862369e+00
[121]  2.610210e+00  1.568666e+00  1.098637e+00  2.028233e+00  2.322498e+00  2.509728e+00  2.028233e+00  1.648714e+00
[129]  3.174085e+00  2.660405e+00  1.987955e+00  3.388010e+00  3.182420e+00  2.272232e+00  2.632031e+00  3.877740e+00
[137]  2.827478e+00  2.541734e+00  1.361015e+00  2.163419e+00  1.435127e+00  1.526104e+00  3.010807e+00  2.850873e+00
[145]  2.708201e+00  2.901594e+00  2.407064e+00  2.912525e+00  3.109259e+00  2.140160e+00  2.442469e+00  3.077507e+00
[153]  2.116347e+00  2.230116e+00  3.148657e+00  2.674295e+00  3.586557e+00  2.760166e+00  2.760166e+00  1.335038e+00
[161]  1.609491e+00  3.063583e+00  2.197324e+00  3.122565e+00  1.648714e+00  3.086682e+00  1.481650e+00  2.424923e+00
[169]  2.398013e+00  3.054192e+00  3.541217e+00  1.252795e+00  2.476664e+00  2.660405e+00  3.049464e+00  3.288624e+00
[177]  2.332255e+00  2.091954e+00  1.887143e+00  2.379662e+00  2.595393e+00  2.862369e+00  3.370971e+00
attr(,"lambda")
[1] 4.102259e-05

The lambda used is 4.102259e-05.

Then we can find the original row of this “outlier” by finding the index of the value in the transformed variable.

maleDf[which(boxCoxMaleSuicideRate == min(boxCoxMaleSuicideRate)), ]

The original suicide rate is 0. Which makes sense because the Box-Cox formula is \(\frac{y^\lambda-1}{\lambda}\) and if \(y\) is 0 and \(\lambda\) is 4.102259e-05, then the result is \(\frac{-1}{4.102259e-05} = -24376.81\)

We will remove this value and draw the histogram again so that we can see a better picture of the distribution of the rest of the group.

hist(boxCoxMaleSuicideRate[boxCoxMaleSuicideRate != min(boxCoxMaleSuicideRate)], breaks=20)

The distribution is not exactly symmetric but is definitely more similar to a normal distribution than the original values.

Now we will apply Box-Cox transformation on the female suicide rate.

femaleDf <- df %>% filter(Sex=="Female")
boxCoxFemaleSuicideRate <- BoxCox(femaleDf$SuicideRate, lambda="auto")
hist(boxCoxFemaleSuicideRate)

Again, the distribution is not exactly symmetric but is definitely more similar to a normal distribution than the original values.

References

• United Nations 2019, World Population Prospects, United Nations, viewed 1 June 2020, <https://population.un.org/wpp/Download/Standard/CSV/>.
• World Health Organization 2018, Suicide rate estimates, age-standardized, World Health Organization, viewed 1 June 2020, <https://apps.who.int/gho/data/node.main.MHSUICIDEASDR?lang=en>.