### Scraping web data in table format

In addition to the `tidyverse` collection of packages, we will also load the `rvest` package which we will allow to scrape data from the web.

```{r}
library(rvest)
library(tidyverse)
```

The goal of this exercise is to scrape the counts of social security number applicants by year in the US, and then clean the data so that we can generate a plot showing the evolution in this variable over time.

The first step is to read the html code from the website we want to scrape using the `read_html()` function. If we want to see the html in text format, we can then use `html_text()`.

```{r}
url <- "https://www.ssa.gov/oact/babynames/numberUSbirths.html"
# Storing the url's HTML code
html_content <- read_html(url)
```


```{r}
# Not very informative
str(html_content)
print(html_content)
# Looking at first 1000 characters
substr(html_content, 1, 1000)
```

To extract all the tables in the html code automatically, we use `html_table()`. Note that it returns a list of data frames which has length 1 here as there is only one table on this website.

```{r}

# Extracting tables in the document
tab <- html_table(html_content, fill = TRUE)

# Check object
str(tab)
# Website only had one table -> list of length 1 containing the dataframe

# Save the dataframe (tibble)
social_security_data <- tab[[1]]
social_security_data <- as_tibble(social_security_data)
social_security_data
```

Now let us clean the data so that we can use it for our analysis. We need to convert the population values into a numeric format, which requires deleting the commas. We will also change the variable names so that it is easier to work with them.

```{r}
## Remove commas and then make variable numeric
social_security_data$Male <-  as.numeric(gsub(",", "", social_security_data$Male))
social_security_data$Female <- as.numeric(gsub(",", "", social_security_data$Female))
social_security_data$Total <- as.numeric(gsub(",", "", social_security_data$Total))

## Rename variables
colnames(social_security_data) <- c('year', 'male', 'female', 'total')
```

Now we can plot to see how the number of people applying for a Social Security Number in the US has increased over time (using `ggplot2` requires the data to be in a tidy long format):

```{r}
# Pivot longer function
social_security_data_long <- pivot_longer(data = social_security_data[, 1:3],
  cols = c("male", "female"), names_to = "sex", values_to = "individuals")

# Creating the plot
ggplot(social_security_data_long) +
  aes(x = year, y = individuals, group = sex, colour = sex) +
  geom_line() + xlab("Year") + ylab("Individuals") +
  scale_colour_manual(name ="Birth Sex", values = c(1, 2), labels = c("Female", "Male"))
```

### Scraping web data in table format: A more advanced example

When there are multiple tables on the website, scraping them becomes a bit more complicated. Let's work through an exemplary scenario: Scraping a table from Wikipedia with a list of the most populated cities in the United States.

```{r}
url <- "https://en.wikipedia.org/wiki/List_of_United_States_cities_by_population"
html <- read_html(url)
tables <- html_table(html, fill=TRUE)
class(tables)
length(tables)
```

The function returns 15 different tables. The option `fill=TRUE` is used because some of the tables appear to have incomplete rows.

One option is to go through these tables in this list manually and select the correct one:

```{r}
tables[5] # position 5 happens to be the correct table (checked 23/10/2023)
```

Alternatively, the table of interest can be identified directly with a right-click and "Inspect" or "Inspect Element" in the browser. Clicking on the relevant part of the page's code then allows to copy identifiers such as the CSS selector.

CSS selector: table.wikitable:nth-child(25) (extracted with Firefox "Inspect Element" on 23/10/2023)

An alternative option is to look at the full source code of the website. In Google Chrome e.g., go to _View_ > _Developer_ > _View Source_. All browsers should have similar options to view the source code of a website. In the source code, search for the text of the page (e.g. _2022 rank_). Right above it you will see: `<table class="wikitable sortable ..." ...>`. This is the CSS selector. Using this selector, however, might still return several tables.

Now that we know what we're looking for, let's use `html_elements()` to identify all the elements of the page that have that CSS class (note that we use a dot before the name of the class because the R function expects CSS notation). Whereas we will find several elements with the CSS selector ".wikitable" which we would have to search subsequently, we will only find one with the very specific selector obtained with inspect element.

```{r}
table_raw <- html_elements(html, css = "tba")
length(table_raw) # now a list of length 1
```

Transforming into a tibble with `html_table()`:

```{r}
data_pop <- html_table(table_raw[[1]], fill=TRUE)
data_pop
```

As in the previous case, we still need to clean the data before we can use it. For this particular example, let's see if this dataset provides evidence in support of [Zipf's law for population ranks](https://en.wikipedia.org/wiki/Zipf%27s_law). Keeping only the columns of interest and transforming into a tibble:

```{r}
data_pop <- data_pop %>% select(c("City", "2022rank", "2022estimate"))
data_pop
```

Renaming and cleaning columns:

```{r}
# Renaming the columns
data_pop <- rename(data_pop, city_name = "City",
                   rank = "2022rank",
                   population = "2022estimate")

# Removing superscripts in the city names
data_pop$city_name <- gsub("\\[.*\\]", "", data_pop$city_name)

# Removes commas and transform population figures into numbers
data_pop$population <- as.numeric(gsub(",", "", data_pop$population))

data_pop
```

Now we're ready to generate the figure:

```{r}
p <- ggplot(data_pop, aes(x=rank, y=population, label=city_name)) +
  geom_point() + geom_text(hjust=-.1, size=3) +
	scale_x_log10("log(rank)") + 
  scale_y_log10("log(population)", labels=scales::comma) +
  theme_minimal()
p
```

These power laws (https://en.wikipedia.org/wiki/Power_law) are remarkably general. For example, have a look how the first 10 million words in 30 Wikipedias (dumps from October 2015): https://en.wikipedia.org/wiki/Zipf%27s_law#/media/File:Zipf_30wiki_en_labels.png

Side note: We can also try to get an idea whether the distribution follows Zipf's law by estimating a log-log linear regression.

```{r}
summary(lm(log(rank) ~ log(population), data = data_pop))
```