Introduction

This is the R notebook for the analyses in the paper:

Init

options(digits = 2)
library(pacman)
p_load(kirkegaard, readxl, rms)

Data

Data were collected from a variety of French and Dutch-language documents, both public and privately in nature. The data covers various subdivisions of Belgium (provinces, arrondissements and municipalities/communes/gemeenten), including districts of cities.

#Brussels
brussels = readxl::read_xlsx("data/Brussels.xlsx") %>% 
  df_legalize_names()

#legalize names in data frame, keep labels
brussels %<>% 
  #mutate
  mutate(
    #as fraction
    North_African_proportion_2016_pct = North_African_proportion_2016_pct/100
  )

#municipals
municipals_muslims = readxl::read_xls("data/Moslims-Marokkaans-Turks-Ander-gemeente-2011.xls", skip = 1) %>% 
  df_legalize_names() %>% 
  mutate(NIS_code = as.numeric(NIS_code))
## Warning in evalq(as.numeric(NIS_code), <environment>): NAs introduced by
## coercion
municipals_income = readxl::read_xls("data/fisc2009_c_nl_tcm325-151868.xls", sheet = "per indiv. aang.", skip = 6) %>% 
  df_legalize_names()

municipals_age = readxl::read_xls("data/Bevolking-leeftijd-gemeente-2014-65plus.xls", skip = 2) %>% 
  df_legalize_names() %>% 
  rename(
    NIS_code = NIS_CODE
  ) %>% 
  mutate(
    NIS_code = NIS_code %>% as.numeric()
  )
## Warning in function_list[[k]](value): NAs introduced by coercion
municipals_density = readxl::read_xls("data/density.xls", sheet = 2) %>% 
  df_legalize_names() %>% 
  filter(!is.na(NIS_code)) %>% 
    mutate(
    NIS_code = NIS_code %>% as.numeric()
  )

#merge by NIS
municipals = left_join(municipals_muslims, municipals_income, by = "NIS_code") %>% 
  left_join(municipals_age, by = "NIS_code") %>% 
  left_join(municipals_density, by = "NIS_code") %>% 
  #mutate
  mutate(
    Muslim_pct = Totaal_moslims,
    median_income = Mediaan_inkomen_per_aangifte,
    mean_income = Gemiddeld_inkomen_per_aangifte,
    income_IQR_ratio = Interkwartiele_co_ffici_nt,
    name = Gemeente,
    population = Totale_bevolking,
    density = Bevolkings_dichtheid_inw_km,
    #age vars
    age_65ao_pct = pct65_Tot,
    age_0_19_pct = (Totale_bevolking_Minder_dan_5_jaar + Van_5_tot_9_jaar_2 + Van_10_tot_14_jaar_2 + Van_15_tot_19_jaar_2)/Tot
  )
## Warning: Column `NIS_code` has different attributes on LHS and RHS of join

Brussels

S factor

Social data usually correlate strongly, meaning that (in latent variable theory) they measure a common underlying construct. This we call S, for general socioeconomic factor. It can be extracted with factor analysis.

#social data
brussels_social = brussels[-c(1:3)]

#impute single missing datapoint
brussels_social %<>% miss_impute()
## Warning in predict.lm(object, newdata, se.fit, scale = 1, type =
## ifelse(type == : prediction from a rank-deficient fit may be misleading
#analyse
S_fa = silence(fa(brussels_social))
S_fa_wtd = silence(fa(brussels_social, weight = brussels$Total_population_of_the_municipality_2014 %>% sqrt))

#there are often warnings due to the small number of cases. In most cases, they are not problematic for the analysis.
S_fa_methods = fa_all_methods(brussels_social, messages = F)

S_fa_methods$loadings %>% wtd.cors()
##        minres pa minchi
## minres      1  1      1
## pa          1  1      1
## minchi      1  1      1
#no variation in loadings by method

S_fa_methods$scores %>% wtd.cors()
##                   regression_minres Thurstone_minres Bartlett_minres
## regression_minres              1.00             1.00            0.99
## Thurstone_minres               1.00             1.00            0.99
## Bartlett_minres                0.99             0.99            1.00
## regression_pa                  1.00             1.00            0.99
## Thurstone_pa                   1.00             1.00            0.99
## Bartlett_pa                    0.99             0.99            1.00
## regression_minchi              0.99             0.99            1.00
## Thurstone_minchi               0.99             0.99            1.00
## tenBerge_minchi                0.99             0.99            1.00
## Bartlett_minchi                0.99             0.99            1.00
##                   regression_pa Thurstone_pa Bartlett_pa regression_minchi
## regression_minres          1.00         1.00        0.99              0.99
## Thurstone_minres           1.00         1.00        0.99              0.99
## Bartlett_minres            0.99         0.99        1.00              1.00
## regression_pa              1.00         1.00        0.99              0.99
## Thurstone_pa               1.00         1.00        0.99              0.99
## Bartlett_pa                0.99         0.99        1.00              1.00
## regression_minchi          0.99         0.99        1.00              1.00
## Thurstone_minchi           0.99         0.99        1.00              1.00
## tenBerge_minchi            0.99         0.99        1.00              1.00
## Bartlett_minchi            0.99         0.99        1.00              1.00
##                   Thurstone_minchi tenBerge_minchi Bartlett_minchi
## regression_minres             0.99            0.99            0.99
## Thurstone_minres              0.99            0.99            0.99
## Bartlett_minres               1.00            1.00            1.00
## regression_pa                 0.99            0.99            0.99
## Thurstone_pa                  0.99            0.99            0.99
## Bartlett_pa                   1.00            1.00            1.00
## regression_minchi             1.00            1.00            1.00
## Thurstone_minchi              1.00            1.00            1.00
## tenBerge_minchi               1.00            1.00            1.00
## Bartlett_minchi               1.00            1.00            1.00
#no variation in scores by method

#loadings
fa_plot_loadings(list(standard = S_fa,
                      weights = S_fa_wtd))

GG_save("figs/Brussels_loadings.png")
#extreme results!

#save scores
brussels$S = S_fa$scores[, 1] * -1

Demographics

Examine relationship to North-African demographics.

#specific S variables
wtd.cors(cbind(brussels$North_African_proportion_2016_pct, brussels_social))[1, -1]
##                                                         Wealth_index_2013 
##                                                                     -0.95 
##                            Mean_wages_compared_with_the_national_mean_pct 
##                                                                     -0.95 
##                 Share_of_households_in_demand_for_social_housing_2012_pct 
##                                                                      0.86 
##                                                Unemployment_rate_2012_pct 
##                                                                      0.96 
##                                                  Median_income_2012_euros 
##                                                                     -0.90 
##          _18_64_years_beneficiaries_of_social_integration_income_2012_pct 
##                                                                      0.89 
##           Juvenile_delinquency_in_the_total_population_under_18_years_pct 
##                                                                      0.91 
##                       Delinquency_in_the_total_population_18_25_years_pct 
##                                                                      0.92 
##                                   Proportion_with_a_university_degree_pct 
##                                                                     -0.80 
##                                            Part_with_low_Birth_weight_pct 
##                                                                      0.44 
##                                                  Life_expectancy_at_birth 
##                                                                     -0.75 
##                    Children_born_in_a_family_with_no_income_from_work_pct 
##                                                                      0.95 
##                                                     Population_growth_pct 
##                                                                      0.81 
##                                                 Houses_median_prize_euros 
##                                                                     -0.67 
##                                                Fraudulent_declaration_pct 
##                                                                      0.67 
##                Rate_of_borrowers_with_at_least_one_non_regularized_credit 
##                                                                      0.90 
##                                            Mothers_20_years_pct_2009_2013 
##                                                                      0.81 
##                               Fetal_infantile_mortality_on_1000_2009_2013 
##                                                                      0.65 
##                        Proportion_with_two_years_late_in_school_2013_2014 
##                                                                      0.92 
## _2013_2014_Special_humanity_designed_for_young_with_mental_deficiency_pct 
##                                                                      0.84 
##       Relative_risk_of_mortality_by_cardiovascular_diseases_men_2009_2013 
##                                                                      0.62 
##            General_humanities_technical_and_professional_humanities_ratio 
##                                                                     -0.85
#scatter
GG_scatter(brussels, "North_African_proportion_2016_pct", "S", case_names = "Name", repel_names = T) +
  scale_x_continuous(labels = scales::percent)