Analysis on Game of Thrones marriage survival

In this notebook, I’ll …

• Download data from the “Song of Ice and Fire API”, https://anapioficeandfire.com/. A web API is a way to download data from the internet. The asoif API is a repository that contains a lot of GOT related data. There is an R package for that. https://github.com/MangoTheCat/GoTr
• Make data rectangling and transformations to ready data for analysis using the tidyverse https://www.tidyverse.org/
• I will create visualizations using the fantastic gameofhrones color palettes https://cran.r-project.org/web/packages/gameofthrones/index.html
• I will create Westeros maps using shapefiles from here: https://www.cartographersguild.com/showthread.php?t=30472
• I will make an animated plot using gganimate https://gganimate.com/
• I will conduct a survival analysis to see if married characters die sooner

First we have to download the data from the web. This API is “paginated”, which means that the data can be downloaded as it was a homepage with a fixed page size (e.g. 50 book characters in a page), and you have access to ‘next’ and ‘previous’ button links. So we need to build a loop to collect data from all characters.

We collected a whooping 2134 characters from the Game of Thrones universe!

In the next block, we transform the data from the awkward hierarchical data structure into a data frame, that is easy to use. We also clean the data to get meaningful variables.

Birth and death in the world of Game of Thrones

Births show a smooth pattern without sudden increases. The number of deaths seems to show a few sudden increases around 135 AC, 210AC, 280AC and most notably around 300AC. These likely coincide with significant event - such as wars or catastrophes - in the history of ASOIF.

Of course the conclusion is not: by 300 AC more people die than was born, because the data is heavily biased by missing values. As we get closer to the events of the books, more and more birth and death dates are documented. Still, it is telling that by 300, more deaths are documented than births…

cumulative_birth <-
    got_df %>% 
    group_by(year = birth_year) %>% 
    summarise(yearly_birth = sum(!is.na(birth_year))) %>% 
    ungroup() %>% 
    transmute(year, birth = cumsum(yearly_birth))

cumulative_death <-
    got_df %>% 
    group_by(year = death_year) %>% 
    summarise(yearly_death = sum(!is.na(death_year))) %>% 
    ungroup() %>% 
    transmute(year, death = cumsum(yearly_death)) %>% 
    drop_na()
    
    tibble(year = -28:300) %>% 
    left_join(cumulative_birth, by = "year") %>% 
    left_join(cumulative_death, by = "year") %>% 
    gather(event, number, -year, na.rm = TRUE) %>% 
        ggplot() +
        aes(x = year, y = number, color = event) +
        geom_line(size = 2) +
        scale_x_continuous(labels = scales::unit_format(unit = "A.C.")) +
        scale_color_got_d(option = "Martell") +
        labs(title = "Cumulative number of births and deaths over time for named characters in ASOIF",
             y = "Cumulative number of events",
             color = "Event")

Let’s plot a map of all the deaths and births on a proper Westeros map!

# Prepare the data for the animation
continents <- readOGR("map", "continents")

## OGR data source with driver: ESRI Shapefile 
## Source: "D:\Documents\GitHub\own projects\got_analysis\map", layer: "continents"
## with 3 features
## It has 2 fields
## Integer64 fields read as strings:  id

westeros <- readOGR("map", "political") # The actual map

## OGR data source with driver: ESRI Shapefile 
## Source: "D:\Documents\GitHub\own projects\got_analysis\map", layer: "political"
## with 12 features
## It has 3 fields
## Integer64 fields read as strings:  id

islands <- readOGR("map", "islands")

## OGR data source with driver: ESRI Shapefile 
## Source: "D:\Documents\GitHub\own projects\got_analysis\map", layer: "islands"
## with 86 features
## It has 3 fields
## Integer64 fields read as strings:  id

locations <- readOGR("map", "locations") # Places data

## OGR data source with driver: ESRI Shapefile 
## Source: "D:\Documents\GitHub\own projects\got_analysis\map", layer: "locations"
## with 247 features
## It has 5 fields
## Integer64 fields read as strings:  id size

places <-
    bind_cols(locations@data, as_tibble(locations@coords)) %>% 
    as_tibble() %>%
    select(id, place = name, 
           place_size = size, place_type = type, 
           lat = coords.x1, long = coords.x2) %>% 
    mutate_if(is.factor, as.character) %>%
    mutate_at(vars(id, place_size), as.integer) %>% 
    drop_na(place)

coordinates  <-
    got_df %>%
    drop_na(name) %>% 
    # Add more precise place names as rows, this will duplicate some rows, but will make
    # it possible to match more place names
    separate_rows(birth_place, sep = ",") %>% 
    separate_rows(death_place, sep = ",") %>% 
    # Add pefixes to birth and death info and join them to this dataset
    left_join(rename_all(places, ~paste0("birth_", .)), by = "birth_place") %>% 
    left_join(rename_all(places, ~paste0("death_", .)), by = "death_place") %>% 
    # Now, fill up the non-recognized coordinates by the recognized ones
    # arrange will make sure that the recognized will be on top when there are more places
    group_by(id) %>% 
    arrange(birth_long, death_long, .by_group = TRUE) %>% 
    fill(birth_long, birth_lat, death_long, death_lat) %>% 
    ungroup() %>% 
    # Keep only unique characters
    distinct(id, .keep_all = TRUE)

# Create a dataset with all notable events by time and place
death_birth <-
    coordinates %>% 
    select(id, name, matches("_year$|_lat$|_long$|_place$")) %>% 
    gather(variable, value, -id, -name, na.rm = TRUE) %>% 
    separate(variable, c("event", "variable"), sep = "_") %>% 
    spread(variable, value, convert = TRUE) %>% 
    drop_na() %>% 
    distinct(name, event, .keep_all = TRUE) %>% 
    mutate(decade = floor(year/10)*10) %>% 
    arrange(year)

event_places <-
    death_birth %>% 
    count(event, decade, lat, long, place) %>% 
    # Need to add decades that has no data for the animation later
    full_join(crossing(decade = seq(-30, 300, 10), 
                       event = c("birth", "death")), 
              by = c("event", "decade")) %>% 
    arrange(decade)

death_places <-
    event_places %>% 
    filter(event == "death") %>% 
    group_by(lat, long, place) %>% 
    summarise(n = sum(n)) %>% 
    drop_na()

# Plot a blank map first
base_map <-
    ggplot() +
    # Adds the continents
    geom_polygon(data = continents, 
                 aes(x = long, y = lat, group = group), 
                 color = "black", fill = "wheat") +
    # Adds the islands
    geom_polygon(data = islands, 
                 aes(x = long, y = lat, group = group),
                 color = "black", fill = "wheat") +
    # Adds the kingdomes of westeros
    geom_polygon(data = westeros,
                 aes(x = long, y = lat, group = group),
                 color = "black", fill = "wheat") +
    coord_map(xlim = c(-10, 75), ylim = c(-15, 50)) +
    theme_map() +
    theme(panel.background = element_rect(fill = "lightblue"))

## Regions defined for each Polygons
## Regions defined for each Polygons
## Regions defined for each Polygons

Let’s plot all the places where each notable death happened!

base_map +
    geom_point(data = death_places, 
               aes(x = lat, y = long, size = n)) +
    scale_color_got(option = "Greyjoy", direction = 1) +
    scale_size(range = c(2, 7)) +
    geom_text_repel(data = death_places, 
              aes(x = lat, y = long, label = place)) +
    labs(title = "Deadliest places in the word of ASOIF",
         size = "All notable deaths") +
    theme(legend.position = "right")

If we also wanted to include a timeline, we could do an animated map of all deaths.

event_map <-
    base_map +
    geom_point(data = filter(event_places, event == "death"),
               aes(x = lat,
                   y = long,
                   size = n)) +
    scale_color_got_d(option = "Targaryen", direction = -1) +
    scale_size(range = c(3, 7)) +
    geom_text(data = filter(event_places, event == "death"),
                    aes(x = lat,
                        y = long,
                        label = place)) +
    labs(title = "Deaths of characters over time in the {closest_state}s",
         color = "Event",
         size = "Notable deaths by decade") +
    theme(legend.position = "right") +
    transition_states(decade) +
    ease_aes("quadratic-in-out") +
    enter_grow() +
    exit_fade()

animate(event_map, nframes = 200, end_pause = 30, fps = 10,
        width = 900, height = 500)

Analysis question: Is there a smaller chance of survival in game of thrones after marriage?

People die left and right in game of thrones, but no occasion is as deadly as a wedding. But what about marriage? How people fare if they survived their wedding? In real life, lifetime expectancy is increased by marriage, which is obviously not the case in Game of Thrones. But exactly how dangerous is to be married in GOT?

We are building a cox regression on got character data with the outcome variable that shows if the character is alive or not at particular ages. We are only looking for main effects of marital status (if the character has a spouse). Disclaimer: This is obviously a very simple minded analysis, and I only did that to show how this kind of stuff can be done in R.

got_df %>% 
    select(is_alive, has_spouse, gender) %>% 
    ftable()

##                     gender Female Male
## is_alive has_spouse                   
## FALSE    FALSE                 52  402
##          TRUE                  50   81
## TRUE     FALSE                 79  144
##          TRUE                  32   31

To explore data, we create a box plot visualization about age by gender by marital status by survival status (based on 458 characters).

got_df %>% 
    drop_na(is_alive, has_spouse, gender, age) %>% 
    ggplot() +
      aes(y = age, x = gender, fill = gender) +
      geom_boxplot() +
      facet_grid(is_alive~has_spouse, labeller = label_both)

Now we create two statistical models (Cox regression) to test the effect of the predictors on the outcome variable (survival status).

survival_df <-
  got_df %>% 
  filter(!is.na(is_alive)) %>% 
  mutate(death_year = if_else(is_alive == TRUE & birth_year >= 200, 301, death_year),
         status = if_else(is_alive == TRUE, 0, 1),
         age = death_year - birth_year) %>% 
  drop_na(death_year)

survival_model1 <- coxph(Surv(age, status) ~ has_spouse, data = survival_df)
survival_model2 <- coxph(Surv(age, status) ~ has_spouse + gender, data = survival_df)
    
AIC(survival_model1, survival_model2)

The models show that being married is a significant health risk in the world of ASOIF. Comparing the models shows that gender does not add much to this finding, therefore we will just visualize the one without gender.

We also calculate odds ratios, and confidence intervals. Odds ratios suggest that it is 1.4 95%CI[1.03 1.88] times more likely to die if a character is married.

At first, there is little difference between characters with and without spouses, and young characters with spouses have a better chance of survival but after about 35 years of age,those who have a spouse die earlier.

Not getting wed is a significant predictor of longer life. Those having a spouse die at a median age of 48 years 95%CI[54 80], while the unmarried die at 66 years 95%CI[41 53]. This is a median difference of 18 years!

summary(survival_model1)

## Call:
## coxph(formula = Surv(age, status) ~ has_spouse, data = survival_df)
## 
##   n= 458, number of events= 172 
##    (331 observations deleted due to missingness)
## 
##                  coef exp(coef) se(coef)     z Pr(>|z|)  
## has_spouseTRUE 0.3285    1.3888   0.1544 2.127   0.0334 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                exp(coef) exp(-coef) lower .95 upper .95
## has_spouseTRUE     1.389       0.72     1.026      1.88
## 
## Concordance= 0.489  (se = 0.021 )
## Likelihood ratio test= 4.52  on 1 df,   p=0.03
## Wald test            = 4.52  on 1 df,   p=0.03
## Score (logrank) test = 4.56  on 1 df,   p=0.03

And finally, let’s plot the predictions for the survival chance over time by marital status.

survfit(Surv(age, status) ~ has_spouse, data = survival_df) %>% 
    ggsurv(.,
           cens.shape = 16, 
           size.est = 1.2) +
    scale_color_got_d(option = "Greyjoy", direction = 1) +
    scale_y_continuous(labels = scales::percent_format()) +
    labs(title = "Survival of ASOIF characters based on marital status",
         subtitle = "At first, there is little difference between characters with and without spouses, and young characters with spouses have\na better chance of survival but after about 35 years of age,those who have a spouse die earlier",
         y = "Estimated % survival", 
         x = "Age") +
    annotate("text",
             label = c("Has spouse", "No spouse"),
             x = c(100, 100),
             y = c(.1, .25),
             color = got(option = "Greyjoy", direction = 1, n = 2)) +
    guides(color = FALSE, linetype = FALSE)

## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

Conclusions

• we looked at significant life events (births and deaths) of notable characters in ASOIF
• Married characters die approximately 18 years sooner than unmarried characters