Time preferences and national IQ

About

This is the R notebook for the paper:

• Kirkegaard, E. O. W. & Karlin, A. (2020). National intelligence is more important for explaining country well-being than time preference and other measured non-cognitive traits. Mankind Quarterly.

Init

#options
options(
  digits = 3
  )

#packages
#you have to use old broom version
#devtools::install_version("broom", "0.5.6")
library(pacman)
p_load(kirkegaard,
       readxl,
       rms,
       haven,
       sf,
       tmap
       )

#theme
theme_set(theme_bw())

Functions

#describe without bug, remove the custom class
describe = function(...) {
  y = psych::describe(...)
  class(y) = "data.frame"
  y
}

#drop unused factor levels
df_drop_fct_levels = function(x, drop_unused = T) {
  #find factors
  x_fct = map_lgl(x, is.factor)
  
  #loop and drop unused
  for (v in names(x[x_fct])) {
    #unused
    if (drop_unused) {
      x[[v]] = x[[v]] %>% fct_drop()
    }
  }
  
  x
}

#test it
assert_that(iris[1:5, ] %>% pull(Species) %>% levels() %>% equals(levels(iris$Species)) %>% all())

## [1] TRUE

assert_that(iris[1:5, ] %>% df_drop_fct_levels() %>% pull(Species) %>% levels() %>% equals("setosa"))

## [1] TRUE

#spatial lag function
#make spatial lag variables
#slow because it recalculates distances every time
#we cant just precompute because the function cannot handle missing data properly
add_spatial_lag = function(x, vars) {
  #some temp var to merge by
  if (nrow(x) == 0) stop("0 rows in `x`", call. = F)
  x[["..tmp"]] = 1:nrow(x)
  
  for (v in vars) {
    #no NA
    tmp = x[!is.na(x[[v]]) & !is.na(x$lat) & !is.na(x$lon), ] %>% st_drop_geometry()
    tmp_var = str_glue("{v}_lag") %>% as.character()
    
    #estimate
    tmp[[tmp_var]] = spatialstatstools::SAC_knsnr(tmp, v, k = 3)
    
    #join
    x[[tmp_var]] = NULL
    x = left_join(x, tmp %>% select(`..tmp`, !!tmp_var), by = "..tmp")
  }
  
  #remove temp var
  x[["..tmp"]] = NULL
  
  x
}

Data

#Wang paper
#https://www.sciencedirect.com/science/article/abs/pii/S0167487015001439
timepref_wang = read_xlsx("data/time_pref.xlsx") %>% 
  df_legalize_names() %>% 
  mutate(ISO = pu_translate(Country))

## No exact match: Bosnia & Herzigovina

## Best fuzzy match found: Bosnia & Herzigovina -> Bosnia & Herzegovina with distance 1.00

#Falk
#https://academic.oup.com/qje/article/133/4/1645/5025666
timepref_falk = read_dta("data/GPS_Dataset/GPS_dataset_country_level/country_v11.dta") %>% 
  df_legalize_names() %>% 
  rename(ISO = isocode) %>% 
  mutate(ISO = str_to_upper(ISO))

#Rieger meta-analysis of time pref
timepref_rieger = readxl::read_xlsx("data/Rieger et al 2020.xlsx") %>% 
  df_legalize_names() %>% 
  #remove aggregates
  filter(Aggregate == 0) %>% 
  mutate(
    ISO = pu_translate(Unit)
  )

#data from cog elite study
#https://rpubs.com/EmilOWK/cognitive_elite_2019
cogelite = read_rds("data/cog_elite_data_out.rds")

#remove old SPI data
cogelite %<>% select(-(SocialProgressIndex:Numberofgloballyrankeduniversities0none550), -SPI_g)


#Becker IQ data
#1.3.3
becker = read_excel("data/NIQ-DATASET-V1.3.3/NIQ-DATA (V1.3.3).xlsx", sheet = "FAV", skip = 1) %>% 
  df_legalize_names()

## New names:
## * `` -> ...1
## * `` -> ...2

names(becker)[1:2] = c("ISO", "Country")
becker %<>% mutate(
  IQb = QNWplusSASplusGEO,
  IQr = R
)

#World bank IQs
worldbank = read_excel("data/hcidataoctober2018.xlsx", sheet = 4) %>% 
  mutate(
    ISO = pu_translate(`Country Name`),
    IQwb = `HUMAN CAPITAL INDEX`
  )

## No exact match: eSwatini

## Best fuzzy match found: eSwatini -> Eswatini with distance 2.00

#SPI data
SPI = read_excel("data/2014-2019-SPI-Public.xlsx", sheet = "2019") %>% 
  #remove stuff in parens first
  #then clean names normally
  set_colnames(names(.) %>% str_replace_all("\\(.*\\)", "")) %>% 
  df_legalize_names() %>% 
  rename(ISO = Code)
SPI_indicators = SPI %>% select(Undernourishment:Percent_of_tertiary_students_enrolled_in_globally_ranked_universities) %>% names()

#remove world average
SPI = SPI[-1, ]

#ecnomics data
#world bank data
#https://data.worldbank.org/indicator/NY.GNP.PCAP.PP.CD
#https://data.worldbank.org/indicator/NY.GDP.PCAP.PP.CD
#these have missing data, which we can impute
GNI = read_excel("data/API_NY.GNP.PCAP.PP.CD_DS2_en_excel_v2_1218005.xls", skip = 3)
GDP = read_excel("data/API_NY.GDP.PCAP.PP.CD_DS2_en_excel_v2_1217845.xls", skip = 3)

#save back, add growth
GNI %<>% mutate(
  GNI_growth = (`2019` / `1990`) %>% log10(),
  GNI_2019 = `2019` %>% log10(),
  ISO = `Country Code`
)

GDP %<>% mutate(
  GDP_growth = (`2019` / `1990`) %>% log10(),
  GDP_2019 = `2019` %>% log10(),
  ISO = `Country Code`
)

#these are based on self-report data gathered by Gallup over the years
#https://news.gallup.com/poll/166211/worldwide-median-household-income-000.aspx
median_income = read_csv("data/median income gallup.csv") %>% 
  mutate(
    median_income = medianPerCapitaIncome %>% log10(),
    median_income_hh = medianHouseholdIncome %>% log10(),
    ISO = pu_translate(country)
  )

## 
## ── Column specification ────────────────────────────────────────────────────────
## cols(
##   country = col_character(),
##   medianHouseholdIncome = col_double(),
##   medianPerCapitaIncome = col_double(),
##   pop2020 = col_double()
## )

#economic outcomes
econ_outcomes = c("median_income", "median_income_hh", "GNI_2019", "GNI_growth", "GDP_2019", "GDP_growth")

#merge
d = full_join(timepref_wang, cogelite, by = "ISO") %>% 
  full_join(timepref_falk, by = "ISO") %>% 
  full_join(timepref_rieger, by = "ISO") %>% 
  left_join(becker %>% select(ISO, IQb, IQr), by = "ISO") %>% 
  left_join(worldbank %>% select(ISO, IQwb), by = "ISO") %>% 
  left_join(median_income %>% select(ISO, median_income, median_income_hh), by = "ISO") %>% 
  left_join(GNI %>% select(ISO, GNI_2019, GNI_growth), by = "ISO") %>% 
  left_join(GDP %>% select(ISO, GDP_2019, GDP_growth), by = "ISO")

assert_that(!any(duplicated(d$ISO)))

## [1] TRUE

#Hong Kong data
d$IQ[d$ISO == "HKG"] = 105.7
d$UN_macroregion[d$ISO == "HKG"] = "Eastern Asia"
d$Names[d$ISO == "HKG"] = "Hong Kong"
d$lat[d$ISO == "HKG"] = 22.321751
d$lon[d$ISO == "HKG"] = 114.169663

#spatial data
world = sf::read_sf("data/TM_WORLD_BORDERS-0.3/TM_WORLD_BORDERS-0.3.shp") %>% 
  rename(ISO = ISO3)

Recode & spatial

#fix variables
d %<>% mutate(
  #rename
  time_pref_wang = Choose_to_wait,
  time_pref_falk = patience,
  time_pref_rieger10 = UP10,
  time_pref_rieger6 = UP6,
  HDI2018 = Human_Development_Index_HDI,
  IQlv = IQ,
  
  #standard error of proportion (approx)
  time_pref_wang_SE = sqrt((time_pref_wang * (1 - time_pref_wang)) / N),
  time_pref_wang_RSE = 1 / time_pref_wang_SE,
  
  popsqrt = population2017 %>% sqrt(),
  
  #region coding
  UN_macroregion = plyr::revalue(UN_region, c("Australia and New Zealand" = "N & W Europe + offshoots",
                          "Europe" = "N & W Europe + offshoots",
                          "Northern Europe" = "N & W Europe + offshoots",
                          "Eastern Africa" = "Africa",
                          "Middle Africa" = "Africa",
                          "Western Africa" = "Africa",
                          "Southern Africa" = "Africa",
                          "Europe" = "North & Western Europe",
                          "Northern America" = "N & W Europe + offshoots",
                          "South America" = "Latin America",
                          "Central America" = "Latin America",
                          "Western Asia" = "MENA",
                          "Northern Africa" = "MENA",
                          "Western Europe" = "N & W Europe + offshoots"
                          ), warn_missing = F)
)

#standardize to this subset
d$time_pref_wang_z = d$time_pref_wang %>% standardize()
d$time_pref_falk_z = d$time_pref_falk %>% standardize()
d$time_pref_rieger10_z = d$time_pref_rieger10 %>% standardize()
d$time_pref_rieger6_z = d$time_pref_rieger6 %>% standardize()

SPI

We impute first, and then split indicators into groups.

#impute with at most 30% missing data
SPI_imputed = SPI %>% 
  #indicators + ISO
  select(ISO, !!SPI_indicators) %>% 
  #impute
  {
    #max NA per case
    (max_na = floor(ncol(.[-1]) / 3))
    
    #impute
    #remove the non-data ISO variable and then add back
    ISO = .$ISO
    cbind(ISO = ISO, miss_impute(.[-1], max_na = max_na))
  } %>% 
  as_tibble()

#looks ok? remove the cases with missing data still
SPI %>% select(!!SPI_indicators) %>% miss_amount()

## cases with missing data  vars with missing data cells with missing data 
##                   0.670                   1.000                   0.206

SPI_imputed %>% miss_amount()

## cases with missing data  vars with missing data cells with missing data 
##                   0.260                   0.980                   0.178

SPI_imputed %>% miss_by_var()

##                                                                   ISO 
##                                                                     0 
##                                                      Undernourishment 
##                                                                    55 
##                                               Maternal_mortality_rate 
##                                                                    42 
##                                                  Child_mortality_rate 
##                                                                    42 
##                                                        Child_stunting 
##                                                                    42 
##                                       Deaths_from_infectious_diseases 
##                                                                    42 
##                               Access_to_at_least_basic_drinking_water 
##                                                                     9 
##                                                 Access_to_piped_water 
##                                                                    10 
##                        Access_to_at_least_basic_sanitation_facilities 
##                                                                    12 
##                                                 Rural_open_defecation 
##                                                                    51 
##                                                 Access_to_electricity 
##                                                                    24 
##                                         Quality_of_electricity_supply 
##                                                                    60 
##                           Household_air_pollution_attributable_deaths 
##                                                                    42 
##                      Access_to_clean_fuels_and_technology_for_cooking 
##                                                                    45 
##                                                         Homicide_rate 
##                                                                    23 
##                                                 Perceived_criminality 
##                                                                    59 
##                                        Political_killings_and_torture 
##                                                                    58 
##                                                        Traffic_deaths 
##                                                                    42 
##                                                   Adult_literacy_rate 
##                                                                    50 
##                                             Primary_school_enrollment 
##                                                                    40 
##                                           Secondary_school_enrollment 
##                                                                    43 
##                                 Gender_parity_in_secondary_enrollment 
##                                                                    40 
##                                           Access_to_quality_education 
##                                                                    59 
##                                        Mobile_telephone_subscriptions 
##                                                                    27 
##                                                        Internet_users 
##                                                                    22 
##                                           Access_to_online_governance 
##                                                                    44 
##                                                      Media_censorship 
##                                                                    59 
##                                                 Life_expectancy_at_60 
##                                                                    42 
##                       Premature_deaths_from_non_communicable_diseases 
##                                                                    42 
##                                          Access_to_essential_services 
##                                                                    42 
##                                          Access_to_quality_healthcare 
##                                                                    59 
##                             Outdoor_air_pollution_attributable_deaths 
##                                                                    42 
##                                              Greenhouse_gas_emissions 
##                                                                    45 
##                                                      Biome_protection 
##                                                                     6 
##                                                      Political_rights 
##                                                                    40 
##                                                 Freedom_of_expression 
##                                                                    59 
##                                                   Freedom_of_religion 
##                                                                    59 
##                                                     Access_to_justice 
##                                                                    59 
##                                             Property_rights_for_women 
##                                                                    59 
##                                                 Vulnerable_employment 
##                                                                    50 
##                                                        Early_marriage 
##                                                                    55 
##                                    Satisfied_demand_for_contraception 
##                                                                    48 
##                                                            Corruption 
##                                                                    55 
##                                       Acceptance_of_gays_and_lesbians 
##                                                                    57 
##                        Discrimination_and_violence_against_minorities 
##                                                                    58 
##                                 Equality_of_political_power_by_gender 
##                                                                    59 
##                 Equality_of_political_power_by_socioeconomic_position 
##                                                                    59 
##                           Equality_of_political_power_by_social_group 
##                                                                    59 
##                                           Years_of_tertiary_schooling 
##                                                                    47 
##                                       Women_s_average_years_in_school 
##                                                                    56 
##                                          Globally_ranked_universities 
##                                                                     1 
## Percent_of_tertiary_students_enrolled_in_globally_ranked_universities 
##                                                                     1

SPI_imputed %<>% miss_filter()
nrow(SPI_imputed)

## [1] 176

SPI_imputed %>% miss_amount()

## cases with missing data  vars with missing data cells with missing data 
##                       0                       0                       0

#factor analyze
#for total set, and each broad subset
(SPI_g = fa(SPI_imputed %>% select(!!SPI_indicators)))

## Factor Analysis using method =  minres
## Call: fa(r = SPI_imputed %>% select(!!SPI_indicators))
## Standardized loadings (pattern matrix) based upon correlation matrix
##                                                                         MR1
## Undernourishment                                                      -0.80
## Maternal_mortality_rate                                               -0.83
## Child_mortality_rate                                                  -0.88
## Child_stunting                                                        -0.88
## Deaths_from_infectious_diseases                                       -0.79
## Access_to_at_least_basic_drinking_water                                0.87
## Access_to_piped_water                                                  0.85
## Access_to_at_least_basic_sanitation_facilities                         0.85
## Rural_open_defecation                                                 -0.64
## Access_to_electricity                                                  0.81
## Quality_of_electricity_supply                                          0.87
## Household_air_pollution_attributable_deaths                           -0.84
## Access_to_clean_fuels_and_technology_for_cooking                       0.85
## Homicide_rate                                                         -0.16
## Perceived_criminality                                                 -0.59
## Political_killings_and_torture                                         0.65
## Traffic_deaths                                                        -0.68
## Adult_literacy_rate                                                    0.78
## Primary_school_enrollment                                              0.66
## Secondary_school_enrollment                                            0.91
## Gender_parity_in_secondary_enrollment                                  0.54
## Access_to_quality_education                                            0.75
## Mobile_telephone_subscriptions                                         0.66
## Internet_users                                                         0.91
## Access_to_online_governance                                            0.77
## Media_censorship                                                       0.45
## Life_expectancy_at_60                                                  0.84
## Premature_deaths_from_non_communicable_diseases                       -0.61
## Access_to_essential_services                                           0.95
## Access_to_quality_healthcare                                           0.83
## Outdoor_air_pollution_attributable_deaths                             -0.19
## Greenhouse_gas_emissions                                              -0.45
## Biome_protection                                                       0.24
## Political_rights                                                       0.62
## Freedom_of_expression                                                  0.46
## Freedom_of_religion                                                    0.32
## Access_to_justice                                                      0.67
## Property_rights_for_women                                              0.68
## Vulnerable_employment                                                 -0.86
## Early_marriage                                                        -0.74
## Satisfied_demand_for_contraception                                     0.66
## Corruption                                                             0.80
## Acceptance_of_gays_and_lesbians                                        0.63
## Discrimination_and_violence_against_minorities                        -0.47
## Equality_of_political_power_by_gender                                  0.51
## Equality_of_political_power_by_socioeconomic_position                  0.48
## Equality_of_political_power_by_social_group                            0.47
## Years_of_tertiary_schooling                                            0.80
## Women_s_average_years_in_school                                        0.87
## Globally_ranked_universities                                           0.30
## Percent_of_tertiary_students_enrolled_in_globally_ranked_universities  0.62
##                                                                          h2
## Undernourishment                                                      0.648
## Maternal_mortality_rate                                               0.696
## Child_mortality_rate                                                  0.776
## Child_stunting                                                        0.774
## Deaths_from_infectious_diseases                                       0.628
## Access_to_at_least_basic_drinking_water                               0.748
## Access_to_piped_water                                                 0.725
## Access_to_at_least_basic_sanitation_facilities                        0.717
## Rural_open_defecation                                                 0.410
## Access_to_electricity                                                 0.654
## Quality_of_electricity_supply                                         0.759
## Household_air_pollution_attributable_deaths                           0.698
## Access_to_clean_fuels_and_technology_for_cooking                      0.724
## Homicide_rate                                                         0.026
## Perceived_criminality                                                 0.350
## Political_killings_and_torture                                        0.426
## Traffic_deaths                                                        0.463
## Adult_literacy_rate                                                   0.603
## Primary_school_enrollment                                             0.435
## Secondary_school_enrollment                                           0.825
## Gender_parity_in_secondary_enrollment                                 0.295
## Access_to_quality_education                                           0.569
## Mobile_telephone_subscriptions                                        0.441
## Internet_users                                                        0.824
## Access_to_online_governance                                           0.594
## Media_censorship                                                      0.204
## Life_expectancy_at_60                                                 0.710
## Premature_deaths_from_non_communicable_diseases                       0.372
## Access_to_essential_services                                          0.904
## Access_to_quality_healthcare                                          0.694
## Outdoor_air_pollution_attributable_deaths                             0.035
## Greenhouse_gas_emissions                                              0.199
## Biome_protection                                                      0.060
## Political_rights                                                      0.380
## Freedom_of_expression                                                 0.215
## Freedom_of_religion                                                   0.105
## Access_to_justice                                                     0.449
## Property_rights_for_women                                             0.456
## Vulnerable_employment                                                 0.746
## Early_marriage                                                        0.554
## Satisfied_demand_for_contraception                                    0.430
## Corruption                                                            0.633
## Acceptance_of_gays_and_lesbians                                       0.395
## Discrimination_and_violence_against_minorities                        0.217
## Equality_of_political_power_by_gender                                 0.260
## Equality_of_political_power_by_socioeconomic_position                 0.229
## Equality_of_political_power_by_social_group                           0.217
## Years_of_tertiary_schooling                                           0.644
## Women_s_average_years_in_school                                       0.763
## Globally_ranked_universities                                          0.091
## Percent_of_tertiary_students_enrolled_in_globally_ranked_universities 0.381
##                                                                          u2 com
## Undernourishment                                                      0.352   1
## Maternal_mortality_rate                                               0.304   1
## Child_mortality_rate                                                  0.224   1
## Child_stunting                                                        0.226   1
## Deaths_from_infectious_diseases                                       0.372   1
## Access_to_at_least_basic_drinking_water                               0.252   1
## Access_to_piped_water                                                 0.275   1
## Access_to_at_least_basic_sanitation_facilities                        0.283   1
## Rural_open_defecation                                                 0.590   1
## Access_to_electricity                                                 0.346   1
## Quality_of_electricity_supply                                         0.241   1
## Household_air_pollution_attributable_deaths                           0.302   1
## Access_to_clean_fuels_and_technology_for_cooking                      0.276   1
## Homicide_rate                                                         0.974   1
## Perceived_criminality                                                 0.650   1
## Political_killings_and_torture                                        0.574   1
## Traffic_deaths                                                        0.537   1
## Adult_literacy_rate                                                   0.397   1
## Primary_school_enrollment                                             0.565   1
## Secondary_school_enrollment                                           0.175   1
## Gender_parity_in_secondary_enrollment                                 0.705   1
## Access_to_quality_education                                           0.431   1
## Mobile_telephone_subscriptions                                        0.559   1
## Internet_users                                                        0.176   1
## Access_to_online_governance                                           0.406   1
## Media_censorship                                                      0.796   1
## Life_expectancy_at_60                                                 0.290   1
## Premature_deaths_from_non_communicable_diseases                       0.628   1
## Access_to_essential_services                                          0.096   1
## Access_to_quality_healthcare                                          0.306   1
## Outdoor_air_pollution_attributable_deaths                             0.965   1
## Greenhouse_gas_emissions                                              0.801   1
## Biome_protection                                                      0.940   1
## Political_rights                                                      0.620   1
## Freedom_of_expression                                                 0.785   1
## Freedom_of_religion                                                   0.895   1
## Access_to_justice                                                     0.551   1
## Property_rights_for_women                                             0.544   1
## Vulnerable_employment                                                 0.254   1
## Early_marriage                                                        0.446   1
## Satisfied_demand_for_contraception                                    0.570   1
## Corruption                                                            0.367   1
## Acceptance_of_gays_and_lesbians                                       0.605   1
## Discrimination_and_violence_against_minorities                        0.783   1
## Equality_of_political_power_by_gender                                 0.740   1
## Equality_of_political_power_by_socioeconomic_position                 0.771   1
## Equality_of_political_power_by_social_group                           0.783   1
## Years_of_tertiary_schooling                                           0.356   1
## Women_s_average_years_in_school                                       0.237   1
## Globally_ranked_universities                                          0.909   1
## Percent_of_tertiary_students_enrolled_in_globally_ranked_universities 0.619   1
## 
##                  MR1
## SS loadings    25.15
## Proportion Var  0.49
## 
## Mean item complexity =  1
## Test of the hypothesis that 1 factor is sufficient.
## 
## The degrees of freedom for the null model are  1275  and the objective function was  75.4 with Chi Square of  11846
## The degrees of freedom for the model are 1224  and the objective function was  37.8 
## 
## The root mean square of the residuals (RMSR) is  0.13 
## The df corrected root mean square of the residuals is  0.13 
## 
## The harmonic number of observations is  176 with the empirical chi square  7737  with prob <  0 
## The total number of observations was  176  with Likelihood Chi Square =  5909  with prob <  0 
## 
## Tucker Lewis Index of factoring reliability =  0.536
## RMSEA index =  0.147  and the 90 % confidence intervals are  0.144 0.152
## BIC =  -420
## Fit based upon off diagonal values = 0.93
## Measures of factor score adequacy             
##                                                    MR1
## Correlation of (regression) scores with factors   1.00
## Multiple R square of scores with factors          0.99
## Minimum correlation of possible factor scores     0.98

(SPI_bhn = fa(SPI_imputed %>% select(Undernourishment:Traffic_deaths)))

## Factor Analysis using method =  minres
## Call: fa(r = SPI_imputed %>% select(Undernourishment:Traffic_deaths))
## Standardized loadings (pattern matrix) based upon correlation matrix
##                                                    MR1    h2   u2 com
## Undernourishment                                  0.83 0.692 0.31   1
## Maternal_mortality_rate                           0.89 0.787 0.21   1
## Child_mortality_rate                              0.92 0.839 0.16   1
## Child_stunting                                    0.87 0.763 0.24   1
## Deaths_from_infectious_diseases                   0.89 0.785 0.21   1
## Access_to_at_least_basic_drinking_water          -0.93 0.866 0.13   1
## Access_to_piped_water                            -0.87 0.761 0.24   1
## Access_to_at_least_basic_sanitation_facilities   -0.93 0.874 0.13   1
## Rural_open_defecation                             0.70 0.494 0.51   1
## Access_to_electricity                            -0.90 0.815 0.18   1
## Quality_of_electricity_supply                    -0.82 0.680 0.32   1
## Household_air_pollution_attributable_deaths       0.90 0.808 0.19   1
## Access_to_clean_fuels_and_technology_for_cooking -0.91 0.826 0.17   1
## Homicide_rate                                     0.14 0.019 0.98   1
## Perceived_criminality                             0.47 0.225 0.77   1
## Political_killings_and_torture                   -0.50 0.253 0.75   1
## Traffic_deaths                                    0.66 0.430 0.57   1
## 
##                  MR1
## SS loadings    10.92
## Proportion Var  0.64
## 
## Mean item complexity =  1
## Test of the hypothesis that 1 factor is sufficient.
## 
## The degrees of freedom for the null model are  136  and the objective function was  20.9 with Chi Square of  3515
## The degrees of freedom for the model are 119  and the objective function was  3.58 
## 
## The root mean square of the residuals (RMSR) is  0.07 
## The df corrected root mean square of the residuals is  0.07 
## 
## The harmonic number of observations is  176 with the empirical chi square  219  with prob <  7.4e-08 
## The total number of observations was  176  with Likelihood Chi Square =  600  with prob <  2.6e-65 
## 
## Tucker Lewis Index of factoring reliability =  0.837
## RMSEA index =  0.151  and the 90 % confidence intervals are  0.14 0.164
## BIC =  -15.1
## Fit based upon off diagonal values = 0.99
## Measures of factor score adequacy             
##                                                    MR1
## Correlation of (regression) scores with factors   0.99
## Multiple R square of scores with factors          0.98
## Minimum correlation of possible factor scores     0.96

(SPI_fw = fa(SPI_imputed %>% select(Adult_literacy_rate:Biome_protection)))

## Factor Analysis using method =  minres
## Call: fa(r = SPI_imputed %>% select(Adult_literacy_rate:Biome_protection))
## Standardized loadings (pattern matrix) based upon correlation matrix
##                                                   MR1    h2    u2 com
## Adult_literacy_rate                              0.76 0.580 0.420   1
## Primary_school_enrollment                        0.65 0.425 0.575   1
## Secondary_school_enrollment                      0.89 0.798 0.202   1
## Gender_parity_in_secondary_enrollment            0.52 0.272 0.728   1
## Access_to_quality_education                      0.76 0.584 0.416   1
## Mobile_telephone_subscriptions                   0.67 0.448 0.552   1
## Internet_users                                   0.91 0.827 0.173   1
## Access_to_online_governance                      0.78 0.604 0.396   1
## Media_censorship                                 0.39 0.151 0.849   1
## Life_expectancy_at_60                            0.86 0.734 0.266   1
## Premature_deaths_from_non_communicable_diseases -0.64 0.410 0.590   1
## Access_to_essential_services                     0.96 0.922 0.078   1
## Access_to_quality_healthcare                     0.84 0.710 0.290   1
## Outdoor_air_pollution_attributable_deaths       -0.20 0.039 0.961   1
## Greenhouse_gas_emissions                        -0.44 0.192 0.808   1
## Biome_protection                                 0.25 0.060 0.940   1
## 
##                 MR1
## SS loadings    7.76
## Proportion Var 0.48
## 
## Mean item complexity =  1
## Test of the hypothesis that 1 factor is sufficient.
## 
## The degrees of freedom for the null model are  120  and the objective function was  14.9 with Chi Square of  2512
## The degrees of freedom for the model are 104  and the objective function was  4.61 
## 
## The root mean square of the residuals (RMSR) is  0.09 
## The df corrected root mean square of the residuals is  0.1 
## 
## The harmonic number of observations is  176 with the empirical chi square  372  with prob <  9.1e-32 
## The total number of observations was  176  with Likelihood Chi Square =  775  with prob <  3.4e-103 
## 
## Tucker Lewis Index of factoring reliability =  0.675
## RMSEA index =  0.191  and the 90 % confidence intervals are  0.18 0.205
## BIC =  238
## Fit based upon off diagonal values = 0.96
## Measures of factor score adequacy             
##                                                    MR1
## Correlation of (regression) scores with factors   0.99
## Multiple R square of scores with factors          0.97
## Minimum correlation of possible factor scores     0.94

(SPI_o = fa(SPI_imputed %>% select(Political_rights:Percent_of_tertiary_students_enrolled_in_globally_ranked_universities)))

## Factor Analysis using method =  minres
## Call: fa(r = SPI_imputed %>% select(Political_rights:Percent_of_tertiary_students_enrolled_in_globally_ranked_universities))
## Standardized loadings (pattern matrix) based upon correlation matrix
##                                                                         MR1
## Political_rights                                                       0.86
## Freedom_of_expression                                                  0.76
## Freedom_of_religion                                                    0.60
## Access_to_justice                                                      0.85
## Property_rights_for_women                                              0.75
## Vulnerable_employment                                                 -0.63
## Early_marriage                                                        -0.57
## Satisfied_demand_for_contraception                                     0.60
## Corruption                                                             0.86
## Acceptance_of_gays_and_lesbians                                        0.69
## Discrimination_and_violence_against_minorities                        -0.58
## Equality_of_political_power_by_gender                                  0.72
## Equality_of_political_power_by_socioeconomic_position                  0.68
## Equality_of_political_power_by_social_group                            0.70
## Years_of_tertiary_schooling                                            0.68
## Women_s_average_years_in_school                                        0.70
## Globally_ranked_universities                                           0.30
## Percent_of_tertiary_students_enrolled_in_globally_ranked_universities  0.58
##                                                                          h2
## Political_rights                                                      0.733
## Freedom_of_expression                                                 0.571
## Freedom_of_religion                                                   0.362
## Access_to_justice                                                     0.717
## Property_rights_for_women                                             0.567
## Vulnerable_employment                                                 0.392
## Early_marriage                                                        0.330
## Satisfied_demand_for_contraception                                    0.360
## Corruption                                                            0.739
## Acceptance_of_gays_and_lesbians                                       0.475
## Discrimination_and_violence_against_minorities                        0.333
## Equality_of_political_power_by_gender                                 0.517
## Equality_of_political_power_by_socioeconomic_position                 0.468
## Equality_of_political_power_by_social_group                           0.495
## Years_of_tertiary_schooling                                           0.465
## Women_s_average_years_in_school                                       0.494
## Globally_ranked_universities                                          0.087
## Percent_of_tertiary_students_enrolled_in_globally_ranked_universities 0.336
##                                                                         u2 com
## Political_rights                                                      0.27   1
## Freedom_of_expression                                                 0.43   1
## Freedom_of_religion                                                   0.64   1
## Access_to_justice                                                     0.28   1
## Property_rights_for_women                                             0.43   1
## Vulnerable_employment                                                 0.61   1
## Early_marriage                                                        0.67   1
## Satisfied_demand_for_contraception                                    0.64   1
## Corruption                                                            0.26   1
## Acceptance_of_gays_and_lesbians                                       0.53   1
## Discrimination_and_violence_against_minorities                        0.67   1
## Equality_of_political_power_by_gender                                 0.48   1
## Equality_of_political_power_by_socioeconomic_position                 0.53   1
## Equality_of_political_power_by_social_group                           0.50   1
## Years_of_tertiary_schooling                                           0.54   1
## Women_s_average_years_in_school                                       0.51   1
## Globally_ranked_universities                                          0.91   1
## Percent_of_tertiary_students_enrolled_in_globally_ranked_universities 0.66   1
## 
##                 MR1
## SS loadings    8.44
## Proportion Var 0.47
## 
## Mean item complexity =  1
## Test of the hypothesis that 1 factor is sufficient.
## 
## The degrees of freedom for the null model are  153  and the objective function was  15.3 with Chi Square of  2570
## The degrees of freedom for the model are 135  and the objective function was  5.99 
## 
## The root mean square of the residuals (RMSR) is  0.13 
## The df corrected root mean square of the residuals is  0.14 
## 
## The harmonic number of observations is  176 with the empirical chi square  925  with prob <  6.8e-118 
## The total number of observations was  176  with Likelihood Chi Square =  1003  with prob <  1.5e-132 
## 
## Tucker Lewis Index of factoring reliability =  0.591
## RMSEA index =  0.191  and the 90 % confidence intervals are  0.181 0.203
## BIC =  305
## Fit based upon off diagonal values = 0.93
## Measures of factor score adequacy             
##                                                    MR1
## Correlation of (regression) scores with factors   0.98
## Multiple R square of scores with factors          0.95
## Minimum correlation of possible factor scores     0.91

#scores
SPI_imputed$SPI_g = SPI_g$scores[, 1] %>% standardize()
SPI_imputed$SPI_bhn = SPI_bhn$scores[, 1] %>% standardize() %>% multiply_by(-1)
SPI_imputed$SPI_fw = SPI_fw$scores[, 1] %>% standardize()
SPI_imputed$SPI_o = SPI_o$scores[, 1] %>% standardize()

#correlations
SPI_imputed %>% select(SPI_g:SPI_o) %>% GG_heatmap(legend_position = c(.4, .7))

#merge to main
#remove existing SPI variables in main
d = d[, setdiff(names(d), names(SPI_imputed[-1]))]
assert_that(!any(duplicated(SPI_imputed$ISO)))

## [1] TRUE

assert_that(!any(duplicated(d$ISO)))

## [1] TRUE

d = full_join(d, SPI_imputed, by = c("ISO" = "ISO"))
assert_that(!any(duplicated(d$ISO)))

## [1] TRUE

Spatial

#join with spatial data
d_geo = full_join(world, d, by = "ISO")
assert_that(!any(duplicated(d_geo$ISO)))

## [1] TRUE

#compute spatial lag variables
d_geo = suppressMessages(add_spatial_lag(d_geo, c("HDI2018", econ_outcomes, names(SPI_imputed[-1]))))
assert_that(!any(duplicated(d_geo$ISO)))

## [1] TRUE

#back to no geodata
d = d_geo %>% st_drop_geometry()

#assert no merge issues
assert_that(!any(duplicated(d$ISO)))

## [1] TRUE

Prep subsets and standardized variants

#standardized data subset
std_sub = d %>% 
  select(
    #basic meta
    Names, ISO, UN_macroregion,
    
    #cognitive
    IQlv, IQr, IQb,
    
    #SPI/well-being
    SPI_g, SPI_g_lag,
    HDI2018, HDI2018_lag,
    !!SPI_indicators, !!(SPI_indicators + "_lag"),
    !!econ_outcomes, !!(econ_outcomes + "_lag"),
    
    #noncognitive
    time_pref_wang, time_pref_falk, patience:trust, time_pref_rieger10, time_pref_rieger6, popsqrt
    ) %>% 
  df_standardize(messages = F, exclude = c("popsqrt"))

#reverse negative indicators
for (v in SPI_indicators) {
  #negative loading? reverse
  if (SPI_g$loadings[, 1][v] < 0) std_sub[[v]] = std_sub[[v]] * -1
}

#subsets
#to keep sample size consistent across models
#the function call removes unused factor levels, which cause errors
std_sub_falk = std_sub %>% filter(!is.na(time_pref_falk)) %>% df_drop_fct_levels()
std_sub_falk_r = std_sub %>% filter(!is.na(time_pref_falk), !is.na(IQr)) %>% df_drop_fct_levels()
std_sub_falk_b = std_sub %>% filter(!is.na(time_pref_falk), !is.na(IQb)) %>% df_drop_fct_levels()

std_sub_wang = std_sub %>% filter(!is.na(time_pref_wang)) %>% df_drop_fct_levels()

std_sub_rieger10 = std_sub %>% filter(!is.na(time_pref_rieger10)) %>% df_drop_fct_levels()
std_sub_rieger6 = std_sub %>% filter(!is.na(time_pref_rieger6)) %>% df_drop_fct_levels()

Descriptive

#Wang data sample sizes
timepref_wang$N %>% describe()

#main variables
d %>% select(
  #IQ variables + HDI
  IQlv, IQr, HDI2018, 
  #time pref
  time_pref_wang, time_pref_falk, time_pref_rieger10, time_pref_rieger6,
  #SPI variables
  SPI_g:SPI_o,
  Undernourishment:Percent_of_tertiary_students_enrolled_in_globally_ranked_universities) %>% 
  describe()

#UN codings
d$UN_continent %>% table2()

d$UN_macroregion %>% table2()

d$UN_region %>% table2()

Restriction of range

How bad is restriction of range? Compare SDs of national IQ and SPI index.

#unrestricted SD
d %>% select(IQ, SPI_g) %>% describe()

#restricted SDs
d %>% filter(!is.na(time_pref_falk)) %>% select(IQ, SPI_g) %>% describe()

d %>% filter(!is.na(time_pref_wang)) %>% select(IQ, SPI_g) %>% describe()

d %>% filter(!is.na(time_pref_rieger10)) %>% select(IQ, SPI_g) %>% describe()

Correlations & Scatterplots

#weighted and unweighted cors
combine_upperlower(
  .upper.tri = d %>% select(IQlv, IQb, IQr, IQwb, SPI_g, HDI2018, time_pref_wang, time_pref_falk, time_pref_rieger10, time_pref_rieger6, risktaking:trust) %>% wtd.cors(),
  .lower.tri = d %>% select(IQlv, IQb, IQr, IQwb, SPI_g, HDI2018, time_pref_wang, time_pref_falk, time_pref_rieger10, time_pref_rieger6, risktaking:trust) %>% wtd.cors(weight = d$popsqrt)
)

##                      IQlv     IQb    IQr   IQwb   SPI_g HDI2018 time_pref_wang
## IQlv                   NA  0.8782  0.980  0.889  0.8004  0.8150         0.6269
## IQb                 0.908      NA  0.878  0.850  0.7262  0.7676         0.6387
## IQr                 0.976  0.9089     NA  0.908  0.8207  0.8327         0.5767
## IQwb                0.878  0.8562  0.902     NA  0.9254  0.9395         0.6314
## SPI_g               0.772  0.7174  0.800  0.918      NA  0.9566         0.6079
## HDI2018             0.824  0.7783  0.843  0.934  0.9622      NA         0.6630
## time_pref_wang      0.612  0.6041  0.526  0.616  0.6185  0.6638             NA
## time_pref_falk      0.612  0.6231  0.626  0.634  0.5633  0.5937         0.5750
## time_pref_rieger10  0.574  0.5890  0.561  0.598  0.5676  0.5839         0.8183
## time_pref_rieger6   0.589  0.6022  0.576  0.596  0.5487  0.5694         0.8122
## risktaking         -0.328 -0.2388 -0.271 -0.200 -0.1401 -0.1150        -0.3184
## posrecip            0.372  0.3006  0.357  0.283  0.1853  0.2137         0.2139
## negrecip            0.200  0.1976  0.212  0.140  0.1022  0.1710        -0.0115
## altruism            0.186  0.0958  0.181  0.108  0.0523  0.0787         0.0623
## trust               0.431  0.3442  0.414  0.229  0.1437  0.2394         0.1199
##                    time_pref_falk time_pref_rieger10 time_pref_rieger6
## IQlv                        0.537            0.55225           0.56577
## IQb                         0.540            0.58763           0.59858
## IQr                         0.550            0.56156           0.57388
## IQwb                        0.601            0.59515           0.59319
## SPI_g                       0.586            0.57299           0.55425
## HDI2018                     0.591            0.57910           0.56695
## time_pref_wang              0.625            0.81879           0.82857
## time_pref_falk                 NA            0.89020           0.89667
## time_pref_rieger10          0.838                 NA           0.95755
## time_pref_rieger6           0.850            0.96405                NA
## risktaking                  0.192            0.00988          -0.00632
## posrecip                    0.127            0.08332           0.07039
## negrecip                    0.199            0.22480           0.23122
## altruism                    0.166            0.06181           0.07814
## trust                       0.247            0.12091           0.12212
##                    risktaking posrecip negrecip altruism   trust
## IQlv                  -0.3896   0.2469   0.2725  -0.0335  0.3818
## IQb                   -0.2978   0.1721   0.3026  -0.1215  0.2959
## IQr                   -0.3575   0.2204   0.2749  -0.0419  0.3716
## IQwb                  -0.2747   0.2289   0.2060  -0.0439  0.2662
## SPI_g                 -0.1915   0.1865   0.1466  -0.0122  0.2511
## HDI2018               -0.1557   0.1830   0.2219  -0.0189  0.3211
## time_pref_wang        -0.2140   0.0547  -0.0341  -0.1673  0.2239
## time_pref_falk         0.2304   0.0155   0.2581  -0.0104  0.1899
## time_pref_rieger10     0.1327   0.0749   0.2232  -0.0464  0.1687
## time_pref_rieger6      0.1087   0.0240   0.2560  -0.0648  0.1581
## risktaking                 NA  -0.2561   0.1928  -0.0152 -0.0615
## posrecip              -0.2149       NA  -0.1536   0.7110  0.3631
## negrecip               0.1398  -0.1693       NA  -0.1320  0.1599
## altruism              -0.0495   0.7804  -0.1103       NA  0.2727
## trust                  0.0278   0.5580   0.1291   0.4703      NA

#heatmap
d %>% 
  select(IQlv, IQb, IQr, IQwb, SPI_g, HDI2018, time_pref_wang, time_pref_falk, time_pref_rieger10,time_pref_rieger6, risktaking:trust) %>% 
  wtd.cors(weight = d$popsqrt) %>% 
  GG_heatmap(reorder_vars = F)

GG_save("figs/heatmap_correlations.png")

#plots
#Falk data
GG_scatter(d, "IQlv", "time_pref_falk", case_names = "Names", color = "UN_macroregion", repel_names = T, text_pos = "tl", weights = "popsqrt") +
  # ggtitle("Time preference (Falk et al) and national IQ") +
  scale_x_continuous("National IQ (L&V 2012)") +
  scale_color_discrete("Region")

## `geom_smooth()` using formula 'y ~ x'

GG_save("figs/IQ_time_pref_falk.png")

## `geom_smooth()` using formula 'y ~ x'

#Wang data
GG_scatter(d, "IQlv", "time_pref_wang", case_names = "Names", color = "UN_macroregion", repel_names = T, text_pos = "tl", weights = "popsqrt") +
  ggtitle("Time preference (Wang et al) and national IQ") +
  scale_x_continuous("National IQ (L&V 2012)") +
  scale_color_discrete("Region")

## `geom_smooth()` using formula 'y ~ x'

GG_save("figs/IQ_time_pref_wang.png")

## `geom_smooth()` using formula 'y ~ x'

#Rieger meta-analysis
#10 datasets
GG_scatter(d, "IQlv", "time_pref_rieger10", case_names = "Names", color = "UN_macroregion", repel_names = T, text_pos = "tl", weights = "popsqrt") +
  ggtitle("Time preference (general) and national IQ") +
  scale_x_continuous("National IQ (L&V 2012)") +
  scale_color_discrete("Region")

## `geom_smooth()` using formula 'y ~ x'

GG_save("figs/IQ_time_pref_rieger10.png")

## `geom_smooth()` using formula 'y ~ x'

#6 datasets
GG_scatter(d, "IQlv", "time_pref_rieger6", case_names = "Names", color = "UN_macroregion", repel_names = T, text_pos = "tl", weights = "popsqrt") +
  ggtitle("Time preference (general) and national IQ") +
  scale_x_continuous("National IQ (L&V 2012)") +
  scale_color_discrete("Region")

## `geom_smooth()` using formula 'y ~ x'

GG_save("figs/IQ_time_pref_rieger6.png")

## `geom_smooth()` using formula 'y ~ x'

Jensen’s method

Do variables with stronger loadings on S also correlate more strongly with X?

#IQ
fa_Jensens_method(SPI_g, d, "IQlv") +
  scale_x_continuous("S loading") +
  scale_y_continuous("Correlation (national IQ x S indicator)") +
  coord_cartesian(xlim = c(0, 1.05), ylim = c(0, 1))

## Using Pearson correlations for the criterion-indicators relationships.

## `geom_smooth()` using formula 'y ~ x'

GG_save("figs/jensen_method_IQlv.png")

## `geom_smooth()` using formula 'y ~ x'

#Time pref Falk
fa_Jensens_method(SPI_g, d, "time_pref_falk") +
  scale_x_continuous("S loading") +
  scale_y_continuous("Correlation (time preference x S indicator)") +
  coord_cartesian(xlim = c(0, 1.05), ylim = c(0, 1))

## Using Pearson correlations for the criterion-indicators relationships.
## `geom_smooth()` using formula 'y ~ x'

GG_save("figs/jensen_method_falk.png")

## `geom_smooth()` using formula 'y ~ x'

#Time pref Rieger10
fa_Jensens_method(SPI_g, d, "time_pref_rieger10") +
  scale_x_continuous("S loading") +
  scale_y_continuous("Correlation (time preference x S indicator)") +
  coord_cartesian(xlim = c(0, 1.05), ylim = c(0, 1))

## Using Pearson correlations for the criterion-indicators relationships.
## `geom_smooth()` using formula 'y ~ x'

GG_save("figs/jensen_method_rieger10.png")

## `geom_smooth()` using formula 'y ~ x'

Models

Main set

#SPI ~ LV + Falk
list(
  ols(SPI_g ~ IQlv, data = std_sub, weights = std_sub$popsqrt),
  ols(SPI_g ~ IQlv, data = std_sub_falk, weights = std_sub_falk$popsqrt),
  ols(SPI_g ~ time_pref_falk, data = std_sub_falk, weights = std_sub_falk$popsqrt),
  ols(SPI_g ~ IQlv + time_pref_falk, data = std_sub_falk, weights = std_sub_falk$popsqrt),
  ols(SPI_g ~ IQlv + time_pref_falk + UN_macroregion, data = std_sub_falk, weights = std_sub_falk$popsqrt),
  #other measures
  ols(SPI_g ~ IQlv + time_pref_falk + risktaking + posrecip + negrecip + altruism + trust, data = std_sub_falk, weights = std_sub_falk$popsqrt),
  ols(SPI_g ~ IQlv + time_pref_falk + risktaking + posrecip + negrecip + altruism + trust + UN_macroregion, data = std_sub_falk, weights = std_sub_falk$popsqrt),
  ols(SPI_g ~ SPI_g_lag + IQlv + time_pref_falk + risktaking + posrecip + negrecip + altruism + trust, data = std_sub_falk, weights = std_sub_falk$popsqrt)
) ->
  model_set_P


#SPI ~ LV + Wang
list(
  ols(SPI_g ~ IQlv, data = std_sub, weights = std_sub$popsqrt),
  ols(SPI_g ~ IQlv, data = std_sub_wang, weights = std_sub_wang$popsqrt),
  ols(SPI_g ~ time_pref_wang, data = std_sub_wang, weights = std_sub_wang$popsqrt),
  ols(SPI_g ~ IQlv + time_pref_wang, data = std_sub_wang, weights = std_sub_wang$popsqrt),
  ols(SPI_g ~ IQlv + time_pref_wang + UN_macroregion, data = std_sub_wang, weights = std_sub_wang$popsqrt),
  ols(SPI_g ~ SPI_g_lag + IQlv + time_pref_wang, data = std_sub_wang, weights = std_sub_wang$popsqrt)
)  ->
  model_set_S1


#SPI ~ LV + Rieger10
list(
  ols(SPI_g ~ IQlv, data = std_sub, weights = std_sub$popsqrt),
  ols(SPI_g ~ IQlv, data = std_sub_rieger10, weights = std_sub_rieger10$popsqrt),
  ols(SPI_g ~ time_pref_rieger10, data = std_sub_rieger10, weights = std_sub_rieger10$popsqrt),
  ols(SPI_g ~ IQlv + time_pref_rieger10, data = std_sub_rieger10, weights = std_sub_rieger10$popsqrt),
  ols(SPI_g ~ IQlv + time_pref_rieger10 + UN_macroregion, data = std_sub_rieger10, weights = std_sub_rieger10$popsqrt),
  ols(SPI_g ~ SPI_g_lag + IQlv + time_pref_rieger10, data = std_sub_rieger10, weights = std_sub_rieger10$popsqrt)
)  ->
  model_set_S2

#HDI ~ LV + Falk
list(
  ols(HDI2018 ~ IQlv, data = std_sub, weights = std_sub$popsqrt),
  ols(HDI2018 ~ IQlv, data = std_sub_falk, weights = std_sub_falk$popsqrt),
  ols(HDI2018 ~ time_pref_falk, data = std_sub_falk, weights = std_sub_falk$popsqrt),
  ols(HDI2018 ~ IQlv + time_pref_falk, data = std_sub_falk, weights = std_sub_falk$popsqrt),
  ols(HDI2018 ~ IQlv + time_pref_falk + UN_macroregion, data = std_sub_falk, weights = std_sub_falk$popsqrt),
  ols(HDI2018 ~ HDI2018_lag + IQlv + time_pref_falk + risktaking + posrecip + negrecip + altruism + trust, data = std_sub_falk, weights = std_sub_falk$popsqrt)
)   ->
  model_set_S3

#SPI ~ Rindermann + Falk
list(
  ols(SPI_g ~ IQr, data = std_sub, weights = std_sub$popsqrt),
  ols(SPI_g ~ IQr, data = std_sub_falk_r, weights = std_sub_falk_r$popsqrt),
  ols(SPI_g ~ time_pref_falk, data = std_sub_falk_r, weights = std_sub_falk_r$popsqrt),
  ols(SPI_g ~ IQr + time_pref_falk, data = std_sub_falk_r, weights = std_sub_falk_r$popsqrt),
  ols(SPI_g ~ IQr + time_pref_falk + UN_macroregion, data = std_sub_falk_r, weights = std_sub_falk_r$popsqrt),
  #other measures
  ols(SPI_g ~ IQr + time_pref_falk + risktaking + posrecip + negrecip + altruism + trust, data = std_sub_falk_r, weights = std_sub_falk_r$popsqrt),
  ols(SPI_g ~ IQr + time_pref_falk + risktaking + posrecip + negrecip + altruism + trust + UN_macroregion, data = std_sub_falk_r, weights = std_sub_falk_r$popsqrt),
  ols(SPI_g ~ SPI_g_lag + IQr + time_pref_falk + risktaking + posrecip + negrecip + altruism + trust, data = std_sub_falk_r, weights = std_sub_falk_r$popsqrt)
) ->
  model_set_S4

#SPI ~ Becker + Falk
list(
  ols(SPI_g ~ IQb, data = std_sub, weights = std_sub$popsqrt),
  ols(SPI_g ~ IQb, data = std_sub_falk_b, weights = std_sub_falk_b$popsqrt),
  ols(SPI_g ~ time_pref_falk, data = std_sub_falk_b, weights = std_sub_falk_b$popsqrt),
  ols(SPI_g ~ IQb + time_pref_falk, data = std_sub_falk_b, weights = std_sub_falk_b$popsqrt),
  ols(SPI_g ~ IQb + time_pref_falk + UN_macroregion, data = std_sub_falk_b, weights = std_sub_falk_b$popsqrt),
  #other measures
  ols(SPI_g ~ IQb + time_pref_falk + risktaking + posrecip + negrecip + altruism + trust, data = std_sub_falk_b, weights = std_sub_falk_b$popsqrt),
  ols(SPI_g ~ IQb + time_pref_falk + risktaking + posrecip + negrecip + altruism + trust + UN_macroregion, data = std_sub_falk_b, weights = std_sub_falk_b$popsqrt),
  ols(SPI_g ~ SPI_g_lag + IQb + time_pref_falk + risktaking + posrecip + negrecip + altruism + trust, data = std_sub_falk_b, weights = std_sub_falk_b$popsqrt)
) ->
  model_set_S5

#SPI ~ LV + Falk / NO weights
list(
  ols(SPI_g ~ IQlv, data = std_sub),
  ols(SPI_g ~ IQlv, data = std_sub_falk),
  ols(SPI_g ~ time_pref_falk, data = std_sub_falk),
  ols(SPI_g ~ IQlv + time_pref_falk, data = std_sub_falk),
  ols(SPI_g ~ IQlv + time_pref_falk + UN_macroregion, data = std_sub_falk),
  #other measures
  ols(SPI_g ~ IQlv + time_pref_falk + risktaking + posrecip + negrecip + altruism + trust, data = std_sub_falk),
  ols(SPI_g ~ IQlv + time_pref_falk + risktaking + posrecip + negrecip + altruism + trust + UN_macroregion, data = std_sub_falk),
  ols(SPI_g ~ SPI_g_lag + IQlv + time_pref_falk + risktaking + posrecip + negrecip + altruism + trust, data = std_sub_falk)
) ->
  model_set_S6

#SPI ~ LV + Falk, robust
#this function is picky with factors
#not included in the main results because added later and requires too much hacking to get to work
#nonetheless shows about the same results
list(
  MASS::rlm(SPI_g ~ IQlv, data = std_sub, weights = std_sub$popsqrt, maxit = 999),
  MASS::rlm(SPI_g ~ IQlv, data = std_sub_falk, weights = std_sub_falk$popsqrt, maxit = 999),
  MASS::rlm(SPI_g ~ time_pref_falk, data = std_sub_falk, weights = std_sub_falk$popsqrt, maxit = 999),
  MASS::rlm(SPI_g ~ IQlv + time_pref_falk, data = std_sub_falk, weights = std_sub_falk$popsqrt, maxit = 999),
  MASS::rlm(SPI_g ~ IQlv + time_pref_falk + UN_macroregion, data = std_sub_falk %>% mutate(UN_macroregion = as.character(UN_macroregion)), weights = std_sub_falk$popsqrt, maxit = 999),
  #other measures
  MASS::rlm(SPI_g ~ IQlv + time_pref_falk + risktaking + posrecip + negrecip + altruism + trust, data = std_sub_falk, weights = std_sub_falk$popsqrt, maxit = 999),
  MASS::rlm(SPI_g ~ IQlv + time_pref_falk + risktaking + posrecip + negrecip + altruism + trust + UN_macroregion, data = std_sub_falk %>% mutate(UN_macroregion = as.character(UN_macroregion)), weights = std_sub_falk$popsqrt, maxit = 999),
  MASS::rlm(SPI_g ~ SPI_g_lag + IQlv + time_pref_falk + risktaking + posrecip + negrecip + altruism + trust, data = std_sub_falk, weights = std_sub_falk$popsqrt, maxit = 999)
) ->
  model_set_S7

#quantitative summary across models
extract_func = function(x) {
  # browser()
  #last model
  xlast = x %>% last() %>% .[[1]]
  
  #model summary
  tibble(
    lag = coef(xlast)[2],
    IQ = coef(xlast)[3],
    time_pref = coef(xlast)[4],
    risktaking = coef(xlast)["risktaking"],
    posrecip = coef(xlast)["posrecip"],
    negrecip = coef(xlast)["negrecip"],
    altruism = coef(xlast)["altruism"],
    trust = coef(xlast)["trust"],
    N = xlast$stats[1],
    R2adj = xlast %>% summary.lm() %>% .$adj.r.squared
  )
}

#list of list of models
ll_models = list(
  #primary set
  model_set_P,
  #supplemental
  model_set_S1,
  model_set_S2,
  model_set_S3,
  model_set_S4,
  model_set_S5,
  model_set_S6
  )

#quantitative summary
map_df(ll_models, extract_func) ->
  param_sum

param_sum

param_sum %>% colMeans(na.rm = T)

##        lag         IQ  time_pref risktaking   posrecip   negrecip   altruism 
##    0.39397    0.41161    0.13547    0.08797   -0.00898   -0.04918    0.00924 
##      trust          N      R2adj 
##   -0.16289   77.42857    0.70340

param_sum %>% describe()

#model accuracy of models
map2_df(ll_models, seq_along(ll_models), function(x, i) {
  #subset to to wanted models
  #note that the lists do not have the same number of models, so we need to use negative indexing
  #models 1-2, are pure national IQ
  #model 3 is pure time pref
  #model 4 is combined national IQ + time pref
  #next to last model is full model with regional dummies
  #last model is full model with spatial lag
  
  xx = x[c(2, 3, 4, length(x)-1, length(x))]
  
  #model adj. R2s
  tibble(
    model_set = i,
    model = c("1. National IQ", "2. Time pref.", "3. National IQ + time pref.", "4. Full model with dummies", "5. Full model with spatial lag"),
    R2adj = map_dbl(xx, ~summary.lm(.) %>% .$adj.r.squared)
  )
}) %>% 
  ggplot(aes(model, R2adj, color = factor(model_set), group = factor(model_set))) +
  geom_line() +
  scale_x_discrete(guide = guide_axis(angle = -7)) +
  scale_y_continuous("Adjusted R²", labels = scales::percent) +
  scale_color_discrete("Model set", labels = c("SPI ~ LV + Falk", "SPI ~ LV + Wang", "SPI ~ LV + Rieger10", "HDI ~ LV + Falk", "SPI ~ R + Falk", "SPI ~ B + Falk", "SPI ~ LV + Falk\nno weights"))

GG_save("figs/model_r2s.png")

#test addition of time pref
map_dbl(ll_models, function(x) {
  #test model 2 vs. 4
  # browser()
  lrtest(x[[2]], x[[4]])$stats[3]
})

## [1] 0.06064 0.01550 0.00619 0.03750 0.14099 0.03452 0.00768

#print model summaries
model_set_P %>% summarize_models(asterisks_only = F)

model_set_S1 %>% summarize_models(asterisks_only = F)

model_set_S2 %>% summarize_models(asterisks_only = F)

model_set_S3 %>% summarize_models(asterisks_only = F)

model_set_S4 %>% summarize_models(asterisks_only = F)

model_set_S5 %>% summarize_models(asterisks_only = F)

model_set_S6 %>% summarize_models(asterisks_only = F)

SPI indicators

Fit specific models for every SPI indicator to assess robustness of result.

#fit function
fit_models = function(y) {
  #prep local subset of data because of bug in rms, we need to drop unused factor levels
  # browser()
  y_sym = as.symbol(y)
  local_subset = std_sub %>% filter(!is.na(!!y_sym)) %>% df_drop_fct_levels()
  local_subset_falk = std_sub %>% filter(!is.na(!!y_sym), !is.na(time_pref_falk)) %>% df_drop_fct_levels()
  
  #fit
  #but with changed formula
  list(
    ols(str_glue("{y} ~ IQlv") %>% as.formula(), data = local_subset, weights = local_subset$popsqrt),
    ols(str_glue("{y} ~ IQlv") %>% as.formula(), data = local_subset_falk, weights = local_subset_falk$popsqrt),
    ols(str_glue("{y} ~ time_pref_falk") %>% as.formula(), data = local_subset_falk, weights = local_subset_falk$popsqrt),
    ols(str_glue("{y} ~ IQlv + time_pref_falk") %>% as.formula(), data = local_subset_falk, weights = local_subset_falk$popsqrt),
    ols(str_glue("{y} ~ IQlv + time_pref_falk + UN_macroregion") %>% as.formula(), data = local_subset_falk, weights = local_subset_falk$popsqrt),
    #other measures
    ols(str_glue("{y} ~ IQlv + time_pref_falk + risktaking + posrecip + negrecip + altruism + trust") %>% as.formula(), data = local_subset_falk, weights = local_subset_falk$popsqrt),
    ols(str_glue("{y} ~ IQlv + time_pref_falk + risktaking + posrecip + negrecip + altruism + trust + UN_macroregion") %>% as.formula(), data = local_subset_falk, weights = local_subset_falk$popsqrt),
    ols(str_glue("{y} ~ {y}_lag + IQlv + time_pref_falk + risktaking + posrecip + negrecip + altruism + trust") %>% as.formula(), data = local_subset_falk, weights = local_subset_falk$popsqrt)
  )
}

#loop and fit for each SPI indicator
SPI_indicator_fits = map(SPI_indicators, fit_models)

#extract summary to long form
SPI_indicator_fits_long = map_df(SPI_indicator_fits, function(x) {
  # browser()
  #summarize
  x %>% 
    summarize_models(beta_digits = 5) %>% 
    #just numbers
    {
      #convert to just betas
      cbind(.[1], map_df(.[-1], function(xx) {
        str_match(xx, "[\\d\\.]+") %>% as.numeric()
      }))
    } %>% 
    mutate(
      outcome = x[[1]] %>% terms() %>% .[[2]] %>% as.character()
    )
})

#print for inspection
SPI_indicator_fits_long

#fix up a bit
SPI_indicator_fits_long2 = SPI_indicator_fits_long %>% 
  {
    colnames(.)[1:9] = c("predictor", "model_" + 1:8)
    .
  } %>% 
  filter(predictor %in% c("IQlv", "time_pref_falk"))

#model 4
#plot simple comparison
SPI_indicator_fits_long2 %>% 
  GG_denhist("model_4", "predictor", vline = median) +
  scale_x_continuous("Standardized beta")

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

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

GG_save("figs/SPI_indicators_beta_dist.png")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

#numbers
describeBy(SPI_indicator_fits_long2$model_4, SPI_indicator_fits_long2$predictor)

## 
##  Descriptive statistics by group 
## group: IQlv
##    vars  n mean   sd median trimmed  mad  min  max range  skew kurtosis   se
## X1    1 51 0.37 0.21    0.4    0.37 0.26 0.01 0.71   0.7 -0.11    -1.23 0.03
## ------------------------------------------------------------ 
## group: time_pref_falk
##    vars  n mean   sd median trimmed  mad min  max range skew kurtosis   se
## X1    1 51 0.19 0.19   0.14    0.16 0.12   0 1.14  1.14 2.65     9.86 0.03

#relationship
SPI_indicator_fits_long2 %>% 
  select(predictor, model_4, outcome) %>% 
  spread(key = predictor, value = model_4) %>% 
  mutate(
    ratio = IQlv / time_pref_falk
  ) %>% 
  {
    describe(.$ratio) %>% print()
    .[-26, ]
  } %>% 
  GG_scatter("IQlv", "time_pref_falk", case_names = "outcome")

##    vars  n mean   sd median trimmed  mad    min max range skew kurtosis   se
## X1    1 51 14.4 53.1   2.99    3.96 3.46 0.0448 373   373 6.08       38 7.44

## `geom_smooth()` using formula 'y ~ x'

#within model comparison
SPI_indicator_fits_long2 %>% 
  #fix factor level
  mutate(
    outcome = factor(outcome %>% str_clean(), levels = SPI_indicators %>% str_clean()) %>% fct_reorder(model_4)
  ) %>% 
  ggplot(aes(model_4, outcome)) +
  geom_point(aes(color = predictor)) +
  # scale_color_discrete(guide = F) +
  scale_x_continuous("Standardized beta", breaks = seq(0, 2, .2))

GG_save("figs/SPI_indicators_beta_by_side.png")

#model 7
#plot simple comparison
SPI_indicator_fits_long2 %>% 
  GG_denhist("model_7", "predictor", vline = median) +
  scale_x_continuous("Standardized beta")

## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

GG_save("figs/SPI_indicators_beta_dist_full.png")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

#numbers
describeBy(SPI_indicator_fits_long2$model_7, SPI_indicator_fits_long2$predictor)

## 
##  Descriptive statistics by group 
## group: IQlv
##    vars  n mean   sd median trimmed  mad  min  max range skew kurtosis   se
## X1    1 51  0.4 0.21   0.39     0.4 0.23 0.01 0.85  0.84 0.16    -0.78 0.03
## ------------------------------------------------------------ 
## group: time_pref_falk
##    vars  n mean   sd median trimmed mad min  max range skew kurtosis   se
## X1    1 51 0.14 0.12   0.11    0.13 0.1   0 0.66  0.66 1.83     4.91 0.02

#relationship
SPI_indicator_fits_long2 %>% 
  select(predictor, model_7, outcome) %>% 
  spread(key = predictor, value = model_7) %>% 
  mutate(
    ratio = IQlv / time_pref_falk
  ) %>% 
  {
    describe(.$ratio) %>% print()
    .[-26, ]
  } %>% 
  GG_scatter("IQlv", "time_pref_falk", case_names = "outcome")

##    vars  n mean   sd median trimmed  mad    min max range skew kurtosis   se
## X1    1 51 11.2 38.8   3.47    4.21 2.43 0.0717 275   275 6.22     39.2 5.43

## `geom_smooth()` using formula 'y ~ x'

#within model comparison
SPI_indicator_fits_long2 %>% 
  #fix factor level
  mutate(
    outcome = factor(outcome %>% str_clean(), levels = SPI_indicators %>% str_clean()) %>% fct_reorder(model_7)
  ) %>% 
  ggplot(aes(model_7, outcome)) +
  geom_point(aes(color = predictor)) +
  scale_color_discrete(labels = c("IQlv", "Time pref.")) +
  scale_x_continuous("Standardized beta", breaks = seq(0, 2, .2))

GG_save("figs/SPI_indicators_beta_by_side_full.png")

Economic outcomes

Since SPI does not include these.

#correlations with other variables
d %>% select(IQlv, SPI_g, !!econ_outcomes) %>% wtd.cors(weight = d$popsqrt)

##                   IQlv  SPI_g median_income median_income_hh GNI_2019
## IQlv             1.000 0.7717       0.78641           0.7344   0.8113
## SPI_g            0.772 1.0000       0.89925           0.8709   0.9337
## median_income    0.786 0.8992       1.00000           0.9689   0.9185
## median_income_hh 0.734 0.8709       0.96886           1.0000   0.9043
## GNI_2019         0.811 0.9337       0.91855           0.9043   1.0000
## GNI_growth       0.380 0.0258      -0.04735          -0.0909   0.0745
## GDP_2019         0.804 0.9327       0.91695           0.9016   0.9991
## GDP_growth       0.384 0.0656      -0.00336          -0.0450   0.0967
##                  GNI_growth GDP_2019 GDP_growth
## IQlv                 0.3803   0.8042    0.38404
## SPI_g                0.0258   0.9327    0.06557
## median_income       -0.0473   0.9169   -0.00336
## median_income_hh    -0.0909   0.9016   -0.04500
## GNI_2019             0.0745   0.9991    0.09675
## GNI_growth           1.0000   0.0766    0.99638
## GDP_2019             0.0766   1.0000    0.10009
## GDP_growth           0.9964   0.1001    1.00000

d %>% select(IQlv, SPI_g, !!econ_outcomes) %>% wtd.cors()

##                   IQlv SPI_g median_income median_income_hh GNI_2019 GNI_growth
## IQlv             1.000 0.800        0.8410           0.7908    0.780     0.2613
## SPI_g            0.800 1.000        0.9009           0.8639    0.916     0.1533
## median_income    0.841 0.901        1.0000           0.9659    0.923     0.0928
## median_income_hh 0.791 0.864        0.9659           1.0000    0.915     0.0798
## GNI_2019         0.780 0.916        0.9229           0.9146    1.000     0.1769
## GNI_growth       0.261 0.153        0.0928           0.0798    0.177     1.0000
## GDP_2019         0.767 0.913        0.9206           0.9109    0.996     0.1819
## GDP_growth       0.244 0.161        0.1462           0.1361    0.177     0.9820
##                  GDP_2019 GDP_growth
## IQlv                0.767      0.244
## SPI_g               0.913      0.161
## median_income       0.921      0.146
## median_income_hh    0.911      0.136
## GNI_2019            0.996      0.177
## GNI_growth          0.182      0.982
## GDP_2019            1.000      0.199
## GDP_growth          0.199      1.000

d %>% select(IQlv, SPI_g, !!econ_outcomes) %>% pairwiseCount()

##                  IQlv SPI_g median_income median_income_hh GNI_2019 GNI_growth
## IQlv              196   174           131              131      170        141
## SPI_g             174   176           130              130      163        133
## median_income     131   130           131              131      124        101
## median_income_hh  131   130           131              131      124        101
## GNI_2019          170   163           124              124      174        142
## GNI_growth        141   133           101              101      142        142
## GDP_2019          172   165           126              126      174        142
## GDP_growth        152   144           109              109      154        142
##                  GDP_2019 GDP_growth
## IQlv                  172        152
## SPI_g                 165        144
## median_income         126        109
## median_income_hh      126        109
## GNI_2019              174        154
## GNI_growth            142        142
## GDP_2019              177        155
## GDP_growth            155        155

#pairwise plots
d %>% select(IQlv, SPI_g, !!econ_outcomes) %>% GGally::ggpairs()

## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2

## Warning: Removed 65 rows containing non-finite values (stat_density).

## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 87 rows containing missing values

## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 130 rows containing missing values

## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 130 rows containing missing values

## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 91 rows containing missing values

## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 120 rows containing missing values

## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 89 rows containing missing values

## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 109 rows containing missing values

## Warning: Removed 87 rows containing missing values (geom_point).

## Warning: Removed 85 rows containing non-finite values (stat_density).

## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 131 rows containing missing values

## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 131 rows containing missing values

## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 98 rows containing missing values

## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 128 rows containing missing values

## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 96 rows containing missing values

## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 117 rows containing missing values

## Warning: Removed 130 rows containing missing values (geom_point).

## Warning: Removed 131 rows containing missing values (geom_point).

## Warning: Removed 130 rows containing non-finite values (stat_density).

## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 130 rows containing missing values

## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 137 rows containing missing values

## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 160 rows containing missing values

## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 135 rows containing missing values

## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 152 rows containing missing values

## Warning: Removed 130 rows containing missing values (geom_point).

## Warning: Removed 131 rows containing missing values (geom_point).

## Warning: Removed 130 rows containing missing values (geom_point).

## Warning: Removed 130 rows containing non-finite values (stat_density).

## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 137 rows containing missing values

## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 160 rows containing missing values

## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 135 rows containing missing values

## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 152 rows containing missing values

## Warning: Removed 91 rows containing missing values (geom_point).

## Warning: Removed 98 rows containing missing values (geom_point).

## Warning: Removed 137 rows containing missing values (geom_point).

## Warning: Removed 137 rows containing missing values (geom_point).

## Warning: Removed 87 rows containing non-finite values (stat_density).

## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 119 rows containing missing values

## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 87 rows containing missing values

## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 107 rows containing missing values

## Warning: Removed 120 rows containing missing values (geom_point).

## Warning: Removed 128 rows containing missing values (geom_point).

## Warning: Removed 160 rows containing missing values (geom_point).

## Warning: Removed 160 rows containing missing values (geom_point).

## Warning: Removed 119 rows containing missing values (geom_point).

## Warning: Removed 119 rows containing non-finite values (stat_density).

## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 119 rows containing missing values

## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 119 rows containing missing values

## Warning: Removed 89 rows containing missing values (geom_point).

## Warning: Removed 96 rows containing missing values (geom_point).

## Warning: Removed 135 rows containing missing values (geom_point).

## Warning: Removed 135 rows containing missing values (geom_point).

## Warning: Removed 87 rows containing missing values (geom_point).

## Warning: Removed 119 rows containing missing values (geom_point).

## Warning: Removed 84 rows containing non-finite values (stat_density).

## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 106 rows containing missing values

## Warning: Removed 109 rows containing missing values (geom_point).

## Warning: Removed 117 rows containing missing values (geom_point).

## Warning: Removed 152 rows containing missing values (geom_point).

## Warning: Removed 152 rows containing missing values (geom_point).

## Warning: Removed 107 rows containing missing values (geom_point).

## Warning: Removed 119 rows containing missing values (geom_point).

## Warning: Removed 106 rows containing missing values (geom_point).

## Warning: Removed 106 rows containing non-finite values (stat_density).

GG_save("figs/economic_pairwise.png")

## Warning: Removed 65 rows containing non-finite values (stat_density).

## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 87 rows containing missing values

## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 130 rows containing missing values

## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 130 rows containing missing values

## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 91 rows containing missing values

## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 120 rows containing missing values

## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 89 rows containing missing values

## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 109 rows containing missing values

## Warning: Removed 87 rows containing missing values (geom_point).

## Warning: Removed 85 rows containing non-finite values (stat_density).

## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 131 rows containing missing values

## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 131 rows containing missing values

## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 98 rows containing missing values

## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 128 rows containing missing values

## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 96 rows containing missing values

## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 117 rows containing missing values

## Warning: Removed 130 rows containing missing values (geom_point).

## Warning: Removed 131 rows containing missing values (geom_point).

## Warning: Removed 130 rows containing non-finite values (stat_density).

## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 130 rows containing missing values

## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 137 rows containing missing values

## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 160 rows containing missing values

## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 135 rows containing missing values

## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 152 rows containing missing values

## Warning: Removed 130 rows containing missing values (geom_point).

## Warning: Removed 131 rows containing missing values (geom_point).

## Warning: Removed 130 rows containing missing values (geom_point).

## Warning: Removed 130 rows containing non-finite values (stat_density).

## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 137 rows containing missing values

## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 160 rows containing missing values

## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 135 rows containing missing values

## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 152 rows containing missing values

## Warning: Removed 91 rows containing missing values (geom_point).

## Warning: Removed 98 rows containing missing values (geom_point).

## Warning: Removed 137 rows containing missing values (geom_point).

## Warning: Removed 137 rows containing missing values (geom_point).

## Warning: Removed 87 rows containing non-finite values (stat_density).

## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 119 rows containing missing values

## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 87 rows containing missing values

## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 107 rows containing missing values

## Warning: Removed 120 rows containing missing values (geom_point).

## Warning: Removed 128 rows containing missing values (geom_point).

## Warning: Removed 160 rows containing missing values (geom_point).

## Warning: Removed 160 rows containing missing values (geom_point).

## Warning: Removed 119 rows containing missing values (geom_point).

## Warning: Removed 119 rows containing non-finite values (stat_density).

## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 119 rows containing missing values

## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 119 rows containing missing values

## Warning: Removed 89 rows containing missing values (geom_point).

## Warning: Removed 96 rows containing missing values (geom_point).

## Warning: Removed 135 rows containing missing values (geom_point).

## Warning: Removed 135 rows containing missing values (geom_point).

## Warning: Removed 87 rows containing missing values (geom_point).

## Warning: Removed 119 rows containing missing values (geom_point).

## Warning: Removed 84 rows containing non-finite values (stat_density).

## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 106 rows containing missing values

## Warning: Removed 109 rows containing missing values (geom_point).

## Warning: Removed 117 rows containing missing values (geom_point).

## Warning: Removed 152 rows containing missing values (geom_point).

## Warning: Removed 152 rows containing missing values (geom_point).

## Warning: Removed 107 rows containing missing values (geom_point).

## Warning: Removed 119 rows containing missing values (geom_point).

## Warning: Removed 106 rows containing missing values (geom_point).

## Warning: Removed 106 rows containing non-finite values (stat_density).

#spearman
d %>% select(IQlv, SPI_g, !!econ_outcomes) %>% cor(use = "pairwise", method = "spearman")

##                   IQlv SPI_g median_income median_income_hh GNI_2019 GNI_growth
## IQlv             1.000 0.814         0.848            0.800    0.788      0.302
## SPI_g            0.814 1.000         0.929            0.887    0.921      0.215
## median_income    0.848 0.929         1.000            0.979    0.952      0.102
## median_income_hh 0.800 0.887         0.979            1.000    0.933      0.062
## GNI_2019         0.788 0.921         0.952            0.933    1.000      0.197
## GNI_growth       0.302 0.215         0.102            0.062    0.197      1.000
## GDP_2019         0.775 0.918         0.950            0.930    0.997      0.196
## GDP_growth       0.287 0.229         0.164            0.131    0.213      0.963
##                  GDP_2019 GDP_growth
## IQlv                0.775      0.287
## SPI_g               0.918      0.229
## median_income       0.950      0.164
## median_income_hh    0.930      0.131
## GNI_2019            0.997      0.213
## GNI_growth          0.196      0.963
## GDP_2019            1.000      0.226
## GDP_growth          0.226      1.000

#refit main models with economic outcomes
econ_models_fit = map(econ_outcomes, fit_models)

#extract summary to long form
econ_fits_long = map_df(econ_models_fit, function(x) {
  # browser()
  #summarize
  x %>% 
    summarize_models(beta_digits = 5) %>% 
    #just numbers
    {
      #convert to just betas
      cbind(.[1], map_df(.[-1], function(xx) {
        str_match(xx, "[\\d\\.]+") %>% as.numeric()
      }))
    } %>% 
    mutate(
      outcome = x[[1]] %>% terms() %>% .[[2]] %>% as.character()
    )
})

#print for inspection
econ_fits_long

#fix up a bit
econ_fits_long2 = econ_fits_long %>% 
  {
    colnames(.)[1:9] = c("predictor", "model_" + 1:8)
    .
  } %>% 
  filter(predictor %in% c("IQlv", "time_pref_falk"))

#numbers
describeBy(econ_fits_long2$model_4, econ_fits_long2$predictor)

## 
##  Descriptive statistics by group 
## group: IQlv
##    vars n mean   sd median trimmed  mad  min  max range  skew kurtosis   se
## X1    1 6 0.54 0.04   0.54    0.54 0.03 0.46 0.58  0.11 -0.81    -0.86 0.02
## ------------------------------------------------------------ 
## group: time_pref_falk
##    vars n mean   sd median trimmed  mad  min  max range skew kurtosis   se
## X1    1 6 0.21 0.08   0.19    0.21 0.08 0.13 0.34  0.21 0.53    -1.53 0.03

describeBy(econ_fits_long2$model_7, econ_fits_long2$predictor)

## 
##  Descriptive statistics by group 
## group: IQlv
##    vars n mean   sd median trimmed  mad  min  max range  skew kurtosis   se
## X1    1 6 0.47 0.18   0.52    0.47 0.19 0.25 0.64   0.4 -0.18    -2.11 0.07
## ------------------------------------------------------------ 
## group: time_pref_falk
##    vars n mean   sd median trimmed  mad  min  max range skew kurtosis   se
## X1    1 6  0.2 0.17   0.15     0.2 0.16 0.01 0.44  0.42 0.34    -1.89 0.07

#relationship
econ_fits_long2 %>% 
  select(predictor, model_7, outcome) %>% 
  spread(key = predictor, value = model_7) %>% 
  mutate(
    ratio = IQlv / time_pref_falk
  ) %>% 
  {
    describe(.$ratio) %>% print()
    .[-26, ]
  } %>% 
  GG_scatter("IQlv", "time_pref_falk", case_names = "outcome")

##    vars n mean   sd median trimmed   mad  min  max range skew kurtosis   se
## X1    1 6 7.74 11.7   1.84    7.74 0.485 1.47 30.9  29.5  1.2   -0.389 4.78

## `geom_smooth()` using formula 'y ~ x'

#within model comparison
econ_fits_long2 %>% 
  #fix factor level
  mutate(
    outcome = factor(outcome %>% str_clean(), levels = econ_outcomes %>% str_clean())
  ) %>% 
  ggplot(aes(model_7, outcome)) +
  geom_point(aes(color = predictor), size = 3) +
  scale_color_discrete(labels = c("IQlv", "Time pref.")) +
  scale_x_continuous("Standardized beta", breaks = seq(0, 2, .1))

GG_save("figs/econ_beta_by_side_full.png")

#direct comparison, model 4
econ_fits_long2 %>% 
  #fix factor level
  mutate(
    outcome = factor(outcome %>% str_clean(), levels = econ_outcomes %>% str_clean()) %>% fct_reorder(model_4)
  ) %>% 
  ggplot(aes(model_4, outcome)) +
  geom_point(aes(color = predictor)) +
  scale_color_discrete(labels = c("IQlv", "Time pref.")) +
  scale_x_continuous("Standardized beta", breaks = seq(0, 2, .1))

GG_save("figs/econ_beta_by_side_direct.png")

Maps

#plot time preference
#max values, so we can set the guide
d_geo$time_pref_rieger10_z %>% describe()

d_geo$time_pref_wang_z %>% describe()

d_geo$time_pref_falk_z %>% describe()

#so -3 to 3 is reasonable

#Falk
ggplot(d_geo) + 
  geom_sf(aes(fill = time_pref_falk_z), lwd = .1) +
  scale_fill_gradient("Time preference\n(Falk et al)", high = "green", low = "red", limits = c(-3, 3)) +
  theme_classic()

GG_save("figs/map_time_pref_falk.png")

#Wang
ggplot(d_geo) + 
  geom_sf(aes(fill = time_pref_wang_z), lwd = .1) +
  scale_fill_gradient("Time preference\n(Wang et al)", high = "green", low = "red", limits = c(-3, 3)) +
  theme_classic()

GG_save("figs/map_time_pref_wang.png")

#Rieger
ggplot(d_geo) + 
  geom_sf(aes(fill = time_pref_rieger10_z), lwd = .1) +
  scale_fill_gradient("Time preference", high = "green", low = "red", limits = c(-3, 3)) +
  theme_classic()

GG_save("figs/map_time_pref_rieger10.png")

#regional coding
map_continents = tm_shape(d_geo) +
  tm_fill("UN_continent", title = "Continents") +
  tm_borders()
map_continents

## Warning: The shape d_geo contains empty units.

tmap_save(map_continents, "figs/continents map.png")

## Warning: The shape d_geo contains empty units.

## Warning: The shape d_geo contains empty units.

## Map saved to /media/luks8tb2/Projects/Time preferences national IQs/figs/continents map.png

## Resolution: 3054 by 1444 pixels

## Size: 10.2 by 4.81 inches (300 dpi)

map_macroregions = tm_shape(d_geo) +
  tm_fill("UN_macroregion", title = "Macroregions") +
  tm_borders()
map_macroregions

## Warning: The shape d_geo contains empty units.

tmap_save(map_macroregions, "figs/macroregions map.png")

## Warning: The shape d_geo contains empty units.

## Warning: The shape d_geo contains empty units.

## Map saved to /media/luks8tb2/Projects/Time preferences national IQs/figs/macroregions map.png

## Resolution: 3054 by 1444 pixels

## Size: 10.2 by 4.81 inches (300 dpi)

map_regions = tm_shape(d_geo) +
  tm_fill("UN_region", title = "Regions") +
  tm_borders()
map_regions

## Warning: The shape d_geo contains empty units.

tmap_save(map_regions, "figs/regions map.png")

## Warning: The shape d_geo contains empty units.

## Warning: The shape d_geo contains empty units.

## Map saved to /media/luks8tb2/Projects/Time preferences national IQs/figs/regions map.png

## Resolution: 3054 by 1444 pixels

## Size: 10.2 by 4.81 inches (300 dpi)

Data print

Some of the main variables

#print data that has data for at least one time pref variable
std_sub %>% 
  select(Names, SPI_g, IQlv, time_pref_falk, time_pref_wang, time_pref_rieger10) %>% 
  filter(!is.na(time_pref_falk) | !is.na(time_pref_wang) | !is.na(time_pref_rieger10))

Meta

write_sessioninfo()

## R version 4.0.3 (2020-10-10)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Linux Mint 19.3
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1
## 
## locale:
##  [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=de_DE.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=de_DE.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=de_DE.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] tmap_3.2              sf_0.9-6              haven_2.3.1          
##  [4] rms_6.0-1             SparseM_1.78          readxl_1.3.1         
##  [7] kirkegaard_2020-11-08 metafor_2.4-0         Matrix_1.2-18        
## [10] psych_2.0.9           magrittr_1.5          assertthat_0.2.1     
## [13] weights_1.0.1         mice_3.11.0           gdata_2.18.0         
## [16] Hmisc_4.4-1           Formula_1.2-4         survival_3.2-7       
## [19] lattice_0.20-41       forcats_0.5.0         stringr_1.4.0        
## [22] dplyr_1.0.2           purrr_0.3.4           readr_1.4.0          
## [25] tidyr_1.1.2           tibble_3.0.4          ggplot2_3.3.2        
## [28] tidyverse_1.3.0       pacman_0.5.1         
## 
## loaded via a namespace (and not attached):
##   [1] backports_1.1.10        lwgeom_0.2-5            plyr_1.8.6             
##   [4] sp_1.4-4                splines_4.0.3           crosstalk_1.1.0.1      
##   [7] leaflet_2.0.3           TH.data_1.0-10          digest_0.6.26          
##  [10] htmltools_0.5.0         fansi_0.4.1             checkmate_2.0.0        
##  [13] cluster_2.1.0           openxlsx_4.2.2          modelr_0.1.8           
##  [16] matrixStats_0.57.0      sandwich_3.0-0          jpeg_0.1-8.1           
##  [19] colorspace_1.4-1        ggrepel_0.8.2           blob_1.2.1             
##  [22] rvest_0.3.6             xfun_0.18               leafem_0.1.3           
##  [25] crayon_1.3.4            jsonlite_1.7.1          zoo_1.8-8              
##  [28] glue_1.4.2              stars_0.4-3             gtable_0.3.0           
##  [31] MatrixModels_0.4-1      car_3.0-10              DEoptimR_1.0-8         
##  [34] abind_1.4-5             VIM_6.0.0               scales_1.1.1           
##  [37] mvtnorm_1.1-1           GGally_2.0.0            DBI_1.1.0              
##  [40] Rcpp_1.0.5              viridisLite_0.3.0       laeken_0.5.1           
##  [43] htmlTable_2.1.0         tmvnsim_1.0-2           units_0.6-7            
##  [46] foreign_0.8-79          vcd_1.4-8               htmlwidgets_1.5.2      
##  [49] httr_1.4.2              RColorBrewer_1.1-2      geosphere_1.5-10       
##  [52] ellipsis_0.3.1          reshape_0.8.8           farver_2.0.3           
##  [55] pkgconfig_2.0.3         XML_3.99-0.5            nnet_7.3-14            
##  [58] dbplyr_1.4.4            labeling_0.4.2          tidyselect_1.1.0       
##  [61] rlang_0.4.8             tmaptools_3.1           multilevel_2.6         
##  [64] munsell_0.5.0           cellranger_1.1.0        tools_4.0.3            
##  [67] cli_2.1.0               generics_0.0.2          ranger_0.12.1          
##  [70] broom_0.5.6             evaluate_0.14           yaml_2.2.1             
##  [73] leafsync_0.1.0          knitr_1.30              fs_1.5.0               
##  [76] zip_2.1.1               robustbase_0.93-6       nlme_3.1-150           
##  [79] quantreg_5.74           xml2_1.3.2              psychometric_2.2       
##  [82] compiler_4.0.3          rstudioapi_0.11         curl_4.3               
##  [85] png_0.1-7               e1071_1.7-4             reprex_0.3.0           
##  [88] stringi_1.5.3           classInt_0.4-3          vctrs_0.3.4            
##  [91] stringdist_0.9.6.3      pillar_1.4.6            lifecycle_0.2.0        
##  [94] lmtest_0.9-38           data.table_1.13.2       conquer_1.0.2          
##  [97] raster_3.3-13           R6_2.4.1                latticeExtra_0.6-29    
## [100] KernSmooth_2.23-18      gridExtra_2.3           rio_0.5.16             
## [103] codetools_0.2-17        polspline_1.1.19        dichromat_2.0-0        
## [106] boot_1.3-25             MASS_7.3-53             gtools_3.8.2           
## [109] spatialstatstools_0.1.0 withr_2.3.0             mnormt_2.0.2           
## [112] multcomp_1.4-14         mgcv_1.8-33             parallel_4.0.3         
## [115] hms_0.5.3               grid_4.0.3              rpart_4.1-15           
## [118] class_7.3-17            rmarkdown_2.5           carData_3.0-4          
## [121] lubridate_1.7.9         base64enc_0.1-3

#data out
d %>% write_rds("data/data_out.rds", compress = "xz")
d %>% write_csv("data/data_out.csv.gz")

#OSF repo update
if (F) {
  library(osfr)
  osf_auth(read_lines("~/.config/osf_token"))
  osf_proj = osf_retrieve_node("https://osf.io/pagz7/")
  osf_upload(osf_proj, path = c("notebook.html", "notebook.Rmd", "data", "figs", "Time preferences national IQs.Rproj"), conflicts = "overwrite", recurse = T)
}