Init

library(kirkegaard)
## Loading required package: tidyverse
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.1     ✔ stringr   1.5.2
## ✔ ggplot2   4.0.0     ✔ tibble    3.3.0
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.1
## ✔ purrr     1.1.0     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
## Loading required package: magrittr
## 
## 
## Attaching package: 'magrittr'
## 
## 
## The following object is masked from 'package:purrr':
## 
##     set_names
## 
## 
## The following object is masked from 'package:tidyr':
## 
##     extract
## 
## 
## Loading required package: weights
## 
## Loading required package: assertthat
## 
## 
## Attaching package: 'assertthat'
## 
## 
## The following object is masked from 'package:tibble':
## 
##     has_name
## 
## 
## Loading required package: psych
## 
## 
## Attaching package: 'psych'
## 
## 
## The following objects are masked from 'package:ggplot2':
## 
##     %+%, alpha
## 
## 
## Loading required package: robustbase
## 
## 
## Attaching package: 'kirkegaard'
## 
## 
## The following object is masked from 'package:psych':
## 
##     rescale
## 
## 
## The following object is masked from 'package:assertthat':
## 
##     are_equal
## 
## 
## The following object is masked from 'package:purrr':
## 
##     is_logical
## 
## 
## The following object is masked from 'package:base':
## 
##     +
load_packages(
  rnaturalearth, rnaturalearthdata,
  sf,
  fixest,
  ggeffects,
  readxl,
  broom
)
## 
## Attaching package: 'rnaturalearthdata'
## 
## The following object is masked from 'package:rnaturalearth':
## 
##     countries110
## 
## Linking to GEOS 3.10.2, GDAL 3.4.1, PROJ 8.2.1; sf_use_s2() is TRUE
theme_set(theme_bw())

options(
    digits = 3
)

#multithreading
#library(future)
#plan(multisession(workers = 8))

Functions

Data

#fertility data
#https://ourworldindata.org/fertility-rate
d_tfr = read_csv("data/children-born-per-woman/children-born-per-woman.csv") %>% 
  rename(
    year = Year,
    ISO3 = Code, 
    TFR = `Fertility rate (period), historical`
  ) %>% 
  select(
    -Entity
  ) %>% 
  miss_filter()
## Rows: 18958 Columns: 4
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (2): Entity, Code
## dbl (2): Year, Fertility rate (period), historical
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#spending data
#https://data-explorer.oecd.org/vis?lc=en&tm=family&pg=0&hc[Programme%20type]=Family&snb=63&vw=tb&df[ds]=dsDisseminateFinalDMZ&df[id]=DSD_SOCX_AGG%40DF_PUB_FAM&df[ag]=OECD.ELS.SPD&df[vs]=&pd=%2C&dq=.A..PT_B1GQ..C%2BK%2B_T.TP51.&ly[rw]=REF_AREA&ly[cl]=TIME_PERIOD&ly[rs]=SPENDING_TYPE&to[TIME_PERIOD]=false&lb=bt
d_spending = read_csv("data/OECD.ELS.SPD,DSD_SOCX_AGG@DF_PUB_FAM,+.A..PT_B1GQ..C+K+_T.TP51..csv") %>% 
  select(
    ISO3 = REF_AREA,
    spending = OBS_VALUE,
    spending_type = `Spending type`,
    year = TIME_PERIOD,
  ) %>% 
  miss_filter()
## Rows: 5478 Columns: 34
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (26): STRUCTURE, STRUCTURE_ID, STRUCTURE_NAME, ACTION, REF_AREA, Referen...
## dbl  (4): TIME_PERIOD, OBS_VALUE, UNIT_MULT, DECIMALS
## lgl  (4): Time period, Observation value, BASE_PER, Base period
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#excel file
d_spending2 = read_excel("data/PF1_1_Public_spending_on_family_benefits.xlsx", sheet = 3) %>% 
  tidyr::fill(Country) %>% 
  mutate(
    ISO3 = pu_translate(Country),
    spending_type = if_else(is.na(spending_type), "Total", spending_type)
  ) %>% 
  select(-Country) %>%
  pivot_longer(
    cols = -c(ISO3, spending_type),
    names_to = "year",
    values_to = "spending"
  ) %>% 
  mutate(
    spending = as.numeric(spending),
    year = as.numeric(year)
  )
## No exact match: Türkiye
## Best fuzzy match found: Türkiye -> Turkije with distance 2.00
## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `spending = as.numeric(spending)`.
## Caused by warning:
## ! NAs introduced by coercion
#get map data for countries
d_world <- ne_countries(scale = "medium", returnclass = "sf")
d_world$ISO3 = d_world$adm0_a3

#join
d = inner_join(
  d_tfr,
  d_spending,
  by = c("ISO3", "year")
)

#mutate
d %<>% mutate(
  spending_type = fct_relevel(spending_type, "Total")
)

#add names back
d$name = pu_translate(d$ISO3, reverse = T)

#join spatial data
d_world_data = left_join(
  d_world,
  d,
  by = "ISO3"
)

#join again with alt dataset
d2 = inner_join(
  d_tfr,
  d_spending2,
  by = c("ISO3", "year")
)

#mutate
d2$name = pu_translate(d2$ISO3, reverse = T)

Analysis

#map
d_world_data %>% 
  #set year to 2019, so rows aren't dropped
  mutate(
    year = if_else(is.na(year), 2019, year),
    spending_type = if_else(is.na(spending_type), "Total", spending_type)
  ) %>% 
  filter(
    year == 2019,
    spending_type == "Total"
  ) %>% 
  ggplot() +
  geom_sf(aes(fill = spending)) +
  scale_fill_viridis_c(option = "plasma") +
  labs(
    title = "Family spending of GDP % in 2019",
    fill = "Family spending of GDP %"
  )

GG_save("figs/family_spending_gdp_2019_map.png")

#scatter plot for 2019
d %>% 
  filter(
    year == 2019,
    spending_type == "Total"
    ) %>% 
  GG_scatter("spending", "TFR", case_names = "name") +
  labs(
    title = "Public spending on family of GDP % vs total fertility rate (TFR) in 2019",
    x = "Family spending of GDP %",
    y = "Total fertility rate"
  )
## `geom_smooth()` using formula = 'y ~ x'

GG_save("figs/family_spending_gdp_vs_tfr_2019_scatter.png")
## `geom_smooth()` using formula = 'y ~ x'
#do regression for each year of data
d_reg = d %>% 
  filter(
    spending_type == "Total"
  ) %>% 
  group_by(year) %>% 
  group_modify(~{
    mod = lm(TFR ~ spending, data = .x)
    
    r_test = cor.test(.x$TFR, .x$spending, use = "pairwise.complete")
    
    tibble(
      intercept = coef(mod)[1],
      slope = coef(mod)[2],
      slope_lower = confint(mod)[2,1],
      slope_upper = confint(mod)[2,2],
      r_squared = summary(mod)$adj.r.squared,
      r = r_test$estimate,
      r_lower = r_test$conf.int[1],
      r_upper = r_test$conf.int[2],
      p_value = r_test$p.value
    )
  })

#plot results
d_reg %>%
  ggplot(aes(year, r)) +
  geom_ribbon(aes(ymin = r_lower, ymax = r_upper), alpha = 0.2) +
  geom_line() +
  labs(
    title = "Correlation between public spending on family of GDP % and total fertility rate over time",
    x = "Year",
    y = "Pearson's r"
  )

GG_save("figs/family_spending_gdp_vs_tfr_correlation_over_time.png")

#slope
d_reg %>%
  ggplot(aes(year, slope)) +
  geom_ribbon(aes(ymin = slope_lower, ymax = slope_upper), alpha = 0.2) +
  geom_line() +
  labs(
    title = "Slope of family spending of GDP % vs total fertility rate over time",
    x = "Year",
    y = "Slope"
  )

#fixed effects
tfr_fixed = fixest::feols(
  TFR ~ spending | ISO3 + year,
  data = d %>% filter(spending_type == "Total")
)

summary(tfr_fixed)
## OLS estimation, Dep. Var.: TFR
## Observations: 1,439
## Fixed-effects: ISO3: 42,  year: 44
## Standard-errors: IID 
##          Estimate Std. Error t value Pr(>|t|) 
## spending   0.0186     0.0138    1.35  0.17763 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 0.207009     Adj. R2: 0.759279
##                  Within R2: 0.001343
#marginal effects
ggeffects::ggpredict(tfr_fixed, terms = "spending [all]") %>%
  plot()

#spending type interaction
tfr_fixed2 = fixest::feols(
  TFR ~ spending : spending_type | ISO3 + year,
  data = d %>% filter(spending_type != "Total") %>% mutate(
    spending_type = spending_type %>% fct_drop()
  )
)

summary(tfr_fixed2)
## OLS estimation, Dep. Var.: TFR
## Observations: 2,864
## Fixed-effects: ISO3: 42,  year: 44
## Standard-errors: IID 
##                                        Estimate Std. Error t value Pr(>|t|) 
## spending:spending_typeIn cash spending  0.00945    0.00753    1.26  0.20933 
## spending:spending_typeIn kind spending  0.01644    0.01179    1.39  0.16333 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 0.202093     Adj. R2: 0.763369
##                  Within R2: 7.568e-4
#with lags
#since pregnancy lasts ~9 months, we can try lags
d_lag = d %>%
  arrange(ISO3, year) %>%
  group_by(ISO3, spending_type) %>%
  mutate(
    spending_lag1 = lag(spending, 1),
    spending_lag2 = lag(spending, 2),
    spending_lag3 = lag(spending, 3)
  ) %>%
  ungroup()

tfr_fixed_lags = fixest::feols(
  TFR ~ spending + spending_lag1 + spending_lag2 + spending_lag3 | ISO3 + year,
  data = d_lag %>% filter(spending_type == "Total")
)
## NOTE: 126 observations removed because of NA values (RHS: 126).
summary(tfr_fixed_lags)
## OLS estimation, Dep. Var.: TFR
## Observations: 1,313
## Fixed-effects: ISO3: 42,  year: 41
## Standard-errors: IID 
##               Estimate Std. Error t value Pr(>|t|) 
## spending      -0.00347     0.0195  -0.178  0.85880 
## spending_lag1  0.01154     0.0235   0.491  0.62335 
## spending_lag2 -0.02529     0.0273  -0.925  0.35496 
## spending_lag3  0.04062     0.0273   1.489  0.13678 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 0.177711     Adj. R2: 0.794053
##                  Within R2: 0.002882

With alt dataset

d2 %>% 
  filter(
    year == 2021,
    spending_type == "Total"
    ) %>% 
  GG_scatter("spending", "TFR", case_names = "name") +
  labs(
    title = "Public spending on family of GDP % vs total fertility rate (TFR) in 2021",
    subtitle = "Using alternative dataset from OECD Excel file, includes tax breaks",
    x = "Family spending of GDP %",
    y = "Total fertility rate"
  )
## `geom_smooth()` using formula = 'y ~ x'

GG_save("figs/family_spending_gdp_vs_tfr_2021_scatter_altdata.png")
## `geom_smooth()` using formula = 'y ~ x'
#fixed effects again
tfr_fixed_alt = fixest::feols(
  TFR ~ spending | ISO3 + year,
  data = d2 %>% filter(spending_type == "Total")
)
## NOTE: 943 observations removed because of NA values (RHS: 943).
summary(tfr_fixed_alt)
## OLS estimation, Dep. Var.: TFR
## Observations: 653
## Fixed-effects: ISO3: 38,  year: 21
## Standard-errors: IID 
##          Estimate Std. Error t value Pr(>|t|) 
## spending   0.0179     0.0114    1.57  0.11707 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 0.110072     Adj. R2: 0.886721
##                  Within R2: 0.00413
#split by type
tfr_fixed2_alt = fixest::feols(
  TFR ~ spending : spending_type | ISO3 + year,
  data = d2 %>% filter(spending_type != "Total") %>% mutate(
    spending_type = spending_type %>% fct_drop()
  )
)
## NOTE: 1,442 observations removed because of NA values (RHS: 1,442).
summary(tfr_fixed2_alt)
## OLS estimation, Dep. Var.: TFR
## Observations: 3,346
## Fixed-effects: ISO3: 38,  year: 42
## Standard-errors: IID 
##                                               Estimate Std. Error t value
## spending:spending_typeCash                     0.00149    0.00581   0.256
## spending:spending_typeServices                 0.00418    0.00901   0.463
## spending:spending_typeTax-breaks for families  0.14030    0.03113   4.506
##                                                 Pr(>|t|)    
## spending:spending_typeCash                    7.9801e-01    
## spending:spending_typeServices                6.4312e-01    
## spending:spending_typeTax-breaks for families 6.8387e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 0.195266     Adj. R2: 0.767416
##                  Within R2: 0.006646
#just tax breaks
tfr_fixed_taxbreaks = fixest::feols(
  TFR ~ spending | ISO3 + year,
  data = d2 %>% filter(spending_type == "Tax-breaks for families")
)
## NOTE: 929 observations removed because of NA values (RHS: 929).
summary(tfr_fixed_taxbreaks)
## OLS estimation, Dep. Var.: TFR
## Observations: 667
## Fixed-effects: ISO3: 37,  year: 21
## Standard-errors: IID 
##          Estimate Std. Error t value   Pr(>|t|)    
## spending    0.213     0.0361     5.9 5.8761e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 0.105971     Adj. R2: 0.896283
##                  Within R2: 0.054149
#marginal effects
#bugged??
# ggeffects::ggpredict(tfr_fixed_taxbreaks, terms = "spending [all]") %>% 
#   plot()

#cross-sectional
d2 %>% 
  filter(
    year == 2021,
    spending_type == "Tax-breaks for families"
  ) %>% 
  GG_scatter("spending", "TFR", case_names = "name") +
  labs(
    title = "Public spending on family of GDP % vs total fertility rate (TFR) in 2021",
    subtitle = "Using alternative dataset from OECD Excel file, tax breaks only",
    x = "Family spending of GDP %, as tax reductions for families",
    y = "Total fertility rate"
  )
## `geom_smooth()` using formula = 'y ~ x'

GG_save("figs/family_spending_gdp_vs_tfr_2021_scatter_altdata_taxbreaks.png")
## `geom_smooth()` using formula = 'y ~ x'
#data check
d2 %>% 
  filter(spending_type == "Tax-breaks for families") %>%
  group_by(ISO3) %>%
  summarise(
    n_obs = sum(!is.na(spending)),
    mean_spending = mean(spending, na.rm=T),
    sd_spending = sd(spending, na.rm=T),
    range = max(spending, na.rm=T) - min(spending, na.rm=T)
  ) %>%
  filter(sd_spending > 0) %>%
  arrange(desc(sd_spending))
## Warning: There were 2 warnings in `summarise()`.
## The first warning was:
## ℹ In argument: `range = max(spending, na.rm = T) - min(spending, na.rm = T)`.
## ℹ In group 7: `ISO3 = "COL"`.
## Caused by warning in `max()`:
## ! no non-missing arguments to max; returning -Inf
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 1 remaining warning.
#residuals
# augment(
#   tfr_fixed2_alt, 
#   data = d2 %>% filter(spending_type != "Total") %>% mutate(spending_type = spending_type %>% fct_drop())
#   ) %>%
#   filter(spending_type == "Tax-breaks for families") %>%
#   ggplot(aes(x = spending, y = .resid)) +
#   geom_point() +
#   geom_smooth()

#lags and leads
#issues with colinearity here
feols(TFR ~ lead(spending, 1) + spending + 
            lag(spending, 1) | ISO3 + year,
      data = d2 %>% filter(spending_type == "Tax-breaks for families"))
## NOTES: 1,009 observations removed because of NA values (RHS: 1,009).
##        1/0 fixed-effect singleton was removed (1 observation).
## The variable 'lag(spending, 1)' has been removed because of collinearity (see
## $collin.var).
## OLS estimation, Dep. Var.: TFR
## Observations: 586
## Fixed-effects: ISO3: 36,  year: 20
## Standard-errors: IID 
##                   Estimate Std. Error t value   Pr(>|t|)    
## lead(spending, 1)    0.115     0.0385    2.99 0.00294375 ** 
## spending             0.157     0.0386    4.05 0.00005876 ***
## ... 1 variable was removed because of collinearity (lag(spending, 1))
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 0.096093     Adj. R2: 0.915766
##                  Within R2: 0.074038
#first differences
d2_tax_fd <- d2 %>%
  filter(spending_type == "Tax-breaks for families") %>%
  arrange(ISO3, year) %>%
  group_by(ISO3) %>%
  mutate(
    d_TFR = TFR - lag(TFR),
    d_spending = spending - lag(spending),
    d_spending_lag1 = lag(d_spending),
    d_spending_lag2 = lag(d_spending, 2),
    d_spending_lead1 = lead(d_spending),
    d_spending_lead2 = lead(d_spending, 2),
  ) %>%
  ungroup()

feols(d_TFR ~ d_spending + d_spending_lag1 + d_spending_lag2 + d_spending_lead1 + d_spending_lead2 | year + ISO3, 
      data = d2_tax_fd)
## NOTE: 1,155 observations removed because of NA values (LHS: 38, RHS: 1,155).
## OLS estimation, Dep. Var.: d_TFR
## Observations: 441
## Fixed-effects: year: 16,  ISO3: 35
## Standard-errors: IID 
##                  Estimate Std. Error t value Pr(>|t|)    
## d_spending        0.03941     0.0220  1.7937 0.073649 .  
## d_spending_lag1   0.03215     0.0221  1.4528 0.147082    
## d_spending_lag2   0.01851     0.0171  1.0825 0.279711    
## d_spending_lead1  0.01942     0.0234  0.8308 0.406605    
## d_spending_lead2  0.00201     0.0216  0.0931 0.925834    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 0.032244     Adj. R2: 0.349239
##                  Within R2: 0.01023
# Safer version with error handling
country_fd_results <- d2_tax_fd %>%
  filter(!is.na(d_TFR), !is.na(d_spending)) %>%
  group_by(ISO3) %>%
  filter(n() >= 10, sd(d_spending, na.rm = TRUE) > 0) %>%  # Exclude zero variation
  do(
    tidy(lm(d_TFR ~ d_spending, data = .), conf.int = TRUE)
  ) %>%
  filter(term == "d_spending") %>%
  arrange(p.value)

country_fd_results
#just the slope within each country
d2_tax_slopes <- d2 %>%
  filter(spending_type == "Tax-breaks for families") %>%
  group_by(ISO3) %>%
  filter(n() >= 10, sd(spending, na.rm = TRUE) > 0) %>%  # Exclude zero variation
  do(
    tidy(lm(TFR ~ spending + year, data = .), conf.int = TRUE)
  ) %>%
  filter(term != "(Intercept)") %>%
  arrange(p.value)
## Warning in qt(a, object$df.residual): NaNs produced
#plot
d2_tax_slopes %>%
  filter(
    term == "spending"
  ) %>% 
  ggplot(aes(x = reorder(ISO3, estimate), y = estimate)) +
  geom_point() +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.2) +
  coord_flip() +
  labs(
    title = "Slopes of TFR vs family spending of GDP % within countries",
    subtitle = "Tax breaks only, linear regresison within country, year as control",
    x = "Country (ISO3)",
    y = "Slope"
  )

GG_save("figs/family_spending_gdp_vs_tfr_slopes_within_countries.png")

#plot SVN
d2 %>%
  filter(
    ISO3 == "SVN",
    spending_type == "Tax-breaks for families"
    ) %>%
  ggplot(aes(spending, TFR)) +
  geom_point() +
  geom_line() +
  geom_text(aes(label = year)) +
  geom_smooth(method = "lm") +
  labs(
    title = "Slovenia (SVN) TFR vs family spending of GDP %, tax breaks only",
    x = "Family spending of GDP %",
    y = "Total fertility rate"
  )
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 21 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 21 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 21 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 21 rows containing missing values or values outside the scale range
## (`geom_text()`).

GG_save("figs/family_spending_gdp_vs_tfr_svn.png")
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 21 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 21 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 21 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 21 rows containing missing values or values outside the scale range
## (`geom_text()`).
#hungary
d2 %>%
  filter(
    ISO3 == "HUN",
    spending_type == "Tax-breaks for families"
  ) %>%
  ggplot(aes(spending, TFR)) +
  geom_point() +
  geom_line() +
  geom_text(aes(label = year)) +
  geom_smooth(method = "lm") +
  labs(
    title = "Hungary (HUN) TFR vs family spending of GDP %, tax breaks only",
    x = "Family spending of GDP %, tax breaks only",
    y = "Total fertility rate"
  )
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 21 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 21 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 21 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 21 rows containing missing values or values outside the scale range
## (`geom_text()`).

GG_save("figs/family_spending_gdp_vs_tfr_hun.png")
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 21 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 21 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 21 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 21 rows containing missing values or values outside the scale range
## (`geom_text()`).
#refit model, but exclude countries with very low variation
tfr_fixed2_alt_restricted = fixest::feols(
  TFR ~ spending | ISO3 + year,
  data = d2 %>% 
    filter(spending_type == "Tax-breaks for families") %>% 
    mutate(
      spending_type = spending_type %>% fct_drop()
    ) %>%
    group_by(ISO3) %>%
    filter(
      n() >= 10,
      sd(spending, na.rm = T) > 0.1
    ) %>%
    ungroup()
)
## NOTE: 356 observations removed because of NA values (RHS: 356).
summary(tfr_fixed2_alt_restricted)
## OLS estimation, Dep. Var.: TFR
## Observations: 274
## Fixed-effects: ISO3: 15,  year: 21
## Standard-errors: IID 
##          Estimate Std. Error t value  Pr(>|t|)    
## spending    0.156     0.0307    5.09 7.105e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 0.076283     Adj. R2: 0.95562 
##                  Within R2: 0.098342

Meta

#versions
write_sessioninfo()
## R version 4.5.1 (2025-06-13)
## Platform: x86_64-pc-linux-gnu
## Running under: Linux Mint 21.1
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0 
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0  LAPACK version 3.10.0
## 
## 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=en_DK.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_DK.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_DK.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: Europe/Brussels
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] broom_1.0.10            readxl_1.4.5            ggeffects_2.3.1        
##  [4] fixest_0.13.2           sf_1.0-21               rnaturalearthdata_1.0.0
##  [7] rnaturalearth_1.1.0     kirkegaard_2025-10-17   robustbase_0.99-6      
## [10] psych_2.5.6             assertthat_0.2.1        weights_1.1.2          
## [13] magrittr_2.0.4          lubridate_1.9.4         forcats_1.0.1          
## [16] stringr_1.5.2           dplyr_1.1.4             purrr_1.1.0            
## [19] readr_2.1.5             tidyr_1.3.1             tibble_3.3.0           
## [22] ggplot2_4.0.0           tidyverse_2.0.0        
## 
## loaded via a namespace (and not attached):
##   [1] RColorBrewer_1.1-3  rstudioapi_0.17.1   jsonlite_2.0.0     
##   [4] shape_1.4.6.1       datawizard_1.3.0    TH.data_1.1-4      
##   [7] estimability_1.5.1  jomo_2.7-6          farver_2.1.2       
##  [10] nloptr_2.2.1        rmarkdown_2.30      ragg_1.5.0         
##  [13] vctrs_0.6.5         minqa_1.2.8         base64enc_0.1-3    
##  [16] htmltools_0.5.8.1   haven_2.5.5         cellranger_1.1.0   
##  [19] Formula_1.2-5       mitml_0.4-5         sass_0.4.10        
##  [22] KernSmooth_2.23-26  bslib_0.9.0         htmlwidgets_1.6.4  
##  [25] sandwich_3.1-1      emmeans_1.11.2-8    zoo_1.8-14         
##  [28] cachem_1.1.0        lifecycle_1.0.4     iterators_1.0.14   
##  [31] pkgconfig_2.0.3     Matrix_1.7-4        R6_2.6.1           
##  [34] fastmap_1.2.0       rbibutils_2.3       digest_0.6.37      
##  [37] numDeriv_2016.8-1.1 colorspace_2.1-2    textshaping_1.0.4  
##  [40] Hmisc_5.2-4         labeling_0.4.3      timechange_0.3.0   
##  [43] gdata_3.0.1         mgcv_1.9-3          compiler_4.5.1     
##  [46] proxy_0.4-27        bit64_4.6.0-1       withr_3.0.2        
##  [49] htmlTable_2.4.3     S7_0.2.0            backports_1.5.0    
##  [52] DBI_1.2.3           pan_1.9             MASS_7.3-65        
##  [55] classInt_0.4-11     gtools_3.9.5        tools_4.5.1        
##  [58] units_1.0-0         foreign_0.8-90      nnet_7.3-20        
##  [61] glue_1.8.0          nlme_3.1-168        stringmagic_1.2.0  
##  [64] grid_4.5.1          stringdist_0.9.15   checkmate_2.3.3    
##  [67] cluster_2.1.8.1     generics_0.1.4      gtable_0.3.6       
##  [70] tzdb_0.5.0          class_7.3-23        data.table_1.17.8  
##  [73] hms_1.1.3           foreach_1.5.2       pillar_1.11.1      
##  [76] vroom_1.6.6         splines_4.5.1       lattice_0.22-7     
##  [79] survival_3.8-3      bit_4.6.0           dreamerr_1.5.0     
##  [82] tidyselect_1.2.1    knitr_1.50          reformulas_0.4.1   
##  [85] gridExtra_2.3       xfun_0.53           DEoptimR_1.1-4     
##  [88] stringi_1.8.7       yaml_2.3.10         boot_1.3-32        
##  [91] evaluate_1.0.5      codetools_0.2-20    archive_1.1.12     
##  [94] cli_3.6.5           rpart_4.1.24        systemfonts_1.3.1  
##  [97] xtable_1.8-4        Rdpack_2.6.4        jquerylib_0.1.4    
## [100] Rcpp_1.1.0          coda_0.19-4.1       parallel_4.5.1     
## [103] lme4_1.1-37         glmnet_4.1-10       viridisLite_0.4.2  
## [106] mvtnorm_1.3-3       scales_1.4.0        e1071_1.7-16       
## [109] insight_1.4.2       crayon_1.5.3        rlang_1.1.6        
## [112] multcomp_1.4-28     mnormt_2.1.1        mice_3.18.0
#write data to file for reuse
d %>% write_rds("data/data_for_reuse.rds")
d2 %>% write_rds("data/data_for_reuse_alt.rds")

#OSF
if (F) {
  library(osfr)
  
  #login
  osf_auth(readr::read_lines("~/.config/osf_token"))
  
  #the project we will use
  osf_proj = osf_retrieve_node("https://osf.io/XXX/")
  
  #upload all files in project
  #overwrite existing (versioning)
  osf_upload(
    osf_proj,
    path = c("data", "figures", "papers", "notebook.Rmd", "notebook.html", "sessions_info.txt"), 
    conflicts = "overwrite"
    )
}