The effect of legislation on sex ratio at birth

This is an example of using EventStudies on long-form data from social sciences. Specifically, we're going to be using data which

# Load the packages:
using TSFrames         # time series package - EventStudies only accepts this, since it has a defined index type
using EventStudies     # event study package
using DataFrames       # data manipulation
using CairoMakie       # plotting

Load the data.

data = EventStudies.load_data("bigdata.csv")

<div><div style = "float: left;"><span>16992×11 DataFrame</span></div><div style = "float: right;"><span style = "font-style: italic;">16982 rows omitted</span></div><div style = "clear: both;"></div></div><div class = "data-frame" style = "overflow-x: scroll;"><table class = "data-frame" style = "margin-bottom: 6px;"><thead><tr class = "header"><th class = "rowNumber" style = "font-weight: bold; text-align: right;">Row</th><th style = "text-align: left;">country</th><th style = "text-align: left;">iso_code</th><th style = "text-align: left;">year</th><th style = "text-align: left;">totpop</th><th style = "text-align: left;">sex_rat</th><th style = "text-align: left;">birthrate</th><th style = "text-align: left;">intyear</th><th style = "text-align: left;">year_es</th><th style = "text-align: left;">intyear_es</th><th style = "text-align: left;">agrigdp_rat</th><th style = "text-align: left;">gender_inequality_index</th></tr><tr class = "subheader headerLastRow"><th class = "rowNumber" style = "font-weight: bold; text-align: right;"></th><th title = "String" style = "text-align: left;">String</th><th title = "InlineStrings.String3" style = "text-align: left;">String3</th><th title = "Dates.Date" style = "text-align: left;">Date</th><th title = "Float64" style = "text-align: left;">Float64</th><th title = "Float64" style = "text-align: left;">Float64</th><th title = "Float64" style = "text-align: left;">Float64</th><th title = "Union{Missing, Dates.Date}" style = "text-align: left;">Date</th><th title = "Int64" style = "text-align: left;">Int64</th><th title = "Union{Missing, Int64}" style = "text-align: left;">Int64?</th><th title = "Union{Missing, Float64}" style = "text-align: left;">Float64?</th><th title = "Union{Missing, Float64}" style = "text-align: left;">Float64?</th></tr></thead><tbody><tr><td class = "rowNumber" style = "font-weight: bold; text-align: right;">1</td><td style = "text-align: left;">Burundi</td><td style = "text-align: left;">BDI</td><td style = "text-align: left;">1950-01-01</td><td style = "text-align: right;">2229.32</td><td style = "text-align: right;">102.5</td><td style = "text-align: right;">51.334</td><td style = "font-style: italic; text-align: left;">missing</td><td style = "text-align: right;">1950</td><td style = "font-style: italic; text-align: right;">missing</td><td style = "font-style: italic; text-align: right;">missing</td><td style = "font-style: italic; text-align: right;">missing</td></tr><tr><td class = "rowNumber" style = "font-weight: bold; text-align: right;">2</td><td style = "text-align: left;">Burundi</td><td style = "text-align: left;">BDI</td><td style = "text-align: left;">1951-01-01</td><td style = "text-align: right;">2278.9</td><td style = "text-align: right;">102.5</td><td style = "text-align: right;">50.678</td><td style = "font-style: italic; text-align: left;">missing</td><td style = "text-align: right;">1951</td><td style = "font-style: italic; text-align: right;">missing</td><td style = "font-style: italic; text-align: right;">missing</td><td style = "font-style: italic; text-align: right;">missing</td></tr><tr><td class = "rowNumber" style = "font-weight: bold; text-align: right;">3</td><td style = "text-align: left;">Burundi</td><td style = "text-align: left;">BDI</td><td style = "text-align: left;">1952-01-01</td><td style = "text-align: right;">2327.59</td><td style = "text-align: right;">102.5</td><td style = "text-align: right;">49.933</td><td style = "font-style: italic; text-align: left;">missing</td><td style = "text-align: right;">1952</td><td style = "font-style: italic; text-align: right;">missing</td><td style = "font-style: italic; text-align: right;">missing</td><td style = "font-style: italic; text-align: right;">missing</td></tr><tr><td class = "rowNumber" style = "font-weight: bold; text-align: right;">4</td><td style = "text-align: left;">Burundi</td><td style = "text-align: left;">BDI</td><td style = "text-align: left;">1953-01-01</td><td style = "text-align: right;">2375.48</td><td style = "text-align: right;">102.5</td><td style = "text-align: right;">49.448</td><td style = "font-style: italic; text-align: left;">missing</td><td style = "text-align: right;">1953</td><td style = "font-style: italic; text-align: right;">missing</td><td style = "font-style: italic; text-align: right;">missing</td><td style = "font-style: italic; text-align: right;">missing</td></tr><tr><td class = "rowNumber" style = "font-weight: bold; text-align: right;">5</td><td style = "text-align: left;">Burundi</td><td style = "text-align: left;">BDI</td><td style = "text-align: left;">1954-01-01</td><td style = "text-align: right;">2422.72</td><td style = "text-align: right;">102.5</td><td style = "text-align: right;">48.883</td><td style = "font-style: italic; text-align: left;">missing</td><td style = "text-align: right;">1954</td><td style = "font-style: italic; text-align: right;">missing</td><td style = "font-style: italic; text-align: right;">missing</td><td style = "font-style: italic; text-align: right;">missing</td></tr><tr><td style = "text-align: right;">&vellip;</td><td style = "text-align: right;">&vellip;</td><td style = "text-align: right;">&vellip;</td><td style = "text-align: right;">&vellip;</td><td style = "text-align: right;">&vellip;</td><td style = "text-align: right;">&vellip;</td><td style = "text-align: right;">&vellip;</td><td style = "text-align: right;">&vellip;</td><td style = "text-align: right;">&vellip;</td><td style = "text-align: right;">&vellip;</td><td style = "text-align: right;">&vellip;</td><td style = "text-align: right;">&vellip;</td></tr><tr><td class = "rowNumber" style = "font-weight: bold; text-align: right;">16988</td><td style = "text-align: left;">Wallis and Futuna Islands</td><td style = "text-align: left;">WLF</td><td style = "text-align: left;">2017-01-01</td><td style = "text-align: right;">12.002</td><td style = "text-align: right;">104.8</td><td style = "text-align: right;">12.713</td><td style = "font-style: italic; text-align: left;">missing</td><td style = "text-align: right;">2017</td><td style = "font-style: italic; text-align: right;">missing</td><td style = "font-style: italic; text-align: right;">missing</td><td style = "font-style: italic; text-align: right;">missing</td></tr><tr><td class = "rowNumber" style = "font-weight: bold; text-align: right;">16989</td><td style = "text-align: left;">Wallis and Futuna Islands</td><td style = "text-align: left;">WLF</td><td style = "text-align: left;">2018-01-01</td><td style = "text-align: right;">11.87</td><td style = "text-align: right;">104.8</td><td style = "text-align: right;">12.381</td><td style = "font-style: italic; text-align: left;">missing</td><td style = "text-align: right;">2018</td><td style = "font-style: italic; text-align: right;">missing</td><td style = "font-style: italic; text-align: right;">missing</td><td style = "font-style: italic; text-align: right;">missing</td></tr><tr><td class = "rowNumber" style = "font-weight: bold; text-align: right;">16990</td><td style = "text-align: left;">Wallis and Futuna Islands</td><td style = "text-align: left;">WLF</td><td style = "text-align: left;">2019-01-01</td><td style = "text-align: right;">11.761</td><td style = "text-align: right;">104.8</td><td style = "text-align: right;">12.189</td><td style = "font-style: italic; text-align: left;">missing</td><td style = "text-align: right;">2019</td><td style = "font-style: italic; text-align: right;">missing</td><td style = "font-style: italic; text-align: right;">missing</td><td style = "font-style: italic; text-align: right;">missing</td></tr><tr><td class = "rowNumber" style = "font-weight: bold; text-align: right;">16991</td><td style = "text-align: left;">Wallis and Futuna Islands</td><td style = "text-align: left;">WLF</td><td style = "text-align: left;">2020-01-01</td><td style = "text-align: right;">11.667</td><td style = "text-align: right;">104.8</td><td style = "text-align: right;">11.997</td><td style = "font-style: italic; text-align: left;">missing</td><td style = "text-align: right;">2020</td><td style = "font-style: italic; text-align: right;">missing</td><td style = "font-style: italic; text-align: right;">missing</td><td style = "font-style: italic; text-align: right;">missing</td></tr><tr><td class = "rowNumber" style = "font-weight: bold; text-align: right;">16992</td><td style = "text-align: left;">Wallis and Futuna Islands</td><td style = "text-align: left;">WLF</td><td style = "text-align: left;">2021-01-01</td><td style = "text-align: right;">11.642</td><td style = "text-align: right;">104.8</td><td style = "text-align: right;">11.816</td><td style = "font-style: italic; text-align: left;">missing</td><td style = "text-align: right;">2021</td><td style = "font-style: italic; text-align: right;">missing</td><td style = "font-style: italic; text-align: right;">missing</td><td style = "font-style: italic; text-align: right;">missing</td></tr></tbody></table></div>

This data is in "long" format, so we need to convert it to "wide" format.

We do this by grouping by country,

country_data = groupby(data, :country)

GroupedDataFrame with 236 groups based on key: country

⋮

filtering to only include countries with events defined,

countries_with_events = filter(sdf -> !(all(ismissing.(sdf.intyear_es))), country_data)

GroupedDataFrame with 15 groups based on key: country

⋮

and finally, creating a TSFrame with the sex_rat variable for each country.

relevant_country_ts = TSFrames.join((TSFrame(DataFrame([:Index => sdf.year_es, Symbol(sdf.country[1]) => sdf.sex_rat])) for sdf in countries_with_events)...; jointype = :JoinInner)

72×15 TSFrame with Int64 Index
 Index  Libya    China    China, Hong Kong SAR  China, Taiwan Province of Chin ⋯
 Int64  Float64  Float64  Float64               Float64                        ⋯
────────────────────────────────────────────────────────────────────────────────
  1950    105.7    106.1                 109.8                            105. ⋯
  1951    105.7    106.1                 109.7                            105.
  1952    105.7    106.1                 109.4                            105.
  1953    105.7    106.1                 109.2                            105.
  1954    105.6    106.1                 108.9                            105. ⋯
  1955    105.6    106.1                 108.9                            105.
  1956    105.6    106.1                 108.3                            105.
  1957    105.6    106.1                 109.0                            105.
   ⋮       ⋮        ⋮              ⋮                           ⋮               ⋱
  2015    105.7    115.0                 107.6                            107. ⋯
  2016    105.7    114.5                 107.6                            107.
  2017    105.7    114.0                 107.6                            107.
  2018    105.7    113.4                 107.7                            107.
  2019    105.7    112.9                 107.7                            107. ⋯
  2020    105.7    112.3                 107.7                            107.
  2021    105.7    111.8                 107.7                            107.
                                                  12 columns and 57 rows omitted

Create a mapping from country to event time, for each event.

Note that this can have multiple events per country.

event_times = [Symbol(sdf.country[1]) => sdf.intyear_es[findfirst(!ismissing, sdf.intyear_es)] for sdf in countries_with_events]

15-element Vector{Pair{Symbol, Int64}}:
                                    :Libya => 1986
                                    :China => 1983
            Symbol("China, Hong Kong SAR") => 1998
 Symbol("China, Taiwan Province of China") => 1974
               Symbol("Republic of Korea") => 1974
                                    :India => 1972
                                 :Pakistan => 1973
                        Symbol("Viet Nam") => 1999
                                  :Armenia => 1991
                               :Azerbaijan => 1991
                                  :Georgia => 1991
              Symbol("Russian Federation") => 1985
                                  :Albania => 1991
                                    :Spain => 1975
                                     :Cuba => 1990

Convert to event time using EventStudies.jl (perform the event study).

eventtime_ts, retcodes = EventStudies.eventstudy(levels_to_returns(relevant_country_ts), event_times, 4)

(9×15 TSFrame with Int64 Index
 Index  Libya       China      China, Hong Kong SAR  China, Taiwan Province of ⋯
 Int64  Float64     Float64    Float64               Float64                   ⋯
────────────────────────────────────────────────────────────────────────────────
    -4   0.0951928  0.0935016             0.371402                        -0.1 ⋯
    -3   0.0951022  0.0934143            -0.0927214                       -0.0
    -2   0.0950119  0.186567              0.0                             -0.0
    -1   0.284495   0.18622               0.46275                          0.3
     0   0.189215   0.278681              0.276625                         0.0 ⋯
     1   0.283153   0.277907              0.0                             -0.2
     2   0.0942063  0.461468              0.0920387                        0.0
     3   0.0        0.459349             -0.276371                         0.1
     4  -0.0942063  0.548447             -0.277136                         0.4 ⋯
                                                              12 columns omitted, EventStudies.EventStatus[EventStudies.Success(), EventStudies.Success(), EventStudies.Success(), EventStudies.Success(), EventStudies.Success(), EventStudies.Success(), EventStudies.Success(), EventStudies.Success(), EventStudies.Success(), EventStudies.Success(), EventStudies.Success(), EventStudies.Success(), EventStudies.Success(), EventStudies.Success(), EventStudies.Success()])

Perform bootstrap inference with a 95% confidence interval

t0, lower, upper = inference(BootstrapInference(), eventtime_ts, 0.95)

([0.08750559901358794, 0.05639795220573044, 0.037076052234521484, 0.14915962548207298, 0.2053603997918548, 0.08716474148950344, 0.30419233442415816, 0.27634298013774483, 0.3610686693636171], [0.012666142999793852, -0.04002900293414792, -0.03142136506890333, 0.04800203081096841, 0.09961761892883336, -0.08669979665881751, 0.13652946052636983, 0.06262873187748971, 0.15695636073234576], [0.16239711651013033, 0.1561040473290407, 0.1361553928509711, 0.25479841834865685, 0.3555889811206953, 0.24829626602792015, 0.5212877911822056, 0.5702681228666144, 0.5940640931676833])

we're technically done with the event study here!

now, it's time to plot plot each variable in the event study

fig, ax, plt = series(
    index(eventtime_ts), Matrix(remap_cumsum(eventtime_ts))';
    labels = names(eventtime_ts),
    color = Makie.resample_cmap(:rainbow_bgyrm_35_85_c71_n256, 15),
    linewidth = 3,
    axis = (
        xlabel = "Years relative to event",
        xticks = WilkinsonTicks(8, k_min = 8, k_max = 8),
        ylabel = "Cumulative change\nin sex ratio (%)",
        title = "Variable: sex ratio at birth",
        subtitle = "Event: some legislation",
        titlealign = :left,
    ),
    figure = (resolution = (1000, 1000),)
)
leg = Legend(fig[2, 1], ax; nbanks = 3, tellwidth = false, tellheight = true)
fig

plot inference

fig, ax, mean_plt = lines(
    index(eventtime_ts), t0;
    label = "Mean",
    linewidth = 3,
    axis = (
        xlabel = "Years relative to event",
        xticks = WilkinsonTicks(8, k_min = 8, k_max = 8),
        ylabel = "Cumulative change\nin sex ratio (%)",
        title = "Variable: sex ratio at birth",
        subtitle = "Event: some legislation",
        titlealign = :left,
    ),
)

band_plt = band!(ax,
    index(eventtime_ts), lower, upper;
    label = "95% CI",
    linewidth = 3,
    color = Makie.wong_colors(0.5)[2],
)

leg = axislegend(ax, position = :lt)

fig

plot inference using errorbars

fig, ax, plt = rangebars(
    index(eventtime_ts),
    lower, upper;
    whiskerwidth = 16,
    label = "95% CI",
    axis = (
        xlabel = "Years relative to event",
        xticks = WilkinsonTicks(8, k_min = 8, k_max = 8),
        ylabel = "Cumulative change\nin sex ratio (%)",
        title = "Variable: sex ratio at birth",
        subtitle = "Event: some legislation",
        titlealign = :left,
    ),
)

scatter!(ax, index(eventtime_ts), t0; label = "Mean value")

Scatter{Tuple{Vector{Point{2, Float32}}}}

small hack for a feature which we should add into Makie

function Makie.legendelements(plot::Rangebars, legend)
    line_points = lift(plot.whiskerwidth) do ww
        if ww ≤ 3
            Point2f[(0.5, 1), (0.5, 0)]
        else
            Point2f[(0.25, 1), (0.75, 1), (NaN, NaN), (0.5, 1), (0.5, 0), (NaN, NaN), (0.25, 0), (0.75, 0)]
        end
    end
    return Makie.LegendElement[
        LineElement(
            points = line_points,
            color = Makie.scalar_lift(plot.color, legend.linecolor),
            linewidth = Makie.scalar_lift(plot.linewidth, legend.linewidth)
        ),
    ]
end

axislegend(ax, position = :lt)

fig

replicate eventstudies.R results

convert to event time

eventtime_ts, retcodes = EventStudies.eventstudy(levels_to_returns(relevant_country_ts), event_times, -2:3)

(6×15 TSFrame with Int64 Index
 Index  Libya      China     China, Hong Kong SAR  China, Taiwan Province of C ⋯
 Int64  Float64    Float64   Float64               Float64                     ⋯
────────────────────────────────────────────────────────────────────────────────
    -2  0.0950119  0.186567             0.0                             -0.094 ⋯
    -1  0.284495   0.18622              0.46275                          0.375
     0  0.189215   0.278681             0.276625                         0.0
     1  0.283153   0.277907             0.0                             -0.281
     2  0.0942063  0.461468             0.0920387                        0.094 ⋯
     3  0.0        0.459349            -0.276371                         0.187
                                                              12 columns omitted, EventStudies.EventStatus[EventStudies.Success(), EventStudies.Success(), EventStudies.Success(), EventStudies.Success(), EventStudies.Success(), EventStudies.Success(), EventStudies.Success(), EventStudies.Success(), EventStudies.Success(), EventStudies.Success(), EventStudies.Success(), EventStudies.Success(), EventStudies.Success(), EventStudies.Success(), EventStudies.Success()])

perform inference with a 95% confidence interval

t0, lower, upper = inference(BootstrapInference(), eventtime_ts, 0.95)

([0.037076052234521484, 0.14915962548207298, 0.2053603997918548, 0.08716474148950344, 0.30419233442415816, 0.27634298013774483], [-0.0315635911707958, 0.03701683035704321, 0.09358506164096725, -0.0749783996219527, 0.13618171162839152, 0.04323068287054096], [0.12989470094049999, 0.2497554490344732, 0.34647341283161553, 0.24032409540363, 0.5209548982132289, 0.5761557157510687])

we're technically done with the event study here!

plot inference

fig, ax, mean_plt = lines(
    index(eventtime_ts), t0;
    label = "Mean",
    linewidth = 3,
    axis = (
        xlabel = "Years relative to event",
        xticks = WilkinsonTicks(8, k_min = 8, k_max = 8),
        ylabel = "Cumulative change\nin sex ratio (%)",
        title = "Variable: sex ratio at birth",
        subtitle = "Event: some legislation",
        titlealign = :left,
    ),
)

band_plt = band!(ax,
    index(eventtime_ts), lower, upper;
    label = "95% CI",
    linewidth = 3,
    color = Makie.wong_colors(0.5)[2],
)

leg = axislegend(ax, position = :lt)

ylims!(ax, -0.75, 0.75) # match R
ylims!(ax, -0.5, 0.5) # match R
fig

This page was generated using Literate.jl.