Analysis of demographics

library(readr)
library(dplyr)
library(scales)
library(tidyr)
library(ggplot2)
library(cowplot)
library(stringr)
library(forcats)
library(ggallin)

theme_set(theme_bw())

Read in data

alpha_threshold <- qnorm(0.975)
countries <- c("France", "Germany", "Spain", "Italy", "Singapore", "USA" )
countries_code <- c("FR", "DE", "ES", "IT", "SG", "US")
colorBlindGrey6 <- c("#0072B2", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#CC79A7")
races <- c("black", "hispanic_latino", "white", "asian", 
           "hawaiian_pacific_islander", "american_indian")  
Races <- rev(c('White', 'Black', 'Hispanic/Latino', 'Asian', 
               'American Indian', 'Hawaiian/Pacific Islander', 'Other/Unknown'))

demo_desc <- read_csv('results/demo_prelim.csv') %>% 
  filter(siteid != 'SITE157', siteid != 'SITE712') %>% 
  mutate(Country = fct_relevel(Country, countries))

## 
## ── Column specification ────────────────────────────────────────────────────────
## cols(
##   .default = col_double(),
##   siteid = col_character(),
##   Country = col_character()
## )
## ℹ Use `spec()` for the full column specifications.

demo_perc <- demo_desc %>% 
  mutate_at(vars(- num_patients_all, - siteid, - Country, - median_age),
            list(~ . / .data$num_patients_all)) %>% 
  mutate(Country = as.factor(Country)) 

demo_perc %>% 
  rename('perc_patients_ever_severe' = num_patients_ever_severe,
         'perc_patients_never_severe' = num_patients_never_severe) %>% 
  write_csv('results/demographics.csv')

Overall statistics

demo_desc %>% 
  summarise(across(c(num_patients_all, female, male, unknown_other_sex), sum)) 

## # A tibble: 1 x 4
##   num_patients_all female  male unknown_other_sex
##              <dbl>  <dbl> <dbl>             <dbl>
## 1            34647  13546 20814               287

demo_desc %>% 
  summarise(across(c(num_patients_all, races, unknown_other_race), sum)) %>% 
  print()

## Note: Using an external vector in selections is ambiguous.
## ℹ Use `all_of(races)` instead of `races` to silence this message.
## ℹ See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
## This message is displayed once per session.

## # A tibble: 1 x 8
##   num_patients_all black hispanic_latino white asian hawaiian_pacifi…
##              <dbl> <dbl>           <dbl> <dbl> <dbl>            <dbl>
## 1            34647  5706            1400  8532   621                9
## # … with 2 more variables: american_indian <dbl>, unknown_other_race <dbl>

unknowns <- demo_perc %>%
  select(siteid, Country, contains('unknown')) %>% 
  pivot_longer(- c(siteid, Country))
unknowns %>% filter(value < 0)

## # A tibble: 6 x 4
##   siteid  Country name                  value
##   <chr>   <fct>   <chr>                 <dbl>
## 1 SITE824 USA     unknown_other_sex  -0.00293
## 2 SITE864 USA     unknown_other_race -0.00443
## 3 SITE295 USA     unknown_other_race -0.0822 
## 4 SITE200 USA     unknown_other_race -0.0145 
## 5 SITE775 USA     unknown_other_race -0.0769 
## 6 SITE923 USA     unknown_other_race -0.107

unknowns %>% filter(value < 1, value > 0.1)

## # A tibble: 9 x 4
##   siteid  Country   name               value
##   <chr>   <fct>     <chr>              <dbl>
## 1 SITE824 USA       unknown_other_race 0.152
## 2 SITE599 USA       unknown_other_race 0.14 
## 3 SITE116 Singapore unknown_other_race 0.378
## 4 SITE798 USA       unknown_other_race 0.173
## 5 SITE510 USA       unknown_other_race 0.367
## 6 SITE533 USA       unknown_other_race 0.148
## 7 SITE432 USA       unknown_other_sex  0.107
## 8 SITE432 USA       unknown_other_race 0.214
## 9 SITE734 USA       unknown_other_race 0.258

# unknowns %>% 
#   ggplot(aes(y = value, x = name)) +
#   ggbeeswarm::geom_beeswarm() +
#   facet_wrap(~ name, scales = 'free_x')

Severity vs. gender/age/race?

demo_gender <- demo_perc %>% filter(unknown_other_sex < 1)
lm(num_patients_ever_severe ~ female, demo_gender) %>% 
  summary()

## 
## Call:
## lm(formula = num_patients_ever_severe ~ female, data = demo_gender)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.55272 -0.04966 -0.00420  0.10740  0.42322 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.00420    0.04687   21.43  < 2e-16 ***
## female      -1.32622    0.13873   -9.56 9.02e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1833 on 39 degrees of freedom
## Multiple R-squared:  0.7009, Adjusted R-squared:  0.6932 
## F-statistic: 91.38 on 1 and 39 DF,  p-value: 9.02e-12

age_perc <- demo_perc %>% 
  filter(unknown_other_age_group < 1) %>% 
  select(Country, num_patients_all, num_patients_ever_severe, contains('to'), ends_with('plus'))

lm(num_patients_ever_severe ~ . - Country - num_patients_all, age_perc) %>% 
  summary()

## 
## Call:
## lm(formula = num_patients_ever_severe ~ . - Country - num_patients_all, 
##     data = age_perc)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.56212 -0.07390  0.03572  0.10181  0.38724 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   -3.731      1.712  -2.179  0.03659 * 
## `00to02`      -4.334      7.509  -0.577  0.56772   
## `03to05`      73.690     33.517   2.199  0.03503 * 
## `06to11`     -18.079     15.938  -1.134  0.26484   
## `12to17`       4.277     15.450   0.277  0.78362   
## `18to25`       1.723      3.225   0.534  0.59675   
## `26to49`       3.803      1.815   2.096  0.04386 * 
## `50to69`       4.476      1.902   2.354  0.02470 * 
## `70to79`       6.186      1.887   3.278  0.00246 **
## `80plus`       2.964      1.638   1.810  0.07943 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.217 on 33 degrees of freedom
## Multiple R-squared:  0.6591, Adjusted R-squared:  0.5661 
## F-statistic: 7.089 on 9 and 33 DF,  p-value: 1.225e-05

lm(num_patients_ever_severe ~ black + white + asian + hispanic_latino + 
     american_indian + hawaiian_pacific_islander, 
   data = filter(demo_perc, unknown_other_race < 1)) %>% 
  summary()

## 
## Call:
## lm(formula = num_patients_ever_severe ~ black + white + asian + 
##     hispanic_latino + american_indian + hawaiian_pacific_islander, 
##     data = filter(demo_perc, unknown_other_race < 1))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.5508 -0.1224  0.1029  0.1451  0.3034 
## 
## Coefficients:
##                           Estimate Std. Error t value Pr(>|t|)
## (Intercept)                 0.2270     0.6028   0.377    0.709
## black                       0.7031     0.6360   1.105    0.279
## white                       0.6638     0.6880   0.965    0.344
## asian                      -0.4146     1.1398  -0.364    0.719
## hispanic_latino            -0.3319     0.6205  -0.535    0.597
## american_indian           -31.0987    33.4578  -0.929    0.361
## hawaiian_pacific_islander -60.7964    71.6772  -0.848    0.404
## 
## Residual standard error: 0.2678 on 26 degrees of freedom
## Multiple R-squared:  0.3673, Adjusted R-squared:  0.2212 
## F-statistic: 2.515 on 6 and 26 DF,  p-value: 0.04703

lm(num_patients_ever_severe ~ median_age, demo_perc) %>% 
  summary()

## 
## Call:
## lm(formula = num_patients_ever_severe ~ median_age, data = demo_perc)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.69883 -0.22095 -0.01098  0.26950  0.45614 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -0.417341   0.452308  -0.923   0.3616  
## median_age   0.015701   0.006709   2.340   0.0242 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3132 on 41 degrees of freedom
## Multiple R-squared:  0.1178, Adjusted R-squared:  0.09633 
## F-statistic: 5.477 on 1 and 41 DF,  p-value: 0.02421

Overall plot (A)

gender_count <- demo_desc %>% 
  # filter(unknown_other_sex/num_patients_all < 1) %>% 
  # mutate(siteid = as.factor(siteid)) %>% 
  # mutate(num_patients_all_gender = female + male) %>% 
  group_by(Country) %>% 
  summarise(num_patients_country = sum(num_patients_all),
            male_count = sum(male),
            female_count = sum(female),
            .groups = 'drop') %>% 
  {.}

showtext::showtext_auto()
female <- intToUtf8(9792)
male <- intToUtf8(9794)

overall_plot <- gender_count %>% 
  mutate(male_count = -male_count) %>% 
  pivot_longer(- c(Country, num_patients_country)) %>% 
  ggplot(aes(y = fct_relevel(Country, rev(countries)), 
             x = value, color = name, 
             fill = Country)) +
  geom_col(position = 'identity') +
  scale_color_manual(values = c('white', 'white')) +
  scale_x_continuous(breaks = c(-5000, -2000, 0, 2000, 5000),
                     labels = c(5000, 2000, 0, 2000, 5000),
                     trans = ssqrt_trans) +
  geom_vline(xintercept = 0, color = 'white') +
  geom_text(data = . %>% filter(name == 'male_count'),
            aes(y = Country, x = value - 22*sqrt(abs(value)), 
                label = abs(value)), color = 'grey40', size = 3) +
    geom_text(data = . %>% filter(name == 'female_count'),
            aes(y = Country, x = value + 19*sqrt(abs(value)), 
                label = abs(value)), color = 'grey40', size = 3) +
  geom_text(data = filter(gender_count, Country == 'France'),
            aes(y = 6, x = female_count - 600), fontface = 2,
            label = female, size = 3, color = 'white') +
  geom_text(data = filter(gender_count, Country == 'France'),
            aes(y = 6, x = - male_count + 600), fontface = 2,
            label = male, size = 3, color = 'white') +
  labs(y = NULL, 
       x = 'Number of male/female patients per country', 
       fill = NULL) +
  scale_fill_manual(values = colorBlindGrey6) +
  theme(legend.position = 'None',
        plot.margin = margin(15, 5.5, 15.5, 55.5),
        axis.ticks.y = element_blank(),
        panel.grid.major.y = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.spacing.y = unit(0, 'pt'),
        strip.background = element_blank(),
        strip.placement = 'inside',
        legend.background = element_blank(),
        legend.key.height = unit(1, 'line'),
        legend.title = element_text(hjust = 0, size = 9)) +
  NULL

Age plot (B)

age_count <- demo_desc %>% 
  filter(unknown_other_age_group/num_patients_all < 1) %>% 
  select(Country, num_patients_all, contains('to'), ends_with('plus')) %>% 
  group_by(Country) %>% 
  summarise(across(.fns = sum), .groups = 'drop') %>% 
  mutate_at(vars(- num_patients_all, - Country),
            list(~ . / .data$num_patients_all)) %>% 
  select(-num_patients_all) %>% 
  pivot_longer(- Country) %>% 
  mutate(name = case_when(
    name == '00to02' ~ '0to2',
    name == '03to05' ~ '3to5',
    name == '06to11' ~ '6to11',
    TRUE ~ name
  )) %>% 
  mutate(name = gsub('to', '–', name) %>% gsub('plus', '+', .))
age_count$name <- factor(age_count$name, levels = unique(age_count$name))

age_plot <- age_count  %>% 
  ggplot(aes(name, value, fill = Country, group = Country)) +
  geom_col(position = 'dodge2') + 
  coord_flip() +
  facet_grid(rows = vars(Country)) +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1),
                     breaks = seq(0, 1, 0.2),
                     expand = expansion(add = c(0,0.02))) +
  scale_size_continuous(guide = FALSE) +
  labs(x = NULL, y = 'Patients proportion per age group') +
  scale_color_manual(values = colorBlindGrey6, drop = FALSE) +
  scale_fill_manual(values = colorBlindGrey6, drop = FALSE) +
  theme(legend.position = 'None',
        axis.text.y = element_text(size = 7),
        panel.grid.major.y = element_blank(),
        plot.margin = margin(15, 5.5, 5.5, 5.5, 'pt')) +
  NULL

Severity vs. age plot (C)

country_codes <- read_csv('https://raw.githubusercontent.com/lukes/ISO-3166-Countries-with-Regional-Codes/master/all/all.csv') %>%
  mutate(name = ifelse(name == 'United States of America', 'USA', name)) %>%
  filter(name %in% countries) %>%
  select(Country = name, coun_code = `alpha-2`) %>%
  mutate(coun_code = as.factor(coun_code) %>% fct_relevel(countries_code))

## 
## ── Column specification ────────────────────────────────────────────────────────
## cols(
##   name = col_character(),
##   `alpha-2` = col_character(),
##   `alpha-3` = col_character(),
##   `country-code` = col_character(),
##   `iso_3166-2` = col_character(),
##   region = col_character(),
##   `sub-region` = col_character(),
##   `intermediate-region` = col_character(),
##   `region-code` = col_character(),
##   `sub-region-code` = col_character(),
##   `intermediate-region-code` = col_character()
## )

age_severity <- demo_perc %>% 
  left_join(country_codes, by = 'Country') %>%
  ggplot(aes(median_age, num_patients_ever_severe, 
             color = coun_code)) +
  geom_point(aes(size = num_patients_all), alpha = 0.8) +
  scale_y_continuous(labels = scales::percent_format(),
                     breaks = seq(0, 1, 0.2)) +
  labs(x = 'Median age estimate', 
       y = 'Severe case proportion', 
       color = NULL) +
  scale_size(trans = "log10", guide = FALSE) +
  scale_color_manual(values = colorBlindGrey6) +
  coord_cartesian(ylim = c(0, 1)) +
  theme(
    # legend.position = c(0.9, 0.7),
    legend.position = 'left',
    legend.background = element_blank(),
    plot.margin = margin(5.5, 5.5, 10.5, 15.5),
    panel.grid.minor = element_blank(),
    legend.key.height = unit(0.7, 'line'),
    legend.key.width = unit(0.4, 'line'),
    legend.title = element_text(hjust = 0, size = 9)) +
  NULL

Race plot (D)

race_count <- demo_desc %>% 
  filter(unknown_other_race/num_patients_all < 1) %>% 
  select(Country, num_patients_all, all_of(races), unknown_other_race) %>% 
  group_by(Country) %>% 
  summarise(across(.fns = sum), .groups = 'drop') %>% 
  mutate_at(vars(- num_patients_all, - Country),
            list(~ . / .data$num_patients_all)) %>% 
  select(-num_patients_all) %>% 
  pivot_longer(- Country) %>% 
  mutate(name = gsub('_', ' ', name) %>% 
           str_to_title(.) %>% 
           gsub('Hispanic ', 'Hispanic/', .) %>% 
           gsub('Hawaiian ', 'Hawaiian/', .) %>% 
           gsub('Unknown Other Race', 'Other/Unknown',. ),
         name = fct_relevel(name, Races))

race_plot <- race_count %>% 
  # remove one site in Italy that reported all white ethnicity
  filter(Country != 'Italy') %>% 
  ggplot(aes(name, value, fill = Country, group = Country)) +
  geom_col(position = 'dodge2') + 
  coord_flip() +
  facet_grid(rows = vars(Country)) +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1),
                     breaks = seq(0, 1, 0.2),
                     expand = expansion(add = c(0,0.02))) +
  scale_size_continuous(guide = FALSE) +
  labs(x = NULL, y = 'Patient proportion per race/ethnicity group') +
  scale_color_manual(values = colorBlindGrey6, drop = FALSE) +
  scale_fill_manual(values = colorBlindGrey6, drop = FALSE) +
  theme(legend.position = 'None',
        axis.text.y = element_text(size = 8),
        panel.grid.major.y = element_blank(),
        axis.title.x = element_text(hjust = .8)) +
  NULL

Composite plot

gender_race <- plot_grid(overall_plot, age_severity, race_plot, ncol = 1,
          rel_heights = c(0.7, 1, 1), 
          labels = c('A', 'B', 'C'), axis = 'r')
abcd <- plot_grid(gender_race, age_plot, ncol = 2,
                  rel_widths = c(1.25, 1 ), labels = c('', 'D'),
                  axis = 'b')
abcd

ggsave('figs/composite_plot.png', abcd, height = 6.5, width = 9)
ggsave('figs/tiffs/fig_2.tiff', abcd, height = 6.5, width = 9, dpi = 300)