From 8495eb70fa91926c8d1ad79d2651ef1f8dd9634d Mon Sep 17 00:00:00 2001
From: lelaboratoire If we set a threshold at least 200 articles a year, we should only consider articles from 1998 on. Number of keynote speakers/fellows across years: Honorees that didn’t receive a gender prediction: Chung-I Wu. In summary, the NA predictions mostly include initials only, hyphenated names and perhaps names with accent marks. Only keep articles from 2002 because few authors had gender predictions before 2002. See 093.summary-stats for more details. The two groups of scientists did not have a significant association with the gender predicted from fore names (P = 0.095815). Interaction terms do not predict The two groups of scientists did not have a significant association with the gender predicted from fore names (P = 0.091898). Interaction terms do not predict Only keep articles from 2002 because few authors had nationality predictions before 2002 (mostly due to missing metadata). See 093.summary-stats for more details. Names grouping: Prepare data frames for later analyses: Compared to the name collection of Pubmed authors, honorees with Celtic/English names are overrepresented while honorees with East Asian names are underrepresented. Interaction terms do not predict A Celtic/English name has 1.2117973 the odds of being selected as an honoree, significantly higher compared to other names (\(\beta_\textrm{Celtic/English} =\) 0.1921, P = 0.00067307). An East Asian name has 0.1891745 the odds of being selected as an honoree, significantly lower than to other names (\(\beta_\textrm{East Asian} =\) -1.6651, P = 5.89e-09). The two groups of scientists did not have a significant association with names predicted to be European (P = 0.37934) or in Other categories (P = 0.43538). It’s difficult to come to a conclusion for other regions with so few data points and the imperfect accuracy of our prediction. There seems to be little difference between the proportion of keynote speakers of African, Arabic, South Asian and Hispanic origin than those in the field. However, just because a nationality isn’t underrepresented against the field doesn’t mean scientists from that nationality are appropriately represented. A Celtic/English name has 1.2141832 the odds of being selected as an honoree, significantly higher compared to other names (\(\beta_\textrm{Celtic/English} =\) 0.19407, P = 0.00054646). An East Asian name has 0.1889996 the odds of being selected as an honoree, significantly lower than to other names (\(\beta_\textrm{East Asian} =\) -1.666, P = 6.0164e-09). The two groups of scientists did not have a significant association with names predicted to be European (P = 0.38583) or in Other categories (P = 0.44544). Sincere thanks to the reviewers, Byron Smith and Katie Pollard, for their detailed suggestion with code. The question of what unit one should use to perform this type of analyses is a difficult one. We present here an alternative analysis that treats names as units instead of honors and authorships. We caution that this approach does not distinguish scientists who were honored 4 times vs. one time and hence may yield a conservative estimate of disparity. Further, different authors may have the same names, and to sum Nonetheless, the finding here is consistent with what we have seen above where East Asian names are underrepresented in the honoree group. In this section, we show that a 10-year lag model results in a similar result for East Asian scientists’ names, even though the effect size is less striking. For example, if we assume that honors accrue 10 years after their most prolific year with respect to authorships, the proportion of honor associated with East Asian name origins in 2019 is still substantially less than the proportion of senior authorships associated with East Asian names in 2009. It’s difficult to come to a conclusion for other regions with so few data points and the imperfect accuracy of our prediction. There seems to be little difference between the proportion of keynote speakers of African, Arabic, South Asian and Hispanic origin than those in the field. However, just because a nationality isn’t underrepresented against the field doesn’t mean scientists from that nationality are appropriately represented. Adapted from Plotting ROC curves
Name origin prediction method performance
-
+## ── Attaching packages ────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──
-## ✓ ggplot2 3.3.2 ✓ purrr 0.3.4
## ✓ tibble 3.0.4 ✓ dplyr 1.0.2
## ✓ tidyr 1.1.2 ✓ stringr 1.4.0
## ✓ readr 1.4.0 ✓ forcats 0.5.0
## ── Conflicts ───────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
+
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
# still need to install caret for the calibration function because tidymodels's
@@ -1769,7 +1769,7 @@
Name origin prediction method performance
mutate(Sensitivity = tpr, Specificity = 1-fpr, dSens = c(abs(diff(1-tpr)), 0)) %>%
ungroup()##
-## ── Column specification ───────────────────────────────────────────────────────────────────────────────
+## ── Column specification ────────────────────────────────────────────────────────
## cols(
## fpr = col_double(),
## tpr = col_double(),
@@ -1811,7 +1811,7 @@
Name origin prediction method performance
mutate(y_true = as.factor(truth)) %>%
select(-truth)
-##
-## ── Column specification ───────────────────────────────────────────────────────────────────────────────
+## ── Column specification ────────────────────────────────────────────────────────
## cols(
## African = col_double(),
## CelticEnglish = col_double(),
@@ -1844,18 +1844,30 @@
Name origin prediction method performance
## ℹ See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
## This message is displayed once per session.
+## calibModelVar bin Percent Lower Upper Count midpoint region
-## 1 prob [0,0.0909] 0.973038 0.9061138 1.043559 777 4.545455 EastAsian
-## 2 prob (0.0909,0.182] 12.715105 10.7555376 14.887108 133 13.636364 EastAsian
-## 3 prob (0.182,0.273] 20.620843 16.9791523 24.652952 93 22.727273 EastAsian
-## 4 prob (0.273,0.364] 29.924242 24.4643714 35.841394 79 31.818182 EastAsian
-## 5 prob (0.364,0.455] 35.897436 29.7515540 42.405681 84 40.909091 EastAsian
-## 6 prob (0.455,0.545] 38.536585 31.8402892 45.569554 79 50.000000 EastAsian
-## 7 prob (0.545,0.636] 45.637584 37.4635833 53.988516 68 59.090909 EastAsian
-## 8 prob (0.636,0.727] 56.953642 48.6544756 64.974492 86 68.181818 EastAsian
-## 9 prob (0.727,0.818] 61.421320 54.2394760 68.253900 121 77.272727 EastAsian
-## 10 prob (0.818,0.909] 71.764706 66.6571343 76.488532 244 86.363636 EastAsian
-## 11 prob (0.909,1] 97.209555 96.8524649 97.536348 8953 95.454545 EastAsian
## calibModelVar bin Percent Lower Upper Count midpoint
+## 1 prob [0,0.0909] 0.973038 0.9061138 1.043559 777 4.545455
+## 2 prob (0.0909,0.182] 12.715105 10.7555376 14.887108 133 13.636364
+## 3 prob (0.182,0.273] 20.620843 16.9791523 24.652952 93 22.727273
+## 4 prob (0.273,0.364] 29.924242 24.4643714 35.841394 79 31.818182
+## 5 prob (0.364,0.455] 35.897436 29.7515540 42.405681 84 40.909091
+## 6 prob (0.455,0.545] 38.536585 31.8402892 45.569554 79 50.000000
+## 7 prob (0.545,0.636] 45.637584 37.4635833 53.988516 68 59.090909
+## 8 prob (0.636,0.727] 56.953642 48.6544756 64.974492 86 68.181818
+## 9 prob (0.727,0.818] 61.421320 54.2394760 68.253900 121 77.272727
+## 10 prob (0.818,0.909] 71.764706 66.6571343 76.488532 244 86.363636
+## 11 prob (0.909,1] 97.209555 96.8524649 97.536348 8953 95.454545
+## region
+## 1 EastAsian
+## 2 EastAsian
+## 3 EastAsian
+## 4 EastAsian
+## 5 EastAsian
+## 6 EastAsian
+## 7 EastAsian
+## 8 EastAsian
+## 9 EastAsian
+## 10 EastAsian
+## 11 EastAsian
fig_3b <- bind_rows(cal_dfs) %>%
recode_region() %>%
ggplot(aes(x = midpoint/100, y = Percent/100, color = fct_relevel(region, as.character(auc_df$region)))) +
diff --git a/docs/092.save-raws-to-Rdata.html b/docs/092.save-raws-to-Rdata.html
index aad4c8d..6a6034e 100644
--- a/docs/092.save-raws-to-Rdata.html
+++ b/docs/092.save-raws-to-Rdata.html
@@ -4300,7 +4300,7 @@
Save raw data into .Rdata to be loaded in by analys
General data read-in
##
-## ── Column specification ───────────────────────────────────────────────────────────────────────────────
+## ── Column specification ────────────────────────────────────────────────────────
## cols(
## fore_name = col_character(),
## last_name = col_character(),
@@ -4334,7 +4334,7 @@
General data read-in
publication_date = ymd(publication_date, truncated = 2)) %>%
filter(year(publication_date) < 2020)##
-## ── Column specification ───────────────────────────────────────────────────────────────────────────────
+## ── Column specification ────────────────────────────────────────────────────────
## cols(
## pmid = col_double(),
## pmcid = col_character(),
@@ -4358,7 +4358,7 @@
General data read-in
left_join(select(all_full_names, - full_name), by = c('fore_name', 'last_name')) %>%
filter(year(year) < 2020, conference != 'PSB') # remove PSB, exclude ISCB Fellows and ISMB speakers in 2020 for now##
-## ── Column specification ───────────────────────────────────────────────────────────────────────────────
+## ── Column specification ────────────────────────────────────────────────────────
## cols(
## year = col_double(),
## full_name = col_character(),
@@ -4371,9 +4371,10 @@
General data read-in
## )
+## # … with 11 variables: year <date>, full_name <chr>, fore_name <chr>,
+## # last_name <chr>, conference <chr>, source <chr>, affiliations <chr>,
+## # afflcountries <chr>, publication_date <date>, fore_name_simple <chr>,
+## # last_name_simple <chr>
## # A tibble: 0 x 11
-## # … with 11 variables: year <date>, full_name <chr>, fore_name <chr>, last_name <chr>,
-## # conference <chr>, source <chr>, affiliations <chr>, afflcountries <chr>, publication_date <date>,
-## # fore_name_simple <chr>, last_name_simple <chr>
large_jours <- articles %>%
count(journal, sort = T) %>%
head(10)
@@ -4384,7 +4385,7 @@
General data read-in
left_join(all_full_names, by = 'full_name')## Warning: Missing column names filled in: 'X1' [1]
##
-## ── Column specification ───────────────────────────────────────────────────────────────────────────────
+## ── Column specification ────────────────────────────────────────────────────────
## cols(
## X1 = col_character(),
## African = col_double(),
@@ -4405,8 +4406,8 @@
General data read-in
corr_authors %>%
count(year, name = 'Number of articles with last authors') %>%
DT::datatable(rownames = F)
corr_authors <- corr_authors %>%
add_count(year, name = 'n_aut_yr') %>%
diff --git a/docs/093.summary-stats.html b/docs/093.summary-stats.html
index 20df43e..54e91c6 100644
--- a/docs/093.summary-stats.html
+++ b/docs/093.summary-stats.html
@@ -4317,8 +4317,8 @@
Honorees
count(fore_name, last_name) %>%
arrange(desc(n)) %>%
DT::datatable()keynotes %>%
select(year, conference) %>%
@@ -4364,7 +4364,7 @@
Authors
Gender analysis
##
-## ── Column specification ───────────────────────────────────────────────────────────────────────────────
+## ── Column specification ────────────────────────────────────────────────────────
## cols(
## fore_name_simple = col_character(),
## n_authors = col_double(),
@@ -4412,7 +4412,8 @@
Gender analysis
pull(pred.asi) %>%
mean(na.rm = T)
-## [1] "Proceeding with surname-only predictions..."
+## Warning in merge_surnames(voter.file, impute.missing = impute.missing): 1305 surnames were not matched.
## Warning in merge_surnames(voter.file, impute.missing = impute.missing): 1305
+## surnames were not matched.
## [1] 0.8174599
Name origin analysis
## ℹ Unknown levels in `f`: OtherCategories
## ℹ Input `region` is `fct_recode(...)`.
-
-
+
+
## Warning: Unknown levels in `f`: OtherCategories
## # A tibble: 2 x 3
## `!is.na(African)` `is.na(fore_name_simple.x)` n
diff --git a/docs/10.visualize-gender.html b/docs/10.visualize-gender.html
index c8c3c08..1f524c3 100644
--- a/docs/10.visualize-gender.html
+++ b/docs/10.visualize-gender.html
@@ -1671,18 +1671,18 @@
Setups
library(lubridate)
library(rnaturalearth)
library(wru)
-source('utils/r-utils.R')
+source("utils/r-utils.R")
theme_set(theme_bw() + theme(legend.title = element_blank()))Load data
load('Rdata/raws.Rdata')
+
+gender_df <- read_tsv("data/gender/genderize.tsv")load("Rdata/raws.Rdata")
alpha_threshold <- qnorm(0.975)
-gender_df <- read_tsv('data/gender/genderize.tsv')
##
-## ── Column specification ───────────────────────────────────────────────────────────────────────────────
+## ── Column specification ────────────────────────────────────────────────────────
## cols(
## fore_name_simple = col_character(),
## n_authors = col_double(),
@@ -1691,11 +1691,11 @@
Load data
## probability_male = col_double()
## )pubmed_gender_df <- corr_authors %>%
- filter(year(year) >= 2002) %>%
- left_join(gender_df, by = 'fore_name_simple')
+ filter(year(year) >= 2002) %>%
+ left_join(gender_df, by = "fore_name_simple")
-iscb_gender_df <- keynotes %>%
- left_join(gender_df, by = 'fore_name_simple')
+iscb_gender_df <- keynotes %>%
+ left_join(gender_df, by = "fore_name_simple")
start_year <- 1993
end_year <- 2019
@@ -1713,53 +1713,56 @@
Prepare data frames for later analyses
iscb_pubmed <- iscb_gender_df %>%
- rename('journal' = conference) %>%
+ rename("journal" = conference) %>%
select(year, journal, probability_male, publication_date) %>%
- mutate(type = 'Keynote speakers/Fellows',
- adjusted_citations = 1) %>%
+ mutate(
+ type = "Keynote speakers/Fellows",
+ adjusted_citations = 1
+ ) %>%
bind_rows(
pubmed_gender_df %>%
select(year, journal, probability_male, publication_date, adjusted_citations) %>%
- mutate(type = 'Pubmed authors')
+ mutate(type = "Pubmed authors")
) %>%
- mutate(probability_female = 1 - probability_male) %>%
- pivot_longer(contains('probability'),
- names_to = 'gender',
- values_to = 'probabilities') %>%
- filter(!is.na(probabilities)) %>%
+ mutate(probability_female = 1 - probability_male) %>%
+ pivot_longer(contains("probability"),
+ names_to = "gender",
+ values_to = "probabilities"
+ ) %>%
+ filter(!is.na(probabilities)) %>%
group_by(type, year, gender) %>%
mutate(
pmc_citations_year = mean(adjusted_citations),
- weight = adjusted_citations/pmc_citations_year,
- weighted_probs = probabilities*weight
+ weight = adjusted_citations / pmc_citations_year,
+ weighted_probs = probabilities * weight
# weight = 1
- )
+ )
-iscb_pubmed_sum <- iscb_pubmed %>%
+iscb_pubmed_sum <- iscb_pubmed %>%
summarise(
# n = n(),
mean_prob = mean(weighted_probs),
# mean_prob = mean(probabilities, na.rm = T),
# sd_prob = sd(probabilities, na.rm = T),
- se_prob = sqrt(var(probabilities) * sum(weight^2)/(sum(weight)^2)),
+ se_prob = sqrt(var(probabilities) * sum(weight^2) / (sum(weight)^2)),
# n = mean(n),
me_prob = alpha_threshold * se_prob,
- .groups = 'drop'
+ .groups = "drop"
)
# https://stats.stackexchange.com/questions/25895/computing-standard-error-in-weighted-mean-estimation
# save(iscb_pubmed, file = 'Rdata/iscb-pubmed_gender.Rdata')
Figures for paper
-
Mean and standard deviation of predicted probabilities
-iscb_pubmed_sum %>% filter(gender == 'probability_male') %>%
- gam_and_ci(df2 = iscb_pubmed %>% filter(gender == 'probability_male'),
- start_y = start_year, end_y = end_year) +
+
iscb_pubmed_sum %>%
+ filter(gender == "probability_male") %>%
+ gam_and_ci(
+ df2 = iscb_pubmed %>% filter(gender == "probability_male"),
+ start_y = start_year, end_y = end_year
+ ) +
theme(legend.position = c(0.88, 0.2))
## `geom_smooth()` using formula 'y ~ s(x, bs = "cs")'
@@ -1801,220 +1814,53 @@ ## Warning: Removed 13171 rows containing non-finite values (stat_smooth).
Mean and standard deviation of predicted probabilities
Hypothesis testing
-
+iscb_lm <- iscb_pubmed %>%
- filter(gender == 'probability_female', !is.na(weighted_probs)) %>%
+
-iscb_lm <- iscb_pubmed %>%
+ filter(gender == "probability_female", !is.na(weighted_probs)) %>%
mutate(type = as.factor(type)) %>%
- mutate(type = relevel(type, ref = 'Pubmed authors'),
- year = as.factor(year))
-
-main_lm <- glm(type ~ year + weighted_probs,
- data = iscb_lm, family = 'binomial')
-
-summary(main_lm)
-##
-## Call:
-## glm(formula = type ~ year + weighted_probs, family = "binomial",
-## data = iscb_lm)
-##
-## Deviance Residuals:
-## Min 1Q Median 3Q Max
-## -0.2521 -0.0675 -0.0602 -0.0536 3.7156
-##
-## Coefficients:
-## Estimate Std. Error z value Pr(>|z|)
-## (Intercept) 1.655e+01 1.385e+03 0.012 0.9905
-## year1994-01-01 8.750e-03 1.959e+03 0.000 1.0000
-## year1995-01-01 -3.840e-02 1.959e+03 0.000 1.0000
-## year1996-01-01 5.925e-03 1.833e+03 0.000 1.0000
-## year1997-01-01 -6.766e-03 1.537e+03 0.000 1.0000
-## year1998-01-01 -1.777e-03 1.549e+03 0.000 1.0000
-## year1999-01-01 8.430e-03 1.510e+03 0.000 1.0000
-## year2000-01-01 -1.929e-03 1.510e+03 0.000 1.0000
-## year2001-01-01 -1.254e-03 1.510e+03 0.000 1.0000
-## year2002-01-01 -2.118e+01 1.385e+03 -0.015 0.9878
-## year2003-01-01 -2.150e+01 1.385e+03 -0.016 0.9876
-## year2004-01-01 -2.191e+01 1.385e+03 -0.016 0.9874
-## year2005-01-01 -2.227e+01 1.385e+03 -0.016 0.9872
-## year2006-01-01 -2.244e+01 1.385e+03 -0.016 0.9871
-## year2007-01-01 -2.231e+01 1.385e+03 -0.016 0.9872
-## year2008-01-01 -2.264e+01 1.385e+03 -0.016 0.9870
-## year2009-01-01 -2.241e+01 1.385e+03 -0.016 0.9871
-## year2010-01-01 -2.277e+01 1.385e+03 -0.016 0.9869
-## year2011-01-01 -2.287e+01 1.385e+03 -0.017 0.9868
-## year2012-01-01 -2.287e+01 1.385e+03 -0.017 0.9868
-## year2013-01-01 -2.291e+01 1.385e+03 -0.017 0.9868
-## year2014-01-01 -2.310e+01 1.385e+03 -0.017 0.9867
-## year2015-01-01 -2.303e+01 1.385e+03 -0.017 0.9867
-## year2016-01-01 -2.280e+01 1.385e+03 -0.016 0.9869
-## year2017-01-01 -2.321e+01 1.385e+03 -0.017 0.9866
-## year2018-01-01 -2.321e+01 1.385e+03 -0.017 0.9866
-## year2019-01-01 -2.346e+01 1.385e+03 -0.017 0.9865
-## weighted_probs 1.548e-01 9.295e-02 1.665 0.0958 .
-## ---
-## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
-##
-## (Dispersion parameter for binomial family taken to be 1)
-##
-## Null deviance: 5692.3 on 154030 degrees of freedom
-## Residual deviance: 4566.8 on 154003 degrees of freedom
-## AIC: 4622.8
-##
-## Number of Fisher Scoring iterations: 15
type
over and above the main effect of gender probability and year.
+scaled_iscb <- iscb_lm
+ mutate(type = type %>% relevel(ref = "Pubmed authors"))
-scaled_iscb <- iscb_lm %>%
+ filter(year(year) >= 2002)
# scaled_iscb$s_prob <- scale(scaled_iscb$weighted_probs, scale = F)
# scaled_iscb$s_year <- scale(scaled_iscb$year, scale = F)
-main_lm <- glm(type ~ year + weighted_probs,
- data = scaled_iscb %>% mutate(year = as.factor(year)),
- family = 'binomial')
-summary(main_lm)
+main_lm <- glm(type ~ year + weighted_probs,
+ data = scaled_iscb, # %>% mutate(year = as.factor(year))
+ family = "binomial"
+)
+
+broom::tidy(main_lm)##
-## Call:
-## glm(formula = type ~ year + weighted_probs, family = "binomial",
-## data = scaled_iscb %>% mutate(year = as.factor(year)))
-##
-## Deviance Residuals:
-## Min 1Q Median 3Q Max
-## -0.2521 -0.0675 -0.0602 -0.0536 3.7156
-##
-## Coefficients:
-## Estimate Std. Error z value Pr(>|z|)
-## (Intercept) 1.655e+01 1.385e+03 0.012 0.9905
-## year1994-01-01 8.750e-03 1.959e+03 0.000 1.0000
-## year1995-01-01 -3.840e-02 1.959e+03 0.000 1.0000
-## year1996-01-01 5.925e-03 1.833e+03 0.000 1.0000
-## year1997-01-01 -6.766e-03 1.537e+03 0.000 1.0000
-## year1998-01-01 -1.777e-03 1.549e+03 0.000 1.0000
-## year1999-01-01 8.430e-03 1.510e+03 0.000 1.0000
-## year2000-01-01 -1.929e-03 1.510e+03 0.000 1.0000
-## year2001-01-01 -1.254e-03 1.510e+03 0.000 1.0000
-## year2002-01-01 -2.118e+01 1.385e+03 -0.015 0.9878
-## year2003-01-01 -2.150e+01 1.385e+03 -0.016 0.9876
-## year2004-01-01 -2.191e+01 1.385e+03 -0.016 0.9874
-## year2005-01-01 -2.227e+01 1.385e+03 -0.016 0.9872
-## year2006-01-01 -2.244e+01 1.385e+03 -0.016 0.9871
-## year2007-01-01 -2.231e+01 1.385e+03 -0.016 0.9872
-## year2008-01-01 -2.264e+01 1.385e+03 -0.016 0.9870
-## year2009-01-01 -2.241e+01 1.385e+03 -0.016 0.9871
-## year2010-01-01 -2.277e+01 1.385e+03 -0.016 0.9869
-## year2011-01-01 -2.287e+01 1.385e+03 -0.017 0.9868
-## year2012-01-01 -2.287e+01 1.385e+03 -0.017 0.9868
-## year2013-01-01 -2.291e+01 1.385e+03 -0.017 0.9868
-## year2014-01-01 -2.310e+01 1.385e+03 -0.017 0.9867
-## year2015-01-01 -2.303e+01 1.385e+03 -0.017 0.9867
-## year2016-01-01 -2.280e+01 1.385e+03 -0.016 0.9869
-## year2017-01-01 -2.321e+01 1.385e+03 -0.017 0.9866
-## year2018-01-01 -2.321e+01 1.385e+03 -0.017 0.9866
-## year2019-01-01 -2.346e+01 1.385e+03 -0.017 0.9865
-## weighted_probs 1.548e-01 9.295e-02 1.665 0.0958 .
-## ---
-## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
-##
-## (Dispersion parameter for binomial family taken to be 1)
-##
-## Null deviance: 5692.3 on 154030 degrees of freedom
-## Residual deviance: 4566.8 on 154003 degrees of freedom
-## AIC: 4622.8
-##
-## Number of Fisher Scoring iterations: 15
## # A tibble: 3 x 5
+## term estimate std.error statistic p.value
+## <chr> <dbl> <dbl> <dbl> <dbl>
+## 1 (Intercept) -2.06 0.478 -4.30 1.67e- 5
+## 2 year -0.000271 0.0000320 -8.47 2.42e-17
+## 3 weighted_probs 0.155 0.0921 1.69 9.19e- 2
-inte_lm <- glm(
# type ~ scale(year, scale = F) * scale(weighted_probs, scale = F),
- # type ~ s_year * s_prob,
+ # type ~ s_year * s_prob,
type ~ year * weighted_probs,
- data = scaled_iscb %>% mutate(year = as.factor(year))
- ,
- family = 'binomial')
-summary(inte_lm)
-##
-## Call:
-## glm(formula = type ~ year * weighted_probs, family = "binomial",
-## data = scaled_iscb %>% mutate(year = as.factor(year)))
-##
-## Deviance Residuals:
-## Min 1Q Median 3Q Max
-## -0.4974 -0.0666 -0.0600 -0.0544 3.7730
-##
-## Coefficients:
-## Estimate Std. Error z value Pr(>|z|)
-## (Intercept) 1.657e+01 2.486e+03 0.007 0.995
-## year1994-01-01 1.785e-05 5.665e+03 0.000 1.000
-## year1995-01-01 1.785e-05 3.020e+03 0.000 1.000
-## year1996-01-01 1.785e-05 2.976e+03 0.000 1.000
-## year1997-01-01 1.785e-05 2.596e+03 0.000 1.000
-## year1998-01-01 1.785e-05 2.592e+03 0.000 1.000
-## year1999-01-01 1.785e-05 2.646e+03 0.000 1.000
-## year2000-01-01 1.785e-05 2.568e+03 0.000 1.000
-## year2001-01-01 1.785e-05 2.567e+03 0.000 1.000
-## year2002-01-01 -2.101e+01 2.486e+03 -0.008 0.993
-## year2003-01-01 -2.147e+01 2.486e+03 -0.009 0.993
-## year2004-01-01 -2.195e+01 2.486e+03 -0.009 0.993
-## year2005-01-01 -2.227e+01 2.486e+03 -0.009 0.993
-## year2006-01-01 -2.237e+01 2.486e+03 -0.009 0.993
-## year2007-01-01 -2.244e+01 2.486e+03 -0.009 0.993
-## year2008-01-01 -2.268e+01 2.486e+03 -0.009 0.993
-## year2009-01-01 -2.241e+01 2.486e+03 -0.009 0.993
-## year2010-01-01 -2.274e+01 2.486e+03 -0.009 0.993
-## year2011-01-01 -2.295e+01 2.486e+03 -0.009 0.993
-## year2012-01-01 -2.288e+01 2.486e+03 -0.009 0.993
-## year2013-01-01 -2.288e+01 2.486e+03 -0.009 0.993
-## year2014-01-01 -2.294e+01 2.486e+03 -0.009 0.993
-## year2015-01-01 -2.291e+01 2.486e+03 -0.009 0.993
-## year2016-01-01 -2.288e+01 2.486e+03 -0.009 0.993
-## year2017-01-01 -2.307e+01 2.486e+03 -0.009 0.993
-## year2018-01-01 -2.332e+01 2.486e+03 -0.009 0.993
-## year2019-01-01 -2.352e+01 2.486e+03 -0.009 0.992
-## weighted_probs 2.470e-04 2.815e+04 0.000 1.000
-## year1994-01-01:weighted_probs -2.470e-04 2.952e+05 0.000 1.000
-## year1995-01-01:weighted_probs -2.470e-04 2.831e+04 0.000 1.000
-## year1996-01-01:weighted_probs -2.470e-04 4.246e+04 0.000 1.000
-## year1997-01-01:weighted_probs -2.470e-04 2.829e+04 0.000 1.000
-## year1998-01-01:weighted_probs -2.470e-04 2.828e+04 0.000 1.000
-## year1999-01-01:weighted_probs -2.470e-04 4.592e+04 0.000 1.000
-## year2000-01-01:weighted_probs -2.470e-04 2.827e+04 0.000 1.000
-## year2001-01-01:weighted_probs -2.470e-04 2.827e+04 0.000 1.000
-## year2002-01-01:weighted_probs -1.987e+00 2.815e+04 0.000 1.000
-## year2003-01-01:weighted_probs -1.160e-01 2.815e+04 0.000 1.000
-## year2004-01-01:weighted_probs 2.712e-01 2.815e+04 0.000 1.000
-## year2005-01-01:weighted_probs 8.687e-02 2.815e+04 0.000 1.000
-## year2006-01-01:weighted_probs -2.989e-01 2.815e+04 0.000 1.000
-## year2007-01-01:weighted_probs 6.183e-01 2.815e+04 0.000 1.000
-## year2008-01-01:weighted_probs 2.592e-01 2.815e+04 0.000 1.000
-## year2009-01-01:weighted_probs 1.071e-01 2.815e+04 0.000 1.000
-## year2010-01-01:weighted_probs -5.974e-02 2.815e+04 0.000 1.000
-## year2011-01-01:weighted_probs 3.763e-01 2.815e+04 0.000 1.000
-## year2012-01-01:weighted_probs 1.388e-01 2.815e+04 0.000 1.000
-## year2013-01-01:weighted_probs 2.515e-02 2.815e+04 0.000 1.000
-## year2014-01-01:weighted_probs -7.596e-01 2.815e+04 0.000 1.000
-## year2015-01-01:weighted_probs -4.676e-01 2.815e+04 0.000 1.000
-## year2016-01-01:weighted_probs 3.395e-01 2.815e+04 0.000 1.000
-## year2017-01-01:weighted_probs -5.012e-01 2.815e+04 0.000 1.000
-## year2018-01-01:weighted_probs 4.535e-01 2.815e+04 0.000 1.000
-## year2019-01-01:weighted_probs 3.118e-01 2.815e+04 0.000 1.000
-##
-## (Dispersion parameter for binomial family taken to be 1)
-##
-## Null deviance: 5692.3 on 154030 degrees of freedom
-## Residual deviance: 4555.1 on 153977 degrees of freedom
-## AIC: 4663.1
-##
-## Number of Fisher Scoring iterations: 15
+ data = scaled_iscb, # %>% mutate(year = as.factor(year))
+ family = "binomial"
+)
+broom::tidy(inte_lm)
+anova(main_lm, inte_lm, test = 'Chisq')
+## # A tibble: 4 x 5
+## term estimate std.error statistic p.value
+## <chr> <dbl> <dbl> <dbl> <dbl>
+## 1 (Intercept) -1.91 0.523 -3.65 2.59e- 4
+## 2 year -0.000281 0.0000351 -7.99 1.30e-15
+## 3 weighted_probs -0.445 0.901 -0.494 6.21e- 1
+## 4 year:weighted_probs 0.0000402 0.0000593 0.679 4.97e- 1
anova(main_lm, inte_lm, test = "Chisq")
-## Analysis of Deviance Table
##
## Model 1: type ~ year + weighted_probs
## Model 2: type ~ year * weighted_probs
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
-## 1 154003 4566.8
-## 2 153977 4555.1 26 11.719 0.9926
-mean(scaled_iscb$year)
-## Warning in mean.default(scaled_iscb$year): argument is not numeric or logical: returning NA
-## [1] NA
-mean(scaled_iscb$weighted_probs)
-## [1] 0.2260815
+## 1 153942 4582.3
+## 2 153941 4581.8 1 0.47034 0.4928
+# inte_lm <- glm(type ~ (year * weighted_probs),
-# data = iscb_lm,
-# family = 'binomial')
+# mean(scaled_iscb$year)
+# mean(scaled_iscb$weighted_probs)
type
over and above the main effect of gender probability and year.sessionInfo()
+## [1] colorspace_2.0-0 ellipsis_0.3.1 class_7.3-17
+## [4] rprojroot_1.3-2 fs_1.5.0 rstudioapi_0.12
+## [7] farver_2.0.3 remotes_2.2.0 DT_0.16
+## [10] prodlim_2019.11.13 fansi_0.4.1 xml2_1.3.2
+## [13] codetools_0.2-16 splines_4.0.3 knitr_1.30
+## [16] pkgload_1.1.0 jsonlite_1.7.1 pROC_1.16.2
+## [19] broom_0.7.2 dbplyr_2.0.0 rgeos_0.5-5
+## [22] compiler_4.0.3 httr_1.4.2 backports_1.2.0
+## [25] assertthat_0.2.1 Matrix_1.2-18 cli_2.1.0
+## [28] htmltools_0.5.0 prettyunits_1.1.1 tools_4.0.3
+## [31] gtable_0.3.0 glue_1.4.2 rnaturalearthdata_0.1.0
+## [34] reshape2_1.4.4 Rcpp_1.0.5 cellranger_1.1.0
+## [37] vctrs_0.3.4 svglite_1.2.3.2 nlme_3.1-149
+## [40] iterators_1.0.13 crosstalk_1.1.0.1 timeDate_3043.102
+## [43] gower_0.2.2 xfun_0.19 ps_1.4.0
+## [46] testthat_3.0.0 rvest_0.3.6 lifecycle_0.2.0
+## [49] devtools_2.3.2 MASS_7.3-53 scales_1.1.1
+## [52] ipred_0.9-9 hms_0.5.3 RColorBrewer_1.1-2
+## [55] yaml_2.2.1 curl_4.3 memoise_1.1.0
+## [58] rpart_4.1-15 stringi_1.5.3 desc_1.2.0
+## [61] foreach_1.5.1 e1071_1.7-4 pkgbuild_1.1.0
+## [64] lava_1.6.8.1 systemfonts_0.3.2 rlang_0.4.8
+## [67] pkgconfig_2.0.3 evaluate_0.14 sf_0.9-6
+## [70] recipes_0.1.15 htmlwidgets_1.5.2 labeling_0.4.2
+## [73] cowplot_1.1.0 tidyselect_1.1.0 processx_3.4.4
+## [76] plyr_1.8.6 magrittr_1.5 R6_2.5.0
+## [79] generics_0.1.0 DBI_1.1.0 mgcv_1.8-33
+## [82] pillar_1.4.6 haven_2.3.1 withr_2.3.0
+## [85] units_0.6-7 survival_3.2-7 sp_1.4-4
+## [88] nnet_7.3-14 modelr_0.1.8 crayon_1.3.4
+## [91] KernSmooth_2.23-17 utf8_1.1.4 rmarkdown_2.5
+## [94] usethis_1.6.3 grid_4.0.3 readxl_1.3.1
+## [97] data.table_1.13.2 callr_3.5.1 ModelMetrics_1.2.2.2
+## [100] reprex_0.3.0 digest_0.6.27 classInt_0.4-3
+## [103] stats4_4.0.3 munsell_0.5.0 viridisLite_0.3.0
+## [106] sessioninfo_1.1.1
## R version 4.0.3 (2020-10-10)
## Platform: x86_64-pc-linux-gnu (64-bit)
@@ -2024,51 +1870,63 @@
Hypothesis testing
## BLAS/LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.8.so
##
## locale:
-## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=en_US.UTF-8
-## [4] LC_COLLATE=en_US.UTF-8 LC_MONETARY=en_US.UTF-8 LC_MESSAGES=C
-## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C LC_ADDRESS=C
-## [10] LC_TELEPHONE=C LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
+## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
+## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
+## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=C
+## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
+## [9] LC_ADDRESS=C LC_TELEPHONE=C
+## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
-## [1] gdtools_0.2.2 wru_0.1-10 rnaturalearth_0.1.0 lubridate_1.7.9.2
-## [5] caret_6.0-86 lattice_0.20-41 forcats_0.5.0 stringr_1.4.0
-## [9] dplyr_1.0.2 purrr_0.3.4 readr_1.4.0 tidyr_1.1.2
+## [1] gdtools_0.2.2 wru_0.1-10 rnaturalearth_0.1.0
+## [4] lubridate_1.7.9.2 caret_6.0-86 lattice_0.20-41
+## [7] forcats_0.5.0 stringr_1.4.0 dplyr_1.0.2
+## [10] purrr_0.3.4 readr_1.4.0 tidyr_1.1.2
## [13] tibble_3.0.4 ggplot2_3.3.2 tidyverse_1.3.0
##
## loaded via a namespace (and not attached):
-## [1] colorspace_2.0-0 ellipsis_0.3.1 class_7.3-17 rprojroot_1.3-2
-## [5] fs_1.5.0 rstudioapi_0.12 farver_2.0.3 remotes_2.2.0
-## [9] DT_0.16 prodlim_2019.11.13 fansi_0.4.1 xml2_1.3.2
-## [13] codetools_0.2-16 splines_4.0.3 knitr_1.30 pkgload_1.1.0
-## [17] jsonlite_1.7.1 pROC_1.16.2 broom_0.7.2 dbplyr_2.0.0
-## [21] rgeos_0.5-5 compiler_4.0.3 httr_1.4.2 backports_1.2.0
-## [25] assertthat_0.2.1 Matrix_1.2-18 cli_2.1.0 htmltools_0.5.0
-## [29] prettyunits_1.1.1 tools_4.0.3 gtable_0.3.0 glue_1.4.2
-## [33] rnaturalearthdata_0.1.0 reshape2_1.4.4 Rcpp_1.0.5 cellranger_1.1.0
-## [37] vctrs_0.3.4 svglite_1.2.3.2 nlme_3.1-149 iterators_1.0.13
-## [41] crosstalk_1.1.0.1 timeDate_3043.102 gower_0.2.2 xfun_0.19
-## [45] ps_1.4.0 testthat_3.0.0 rvest_0.3.6 lifecycle_0.2.0
-## [49] devtools_2.3.2 MASS_7.3-53 scales_1.1.1 ipred_0.9-9
-## [53] hms_0.5.3 RColorBrewer_1.1-2 yaml_2.2.1 curl_4.3
-## [57] memoise_1.1.0 rpart_4.1-15 stringi_1.5.3 desc_1.2.0
-## [61] foreach_1.5.1 e1071_1.7-4 pkgbuild_1.1.0 lava_1.6.8.1
-## [65] systemfonts_0.3.2 rlang_0.4.8 pkgconfig_2.0.3 evaluate_0.14
-## [69] sf_0.9-6 recipes_0.1.15 htmlwidgets_1.5.2 labeling_0.4.2
-## [73] cowplot_1.1.0 tidyselect_1.1.0 processx_3.4.4 plyr_1.8.6
-## [77] magrittr_1.5 R6_2.5.0 generics_0.1.0 DBI_1.1.0
-## [81] mgcv_1.8-33 pillar_1.4.6 haven_2.3.1 withr_2.3.0
-## [85] units_0.6-7 survival_3.2-7 sp_1.4-4 nnet_7.3-14
-## [89] modelr_0.1.8 crayon_1.3.4 KernSmooth_2.23-17 utf8_1.1.4
-## [93] rmarkdown_2.5 usethis_1.6.3 grid_4.0.3 readxl_1.3.1
-## [97] data.table_1.13.2 callr_3.5.1 ModelMetrics_1.2.2.2 reprex_0.3.0
-## [101] digest_0.6.27 classInt_0.4-3 stats4_4.0.3 munsell_0.5.0
-## [105] viridisLite_0.3.0 sessioninfo_1.1.1Representation analysis of name origins
+region_cols <- c("#ffffb3", "#fccde5", "#b3de69", "#fdb462", "#80b1d3", "#8dd3c7", "#bebada", "#fb8072", "#bc80bd", "#ccebc5")
alpha_threshold <- qnorm(0.975)
-load('Rdata/raws.Rdata')
+load("Rdata/raws.Rdata")
pubmed_nat_df <- corr_authors %>%
- filter(year(year) >= 2002) %>%
- left_join(nationalize_df, by = c('fore_name', 'last_name')) %>%
+ filter(year(year) >= 2002) %>%
+ left_join(nationalize_df, by = c("fore_name", "last_name")) %>%
group_by(pmid, journal, publication_date, year, adjusted_citations) %>%
summarise_at(vars(African:SouthAsian), mean, na.rm = T) %>%
ungroup()
iscb_nat_df <- keynotes %>%
- left_join(nationalize_df, by = c('fore_name', 'last_name'))
+ left_join(nationalize_df, by = c("fore_name", "last_name"))
start_year <- 1992
end_year <- 2019
@@ -1686,36 +1686,40 @@
Representation analysis of name origins
my_confs <- unique(iscb_nat_df$conference)
n_jours <- length(my_jours)
n_confs <- length(my_confs)
-region_levels <- paste(c('Celtic/English', 'European', 'East Asian', 'Hispanic', 'South Asian', 'Arabic', 'Hebrew', 'African', 'Nordic', 'Greek'), 'names')
+region_levels <- paste(c("Celtic/English", "European", "East Asian", "Hispanic", "South Asian", "Arabic", "Hebrew", "African", "Nordic", "Greek"), "names")
region_levels_let <- toupper(letters[1:8])
-region_cols <- c('#ffffb3', '#fccde5', '#b3de69', '#fdb462', '#80b1d3', '#8dd3c7', '#bebada', '#fb8072', '#bc80bd', '#ccebc5')
-our_country_map <- read_tsv('https://raw.githubusercontent.com/greenelab/wiki-nationality-estimate/7c22d0a5f661ce5aeb785215095deda40973ff17/data/country_to_region_NamePrism.tsv') %>%
- rename('region' = Region) %>%
+
our_country_map <- read_tsv("https://raw.githubusercontent.com/greenelab/wiki-nationality-estimate/7c22d0a5f661ce5aeb785215095deda40973ff17/data/country_to_region_NamePrism.tsv") %>%
+ rename("region" = Region) %>%
recode_region()
##
-## ── Column specification ───────────────────────────────────────────────────────────────────────────────
+## ── Column specification ────────────────────────────────────────────────────────
## cols(
## Country = col_character(),
## Region = col_character()
## )
+ geom_sf(aes(fill = fct_relevel(region, region_levels))) +
+ coord_sf(crs = "+proj=eqearth +wktext") +
+ scale_fill_manual(
+ values = region_cols,
+ na.translate = FALSE
+ ) +
+ theme(
+ panel.background = element_rect(fill = "azure"),
+ legend.title = element_blank(),
+ legend.position = "bottom",
+ panel.border = element_rect(fill = NA)
+ ))my_world <- world %>%
- select(- geometry) %>%
- rename(Country = 'name') %>%
- left_join(our_country_map, by = 'Country')
+ select(-geometry) %>%
+ rename(Country = "name") %>%
+ left_join(our_country_map, by = "Country")
(gworld <- ggplot(data = my_world) +
- geom_sf(aes(fill = fct_relevel(region, region_levels))) +
- coord_sf(crs = "+proj=eqearth +wktext") +
- scale_fill_manual(values = region_cols,
- na.translate = FALSE) +
- theme(panel.background = element_rect(fill = "azure"),
- legend.title = element_blank(),
- legend.position = 'bottom',
- panel.border = element_rect(fill = NA)))
+ggsave('figs/2020-01-31_groupings.png', gworld, width = 7.2, height = 4.3)
-ggsave('figs/2020-01-31_groupings.svg', gworld, width = 7.2, height = 4.3)
ggsave("figs/2020-01-31_groupings.png", gworld, width = 7.2, height = 4.3)
+ggsave("figs/2020-01-31_groupings.svg", gworld, width = 7.2, height = 4.3)
Descriptive statistics
Descriptive statistics
+ filter(region != "OtherCategories")
iscb_pubmed_oth <- iscb_nat_df %>%
- rename('journal' = conference) %>%
+ rename("journal" = conference) %>%
select(year, journal, African:SouthAsian, publication_date) %>%
- mutate(type = 'Keynote speakers/Fellows',
- adjusted_citations = 1,
- pmid = -9999) %>%
+ mutate(
+ type = "Keynote speakers/Fellows",
+ adjusted_citations = 1,
+ pmid = -9999
+ ) %>%
bind_rows(
pubmed_nat_df %>%
select(pmid, year, journal, African:SouthAsian, publication_date, adjusted_citations) %>%
- mutate(type = 'Pubmed authors')
+ mutate(type = "Pubmed authors")
) %>%
mutate(OtherCategories = SouthAsian + Hispanic + Jewish + Muslim + Nordic + Greek + African) %>%
pivot_longer(c(African:SouthAsian, OtherCategories),
- names_to = 'region',
- values_to = 'probabilities') %>%
- filter(!is.na(probabilities)) %>%
+ names_to = "region",
+ values_to = "probabilities"
+ ) %>%
+ filter(!is.na(probabilities)) %>%
group_by(type, year, region) %>%
mutate(
pmc_citations_year = mean(adjusted_citations),
- weight = adjusted_citations/pmc_citations_year,
- weighted_probs = probabilities*weight
- )
+ weight = adjusted_citations / pmc_citations_year,
+ weighted_probs = probabilities * weight
+ )
-iscb_pubmed_sum_oth <- iscb_pubmed_oth %>%
+iscb_pubmed_sum_oth <- iscb_pubmed_oth %>%
summarise(
mean_prob = mean(weighted_probs),
- se_prob = sqrt(var(probabilities) * sum(weight^2)/(sum(weight)^2)),
+ se_prob = sqrt(var(probabilities) * sum(weight^2) / (sum(weight)^2)),
me_prob = alpha_threshold * se_prob,
- .groups = 'drop'
+ .groups = "drop"
)
iscb_pubmed_sum <- iscb_pubmed_sum_oth %>%
- filter(region != 'OtherCategories')
Prepare data frames for analysis
By conference keynotes/fellows
-i <- 0
-iscb_nat <- vector('list', length = n_confs)
+iscb_nat <- vector("list", length = n_confs)
-for (conf in my_confs){
+for (conf in my_confs) {
i <- i + 1
iscb_nat[[i]] <- iscb_pubmed_oth %>%
- filter(region != 'OtherCategories', type != 'Pubmed authors' & journal == conf) %>%
+ filter(region != "OtherCategories", type != "Pubmed authors" & journal == conf) %>%
group_by(type, year, region, journal) %>%
- summarise(mean_prob = mean(weighted_probs), .groups = 'drop') %>%
- {.}
+ summarise(mean_prob = mean(weighted_probs), .groups = "drop")
}
+save(my_world, iscb_pubmed_oth, iscb_nat, file = 'Rdata/iscb-pubmed_nat.Rdata')
save(my_world, iscb_pubmed_oth, iscb_nat, file = "Rdata/iscb-pubmed_nat.Rdata")
Figures for paper
Figure 4
+fig_4 <- cowplot::plot_grid(fig_4a, fig_4b, labels = "AUTO", ncol = 1, rel_heights = c(1.3, 1))
fig_4a <- iscb_pubmed_sum %>%
- filter(year < '2020-01-01') %>%
- region_breakdown('main', region_levels, fct_rev(type)) +
+ filter(year < "2020-01-01") %>%
+ region_breakdown("main", region_levels, fct_rev(type)) +
guides(fill = guide_legend(nrow = 2)) +
- theme(legend.position = 'bottom')
+ theme(legend.position = "bottom")
-large_regions <- c('CelticEnglish', 'EastAsian', 'European', 'OtherCategories')
+large_regions <- c("CelticEnglish", "EastAsian", "European", "OtherCategories")
## Mean and standard deviation of predicted probabilities:
fig_4b <- iscb_pubmed_sum_oth %>%
filter(region %in% large_regions) %>%
recode_region() %>%
gam_and_ci(
- df2 = iscb_pubmed_oth %>%
+ df2 = iscb_pubmed_oth %>%
filter(region %in% large_regions) %>%
recode_region(),
- start_y = start_year, end_y = end_year) +
- theme(legend.position = c(0.88, 0.83),
- panel.grid.minor = element_blank(),
- legend.margin = margin(-0.5, 0, 0, 0, unit='cm'),
- legend.text = element_text(size = 7)) +
+ start_y = start_year, end_y = end_year
+ ) +
+ theme(
+ legend.position = c(0.88, 0.83),
+ panel.grid.minor = element_blank(),
+ legend.margin = margin(-0.5, 0, 0, 0, unit = "cm"),
+ legend.text = element_text(size = 7)
+ ) +
facet_wrap(vars(fct_relevel(region, large_regions)), nrow = 1)
-fig_4 <- cowplot::plot_grid(fig_4a, fig_4b, labels = 'AUTO', ncol = 1, rel_heights = c(1.3,1))
## `geom_smooth()` using formula 'y ~ s(x, bs = "cs")'
-fig_4
+ggsave('figs/region_breakdown.png', fig_4, width = 6.7, height = 5.5)
-ggsave('figs/region_breakdown.svg', fig_4, width = 6.7, height = 5.5)
ggsave("figs/region_breakdown.png", fig_4, width = 6.7, height = 5.5)
+ggsave("figs/region_breakdown.svg", fig_4, width = 6.7, height = 5.5)
Hypothesis testing
-iscb_lm <- iscb_pubmed_oth %>%
- ungroup() %>%
+
iscb_lm <- iscb_pubmed_oth %>%
+ ungroup() %>%
mutate(
# year = c(scale(year(year))),
- year = as.factor(year),
- type = as.factor(type) %>% relevel(ref = 'Pubmed authors'))
-main_lm <- function(regioni){
- glm(type ~ year + weighted_probs,
- data = iscb_lm %>%
- filter(region == regioni, !is.na(probabilities)) ,
- family = 'binomial')
+ # year = as.factor(year),
+ type = as.factor(type) %>% relevel(ref = "Pubmed authors")
+ )
+main_lm <- function(regioni) {
+ glm(type ~ year + weighted_probs,
+ data = iscb_lm %>%
+ filter(region == regioni, !is.na(probabilities), year(year) >= 2002),
+ family = "binomial"
+ )
}
-inte_lm <- function(regioni){
- glm(type ~ year * weighted_probs,
- data = iscb_lm %>%
- filter(region == regioni, !is.na(weighted_probs)),
- family = 'binomial')
+inte_lm <- function(regioni) {
+ glm(type ~ year * weighted_probs,
+ data = iscb_lm %>%
+ filter(region == regioni, !is.na(weighted_probs), year(year) >= 2002),
+ family = "binomial"
+ )
}
main_list <- lapply(large_regions, main_lm)
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+lapply(main_list, broom::tidy)names(main_list) <- large_regions
-lapply(main_list, summary)
+## # A tibble: 3 x 5
+## term estimate std.error statistic p.value
+## <chr> <dbl> <dbl> <dbl> <dbl>
+## 1 (Intercept) -2.06 0.479 -4.30 1.73e- 5
+## 2 year -0.000273 0.0000319 -8.57 1.03e-17
+## 3 weighted_probs 0.0724 0.0948 0.763 4.45e- 1
## $CelticEnglish
-##
-## Call:
-## glm(formula = type ~ year + weighted_probs, family = "binomial",
-## data = iscb_lm %>% filter(region == regioni, !is.na(probabilities)))
-##
-## Deviance Residuals:
-## Min 1Q Median 3Q Max
-## -0.6215 -0.0651 -0.0580 -0.0513 3.7315
-##
-## Coefficients:
-## Estimate Std. Error z value Pr(>|z|)
-## (Intercept) 1.648e+01 1.385e+03 0.012 0.990508
-## year1994-01-01 -6.727e-02 1.959e+03 0.000 0.999973
-## year1995-01-01 -9.883e-02 1.959e+03 0.000 0.999960
-## year1996-01-01 -1.957e-02 1.832e+03 0.000 0.999991
-## year1997-01-01 7.983e-03 1.536e+03 0.000 0.999996
-## year1998-01-01 2.415e-02 1.548e+03 0.000 0.999988
-## year1999-01-01 1.437e-02 1.509e+03 0.000 0.999992
-## year2000-01-01 3.791e-02 1.509e+03 0.000 0.999980
-## year2001-01-01 1.267e-02 1.509e+03 0.000 0.999993
-## year2002-01-01 -2.118e+01 1.385e+03 -0.015 0.987800
-## year2003-01-01 -2.150e+01 1.385e+03 -0.016 0.987612
-## year2004-01-01 -2.191e+01 1.385e+03 -0.016 0.987380
-## year2005-01-01 -2.227e+01 1.385e+03 -0.016 0.987172
-## year2006-01-01 -2.244e+01 1.385e+03 -0.016 0.987073
-## year2007-01-01 -2.231e+01 1.385e+03 -0.016 0.987147
-## year2008-01-01 -2.263e+01 1.385e+03 -0.016 0.986961
-## year2009-01-01 -2.240e+01 1.385e+03 -0.016 0.987095
-## year2010-01-01 -2.277e+01 1.385e+03 -0.016 0.986881
-## year2011-01-01 -2.287e+01 1.385e+03 -0.017 0.986825
-## year2012-01-01 -2.287e+01 1.385e+03 -0.017 0.986825
-## year2013-01-01 -2.285e+01 1.385e+03 -0.017 0.986835
-## year2014-01-01 -2.310e+01 1.385e+03 -0.017 0.986692
-## year2015-01-01 -2.303e+01 1.385e+03 -0.017 0.986735
-## year2016-01-01 -2.279e+01 1.385e+03 -0.016 0.986870
-## year2017-01-01 -2.320e+01 1.385e+03 -0.017 0.986635
-## year2018-01-01 -2.319e+01 1.385e+03 -0.017 0.986640
-## year2019-01-01 -2.344e+01 1.385e+03 -0.017 0.986498
-## weighted_probs 1.921e-01 5.650e-02 3.400 0.000673 ***
-## ---
-## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
-##
-## (Dispersion parameter for binomial family taken to be 1)
-##
-## Null deviance: 5755.8 on 163974 degrees of freedom
-## Residual deviance: 4611.5 on 163947 degrees of freedom
-## AIC: 4667.5
-##
-## Number of Fisher Scoring iterations: 15
-##
+## # A tibble: 3 x 5
+## term estimate std.error statistic p.value
+## <chr> <dbl> <dbl> <dbl> <dbl>
+## 1 (Intercept) -2.17 0.481 -4.52 6.06e- 6
+## 2 year -0.000268 0.0000320 -8.37 5.64e-17
+## 3 weighted_probs 0.194 0.0561 3.46 5.46e- 4
##
## $EastAsian
-##
-## Call:
-## glm(formula = type ~ year + weighted_probs, family = "binomial",
-## data = iscb_lm %>% filter(region == regioni, !is.na(probabilities)))
-##
-## Deviance Residuals:
-## Min 1Q Median 3Q Max
-## -0.1463 -0.0682 -0.0594 -0.0507 4.0267
-##
-## Coefficients:
-## Estimate Std. Error z value Pr(>|z|)
-## (Intercept) 1.658e+01 1.385e+03 0.012 0.990
-## year1994-01-01 6.665e-02 1.959e+03 0.000 1.000
-## year1995-01-01 -1.257e-02 1.959e+03 0.000 1.000
-## year1996-01-01 -2.152e-03 1.833e+03 0.000 1.000
-## year1997-01-01 6.364e-03 1.537e+03 0.000 1.000
-## year1998-01-01 -4.470e-03 1.549e+03 0.000 1.000
-## year1999-01-01 7.432e-03 1.510e+03 0.000 1.000
-## year2000-01-01 6.429e-01 1.502e+03 0.000 1.000
-## year2001-01-01 -5.036e-03 1.510e+03 0.000 1.000
-## year2002-01-01 -2.111e+01 1.385e+03 -0.015 0.988
-## year2003-01-01 -2.142e+01 1.385e+03 -0.015 0.988
-## year2004-01-01 -2.181e+01 1.385e+03 -0.016 0.987
-## year2005-01-01 -2.217e+01 1.385e+03 -0.016 0.987
-## year2006-01-01 -2.233e+01 1.385e+03 -0.016 0.987
-## year2007-01-01 -2.221e+01 1.385e+03 -0.016 0.987
-## year2008-01-01 -2.252e+01 1.385e+03 -0.016 0.987
-## year2009-01-01 -2.228e+01 1.385e+03 -0.016 0.987
-## year2010-01-01 -2.264e+01 1.385e+03 -0.016 0.987
-## year2011-01-01 -2.273e+01 1.385e+03 -0.016 0.987
-## year2012-01-01 -2.272e+01 1.385e+03 -0.016 0.987
-## year2013-01-01 -2.268e+01 1.385e+03 -0.016 0.987
-## year2014-01-01 -2.292e+01 1.385e+03 -0.017 0.987
-## year2015-01-01 -2.285e+01 1.385e+03 -0.016 0.987
-## year2016-01-01 -2.260e+01 1.385e+03 -0.016 0.987
-## year2017-01-01 -2.299e+01 1.385e+03 -0.017 0.987
-## year2018-01-01 -2.297e+01 1.385e+03 -0.017 0.987
-## year2019-01-01 -2.320e+01 1.385e+03 -0.017 0.987
-## weighted_probs -1.665e+00 2.861e-01 -5.820 5.89e-09 ***
-## ---
-## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
-##
-## (Dispersion parameter for binomial family taken to be 1)
-##
-## Null deviance: 5755.8 on 163974 degrees of freedom
-## Residual deviance: 4560.8 on 163947 degrees of freedom
-## AIC: 4616.8
-##
-## Number of Fisher Scoring iterations: 15
-##
+## # A tibble: 3 x 5
+## term estimate std.error statistic p.value
+## <chr> <dbl> <dbl> <dbl> <dbl>
+## 1 (Intercept) -2.30 0.480 -4.80 1.58e- 6
+## 2 year -0.000242 0.0000321 -7.54 4.82e-14
+## 3 weighted_probs -1.67 0.286 -5.82 6.02e- 9
##
## $European
-##
-## Call:
-## glm(formula = type ~ year + weighted_probs, family = "binomial",
-## data = iscb_lm %>% filter(region == regioni, !is.na(probabilities)))
-##
-## Deviance Residuals:
-## Min 1Q Median 3Q Max
-## -0.1960 -0.0638 -0.0588 -0.0509 3.7279
-##
-## Coefficients:
-## Estimate Std. Error z value Pr(>|z|)
-## (Intercept) 1.655e+01 1.385e+03 0.012 0.990
-## year1994-01-01 1.671e-02 1.959e+03 0.000 1.000
-## year1995-01-01 1.751e-02 1.959e+03 0.000 1.000
-## year1996-01-01 -5.959e-03 1.833e+03 0.000 1.000
-## year1997-01-01 -1.572e-03 1.537e+03 0.000 1.000
-## year1998-01-01 -4.096e-03 1.549e+03 0.000 1.000
-## year1999-01-01 -5.232e-03 1.510e+03 0.000 1.000
-## year2000-01-01 -1.138e-03 1.510e+03 0.000 1.000
-## year2001-01-01 -1.393e-02 1.510e+03 0.000 1.000
-## year2002-01-01 -2.120e+01 1.385e+03 -0.015 0.988
-## year2003-01-01 -2.153e+01 1.385e+03 -0.016 0.988
-## year2004-01-01 -2.194e+01 1.385e+03 -0.016 0.987
-## year2005-01-01 -2.230e+01 1.385e+03 -0.016 0.987
-## year2006-01-01 -2.248e+01 1.385e+03 -0.016 0.987
-## year2007-01-01 -2.235e+01 1.385e+03 -0.016 0.987
-## year2008-01-01 -2.267e+01 1.385e+03 -0.016 0.987
-## year2009-01-01 -2.244e+01 1.385e+03 -0.016 0.987
-## year2010-01-01 -2.281e+01 1.385e+03 -0.016 0.987
-## year2011-01-01 -2.291e+01 1.385e+03 -0.017 0.987
-## year2012-01-01 -2.291e+01 1.385e+03 -0.017 0.987
-## year2013-01-01 -2.290e+01 1.385e+03 -0.017 0.987
-## year2014-01-01 -2.315e+01 1.385e+03 -0.017 0.987
-## year2015-01-01 -2.307e+01 1.385e+03 -0.017 0.987
-## year2016-01-01 -2.284e+01 1.385e+03 -0.016 0.987
-## year2017-01-01 -2.325e+01 1.385e+03 -0.017 0.987
-## year2018-01-01 -2.325e+01 1.385e+03 -0.017 0.987
-## year2019-01-01 -2.350e+01 1.385e+03 -0.017 0.986
-## weighted_probs 7.241e-02 8.237e-02 0.879 0.379
-##
-## (Dispersion parameter for binomial family taken to be 1)
-##
-## Null deviance: 5755.8 on 163974 degrees of freedom
-## Residual deviance: 4618.3 on 163947 degrees of freedom
-## AIC: 4674.3
-##
-## Number of Fisher Scoring iterations: 15
-##
+## # A tibble: 3 x 5
+## term estimate std.error statistic p.value
+## <chr> <dbl> <dbl> <dbl> <dbl>
+## 1 (Intercept) -2.08 0.480 -4.32 1.53e- 5
+## 2 year -0.000272 0.0000319 -8.53 1.46e-17
+## 3 weighted_probs 0.0713 0.0822 0.867 3.86e- 1
##
## $OtherCategories
-##
-## Call:
-## glm(formula = type ~ year + weighted_probs, family = "binomial",
-## data = iscb_lm %>% filter(region == regioni, !is.na(probabilities)))
-##
-## Deviance Residuals:
-## Min 1Q Median 3Q Max
-## -0.1628 -0.0636 -0.0589 -0.0508 3.7283
-##
-## Coefficients:
-## Estimate Std. Error z value Pr(>|z|)
-## (Intercept) 1.655e+01 1.385e+03 0.012 0.990
-## year1994-01-01 1.077e-02 1.959e+03 0.000 1.000
-## year1995-01-01 1.868e-02 1.959e+03 0.000 1.000
-## year1996-01-01 1.373e-02 1.833e+03 0.000 1.000
-## year1997-01-01 -1.182e-03 1.537e+03 0.000 1.000
-## year1998-01-01 -5.156e-03 1.549e+03 0.000 1.000
-## year1999-01-01 3.844e-04 1.510e+03 0.000 1.000
-## year2000-01-01 4.700e-03 1.510e+03 0.000 1.000
-## year2001-01-01 9.411e-03 1.510e+03 0.000 1.000
-## year2002-01-01 -2.119e+01 1.385e+03 -0.015 0.988
-## year2003-01-01 -2.152e+01 1.385e+03 -0.016 0.988
-## year2004-01-01 -2.193e+01 1.385e+03 -0.016 0.987
-## year2005-01-01 -2.230e+01 1.385e+03 -0.016 0.987
-## year2006-01-01 -2.247e+01 1.385e+03 -0.016 0.987
-## year2007-01-01 -2.234e+01 1.385e+03 -0.016 0.987
-## year2008-01-01 -2.267e+01 1.385e+03 -0.016 0.987
-## year2009-01-01 -2.244e+01 1.385e+03 -0.016 0.987
-## year2010-01-01 -2.280e+01 1.385e+03 -0.016 0.987
-## year2011-01-01 -2.291e+01 1.385e+03 -0.017 0.987
-## year2012-01-01 -2.291e+01 1.385e+03 -0.017 0.987
-## year2013-01-01 -2.290e+01 1.385e+03 -0.017 0.987
-## year2014-01-01 -2.314e+01 1.385e+03 -0.017 0.987
-## year2015-01-01 -2.307e+01 1.385e+03 -0.017 0.987
-## year2016-01-01 -2.284e+01 1.385e+03 -0.016 0.987
-## year2017-01-01 -2.325e+01 1.385e+03 -0.017 0.987
-## year2018-01-01 -2.325e+01 1.385e+03 -0.017 0.987
-## year2019-01-01 -2.350e+01 1.385e+03 -0.017 0.986
-## weighted_probs 7.448e-02 9.548e-02 0.780 0.435
-##
-## (Dispersion parameter for binomial family taken to be 1)
-##
-## Null deviance: 5755.8 on 163974 degrees of freedom
-## Residual deviance: 4618.4 on 163947 degrees of freedom
-## AIC: 4674.4
-##
-## Number of Fisher Scoring iterations: 15
inte_list <- lapply(large_regions, inte_lm)
-## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
+lapply(inte_list, summary)
lapply(inte_list, broom::tidy)
-## [[1]]
-##
-## Call:
-## glm(formula = type ~ year * weighted_probs, family = "binomial",
-## data = iscb_lm %>% filter(region == regioni, !is.na(weighted_probs)))
-##
-## Deviance Residuals:
-## Min 1Q Median 3Q Max
-## -0.5406 -0.0661 -0.0589 -0.0514 3.7625
-##
-## Coefficients:
-## Estimate Std. Error z value Pr(>|z|)
-## (Intercept) 1.657e+01 2.261e+03 0.007 0.994
-## year1994-01-01 1.551e-07 1.105e+04 0.000 1.000
-## year1995-01-01 1.703e-05 3.833e+05 0.000 1.000
-## year1996-01-01 1.549e-07 3.070e+03 0.000 1.000
-## year1997-01-01 1.551e-07 2.483e+03 0.000 1.000
-## year1998-01-01 1.547e-07 2.445e+03 0.000 1.000
-## year1999-01-01 1.550e-07 2.418e+03 0.000 1.000
-## year2000-01-01 1.547e-07 2.391e+03 0.000 1.000
-## year2001-01-01 1.548e-07 2.417e+03 0.000 1.000
-## year2002-01-01 -2.120e+01 2.261e+03 -0.009 0.993
-## year2003-01-01 -2.155e+01 2.261e+03 -0.010 0.992
-## year2004-01-01 -2.202e+01 2.261e+03 -0.010 0.992
-## year2005-01-01 -2.241e+01 2.261e+03 -0.010 0.992
-## year2006-01-01 -2.251e+01 2.261e+03 -0.010 0.992
-## year2007-01-01 -2.248e+01 2.261e+03 -0.010 0.992
-## year2008-01-01 -2.269e+01 2.261e+03 -0.010 0.992
-## year2009-01-01 -2.251e+01 2.261e+03 -0.010 0.992
-## year2010-01-01 -2.283e+01 2.261e+03 -0.010 0.992
-## year2011-01-01 -2.295e+01 2.261e+03 -0.010 0.992
-## year2012-01-01 -2.290e+01 2.261e+03 -0.010 0.992
-## year2013-01-01 -2.292e+01 2.261e+03 -0.010 0.992
-## year2014-01-01 -2.320e+01 2.261e+03 -0.010 0.992
-## year2015-01-01 -2.315e+01 2.261e+03 -0.010 0.992
-## year2016-01-01 -2.290e+01 2.261e+03 -0.010 0.992
-## year2017-01-01 -2.333e+01 2.261e+03 -0.010 0.992
-## year2018-01-01 -2.335e+01 2.261e+03 -0.010 0.992
-## year2019-01-01 -2.347e+01 2.261e+03 -0.010 0.992
-## weighted_probs 6.439e-07 3.794e+03 0.000 1.000
-## year1994-01-01:weighted_probs -6.445e-07 1.374e+04 0.000 1.000
-## year1995-01-01:weighted_probs -1.794e-05 3.928e+05 0.000 1.000
-## year1996-01-01:weighted_probs -6.439e-07 4.804e+03 0.000 1.000
-## year1997-01-01:weighted_probs -6.445e-07 4.207e+03 0.000 1.000
-## year1998-01-01:weighted_probs -6.440e-07 4.197e+03 0.000 1.000
-## year1999-01-01:weighted_probs -6.441e-07 4.092e+03 0.000 1.000
-## year2000-01-01:weighted_probs -6.438e-07 4.207e+03 0.000 1.000
-## year2001-01-01:weighted_probs -6.437e-07 4.076e+03 0.000 1.000
-## year2002-01-01:weighted_probs 1.041e-02 3.794e+03 0.000 1.000
-## year2003-01-01:weighted_probs 8.070e-02 3.794e+03 0.000 1.000
-## year2004-01-01:weighted_probs 2.619e-01 3.794e+03 0.000 1.000
-## year2005-01-01:weighted_probs 3.303e-01 3.794e+03 0.000 1.000
-## year2006-01-01:weighted_probs 1.454e-01 3.794e+03 0.000 1.000
-## year2007-01-01:weighted_probs 3.719e-01 3.794e+03 0.000 1.000
-## year2008-01-01:weighted_probs 9.020e-02 3.794e+03 0.000 1.000
-## year2009-01-01:weighted_probs 2.450e-01 3.794e+03 0.000 1.000
-## year2010-01-01:weighted_probs 1.140e-01 3.794e+03 0.000 1.000
-## year2011-01-01:weighted_probs 1.726e-01 3.794e+03 0.000 1.000
-## year2012-01-01:weighted_probs -5.235e-02 3.794e+03 0.000 1.000
-## year2013-01-01:weighted_probs 1.263e-01 3.794e+03 0.000 1.000
-## year2014-01-01:weighted_probs 2.076e-01 3.794e+03 0.000 1.000
-## year2015-01-01:weighted_probs 3.125e-01 3.794e+03 0.000 1.000
-## year2016-01-01:weighted_probs 2.667e-01 3.794e+03 0.000 1.000
-## year2017-01-01:weighted_probs 3.387e-01 3.794e+03 0.000 1.000
-## year2018-01-01:weighted_probs 4.747e-01 3.794e+03 0.000 1.000
-## year2019-01-01:weighted_probs -1.915e-01 3.794e+03 0.000 1.000
-##
-## (Dispersion parameter for binomial family taken to be 1)
-##
-## Null deviance: 5755.8 on 163974 degrees of freedom
-## Residual deviance: 4607.9 on 163921 degrees of freedom
-## AIC: 4715.9
-##
-## Number of Fisher Scoring iterations: 15
-##
+## # A tibble: 4 x 5
+## term estimate std.error statistic p.value
+## <chr> <dbl> <dbl> <dbl> <dbl>
+## 1 (Intercept) -2.11 0.507 -4.17 3.06e- 5
+## 2 year -0.000272 0.0000338 -8.05 8.48e-16
+## 3 weighted_probs 0.00264 0.525 0.00502 9.96e- 1
+## 4 year:weighted_probs 0.0000131 0.0000353 0.370 7.11e- 1
##
## [[2]]
-##
-## Call:
-## glm(formula = type ~ year * weighted_probs, family = "binomial",
-## data = iscb_lm %>% filter(region == regioni, !is.na(weighted_probs)))
-##
-## Deviance Residuals:
-## Min 1Q Median 3Q Max
-## -0.1634 -0.0704 -0.0618 -0.0518 4.1800
-##
-## Coefficients:
-## Estimate Std. Error z value Pr(>|z|)
-## (Intercept) 1.657e+01 1.839e+03 0.009 0.993
-## year1994-01-01 5.534e-06 2.678e+03 0.000 1.000
-## year1995-01-01 5.534e-06 3.525e+03 0.000 1.000
-## year1996-01-01 5.534e-06 2.371e+03 0.000 1.000
-## year1997-01-01 5.534e-06 2.011e+03 0.000 1.000
-## year1998-01-01 5.534e-06 2.058e+03 0.000 1.000
-## year1999-01-01 5.534e-06 1.966e+03 0.000 1.000
-## year2000-01-01 5.534e-06 1.967e+03 0.000 1.000
-## year2001-01-01 5.534e-06 2.002e+03 0.000 1.000
-## year2002-01-01 -2.113e+01 1.839e+03 -0.011 0.991
-## year2003-01-01 -2.147e+01 1.839e+03 -0.012 0.991
-## year2004-01-01 -2.165e+01 1.839e+03 -0.012 0.991
-## year2005-01-01 -2.218e+01 1.839e+03 -0.012 0.990
-## year2006-01-01 -2.222e+01 1.839e+03 -0.012 0.990
-## year2007-01-01 -2.229e+01 1.839e+03 -0.012 0.990
-## year2008-01-01 -2.272e+01 1.839e+03 -0.012 0.990
-## year2009-01-01 -2.210e+01 1.839e+03 -0.012 0.990
-## year2010-01-01 -2.249e+01 1.839e+03 -0.012 0.990
-## year2011-01-01 -2.272e+01 1.839e+03 -0.012 0.990
-## year2012-01-01 -2.258e+01 1.839e+03 -0.012 0.990
-## year2013-01-01 -2.291e+01 1.839e+03 -0.012 0.990
-## year2014-01-01 -2.265e+01 1.839e+03 -0.012 0.990
-## year2015-01-01 -2.271e+01 1.839e+03 -0.012 0.990
-## year2016-01-01 -2.250e+01 1.839e+03 -0.012 0.990
-## year2017-01-01 -2.304e+01 1.839e+03 -0.013 0.990
-## year2018-01-01 -2.283e+01 1.839e+03 -0.012 0.990
-## year2019-01-01 -2.312e+01 1.839e+03 -0.013 0.990
-## weighted_probs 5.664e-04 1.271e+05 0.000 1.000
-## year1994-01-01:weighted_probs -5.664e-04 1.303e+05 0.000 1.000
-## year1995-01-01:weighted_probs -5.665e-04 1.316e+06 0.000 1.000
-## year1996-01-01:weighted_probs -5.664e-04 1.673e+05 0.000 1.000
-## year1997-01-01:weighted_probs -5.664e-04 1.320e+05 0.000 1.000
-## year1998-01-01:weighted_probs -5.664e-04 1.552e+05 0.000 1.000
-## year1999-01-01:weighted_probs -5.664e-04 1.297e+05 0.000 1.000
-## year2000-01-01:weighted_probs -5.664e-04 1.271e+05 0.000 1.000
-## year2001-01-01:weighted_probs -5.664e-04 1.496e+05 0.000 1.000
-## year2002-01-01:weighted_probs -9.616e-01 1.271e+05 0.000 1.000
-## year2003-01-01:weighted_probs -4.927e-01 1.271e+05 0.000 1.000
-## year2004-01-01:weighted_probs -1.039e+01 1.271e+05 0.000 1.000
-## year2005-01-01:weighted_probs -1.171e+00 1.271e+05 0.000 1.000
-## year2006-01-01:weighted_probs -4.484e+00 1.271e+05 0.000 1.000
-## year2007-01-01:weighted_probs -4.192e-01 1.271e+05 0.000 1.000
-## year2008-01-01:weighted_probs 2.439e-01 1.271e+05 0.000 1.000
-## year2009-01-01:weighted_probs -1.025e+01 1.271e+05 0.000 1.000
-## year2010-01-01:weighted_probs -7.067e+00 1.271e+05 0.000 1.000
-## year2011-01-01:weighted_probs -1.600e+00 1.271e+05 0.000 1.000
-## year2012-01-01:weighted_probs -5.715e+00 1.271e+05 0.000 1.000
-## year2013-01-01:weighted_probs 7.159e-02 1.271e+05 0.000 1.000
-## year2014-01-01:weighted_probs -2.094e+01 1.271e+05 0.000 1.000
-## year2015-01-01:weighted_probs -4.384e+00 1.271e+05 0.000 1.000
-## year2016-01-01:weighted_probs -3.331e+00 1.271e+05 0.000 1.000
-## year2017-01-01:weighted_probs -1.043e+00 1.271e+05 0.000 1.000
-## year2018-01-01:weighted_probs -4.340e+00 1.271e+05 0.000 1.000
-## year2019-01-01:weighted_probs -2.436e+00 1.271e+05 0.000 1.000
-##
-## (Dispersion parameter for binomial family taken to be 1)
-##
-## Null deviance: 5755.8 on 163974 degrees of freedom
-## Residual deviance: 4516.4 on 163921 degrees of freedom
-## AIC: 4624.4
-##
-## Number of Fisher Scoring iterations: 15
-##
+## # A tibble: 4 x 5
+## term estimate std.error statistic p.value
+## <chr> <dbl> <dbl> <dbl> <dbl>
+## 1 (Intercept) -2.53 0.506 -5.01 5.54e- 7
+## 2 year -0.000227 0.0000338 -6.71 1.94e-11
+## 3 weighted_probs 1.96 2.45 0.800 4.24e- 1
+## 4 year:weighted_probs -0.000239 0.000164 -1.46 1.44e- 1
##
## [[3]]
-##
-## Call:
-## glm(formula = type ~ year * weighted_probs, family = "binomial",
-## data = iscb_lm %>% filter(region == regioni, !is.na(weighted_probs)))
-##
-## Deviance Residuals:
-## Min 1Q Median 3Q Max
-## -0.4562 -0.0634 -0.0575 -0.0511 3.7671
-##
-## Coefficients:
-## Estimate Std. Error z value Pr(>|z|)
-## (Intercept) 1.657e+01 1.727e+03 0.010 0.992
-## year1994-01-01 7.478e-08 2.308e+04 0.000 1.000
-## year1995-01-01 -5.512e-08 5.655e+04 0.000 1.000
-## year1996-01-01 8.130e-08 2.388e+03 0.000 1.000
-## year1997-01-01 8.129e-08 1.961e+03 0.000 1.000
-## year1998-01-01 8.110e-08 1.966e+03 0.000 1.000
-## year1999-01-01 8.127e-08 1.905e+03 0.000 1.000
-## year2000-01-01 8.109e-08 1.885e+03 0.000 1.000
-## year2001-01-01 8.096e-08 1.945e+03 0.000 1.000
-## year2002-01-01 -2.120e+01 1.727e+03 -0.012 0.990
-## year2003-01-01 -2.155e+01 1.727e+03 -0.012 0.990
-## year2004-01-01 -2.190e+01 1.727e+03 -0.013 0.990
-## year2005-01-01 -2.232e+01 1.727e+03 -0.013 0.990
-## year2006-01-01 -2.258e+01 1.727e+03 -0.013 0.990
-## year2007-01-01 -2.217e+01 1.727e+03 -0.013 0.990
-## year2008-01-01 -2.247e+01 1.727e+03 -0.013 0.990
-## year2009-01-01 -2.245e+01 1.727e+03 -0.013 0.990
-## year2010-01-01 -2.278e+01 1.727e+03 -0.013 0.989
-## year2011-01-01 -2.297e+01 1.727e+03 -0.013 0.989
-## year2012-01-01 -2.298e+01 1.727e+03 -0.013 0.989
-## year2013-01-01 -2.286e+01 1.727e+03 -0.013 0.989
-## year2014-01-01 -2.304e+01 1.727e+03 -0.013 0.989
-## year2015-01-01 -2.312e+01 1.727e+03 -0.013 0.989
-## year2016-01-01 -2.291e+01 1.727e+03 -0.013 0.989
-## year2017-01-01 -2.321e+01 1.727e+03 -0.013 0.989
-## year2018-01-01 -2.332e+01 1.727e+03 -0.014 0.989
-## year2019-01-01 -2.366e+01 1.727e+03 -0.014 0.989
-## weighted_probs 1.938e-06 4.121e+03 0.000 1.000
-## year1994-01-01:weighted_probs -1.570e-06 1.402e+06 0.000 1.000
-## year1995-01-01:weighted_probs 2.381e-05 1.069e+07 0.000 1.000
-## year1996-01-01:weighted_probs -1.939e-06 5.341e+03 0.000 1.000
-## year1997-01-01:weighted_probs -1.938e-06 4.763e+03 0.000 1.000
-## year1998-01-01:weighted_probs -1.938e-06 4.612e+03 0.000 1.000
-## year1999-01-01:weighted_probs -1.938e-06 4.443e+03 0.000 1.000
-## year2000-01-01:weighted_probs -1.938e-06 4.469e+03 0.000 1.000
-## year2001-01-01:weighted_probs -1.938e-06 4.384e+03 0.000 1.000
-## year2002-01-01:weighted_probs 3.218e-02 4.121e+03 0.000 1.000
-## year2003-01-01:weighted_probs 8.794e-02 4.121e+03 0.000 1.000
-## year2004-01-01:weighted_probs -9.221e-02 4.121e+03 0.000 1.000
-## year2005-01-01:weighted_probs 5.542e-02 4.121e+03 0.000 1.000
-## year2006-01-01:weighted_probs 2.826e-01 4.121e+03 0.000 1.000
-## year2007-01-01:weighted_probs -7.673e-01 4.121e+03 0.000 1.000
-## year2008-01-01:weighted_probs -8.470e-01 4.121e+03 0.000 1.000
-## year2009-01-01:weighted_probs 3.223e-02 4.121e+03 0.000 1.000
-## year2010-01-01:weighted_probs -8.737e-02 4.121e+03 0.000 1.000
-## year2011-01-01:weighted_probs 1.924e-01 4.121e+03 0.000 1.000
-## year2012-01-01:weighted_probs 2.004e-01 4.121e+03 0.000 1.000
-## year2013-01-01:weighted_probs -1.344e-01 4.121e+03 0.000 1.000
-## year2014-01-01:weighted_probs -4.055e-01 4.121e+03 0.000 1.000
-## year2015-01-01:weighted_probs 1.523e-01 4.121e+03 0.000 1.000
-## year2016-01-01:weighted_probs 2.242e-01 4.121e+03 0.000 1.000
-## year2017-01-01:weighted_probs -1.849e-01 4.121e+03 0.000 1.000
-## year2018-01-01:weighted_probs 2.411e-01 4.121e+03 0.000 1.000
-## year2019-01-01:weighted_probs 4.947e-01 4.121e+03 0.000 1.000
-##
-## (Dispersion parameter for binomial family taken to be 1)
-##
-## Null deviance: 5755.8 on 163974 degrees of freedom
-## Residual deviance: 4609.4 on 163921 degrees of freedom
-## AIC: 4717.4
-##
-## Number of Fisher Scoring iterations: 15
-##
+## # A tibble: 4 x 5
+## term estimate std.error statistic p.value
+## <chr> <dbl> <dbl> <dbl> <dbl>
+## 1 (Intercept) -1.87 0.538 -3.47 5.20e- 4
+## 2 year -0.000286 0.0000358 -8.01 1.19e-15
+## 3 weighted_probs -0.587 0.779 -0.753 4.51e- 1
+## 4 year:weighted_probs 0.0000441 0.0000511 0.862 3.89e- 1
##
## [[4]]
-##
-## Call:
-## glm(formula = type ~ year * weighted_probs, family = "binomial",
-## data = iscb_lm %>% filter(region == regioni, !is.na(weighted_probs)))
-##
-## Deviance Residuals:
-## Min 1Q Median 3Q Max
-## -0.3124 -0.0638 -0.0592 -0.0506 3.7349
-##
-## Coefficients:
-## Estimate Std. Error z value Pr(>|z|)
-## (Intercept) 1.657e+01 2.111e+03 0.008 0.994
-## year1994-01-01 7.554e-07 3.745e+03 0.000 1.000
-## year1995-01-01 7.562e-07 6.822e+03 0.000 1.000
-## year1996-01-01 7.560e-07 2.790e+03 0.000 1.000
-## year1997-01-01 7.560e-07 2.296e+03 0.000 1.000
-## year1998-01-01 7.562e-07 2.334e+03 0.000 1.000
-## year1999-01-01 7.563e-07 2.246e+03 0.000 1.000
-## year2000-01-01 7.562e-07 2.256e+03 0.000 1.000
-## year2001-01-01 7.561e-07 2.238e+03 0.000 1.000
-## year2002-01-01 -2.124e+01 2.111e+03 -0.010 0.992
-## year2003-01-01 -2.147e+01 2.111e+03 -0.010 0.992
-## year2004-01-01 -2.195e+01 2.111e+03 -0.010 0.992
-## year2005-01-01 -2.220e+01 2.111e+03 -0.011 0.992
-## year2006-01-01 -2.232e+01 2.111e+03 -0.011 0.992
-## year2007-01-01 -2.224e+01 2.111e+03 -0.011 0.992
-## year2008-01-01 -2.272e+01 2.111e+03 -0.011 0.991
-## year2009-01-01 -2.248e+01 2.111e+03 -0.011 0.992
-## year2010-01-01 -2.289e+01 2.111e+03 -0.011 0.991
-## year2011-01-01 -2.288e+01 2.111e+03 -0.011 0.991
-## year2012-01-01 -2.293e+01 2.111e+03 -0.011 0.991
-## year2013-01-01 -2.287e+01 2.111e+03 -0.011 0.991
-## year2014-01-01 -2.322e+01 2.111e+03 -0.011 0.991
-## year2015-01-01 -2.311e+01 2.111e+03 -0.011 0.991
-## year2016-01-01 -2.286e+01 2.111e+03 -0.011 0.991
-## year2017-01-01 -2.333e+01 2.111e+03 -0.011 0.991
-## year2018-01-01 -2.335e+01 2.111e+03 -0.011 0.991
-## year2019-01-01 -2.354e+01 2.111e+03 -0.011 0.991
-## weighted_probs 2.123e-06 5.917e+03 0.000 1.000
-## year1994-01-01:weighted_probs -2.119e-06 2.323e+04 0.000 1.000
-## year1995-01-01:weighted_probs -2.134e-06 3.764e+05 0.000 1.000
-## year1996-01-01:weighted_probs -2.121e-06 1.749e+04 0.000 1.000
-## year1997-01-01:weighted_probs -2.123e-06 6.287e+03 0.000 1.000
-## year1998-01-01:weighted_probs -2.124e-06 6.281e+03 0.000 1.000
-## year1999-01-01:weighted_probs -2.124e-06 6.185e+03 0.000 1.000
-## year2000-01-01:weighted_probs -2.124e-06 6.436e+03 0.000 1.000
-## year2001-01-01:weighted_probs -2.123e-06 6.679e+03 0.000 1.000
-## year2002-01-01:weighted_probs 1.885e-01 5.917e+03 0.000 1.000
-## year2003-01-01:weighted_probs -2.208e-01 5.917e+03 0.000 1.000
-## year2004-01-01:weighted_probs 6.477e-02 5.917e+03 0.000 1.000
-## year2005-01-01:weighted_probs -4.335e-01 5.917e+03 0.000 1.000
-## year2006-01-01:weighted_probs -7.937e-01 5.917e+03 0.000 1.000
-## year2007-01-01:weighted_probs -5.026e-01 5.917e+03 0.000 1.000
-## year2008-01-01:weighted_probs 1.947e-01 5.917e+03 0.000 1.000
-## year2009-01-01:weighted_probs 1.551e-01 5.917e+03 0.000 1.000
-## year2010-01-01:weighted_probs 2.730e-01 5.917e+03 0.000 1.000
-## year2011-01-01:weighted_probs -9.967e-02 5.917e+03 0.000 1.000
-## year2012-01-01:weighted_probs 6.868e-02 5.917e+03 0.000 1.000
-## year2013-01-01:weighted_probs -1.147e-01 5.917e+03 0.000 1.000
-## year2014-01-01:weighted_probs 2.779e-01 5.917e+03 0.000 1.000
-## year2015-01-01:weighted_probs 1.308e-01 5.917e+03 0.000 1.000
-## year2016-01-01:weighted_probs 9.547e-02 5.917e+03 0.000 1.000
-## year2017-01-01:weighted_probs 2.577e-01 5.917e+03 0.000 1.000
-## year2018-01-01:weighted_probs 3.356e-01 5.917e+03 0.000 1.000
-## year2019-01-01:weighted_probs 1.585e-01 5.917e+03 0.000 1.000
-##
-## (Dispersion parameter for binomial family taken to be 1)
-##
-## Null deviance: 5755.8 on 163974 degrees of freedom
-## Residual deviance: 4613.5 on 163921 degrees of freedom
-## AIC: 4721.5
-##
-## Number of Fisher Scoring iterations: 15
+for (i in 1:4){
- print(anova(main_list[[i]], inte_list[[i]], test = 'Chisq'))
+## # A tibble: 4 x 5
+## term estimate std.error statistic p.value
+## <chr> <dbl> <dbl> <dbl> <dbl>
+## 1 (Intercept) -1.80 0.527 -3.42 6.25e- 4
+## 2 year -0.000290 0.0000351 -8.29 1.17e-16
+## 3 weighted_probs -0.880 0.849 -1.04 3.00e- 1
+## 4 year:weighted_probs 0.0000628 0.0000544 1.15 2.48e- 1
for (i in 1:4) {
+ print(anova(main_list[[i]], inte_list[[i]], test = "Chisq"))
}
+## 1 163886 4634.3
+## 2 163885 4633.1 1 1.2103 0.2713
## Analysis of Deviance Table
##
## Model 1: type ~ year + weighted_probs
## Model 2: type ~ year * weighted_probs
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
-## 1 163947 4611.5
-## 2 163921 4607.9 26 3.6735 1
+## 1 163886 4627.2
+## 2 163885 4627.1 1 0.14187 0.7064
## Analysis of Deviance Table
##
## Model 1: type ~ year + weighted_probs
## Model 2: type ~ year * weighted_probs
-## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
-## 1 163947 4560.8
-## 2 163921 4516.4 26 44.388 0.01373 *
-## ---
-## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
+## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
+## 1 163886 4576.6
+## 2 163885 4574.6 1 2.0829 0.149
## Analysis of Deviance Table
##
## Model 1: type ~ year + weighted_probs
## Model 2: type ~ year * weighted_probs
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
-## 1 163947 4618.3
-## 2 163921 4609.4 26 8.8572 0.9993
+## 1 163886 4634.1
+## 2 163885 4633.4 1 0.74565 0.3879
## Analysis of Deviance Table
##
## Model 1: type ~ year + weighted_probs
## Model 2: type ~ year * weighted_probs
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
-## 1 163947 4618.4
-## 2 163921 4613.5 26 4.9708 1
type
over and above the main effect of name origin probability and year (p > 0.01).Conclusion
-Supplementary Figure S5
-df2 <- iscb_pubmed_oth %>%
- filter(region != 'OtherCategories') %>%
- recode_region()
+
Alternative approach
+adjusted_citations
across them may not be optimal.
+keynotes_post_2002 <- keynotes %>%
+ filter(year(year) >= 2002) %>%
+ separate_rows(afflcountries, sep = ",") %>%
+ filter(afflcountries == "United States") %>%
+ group_by(fore_name_simple, last_name_simple) %>%
+ summarise_at("year", n_distinct, na.rm = T)
+
+# nationalize_df was not unique, so the left join to corr_authors resulted
+# in (mostly) duplicate rows.
+# FIXME: I was getting occasional crashes on this line, and it's slow.
+# TRANG: fixed on June 3, 2021 using distinct().
+# Also, the duplication was intentional.
+# Please see our extensive discussion on the merge/join on full names vs
+# fore_name and last_name here:
+# <https://github.com/greenelab/iscb-diversity/issues/6>
+
+distinct_nationalize_df <- nationalize_df %>%
+ distinct(fore_name_simple, last_name_simple, .keep_all = TRUE)
+
+# Calculate sum of adjusted citations for all publications for a first-name/
+# last-name pair in db since 2002
+# where the author countries include US.
+authors <- corr_authors %>%
+ filter(year(year) >= 2002) %>%
+ separate_rows(countries, sep = ",") %>%
+ filter(countries == "US") %>%
+ group_by(fore_name_simple, last_name_simple) %>%
+ summarise_at(vars(adjusted_citations), sum, na.rm = T) %>%
+ left_join(
+ keynotes_post_2002[c("fore_name_simple", "last_name_simple", "year")],
+ by = c("fore_name_simple", "last_name_simple")
+ ) %>%
+ left_join(distinct_nationalize_df, by = c("fore_name_simple", "last_name_simple")) %>%
+ mutate(OtherCategories = SouthAsian + Hispanic + Jewish + Muslim + Nordic + Greek + African)
+
+for (large_region in large_regions) {
+ glm(
+ as.formula(paste("honoree ~ adjusted_citations +", large_region)),
+ data = authors %>% mutate(honoree = !is.na(year)),
+ family = "binomial",
+ control = list(epsilon = 1e-12, maxit = 55, trace = FALSE)
+ ) %>%
+ broom::tidy() %>%
+ print()
+}
+
+
+
+
+
+
+## # A tibble: 3 x 5
+## term estimate std.error statistic p.value
+## <chr> <dbl> <dbl> <dbl> <dbl>
+## 1 (Intercept) -5.79 0.160 -36.2 2.15e-286
+## 2 adjusted_citations 0.0555 0.00426 13.0 8.68e- 39
+## 3 CelticEnglish 0.453 0.272 1.66 9.65e- 2
+## # A tibble: 3 x 5
+## term estimate std.error statistic p.value
+## <chr> <dbl> <dbl> <dbl> <dbl>
+## 1 (Intercept) -5.42 0.130 -41.5 0.
+## 2 adjusted_citations 0.0566 0.00446 12.7 7.36e-37
+## 3 EastAsian -3.18 1.02 -3.10 1.91e- 3
+## # A tibble: 3 x 5
+## term estimate std.error statistic p.value
+## <chr> <dbl> <dbl> <dbl> <dbl>
+## 1 (Intercept) -5.75 0.154 -37.3 6.30e-304
+## 2 adjusted_citations 0.0560 0.00423 13.2 5.74e- 40
+## 3 European 0.393 0.294 1.33 1.82e- 1
+## # A tibble: 3 x 5
+## term estimate std.error statistic p.value
+## <chr> <dbl> <dbl> <dbl> <dbl>
+## 1 (Intercept) -5.76 0.163 -35.2 6.52e-272
+## 2 adjusted_citations 0.0563 0.00425 13.3 3.70e- 40
+## 3 OtherCategories 0.357 0.299 1.19 2.33e- 1
Time lag
+
-year_lag <- period(10, "years")
+iscb_pubmed_oth_lag <- iscb_nat_df %>%
+ rename("journal" = conference) %>%
+ select(year, journal, African:SouthAsian, publication_date) %>%
+ mutate(
+ type = "Keynote speakers/Fellows",
+ adjusted_citations = 1,
+ pmid = -9999
+ ) %>%
+ bind_rows(
+ pubmed_nat_df %>%
+ select(pmid, year, journal, African:SouthAsian, publication_date, adjusted_citations) %>%
+ mutate(type = "Pubmed authors", year = year + year_lag)
+ ) %>%
+ mutate(OtherCategories = SouthAsian + Hispanic + Jewish + Muslim + Nordic + Greek + African) %>%
+ pivot_longer(c(African:SouthAsian, OtherCategories),
+ names_to = "region",
+ values_to = "probabilities"
+ ) %>%
+ filter(!is.na(probabilities), year(year) >= 2002) %>%
+ group_by(type, year, region) %>%
+ mutate(
+ pmc_citations_year = mean(adjusted_citations),
+ weight = adjusted_citations / pmc_citations_year,
+ weighted_probs = probabilities * weight
+ )
-fig_s5 <- iscb_pubmed_sum %>%
- recode_region() %>%
- gam_and_ci(df2 = df2,
- start_y = start_year, end_y = end_year) +
- theme(legend.position = c(0.8, 0.1)) +
- facet_wrap(vars(fct_relevel(region, region_levels)), ncol = 3)
-fig_s5
-
-## `geom_smooth()` using formula 'y ~ s(x, bs = "cs")'
-ggsave('figs/fig_s5.png', fig_s5, width = 6, height = 6)
-## `geom_smooth()` using formula 'y ~ s(x, bs = "cs")'
-ggsave('figs/fig_s5.svg', fig_s5, width = 6, height = 6)
+iscb_lm_lag <- iscb_pubmed_oth_lag %>%
+ ungroup() %>%
+ mutate(type = as.factor(type) %>% relevel(ref = "Pubmed authors"))
+
+main_lm <- function(regioni) {
+ glm(type ~ year + weighted_probs,
+ data = iscb_lm_lag %>%
+ filter(region == regioni, !is.na(weighted_probs)),
+ family = "binomial"
+ )
+}
+
+main_list <- lapply(large_regions, main_lm)
+names(main_list) <- large_regions
+lapply(main_list, broom::tidy)
+## `geom_smooth()` using formula 'y ~ s(x, bs = "cs")'
+## $CelticEnglish
+## # A tibble: 3 x 5
+## term estimate std.error statistic p.value
+## <chr> <dbl> <dbl> <dbl> <dbl>
+## 1 (Intercept) 29.1 1.18 24.7 2.46e-134
+## 2 year -0.00210 0.0000749 -28.0 2.09e-172
+## 3 weighted_probs 0.0199 0.102 0.196 8.45e- 1
+##
+## $EastAsian
+## # A tibble: 3 x 5
+## term estimate std.error statistic p.value
+## <chr> <dbl> <dbl> <dbl> <dbl>
+## 1 (Intercept) 29.0 1.18 24.6 5.89e-134
+## 2 year -0.00209 0.0000749 -27.8 1.41e-170
+## 3 weighted_probs -0.708 0.292 -2.43 1.52e- 2
+##
+## $European
+## # A tibble: 3 x 5
+## term estimate std.error statistic p.value
+## <chr> <dbl> <dbl> <dbl> <dbl>
+## 1 (Intercept) 29.1 1.18 24.7 4.84e-135
+## 2 year -0.00210 0.0000748 -28.0 7.61e-173
+## 3 weighted_probs 0.00832 0.106 0.0788 9.37e- 1
+##
+## $OtherCategories
+## # A tibble: 3 x 5
+## term estimate std.error statistic p.value
+## <chr> <dbl> <dbl> <dbl> <dbl>
+## 1 (Intercept) 29.1 1.18 24.7 3.41e-135
+## 2 year -0.00210 0.0000748 -28.1 3.25e-173
+## 3 weighted_probs 0.170 0.106 1.60 1.09e- 1
Supplementary Figure S5
+# df2 <- iscb_pubmed_oth %>%
+# filter(region != "OtherCategories") %>%
+# recode_region()
+#
+# fig_s5 <- iscb_pubmed_sum %>%
+# recode_region() %>%
+# gam_and_ci(
+# df2 = df2,
+# start_y = start_year, end_y = end_year
+# ) +
+# theme(legend.position = c(0.8, 0.1)) +
+# facet_wrap(vars(fct_relevel(region, region_levels)), ncol = 3)
+# fig_s5
+# ggsave("figs/fig_s5.png", fig_s5, width = 6, height = 6)
+# ggsave("figs/fig_s5.svg", fig_s5, width = 6, height = 6)
sessionInfo()
-## R version 4.0.3 (2020-10-10)
## Platform: x86_64-pc-linux-gnu (64-bit)
@@ -2411,52 +2140,63 @@
Supplementary Figure S5
## BLAS/LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.8.so
##
## locale:
-## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=en_US.UTF-8
-## [4] LC_COLLATE=en_US.UTF-8 LC_MONETARY=en_US.UTF-8 LC_MESSAGES=C
-## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C LC_ADDRESS=C
-## [10] LC_TELEPHONE=C LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
+## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
+## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
+## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=C
+## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
+## [9] LC_ADDRESS=C LC_TELEPHONE=C
+## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
-## [1] gdtools_0.2.2 wru_0.1-10 rnaturalearth_0.1.0 lubridate_1.7.9.2
-## [5] caret_6.0-86 lattice_0.20-41 forcats_0.5.0 stringr_1.4.0
-## [9] dplyr_1.0.2 purrr_0.3.4 readr_1.4.0 tidyr_1.1.2
+## [1] gdtools_0.2.2 wru_0.1-10 rnaturalearth_0.1.0
+## [4] lubridate_1.7.9.2 caret_6.0-86 lattice_0.20-41
+## [7] forcats_0.5.0 stringr_1.4.0 dplyr_1.0.2
+## [10] purrr_0.3.4 readr_1.4.0 tidyr_1.1.2
## [13] tibble_3.0.4 ggplot2_3.3.2 tidyverse_1.3.0
##
## loaded via a namespace (and not attached):
-## [1] colorspace_2.0-0 ellipsis_0.3.1 class_7.3-17 rprojroot_1.3-2
-## [5] fs_1.5.0 rstudioapi_0.12 farver_2.0.3 remotes_2.2.0
-## [9] DT_0.16 prodlim_2019.11.13 fansi_0.4.1 xml2_1.3.2
-## [13] codetools_0.2-16 splines_4.0.3 knitr_1.30 pkgload_1.1.0
-## [17] jsonlite_1.7.1 pROC_1.16.2 broom_0.7.2 dbplyr_2.0.0
-## [21] rgeos_0.5-5 compiler_4.0.3 httr_1.4.2 backports_1.2.0
-## [25] assertthat_0.2.1 Matrix_1.2-18 cli_2.1.0 htmltools_0.5.0
-## [29] prettyunits_1.1.1 tools_4.0.3 gtable_0.3.0 glue_1.4.2
-## [33] rnaturalearthdata_0.1.0 reshape2_1.4.4 Rcpp_1.0.5 cellranger_1.1.0
-## [37] vctrs_0.3.4 svglite_1.2.3.2 nlme_3.1-149 iterators_1.0.13
-## [41] crosstalk_1.1.0.1 timeDate_3043.102 gower_0.2.2 xfun_0.19
-## [45] ps_1.4.0 testthat_3.0.0 rvest_0.3.6 lifecycle_0.2.0
-## [49] devtools_2.3.2 MASS_7.3-53 scales_1.1.1 ipred_0.9-9
-## [53] hms_0.5.3 RColorBrewer_1.1-2 yaml_2.2.1 curl_4.3
-## [57] memoise_1.1.0 rpart_4.1-15 stringi_1.5.3 desc_1.2.0
-## [61] foreach_1.5.1 e1071_1.7-4 pkgbuild_1.1.0 lava_1.6.8.1
-## [65] systemfonts_0.3.2 rlang_0.4.8 pkgconfig_2.0.3 evaluate_0.14
-## [69] sf_0.9-6 recipes_0.1.15 htmlwidgets_1.5.2 labeling_0.4.2
-## [73] cowplot_1.1.0 tidyselect_1.1.0 processx_3.4.4 plyr_1.8.6
-## [77] magrittr_1.5 R6_2.5.0 generics_0.1.0 DBI_1.1.0
-## [81] mgcv_1.8-33 pillar_1.4.6 haven_2.3.1 withr_2.3.0
-## [85] units_0.6-7 survival_3.2-7 sp_1.4-4 nnet_7.3-14
-## [89] modelr_0.1.8 crayon_1.3.4 KernSmooth_2.23-17 utf8_1.1.4
-## [93] rmarkdown_2.5 usethis_1.6.3 grid_4.0.3 readxl_1.3.1
-## [97] data.table_1.13.2 callr_3.5.1 ModelMetrics_1.2.2.2 reprex_0.3.0
-## [101] digest_0.6.27 classInt_0.4-3 stats4_4.0.3 munsell_0.5.0
-## [105] viridisLite_0.3.0 sessioninfo_1.1.1Country enrichment table
formatPercentage('Author proportion', 2) %>%
formatRound(c('Observed', 'Expected', 'Observed - Expected',
'Enrichment', 'Log2(enrichment)'), 1)
-
-
+
+
Compute enrichment from proportion comparisons
epitools::riskratio()
.Log enrichment figure
## BLAS/LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.8.so
##
## locale:
-## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=en_US.UTF-8
-## [4] LC_COLLATE=en_US.UTF-8 LC_MONETARY=en_US.UTF-8 LC_MESSAGES=C
-## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C LC_ADDRESS=C
-## [10] LC_TELEPHONE=C LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
+## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
+## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
+## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=C
+## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
+## [9] LC_ADDRESS=C LC_TELEPHONE=C
+## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
-## [1] DT_0.16 epitools_0.5-10.1 gdtools_0.2.2 wru_0.1-10
-## [5] rnaturalearth_0.1.0 lubridate_1.7.9.2 caret_6.0-86 lattice_0.20-41
-## [9] forcats_0.5.0 stringr_1.4.0 dplyr_1.0.2 purrr_0.3.4
-## [13] readr_1.4.0 tidyr_1.1.2 tibble_3.0.4 ggplot2_3.3.2
-## [17] tidyverse_1.3.0
+## [1] DT_0.16 epitools_0.5-10.1 gdtools_0.2.2
+## [4] wru_0.1-10 rnaturalearth_0.1.0 lubridate_1.7.9.2
+## [7] caret_6.0-86 lattice_0.20-41 forcats_0.5.0
+## [10] stringr_1.4.0 dplyr_1.0.2 purrr_0.3.4
+## [13] readr_1.4.0 tidyr_1.1.2 tibble_3.0.4
+## [16] ggplot2_3.3.2 tidyverse_1.3.0
##
## loaded via a namespace (and not attached):
-## [1] colorspace_2.0-0 ellipsis_0.3.1 class_7.3-17 rprojroot_1.3-2
-## [5] fs_1.5.0 rstudioapi_0.12 farver_2.0.3 remotes_2.2.0
-## [9] prodlim_2019.11.13 fansi_0.4.1 xml2_1.3.2 codetools_0.2-16
-## [13] splines_4.0.3 knitr_1.30 pkgload_1.1.0 jsonlite_1.7.1
-## [17] pROC_1.16.2 broom_0.7.2 dbplyr_2.0.0 rgeos_0.5-5
-## [21] compiler_4.0.3 httr_1.4.2 backports_1.2.0 assertthat_0.2.1
-## [25] Matrix_1.2-18 cli_2.1.0 htmltools_0.5.0 prettyunits_1.1.1
-## [29] tools_4.0.3 gtable_0.3.0 glue_1.4.2 rnaturalearthdata_0.1.0
-## [33] reshape2_1.4.4 Rcpp_1.0.5 cellranger_1.1.0 vctrs_0.3.4
-## [37] svglite_1.2.3.2 nlme_3.1-149 iterators_1.0.13 crosstalk_1.1.0.1
-## [41] timeDate_3043.102 gower_0.2.2 xfun_0.19 ps_1.4.0
-## [45] testthat_3.0.0 rvest_0.3.6 lifecycle_0.2.0 devtools_2.3.2
-## [49] MASS_7.3-53 scales_1.1.1 ipred_0.9-9 hms_0.5.3
-## [53] RColorBrewer_1.1-2 yaml_2.2.1 curl_4.3 memoise_1.1.0
-## [57] rpart_4.1-15 stringi_1.5.3 desc_1.2.0 foreach_1.5.1
-## [61] e1071_1.7-4 pkgbuild_1.1.0 lava_1.6.8.1 systemfonts_0.3.2
-## [65] rlang_0.4.8 pkgconfig_2.0.3 evaluate_0.14 sf_0.9-6
-## [69] recipes_0.1.15 htmlwidgets_1.5.2 labeling_0.4.2 cowplot_1.1.0
-## [73] tidyselect_1.1.0 processx_3.4.4 plyr_1.8.6 magrittr_1.5
-## [77] R6_2.5.0 generics_0.1.0 DBI_1.1.0 mgcv_1.8-33
-## [81] pillar_1.4.6 haven_2.3.1 withr_2.3.0 units_0.6-7
-## [85] survival_3.2-7 sp_1.4-4 nnet_7.3-14 modelr_0.1.8
-## [89] crayon_1.3.4 KernSmooth_2.23-17 utf8_1.1.4 rmarkdown_2.5
-## [93] usethis_1.6.3 grid_4.0.3 readxl_1.3.1 data.table_1.13.2
-## [97] callr_3.5.1 ModelMetrics_1.2.2.2 reprex_0.3.0 digest_0.6.27
-## [101] classInt_0.4-3 stats4_4.0.3 munsell_0.5.0 viridisLite_0.3.0
-## [105] sessioninfo_1.1.1
+## [1] colorspace_2.0-0 ellipsis_0.3.1 class_7.3-17
+## [4] rprojroot_1.3-2 fs_1.5.0 rstudioapi_0.12
+## [7] farver_2.0.3 remotes_2.2.0 prodlim_2019.11.13
+## [10] fansi_0.4.1 xml2_1.3.2 codetools_0.2-16
+## [13] splines_4.0.3 knitr_1.30 pkgload_1.1.0
+## [16] jsonlite_1.7.1 pROC_1.16.2 broom_0.7.2
+## [19] dbplyr_2.0.0 rgeos_0.5-5 compiler_4.0.3
+## [22] httr_1.4.2 backports_1.2.0 assertthat_0.2.1
+## [25] Matrix_1.2-18 cli_2.1.0 htmltools_0.5.0
+## [28] prettyunits_1.1.1 tools_4.0.3 gtable_0.3.0
+## [31] glue_1.4.2 rnaturalearthdata_0.1.0 reshape2_1.4.4
+## [34] Rcpp_1.0.5 cellranger_1.1.0 vctrs_0.3.4
+## [37] svglite_1.2.3.2 nlme_3.1-149 iterators_1.0.13
+## [40] crosstalk_1.1.0.1 timeDate_3043.102 gower_0.2.2
+## [43] xfun_0.19 ps_1.4.0 testthat_3.0.0
+## [46] rvest_0.3.6 lifecycle_0.2.0 devtools_2.3.2
+## [49] MASS_7.3-53 scales_1.1.1 ipred_0.9-9
+## [52] hms_0.5.3 RColorBrewer_1.1-2 yaml_2.2.1
+## [55] curl_4.3 memoise_1.1.0 rpart_4.1-15
+## [58] stringi_1.5.3 desc_1.2.0 foreach_1.5.1
+## [61] e1071_1.7-4 pkgbuild_1.1.0 lava_1.6.8.1
+## [64] systemfonts_0.3.2 rlang_0.4.8 pkgconfig_2.0.3
+## [67] evaluate_0.14 sf_0.9-6 recipes_0.1.15
+## [70] htmlwidgets_1.5.2 labeling_0.4.2 cowplot_1.1.0
+## [73] tidyselect_1.1.0 processx_3.4.4 plyr_1.8.6
+## [76] magrittr_1.5 R6_2.5.0 generics_0.1.0
+## [79] DBI_1.1.0 mgcv_1.8-33 pillar_1.4.6
+## [82] haven_2.3.1 withr_2.3.0 units_0.6-7
+## [85] survival_3.2-7 sp_1.4-4 nnet_7.3-14
+## [88] modelr_0.1.8 crayon_1.3.4 KernSmooth_2.23-17
+## [91] utf8_1.1.4 rmarkdown_2.5 usethis_1.6.3
+## [94] grid_4.0.3 readxl_1.3.1 data.table_1.13.2
+## [97] callr_3.5.1 ModelMetrics_1.2.2.2 reprex_0.3.0
+## [100] digest_0.6.27 classInt_0.4-3 stats4_4.0.3
+## [103] munsell_0.5.0 viridisLite_0.3.0 sessioninfo_1.1.1
Race/ethnicity predictions
rename('surname' = last_name_simple) %>%
predict_race(surname.only = T, impute.missing = F)
-## [1] "Proceeding with surname-only predictions..."
+## Warning in merge_surnames(voter.file, impute.missing = impute.missing): 5166 surnames were not matched.
## Warning in merge_surnames(voter.file, impute.missing = impute.missing): 5166
+## surnames were not matched.
pubmed_us_race <- pubmed_race_pmids %>%
group_by(pmid, journal, publication_date, year, adjusted_citations) %>%
summarise_at(vars(contains('pred.')), mean, na.rm = T, .groups = 'drop') %>%
@@ -1684,7 +1685,8 @@
Race/ethnicity predictions
rename('surname' = last_name_simple) %>%
predict_race(surname.only = T, impute.missing = F)
-## [1] "Proceeding with surname-only predictions..."
+## Warning in merge_surnames(voter.file, impute.missing = impute.missing): 100 surnames were not matched.
## Warning in merge_surnames(voter.file, impute.missing = impute.missing): 100
+## surnames were not matched.
my_jours <- unique(pubmed_us_race$journal)
my_confs <- unique(iscb_us_race$conference)
n_jours <- length(my_jours)
diff --git a/docs/14.us-name-origin.html b/docs/14.us-name-origin.html
index f88475f..c71358d 100644
--- a/docs/14.us-name-origin.html
+++ b/docs/14.us-name-origin.html
@@ -1666,23 +1666,23 @@
Representation analysis of name origin in the USOnly keep articles from 2002 because few authors had nationality predictions before 2002 (mostly due to missing metadata). See 093.summary-stats for more details.
load('Rdata/raws.Rdata')
+load("Rdata/raws.Rdata")
alpha_threshold <- qnorm(0.975)
-pubmed_nat_df <- corr_authors %>%
- filter(year(year) >= 2002) %>%
- separate_rows(countries, sep = ',') %>%
- filter(countries == 'US') %>%
- left_join(nationalize_df, by = c('fore_name', 'last_name')) %>%
+pubmed_nat_df <- corr_authors %>%
+ filter(year(year) >= 2002) %>%
+ separate_rows(countries, sep = ",") %>%
+ filter(countries == "US") %>%
+ left_join(nationalize_df, by = c("fore_name", "last_name")) %>%
group_by(pmid, journal, publication_date, year, adjusted_citations) %>%
summarise_at(vars(African:SouthAsian), mean, na.rm = T) %>%
ungroup()
iscb_nat_df <- keynotes %>%
- separate_rows(afflcountries, sep = '\\|') %>%
- filter(afflcountries == 'United States') %>%
- left_join(nationalize_df, by = c('fore_name', 'last_name'))
+ separate_rows(afflcountries, sep = "\\|") %>%
+ filter(afflcountries == "United States") %>%
+ left_join(nationalize_df, by = c("fore_name", "last_name"))
start_year <- 1992
end_year <- 2019
@@ -1691,9 +1691,9 @@ Representation analysis of name origin in the US
+region_cols <- c("#ffffb3", "#fccde5", "#b3de69", "#fdb462", "#80b1d3", "#8dd3c7", "#bebada", "#fb8072", "#bc80bd", "#ccebc5")
Prepare data frames for later analyses:
@@ -1703,658 +1703,233 @@iscb_pubmed_oth <- iscb_nat_df %>%
- rename('journal' = conference) %>%
+ rename("journal" = conference) %>%
select(year, journal, African:SouthAsian, publication_date) %>%
- mutate(type = 'Keynote speakers/Fellows',
- adjusted_citations = 1) %>%
+ mutate(
+ type = "Keynote speakers/Fellows",
+ adjusted_citations = 1
+ ) %>%
bind_rows(
pubmed_nat_df %>%
select(year, journal, African:SouthAsian, publication_date, adjusted_citations) %>%
- mutate(type = 'Pubmed authors')
+ mutate(type = "Pubmed authors")
) %>%
mutate(OtherCategories = SouthAsian + Hispanic + Jewish + Muslim + Nordic + Greek + African) %>%
pivot_longer(c(African:SouthAsian, OtherCategories),
- names_to = 'region',
- values_to = 'probabilities') %>%
- filter(!is.na(probabilities)) %>%
+ names_to = "region",
+ values_to = "probabilities"
+ ) %>%
+ filter(!is.na(probabilities)) %>%
group_by(type, year, region) %>%
mutate(
pmc_citations_year = mean(adjusted_citations),
- weight = adjusted_citations/pmc_citations_year,
- weighted_probs = probabilities*weight
- )
+ weight = adjusted_citations / pmc_citations_year,
+ weighted_probs = probabilities * weight
+ )
-iscb_pubmed_sum_oth <- iscb_pubmed_oth %>%
+iscb_pubmed_sum_oth <- iscb_pubmed_oth %>%
summarise(
mean_prob = mean(weighted_probs),
- se_prob = sqrt(var(probabilities) * sum(weight^2)/(sum(weight)^2)),
+ se_prob = sqrt(var(probabilities) * sum(weight^2) / (sum(weight)^2)),
me_prob = alpha_threshold * se_prob,
- .groups = 'drop'
+ .groups = "drop"
)
iscb_pubmed_sum <- iscb_pubmed_sum_oth %>%
- filter(region != 'OtherCategories')
+ filter(region != "OtherCategories")
fig_us_name_origina <- iscb_pubmed_sum %>%
- filter(year < '2020-01-01') %>%
- region_breakdown('main', region_levels, fct_rev(type)) +
+ filter(year < "2020-01-01") %>%
+ region_breakdown("main", region_levels, fct_rev(type)) +
guides(fill = guide_legend(nrow = 2))
-large_regions <- c('CelticEnglish', 'EastAsian', 'European', 'OtherCategories')
+large_regions <- c("CelticEnglish", "EastAsian", "European", "OtherCategories")
## Mean and standard deviation of predicted probabilities:
fig_us_name_originb <- iscb_pubmed_sum_oth %>%
filter(region %in% large_regions) %>%
recode_region() %>%
- gam_and_ci(df2 = iscb_pubmed_oth %>%
- filter(region %in% large_regions) %>%
- recode_region(),
- start_y = start_year, end_y = end_year) +
- theme(legend.position = c(0.88, 0.83),
- panel.grid.minor = element_blank(),
- legend.margin = margin(-0.5, 0, 0, 0, unit='cm'),
- legend.text = element_text(size = 6)) +
+ gam_and_ci(
+ df2 = iscb_pubmed_oth %>%
+ filter(region %in% large_regions) %>%
+ recode_region(),
+ start_y = start_year, end_y = end_year
+ ) +
+ theme(
+ legend.position = c(0.88, 0.83),
+ panel.grid.minor = element_blank(),
+ legend.margin = margin(-0.5, 0, 0, 0, unit = "cm"),
+ legend.text = element_text(size = 6)
+ ) +
facet_wrap(vars(fct_relevel(region, large_regions)), nrow = 1)
-fig_us_name_origin <- cowplot::plot_grid(fig_us_name_origina, fig_us_name_originb, labels = 'AUTO', ncol = 1, rel_heights = c(1.3,1))
+fig_us_name_origin <- cowplot::plot_grid(fig_us_name_origina, fig_us_name_originb, labels = "AUTO", ncol = 1, rel_heights = c(1.3, 1))
## `geom_smooth()` using formula 'y ~ s(x, bs = "cs")'
fig_us_name_origin
-ggsave('figs/us_name_origin.png', fig_us_name_origin, width = 6.5, height = 5.5)
-ggsave('figs/us_name_origin.svg', fig_us_name_origin, width = 6.5, height = 5.5)
+ggsave("figs/us_name_origin.png", fig_us_name_origin, width = 6.5, height = 5.5)
+ggsave("figs/us_name_origin.svg", fig_us_name_origin, width = 6.5, height = 5.5)
iscb_lm <- iscb_pubmed_oth %>%
- ungroup() %>%
+iscb_lm <- iscb_pubmed_oth %>%
+ ungroup() %>%
mutate(
# year = c(scale(year)),
- year = as.factor(year),
- type = relevel(as.factor(type), ref = 'Pubmed authors'))
-main_lm <- function(regioni){
- glm(type ~ year + weighted_probs,
- data = iscb_lm %>%
- filter(region == regioni, !is.na(weighted_probs)),
- family = 'binomial')
+ # year = as.factor(year),
+ type = relevel(as.factor(type), ref = "Pubmed authors")
+ )
+main_lm <- function(regioni) {
+ glm(type ~ year + weighted_probs,
+ data = iscb_lm %>%
+ filter(region == regioni, !is.na(weighted_probs), year(year) >= 2002),
+ family = "binomial"
+ )
}
-inte_lm <- function(regioni){
- glm(type ~ weighted_probs*year,
- data = iscb_lm %>%
- filter(region == regioni, !is.na(weighted_probs)),
- family = 'binomial')
+inte_lm <- function(regioni) {
+ glm(type ~ weighted_probs * year,
+ data = iscb_lm %>%
+ filter(region == regioni, !is.na(weighted_probs), year(year) >= 2002),
+ family = "binomial"
+ )
}
main_list <- lapply(large_regions, main_lm)
names(main_list) <- large_regions
-lapply(main_list, summary)
+lapply(main_list, broom::tidy)
## $CelticEnglish
-##
-## Call:
-## glm(formula = type ~ year + weighted_probs, family = "binomial",
-## data = iscb_lm %>% filter(region == regioni, !is.na(weighted_probs)))
-##
-## Deviance Residuals:
-## Min 1Q Median 3Q Max
-## -0.5719 -0.0867 -0.0759 -0.0676 3.4967
-##
-## Coefficients:
-## Estimate Std. Error z value Pr(>|z|)
-## (Intercept) 1.655e+01 1.385e+03 0.012 0.990
-## year1994-01-01 -1.409e-02 1.959e+03 0.000 1.000
-## year1995-01-01 -2.093e-02 2.771e+03 0.000 1.000
-## year1996-01-01 -7.601e-03 1.959e+03 0.000 1.000
-## year1997-01-01 1.080e-02 2.190e+03 0.000 1.000
-## year1998-01-01 -1.349e-02 2.771e+03 0.000 1.000
-## year1999-01-01 -4.707e-03 1.656e+03 0.000 1.000
-## year2000-01-01 1.478e-03 1.580e+03 0.000 1.000
-## year2001-01-01 -1.299e-03 1.563e+03 0.000 1.000
-## year2002-01-01 -1.903e+01 1.385e+03 -0.014 0.989
-## year2003-01-01 -1.992e+01 1.385e+03 -0.014 0.989
-## year2004-01-01 -1.946e+01 1.385e+03 -0.014 0.989
-## year2005-01-01 -1.970e+01 1.385e+03 -0.014 0.989
-## year2006-01-01 -2.004e+01 1.385e+03 -0.014 0.988
-## year2007-01-01 -1.967e+01 1.385e+03 -0.014 0.989
-## year2008-01-01 -2.009e+01 1.385e+03 -0.015 0.988
-## year2009-01-01 -1.939e+01 1.385e+03 -0.014 0.989
-## year2010-01-01 -1.941e+01 1.385e+03 -0.014 0.989
-## year2011-01-01 -2.010e+01 1.385e+03 -0.015 0.988
-## year2012-01-01 -1.993e+01 1.385e+03 -0.014 0.989
-## year2013-01-01 -2.030e+01 1.385e+03 -0.015 0.988
-## year2014-01-01 -2.202e+01 1.385e+03 -0.016 0.987
-## year2015-01-01 -2.266e+01 1.385e+03 -0.016 0.987
-## year2016-01-01 -2.215e+01 1.385e+03 -0.016 0.987
-## year2017-01-01 -2.239e+01 1.385e+03 -0.016 0.987
-## year2018-01-01 -2.263e+01 1.385e+03 -0.016 0.987
-## year2019-01-01 -2.264e+01 1.385e+03 -0.016 0.987
-## weighted_probs 4.104e-02 8.128e-02 0.505 0.614
-##
-## (Dispersion parameter for binomial family taken to be 1)
-##
-## Null deviance: 2725.4 on 26441 degrees of freedom
-## Residual deviance: 1994.0 on 26414 degrees of freedom
-## AIC: 2050
-##
-## Number of Fisher Scoring iterations: 15
-##
+## # A tibble: 3 x 5
+## term estimate std.error statistic p.value
+## <chr> <dbl> <dbl> <dbl> <dbl>
+## 1 (Intercept) 4.91 0.523 9.39 6.14e-21
+## 2 year -0.000616 0.0000346 -17.8 6.39e-71
+## 3 weighted_probs 0.0434 0.0837 0.519 6.04e- 1
##
## $EastAsian
-##
-## Call:
-## glm(formula = type ~ year + weighted_probs, family = "binomial",
-## data = iscb_lm %>% filter(region == regioni, !is.na(weighted_probs)))
-##
-## Deviance Residuals:
-## Min 1Q Median 3Q Max
-## -0.4162 -0.0942 -0.0743 -0.0721 3.8535
-##
-## Coefficients:
-## Estimate Std. Error z value Pr(>|z|)
-## (Intercept) 1.658e+01 1.385e+03 0.012 0.990450
-## year1994-01-01 6.836e-02 1.959e+03 0.000 0.999972
-## year1995-01-01 -1.480e-02 2.771e+03 0.000 0.999996
-## year1996-01-01 -1.319e-02 1.959e+03 0.000 0.999995
-## year1997-01-01 3.151e-02 2.190e+03 0.000 0.999989
-## year1998-01-01 9.420e-03 2.771e+03 0.000 0.999997
-## year1999-01-01 2.774e-04 1.656e+03 0.000 1.000000
-## year2000-01-01 1.432e-03 1.580e+03 0.000 0.999999
-## year2001-01-01 -2.079e-03 1.563e+03 0.000 0.999999
-## year2002-01-01 -1.899e+01 1.385e+03 -0.014 0.989066
-## year2003-01-01 -1.988e+01 1.385e+03 -0.014 0.988552
-## year2004-01-01 -1.941e+01 1.385e+03 -0.014 0.988823
-## year2005-01-01 -1.961e+01 1.385e+03 -0.014 0.988706
-## year2006-01-01 -1.995e+01 1.385e+03 -0.014 0.988512
-## year2007-01-01 -1.960e+01 1.385e+03 -0.014 0.988711
-## year2008-01-01 -2.000e+01 1.385e+03 -0.014 0.988480
-## year2009-01-01 -1.931e+01 1.385e+03 -0.014 0.988877
-## year2010-01-01 -1.932e+01 1.385e+03 -0.014 0.988874
-## year2011-01-01 -2.003e+01 1.385e+03 -0.014 0.988463
-## year2012-01-01 -1.982e+01 1.385e+03 -0.014 0.988583
-## year2013-01-01 -2.017e+01 1.385e+03 -0.015 0.988382
-## year2014-01-01 -2.187e+01 1.385e+03 -0.016 0.987402
-## year2015-01-01 -2.252e+01 1.385e+03 -0.016 0.987033
-## year2016-01-01 -2.199e+01 1.385e+03 -0.016 0.987335
-## year2017-01-01 -2.223e+01 1.385e+03 -0.016 0.987197
-## year2018-01-01 -2.247e+01 1.385e+03 -0.016 0.987061
-## year2019-01-01 -2.248e+01 1.385e+03 -0.016 0.987055
-## weighted_probs -1.705e+00 5.005e-01 -3.407 0.000657 ***
-## ---
-## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
-##
-## (Dispersion parameter for binomial family taken to be 1)
-##
-## Null deviance: 2725.4 on 26441 degrees of freedom
-## Residual deviance: 1972.6 on 26414 degrees of freedom
-## AIC: 2028.6
-##
-## Number of Fisher Scoring iterations: 15
-##
+## # A tibble: 3 x 5
+## term estimate std.error statistic p.value
+## <chr> <dbl> <dbl> <dbl> <dbl>
+## 1 (Intercept) 4.75 0.520 9.13 6.68e-20
+## 2 year -0.000595 0.0000347 -17.2 5.24e-66
+## 3 weighted_probs -1.77 0.501 -3.52 4.28e- 4
##
## $European
-##
-## Call:
-## glm(formula = type ~ year + weighted_probs, family = "binomial",
-## data = iscb_lm %>% filter(region == regioni, !is.na(weighted_probs)))
-##
-## Deviance Residuals:
-## Min 1Q Median 3Q Max
-## -0.6484 -0.0877 -0.0748 -0.0665 3.5052
-##
-## Coefficients:
-## Estimate Std. Error z value Pr(>|z|)
-## (Intercept) 1.652e+01 1.385e+03 0.012 0.990
-## year1994-01-01 3.895e-02 1.959e+03 0.000 1.000
-## year1995-01-01 4.088e-02 2.771e+03 0.000 1.000
-## year1996-01-01 -6.996e-03 1.959e+03 0.000 1.000
-## year1997-01-01 1.329e-02 2.190e+03 0.000 1.000
-## year1998-01-01 2.411e-02 2.771e+03 0.000 1.000
-## year1999-01-01 -1.664e-02 1.656e+03 0.000 1.000
-## year2000-01-01 -9.510e-03 1.579e+03 0.000 1.000
-## year2001-01-01 -6.537e-03 1.563e+03 0.000 1.000
-## year2002-01-01 -1.903e+01 1.385e+03 -0.014 0.989
-## year2003-01-01 -1.994e+01 1.385e+03 -0.014 0.989
-## year2004-01-01 -1.947e+01 1.385e+03 -0.014 0.989
-## year2005-01-01 -1.971e+01 1.385e+03 -0.014 0.989
-## year2006-01-01 -2.005e+01 1.385e+03 -0.014 0.988
-## year2007-01-01 -1.969e+01 1.385e+03 -0.014 0.989
-## year2008-01-01 -2.010e+01 1.385e+03 -0.015 0.988
-## year2009-01-01 -1.941e+01 1.385e+03 -0.014 0.989
-## year2010-01-01 -1.942e+01 1.385e+03 -0.014 0.989
-## year2011-01-01 -2.011e+01 1.385e+03 -0.015 0.988
-## year2012-01-01 -1.994e+01 1.385e+03 -0.014 0.989
-## year2013-01-01 -2.031e+01 1.385e+03 -0.015 0.988
-## year2014-01-01 -2.203e+01 1.385e+03 -0.016 0.987
-## year2015-01-01 -2.267e+01 1.385e+03 -0.016 0.987
-## year2016-01-01 -2.216e+01 1.385e+03 -0.016 0.987
-## year2017-01-01 -2.240e+01 1.385e+03 -0.016 0.987
-## year2018-01-01 -2.264e+01 1.385e+03 -0.016 0.987
-## year2019-01-01 -2.265e+01 1.385e+03 -0.016 0.987
-## weighted_probs 1.719e-01 1.118e-01 1.537 0.124
-##
-## (Dispersion parameter for binomial family taken to be 1)
-##
-## Null deviance: 2725.4 on 26441 degrees of freedom
-## Residual deviance: 1992.2 on 26414 degrees of freedom
-## AIC: 2048.2
-##
-## Number of Fisher Scoring iterations: 15
-##
+## # A tibble: 3 x 5
+## term estimate std.error statistic p.value
+## <chr> <dbl> <dbl> <dbl> <dbl>
+## 1 (Intercept) 4.88 0.521 9.36 7.82e-21
+## 2 year -0.000616 0.0000345 -17.8 3.48e-71
+## 3 weighted_probs 0.173 0.109 1.58 1.14e- 1
##
## $OtherCategories
-##
-## Call:
-## glm(formula = type ~ year + weighted_probs, family = "binomial",
-## data = iscb_lm %>% filter(region == regioni, !is.na(weighted_probs)))
-##
-## Deviance Residuals:
-## Min 1Q Median 3Q Max
-## -0.4441 -0.0868 -0.0759 -0.0675 3.4965
-##
-## Coefficients:
-## Estimate Std. Error z value Pr(>|z|)
-## (Intercept) 1.655e+01 1.385e+03 0.012 0.990
-## year1994-01-01 6.888e-03 1.959e+03 0.000 1.000
-## year1995-01-01 1.202e-02 2.771e+03 0.000 1.000
-## year1996-01-01 1.050e-02 1.959e+03 0.000 1.000
-## year1997-01-01 -1.567e-02 2.190e+03 0.000 1.000
-## year1998-01-01 8.719e-03 2.771e+03 0.000 1.000
-## year1999-01-01 1.000e-02 1.656e+03 0.000 1.000
-## year2000-01-01 9.243e-04 1.580e+03 0.000 1.000
-## year2001-01-01 3.158e-03 1.563e+03 0.000 1.000
-## year2002-01-01 -1.902e+01 1.385e+03 -0.014 0.989
-## year2003-01-01 -1.992e+01 1.385e+03 -0.014 0.989
-## year2004-01-01 -1.946e+01 1.385e+03 -0.014 0.989
-## year2005-01-01 -1.970e+01 1.385e+03 -0.014 0.989
-## year2006-01-01 -2.004e+01 1.385e+03 -0.014 0.988
-## year2007-01-01 -1.968e+01 1.385e+03 -0.014 0.989
-## year2008-01-01 -2.010e+01 1.385e+03 -0.015 0.988
-## year2009-01-01 -1.940e+01 1.385e+03 -0.014 0.989
-## year2010-01-01 -1.941e+01 1.385e+03 -0.014 0.989
-## year2011-01-01 -2.010e+01 1.385e+03 -0.015 0.988
-## year2012-01-01 -1.993e+01 1.385e+03 -0.014 0.989
-## year2013-01-01 -2.031e+01 1.385e+03 -0.015 0.988
-## year2014-01-01 -2.203e+01 1.385e+03 -0.016 0.987
-## year2015-01-01 -2.267e+01 1.385e+03 -0.016 0.987
-## year2016-01-01 -2.216e+01 1.385e+03 -0.016 0.987
-## year2017-01-01 -2.240e+01 1.385e+03 -0.016 0.987
-## year2018-01-01 -2.264e+01 1.385e+03 -0.016 0.987
-## year2019-01-01 -2.265e+01 1.385e+03 -0.016 0.987
-## weighted_probs 4.746e-02 1.309e-01 0.363 0.717
-##
-## (Dispersion parameter for binomial family taken to be 1)
-##
-## Null deviance: 2725.4 on 26441 degrees of freedom
-## Residual deviance: 1994.1 on 26414 degrees of freedom
-## AIC: 2050.1
-##
-## Number of Fisher Scoring iterations: 15
-inte_list <- lapply(large_regions, inte_lm)
-## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
-lapply(inte_list, summary)
+## # A tibble: 3 x 5
+## term estimate std.error statistic p.value
+## <chr> <dbl> <dbl> <dbl> <dbl>
+## 1 (Intercept) 4.94 0.519 9.52 1.79e-21
+## 2 year -0.000618 0.0000345 -17.9 1.17e-71
+## 3 weighted_probs 0.0518 0.136 0.381 7.03e- 1
+inte_list <- lapply(large_regions, inte_lm)
+lapply(inte_list, broom::tidy)
## [[1]]
-##
-## Call:
-## glm(formula = type ~ weighted_probs * year, family = "binomial",
-## data = iscb_lm %>% filter(region == regioni, !is.na(weighted_probs)))
-##
-## Deviance Residuals:
-## Min 1Q Median 3Q Max
-## -0.9159 -0.0866 -0.0745 -0.0664 3.5838
-##
-## Coefficients: (2 not defined because of singularities)
-## Estimate Std. Error z value Pr(>|z|)
-## (Intercept) 1.657e+01 2.261e+03 0.007 0.994
-## weighted_probs 3.033e-07 3.794e+03 0.000 1.000
-## year1994-01-01 1.785e-07 1.105e+04 0.000 1.000
-## year1995-01-01 -1.181e-07 3.375e+03 0.000 1.000
-## year1996-01-01 1.791e-07 3.361e+03 0.000 1.000
-## year1997-01-01 1.793e-07 3.340e+03 0.000 1.000
-## year1998-01-01 -6.301e-08 3.035e+03 0.000 1.000
-## year1999-01-01 1.792e-07 2.836e+03 0.000 1.000
-## year2000-01-01 1.792e-07 2.596e+03 0.000 1.000
-## year2001-01-01 1.791e-07 2.567e+03 0.000 1.000
-## year2002-01-01 -1.883e+01 2.261e+03 -0.008 0.993
-## year2003-01-01 -2.016e+01 2.261e+03 -0.009 0.993
-## year2004-01-01 -1.951e+01 2.261e+03 -0.009 0.993
-## year2005-01-01 -1.974e+01 2.261e+03 -0.009 0.993
-## year2006-01-01 -2.010e+01 2.261e+03 -0.009 0.993
-## year2007-01-01 -2.004e+01 2.261e+03 -0.009 0.993
-## year2008-01-01 -2.013e+01 2.261e+03 -0.009 0.993
-## year2009-01-01 -1.952e+01 2.261e+03 -0.009 0.993
-## year2010-01-01 -1.939e+01 2.261e+03 -0.009 0.993
-## year2011-01-01 -2.007e+01 2.261e+03 -0.009 0.993
-## year2012-01-01 -1.977e+01 2.261e+03 -0.009 0.993
-## year2013-01-01 -2.029e+01 2.261e+03 -0.009 0.993
-## year2014-01-01 -2.215e+01 2.261e+03 -0.010 0.992
-## year2015-01-01 -2.270e+01 2.261e+03 -0.010 0.992
-## year2016-01-01 -2.222e+01 2.261e+03 -0.010 0.992
-## year2017-01-01 -2.232e+01 2.261e+03 -0.010 0.992
-## year2018-01-01 -2.272e+01 2.261e+03 -0.010 0.992
-## year2019-01-01 -2.255e+01 2.261e+03 -0.010 0.992
-## weighted_probs:year1994-01-01 -3.027e-07 1.374e+04 0.000 1.000
-## weighted_probs:year1995-01-01 NA NA NA NA
-## weighted_probs:year1996-01-01 -3.035e-07 4.926e+03 0.000 1.000
-## weighted_probs:year1997-01-01 -3.037e-07 9.415e+03 0.000 1.000
-## weighted_probs:year1998-01-01 NA NA NA NA
-## weighted_probs:year1999-01-01 -3.035e-07 4.532e+03 0.000 1.000
-## weighted_probs:year2000-01-01 -3.035e-07 4.466e+03 0.000 1.000
-## weighted_probs:year2001-01-01 -3.035e-07 4.262e+03 0.000 1.000
-## weighted_probs:year2002-01-01 -6.572e-01 3.794e+03 0.000 1.000
-## weighted_probs:year2003-01-01 5.211e-01 3.794e+03 0.000 1.000
-## weighted_probs:year2004-01-01 1.002e-01 3.794e+03 0.000 1.000
-## weighted_probs:year2005-01-01 1.119e-01 3.794e+03 0.000 1.000
-## weighted_probs:year2006-01-01 1.729e-01 3.794e+03 0.000 1.000
-## weighted_probs:year2007-01-01 7.828e-01 3.794e+03 0.000 1.000
-## weighted_probs:year2008-01-01 9.151e-02 3.794e+03 0.000 1.000
-## weighted_probs:year2009-01-01 3.593e-01 3.794e+03 0.000 1.000
-## weighted_probs:year2010-01-01 -5.646e-02 3.794e+03 0.000 1.000
-## weighted_probs:year2011-01-01 -9.107e-02 3.794e+03 0.000 1.000
-## weighted_probs:year2012-01-01 -6.291e-01 3.794e+03 0.000 1.000
-## weighted_probs:year2013-01-01 -6.333e-02 3.794e+03 0.000 1.000
-## weighted_probs:year2014-01-01 3.413e-01 3.794e+03 0.000 1.000
-## weighted_probs:year2015-01-01 9.374e-02 3.794e+03 0.000 1.000
-## weighted_probs:year2016-01-01 2.007e-01 3.794e+03 0.000 1.000
-## weighted_probs:year2017-01-01 -3.198e-01 3.794e+03 0.000 1.000
-## weighted_probs:year2018-01-01 2.793e-01 3.794e+03 0.000 1.000
-## weighted_probs:year2019-01-01 -4.811e-01 3.794e+03 0.000 1.000
-##
-## (Dispersion parameter for binomial family taken to be 1)
-##
-## Null deviance: 2725.4 on 26441 degrees of freedom
-## Residual deviance: 1985.6 on 26390 degrees of freedom
-## AIC: 2089.6
-##
-## Number of Fisher Scoring iterations: 15
-##
+## # A tibble: 4 x 5
+## term estimate std.error statistic p.value
+## <chr> <dbl> <dbl> <dbl> <dbl>
+## 1 (Intercept) 5.01 0.563 8.89 5.91e-19
+## 2 weighted_probs -0.228 0.591 -0.386 7.00e- 1
+## 3 year -0.000623 0.0000377 -16.5 2.45e-61
+## 4 weighted_probs:year 0.0000201 0.0000425 0.473 6.36e- 1
##
## [[2]]
-##
-## Call:
-## glm(formula = type ~ weighted_probs * year, family = "binomial",
-## data = iscb_lm %>% filter(region == regioni, !is.na(weighted_probs)))
-##
-## Deviance Residuals:
-## Min 1Q Median 3Q Max
-## -0.4447 -0.0949 -0.0758 -0.0713 3.7045
-##
-## Coefficients: (2 not defined because of singularities)
-## Estimate Std. Error z value Pr(>|z|)
-## (Intercept) 1.657e+01 1.839e+03 0.009 0.993
-## weighted_probs 9.566e-05 1.271e+05 0.000 1.000
-## year1994-01-01 8.884e-07 2.678e+03 0.000 1.000
-## year1995-01-01 8.009e-07 2.979e+03 0.000 1.000
-## year1996-01-01 8.886e-07 4.288e+03 0.000 1.000
-## year1997-01-01 8.884e-07 3.087e+03 0.000 1.000
-## year1998-01-01 -5.577e-07 2.861e+03 0.000 1.000
-## year1999-01-01 8.883e-07 2.141e+03 0.000 1.000
-## year2000-01-01 8.883e-07 2.228e+03 0.000 1.000
-## year2001-01-01 8.884e-07 2.102e+03 0.000 1.000
-## year2002-01-01 -1.903e+01 1.839e+03 -0.010 0.992
-## year2003-01-01 -1.977e+01 1.839e+03 -0.011 0.991
-## year2004-01-01 -1.930e+01 1.839e+03 -0.010 0.992
-## year2005-01-01 -1.935e+01 1.839e+03 -0.011 0.992
-## year2006-01-01 -1.974e+01 1.839e+03 -0.011 0.991
-## year2007-01-01 -1.967e+01 1.839e+03 -0.011 0.991
-## year2008-01-01 -2.008e+01 1.839e+03 -0.011 0.991
-## year2009-01-01 -1.918e+01 1.839e+03 -0.010 0.992
-## year2010-01-01 -1.921e+01 1.839e+03 -0.010 0.992
-## year2011-01-01 -1.986e+01 1.839e+03 -0.011 0.991
-## year2012-01-01 -1.966e+01 1.839e+03 -0.011 0.991
-## year2013-01-01 -2.022e+01 1.839e+03 -0.011 0.991
-## year2014-01-01 -2.148e+01 1.839e+03 -0.012 0.991
-## year2015-01-01 -2.239e+01 1.839e+03 -0.012 0.990
-## year2016-01-01 -2.195e+01 1.839e+03 -0.012 0.990
-## year2017-01-01 -2.230e+01 1.839e+03 -0.012 0.990
-## year2018-01-01 -2.241e+01 1.839e+03 -0.012 0.990
-## year2019-01-01 -2.251e+01 1.839e+03 -0.012 0.990
-## weighted_probs:year1994-01-01 -9.567e-05 1.303e+05 0.000 1.000
-## weighted_probs:year1995-01-01 NA NA NA NA
-## weighted_probs:year1996-01-01 -9.576e-05 1.950e+06 0.000 1.000
-## weighted_probs:year1997-01-01 -9.567e-05 1.430e+05 0.000 1.000
-## weighted_probs:year1998-01-01 NA NA NA NA
-## weighted_probs:year1999-01-01 -9.567e-05 1.424e+05 0.000 1.000
-## weighted_probs:year2000-01-01 -9.565e-05 1.597e+05 0.000 1.000
-## weighted_probs:year2001-01-01 -9.566e-05 1.536e+05 0.000 1.000
-## weighted_probs:year2002-01-01 1.380e-01 1.271e+05 0.000 1.000
-## weighted_probs:year2003-01-01 -9.030e+00 1.271e+05 0.000 1.000
-## weighted_probs:year2004-01-01 -8.390e+00 1.271e+05 0.000 1.000
-## weighted_probs:year2005-01-01 -2.229e+01 1.271e+05 0.000 1.000
-## weighted_probs:year2006-01-01 -1.632e+01 1.271e+05 0.000 1.000
-## weighted_probs:year2007-01-01 -9.985e-02 1.271e+05 0.000 1.000
-## weighted_probs:year2008-01-01 -1.389e-01 1.271e+05 0.000 1.000
-## weighted_probs:year2009-01-01 -9.254e+00 1.271e+05 0.000 1.000
-## weighted_probs:year2010-01-01 -5.983e+00 1.271e+05 0.000 1.000
-## weighted_probs:year2011-01-01 -1.351e+01 1.271e+05 0.000 1.000
-## weighted_probs:year2012-01-01 -9.074e+00 1.271e+05 0.000 1.000
-## weighted_probs:year2013-01-01 -6.626e-01 1.271e+05 0.000 1.000
-## weighted_probs:year2014-01-01 -4.322e+01 1.271e+05 0.000 1.000
-## weighted_probs:year2015-01-01 -5.494e+00 1.271e+05 0.000 1.000
-## weighted_probs:year2016-01-01 -2.295e+00 1.271e+05 0.000 1.000
-## weighted_probs:year2017-01-01 -6.581e-01 1.271e+05 0.000 1.000
-## weighted_probs:year2018-01-01 -2.467e+00 1.271e+05 0.000 1.000
-## weighted_probs:year2019-01-01 -1.020e+00 1.271e+05 0.000 1.000
-##
-## (Dispersion parameter for binomial family taken to be 1)
-##
-## Null deviance: 2725.4 on 26441 degrees of freedom
-## Residual deviance: 1953.1 on 26390 degrees of freedom
-## AIC: 2057.1
-##
-## Number of Fisher Scoring iterations: 15
-##
+## # A tibble: 4 x 5
+## term estimate std.error statistic p.value
+## <chr> <dbl> <dbl> <dbl> <dbl>
+## 1 (Intercept) 4.70 0.540 8.72 2.89e-18
+## 2 weighted_probs -0.524 4.08 -0.128 8.98e- 1
+## 3 year -0.000592 0.0000359 -16.5 6.46e-61
+## 4 weighted_probs:year -0.0000794 0.000261 -0.304 7.61e- 1
##
## [[3]]
-##
-## Call:
-## glm(formula = type ~ weighted_probs * year, family = "binomial",
-## data = iscb_lm %>% filter(region == regioni, !is.na(weighted_probs)))
-##
-## Deviance Residuals:
-## Min 1Q Median 3Q Max
-## -0.7432 -0.0882 -0.0752 -0.0649 3.5346
-##
-## Coefficients: (2 not defined because of singularities)
-## Estimate Std. Error z value Pr(>|z|)
-## (Intercept) 1.657e+01 1.727e+03 0.010 0.992
-## weighted_probs -2.147e-07 4.121e+03 0.000 1.000
-## year1994-01-01 -6.953e-08 2.308e+04 0.000 1.000
-## year1995-01-01 -7.042e-08 2.949e+03 0.000 1.000
-## year1996-01-01 -7.179e-08 2.437e+03 0.000 1.000
-## year1997-01-01 -7.194e-08 4.360e+03 0.000 1.000
-## year1998-01-01 -4.958e-08 2.837e+03 0.000 1.000
-## year1999-01-01 -7.177e-08 2.161e+03 0.000 1.000
-## year2000-01-01 -7.156e-08 2.008e+03 0.000 1.000
-## year2001-01-01 -7.128e-08 2.002e+03 0.000 1.000
-## year2002-01-01 -1.925e+01 1.727e+03 -0.011 0.991
-## year2003-01-01 -1.998e+01 1.727e+03 -0.012 0.991
-## year2004-01-01 -1.946e+01 1.727e+03 -0.011 0.991
-## year2005-01-01 -1.973e+01 1.727e+03 -0.011 0.991
-## year2006-01-01 -2.016e+01 1.727e+03 -0.012 0.991
-## year2007-01-01 -1.947e+01 1.727e+03 -0.011 0.991
-## year2008-01-01 -2.003e+01 1.727e+03 -0.012 0.991
-## year2009-01-01 -1.937e+01 1.727e+03 -0.011 0.991
-## year2010-01-01 -1.949e+01 1.727e+03 -0.011 0.991
-## year2011-01-01 -2.032e+01 1.727e+03 -0.012 0.991
-## year2012-01-01 -2.020e+01 1.727e+03 -0.012 0.991
-## year2013-01-01 -2.033e+01 1.727e+03 -0.012 0.991
-## year2014-01-01 -2.200e+01 1.727e+03 -0.013 0.990
-## year2015-01-01 -2.273e+01 1.727e+03 -0.013 0.990
-## year2016-01-01 -2.224e+01 1.727e+03 -0.013 0.990
-## year2017-01-01 -2.244e+01 1.727e+03 -0.013 0.990
-## year2018-01-01 -2.276e+01 1.727e+03 -0.013 0.989
-## year2019-01-01 -2.281e+01 1.727e+03 -0.013 0.989
-## weighted_probs:year1994-01-01 8.770e-08 1.402e+06 0.000 1.000
-## weighted_probs:year1995-01-01 NA NA NA NA
-## weighted_probs:year1996-01-01 2.151e-07 5.382e+03 0.000 1.000
-## weighted_probs:year1997-01-01 2.173e-07 2.223e+04 0.000 1.000
-## weighted_probs:year1998-01-01 NA NA NA NA
-## weighted_probs:year1999-01-01 2.147e-07 4.917e+03 0.000 1.000
-## weighted_probs:year2000-01-01 2.148e-07 4.695e+03 0.000 1.000
-## weighted_probs:year2001-01-01 2.143e-07 4.808e+03 0.000 1.000
-## weighted_probs:year2002-01-01 7.054e-01 4.121e+03 0.000 1.000
-## weighted_probs:year2003-01-01 1.685e-01 4.121e+03 0.000 1.000
-## weighted_probs:year2004-01-01 -2.523e-02 4.121e+03 0.000 1.000
-## weighted_probs:year2005-01-01 8.222e-02 4.121e+03 0.000 1.000
-## weighted_probs:year2006-01-01 3.867e-01 4.121e+03 0.000 1.000
-## weighted_probs:year2007-01-01 -1.298e+00 4.121e+03 0.000 1.000
-## weighted_probs:year2008-01-01 -3.182e-01 4.121e+03 0.000 1.000
-## weighted_probs:year2009-01-01 -8.382e-02 4.121e+03 0.000 1.000
-## weighted_probs:year2010-01-01 2.568e-01 4.121e+03 0.000 1.000
-## weighted_probs:year2011-01-01 6.134e-01 4.121e+03 0.000 1.000
-## weighted_probs:year2012-01-01 7.233e-01 4.121e+03 0.000 1.000
-## weighted_probs:year2013-01-01 9.890e-02 4.121e+03 0.000 1.000
-## weighted_probs:year2014-01-01 -1.283e-01 4.121e+03 0.000 1.000
-## weighted_probs:year2015-01-01 2.281e-01 4.121e+03 0.000 1.000
-## weighted_probs:year2016-01-01 3.110e-01 4.121e+03 0.000 1.000
-## weighted_probs:year2017-01-01 1.407e-01 4.121e+03 0.000 1.000
-## weighted_probs:year2018-01-01 4.559e-01 4.121e+03 0.000 1.000
-## weighted_probs:year2019-01-01 5.439e-01 4.121e+03 0.000 1.000
-##
-## (Dispersion parameter for binomial family taken to be 1)
-##
-## Null deviance: 2725.4 on 26441 degrees of freedom
-## Residual deviance: 1985.0 on 26390 degrees of freedom
-## AIC: 2089
-##
-## Number of Fisher Scoring iterations: 15
-##
+## # A tibble: 4 x 5
+## term estimate std.error statistic p.value
+## <chr> <dbl> <dbl> <dbl> <dbl>
+## 1 (Intercept) 5.07 0.573 8.85 9.11e-19
+## 2 weighted_probs -0.462 0.821 -0.563 5.73e- 1
+## 3 year -0.000629 0.0000382 -16.5 7.11e-61
+## 4 weighted_probs:year 0.0000437 0.0000549 0.796 4.26e- 1
##
## [[4]]
-##
-## Call:
-## glm(formula = type ~ weighted_probs * year, family = "binomial",
-## data = iscb_lm %>% filter(region == regioni, !is.na(weighted_probs)))
-##
-## Deviance Residuals:
-## Min 1Q Median 3Q Max
-## -1.1945 -0.0876 -0.0720 -0.0671 3.5058
-##
-## Coefficients: (2 not defined because of singularities)
-## Estimate Std. Error z value Pr(>|z|)
-## (Intercept) 1.657e+01 2.111e+03 0.008 0.994
-## weighted_probs -1.127e-07 5.917e+03 0.000 1.000
-## year1994-01-01 -6.732e-08 3.745e+03 0.000 1.000
-## year1995-01-01 -6.591e-08 3.153e+03 0.000 1.000
-## year1996-01-01 -6.759e-08 2.965e+03 0.000 1.000
-## year1997-01-01 -6.714e-08 4.337e+03 0.000 1.000
-## year1998-01-01 -5.799e-08 2.978e+03 0.000 1.000
-## year1999-01-01 -6.742e-08 2.720e+03 0.000 1.000
-## year2000-01-01 -6.735e-08 2.449e+03 0.000 1.000
-## year2001-01-01 -6.752e-08 2.341e+03 0.000 1.000
-## year2002-01-01 -1.925e+01 2.111e+03 -0.009 0.993
-## year2003-01-01 -1.956e+01 2.111e+03 -0.009 0.993
-## year2004-01-01 -1.936e+01 2.111e+03 -0.009 0.993
-## year2005-01-01 -1.976e+01 2.111e+03 -0.009 0.993
-## year2006-01-01 -1.996e+01 2.111e+03 -0.009 0.992
-## year2007-01-01 -1.939e+01 2.111e+03 -0.009 0.993
-## year2008-01-01 -2.014e+01 2.111e+03 -0.010 0.992
-## year2009-01-01 -1.941e+01 2.111e+03 -0.009 0.993
-## year2010-01-01 -1.963e+01 2.111e+03 -0.009 0.993
-## year2011-01-01 -1.998e+01 2.111e+03 -0.009 0.992
-## year2012-01-01 -1.991e+01 2.111e+03 -0.009 0.992
-## year2013-01-01 -2.041e+01 2.111e+03 -0.010 0.992
-## year2014-01-01 -2.210e+01 2.111e+03 -0.010 0.992
-## year2015-01-01 -2.271e+01 2.111e+03 -0.011 0.991
-## year2016-01-01 -2.213e+01 2.111e+03 -0.010 0.992
-## year2017-01-01 -2.252e+01 2.111e+03 -0.011 0.991
-## year2018-01-01 -2.261e+01 2.111e+03 -0.011 0.991
-## year2019-01-01 -2.270e+01 2.111e+03 -0.011 0.991
-## weighted_probs:year1994-01-01 1.113e-07 2.323e+04 0.000 1.000
-## weighted_probs:year1995-01-01 NA NA NA NA
-## weighted_probs:year1996-01-01 1.154e-07 3.364e+04 0.000 1.000
-## weighted_probs:year1997-01-01 1.126e-07 8.179e+03 0.000 1.000
-## weighted_probs:year1998-01-01 NA NA NA NA
-## weighted_probs:year1999-01-01 1.121e-07 2.601e+04 0.000 1.000
-## weighted_probs:year2000-01-01 1.127e-07 7.110e+03 0.000 1.000
-## weighted_probs:year2001-01-01 1.125e-07 6.870e+03 0.000 1.000
-## weighted_probs:year2002-01-01 9.966e-01 5.917e+03 0.000 1.000
-## weighted_probs:year2003-01-01 -2.696e+00 5.917e+03 0.000 1.000
-## weighted_probs:year2004-01-01 -4.243e-01 5.917e+03 0.000 1.000
-## weighted_probs:year2005-01-01 1.989e-01 5.917e+03 0.000 1.000
-## weighted_probs:year2006-01-01 -3.199e-01 5.917e+03 0.000 1.000
-## weighted_probs:year2007-01-01 -1.480e+00 5.917e+03 0.000 1.000
-## weighted_probs:year2008-01-01 1.422e-01 5.917e+03 0.000 1.000
-## weighted_probs:year2009-01-01 5.110e-02 5.917e+03 0.000 1.000
-## weighted_probs:year2010-01-01 8.824e-01 5.917e+03 0.000 1.000
-## weighted_probs:year2011-01-01 -5.120e-01 5.917e+03 0.000 1.000
-## weighted_probs:year2012-01-01 -7.207e-02 5.917e+03 0.000 1.000
-## weighted_probs:year2013-01-01 3.531e-01 5.917e+03 0.000 1.000
-## weighted_probs:year2014-01-01 2.293e-01 5.917e+03 0.000 1.000
-## weighted_probs:year2015-01-01 1.584e-01 5.917e+03 0.000 1.000
-## weighted_probs:year2016-01-01 -1.152e-01 5.917e+03 0.000 1.000
-## weighted_probs:year2017-01-01 3.618e-01 5.917e+03 0.000 1.000
-## weighted_probs:year2018-01-01 -1.096e-01 5.917e+03 0.000 1.000
-## weighted_probs:year2019-01-01 1.475e-01 5.917e+03 0.000 1.000
-##
-## (Dispersion parameter for binomial family taken to be 1)
-##
-## Null deviance: 2725.4 on 26441 degrees of freedom
-## Residual deviance: 1984.7 on 26390 degrees of freedom
-## AIC: 2088.7
-##
-## Number of Fisher Scoring iterations: 15
-for (i in 1:4){
- print(anova(main_list[[i]], inte_list[[i]], test = 'Chisq'))
+## # A tibble: 4 x 5
+## term estimate std.error statistic p.value
+## <chr> <dbl> <dbl> <dbl> <dbl>
+## 1 (Intercept) 5.04 0.587 8.58 9.33e-18
+## 2 weighted_probs -0.320 1.04 -0.308 7.58e- 1
+## 3 year -0.000625 0.0000393 -15.9 6.63e-57
+## 4 weighted_probs:year 0.0000251 0.0000690 0.364 7.16e- 1
+for (i in 1:4) {
+ print(anova(main_list[[i]], inte_list[[i]], test = "Chisq"))
}
## Analysis of Deviance Table
##
## Model 1: type ~ year + weighted_probs
## Model 2: type ~ weighted_probs * year
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
-## 1 26414 1994.0
-## 2 26390 1985.6 24 8.3858 0.9986
+## 1 26398 2066.7
+## 2 26397 2066.5 1 0.2267 0.634
## Analysis of Deviance Table
##
## Model 1: type ~ year + weighted_probs
## Model 2: type ~ weighted_probs * year
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
-## 1 26414 1972.6
-## 2 26390 1953.1 24 19.504 0.7246
+## 1 26398 2043.8
+## 2 26397 2043.7 1 0.089951 0.7642
## Analysis of Deviance Table
##
## Model 1: type ~ year + weighted_probs
## Model 2: type ~ weighted_probs * year
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
-## 1 26414 1992.2
-## 2 26390 1985.0 24 7.1976 0.9996
+## 1 26398 2064.9
+## 2 26397 2064.2 1 0.63868 0.4242
## Analysis of Deviance Table
##
## Model 1: type ~ year + weighted_probs
## Model 2: type ~ weighted_probs * year
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
-## 1 26414 1994.1
-## 2 26390 1984.7 24 9.3931 0.9966
+## 1 26398 2066.8
+## 2 26397 2066.7 1 0.13195 0.7164
Interaction terms do not predict type
over and above the main effect of name origin probability and year (p > 0.01).
An East Asian name has 0.1817074 the odds of being selected as an honoree, significantly lower compared to other names (\(\beta_\textrm{East Asian} =\) -1.7054, P = 0.00065666). The two groups of scientists did not have a significant association with names predicted to be Celtic/English (P = 0.61364), European (P = 0.12421), or in Other categories (P = 0.71685).
+An East Asian name has 0.1709525 the odds of being selected as an honoree, significantly lower compared to other names (\(\beta_\textrm{East Asian} =\) -1.7664, P = 0.00042758). The two groups of scientists did not have a significant association with names predicted to be Celtic/English (P = 0.60355), European (P = 0.11373), or in Other categories (P = 0.70348).
It’s difficult to come to a conclusion for other regions with so few data points and the imperfect accuracy of our prediction. There seems to be little difference between the proportion of keynote speakers of African, Arabic, South Asian and Hispanic origin than those in the field. However, just because a nationality isn’t underrepresented against the field doesn’t mean scientists from that nationality are appropriately represented.
-df2 <- iscb_pubmed_oth %>%
- filter(region != 'OtherCategories') %>%
+df2 <- iscb_pubmed_oth %>%
+ filter(region != "OtherCategories") %>%
recode_region()
fig_s7 <- iscb_pubmed_sum %>%
recode_region() %>%
- gam_and_ci(df2 = df2,
- start_y = start_year, end_y = end_year) +
+ gam_and_ci(
+ df2 = df2,
+ start_y = start_year, end_y = end_year
+ ) +
theme(legend.position = c(0.8, 0.1)) +
facet_wrap(vars(fct_relevel(region, region_levels)), ncol = 3)
-
+
fig_s7
## `geom_smooth()` using formula 'y ~ s(x, bs = "cs")'
-ggsave('figs/fig_s7.png', fig_s7, width = 6, height = 6)
+ggsave("figs/fig_s7.png", fig_s7, width = 6, height = 6)
## `geom_smooth()` using formula 'y ~ s(x, bs = "cs")'
-ggsave('figs/fig_s7.svg', fig_s7, width = 6, height = 6)
+ggsave("figs/fig_s7.svg", fig_s7, width = 6, height = 6)
## `geom_smooth()` using formula 'y ~ s(x, bs = "cs")'
sessionInfo()
## R version 4.0.3 (2020-10-10)
@@ -2365,52 +1940,64 @@ Supplementary Figure S7
## BLAS/LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.8.so
##
## locale:
-## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=en_US.UTF-8
-## [4] LC_COLLATE=en_US.UTF-8 LC_MONETARY=en_US.UTF-8 LC_MESSAGES=C
-## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C LC_ADDRESS=C
-## [10] LC_TELEPHONE=C LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
+## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
+## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
+## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=C
+## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
+## [9] LC_ADDRESS=C LC_TELEPHONE=C
+## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
-## [1] broom_0.7.2 DT_0.16 epitools_0.5-10.1 gdtools_0.2.2
-## [5] wru_0.1-10 rnaturalearth_0.1.0 lubridate_1.7.9.2 caret_6.0-86
-## [9] lattice_0.20-41 forcats_0.5.0 stringr_1.4.0 dplyr_1.0.2
-## [13] purrr_0.3.4 readr_1.4.0 tidyr_1.1.2 tibble_3.0.4
-## [17] ggplot2_3.3.2 tidyverse_1.3.0
+## [1] broom_0.7.2 DT_0.16 epitools_0.5-10.1
+## [4] gdtools_0.2.2 wru_0.1-10 rnaturalearth_0.1.0
+## [7] lubridate_1.7.9.2 caret_6.0-86 lattice_0.20-41
+## [10] forcats_0.5.0 stringr_1.4.0 dplyr_1.0.2
+## [13] purrr_0.3.4 readr_1.4.0 tidyr_1.1.2
+## [16] tibble_3.0.4 ggplot2_3.3.2 tidyverse_1.3.0
##
## loaded via a namespace (and not attached):
-## [1] colorspace_2.0-0 ellipsis_0.3.1 class_7.3-17 rprojroot_1.3-2
-## [5] fs_1.5.0 rstudioapi_0.12 farver_2.0.3 remotes_2.2.0
-## [9] prodlim_2019.11.13 fansi_0.4.1 xml2_1.3.2 codetools_0.2-16
-## [13] splines_4.0.3 knitr_1.30 pkgload_1.1.0 jsonlite_1.7.1
-## [17] pROC_1.16.2 dbplyr_2.0.0 rgeos_0.5-5 compiler_4.0.3
-## [21] httr_1.4.2 backports_1.2.0 assertthat_0.2.1 Matrix_1.2-18
-## [25] cli_2.1.0 htmltools_0.5.0 prettyunits_1.1.1 tools_4.0.3
-## [29] gtable_0.3.0 glue_1.4.2 rnaturalearthdata_0.1.0 reshape2_1.4.4
-## [33] Rcpp_1.0.5 cellranger_1.1.0 vctrs_0.3.4 svglite_1.2.3.2
-## [37] nlme_3.1-149 iterators_1.0.13 crosstalk_1.1.0.1 timeDate_3043.102
-## [41] gower_0.2.2 xfun_0.19 ps_1.4.0 testthat_3.0.0
-## [45] rvest_0.3.6 lifecycle_0.2.0 devtools_2.3.2 MASS_7.3-53
-## [49] scales_1.1.1 ipred_0.9-9 hms_0.5.3 RColorBrewer_1.1-2
-## [53] yaml_2.2.1 curl_4.3 memoise_1.1.0 rpart_4.1-15
-## [57] stringi_1.5.3 desc_1.2.0 foreach_1.5.1 e1071_1.7-4
-## [61] pkgbuild_1.1.0 lava_1.6.8.1 systemfonts_0.3.2 rlang_0.4.8
-## [65] pkgconfig_2.0.3 evaluate_0.14 sf_0.9-6 recipes_0.1.15
-## [69] htmlwidgets_1.5.2 labeling_0.4.2 cowplot_1.1.0 tidyselect_1.1.0
-## [73] processx_3.4.4 plyr_1.8.6 magrittr_1.5 R6_2.5.0
-## [77] generics_0.1.0 DBI_1.1.0 mgcv_1.8-33 pillar_1.4.6
-## [81] haven_2.3.1 withr_2.3.0 units_0.6-7 survival_3.2-7
-## [85] sp_1.4-4 nnet_7.3-14 modelr_0.1.8 crayon_1.3.4
-## [89] KernSmooth_2.23-17 utf8_1.1.4 rmarkdown_2.5 usethis_1.6.3
-## [93] grid_4.0.3 readxl_1.3.1 data.table_1.13.2 callr_3.5.1
-## [97] ModelMetrics_1.2.2.2 reprex_0.3.0 digest_0.6.27 classInt_0.4-3
-## [101] stats4_4.0.3 munsell_0.5.0 viridisLite_0.3.0 sessioninfo_1.1.1
+## [1] colorspace_2.0-0 ellipsis_0.3.1 class_7.3-17
+## [4] rprojroot_1.3-2 fs_1.5.0 rstudioapi_0.12
+## [7] farver_2.0.3 remotes_2.2.0 prodlim_2019.11.13
+## [10] fansi_0.4.1 xml2_1.3.2 codetools_0.2-16
+## [13] splines_4.0.3 knitr_1.30 pkgload_1.1.0
+## [16] jsonlite_1.7.1 pROC_1.16.2 dbplyr_2.0.0
+## [19] rgeos_0.5-5 compiler_4.0.3 httr_1.4.2
+## [22] backports_1.2.0 assertthat_0.2.1 Matrix_1.2-18
+## [25] cli_2.1.0 htmltools_0.5.0 prettyunits_1.1.1
+## [28] tools_4.0.3 gtable_0.3.0 glue_1.4.2
+## [31] rnaturalearthdata_0.1.0 reshape2_1.4.4 Rcpp_1.0.5
+## [34] cellranger_1.1.0 vctrs_0.3.4 svglite_1.2.3.2
+## [37] nlme_3.1-149 iterators_1.0.13 crosstalk_1.1.0.1
+## [40] timeDate_3043.102 gower_0.2.2 xfun_0.19
+## [43] ps_1.4.0 testthat_3.0.0 rvest_0.3.6
+## [46] lifecycle_0.2.0 devtools_2.3.2 MASS_7.3-53
+## [49] scales_1.1.1 ipred_0.9-9 hms_0.5.3
+## [52] RColorBrewer_1.1-2 yaml_2.2.1 curl_4.3
+## [55] memoise_1.1.0 rpart_4.1-15 stringi_1.5.3
+## [58] desc_1.2.0 foreach_1.5.1 e1071_1.7-4
+## [61] pkgbuild_1.1.0 lava_1.6.8.1 systemfonts_0.3.2
+## [64] rlang_0.4.8 pkgconfig_2.0.3 evaluate_0.14
+## [67] sf_0.9-6 recipes_0.1.15 htmlwidgets_1.5.2
+## [70] labeling_0.4.2 cowplot_1.1.0 tidyselect_1.1.0
+## [73] processx_3.4.4 plyr_1.8.6 magrittr_1.5
+## [76] R6_2.5.0 generics_0.1.0 DBI_1.1.0
+## [79] mgcv_1.8-33 pillar_1.4.6 haven_2.3.1
+## [82] withr_2.3.0 units_0.6-7 survival_3.2-7
+## [85] sp_1.4-4 nnet_7.3-14 modelr_0.1.8
+## [88] crayon_1.3.4 KernSmooth_2.23-17 utf8_1.1.4
+## [91] rmarkdown_2.5 usethis_1.6.3 grid_4.0.3
+## [94] readxl_1.3.1 data.table_1.13.2 callr_3.5.1
+## [97] ModelMetrics_1.2.2.2 reprex_0.3.0 digest_0.6.27
+## [100] classInt_0.4-3 stats4_4.0.3 munsell_0.5.0
+## [103] viridisLite_0.3.0 sessioninfo_1.1.1
Rxz}OfHbGqEoWBN29j}9R%Bz@^Ljz 49o?CmM_E`tEE8B2qi4qpnsln%rAbLiiR2jS70e*876M8}`zTn2 zPoFN1Y&C_-i&X}!% )WU+Zm8M#xy)p5Nz!j$~ACm4F8XB5!6!`_Lmz Hn7Zu)QBk?1cYHuq>^1g(CX^yy=@dMEG%_& z-_?zPEV0?+`K bW@%|@^~|YNRU4dn X<{ zw71V;2FzQZ1 s@~fATnSHX06P-)IuC3vW`2 z6V_OwWDPx^E|aX1JN&aV`x(#hhd;Y`Dd}<>^;XMw)fD1qZ+?A^jKa*q5-;W<%*qPo zPXyvAC^INZW5C>ihUd?p=Nym&&4Q)ljyy$Vc7z>A5h_?hGXPAnW`V-oXPadSLq{}S zQNuz(A#PEFD$l-ar?!%jEiW(Y6nuvy;5s*21QdzfVp!oihiegwfc8uckdz6oqwGS@ zT@awc!oslQI4E)^Uty19Jy!yCLOsBQWVMfe{rVMkx|DE*-O8Sf=86V&j&`KMGwzGO z=f*k=%Ux7x>OMC# ;OpEbiaI<&6lG^;lkKKaFQ(9Q>vU>wOJg8U-wMdA5XsYrA|6QrQP-yT6=Q{J z2uJtDIN_pF(Xwy|GA}^cf&C+nDJ!eOn$JIa6diU?NVv}KZ>T7C(pJmYo z^{b+mILLMRX!vcAj;dAdj~XmTU97F{340>=)}`ccG*bB^3Oq&G;pXbP4pO;MX{Rz1 zFgPmClZ@|J6WmoVv>%sHX-p6oM}k9n9wSR1NxbXblVE5}%_^+UJXY@$6A^)mHxJsg zQC{Y+*p@PaIh@);+NO;J>;gjBhZ1o&vZ|^|G^WL}sKsYvX*i+&eDx-h(c+|}q>YuT zjs23@$fL)PziZ3^5r2E%;lK}&Q>b^5F}{njsnniK{2&}R`6=qOBzd@eJ1Wa=rDG($ z%qgYwUmkwj30b<09c3unW1XBz1W7!S4<*Wc>a#94yG+2JAKlH~nW_*8K9Bf4;*ow~ zSLGIor-OTwDnS9N>bo8Z{D9PNOC;oFZ*iIa`H-@WGVVokH=Qs5YBI1n2OHZHVMj~w zPXK7}+(2T e2OK;?o&u@HzsCP$UJwSFrW zKN0JCuDY0|QvUq{Eom}ISba~A8+gVh07BHha72W}$HzxTc0j$va6_)X0Zko8qY`{y zIMVkLFQYFfSA!86AzFXNzmJZV+7aX}x~r?JOZ;HTc&&d24jwp%Tsguy&YbBa8kb03 z_j#+pWAx|vf8Uw)V|0AW*IYA!7)T Jm79yEn@Aj$!DLkz zYhS)x^T%LhT)!7~BrRLpDPX5KeC4(NINT48PWxc6dI2s;SY$gqP5VHh{^ue1Wv+_A zM+MSHScHlnK|bap3z0x}GBS#NZ_ZcO$`M}oUfTAT3iR48s|fLF!JD?WFefNoe5iJL zeRtR`M)fcT;lZXv)r 8b~}?3PH1;;UgylH}H`Y_i^z zRKpMU)KYJM_|^&j1b+U1hOqOK@(re^3FIWN>BE>v5-&$E8prIVUL!i^hCtCEVg@{V z NJ5LN=7m~QtxZd%+5MEYUpJ=Z5 tJ+H4&;x)D zT`*?G3)5e8ZBR@JImI-ZoQi^RjSysDV6elz!Vpzphy7F9y^`j9G+Q+n>zbM((~$2* zyTYvlk-0)b_E~hp>(JBFL%f{HOf@i=J2v&BzItT7@4tTi>bxPWl`-?|Q_u6Hq+=z~ zPqMPI5Q~X-;KD_9He$!qj KKkRU#w(xau JTM#G=m}+_Y{T1hyUai z;O{@Us+FpW93+g 8pp;gyYbO~Q8OMneu7GW~bdsCf3OM^{ zAT`+OozKxkY54iyc!KE>!D&CArzkRSe*2=wfhW1yWo7=xlJPyV-MxK%aZl*eY8o2C z2RL7ql?j@Eeog!yhLv8yhmDIX%rNx#We)G(b90cBx5&_qBkx|gaN!GY6;fPh=jA08 zZSuzv#N#E<#G8BDf*8arLwf7(I|3aeu)K5}(_9j*1*8T0M0e=W0z?*~ZLF7W#Gw?D zfi}_ecEZrma4cw#Kz(DQL819$$|QUi`O4nO?(Ys%W_+(v&C>yY#r3Yx>jYrfHGqej z$BOkdl%@ql|H;lnJlW|nn|r FgZj>+AdI(Vjhf4oVIV4S7ORf{~KL z*${$5U2X05UAsV$_=Z0q6M>}!SPUkJI0wfqe=e!~k-g$Wx>6ceHvzS4w%$5+Sk>xu z0AC+I2X?;rdWRA@(YAtOLDX%YoDq&A;y*p0cdcO@i$^Jn8eY;D6dLFP(25tM85kHg z5ND8jcq Yll(ZnYx=!^%I^i@fgVSnE6&)TXi#7riGcm{h x4B>XZ=}0j7Pc>@1t4ZB^pco15E8h< zoFM-li~w`t5fLoj${Ff;mCKf0S2wXL#j^1SHr3bnRJawN|5D*Hl|*^$?AcxjFpuxW zx17pH^@1$Doe^|N^6PVBf9LTw$6o#{G9iP*aE?qSKIrf7pH7@tZFO1c2>;NJ+KR1A z9+d%M8NJ1W8|wEKO|E=g*(}e=s4s&2Su0@>3i^_Q5<_3Yk^k-wpVbWA`5Af!uyT zv^?`D@HRGG`O@UJL!S;;81OFqvY08wja1y4G0hNu8Ug@yfr+2!V6=O7z?)9cS}$Do z?~IeMh3>oluoPQh$#DAf?{KVrAmV3y9w>&@3!kjQT&r~HQn6PbEVztPIy&7bwYd%n zL|ufR>yTu+I)}yKT7N3ePp0 gq5uW}PKwMjR4RVbs(`4X7vyBY;1eg->LXr{?Ck zi~L2$_Rybwg{hbHg{c-m^Q@?7#q;OC;f(xMKGTRP_Mtp~f7j!Kl#<7pndv~FrK$!V zKK*R@0X*BpgUtu?)x!d(QS?oFd+4QE*J^=6U^A%K7}xGV#;*%whoeJ)dKaPqsQe1T zG%DT*<3#KlBCYDXzP3^f?@07Vtw+BqSCl3 y+TE_ud@!X5-+<+&27f`uKrEhZG(>V`B-eoi-L-RLsh1^MO4rrFA4WoZ9Tg>hO_W zWLwvF%j=V)07e fElN z)CjM8u7^D{a^Jp{TzXO1#&2|lH2F(NJtNe$y2!9Y!CoA0U#KLdq@=k1z$G87qY$%9 z^~{g-ZdMRN>;#{IcPB39va%Hv%h^In96+QDsBS~H4!+!gsdG=4ah!?CkS@#xgl?JS z7x8iXLu&_vX@@6G1Y*XzY%Di#?e__LJ0L-R^#nET;-%nU4hV!unt1GOhXZBmzxtA@ zJR;xP%d6qH^NyX=s*TY1*RRvxTuC1oj_5g2xA)D!?W#MJcYFM3e2Xcmiy~+ EHDkEG=k~? zz&^c-sO-MJG(1hRAW2@NBZ?vkCGO=`qZinQ;k-lNzS*v?-nuHS*je`8I7pObZWC}E zcXby>xculR0_dM^P)5>OWo2awS!JcAU+&o_b&L^06E(B+^{FEDsRLz}sC?v8L&|cL z_Y`zlMO_C1VSgU_4w$0k&7**TfS8yVB1`}(Gts4Tc6QztzojIB3&OH%@C}{s(|b4E zmXsg)x;j*;b5A$q5T(N+Pv>YA9;|ccqpv5tDc?3UD2i>LUx;9{c>{eSUf(E;>2A;S zY3u3(ogdkPQhSxc>L}>wZB%*<>rCxZ^6Os64sIPJ?mZfq dTHI4shW|Z@?_Ep=IQB%(%NhJ xgWb~8Yx1{bZmY6zjJn7I>eUNyLH)`} zD {{e2RtL{s=E$`n;-yz6#ttwY$Ib6A-PJ|%#I6Ze$ z3=kq|@+hXmP6$PVPLrq6^{|TFYG22g=@oYF+{ryQ#`as1ZSL|YfxtYVB{CSe_m0Z- zty{Lp?RTvZgPa|*VFdn)?;9APjZzIFng9IiJxZgtbQQa+QNU*|Q%2Zzqocw`?|j8P zmif#-e}!Ql96=fCU&KwopWpB^>Xz(e_Vx3#B|`f@>(=dn){y~IBFmL_`t)f|9eB+N zuq?HLP@T`yO+`WV)1ee4zUSI28^M7p+3OLc2I zw9-g@PAGTL&nu6#mJWg4T|>bQBuyCi{(Hw?t);IN(!ku+Tz<=*(fiGn#i@@B^*q0c zjBBZ*0?`&rwu|~}x+7;k)UAT9)gZ4`rn6_uaen<~B=fMqz&fVPNZo`3rx+L#FdtvP zUfH^*diS5REC1R{okIq`Amg4A<2#@nBB%>8ymo^63sRixoaM+5xPEZ|v_;k>O0j1i zhlLG(|9%R0lRqFMEBl9vhc}or$br|p6D1vJIUh 0g~HLFmUEEG?CkmI{k- z|9zmFCyt4=x0X5oLNc``!iy=)-^Y*tP#|Wm#LpGR9<@uNgic2Vv@O8zLn#TiF0prN z4EobvI=TwLOMJU>x{ed{6|!I1&2N%AI? ?MIMdHo0b!^;yZE?!Zu8&M)=6C{?AFuMi+7tWfnfTp^92(-`v`1>s5VcWcn(Jo zsMuXQcV<0xG?7WR1wkh13xmlPGBSV{PzJEtO}*Ys0fR3mG9oK$8Ynpo;3e_3LKs#Z zgO;s}Q5*o#WVVmWXV&RxTZE1N =^7r`vKl =zOjuaeg!nwd~<@zE%5tgmHx{EIviYluAEun!l(@_>p!w-JUF98iS^ z(7hxA#uaWe^}=N!aU)}6@CIjD|0v%#kUYTY(~EnS>*Vkh$%AUlH>pQTwiJj-6(Sa| z@EYss^=9Q&V#77rIRof=t}mG}z9R-3504wdUzVZg$dEsdj*gCr$wDFzf~T-i>CoTd zi1ibb$zl~3_d-F0nv8;i0$F(J2kPPh@V&^?c3MPE8tnQ}S!q#GQ4i#Ki-IncAVWSo z57jka(zdc1M?H |xQN=Z&YaH^{4=pXA;a|<00Qcg zLy`|^`0D_#hLq&uOcZ-hv#{LH{I}=oqyn2Q&v>pbf>}r%L8MTEsL94Am#Vs*kw(Dk z4DvtfYq (JFvWO9-6%uS`OhVq &E{=mbA(A}>^+kQA7HLgE#8bubDYDZSS?fL=GCiD5A4s=!#XZp z1A(d|#>#(A*mz%sfz6T%qx_og=o=dy{ZvnfI3+D;gW2!S1egr4qo}*%buWMEDv=8p z 4sclA5KaJo-E~H#aU{eBY$m;#AZD z7sfCSS$}_j6yB*8D>l3YZkU(86xi#|3tXJEQWWB{-(_9x-Tto(Z(=H9@XXi=hIrb; zQ97_B%jdQE+*zoUQ1*!N9*%z$JkVOzEjRsbPlVIpN5pA@nH;fHFXrV<1lb&g)s7O9 zSgJ!_8I0lxAARwW{s>?+N eyIQfu4NvDx#d!98y`FM#JKtyF!#A;RT;Zsa5w 4i!-MjNTZb>71g?o0D<2xs_hkN6L@L zvltY{POg*vr+&|Q3D+c`d;Ogs+&)0*va*q`flM2w=u-{1(KASFYw868LW2~DcO3Xj z7Iq*!yD)D&J_H}mix)56CXI;x-5j_5Q3=Gu+x00%@KRR0JOdI1s>;dHQ4Ys7c5gX( zc~_V_3?U*SyVhRs`pr`~{(8Z(NSKeW%`TQ%0Rg_s_Q!ph d&Er%jTlU)RR)h4M4`o_#;0%m#|tp^1Ef8d z)*dh?-;p(mkq~`wBWVdr;BoJ4qYthOIZ~t5MJkzl {6pIzjV Q5YD3B8bel-K(=dN=dMgO`c={Im!+5sh#AoE3@v9bxlz{CK1koFDBRi=$#-PJ)uz zynjQc!2KAx#vTJXSsCC6^5v`eur(wkE?K$XGc#&mH&1|x9JD~O@9LDQ%Y^%>tagTK z_H{z4YxOxIeM=#B7Al=|vtg9D{b(g!Pe*US3+6B;6rYn%I>F7X53@bEftd49Gb_BB z2G733or@cy#mDTdwciT*%A+s9O4J0U?%5v3$@9g!0Lni&vGe3lRWpcGY;za-^7Cgq z^)6qAWbR~PN@SlQC#9Yn<_w98<5s-9Vc7J)Bp;U1^63GsUL38H@DF_kS;hDcfj?Da zv;xS-ZgXO|>_29mShE&IoIU0ay8;?Q^_qui1ZE3!CCpgpMCv5OQ;> zz)#n;{WM0GPowu04wGU%qCE#0?~o)s5qg#GUD5h~;{nMW1^AlU8$HOyI!htbq!B-= zWMyRJ_d^+?%3Fp!@n5c^j}MI~Y;HPH+tSbc8G)AQRbjxjU8|F36r z9Pd<#uBofbJf89d23@FQs#)ZWXjf!cq!n}T%wMpL`qPQaTXd74s;c^@#fTb1nUR;W zBQv7Uu^!&>)a_fHx`=u`4ml7qV~{~ZW2&2C1e@%#E%AN5=U49V89YHCM)Dx!{f`Vq z?g|)+3&~EFe9^GLUJgw`*HBEukk{%D20TS-Wyq+1{2ALMD=YQMs|CHwDSxkYdmoVv zjfv?&t3v7sbd`$P76qc;jMnKB1QM_Y0MCS>-w0~Yt>k}w|9;GIlHI!##S76?{vv1> zS@rVV7@GLFJ!lb^K{JJq2-$v@QdJEJ2&jiOpx*Nc9UZqfT%%9WXK=Qqg3<_V4g&u{ zNi l&$b3uqef*t7c?gB)ZWs< eX3zFfkJSSh%XHEr?7PS=6pvxq>Z$elhBU2Tc`8QR$-N1Z|i1 zcR@PLHb05bC5Xl( cl{ z<@cLyhVE0+W~_f-+ae8aghUCT66ll18#J_|lg^*>Des8JD`QAV9Hr5Rjg*RZ1sI|P z%+he-#Q5Y7kgO!)0KH&U^ Zbap~9mGm`>zeRlY%U5ezhW{D2 z@{}ETM>&FMN&`8E$%q%Se^$^+BxGJ&{JFBp1ypPB&=Fs8gJ0H~%TxS$K9c_r5I1^p zaU_WXpi3wHSmU9XHY)v7>qfVbCM$af2UD7{{<-Uo>i)O~vKsMgo2$$Q%<}VDNUuIt z^s7@<*~zGWefN$X?Z@sWJ>e*SUx5})`9s&|qziR2+wS_3b0o@TA7J*oeatu5Tblg1 z_wnExrl&%?76uqT#w)Sb4t~0J%FJ|nP0zrfI;x;|IBM9+Iez9W#%KGMyg#%5XQ(8q z?!Ciw&mM9rtAJw+t-7c;+&^NFU(BMpzX%?<#KhNEPW+}L2Grxb%ZWq0=F>Esm7B|w z1SvW8*X6%se%s9a jsITQ1f(RJ=( z10YN1z# %;#OFEy^D^v+PB0EcL>w&A@#r#% z45Sxbe%WJ k `ea*qN&;+t{d-#YgG 8 aE6*;PY~FChT;9;n5p_*C^NzeH!VbKR CzTg3#Oi@vh?B1PU692~OwwL)yn#09CYY)U%?_9`c3$FwBf|n!I5bT0J5_@c_ z;0({n$&s3ktXefeuK^mN-_&brXn=}y`!vj%`^pH8IP_*5Pl6-?1M$Q_Re~>|vh#Z? ztzL#c@?=WL+ipta~BTtW8H&!7hH{bVXs|ESwZ`OYV8)({&2G6$e@x+pGk@UET8~ zzD^)qekJ=XNk}dFw04tKQ)%f<;qpHJS*|3+^-2?kRrbv1ndDO!yzs)ha>=wrPf2Q@ z-!&j5(3NLrT1wk3`w-TD-}S882o&tUB6bOWUQpo|HN$*i&78%o=J(wu+WQ1z0`;Gl zm&aOg(!^3RCxPvP0V%;3#zrq>#NabnFCGu7k%j~QW3<5V37 M~(Gea*zA^=Q(=a?r+q@iC7pLyuCp(5YBkBnO=X*wZV0OKg({J=*yr(!Q&U3O? z=ok@C-qZ|x!y+1W>oVpDNuPZ$w?X+6kB|SdAue~hCR+L{I&lvM`1pcXLf-BpoE7Z7 zc^){_z#v|a2i=O}4gwuQAV^V0I~gEw<)0Tgnb&^$N%*yE*BlIvqm|$&B_(=7DdN{* zopz}}vi#^I$g}_7Vjz(ae;_8q^m oxqJ zhn<}<34fOMSwFM+3`MxGsK}>qyGh%ZFQ=IJcQc5cbvYL5O RHCGV04gz`J)gGZs;U&n&^hi-4HJ- z!o-yP^5r~UvOtt7_clbGX657z&^~WsLcj=aA$@~(8m3 ?_a-uJymsNW3_4H zZ{Q^ku0(M#HNpSJTlNDxw(X*pkb+&w#KZ(aF%Ffj+;PxkkgLR(YA*=~g@CqrLb1 einRAI$+U&F!*L;yo|v&M~FU%^=an0n{-Ey z756l}iz=txBby>k=PvMjWo`T1=J&G0Gf(RY{_S;LG~uotkAv^MC)nMp=% ^(XB#Cc=$o!bqj0ksE(rkC8s5}JNYy}MCHdr8!(EUv`i !rMa!W5=UdIu~hLKj+NxPDybfKwXyAq=(dw)nfY-8s<$>|g; zM7!rvEges18Esh~eN&-`y@2nqU!?KVxqN??Q brzpZhk8hQah^ee6=eXcXHdM<@=6J%Su>+>mczZB=EIgTW6;v2Q?=vyKUxD zaZd%b_RP%8f%W;~@%(X2l!9PWOX3iSMX_&rcuSAZlvK9b>W&sZZjaS |aQ z EpT6nq;Smt)W*7p}mcLoqg{ z8{n<-N=Gvx!7MExIih&+Glu~X@@#EQO=9AqK{t0dU&Shp13M;P2*ou5`|fLo
EzPw#uVw~F1{&yZR@O>N160`%AI;^Kc`o3E)#9>d$&99Xm!bZoms sCJQ?BIKu<2-yTrNGM96KN=6)Pw(z~Dd7juCa zICzQE?uX$cpzVS>qzNE^=mE^kbcLquI{tFWxoSm!@i6hd46qj5kqC~8suk_-50yJ^ zkf-%jdA84~&swlYlyq-uYN}(DJP*8Oa>vKT5zy@2{iTw=j)0+~Xw}VvXGre&si^Er z0V=7gCoP;O38s8BUCYVJZ0a0Tk*NoE0MR~{J5C%ih*5r0$0F;6J{;`zyOi+79)U@n zI#m^sXlP+2P=;ycQjt#F^NN?Y@7gvU19eFoP@0^4Z!JbMX!Fn1l*Zv|>Vu?i%yl(2 zzCYh~D17pZT ef-XMI}NU|E#! zym;&}+g5}A9Xz9WMI8Mwg%E>fk 4xbKly7 zE-O=!urT@5Z@LpH#~%fyW^k8&Y*9BV?&;#;|B>K%VGCnWD&FW>5$9Wb^V*uVR6+Vy zRvAF>vX~pMZ*Te(EG0FI7o(kIc;b{SHRx2<^|-eTmW?MAde*0?fnom}V{aW*WgGPk zA_j^|i%JPdgCLC{Aso5{kq+tZZj?^x4v_{C=?0Z U+C98Q2Vc@rN1XS4BW4e^h_=mO3L z*rZ`kVvEdVg*jnqK4dh`Lq&sB1B1f@Vm1Mb>G^U3AKKFHSMQ_#Kdc;zz^ dk+(+XE4a89%%Pi&}2f?QW^vj-*C?GLtj+92u$bSWtSmVKN_pJU+F5NrSh1Ad3r z#}p&+@R=$CLGm%@>({U@s{dW2P~ +l^*` zIAbi<3Vp(#s*5HC r1ty6Ik@Yl#cj&-~!CpjPn985{4oh3V909D=y;nsae*;3g)rx z{RH}F#=rF0g+{`{!b_7hJTG11HpIzjiFh2$fT>rjFngu*pAZB-pBp}3Afx}AIV&<6 zC28h-5ze_qWl{hGeixW|kP8dYNLUd9q)rf0*sQ}``{gx_9`_N&XNaW$-(S>5E--GJ zKz&31)Xfjs9Yu--uJ^J*f>Kf LN5xAj7%{K$Ymn4Y@0#8 zwU$59)hr8+IeB?O1kxvDY`g*vy}(?!N8Ia?>L8H5yMENZf%+cAWGp*&FFY VAmtf z2K|rqScyLB`%i6;K&qkj=FQU>3D3jQG<(iIAxvHhtlOZtAi`*icq%%&0`(0PpL)}4 z{#Dwp&rz1wvy&`fzbF3;t~_v11Z`O%PZPXOZ&HXnrNRl^gl8T&{pi?O5?&VtC_OhE zQe;5T_l@>VHgJxe3%i@KP?*7ra0IYsBeXJ|of7q$TSsRg>4Q4{D&xt85r}!!QMXUr zKu3-7ui|E9Ev+5HLvP9W$maqrmOF?L*xA{sqKu&{g^IZTPZ&?Mggzb*MoM}=&A$B7 zX_L$=7GgN6!Uyz InP8?E4po2-x;F`sO8$?X6|W1@+rTraFK4MrfdiX% z*q?52bASfHWN6|8iFHRPdDBe|bS+zOfxm`m?N6R1V`lilEe54K#7J-`lMC!(ztpr# zl%TF}H-8mZRmLvjs)sx)C@5YH-uM|P%-{9+2?+_eKsW_EAn 4qt zS@J|^F3LwVyBxW&g 7L;8}I^4+B^?3 P?zWq$Y~ zkMtAMrxqF%0l}uG=O_Z5^fWZLXZq1=w2_~(d{^v0LDjaO?f`BzQdHp;F9e@L8k;O* z_^-%|YMNWNxkaMpwK-Or^AC|a25#J~1^LeK3MWNko&pSzJoJ~d@#H@C#py^RBLl-< zp5!YXHUcmO1?JXfK0^^;$ 8Bg95lYdC=v>5Rw#p`tHeqC525&~PXj zon`m;;2emf2f~?~e*O^)A}DSrAOTTzySUOO`e~w=7%(Vs6e^xIJ-PfAS0HEst79ZK z^18k}Ll(8*sKUkG^Mt!;(V vjSa=c)*!l_q7>Ol1OzdYV;_*dmlKViSGck-A*)OlxcT>6%A_2rkR65GMlC(h z2xs1-2#oN-0s^b~0rW_Pd3keHVv>?w|I$G
!T_=4A*aJqV Ht4<|&CSneNs#hrOhaN&I2y;7**txd zIX@^OqGi9s+MJzDMQ&3&^{;ae?}YOlDUS3s{!<%3t|QX(x_i_Nr!9Opd~@BjN;iZc z2$pGi5}b^1qWEZ iXn+aaruq5ybh9dZjgaN`W#4x zj>B&1e7C6d1PKrV{8_A!< b73>pL=fXsHYL&GmSTeii>VA!+LZR2% z+6uOViw#Rp wR&_2JvN{smjNF&ZN{Whyo%BUWycw|%oyIOPTzLV@xwkH2v{e_r zy)~|*B5ijsMB`_;$4ozfhCrDC!`~b3GZ)Hp?RYLU rcpYlOSdy+ZxhEpF=M06?r)#k&6+j1|J}ioevEO0)J5Dbe0;+pxCM&l z6O0>|ccZQ&1!+K?Lku%?XV%w|;wmCCX7s%4oIkL9Z#a+o@% 0vi3|+W?K|xV7%Kvk9ZyOx{GLHwJ{q*FP1u4_TzE^dpgb=x zFrZG;^2_i4 $IBIjQFEjNfh`pMl^F7~O zTdv$jV)0(hAX4825DTXtko&CHG{+eJXdp{~%wI$qEMF8NQ4ck<{QqKQqg`aKsq+VN zzme38z$Zojj|20ng$fTo6JCFaX~j$cf@PJ-iE_5K0gs7%Ui6k#Kmd>euo#%L5rDM^ zx_a2kn7>m9c %5Y}`cgnCoWlhkvvBb*=gFCCxv^>|9WX%>Y>^1HX^k z4}dTm)WgqhD=I2LqPYYdmW8=FNIe04hOH37iQxItj)fwG?%6X?w0r@pWmgv=EUBQj zgnlbt7^*ef(|zL@iMPtikEUd =9$|PrL0&gR;R9 zR{4%jgHl)>k-HHo@|Hl+5O|Wo{0oL~SWIy<+FHOC1etDVNe^PW7EH<$Pl#^c0BEw@ zpof6ttOKbSg##qT3YW+N4^iR^TLTLVMr@c{fQ%4?6WC=A&Pedo>KhoqYzOL2VsGU& zYGN%#1qEQ`<7Xa`o1cz2Uqa9ed`8 yZJVp#_a!GaJ_waEbuo zTLUUnED1HZ0W5Hc=(xC}F0D%&P|0-f1I( yS)zCBtM6B@EXDqb2BzdY!TuZCZQIchWjiTS@Q(HtD z2qQ<3rw|_We+!J*%+FWl7`Em4P~Gli`OrXKd&b7`f;>-$kW%zx)C%Fo^SEnPk#rz?-MOMgAFt4N3r2 zD<~i*K@@A){xrZj+s7uT?#J`?xl6$gN^A%Y9~k)QNGdmPgECKwe9XKWyuPUH17!vK zxI15(G6OPk>t|p0A~!_Xuz`TAmlYTwFwnE4dAC{X-0}cOYq?HLM8r>kl$5+)^U2GL z6vycumwKqTH!RetW+k_H8bSn@7L%FR>7&kEa5y*#jBVw_xS YsW zVWGZ=2?3v*(>8+Y4l?0Tq=G-(#sxm|^v57oska8cFw+rnN+{LC#kSKqOwUngauf>G zLnpm60&=2E NHG75iA1&?#!Dsr zcyhVz4=dEaW>5My*;{w8EW@ASh1|VYQ3Q187#R)1kaHULi0^I(n!o*bc#d#*1YjzK z^8YD$mOL^~3!b^`1p`Njw;t;QoA;^(M~no-w=v+~Efv$$BthQp5DBbSG*Sd=dsv%b zafWR}1BlW;^wQ4iITcL@`8;|y500R{IH*Fu(o5YV-3jUCm|>;_O9cM4;7G8Ez6iDSai z5(}NcOcbgv9Jn|lb12)wW={a)X #B2o8M;@0-4+1 zj|H0p58FwgCK_k9_D)VMYBh^EkJF&%?|QBX;PuuKYp{P{s4=Q7Jh!)G@7+M{E UWc1#erZSxM zmg~N6+&9DPLHR+4PQ)uXuE$nI53KoiM9P@R2=!h3F_~&gU%4UVhJO5vbJ%_xpZ9WI zx(jl-mS;G{v=vE=3&Swe`r14&ZlDtB2G5g*%4C)&UY?OIlqgu;!HQ`6V_iGdy|f)F zS@mut+^ErY%f2COzo9Bp#;#I3j{@Mlsylpy&v~yM`Yu~g&p(C#A?BMbjGKRB{^{TK zhZAeJ$(RXzf%8^WRK)2A+8ha?>0|ph{csEKT9{kXGqmDPS#K7Pq5{Sl7!daI^73`W zdjYGiXO#>pB|f~@h@mhApS^sZ9hCdYM^wqMQ~F+=z1-XMgxxT*+WOCd3#=>?8(|X; z<{X4^V_y{MZF_CFQ+Qh|fGGvV6i&p<7T^%?1-w)ap_N~?R)Kx5+pJu*whq-dnLB~D zo7TJf=>ADX*Rxw`I+GThuB&%S9~k7#PR*5&IXIee5Vp+>8I&>ya}i-|4uAi%<5ToP zqm^iJf^kH*aw6fyh67LAs5;}r(g3(Dc8)JeFpq>k)2zN*&-)^}^Fj<0vReLO=U};T z1{lFw#fMkyDV<4d-_oBdh3!`a#)XSEsM f;<)e$ zTu}pagMj!<>svZ-IomhG(s-cn3}g zZbogGPe%qbr|=g>OKMPH^TV)}O~u~dCc!NLLX(F+<2(28%p8eW4oc)% %TKSt 12`w7hi0v^!6AjpW6CvakMD5K8g a3^b{44r36lw#6*I6 zO@44;AXySlML^@F-T`8_o^%jC8=rKhbw2+Ph*9u)05M26gDt!W7Y*GxWJT%c%Oqs0 z;W%kY4W(Z28WOZ@#sav~s&)$SSA%gL2s4n+=(kdT9z^7@M1?~hVFX^1AAcBbPQbGW zo)n3!M1D_?Z9BN*Gv7eGw5H}Xl(R1ytkAJ+og$HMzWkIgO}1zRHv+tv1;cLs(2|D} zg+B^KE-n8UpUO#_E)Wv*xW78Qc|xM*URv}JIM&TINA48ZZ74s3F8jvJt+z@_z$kaE zqZ}b9Z~`Ip$O;fnk)83o0-b-RZ{dq;B0uU4NIL?IJTm%~)zzgCus^7k=!wtVfPFLE z1W2G&I5hCeA&}c7>_MEL^~V?t+wbnC6tg0tA6)?P3*P{qD>*rAo*fX^xB)&T#JvIF z>HnME0L&5T+4 e_xxKLFs^a-$gp+@ zK)9Dmj)#5|o;K3kW1Zj>teqf^k$B+Y$jzG6(@jnn6XaSg4*0h+QN!YmExGnN+&WS1 z-#|UYzhEc0-!_^oLtPpmn7%xi*o0pWM^vI9b|+VXF04LZk|L5MxaR0HjS2i?V>z@m zbU9B;b36yt#@0CGdCTzE+4_UIv%OLM9!%<3Vp**?Ir &PyQv$L|SAG?X?3<_Iw1YlO=R*w)tZSA2*gO|cH38Q`1-Rz!M;tHmJ z+f#yyZfo~QN45@B0z|r?w+3>89UFqkwE;<2#}5>=de{VyWWlgsT3R)xF_Z@;UAiB6 zTZwphf6G@+5VwRH;Vm`Idb5EWrtO3$4RcIz`QkatIziP)haR@P1!(oPdHObZez?<0 zP(*r+;&K1avfR9V`P>IMm^6tRSB3@XAp_@|329;RJ7j^)wkp}_O~}j^XHw5PzxNYK zT^6FpKqp)*FnNN`fjnXvn((#af?rsW&
c4`x6)KBxq-v&J5##?P zULb!0xWQKi!-`OPj($VJGXQ9w!i~CB^#>zN+tJ@Vm|&^2ban 5yYA;12|2 zw573Y10U9cU)kh;Hs-}!Us}wZ22ijxQSUa8{_Q`0#%h0*6qdU`pAlDVU#)P=va&nD z$Vn{RV{BrgT550~8>)Q>>5qVGhI+*Fd3NGsucHN2K^R)XxQ91zrcR#ATQ`I?YoVrM z&6C!wiO~ga7PDD==Gr-w|NPc_S*9*x-h){DtsH->jSmpqQXb<=d%WObnOky$s(7>y z9^1EcF#r*-1CSlNiIef;dcwluqRl-s(0 H@H5Qezsc!e0>ggi>%lLT1?)ZLZon1%`Gcq@? zw_Cx%h9N(QJO!#Vl4~f4%r<)6SFf(ggnWO<=xCAq85Sw#hi4*BekDj`>8NP|hEh>G zHd=H8^_3MPdkDy{poGKihbYMpYLpO4=!c{=x(tLL!)6VmNDzN_oK{y?$;a<{JqOT+ zdnV8E1#{p!OqBqgNBn PUhh|r#25bR-OnLa2 z*M)5_Jv-a4i;7g-x?8N3!&*gAu>w#-FhoN+O=T~^7Gq^h_8?6RAp3S8yEH*iD<_AH z!<#O;Z`i^WyzfUwY5ANy-~PW-;iI%tJ^oEwDBcpNU!M2_AV# 3u~9_A6+rnUDHhV&0o)A)y4gTta?waF=pSZ1#2y0^202y; zu(rfwqqyt?ZUhVOlk-_Lp!x2e0G=TSydbh{-`F7Y&M5wfu2Dhl+l)-Ehg9E+3J&r_ z(1s((??esgG}bI8nsUoZiF7v=4z9dJijI)b2~zh&c_MRwO@WF+>g&tjVxK_b>(f$i zX?LCthlZu*6vHRuBS L+jzz`CTr zc7k*qyGsWk+z@)DW0&Vk{H0H)aQl(ca?=q;astTB0(9B6ToTK%euZ|o$n(PVKfDWc z%fC4LTZH1L5RM7;J_FaiC_cPVT%C*DkqZB^j7$fFNCfazHe83D-0P7*aQ|`j|1XlP z09v1ak|uCsIR5|T$(l60xxOpcJSqg6BstE3@HHh|O&?3nM6Skm*o@Z?T#A690|s-C zSn{455Ydcr!H=upXFojw -*5GSMclNk^*w@CQ^=PS*R$O7Fyjij#aLqd`kOA zkwq4r$fwPZCQ@NUpghM(Zvfe46|lyiKhJ}@hj*qJ@zlcb#pcw9BYO(w2%Xv}qYqu7 zx(`|o2mmdwK$1n y)+L;@roexvD9wisaHosw zKplPeZeDkyHwmoc;pG!0qbFb((8{m9W_!7Hb@|bG#R m4<<~mVV7N0-uFsBLW#zPPKF^kmxzsG`s^qVAdE)%)?QU zv|P;kK3ev!jWJj9t}lu^o#NYExC$T$Lc-+8XbO1RYr350pFDCTDhesk^bEb|HV}JU zD*ExG`#;V5uh7{f=E{*ST?#U|hytr(LCKPtm=oRpL!0z(NU 8im{Mr|4r0aTS|{o#0ZU8%NeOCzkwZR zb=8FYbZX1v2qfgxmhOH&E&7kL{^mQGl>zQ_jhuYf$=%tpN1MvxMOoU(r^t>nn=F~U zA-ZVpEQy@gCJ!WvN1oK>_g42eJJ;Whl(ZYchV@!ucRKxgA--+!*pQxu<;wF)`*_}l zD>m%6T?+)&AKJDBNtw~ZiGD6=opxy0JNo{Tb>P1Q#6Z7X6>`8t5nxmUGClY;&APO) zzw^u?5*UT{BXAj_BpwDL !)aNgQ73FX1mp^nx z6E+SpYJ1l0;x{;=ntsoHtF((x?B4hhe60S}MJP#`=B3SL^42I>p b-^dd!C^ij0i#PY>BA60Us%V zKJargnzT8`Hi$7?{Rh8tJ|Yw*Pn(&ss7y8bf|;O2iba#fQZCg>AwJAfEl;J?nTenC z=Q!diW8UA3CaoYJ4Lh;&`Ie>BlanB`(2AnN=h>QI3 gTTxnrX}%p3eVHq$*e2z9zaO{)ysPOv0KJ>vkD>mSTE$$FWDy zlnHOrAvU%K>t^Sg!Nusjm1^=KMo!bBodegpV%*@o-Ev`UECws<)mSAzDS2G_)YbvX z$!^!#R-Ggb7W&N(rg1YxpxwX;RA-0~&txfXo-b2N!AKWqw?o&lADz- YhtL&gXT!GKAnub(iY_iUBOm`QNhJH8|HQ%B0=>J3XRPWP=J9>1 z;I^fQK?pccPi`MyTwX7_bRU{O`4Ox)ZOs}btnZY9k^QESqd<7Jj1^auNprhRW2W!v z14+U0)N!?g`y#y98~Nvp^fN~f5g$j|vbky%_ByCNij>&xwqMDjwYgiGn3Rn(pr&u9 zy`Aw(c(!p!IL5pCLcRNW$*m;eW&7gVJ@t-?hkR3H_`P^As2&11MRFkG~gK0?^Rd zXIL{NgI{uJq*Cf&SV0_he+BeCc_VoZ7JN&yPd3r%w!;39;LZMLAOzZ|{nc#S~*on#W>f zzmgU#FT-qAx-ta=C3KGEV`eykK1D?Op1q;1dRE&328*NauS~9JXc?s3CA)9wrx -=19*YS2JwM_NBN0N^iTB zb+g{cU1}|l{)L^YIua _JOluraH_li(I-d)YVRy9PL?{49U~!y=b{)^k z#hbwLR{M%I>l0zrbS(z$L~RL#PPH;!oY+dEXeRDmpC3h;%z4B3)ssty_kLw-urqv2 zV8jg`D_(CF?Dwp1p(3qT-fhX-yYGjQF18;Z5%p$s~+;|+d4Q4XkWCcer7 LKwg1%uRNFZcQMQu_A)v_v zgQcvl?jJB6US1wG0*8r$vc2xzq{ZPF53 |}( z5Q$Sa7i>aGh350S{BL-Em%V0F#miW7+e{G 8{(=b3+R xC>$Bts9tMY<$5AP1Jit{eX )Af!d!S54j3ndEE~VRp&5+MXY)1LIk+j27xHJCG5~Z-1;eLnT34r_AY@4} z1E(+;F?pG8!85~@gt=*I#>qwi`+{UJ5avCRt11hA`(5mzY~ca%x3@7^(!?uINpf3tlmi}KJdp8uwct$)(|_gw+tbgW51fN(%mhW z$r-fTKh6`63NZS;fcv#(m5wp6KMH+ACNwr6JCrCY#1k9hKX1zrfeJ9HDj|SHVKw=h z_Xn+V`m>(9YK_yD)-To><62g}aitK;)DEw#TaChUzj2VbekjB~ynepcFq#qFF%Z3a zFNRJ!ulIHS;P9$zc3=1Fyg$|6u=I9p_Av25j0}OapNNo}+H1@71QTf)k>C1{RfGmv zSrt}KXoncs?!o3pI$h&DfG>t1PY_IhOXKcwCF{92R2`eG>fW`}dbf>`DM0vV7fXpK z#v?YH@3T_ni-UheCu&S|f5ws6qFJ{GmF$s ?JaL8jaGDVP zZ3%^`8Wifk9Fc`!WZs2N?VNn=&gYJIw(c`>Rmv3I;(f?p9o?PMed^J4^<#(|un;Z_ zLT=2&2A}gT1vWbe$L2((ECo)>_G7N|Pd7Z3D`Nw`b`x(NovnrXjsMOreI{LQ{tSMl zrF@Z6%DZ(t3!_q!F@ofUr%cKt%Q44F(PFwXL~lTg`J*EU)%S`A$sf|-GHmi7xW(P> zm~|}~=FzoB58 ~OZ*q+&niL_pD zb}urI-vBfc LYqouVJ&{0}&->^?6s9c Y>{P`8f=&0xepx`yQz;cxi$5$B_f+$V@RsC}U|BRGozvt_O<2KMa_a8C*<73tA?L zhQ95WE4shAa1IpPAE{6CSYR>CXLppT^$=6_QwR~IfIf#IEtzZa@N}KzVW1DOjQ{XT z8P|qE?UY_+0+ECN#DkZmg5<=hF)4+^_lL4EWh9wWt#z%SPDfR8M>!^GJ49H^)Y$h% z1^N7(?EYdG!>u9wrcL98)%*kYLg|P1PAz-WE*8&KB3^Sbu)6;BdgUT18Y?p4e7}Rn z;=0gLNl`cB!!1i^ECbIONT7%072JM6z5MSz9#}A|xBP*2p@_$e(PsN7Yul2f*-9vn z4hEKhx^s>vNa0fi{pFJA6?<4h- ;S9v;^}4HvE)L^^}+XDW#&8CqTL zNydV~u{oR36vm_pWBV5y7I94rCEO{Bm4kB*XXlr_chP+X3)Pul7j!9@vb-Ia3Z;TP z(#{y#Z9F2*s+#Kp%}#E@aW5K2hKK8m`YTkp6 RLOv zLh%^_+S)wVe!0rQ1CS=tL(3|qo+)ZHzlDv3y*DzLV_EHe0;&58PyJ<&&o7H{zji0E z^&6E|OcZV?7Iat?aFz`dj*s5tbbl*@FsW95_J+u^ij<-xEqB|(RrG0Gf!?mOk}P~P z=1E*In)Xi-qsH;0Sn}IE$GZK 2Swbbo*g4Khrs5czW)RkiyC$zd0o0zsMZra!8&dX9_eCLXBJI~ zPo6rfFzEw;4kGMw!cY}`G@9k^mHnz2?$@vxNI9taH%-sk*>SBK_}BO7RJf{NYiykU zby)OAG$DbLgocv5SS&ZZJz73bknT;V-i&8#do${f%7djOZYl3UR&&OOB3`Jo0V|PN zH8tZ4HbCNU#;RFh2e(w+yg4o0e>PcR* So)%;&%Zse%b!L9?1nriw` z`(KGXsb@;L`YgT*+0hch6vzGXce<0IS_9R6g=_P!(@RhU-}ooG2Cb3Z6LZRIASpzl zCFdHczG(kZ;J0vcaoOr1WUeM@ZfdJr`^s` A3j}(Mwe&R&v+ZH4-s+j|aFo zpRHTlR#AqgKhK&<+OZ|AgJ6yt^{;?aG+bY?oht7+A1?AZfJUd4pF}YTV8(Ozz@$vB zH`b30 t Z6>B S6n?K)8J1iNn8 z(rdc)zAGK7;7GGH&Es6vAi>;=c}>%e!~yKaR5kqMQ@+Y)(GRl(QgrL)O(Ct6#bei& z(V?Wm3ccZ(eP2AXmx!l-(}uTkNit)G6z4XK0Z26l6E0iUDl|=D?F3ve>D9<$D1{@V zgQu_9x(VX@`E?q0m7PJly0*F+j|mRmKRn6H%)%~^PHb9K9~ 2v0$#9Mau-?VnVMF=Gt!+d ze|Y=G^>X8^oP0lh%qj+6Npk!@s=^4s>yXizASjSPK&(;?VBva!{qI#BR|@qSg&{=- zlV3=ltwloFoh^h1@l2)j<+*hnUM-ek;3TD_TWf-XguX KR+zq z44lhPf<04HFE6_GHPia^D&4LA(&tQ=+t#i#=Jafo)J}55H7z;}ixaweG*ZtNnVVK4 zG#wVUh{ut|#B0hJV?2p+g={$*cAILEucSF$bJmi=P7T++xRB{Knt+1?pso*7Dc$@B zmwu`vxQQGpUs)RZiJ>0{iqT+&JZ7U{a(%72TbDCxIrsApO262E=nu0{2iUE;yQ~2P z+&{Yrr4rImGfinY(5X!yoAmg$81kBiC^EKqtGk2D_;DN4?>*w(qiqrDjJ8SqESEc< zE>!J)o!@ZD=&~+&t@cx_n^A@Qv`L;w=C)1*jbw&yjN&%UifeVQt5fIS4;~wp2g|ubf{^?$r{z# zEtRdE1V4Pte5@H?+#+|aCR##P#fZH-qokw+so8{358z$FBR9{Jj2M|pp! 2n%ab(tGh>|Ao<35{-vuwCRT(p#Xv_4 zfv s367LUrd|o%$X D;1MA{tfupUDh>XkdgZ{|-SPHdr+lD^?)+Vx&HM@qA(Yi&&BWMk7d zXPURnQrW=4!sYj`kT1}NkzHCJ?s0Pe*+bVWd3RV0;qC^dV4qv?eW*FhLoSB6m{ayb z0Dh>bp#j{lJUQs?RjX|%XRx#I(gnh^`VEZB{`mdHl9QB_>_MvySpG^F3`VaF3#H(( zvW@mv!GryC-_(-f+G nB?w{R;Vl8?%XVzEsK}ti0fK)sHrC*|12Otffsb+ zr(Kpw_os9>Q8YP0=-|m?f5vg6OGI)Lg`v^O5a%@ls78&^!S0kCDk9|lN*00oem143 z(-tP64wTxkmTcC0C{F&=1>F%fwVs5DWSJ5vnachW|4Fk-rsq2lhhnKAj|2$dfMpOd zf;$49XvE=nYRl9KXBq#ADKm>1hrgd%w6ET6y^5MBcufPFS%-z?8XHBl0O{)vJ@{s4 zTUWhoT2e~+B#HL6FAvdCa75_pKj7vaVS1`e*{lq;-12WD{^aF#13gfIGD=eHmq9UZ zgkt(%9+3cWji7ATdrTXZN#fpGQh#k4xO`;mkqFx4UZsWm8-BNh)av_xnXzyrO*E$4 zBKCj3A4m@APG3oQVf=W*{{Bv&IrzmvvfAC-8`igoavIp&7$f?z2?gtOE%Tkcq&+V| zI(^&hCI0g=IdaIj9GXdUiVc`P$**f#S+*LLc1UWOd(N|g@84ursl`&RWyV=z`Y16q z5E_8YI`bK~;a^|*r Z;fts5O0cX~w?TX&)F z9NHe2F5j2NQh)U0I?6YkiZ*Qu^K0`=6&DvbXI7s>IV6U-6!<=GrJU>lAiTCCjebe{ z)0%h?+^%8- zklgX!h=bDWCm9nBC!Oy~k)?`sYt$`U2nX18Y2ZGZvQCqx!*8_R4*p@hg2dzPY8c8L zuI*-`{4&+<9gg14S-XLWNgCLk-q<+nG5T04ACGkL7v (l2hvwbFTZC- z3&JF^b>&sU?25Yg7}{9N=c{Y?XBu?@>f|7yA^tS1Z)$&?t>(R|#!~iUaC$HwO@cD+ z$JJk@`&%iW#0G#0qf?bv%9IZ zsjakWRL=rAd7#)!3#;Da^f&-=mV1py(E9pz ?@uV1+=Goy;~E*t(XE4e4-{BN>hfri9BuNf@+oMewQLAisQ!z}!zYD)0BCqdS>p z`t{7@%U|buFoFb%kiiu@H?-X+8u&bIB`IJtj#`Zoqkup=(zoxV1jWRcBYj0d%~EBz z64{~y`WrmLFrBWEGs4hQ7oiDfA+@NEN4IYT03<8jNz7;P>sQ_Sd$N;X@A9V3)USM+ zjiMpEe+L}IsoeJY_4QLa#gMuQvzRoU7h*W4>;^$s(8M$`K7I;<4djHD@^Xl{gh}K} zOG^*&@ga=YR|PB{`H&R^8^SGbGA1TV5D+4j+gEsH0QiDUra*y?isw2?UDet6lUrAN z^+H;!g^Ry2Y?HDwhIM|caj)*1+n-+9O XVYqm}C@zxbh#|JG#OXvD@*LC9FK z1~C7nfL8#-zs#0}$c=gL-nvH+8QW`Z74TdMSN#>zyN5lvb^xB!aY%{qrcUfJC>?|0 zen|B^#8-g|UW@|AzWNx9LtvCZI_hcf|9gSV5EPgYJ+o2wsqN1gBwc`-;^#TBcL1DE zq}>F{>95MgpKJZmat5pQyYR_DN7*2f5$NK1jkIDj*+HW{`kB_&pc?nSV(_#HvoiNR zMM+ct>an5T-f!Ahn!XCL>GD|xLo0OlHONvo&7$Vf%y;Klmm&ed$xkPDP`>>dKqLJ@ zxsKwEg+1+=0llkga+pDBUWb_(KsG%h!Gj%cV3ogRGRWA~O%fa9uHsUmP8ayuAIku3 zZEY=%K}8S@Yrb=yYoDVsV0=JiEUB0?c; 7g(|a2Zo>i7Pzx4(Z zLsd(~4LOXZv+?V{svI>sZZom7kMfkEYn-h&GW00smptfX$dbrV=pp2MMT6Dsco!v- zc(Vk>3pbZ?GL!rEq{~_%DDX*r)NOm%1H*2TPRm!Zj*>hxp2;?J^|9@@^Di#wF 1l&pv#gc0LoF0l~5MUb`3mVAb6jx~_?1F3bJ zOtIhq@_syF-s+ev(#77#cV-5{w$wcm9*2;wAn--6M~ITZ5D_w;vWkkYAPk4v*C!nA zLF8W@PFOB2azp$ioQ6AOuz!P`NFBlx#&!T&f_RcN<3uPlDglm>E^^7f l~V>DN?Ip1fkUi4LX6s<06 WMOm&Ns_iol~O K$-j6lsQ*)Ce6^C<9kRzZ!j0$C^(-YWhWban-M!)T z;dgum!XMrL0Wu@mG`=1T(DcBG8JiTfsXm8hmjwK2NZZ-ysMu|!3(ONj60Ekq)lO}K zm#4G|QYlqnbe@=9`8G_B13T{HKtu^humsz|t oq#l$W$z2}yB|OPP79 zr7D&jtIrGyTfUUW?;E@UVEJ;yj8>)LtUxN^r&~osUbzYd)}?vhVcL?(slgJ3P(^fo zUmbjaOwLmuDVWDn$i8(Hgbu3(#rbc6AAHwQrqFFxwfTQW+&~$4b*}Y8#dF!TbZq^z zabZg)dy*Y;5!#XV20D~aB~kNng7a#~OM3$kgIINjS5$aTAi%n@ezxut4!MD#m+W&y zO^`Y)(7G^zEed*XgVKv7{;PzS&MasB5Js4kn0T8EPZ(+4Qdf*-^}IZ{Y&dDvV}YrB za7#{Gw-K)&!`v4^F=yv$I8E#j!&-rEqM?yus8d#0Xb## #NzS~lpg0X!wx{#)0-K#hPls56?LP2z2Z zd`gCb=}`WQk(qHGdL9eQd9=o_=tI7YoHn5IYj}FC=B M{W8TuA!jL z?UkU`xTAToVXrA!{eE(q=ulfurCXa `HL=in #Uf4y*bw#6z_ij| zX{j~y_yKwgz*J+CYY?GP)tswXK ^8Bc|TE|7dp@Oa- zCSNF*$Bv$rRo^5&Se&a&aN@oQ(k}4xB5osNUNas;q4eF8($L%GhC2?bz`T`zG}Hj6 zVSoGfvW{gkil_A6@njAsN`v($Ns0kyp@5ZLB1R;=yRUOIB_=9QfB^6Foj|#U(VFCS zE6{9*(?TTnSToF9n-zh0iD=>~@aUrUIDwa=7>FloS-?9g(C1+O!1%<(m`jExC?5Ox zICjjUFJP#J8IA+YgFmmhhce)+&aje{B6R!YpzJZQ0_$Sl!zC`BZ*#`rmMx>o_-zIY zcM8>0 GAV #uddV|z%z82AMmJ4maC z!98}QWC2%q=#6; pXXuOd^YXrrv{wm(gKJL5sPXDRqVz`VR~g&5>D~JA zPi6#VN+Jef+)rgw5(*T9r+iZ0t^Sg;sZKVv$0Ah)`$6yqPbx?{+YBYXCgwoIon_ OM_uz_dt8GM~>heTWT^$UgGR z%gsg3J-(fUhO$EU=IimvUt`Uu6QLTWml{*XEW;^0BaE?CM&o*wsnwHGjGVZF98XUC zus%mil*X*8$Ed<>uaNIb-}=DBhLUf?=h&AJ#>L#$fuO@cjQ8QGRuw^>uAV3;9?k~* zcgXHFA)#I5I&*+u+6!kiggdxn(auGN3OCr*(0F&~@e4ZpuXOOOBWRUg<1MHVDijR6 zD3aa`;HD5cfl0G4@sO8?0YCV<^r96P6iRIvwng1A-!dWz^j@UD>cAILs91U+4>Rf8 zsYF2h3n3^KKpI2(O;J{s{y?TULFgh}cMuox*)HBA(*w8-D)x085dVTpRZwn72$L7D zqqPL4Wi-EyJSZsKSY4}RjHBNj-f4S&uj#o(C6B{emht$SR$q?t`zm!thm!c#R>m}q zh^~eCYH|x)qeP6S9G7a%t*y`RU6uP_tNA1~MUC9y%F}tu%(!Y+J^v6lr)Fvsd^U?r zj@W{WQi*&g{(~b=PDEVGTjZmDo0 f O!oEi^Qe1mGya%XMHQ;1HI;>8z`xB_S1&trQ{AtK}TXs~`TY2k+rfJBYu&&jV zh0Jyre=0RY-xO))z#?)G3wVji`KWJ;awq+-7J$Lw9XNM?jKqik(Z=h0ni?#fUr%qN zi55OIA;9;u_MRl=+*GWee_#K*KXaBZuL`S2f0_L1w>TOwPfEtqFB!YXQTjuNeY9Jo zAm;Lt4?GW$_uMe{f*^}qG;bhx=!HbN1l3cYN#60=``-naDJGmL_#S(0TSzuxRr-kJ zejvik*bNrYiDnA={Nb cAEP$f8ZPI;6Xd!UVddeW JB{}UDgCC=FpxwE4YdLo3IBzFW#o1Q=>cE zxnJqdbSrKxQnPkp)RDVj>>ab|GlTk;$#=%XDpF-Vn(m9oA1}{qUX!-R7L+^)GV^7T zOi$*pQ<(WeG=p1cSkcRAJ^?o2V zx2-fHEBEf1kpS2cFXn^BVNjji?SEEqKW?RUJ6C%=8Z`n6N31&KNef2?3?6XhJSNkF zI3xAcIO0L`mlUtsgf<7bcU9j_;vA=0xzd!5B`;&)Fk}v@=g&K<5%@955*b$wWOFX( z{f_;3c*S&Q)~eqir)SklnHVFg$- 1g6jsT~HYe)GULSq( zVv+RD@Ge&UU8!nvq3m{f8dg1 o1(cO;}|X; z#HSLA=&-H0K^zB&6t(@4=uN$Kp9X5DUdgQFEY-I$v&nXdG-Uv)86=F5;<6+nkH-#X z1~4z$)F*7g_~@@+Wym7{KA S1PFBd48FJm8y+7aY_5Gc5o%8$mI@j~O&%4v`8qeqB zaev&$ZBVgADDJ0@ROO%|-`c)^zqL?q7 hu^L+M2zk`+_5EBkDvzj1NRUzW=)>$8u02gnf|uU-5)!fH9u72Lgd2XC zBLLTw>+O4U@RZjww0n^xA3W&)B#}E9m2TkuxHk&-3}1xT?aIAzZ~wmE>o)@_Nj &$C?=PX-fOSB8}^pgS>&Fa zjCF4my&TbE=X~+go|s8mxw`jHkKNTe63oM)ReUi%ND;P zlWG<#b&`VEckGm=rB2X<1llb?;K4e&yIzQc34vwKD~$co$MTBlALMu6H8*QDdWTEl z0LsYE2RanCAuTtgR9kE$+Aayehi2cQlLC^GlC76mwxRd|?L}5$la-Q;=C*^8`0w)X z12m-Yo?ma!h8VV>FRkyi5)$vt3rKn#>Tf|WX8K5!@x_DB*SP%1d2U@~!I}y9lkeU9 z{Kfj~8^6sh0`9((wWlPL(VA&^=yCI*l#w5r7qn$TP;PZoi9g#$15T5g8gH@k;wCN$LZy_kN})m_$akxUzlO-FnI+c z(tk`lY+xY%m5Ky#@UXj-1YqD06> gRJ`DIhzj(89KY%S`86QZ35yTs+e>I!i}B#Q`>%aWGE)L zI8 7fQ`B>7}Zgw2dDNU^x3&HC4ODNepy$Fekgi z1Ehf!uqQG@-Tt<@2-1u1Ev-VV1X)xzp5?Sw!2=Gn9DVPFgU#ExGhAuJszQgBri2d~ z*>#<3QHyf<;3CN@^=pcDQjIcQ@_8=04|yrl`?j$KKe(E( zvnNoPo&C;L8_)zFmFpj)+j0Dw3G|RSsAo?*NyeyUQXLlk%u{geEt5X9AW0Ilau4{# z{AH<2rSB(MHM+ a`g<0zFk;?vbX`Edjer*yj!!om8J57>`z5tDcYtFCq1EG9mYWRR_$gOt& z=D1xG--|e_=hIHILU{J7LXlJ>`;nSoLMq4TxD8IcCiFGBC;9pK2syP!I{nTE3(}{A zU6~!4=J%i1^#mTcwXLbG&Ge^Gr>$`=Z{q9QpC(j_rOt9S2JxO45=xW}XVnELy6>~* zzSP!Ywy`(f>poK6I~!*@(v?V1r;wl6YbT!&T*~V|W@|-Xbn3KO)F$%j0^!|z@3TKj zPdIsrz9#pU(Ue5^+S1EIFY+V{ ~^?Icx|3)qG>cdmzUkXj(NJv z_<~0obnup6KJ7Hz!Ya!{Ta%{(T2?6=xUzBY0 x+Ej37r!Pg)hOz3JAN>>IdgH{W9f8N!M&54G56NjUpAbHnZ!x&C z`1RLhb*tOhTo@gxHstVnM#>Lg`4?Duol0$i@{cDYVxr@M53Wqot^AQ@Y&2e0Jhb#v z;K}`o!YUsFufY47WrBl?enm_T7af*Y+$M*lmbT&p*4R rrCY)zceu)tt!czVQkCx3b^O^Dz `QsCA@#i_OjG1{_@|~7{i@(I&LPrlPtIM zP-)sQW5ufaWa@eDeN0?_kBXh2$%eF dZvk3B zIzs2M@zui6_40FP+s}2@R#p_#j%*@C;uzPhmE1D6Vz1J|GJDR_?`W;<5$Xhr8*@W- zD41NFoU%7y6aynzA2d~O*}cOILGt@`_&Z+?x_TrgCT2oLM~ANQ66{^TXQW_-phiE# zSqM!+f+V(klFpwH*_`Gp?0Fg-?0`WKaX6q412J;Y*Pfx9-UxLMSzW-ap0BDL_jH-& z@+%Ww_snm1wrIKgdGSj6sh$Pf$012=j(5MB+f0{q_5P~7wwH U3^cY`0VVUt5WEjV7P>F0~B<=%X@=CHxUocu^S=k!wy$P zY0%;!1_|Z9F@o=ejjq*3+Tok`?zGrXUW3$+I%&)mR|ejvt@{N N2|h0h;9pTDxz zy$If9f!>p3Ps8DZcEVyeD~BqB*1py6bvJ5QikqNhgf E0cpB7UXgfs}KS z7K;ZIvJ&QWPLWmD{xyUlsTXa^@&BCNk{WvtL9MVe#p!$D?qKR|zT^%o8eNliqeYGd z`Nb6(Krm8~cpiSdkoMZqqY0n~jwlHyr Y{JU!8@YL6VQjYRjorQq*H+Vd-z65;%6vP=xRn>;JovN^7{n1%73MWR(D1Mmj zq@>Io>Lc2pv(dMP_4j6EWhr2$OwFqy^2`#snzQo##b4gH+pL6)MMw4blrV8WVK_#2 zrwzE$5}R=45;c^}`AL)9S&MT=p7SKirak@GD&!un6m~>cg+Yq#ID4+Eg-fQy*$2;# zkjZ#lYE(rx$fej2^xM_-X`i^!2V~Ea=RV#4T38rsSE22h_ z@}1^@-ajin?$0T%)K>KLOyT3=Tk7ox$j)JJ-+(*MHr63{9-lGGyWd&j%-b!lh2|UK zU7?lU4~b*xXnJY4ndin?OSYuYq%mld*_oZH2XV5m(&sd5bz!vCWBD5}nr<|HFTTC1 z%|mM~S`adospwZj{K~mK)5tFfsIHp++s73WKWnDmUz|K`HACJws}@z!e5I&4vC`(t zg%}BoJOgL>+&T+E37ZU~-}Ocwj&~_ZF3-w-9}0-)3DZss=RINmdUqt-+01{W%w@Aa z9ZP8YC}O;q?91TuL{0E;P0z>gt}<+PiM1?tD=Dt~=|o9t73%y_v0m@7G#(f7p`7dV z&1HJs_4oE+X|nP?CkUNY=mY{)^Wr|U5-O1ibffBr9dnLn(2$sK&NpeSY7v^gu-Dbs z*Ei&Ehm!3wy4#+r+1~3O1!CNZTCk+xI8~3lnX^iOkfcRtwp{`HXA*sHM2@tUv9hl2 z0sb3X%A59)1r=12O#5t9`)nF`NK=PS-J@)7Ya70PlAM)4N>rSQt8iQ?!gQoAia+|j zWxh-2D bWl_F7WuhNK5+(P(Hd0ijPQyr^7f#CpJdPITxS7RlBRJYw{S9%kis_V NILrYuAw`&BaC&3y21yaqYf!Y}_P|1e^gxm4Q^*y==*V;AkwQ){M-=k9drg0q zylOWZ6Ph@wcR5dg_cpqOhOR5QhN8~-^jTk9Hu_3eHs4Y@S|7^WDpah*g N~v+>9O38%TD6Ny%OZKle9? sxAcp?)d(-AMwjs}!eld%zKnjbOn3lXU(`+x_$%zB%T( zu-G9dvd^CJNkmnGRZ^&IL2m(P7QY}!53NjAdP7QT@IxSc7MFkwVTKbM%*)L!EMh^B z2g8#R`W)UX^(J$j;K`1LqcPH0`?0Z>?UZ=SbY-+pS-r*B>x(55Itgta=i2Mm3S-Bf zgoqUEwLSl98>!Y_)N5@isYS_Og)fBnrB=R*OZW2^3}Hug!-6xkukePl(H!6odvB zW p4z#awuSh0)8{DN+OdKJ9BWB1s3j*thP6hZs8DIP<+a^SpN zo&lsW1H+VtiGi)aqdM~e4+AqoG8@H`*SD-5wB>K!EC&NY+SBE^6|?7r;KTT|^z;m! zaagt#&7w02Pry;kwK*6e4ePaB5#8v}Ow9K(l>5G2>2Uvf@%-64I!6<%3RR9xz1Ej_ zWYiN{s>bx|nI?ZY&r7Xh*|bo_KoOCW82 KV7vF0KSatk{uV8!V{O7W+W~R9pK{n83^TNo1 zRQo>Q0CfCfPh~u27yNiYH}wQHpgthX9~QIgfoPrewQCX*66gY95c>V?9ipp}8|W65 zgs`Se9(n4fN6_A~A;%R6PZC15(+~)*s4&iIj{x&LEUe)l!Lu(kWGnbEJ*CLe=#b (AAc!z0lp!xj%>z+=*xj)tVS?>rZ=2JMD5ifS~C$D~z5f zPzwpBG|SfSQA^HExLMzH$1ljoJ2%j%$Tyyr8VrnJf6kpbvi!ldoYavQ*$yaIobpW5 zOP)+V=d}@tcRp48IosZFFpq)UJA*LLI90vgh5+K)L~GMT>FlNrd- zu4pKdyt8;Zg2%3(^R4DFfs2t+9^D0PZDGu#vacT4c^7(tryBiXn}h3}@d5a6`5Ceu z`B=yMu?{n|LsC7io4HlT`F?DyE(s>oaWpnoUg_Wos4{i#%5X`3Z&lbesdGP}{@TJT zcYQR8NDt2T$sPLaQYSaHnTLtrl}t-(rYEIcoiR+F#`2BtCbi`AY{487lw=Pm_&JqI zdp|@mRHmM8w-QCm#?j+3Uqyj%<-S!N9 (SdL zyBtuniS_3cTnA !Ax(S`g5f&ai_dTN#(s`p$}yD zqa$`;-=`oG+p@h!G!-q@OgPU?{{8OHN6-*M9ga=@GPJd#XOm(mCbuT8z0TO-w{OdS zASwDXB*h+w`P9|UQCfppBS>0JAVlp-)_izz{FN~ z%MnfA*v*8*vv25Ku5GOPkowvnm|`TT`9_Rh^q_uBxI>r3@mF*l8imwyHwHE6F)&HG za7eI``3?<-?c?-w(&O96l=MXG5|}!4>r(S~s)UCRIpfni! @ zn}{3o_fjM@)jfWXC1c*y@}oy`@vA}EDUnVaeBB#EX3`PVOifveER;4c$g6_pPOWPR z@qSbdNh=m0zp!(fg3uqAkwFQ9sW-PD-0tx`iC3jOO;@_!w{W$>`ju;w(>+5$yGtuG z@2vecZ(6mAm1voC;TC ~7bc1VNfCo{-Cm5l?ITDf{8605MyTh#527 zEKQrXg>386_dEL-c5jDq$fyj91m>`}2ARp+a!L*RVoxtMtlSVyKF=EVO<%(NlL#l< zi=gP|_bX2Bqx5e*WyB`9^UE(Oh{EkVwaoIiyWt$6YClmrGG+r+DxAJ#EWhoF- him*o>|*=7LJ?f+&AK)3Ovpz;qpeNZmIHBYwg2(`s$SS zEy+i<8Lsp{*i9({rYzhE1pRMTXdQ)~Wb;3no5==LcxYcc+i+4-Y!t@=`V>(QoypCr zZ#Eovw3kY^!{8HTjJc3rJcL|)eZRty63+CbS_nTv<(HaiuUjsp7(Fs3&CY%cIy0zL zpLLv>gC; F)`^|)&Ko>uWCBuR{R{8 zks7rHVgBtWHffQ0It)fLtELCw9m#y(@MLY=F_+`z<;96cf6vOsrl-`^zOm%_&IRZA z*nDS7e<-a`lVITdbW~=5P;;LNjOUqhry_fKkI3YW-LC-lK`RIJ3#SC3E$nMfE5*rR z$wB>qOgfmQ2_O02i{q}0UK4KY_x*gH8tGGOH;QT(NgoN8?7vrYuWEa^Tt|v!>x~4T z<*LK~fqp@kN1{aeN64(K4sG)$<3wTs`~LHfsBIL)_mMad4~cA4K30OgHhJwtEHq^_ z>$rh2hGJo(Ah?__d5%7;`=~t`2hp0>nVFx>g%Ex$Z!!hFBSiCknTgB~rK`mbYPd!& z6T$|Z63DWAnUlH?EG9%VE*M>IY~YCB)-i-`P$x3pt}Z2VB(hIbxRjT5|23pk1Pq8U zGovP+1SH6d+(el_UNzFp+NSqEK@t2W6;cnCpCr3+g|hCxRq@>KYGxDRU;6iwajIqw z?z$u4@o--(#qdel&EC!i%M@!a#HYmSL2QBvB)~+R5<#P*5Wf2B`F2E36Zk0N=UFz9 zKJH;{_&|}mt>Yj0r&5^pzW%FFe*pHsY)v^ce2cO-C7}yX_FUojN~2HgC}sRu6O_n7 z(D=K6$YfNG!GX{AJR#7t3SkWz@SR_tggie4g{es1^9hGJMV}6-%d;c={=Fz~j_%~D z%jvsYXO5JP-LGtUAOJ=u>{LFI%Z?#>w2X`wcxaK;a{KM37hpNM%g;}ETL93ze=nI( z=kEQx&ur&=aog*Rac3l!(hwv?3aFRchJO5jd0D-`zNG9fB7hN>v;i Yl+TU;7wXl42aNiy#H&f-a$t{rfzeV{PW>ud)sZi&}`~nse zp5n*o;MFNbS0A{g6;I^v<(o%&j+#&LS}Nc-bJgb+m1U?!rBaRT-Wq&<{GOo?$s3R= zf2Le31Kz(I+#*8k-oIDLhr~AMt+dnoRzEbVFLnqCnjxO@zZdyi^fY?xmNE*zG4#Hl z^AeP*YcX)B*|7c1tZRqWEv$;MtKFlE#}-TXui=e!ShKBPrKevuNxZtTVB|Je=#0?w zpVyMn8Wjzu<;b~!J|^^IOG=2mLz_QhsYv^hn^xAg4m*>FOFhp{7Fr+Dvf&tOSyz0N zBovk8gD!`4ei!7+iJ6l>;a>P$VXpk?xZZ)Uc;zp-*}4pS?_Wm1?x}@YC&k@&GP&0y zc=;pJb8> 9vV z+I|CoH%iiNSG$m)8cssGw_<>gE+C=Q{J5b=7BoN8y8SyFjvYHTzqV`AdiuYx0AeTm znw1N9%YXc^h*47@KPcsvJ+f?WcOulsPPiY@-m=UV(-bJkLV%yK1pvI()_G`E!S`qq zV8U7HQK;ObUFrpZ(C(`2%l`>(xO)Vy3(b`BzWp$RAnIqoAfimhbjM~7QHv?v+CIU0 z6Aq;vHlmr_NO+mDa_{ufoKdv>DYvBkE<^0)Z- BmFQ!smnVZd z%zFm z4Pt>EG0rhV^9SRtLjsp%oY(@$Al|iraIuaK+0dUb9LK>BtfV5*-&=veG%QO&Z{tN< za0Ehy4X{c$L+n!+phCYMVyb^VSo?v*yjbFt+LxAu(bL!Q1#xk4#P7ga3+#)RQMJ>8 zQj9rg{?M9D;%!x&$jMS0r&C$qN$F0^1-{)9Fw}o|gU(`L 6e=fG4Q~A+#_u0jXcn4+j8oL9k+TKz&FQ#myIjwjV zoyf-CO}apR Nb85FF*RYTZ*M!y$@xmmtGSl>-02n#YQ?#Apm9k*A`PgA)0Zhj1K3n>&tn^E}UxS*42mE&WH=PuR2t35lfy_$5 zXVT)B%S@N_`QsI)Bu6gD@EAb8;TjNqsYo!RlFp4kHf3I?Rs}Uu*!l9Z4xtrCO+K?9 z_q%$_Mw(#G6CBa;wM#`eq107VBft2Ee>t>mV^~t?O9{#2%@4_(s(zq5VjF~9($Sy? zC=f17OCG(f%MWj>|B#ecZNabu7r+rdHmn)Kt#zAkQk|)IesL>ZbEveMO4bu_-nc?i zLD)-7Ou?)7i F=uqrXdUVfWY+uktOlbi%vp=>;hGo`k*U7!(534A6A{!kZ@bz0p z6%fdh_e;(TEbysQOTJv1aZkq0LYw_u**yt+^>cG t~(~3lFCz0UK6e-!~ O3D 9xDMgL=W( zbmHrHzWF`VnJIIIaNlXWpcF;wZKt=Iq-t-LDy>dM17`(=sEkdl=Ze0&f~GIy-Zv zKRsIZ0VQ{TJ%|YH$vGQ#k29!igsO7jhTlGWRs<6aq#k(5;bYpG ? z@sM{{k>4fVCG{_80u_4ULfEt`Dy2>_sj^J&n_JT1d*