European Parenting Leave Policies

European Parenting Leave Policies

This week we’re exploring European Parenting Leave Policies. The European Parenting Leave Policies (EPLP) Dataset provides harmonised data on maternity, co-parent, paid parental, and job-protected leave regulations across 21 European countries from 1970 to 2024.

The dataset enables quantitative analyses of policy trends, cross-national differences, and the effects of major reforms – for researchers, policymakers, and others interested in family policy.

Given the variety of parental leave schemes across countries, the dataset considers three different dimensions of parental leave duration for each country, if applicable. Dimension 1 (par1) identifies the paid parental leave scheme with the longest possible duration. Dimension 2 (par2) identifies the paid parental leave duration with the highest monthly flat rate payment. Dimension 3 (par3) identifies the duration with the highest replacement rate.

Values that are missing are represented as NA. Some values are missing because they are not applicable. These values are encoded as "Not applicable" for character vectors, and -98 for numeric variables.

• Which countries were the first to implement co-parent leave policies?
• Has parenting leave decreased in any countries?

Cite the dataset as: S. Spitzer et al., “The European Parenting Leave Policies (EPLP) Dataset”. Zenodo, Nov. 19, 2025. doi: 10.5281/zenodo.17648712.

Thank you to Nicola Rennie for curating this week’s dataset.

library(tidyverse)
library(sf)
library(rnaturalearth)
library(gganimate)
library(magick)
library(praise)
library(broom)

eplp <- readr::read_csv('https://raw.githubusercontent.com/rfordatascience/tidytuesday/main/data/2026/2026-06-02/eplp.csv') |>
  mutate(across(where(is.numeric), ~ na_if(.x, -98)))  # changes all of the -98 to NA

Job-protected leave

Using the variable job-protected leave duration for birth mothers (maximum duration), we consider changes across time and country. The variable measures the maximum length of time, after childbirth, that a birth mother is legally entitled to be absent from paid work while her job (or an equivalent position) is protected.

europe <- ne_countries(continent = "europe", returnclass = "sf")
yrs <- sort(unique(eplp$year))

p <- europe |> 
  mutate(country2 = ifelse(iso_a2_eh == "GB", "UK", iso_a2_eh)) |>
  tidyr::expand_grid(year = yrs) |>
  st_as_sf() |>
  left_join(eplp, by = c("country2" = "country", "year")) |>
  ggplot() + 
  geom_sf(aes(fill = jp_ld_m)) + 
  coord_sf(xlim = c(-10, 45), ylim = c(27, 82)) + 
  theme_minimal() + 
  scale_fill_gradient2(
    low = "#91BFDB", mid = "#FFFFBF", high = "#FC8D59",
    name = "job protected \nleave (weeks)",
    midpoint = median(eplp$jp_ld_m, na.rm = TRUE)
  ) + 
  transition_states(year, transition_length = 0, state_length = 4) + 
  labs(title = "Year: {closest_state}")

animate(p)

Over time most countries have increased their number of weeks of job protection.

Predicting protected leave weeks

We can use a linear model and the following variables - par1_ld (paid parental leave duration ) and year and par1_for_whom (who is eligible for paid parental leave) - to predict jd_ld_m (job-protected leave duration for birth mothers ).

Caveat: from the scatterplot (and the residual plot below), it doesn’t seem like a linear model is necessarily the best model to fit here. That is, jd_ld_m is not really a linear model of the variables.

Running models one at a time to investigate the relationship between part1_ld and par1_for_whom in predicting jp_ld_m. We are only interested in predicting for those countries that actually have policies (i.e., where part1_ld \(\ne 0\)). Note that par1_for_whom is either mothers or either.

eplp1 <- eplp |>
  filter(par1_ld != 0)

eplp1 |>
  ggplot(aes(x = par1_ld, y = jp_ld_m, color = year)) +
  geom_point() + 
  geom_smooth(method = "lm", se = FALSE, color = "red") + 
  facet_wrap( ~ par1_for_whom)

Again, with the caveat that the relationship is not truly linear, it seems like both par1_ld and par1_for_whom are important in predicting jp_ld_m. Additionally, there seems to be moderate interaction between par1_ld and par1_for_whom.

lm(jp_ld_m ~ par1_ld, data = eplp1) |>
  tidy()

# A tibble: 2 × 5
  term        estimate std.error statistic   p.value
  <chr>          <dbl>     <dbl>     <dbl>     <dbl>
1 (Intercept)   46.5      2.15        21.6 3.07e- 79
2 par1_ld        0.754    0.0281      26.8 7.65e-109

lm(jp_ld_m ~ par1_ld + par1_for_whom, data = eplp1) |>
  tidy()

# A tibble: 3 × 5
  term                 estimate std.error statistic   p.value
  <chr>                   <dbl>     <dbl>     <dbl>     <dbl>
1 (Intercept)            50.9      2.30       22.1  3.34e- 82
2 par1_ld                 0.742    0.0277     26.8  1.75e-108
3 par1_for_whommothers  -14.4      2.94       -4.91 1.17e-  6

lm(jp_ld_m ~ par1_ld * par1_for_whom, data = eplp1) |>
  tidy()

# A tibble: 4 × 5
  term                         estimate std.error statistic  p.value
  <chr>                           <dbl>     <dbl>     <dbl>    <dbl>
1 (Intercept)                    52.5      2.41       21.8  3.28e-80
2 par1_ld                         0.717    0.0300     23.9  3.06e-92
3 par1_for_whommothers          -23.6      5.22       -4.53 6.93e- 6
4 par1_ld:par1_for_whommothers    0.165    0.0774      2.13 3.33e- 2

To check whether it was reasonable to run a linear model, we produce a residual plot. As expected, the residuals are not particularly scattered (normally) around the y = 0 line. So any inference drawn on the model should be done so with a grain of salt.

lm(jp_ld_m ~ par1_ld * par1_for_whom, data = eplp1) |>
  augment() |>
  ggplot(aes(x = .fitted, y = .resid)) + 
  geom_point() + 
  geom_hline(yintercept = 0)

praise()

[1] "You are pioneering!"