MATH2349 Semester 1, 2020

Required packages

library(readxl)
library(tidyr)
library(dplyr)
library(magrittr)
library(stringr)
library(editrules)

Executive Summary

The aim of this assignment is to pre-process crime data and gambling expenditure data in Victoria to prepare it for future analysis on the association between crime and gambling expenditure.

The two data sets, recorded_offences.xlsx and gambling_expenditure.xls were imported and read. Both data sets have been done string manipulations for strings before merging with a common variable (Local Government Area). To understand the merged data set, inspection of its structure and its data types were conducted. Then, some nominal variables were factorised while no ordinal variables in the data set.

The data was untidy, hence for reshaping the data into a tidy format, column names were renamed, and irrelevant variable was simply dropped. Two new columns were created from the existing ones. The new ranking variables were sorted by year.

Several missing values were identified in the data set. The main cause of missing values was figured out through examining the data with NA values. It was originally recorded as unknown geographical location, so we concluded that it is reasonable to delete from the data set. Accordingly, special values and obvious errors (i.e. inconsistencies) were checked but no its values and errors were not found.

Finally, removal of outliers and transformation were performed to try to reduce the effects of outliers. The capping method was used to handle the outliers. We capped outlier values with the closest 5th percentile. On the heavily right skewed, data, a logarithmic transformation via the natural logarithm was applied to the variable to reduce the skewness before capping. This resulted in a more normal distribution and eliminated much of the perceived outliers.

Data

Crime data set

The main dataset for this assignment is about latest crime taken from Crime Statistics Agency https://www.crimestatistics.vic.gov.au/crime-statistics/latest-crime-data/download-data-0 containing 870 observations from Local government areas about recorded offences rate per 100,000 population by police region from 2010 to 2019. The variables for recorded_offences.xlsx data set consists of:

Variable	Description
Year	The year that data was collected. (quantitative variable)
Police Region	Police regions (levles will be changed to: North West Metro, Eastern, Southern Metro, Western and where Total2 is for unknown geographical location) (qualitative variable)
Local Government Area	79 Local government areas in Victoria and there are totals per police region. (qualitative variable)
Offence Count	Total number of offences (quantitative variable)
Rate per 100,000 population	Number of offences that been conducted in the area per 100,000 population (quantitative variable)

The setwd() function sets the working directory to the folder “Assignment 2”.
readxl package is used to read in the data.

setwd("C:/Users/lisal/Desktop/Data Wrangling/Assignment 2")
offences <- read_excel("recorded_offences.xlsx", sheet = 2)

head() is used to show the first 6 observation of the data. I have subsetted the columns so that each column can be observed.

head(offences)

The variables are renamed by names().

names(offences)[3] <- "LGA"
names(offences)[4] <- "Count"
names(offences)[5] <- "Rate"

Year was filtered after 2011 using filter() due to data unavailability of the year, 2010 in the second data set.

offences2 <- offences %>%  filter( Year > 2010 )

LGA observations are edited by recode() to apply with the same data criteria as another data set, gambling which is explained next.

offences3 <- offences2 %>% mutate(
  LGA = dplyr::recode(LGA, "Nillumbik"="Whittlesea"),
  LGA = recode(LGA, "Ararat"="Northern Grampians"),
  LGA = recode(LGA, "Queenscliffe"="Greater Geelong"),
  LGA = recode(LGA, "Corangamite"="Colac-Otway"),
  LGA = recode(LGA,"Hepburn"="Moorabool"),
  LGA = recode(LGA, "Mount Alexander"="Central Goldfields"),
  LGA = recode(LGA, "Mansfield"="Mitchell"),
  LGA = recode(LGA, "Murrindindi"="Mitchell"),
  LGA = recode(LGA, "Towong"="Alpine"),
  LGA = recode(LGA, "Moira"="Benalla"), 
  LGA = recode(LGA, "Strathbogie"="Benalla"),
  LGA = recode(LGA, "Gannawarra"="Campaspe"),
  LGA = recode(LGA, "Southern Grampians"="Glenelg"))

Gambling Expenditure data set

The second dataset is scraped from the Victorian Commission for Gambling and Liquor Regulation https://www.vcglr.vic.gov.au/resources/data-and-research/gambling-data/gaming-expenditure-local-area. Its variables contain information on total gaming expenditure by local government area. The gambling_expenditure.xls data set consists of below mainly 3 variables:

Variable	Description
LGA Name	57 Local government area in Victoria (qualitative variable)
Region	Gaming venues are classifie one of two regions, C - country, M - metro (qualitative variable, nominal)
Expenditure	Amount of money lost by gaming patrons (quantitative variable)

Note that we can see that the common variable of the two datasets can be Local Government Area (LGA). Having referred to the data clarification of gambling, some local government areas were amalgamated in this data set and stored as a single LGA. For example, CITY OF WHITTLESEA consists of City of Whittlesea and Shire of Nillumbik. Accordingly, the data manipulation in crime has undertaken earlier.
The following is a list of amalgamated LGA’s in gambling:

Local Government Area	Consists of
CITY OF WHITTLESEA	City of Whittlesea and Shire of Nillumbik
SHIRE OF NORTHERN GRAMPIANS	Rural City of Ararat and Shire of Northern Grampians
CITY OF GREATER GEELONG	Borough of Queenscliffe and City of Greater Geelong
SHIRE OF COLAC-OTWAY	Shire of Corangamite and Shire of Colac-Otway
SHIRE OF MOORABOOL	Shire of Hepburn and Shire of Moorabool
SHIRE OF CENTRAL GOLDFIELDS	Shire of Central Goldfields and Shire of Mount Alexander
SHIRE OF MITCHELL	Shire of Mansfield, Shire of Murrindindi and Shire of Mitchell
SHIRE OF ALPINE	Shire of Towong and Shire of Alpine
RURAL CITY OF BENALLA	Shire of Moira, Shire of Strathbogie, and Rural City of Benalla
SHIRE OF CAMPASPE	Shire of Gannawarra and Shire of Campaspe
SHIRE OF GLENELG	Shire of Glenelg and Shire of Southern Grampians

The data discrepancy is also seen by comparing to the main data set, crime, which will be dealt with in the Scan section.
Import and read using the readxl package.

setwd("C:/Users/lisal/Desktop/Data Wrangling/Assignment 2")
gambling <- read_excel("gambling_expenditure.xls", sheet = 3, range = cell_rows(9:66))

From below, it is seen that some LGA Name observations might be better to be edited such as removing the strings “City of”, “Shire of” or “Rural City of” in front of LGA names. This is so that the LGA Name will be able to match with the LGA in the crime so that the two datasets may be joined in the next step:

head(gambling[,1:3])

head(gambling[,4:6])

head(gambling[,7:9])

head(gambling[,10:11])

Pattern matching functions to detect, subset and replace strings under stringr package, str_detect() and str_replace() are used.
mutate() and ifelse() from base R allow users to conduct a logical test across a single variable (or vector), and then populate the fields of a new variable depending on the outcome of the tests.

gambling2 <- gambling %>% 
  mutate(
    `LGA Name` =  
      ifelse(str_detect(`LGA Name`, regex("Shire of", ignore_case = TRUE)), 
             str_replace(`LGA Name`, pattern= regex("Shire of ", ignore_case = TRUE), replacement=""),
      ifelse(str_detect(`LGA Name`, regex("Rural city of ", ignore_case = TRUE)), 
             str_replace(`LGA Name`, pattern= regex("Rural city of ", ignore_case = TRUE), replacement=""),
      ifelse(str_detect(`LGA Name`, regex("City of", ignore_case = TRUE)), 
             str_replace(`LGA Name`, pattern= regex("City of ", ignore_case = TRUE), replacement=""),
                                        "Error"))))

In adition, those almagated oberservations were originally recoded in capital letters so that they are converted to the same format as others.

gambling3 <- gambling2 %>% mutate(
  `LGA Name` = dplyr::recode(`LGA Name`, "WHITTLESEA"="Whittlesea"),
  `LGA Name` = recode(`LGA Name`, "NORTHERN GRAMPIANS"="Northern Grampians"),
  `LGA Name` = recode(`LGA Name`, "GREATER GEELONG"="Greater Geelong"),
  `LGA Name` = recode(`LGA Name`, "COLAC-OTWAY"="Colac-Otway"),
  `LGA Name` = recode(`LGA Name`, "MOORABOOL"="Moorabool"),
  `LGA Name` = recode(`LGA Name`, "CENTRAL GOLDFIELDS"="Central Goldfields"),
  `LGA Name` = recode(`LGA Name`, "MITCHELL"="Mitchell"),
  `LGA Name` = recode(`LGA Name`, "ALPINE"="Alpine"),
  `LGA Name` = recode(`LGA Name`, "BENALLA"="Benalla"), 
  `LGA Name` = recode(`LGA Name`, "CAMPASPE"="Campaspe"),
  `LGA Name` = recode(`LGA Name`, "GLENELG"="Glenelg"))

To combine offences3 and gambling3 into a single tibble, we need to use dplyr::left_join(). The two datasets are joined as following by the common key, LGA, to form a dataset with 15 variables: (requirement #1)

data <- left_join(offences3, gambling3, by = c("LGA" = "LGA Name"))
head(data)

Understand

In this section, we will show that we have fulfilled the requirements #2-4.
The dimensions are checked using dim() function after combined.

dim(data)

## [1] 783  15

Data types are checked using str(). From the following output, the data set includse multiple data types (requirement #2).

str(data)

## tibble [783 x 15] (S3: tbl_df/tbl/data.frame)
##  $ Year                             : num [1:783] 2019 2019 2019 2019 2019 ...
##  $ Police Region                    : chr [1:783] "1 North West Metro" "1 North West Metro" "1 North West Metro" "1 North West Metro" ...
##  $ LGA                              : chr [1:783] "Banyule" "Brimbank" "Darebin" "Hobsons Bay" ...
##  $ Count                            : num [1:783] 9524 19784 14532 6224 21445 ...
##  $ Rate                             : num [1:783] 7247 9410 8855 6377 9204 ...
##  $ Region                           : chr [1:783] "M" "M" "M" "M" ...
##  $ Expenditure
## 1 Jul 10 - 30 Jun 11: num [1:783] 5.84e+07 1.39e+08 8.90e+07 5.14e+07 1.03e+08 ...
##  $ Expenditure
## 1 Jul 11 - 30 Jun 12: num [1:783] 5.78e+07 1.46e+08 8.93e+07 5.26e+07 1.04e+08 ...
##  $ Expenditure
## 1 Jul 12 - 30 Jun 13: num [1:783] 5.53e+07 1.38e+08 8.24e+07 4.79e+07 9.88e+07 ...
##  $ Expenditure
## 1 Jul 13 - 30 Jun 14: num [1:783] 5.50e+07 1.39e+08 8.24e+07 4.72e+07 1.02e+08 ...
##  $ Expenditure
## 1 Jul 14 - 30 Jun 15: num [1:783] 5.45e+07 1.42e+08 8.39e+07 4.64e+07 1.05e+08 ...
##  $ Expenditure
## 1 Jul 15 - 30 Jun 16: num [1:783] 5.60e+07 1.43e+08 8.43e+07 4.68e+07 1.06e+08 ...
##  $ Expenditure
## 1 Jul 16 - 30 Jun 17: num [1:783] 5.58e+07 1.34e+08 8.11e+07 4.69e+07 1.06e+08 ...
##  $ Expenditure
## 1 Jul 17 - 30 Jun 18: num [1:783] 5.85e+07 1.40e+08 8.21e+07 4.74e+07 1.10e+08 ...
##  $ Expenditure
## 1 Jul 18 - 30 Jun 19: num [1:783] 5.78e+07 1.43e+08 8.16e+07 4.70e+07 1.12e+08 ...

After checking the structure of the data types, we can see nominal variable is not in the correct format which needs to be factorised. There are no ordinal variables in the data set.
Police Region and Region are read in as character formats, we factorised these variables using factor() and defined proper levels using levels() argument. (requirement #3, #4)

data$`Police Region` <- 
  factor(data$`Police Region`,
             levels = c("1 North West Metro", "2 Eastern", "3 Southern Metro", "4 Western", 
                        "Justice Institutions and Immigration Facilities", "Unincorporated Vic" ),
             labels = c("North West Metro", "Eastern", "Southern Metro", "Western", "Total2", "Total2"))
levels(data$`Police Region`)

## [1] "North West Metro" "Eastern"          "Southern Metro"   "Western"         
## [5] "Total2"

data$Region <- factor(data$Region,
                          levels = c("C", "M"),
                          labels = c("country", "metro"))
levels(data$Region)

## [1] "country" "metro"

Now all variables seem to have applied proper data type conversion.

Tidy & Manipulate Data I

The second data set, gambling does not conform the tidy data principals and is untidy (requirement #5).
It has set of columns, Expenditure Jul 10 - 30 Jun 11, Expenditure Jul 11 - 30 Jun 12…etc, which are broken down by Expenditure and Year. In short, these columns’ headers represent values, not variable, and multiple variables stored in one column. Therefore, the gambling dataset violates the tidy data principles and is untidy.
To tidy the dataset, firstly, those columns are renamed to each year.

names(data)[7:15] <- c(2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019)

We need to gather those columns into a new pair of variables by gather().
The name of the variable whose values form the column names. In this case, the key argument is year2.
The name of variable whose values are spread over the cells. The value argument, in this case, is the amount of expenditure.
The reason why we created another year variable, year2 is becuase the cartesian product of two sets of year2 and Expenditure was made when gathered. In fact, 7047 observations were generated in the dataset after gathering.

data2 <- data %>% gather(7:15, key = "year2", value = "Expenditure", convert = TRUE)
head(data2)

To eliminate the cartesion product, we need to use filter() to find where Year variables are matched with year2, and then year2 was removed. The final tidy data is shown using head():

data3 <- data2 %>% filter(Year==year2)
data3 <- data3[-c(7)]
head(data3)

This form is tidy. There is one variable in each column, and each row represents a year’s observations by LGA.

Tidy & Manipulate Data II

The mutate_at function is used to add new variables to an existing data, and at affects variables selected with vars().
Ranking variables of Rate, Expenditure split by year were created (requirement #6).
By default, min_rank() assigns 1 to the smallest value and high number to the largest value. In case, we assigned rank 1 to the largest value of a variable, used min_rank(desc(.))

data4 <- data3 %>% 
  group_by(Year) %>%
  mutate_at(vars(`Rate`, Expenditure), funs(Rank=min_rank(desc(.))))
head(data4)[8:9]

dim(data4)

## [1] 783   9

We can see that new variables are added.

Scan I

First, missing values, special values and obvious errors such as inconsistencies are checked, and then these values are handled with proper reasonings (requirement #7).
sum() and is.na() functions were used as a function inside sapply() for calculating the total missing values for each column:

sapply(data4, function(x) sum(is.na(x)))

##             Year    Police Region              LGA            Count 
##                0                0                0                0 
##             Rate           Region      Expenditure        Rate_Rank 
##               36              153              153               36 
## Expenditure_Rank 
##              153

36 observations in the Rate variable are missing values. The next output shows that the details of the first 10 observations of Rate missing variable. Having examined, it is figured out that Police Region is “Total2”, which is unknown geographical location from the data definition.
Hence, Rate variable (originally Rate per 100,000 population) as well as Region are not available. These observations are simply excluded from the dataset.

data4 %>% filter(is.na(Rate))

153 observations of Region and Expenditure are counted as missing values. As noted in the second dataset definition, there is the data discrepancy / unavailability in LGA variable. Therefore, we concluded that it is a major cause of NA values in these variables, and it is reasonable to delete these observation from the dataset.
To handle missing values, na.omit() is used simply to omit the records which have NAs from the analysis.

data5 <- na.omit(data4)
dim(data5)

## [1] 630   9

In the same way, special values were checked whether every numerical column have infinite or NaN values using a function called is.special inside sapply() function which allows to apply to the entire dataset.

is.special <- function(x){ if (is.numeric(x)) (is.infinite(x) | is.nan(x)) } 
sapply(data4, function(x) sum(is.special(x)))

##             Year    Police Region              LGA            Count 
##                0                0                0                0 
##             Rate           Region      Expenditure        Rate_Rank 
##                0                0                0                0 
## Expenditure_Rank 
##                0

Using the editrules package, the rules defined below are to make sure that numerical value must be non negative.
The rule set can be read using editfile function to scan for inconsistencies as follows:

Rules <- editset(expression(
  Year > 0,
  Count > 0,
  Rate >= 0,
  Expenditure >= 0,
  Rate_Rank > 0,
  Expenditure_Rank > 0
))

Rules

## 
## Edit set:
## num1 : 0 < Year
## num2 : 0 < Count
## num3 : 0 <= Rate
## num4 : 0 <= Expenditure
## num5 : 0 < Rate_Rank
## num6 : 0 < Expenditure_Rank

summary(violatedEdits(Rules,data5))

## No violations detected, 0 checks evaluated to NA

## NULL

Scan II

In this section, we will scan all numeric variables (Count, Rate, Expenditure variables) except for the Year variable (requirements #8).
For choosing the appropriate method of outliers detection, we will first check the normality of the distribution using hist().

par(mfrow=c(1,3))

hist(data5$Count, main = "Number of Offences", xlab = "Number")
hist(data5$Rate, main = "Rate per 100,000 population", xlab = "Rate")
hist(data5$Expenditure, main = "Gambling Expenditure", xlab = "Expenditure")

All 3 distributions are right-skewed so that we will use the Tukey’s method of outlier detection as it does not require the assumption of normality.
We can use boxplot() function (under Base graphics) to get the box plot of all variables:

par(mfrow=c(1,3))

data5$Count %>% boxplot(main = "Number of Offences", ylab="Number", col = "lightgrey")
data5$Rate %>% boxplot(main = "Rate per 100,000 population", ylab="Rate", col = "lightgrey")
data5$Expenditure %>% boxplot(main = "Gambling Expenditure", ylab="Expenditure", col = "lightgrey")

According to the Tukey’s method, the Count variable and the Rate variable seem to have many outliers, while any outliers are not spotted in the Expenditure variable.
To understand more outliers, we put a label on the Count variable’s outliers. The following code finds all local government areas that had outliers in the Count variable:

Out_Count <- boxplot.stats(data5$Count)$out

data5 %>% filter(Count %in% Out_Count)

Interestingly, the local government areas which had one of the highest offence counts had also one of the highest Rate_Rank values as well as relativily high Expenditure_Rank values.
Besides, we cannot see any univariate outliers in Expenditure variable, thus we will examine multivariate outliers next.
The plot() function will be used to get the scatter plot and detect multivariate outliers in all 3 variables in order:

data5 %>% plot(Count ~ Rate, data = ., ylab="Count", xlab="Rate", main="Count by Rate")

data5 %>% plot(Expenditure ~ Count, data = ., ylab="Expenditure", xlab="Count", main="Expenditure by Count")

data5 %>% plot(Expenditure ~ Rate, data = ., ylab="Expenditure", xlab="Rate", main="Expenditure by Rate")

According to the scatter plots, there are some possible outliers on the higher right (the data points where both Count and Rate are big) and right hand side of the scatter (the data points with larger Count and middle high Expenditure values) and right hand side of the scatter (the data points with larger Rate and middle high Expenditure values).
For this dataset, we considered that transformation may be more suitable to handle outliers since the distributions are extremely right-skewed. However, as following to the minimum requirement #8, capping will be undertaken in this section by replacing the outliers with the nearest neighboughrs that are not outliers (later in the Transform section another technique will be explored).
In order to cap the outliers we can use a user-defined function as follows (taken from: Stackoverflow https://stackoverflow.com/questions/13339685/how-to-replace-outliers-with-the-5th-and-95th-percentile-values-in-r?utm_medium=organic&utm_source=google_rich_qa&utm_campaign=google_rich_qa):

cap <- function(x){ 
  quantiles <- quantile( x, c(.05, 0.25, 0.75, .95 ) )
  x[ x < quantiles[2] - 1.5*IQR(x) ] <- quantiles[1]
  x[ x > quantiles[3] + 1.5*IQR(x) ] <- quantiles[4]     
  x 
  }

We can apply this function to a subset of the dataset using sapply function.
Check descriptive statistics.

data5_sub <- data5 %>% select(Year, Count, Rate, Expenditure)
data5_capped <- sapply(data5_sub, FUN = cap)
summary(data5_sub)

##       Year          Count            Rate        Expenditure       
##  Min.   :2011   Min.   :   80   Min.   : 2665   Min.   :  2304602  
##  1st Qu.:2013   1st Qu.: 1485   1st Qu.: 5256   1st Qu.: 10575644  
##  Median :2015   Median : 4728   Median : 7269   Median : 23605461  
##  Mean   :2015   Mean   : 6758   Mean   : 7795   Mean   : 42376482  
##  3rd Qu.:2017   3rd Qu.:10311   3rd Qu.: 9345   3rd Qu.: 67532057  
##  Max.   :2019   Max.   :37858   Max.   :32309   Max.   :145619090

Check summary statistics again.

summary(data5_capped)

##       Year          Count            Rate        Expenditure       
##  Min.   :2011   Min.   :   80   Min.   : 2665   Min.   :  2304602  
##  1st Qu.:2013   1st Qu.: 1485   1st Qu.: 5256   1st Qu.: 10575644  
##  Median :2015   Median : 4728   Median : 7269   Median : 23605461  
##  Mean   :2015   Mean   : 6515   Mean   : 7567   Mean   : 42376482  
##  3rd Qu.:2017   3rd Qu.:10311   3rd Qu.: 9345   3rd Qu.: 67532057  
##  Max.   :2019   Max.   :22963   Max.   :15375   Max.   :145619090

We can see that the Count and Rate variable’s values are capped.

Transform

In this section, the Rate variable is transformed (requirement #9). The purpose of this transformation is to decrease the skewness and convert the distribution into a normal distribution.
Having applied different transformations on the same data and selected the Mathematical Operations which transformed the distribution into much more normal.
To illustrate the logarithmic transformation via natural logarithm, let’s revisit the Rate data.

hist(data5$Rate, main = "Rate - Original", xlab = "Rate")

From the histogram, we observe that Rate has a right-skewed distribution. By applying a logarithmic transformation, the Rate distribution would be more symmetrical.
The logarithmic transformation using the log() function as follows:

ln_Rate <- log(data5$Rate) 
hist(ln_Rate, main = "Rate - Transformed")

As seen from the histograms, the ln transformation reduced the skewness and it improved the symmetry of the Rate distribution.

References:

[1] Recorded Offences data. Source: https://www.crimestatistics.vic.gov.au/crime-statistics/latest-crime-data/download-data-0. Last access: 29 June 2020.

[2] Gambling Expenditure data. Source: https://www.vcglr.vic.gov.au/resources/data-and-research/gambling-data/gaming-expenditure-local-area Last access: 29 June 2020.

[3] Stackoverflow, Cap function. Source: https://stackoverflow.com/questions/13339685/how-to-replace-outliers-with-the-5th-and-95th-percentile-values-in-r?utm_medium=organic&utm_source=google_rich_qa&utm_campaign=google_rich_qa Last access: 29 June 2020.

[4] Mutate multiple columns, Reference. https://dplyr.tidyverse.org/reference/mutate_all.html Last access: 29 June 2020.

[5] MATH2349 Data Wrangling course website, Module Notes. http://rare-phoenix-161610.appspot.com/secured/index.html Last access: 29 June 2020.