---
title: 'Data Visualisation Challenge: 3D scatter plots with **plotly**'
author: "Lewis Rendell"
date: "13th March 2019"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## Preparing the data

We load the **rvest** package, for `scraping' (extracting) data from web pages, and the **dplyr** package, for reshaping and transforming the data.

```{r dataPackages, message=FALSE, warning=FALSE}
library(rvest)
library(dplyr)
```

We scrape the data from [the page](https://commons.wikimedia.org/wiki/File:Development_of_a_rational_scale_to_assess_the_harm_of_drugs_of_potential_misuse_%28physical_harm_and_dependence,_NA_free_means%29.svg) on the Wikimedia Commons website where it is displayed.
To do this, we require the URL of the web page, and the 'XPath' of the element we wish to scrape (i.e. the table of data).
A guide of how to obtain this, and how to use **rvest** to conduct the scraping, is provided in [this RBloggers post by Cory Nissen](https://www.r-bloggers.com/using-rvest-to-scrape-an-html-table/).

```{r dataScraping}
# URL of web page containing table of data:
url <- "https://commons.wikimedia.org/wiki/File:Development_of_a_rational_scale_to_assess_the_harm_of_drugs_of_potential_misuse_(physical_harm_and_dependence,_NA_free_means).svg"
# XPath of table of data:
xpath <- '//*[@id="mw-imagepage-content"]/div/table[3]'

# Scrape this table and store as a data frame.
dataRaw <- url %>%
  read_html() %>%                 # Read the web page with the given URL.
  html_node(xpath = xpath) %>%    # Extract the element with the given XPath.
  html_table(header = TRUE)       # Convert to a data frame (with the first row used as the header).
```

In the original data frame, each drug is scored on a scale from 0 to 3 for nine parameters, which are further grouped into three categories. For each of these three categories, the mean score is provided.
The three axes in our 3D scatter plot will correspond to the mean scores in each of these three categories.
These are:

- Mean physical harm
- Mean dependence
- Mean social harm

For the first of these, we take the mean value in which `NA` values are interpreted as zero (some drugs have `NA` scores for 'intravenous harm', as they cannot be injected intravenously).

To make it easier to plot the data, we will select only these three columns (together with the drug names), and rename them so they are easier to reference.

```{r dataWrangling}
# Create a named character vector of the variables to keep, where the names of this vector are the new names to be used.
variableNames <- names(dataRaw)
variablesToKeep <- c(drug = variableNames[1],
                     physicalHarm = variableNames[5],
                     dependence = variableNames[9],
                     socialHarm = variableNames[13])

# Create a new data frame containing only these variables.
dataProcessed <- dataRaw %>% 
  select(variablesToKeep)
```

This gives us the following data frame:

```{r dataFrame, echo=FALSE, warning=FALSE, results='asis'}
library(knitr)
kable(dataProcessed)
```

## Producing a 3D scatter plot

The **graphics** package, which comes included with R when installed, includes the `persp` function for creating 3D surface plots.
For 3D scatter plots, functionality is provided by the **3dscatterplot** and **rgl** package (see [examples here](https://stackoverflow.com/questions/40359343/how-to-use-persp-with-r-to-graph-three-variables-in-a-dataset)).

For our 3D scatter plot we will use the **plotly** package.
This produces an interactive plot that can be rotated and zoomed.

```{r plottingPackages, message=FALSE, warning=FALSE}
library(plotly)
```

[The Plotly web page on 3D scatter plots](https://plot.ly/r/3d-scatter-plots/) gives several examples of the syntax used to generate 3D scatter plots. Note that this syntax differs from that used by the **ggplot2** package.

```{r plot1, message=FALSE, warning=FALSE}
plot_ly(dataProcessed,                                          # Plot data frame.
        x = ~physicalHarm,                                      # Choose variables for x, y, z axes.
        y = ~dependence,
        z = ~socialHarm,
        text = ~drug) %>%                                       # Choose variable for hover text.
  add_markers(hoverinfo = "text") %>%                           # Add points with hover text.
  layout(scene = list(xaxis = list(title = "Physical harm",     # Set axis labels and ranges.
                                   range = c(0, 3)),
                      yaxis = list(title = "Dependence",
                                   range = c(0, 3)),
                      zaxis = list(title = "Social harm",
                                   range = c(0, 3))))
```

We may wish to visualise the relationship between these three variables and the legal status of a drug, to see whether drugs considered more dangerous in law are indeed more dangerous by these measures.
For this purpose, a CSV file has been compiled, containing the legal status of each drug on the list.

All the illegal drugs on the list fall into one of three classes according to the [Misuse of Drugs Act 1971](https://en.wikipedia.org/wiki/Drugs_controlled_by_the_UK_Misuse_of_Drugs_Act): 'Class A', 'Class B' and 'Class C'.
Of the remaining drugs, some can be sold legally, but not for human consumption ([solvents under the Psychoactive Substances Act 2016](https://www.gov.uk/government/publications/psychoactive-substances-act-guidance-for-retailers/psychoactive-substances-act-2016-guidance-for-retailers), and [alkyl nitrites under the Medicines Act 1968](https://www.talktofrank.com/drug/poppers#the-law)) -- these are recorded as 'Legal with restrictions'.
Alcohol and tobacco can be sold for human consumption (subject to age restrictions) and are therefore recorded as 'Legal'.

We read this CSV file and join the data to our existing data frame.

```{r legalityData}
# Read CSV file containing legal status of each drug, storing as a data frame.
drugLegality <- read.csv("drugClasses.csv", stringsAsFactors = FALSE)
# Save the column storing these legal statuses as an ordered factor.
legalityLevels <- c("Class A", "Class B", "Class C", "Legal with restrictions", "Legal")
drugLegality[["legalStatus"]] <- factor(drugLegality[["legalStatus"]],
                                        levels = legalityLevels)
# Join this data frame to the original data frame.
dataWithLegality <- left_join(dataProcessed, drugLegality, by = "drug")
```

This gives us the following data frame:

```{r dataFrameWithLegalStatuses, echo=FALSE, warning=FALSE, results='asis'}
kable(dataWithLegality)
```

Colouring the points according to their legal status, we obtain the following plot:

```{r plot2}
# Define a colour palette: reds/oranges for illegal drugs, blues for legal drugs.
palette <- c("darkred", "red", "orange", "deepskyblue", "blue")

plot_ly(dataWithLegality,
        x = ~physicalHarm,
        y = ~dependence,
        z = ~socialHarm,
        text = ~drug,
        color = ~legalStatus,
        colors = palette) %>% 
  add_markers(hoverinfo = "text") %>%
  layout(scene = list(xaxis = list(title = "Physical harm",
                                   range = c(0, 3)),
                      yaxis = list(title = "Dependence",
                                   range = c(0, 3)),
                      zaxis = list(title = "Social harm",
                                   range = c(0, 3))))
```


