Skip to content

Instantly share code, notes, and snippets.

@jseabold

jseabold/Difference_in_Differences_Differences.ipynb

Last active March 14, 2021 18:13
Show Gist options
  • Star 2 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save jseabold/42b39c4f3c5b7d6b1f73ac7749db2f3b to your computer and use it in GitHub Desktop.
Save jseabold/42b39c4f3c5b7d6b1f73ac7749db2f3b to your computer and use it in GitHub Desktop.
Download ZIP
Working out some differences between R and Python for Causal Inference: The Mixtape
Display the source blob
Display the rendered blob
Raw
Difference_in_Differences_Differences.ipynb
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Difference-in-Differences\n",
"### R and Python Differences\n",
"#### Causal Inference: The Mixtape, Chapter 9"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"R code used for replication. Slightly modified, but not substantively.\n",
"\n",
"Used the `rocker/tidyverse:latest` (on 2021/3/11) to replicate."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Package Install and Setup"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```R\n",
"install.packages(\"pkgapi\")\n",
"install.packages(\"tidyverse\")\n",
"install.packages(\"cli\")\n",
"install.packages(\"haven\")\n",
"install.packages(\"rmarkdown\")\n",
"install.packages(\"learnr\")\n",
"install.packages(\"haven\")\n",
"install.packages(\"stargazer\")\n",
"install.packages(\"estimatr\")\n",
"# This chapter only\n",
"install.packages(\"lfe\")\n",
"install.packages(\"bacondecomp\")\n",
"\n",
"library(learnr)\n",
"library(haven)\n",
"library(tidyverse)\n",
"library(stargazer)\n",
"library(estimatr)\n",
"# This chapter only\n",
"library(lfe)\n",
"library(bacondecomp)\n",
"\n",
"# Make tibbles readable in a dark theme\n",
"options(pillar.subtle = FALSE)\n",
"\n",
"read_data <- function(df) {\n",
" full_path <- paste0(\"https://raw.github.com/scunning1975/mixtape/master/\", df)\n",
" return(haven::read_dta(full_path))\n",
"}\n",
"\n",
"abortion <- read_data(\"abortion.dta\") %>%\n",
" mutate(\n",
" repeal = as_factor(repeal),\n",
" year = as_factor(year),\n",
" fip = as_factor(fip),\n",
" fa = as_factor(fa),\n",
" )\n",
"\n",
"\n",
"# Filter the data to rule out different NULL handling\n",
"\n",
"filtered.data <- abortion %>% filter(bf15 == 1) %>% drop_na(lnr)\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Apples to apples?"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```R\n",
"> filtered.data\n",
"# A tibble: 737 x 39\n",
" fip age race year sex totcase totpop rate totrate id ir crack\n",
" <fct> <dbl> <dbl> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>\n",
" 1 1 15 2 1985 2 5683 106187 6528. 5352. 14 101. 0.217\n",
" 2 1 15 2 1986 2 5344 106831 6351. 5002. 14 112. 0.277\n",
" 3 1 15 2 1987 2 4983 106496 5759. 4679 14 123. 0.211\n",
" 4 1 15 2 1988 2 5276 105238 6140. 5013. 14 134. 0.560\n",
" 5 1 15 2 1989 2 5692 102956 5952. 5529. 14 146. 0.722\n",
" 6 1 15 2 1990 2 5735 100568 6182. 5703. 14 157. 0.509\n",
" 7 1 15 2 1991 2 5587 100595 5189 5554 14 177. 0.858\n",
" 8 1 15 2 1992 2 4596 102094 4228. 4502. 14 197. 0.885\n",
" 9 1 15 2 1993 2 4217 103057 3788. 4092. 14 217. 0.254\n",
"10 1 15 2 1994 2 4173 104922 3890. 3977. 14 238. 0.123\n",
"# … with 727 more rows, and 27 more variables: alcohol <dbl>, income <dbl>,\n",
"# ur <dbl>, poverty <dbl>, black <dbl>, perc1519 <dbl>, perc <dbl>,\n",
"# repeal <fct>, aids <dbl>, aidscapita <dbl>, ac <dbl>, acc <dbl>,\n",
"# trend <dbl>, tsq <dbl>, blk <dbl>, wht <dbl>, male <dbl>, female <dbl>,\n",
"# lnr <dbl>, t <dbl>, younger <dbl>, fa <fct>, pi <dbl>, wm15 <dbl>,\n",
"# wf15 <dbl>, bm15 <dbl>, bf15 <dbl>\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Unweighted Model"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```R\n",
"# Unweighted model to rule out different weighting\n",
"easy.reg <- filtered.data %>%\n",
" lm_robust(lnr ~ repeal*year + fip + acc + ir + pi + alcohol+ crack + poverty+ income+ ur, data = .)\n",
"easy.reg\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```\n",
"1 coefficient not defined because the design matrix is rank deficient\n",
"\n",
" Estimate Std. Error t value Pr(>|t|)\n",
"(Intercept) 7.303800e+00 4.667743e-01 15.64738957 4.728204e-47\n",
"repeal1 -1.966181e+00 3.124754e-01 -6.29227569 5.760881e-10\n",
"year1986 -2.291998e-02 7.306019e-02 -0.31371361 7.538395e-01\n",
"year1987 -1.612419e-01 7.332252e-02 -2.19907763 2.822498e-02\n",
"year1988 -6.897324e-02 8.666039e-02 -0.79590268 4.263801e-01\n",
"year1989 -1.013940e-01 9.994741e-02 -1.01447384 3.107353e-01\n",
"year1990 -2.152838e-01 1.074209e-01 -2.00411365 4.547359e-02\n",
"year1991 -2.628837e-01 1.040223e-01 -2.52718712 1.173477e-02\n",
"year1992 -4.660270e-01 1.207499e-01 -3.85944019 1.250260e-04\n",
"year1993 -5.735496e-01 1.291640e-01 -4.44047401 1.055185e-05\n",
"year1994 -6.423054e-01 1.544166e-01 -4.15956128 3.618181e-05\n",
"year1995 -8.024763e-01 1.707130e-01 -4.70073370 3.168392e-06\n",
"year1996 -1.118392e+00 1.885404e-01 -5.93184165 4.875888e-09\n",
"year1997 -1.132846e+00 2.021089e-01 -5.60512476 3.080736e-08\n",
"year1998 -1.099171e+00 2.331454e-01 -4.71453061 2.967756e-06\n",
"year1999 -1.132353e+00 2.478001e-01 -4.56962383 5.851220e-06\n",
"year2000 -1.338859e+00 2.800623e-01 -4.78057403 2.164913e-06\n",
"fip2 NA NA NA NA\n",
"fip4 -5.557975e-01 1.461032e-01 -3.80414351 1.557604e-04\n",
"fip5 1.027124e-01 6.309669e-02 1.62785673 1.040413e-01\n",
"fip6 8.586731e-01 1.617385e-01 5.30902162 1.516100e-07\n",
"fip8 -6.853199e-01 1.861681e-01 -3.68118925 2.514795e-04\n",
"fip9 -1.086778e+00 2.819098e-01 -3.85505502 1.272369e-04\n",
"fip10 -4.074209e-01 2.079579e-01 -1.95915148 5.052345e-02\n",
"fip11 -1.703412e+00 4.406162e-01 -3.86597578 1.217981e-04\n",
"fip12 -8.086461e-01 1.806336e-01 -4.47671990 8.955628e-06\n",
"fip13 -4.839723e-01 1.468002e-01 -3.29680941 1.031304e-03\n",
"fip15 1.396573e-01 2.268985e-01 0.61550590 5.384368e-01\n",
"fip16 -1.528903e+00 2.521568e-01 -6.06330211 2.265139e-09\n",
"fip17 -5.184233e-01 1.759054e-01 -2.94717116 3.322231e-03\n",
"fip18 3.096018e-02 8.715098e-02 0.35524765 7.225197e-01\n",
"fip19 -3.016982e-02 9.158499e-02 -0.32941878 7.419456e-01\n",
"fip20 2.988341e-01 9.511516e-02 3.14181322 1.755574e-03\n",
"fip21 7.057992e-02 6.934636e-02 1.01778839 3.091584e-01\n",
"fip22 -6.455132e-01 1.193869e-01 -5.40690111 9.026276e-08\n",
"fip23 -2.186881e+00 1.435142e-01 -15.23808608 4.892795e-45\n",
"fip24 -1.041203e+00 1.987084e-01 -5.23985391 2.176488e-07\n",
"fip25 -1.315879e+00 2.321281e-01 -5.66876274 2.166405e-08\n",
"fip26 -6.023982e-01 1.356102e-01 -4.44213045 1.047332e-05\n",
"fip27 -1.471980e-01 1.454018e-01 -1.01235357 3.117469e-01\n",
"fip28 -8.205431e-02 8.535163e-02 -0.96136782 3.367258e-01\n",
"fip29 6.320215e-02 9.472608e-02 0.66720958 5.048758e-01\n",
"fip30 -1.361639e+00 2.043608e-01 -6.66291931 5.740793e-11\n",
"fip31 -5.806144e-03 1.236614e-01 -0.04695194 9.625660e-01\n",
"fip32 -1.602789e+00 3.644174e-01 -4.39822238 1.275664e-05\n",
"fip33 -3.519712e+00 4.081061e-01 -8.62450144 4.940975e-17\n",
"fip34 -1.781443e+00 2.500967e-01 -7.12301649 2.817756e-12\n",
"fip35 -8.645122e-01 1.464854e-01 -5.90169576 5.801428e-09\n",
"fip36 -6.600649e-02 2.209668e-01 -0.29871677 7.652519e-01\n",
"fip37 4.399084e-02 7.480454e-02 0.58807705 5.566855e-01\n",
"fip38 -2.012185e+00 2.122883e-01 -9.47854911 4.770141e-20\n",
"fip39 -1.985972e-01 1.002784e-01 -1.98045899 4.807447e-02\n",
"fip40 4.134393e-01 7.536650e-02 5.48571726 5.910175e-08\n",
"fip41 -6.700750e-01 1.660109e-01 -4.03633198 6.078117e-05\n",
"fip42 -1.478025e-01 1.296681e-01 -1.13985236 2.547690e-01\n",
"fip44 -5.482892e-01 1.684886e-01 -3.25416187 1.196408e-03\n",
"fip45 -1.016353e+00 1.326313e-01 -7.66298865 6.652714e-14\n",
"fip46 -1.365188e+00 1.994112e-01 -6.84609728 1.763767e-11\n",
"fip47 1.856564e-01 8.258824e-02 2.24797574 2.491316e-02\n",
"fip48 -4.381164e-01 1.147762e-01 -3.81713794 1.479553e-04\n",
"fip49 -5.897980e-01 2.589377e-01 -2.27776027 2.306571e-02\n",
"fip50 -1.942781e+00 2.162897e-01 -8.98231023 2.861883e-18\n",
"fip51 -6.568583e-01 1.360757e-01 -4.82715226 1.729195e-06\n",
"fip53 1.296103e+00 2.023888e-01 6.40402747 2.907860e-10\n",
"fip54 -4.509700e-01 1.495571e-01 -3.01537041 2.667053e-03\n",
"fip55 -1.693410e-01 1.701435e-01 -0.99528318 3.199701e-01\n",
"fip56 -1.868125e+00 2.152858e-01 -8.67741899 3.260374e-17\n",
"acc 1.637227e-03 7.074462e-04 2.31427708 2.096434e-02\n",
"ir -1.754211e-05 7.466289e-05 -0.23495086 8.143211e-01\n",
"pi -1.573799e-02 4.520409e-02 -0.34815403 7.278376e-01\n",
"alcohol 4.583269e-01 1.251954e-01 3.66089329 2.718228e-04\n",
"crack -6.620465e-02 2.941035e-02 -2.25106663 2.471565e-02\n",
"poverty -4.853091e-03 8.076517e-03 -0.60088909 5.481240e-01\n",
"income 5.269637e-05 1.882886e-05 2.79870198 5.283445e-03\n",
"ur 3.894939e-04 1.754596e-02 0.02219850 9.822964e-01\n",
"repeal1:year1986 -5.109105e-03 2.492869e-01 -0.02049488 9.836549e-01\n",
"repeal1:year1987 -2.788920e-02 2.869025e-01 -0.09720793 9.225913e-01\n",
"repeal1:year1988 -2.724267e-01 4.003103e-01 -0.68053886 4.964064e-01\n",
"repeal1:year1989 -2.383648e-01 3.509358e-01 -0.67922628 4.972370e-01\n",
"repeal1:year1990 -1.571047e-01 3.073729e-01 -0.51112080 6.094406e-01\n",
"repeal1:year1991 -8.983659e-02 1.963766e-01 -0.45747104 6.474859e-01\n",
"repeal1:year1992 2.352184e-02 2.240274e-01 0.10499538 9.164120e-01\n",
"repeal1:year1993 1.968860e-01 3.170595e-01 0.62097489 5.348345e-01\n",
"repeal1:year1994 2.779483e-01 3.034955e-01 0.91582356 3.601001e-01\n",
"repeal1:year1995 1.521836e-01 2.727116e-01 0.55803869 5.770107e-01\n",
"repeal1:year1996 -1.403581e-01 2.331023e-01 -0.60213077 5.472977e-01\n",
"repeal1:year1997 6.981847e-02 2.060730e-01 0.33880452 7.348668e-01\n",
"repeal1:year1998 -3.284733e-01 2.749362e-01 -1.19472523 2.326315e-01\n",
"repeal1:year1999 -2.417105e-01 3.527950e-01 -0.68513022 4.935067e-01\n",
"repeal1:year2000 -8.237028e-02 2.471748e-01 -0.33324705 7.390556e-01\n",
" CI Lower CI Upper DF\n",
"(Intercept) 6.387227e+00 8.220373e+00 648\n",
"repeal1 -2.579768e+00 -1.352595e+00 648\n",
"year1986 -1.663833e-01 1.205433e-01 648\n",
"year1987 -3.052203e-01 -1.726349e-02 648\n",
"year1988 -2.391423e-01 1.011958e-01 648\n",
"year1989 -2.976539e-01 9.486586e-02 648\n",
"year1990 -4.262189e-01 -4.348617e-03 648\n",
"year1991 -4.671451e-01 -5.862232e-02 648\n",
"year1992 -7.031353e-01 -2.289187e-01 648\n",
"year1993 -8.271802e-01 -3.199190e-01 648\n",
"year1994 -9.455228e-01 -3.390880e-01 648\n",
"year1995 -1.137694e+00 -4.672589e-01 648\n",
"year1996 -1.488615e+00 -7.481678e-01 648\n",
"year1997 -1.529713e+00 -7.359783e-01 648\n",
"year1998 -1.556983e+00 -6.413595e-01 648\n",
"year1999 -1.618941e+00 -6.457652e-01 648\n",
"year2000 -1.888798e+00 -7.889194e-01 648\n",
"fip2 NA NA NA\n",
"fip4 -8.426904e-01 -2.689047e-01 648\n",
"fip5 -2.118628e-02 2.266110e-01 648\n",
"fip6 5.410783e-01 1.176268e+00 648\n",
"fip8 -1.050885e+00 -3.197544e-01 648\n",
"fip9 -1.640345e+00 -5.332107e-01 648\n",
"fip10 -8.157736e-01 9.316821e-04 648\n",
"fip11 -2.568619e+00 -8.382036e-01 648\n",
"fip12 -1.163344e+00 -4.539482e-01 648\n",
"fip13 -7.722338e-01 -1.957108e-01 648\n",
"fip15 -3.058877e-01 5.852024e-01 648\n",
"fip16 -2.024046e+00 -1.033760e+00 648\n",
"fip17 -8.638368e-01 -1.730099e-01 648\n",
"fip18 -1.401722e-01 2.020926e-01 648\n",
"fip19 -2.100090e-01 1.496694e-01 648\n",
"fip20 1.120629e-01 4.856052e-01 648\n",
"fip21 -6.559078e-02 2.067506e-01 648\n",
"fip22 -8.799451e-01 -4.110813e-01 648\n",
"fip23 -2.468690e+00 -1.905072e+00 648\n",
"fip24 -1.431393e+00 -6.510129e-01 648\n",
"fip25 -1.771693e+00 -8.600649e-01 648\n",
"fip26 -8.686867e-01 -3.361097e-01 648\n",
"fip27 -4.327135e-01 1.383175e-01 648\n",
"fip28 -2.496535e-01 8.554485e-02 648\n",
"fip29 -1.228050e-01 2.492093e-01 648\n",
"fip30 -1.762929e+00 -9.603500e-01 648\n",
"fip31 -2.486316e-01 2.370194e-01 648\n",
"fip32 -2.318371e+00 -8.872073e-01 648\n",
"fip33 -4.321082e+00 -2.718342e+00 648\n",
"fip34 -2.272541e+00 -1.290345e+00 648\n",
"fip35 -1.152156e+00 -5.768689e-01 648\n",
"fip36 -4.999039e-01 3.678909e-01 648\n",
"fip37 -1.028977e-01 1.908794e-01 648\n",
"fip38 -2.429041e+00 -1.595329e+00 648\n",
"fip39 -3.955069e-01 -1.687420e-03 648\n",
"fip40 2.654473e-01 5.614314e-01 648\n",
"fip41 -9.960592e-01 -3.440908e-01 648\n",
"fip42 -4.024229e-01 1.068179e-01 648\n",
"fip44 -8.791388e-01 -2.174396e-01 648\n",
"fip45 -1.276792e+00 -7.559134e-01 648\n",
"fip46 -1.756759e+00 -9.736183e-01 648\n",
"fip47 2.348348e-02 3.478292e-01 648\n",
"fip48 -6.634945e-01 -2.127383e-01 648\n",
"fip49 -1.098256e+00 -8.133975e-02 648\n",
"fip50 -2.367494e+00 -1.518068e+00 648\n",
"fip51 -9.240609e-01 -3.896557e-01 648\n",
"fip53 8.986864e-01 1.693520e+00 648\n",
"fip54 -7.446450e-01 -1.572950e-01 648\n",
"fip55 -5.034401e-01 1.647582e-01 648\n",
"fip56 -2.290867e+00 -1.445383e+00 648\n",
"acc 2.480628e-04 3.026390e-03 648\n",
"ir -1.641525e-04 1.290683e-04 648\n",
"pi -1.045022e-01 7.302620e-02 648\n",
"alcohol 2.124893e-01 7.041645e-01 648\n",
"crack -1.239557e-01 -8.453563e-03 648\n",
"poverty -2.071239e-02 1.100621e-02 648\n",
"income 1.572342e-05 8.966932e-05 648\n",
"ur -3.406430e-02 3.484329e-02 648\n",
"repeal1:year1986 -4.946168e-01 4.843986e-01 648\n",
"repeal1:year1987 -5.912600e-01 5.354816e-01 648\n",
"repeal1:year1988 -1.058489e+00 5.136352e-01 648\n",
"repeal1:year1989 -9.274734e-01 4.507438e-01 648\n",
"repeal1:year1990 -7.606719e-01 4.464625e-01 648\n",
"repeal1:year1991 -4.754478e-01 2.957746e-01 648\n",
"repeal1:year1992 -4.163855e-01 4.634291e-01 648\n",
"repeal1:year1993 -4.257020e-01 8.194740e-01 648\n",
"repeal1:year1994 -3.180050e-01 8.739017e-01 648\n",
"repeal1:year1995 -3.833214e-01 6.876886e-01 648\n",
"repeal1:year1996 -5.980851e-01 3.173690e-01 648\n",
"repeal1:year1997 -3.348331e-01 4.744700e-01 648\n",
"repeal1:year1998 -8.683468e-01 2.114002e-01 648\n",
"repeal1:year1999 -9.344700e-01 4.510489e-01 648\n",
"repeal1:year2000 -5.677306e-01 4.029900e-01 648\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Weighted Model"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```R\n",
"\n",
"# Weighted model\n",
"reg <- filtered.data %>%\n",
" lm_robust(lnr ~ repeal*year + fip + acc + ir + pi + alcohol+ crack + poverty+ income+ ur,\n",
" data = ., weights = totpop, clusters = fip)\n",
"reg\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```\n",
"> reg\n",
"1 coefficient not defined because the design matrix is rank deficient\n",
"\n",
" Estimate Std. Error t value Pr(>|t|)\n",
"(Intercept) 7.475343e+00 1.658847e+00 4.50634752 1.369649e-03\n",
"repeal1 -1.749250e+00 6.401727e-01 -2.73246545 5.227752e-02\n",
"year1986 -2.082737e-02 6.663089e-02 -0.31257823 7.601400e-01\n",
"year1987 -2.864665e-01 1.278328e-01 -2.24094668 4.644352e-02\n",
"year1988 -3.494648e-01 1.749727e-01 -1.99725391 7.145663e-02\n",
"year1989 -4.072200e-01 2.413305e-01 -1.68739534 1.206377e-01\n",
"year1990 -4.864117e-01 3.182625e-01 -1.52833529 1.569205e-01\n",
"year1991 -5.353034e-01 3.258975e-01 -1.64255114 1.296883e-01\n",
"year1992 -7.866236e-01 3.994566e-01 -1.96923425 7.870726e-02\n",
"year1993 -9.767166e-01 4.270296e-01 -2.28723402 4.605464e-02\n",
"year1994 -1.064367e+00 4.776021e-01 -2.22856369 5.183896e-02\n",
"year1995 -1.356245e+00 5.388305e-01 -2.51701635 3.155354e-02\n",
"year1996 -1.520062e+00 6.016796e-01 -2.52636483 3.183765e-02\n",
"year1997 -1.571846e+00 6.644358e-01 -2.36568539 4.103017e-02\n",
"year1998 -1.545129e+00 7.509867e-01 -2.05746564 6.887182e-02\n",
"year1999 -1.587414e+00 8.064132e-01 -1.96848660 7.994068e-02\n",
"year2000 -1.672572e+00 9.065627e-01 -1.84495968 9.624871e-02\n",
"fip2 -8.647049e-02 4.026150e-01 -0.21477215 8.350741e-01\n",
"fip4 -9.589386e-01 7.448868e-01 -1.28736141 2.500427e-01\n",
"fip5 8.603103e-02 8.967369e-02 0.95937871 3.665139e-01\n",
"fip6 4.878394e-01 3.649714e-01 1.33665097 2.295107e-01\n",
"fip8 -8.246050e-01 6.117645e-01 -1.34791248 2.189549e-01\n",
"fip9 -1.680892e+00 1.107248e+00 -1.51808057 1.868691e-01\n",
"fip10 -7.457544e-01 7.373248e-01 -1.01143264 3.556510e-01\n",
"fip11 -2.112355e+00 1.901975e+00 -1.11061100 3.145083e-01\n",
"fip12 -1.094748e+00 7.681499e-01 -1.42517446 2.122327e-01\n",
"fip13 -7.215762e-01 3.822429e-01 -1.88774251 1.297186e-01\n",
"fip15 NA NA NA NA\n",
"fip16 -1.793719e+00 2.492800e-01 -7.19560060 1.219555e-05\n",
"fip17 -7.228205e-01 6.786144e-01 -1.06514171 3.273357e-01\n",
"fip18 -1.431537e-01 2.653865e-01 -0.53941586 6.039826e-01\n",
"fip19 -1.497970e-01 2.556753e-01 -0.58588739 5.699621e-01\n",
"fip20 8.004447e-02 2.914899e-01 0.27460466 7.891561e-01\n",
"fip21 -1.103709e-01 6.089921e-02 -1.81235286 9.582570e-02\n",
"fip22 -7.439527e-01 3.869440e-01 -1.92263653 8.984780e-02\n",
"fip23 -2.021225e+00 1.625117e-01 -12.43741167 1.306890e-07\n",
"fip24 -1.474562e+00 7.678151e-01 -1.92046438 1.158568e-01\n",
"fip25 -1.894324e+00 8.342228e-01 -2.27076525 7.845338e-02\n",
"fip26 -7.679939e-01 5.243854e-01 -1.46455997 1.965299e-01\n",
"fip27 -3.203199e-01 5.852182e-01 -0.54735115 6.027069e-01\n",
"fip28 -2.758883e-02 1.573867e-01 -0.17529327 8.638283e-01\n",
"fip29 -1.669216e-01 4.265189e-01 -0.39135812 7.125428e-01\n",
"fip30 -1.266169e+00 2.174461e-01 -5.82290806 7.359898e-05\n",
"fip31 -1.175657e-01 4.240238e-01 -0.27726195 7.881407e-01\n",
"fip32 -1.964808e+00 1.516241e+00 -1.29584203 2.392429e-01\n",
"fip33 -3.590973e+00 1.250231e+00 -2.87224786 2.053894e-02\n",
"fip34 -2.275554e+00 9.861880e-01 -2.30742412 7.003805e-02\n",
"fip35 -1.181403e+00 4.553529e-01 -2.59447868 3.046693e-02\n",
"fip36 -5.981228e-01 2.868533e-01 -2.08511687 6.521830e-02\n",
"fip37 -8.849572e-02 1.815960e-01 -0.48732211 6.372454e-01\n",
"fip38 -1.826491e+00 1.494071e-01 -12.22493041 4.554158e-08\n",
"fip39 -4.892242e-01 3.695991e-01 -1.32366174 2.401485e-01\n",
"fip40 3.048617e-01 7.968022e-02 3.82606433 2.685852e-03\n",
"fip41 -1.122812e+00 5.167108e-01 -2.17299932 7.686070e-02\n",
"fip42 -6.310678e-01 4.975957e-01 -1.26823385 2.539386e-01\n",
"fip44 -1.209623e+00 7.927586e-01 -1.52584082 2.028076e-01\n",
"fip45 -1.110900e+00 2.731853e-01 -4.06647077 4.669939e-03\n",
"fip46 -2.028952e+00 9.538540e-01 -2.12710941 7.665033e-02\n",
"fip47 5.924823e-02 1.439817e-01 0.41149819 6.881190e-01\n",
"fip48 -6.054829e-01 4.816200e-01 -1.25717959 2.583612e-01\n",
"fip49 -1.266724e+00 2.250145e-01 -5.62951974 6.993460e-05\n",
"fip50 -1.997358e+00 2.975760e-01 -6.71209329 1.029781e-04\n",
"fip51 -9.493606e-01 4.788519e-01 -1.98257677 8.834780e-02\n",
"fip53 7.854260e-01 1.193181e-01 6.58262411 3.269844e-05\n",
"fip54 -5.213726e-01 1.350920e-01 -3.85938843 4.060009e-03\n",
"fip55 -4.836776e-01 8.021240e-01 -0.60299604 5.722683e-01\n",
"fip56 -2.192106e+00 8.127714e-01 -2.69707612 3.897138e-02\n",
"acc 2.877740e-03 1.782471e-03 1.61446697 1.621803e-01\n",
"ir 4.438901e-04 5.582632e-04 0.79512697 4.557450e-01\n",
"pi -3.902207e-02 7.056469e-02 -0.55299710 5.908399e-01\n",
"alcohol 4.468421e-01 4.484692e-01 0.99637190 3.416850e-01\n",
"crack 5.282875e-02 3.701847e-02 1.42709140 2.088972e-01\n",
"poverty -2.760983e-03 1.496253e-02 -0.18452648 8.573763e-01\n",
"income 5.847156e-05 6.126983e-05 0.95432877 3.654394e-01\n",
"ur -2.782371e-02 3.877280e-02 -0.71760897 4.928258e-01\n",
"repeal1:year1986 -2.590538e-01 8.090168e-02 -3.20208207 1.487582e-01\n",
"repeal1:year1987 -3.417380e-01 2.317637e-01 -1.47451048 3.384673e-01\n",
"repeal1:year1988 -6.105532e-01 2.908114e-01 -2.09948158 2.336005e-01\n",
"repeal1:year1989 -7.850003e-01 3.244798e-01 -2.41925765 1.954591e-01\n",
"repeal1:year1990 -6.316315e-01 2.111608e-01 -2.99123431 1.471340e-01\n",
"repeal1:year1991 -5.526026e-01 1.802310e-01 -3.06608002 1.408578e-01\n",
"repeal1:year1992 -4.424166e-01 1.942230e-01 -2.27787968 1.998834e-01\n",
"repeal1:year1993 -3.060626e-01 2.702939e-01 -1.13233268 4.094116e-01\n",
"repeal1:year1994 -1.178897e-01 3.434297e-01 -0.34327180 7.736276e-01\n",
"repeal1:year1995 2.102913e-02 3.656786e-01 0.05750715 9.610165e-01\n",
"repeal1:year1996 -1.235341e-01 3.907048e-01 -0.31618282 7.915766e-01\n",
"repeal1:year1997 2.085286e-02 4.940011e-01 0.04221217 9.714596e-01\n",
"repeal1:year1998 -3.605240e-02 6.248853e-01 -0.05769442 9.609411e-01\n",
"repeal1:year1999 1.470912e-02 6.353063e-01 0.02315280 9.842647e-01\n",
"repeal1:year2000 4.145081e-02 6.591620e-01 0.06288410 9.572352e-01\n",
" CI Lower CI Upper DF\n",
"(Intercept) 3.739353e+00 11.2113328136 9.269811\n",
"repeal1 -3.526256e+00 0.0277568459 4.002271\n",
"year1986 -1.665336e-01 0.1248789018 11.618802\n",
"year1987 -5.675546e-01 -0.0053783449 11.087267\n",
"year1988 -7.351575e-01 0.0362278224 10.865956\n",
"year1989 -9.407117e-01 0.1262716370 10.620210\n",
"year1990 -1.193959e+00 0.2211356482 10.168063\n",
"year1991 -1.255635e+00 0.1850283716 10.632544\n",
"year1992 -1.682847e+00 0.1095998262 9.515690\n",
"year1993 -1.932526e+00 -0.0209073505 9.676598\n",
"year1994 -2.139160e+00 0.0104272276 9.319334\n",
"year1995 -2.564579e+00 -0.1479113085 9.548225\n",
"year1996 -2.876043e+00 -0.1640811935 9.228173\n",
"year1997 -3.064791e+00 -0.0789008395 9.418224\n",
"year1998 -3.236540e+00 0.1462817152 9.267253\n",
"year1999 -3.406471e+00 0.2316442438 9.171127\n",
"year2000 -3.705677e+00 0.3605334179 9.544036\n",
"fip2 -1.007452e+00 0.8345109768 8.389543\n",
"fip4 -2.828312e+00 0.9104344031 5.437454\n",
"fip5 -1.221610e-01 0.2942230414 7.701362\n",
"fip6 -4.039180e-01 1.3795967862 6.036167\n",
"fip8 -2.265978e+00 0.6167685390 7.127258\n",
"fip9 -4.489188e+00 1.1274045164 5.234812\n",
"fip10 -2.608255e+00 1.1167466909 5.310508\n",
"fip11 -6.917796e+00 2.6930870276 5.306547\n",
"fip12 -3.056798e+00 0.8673023704 5.108451\n",
"fip13 -1.769084e+00 0.3259318340 4.137528\n",
"fip15 NA NA NA\n",
"fip16 -2.338059e+00 -1.2493797360 11.765097\n",
"fip17 -2.378700e+00 0.9330588456 6.069991\n",
"fip18 -7.529033e-01 0.4665958984 8.171839\n",
"fip19 -7.136306e-01 0.4140366989 10.827461\n",
"fip20 -5.685800e-01 0.7286689580 10.098234\n",
"fip21 -2.435392e-01 0.0227974816 11.621898\n",
"fip22 -1.632441e+00 0.1445351786 8.201777\n",
"fip23 -2.381069e+00 -1.6613816247 10.484407\n",
"fip24 -3.479520e+00 0.5303971176 4.753193\n",
"fip25 -4.114918e+00 0.3262694423 4.485397\n",
"fip26 -2.071779e+00 0.5357916718 5.630993\n",
"fip27 -1.731189e+00 1.0905492203 6.393770\n",
"fip28 -3.712450e-01 0.3160673336 11.771428\n",
"fip29 -1.282549e+00 0.9487052861 4.728274\n",
"fip30 -1.738575e+00 -0.7937624600 12.321694\n",
"fip31 -1.083928e+00 0.8487965241 8.582852\n",
"fip32 -5.608744e+00 1.6791268530 6.484510\n",
"fip33 -6.468934e+00 -0.7130126541 8.081894\n",
"fip34 -4.823955e+00 0.2728472578 4.914493\n",
"fip35 -2.221598e+00 -0.1412092259 8.459797\n",
"fip36 -1.242248e+00 0.0460022436 9.460424\n",
"fip37 -4.969249e-01 0.3199334401 9.355841\n",
"fip38 -2.152534e+00 -1.5004480557 11.832202\n",
"fip39 -1.424956e+00 0.4465074463 5.267956\n",
"fip40 1.300365e-01 0.4796868581 11.290954\n",
"fip41 -2.415431e+00 0.1698059702 5.501889\n",
"fip42 -1.863600e+00 0.6014647486 5.713400\n",
"fip44 -3.423929e+00 1.0046817783 3.940123\n",
"fip45 -1.755568e+00 -0.4662323019 7.071017\n",
"fip46 -4.352295e+00 0.2943914690 6.115606\n",
"fip47 -2.552283e-01 0.3737247851 11.741645\n",
"fip48 -1.803251e+00 0.5922853282 5.625406\n",
"fip49 -1.750785e+00 -0.7826622609 13.565768\n",
"fip50 -2.674212e+00 -1.3205038457 8.689024\n",
"fip51 -2.084470e+00 0.1857485202 6.915689\n",
"fip53 5.240517e-01 1.0468003397 11.445951\n",
"fip54 -8.283047e-01 -0.2144405329 8.750633\n",
"fip55 -2.532199e+00 1.5648433398 5.111069\n",
"fip56 -4.225740e+00 -0.1584722842 5.498038\n",
"acc -1.585233e-03 0.0073407122 5.484059\n",
"ir -9.098678e-04 0.0017976480 6.232597\n",
"pi -1.934910e-01 0.1154468277 11.514779\n",
"alcohol -5.472135e-01 1.4408977472 10.400887\n",
"crack -4.033155e-02 0.1459890456 5.382560\n",
"poverty -3.620879e-02 0.0306868202 9.764713\n",
"income -8.065125e-05 0.0001975944 8.784182\n",
"ur -1.167730e-01 0.0611255680 8.245333\n",
"repeal1:year1986 -8.917936e-01 0.3736859319 1.268677\n",
"repeal1:year1987 -2.115507e+00 1.4320311251 1.284685\n",
"repeal1:year1988 -2.765610e+00 1.5445032932 1.309373\n",
"repeal1:year1989 -3.099724e+00 1.5297234582 1.339972\n",
"repeal1:year1990 -2.070152e+00 0.8068887875 1.379284\n",
"repeal1:year1991 -1.759082e+00 0.6538768744 1.394950\n",
"repeal1:year1992 -1.707836e+00 0.8230030707 1.419972\n",
"repeal1:year1993 -1.988210e+00 1.3760845940 1.464797\n",
"repeal1:year1994 -2.207596e+00 1.9718163483 1.488047\n",
"repeal1:year1995 -2.274876e+00 2.3169340772 1.455928\n",
"repeal1:year1996 -2.614552e+00 2.3674842206 1.440751\n",
"repeal1:year1997 -3.176038e+00 3.2177436826 1.426379\n",
"repeal1:year1998 -4.010702e+00 3.9385967330 1.443069\n",
"repeal1:year1999 -3.904647e+00 3.9340649263 1.473722\n",
"repeal1:year2000 -3.962451e+00 4.0453522150 1.489881\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Python Setup"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"'0.12.2'"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"import numpy as np\n",
"import pandas as pd\n",
"import plotnine as p\n",
"import statsmodels.api as sm\n",
"import statsmodels.formula.api as smf\n",
"\n",
"sm.__version__"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Undo some personal preference options in my environment."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"pd.options.display.show_dimensions = True\n",
"pd.options.display.precision = 8\n",
"pd.options.display.float_format = lambda x: f\"{x:.8f}\"\n",
"np.set_printoptions(precision=8)"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"abortion_data = (\n",
" pd.read_stata(\"https://raw.github.com/scunning1975/mixtape/master/abortion.dta\")\n",
")"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"abortion_bf15 = (\n",
" abortion_data\n",
" .pipe(lambda dta: dta.loc[dta.bf15 == 1])\n",
" .dropna(subset=[\"lnr\"])\n",
")"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>fip</th>\n",
" <th>age</th>\n",
" <th>race</th>\n",
" <th>year</th>\n",
" <th>sex</th>\n",
" <th>totcase</th>\n",
" <th>totpop</th>\n",
" <th>rate</th>\n",
" <th>totrate</th>\n",
" <th>id</th>\n",
" <th>...</th>\n",
" <th>female</th>\n",
" <th>lnr</th>\n",
" <th>t</th>\n",
" <th>younger</th>\n",
" <th>fa</th>\n",
" <th>pi</th>\n",
" <th>wm15</th>\n",
" <th>wf15</th>\n",
" <th>bm15</th>\n",
" <th>bf15</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>19</th>\n",
" <td>1.00000000</td>\n",
" <td>15.00000000</td>\n",
" <td>2.00000000</td>\n",
" <td>1985.00000000</td>\n",
" <td>2</td>\n",
" <td>5683.00000000</td>\n",
" <td>106187</td>\n",
" <td>6527.50000000</td>\n",
" <td>5351.89990234</td>\n",
" <td>14.00000000</td>\n",
" <td>...</td>\n",
" <td>1.00000000</td>\n",
" <td>8.78377914</td>\n",
" <td>1.00000000</td>\n",
" <td>1.00000000</td>\n",
" <td>1.00000000</td>\n",
" <td>0.00000000</td>\n",
" <td>0.00000000</td>\n",
" <td>0.00000000</td>\n",
" <td>0.00000000</td>\n",
" <td>1.00000000</td>\n",
" </tr>\n",
" <tr>\n",
" <th>39</th>\n",
" <td>1.00000000</td>\n",
" <td>15.00000000</td>\n",
" <td>2.00000000</td>\n",
" <td>1986.00000000</td>\n",
" <td>2</td>\n",
" <td>5344.00000000</td>\n",
" <td>106831</td>\n",
" <td>6351.20019531</td>\n",
" <td>5002.29980469</td>\n",
" <td>14.00000000</td>\n",
" <td>...</td>\n",
" <td>1.00000000</td>\n",
" <td>8.75639915</td>\n",
" <td>2.00000000</td>\n",
" <td>1.00000000</td>\n",
" <td>1.00000000</td>\n",
" <td>0.00000000</td>\n",
" <td>0.00000000</td>\n",
" <td>0.00000000</td>\n",
" <td>0.00000000</td>\n",
" <td>1.00000000</td>\n",
" </tr>\n",
" <tr>\n",
" <th>71</th>\n",
" <td>1.00000000</td>\n",
" <td>15.00000000</td>\n",
" <td>2.00000000</td>\n",
" <td>1987.00000000</td>\n",
" <td>2</td>\n",
" <td>4983.00000000</td>\n",
" <td>106496</td>\n",
" <td>5759.10009766</td>\n",
" <td>4679.00000000</td>\n",
" <td>14.00000000</td>\n",
" <td>...</td>\n",
" <td>1.00000000</td>\n",
" <td>8.65853691</td>\n",
" <td>3.00000000</td>\n",
" <td>1.00000000</td>\n",
" <td>1.00000000</td>\n",
" <td>1.00000000</td>\n",
" <td>0.00000000</td>\n",
" <td>0.00000000</td>\n",
" <td>0.00000000</td>\n",
" <td>1.00000000</td>\n",
" </tr>\n",
" <tr>\n",
" <th>89</th>\n",
" <td>1.00000000</td>\n",
" <td>15.00000000</td>\n",
" <td>2.00000000</td>\n",
" <td>1988.00000000</td>\n",
" <td>2</td>\n",
" <td>5276.00000000</td>\n",
" <td>105238</td>\n",
" <td>6139.60009766</td>\n",
" <td>5013.39990234</td>\n",
" <td>14.00000000</td>\n",
" <td>...</td>\n",
" <td>1.00000000</td>\n",
" <td>8.72251511</td>\n",
" <td>4.00000000</td>\n",
" <td>1.00000000</td>\n",
" <td>1.00000000</td>\n",
" <td>1.00000000</td>\n",
" <td>0.00000000</td>\n",
" <td>0.00000000</td>\n",
" <td>0.00000000</td>\n",
" <td>1.00000000</td>\n",
" </tr>\n",
" <tr>\n",
" <th>106</th>\n",
" <td>1.00000000</td>\n",
" <td>15.00000000</td>\n",
" <td>2.00000000</td>\n",
" <td>1989.00000000</td>\n",
" <td>2</td>\n",
" <td>5692.00000000</td>\n",
" <td>102956</td>\n",
" <td>5951.50000000</td>\n",
" <td>5528.60009766</td>\n",
" <td>14.00000000</td>\n",
" <td>...</td>\n",
" <td>1.00000000</td>\n",
" <td>8.69139862</td>\n",
" <td>5.00000000</td>\n",
" <td>1.00000000</td>\n",
" <td>1.00000000</td>\n",
" <td>1.00000000</td>\n",
" <td>0.00000000</td>\n",
" <td>0.00000000</td>\n",
" <td>0.00000000</td>\n",
" <td>1.00000000</td>\n",
" </tr>\n",
" <tr>\n",
" <th>...</th>\n",
" <td>...</td>\n",
" <td>...</td>\n",
" <td>...</td>\n",
" <td>...</td>\n",
" <td>...</td>\n",
" <td>...</td>\n",
" <td>...</td>\n",
" <td>...</td>\n",
" <td>...</td>\n",
" <td>...</td>\n",
" <td>...</td>\n",
" <td>...</td>\n",
" <td>...</td>\n",
" <td>...</td>\n",
" <td>...</td>\n",
" <td>...</td>\n",
" <td>...</td>\n",
" <td>...</td>\n",
" <td>...</td>\n",
" <td>...</td>\n",
" <td>...</td>\n",
" </tr>\n",
" <tr>\n",
" <th>19380</th>\n",
" <td>56.00000000</td>\n",
" <td>15.00000000</td>\n",
" <td>2.00000000</td>\n",
" <td>1992.00000000</td>\n",
" <td>2</td>\n",
" <td>5.00000000</td>\n",
" <td>375</td>\n",
" <td>1242.19995117</td>\n",
" <td>1333.30004883</td>\n",
" <td>1214.00000000</td>\n",
" <td>...</td>\n",
" <td>1.00000000</td>\n",
" <td>7.12463903</td>\n",
" <td>8.00000000</td>\n",
" <td>1.00000000</td>\n",
" <td>101.00000000</td>\n",
" <td>1.00000000</td>\n",
" <td>0.00000000</td>\n",
" <td>0.00000000</td>\n",
" <td>0.00000000</td>\n",
" <td>1.00000000</td>\n",
" </tr>\n",
" <tr>\n",
" <th>19401</th>\n",
" <td>56.00000000</td>\n",
" <td>15.00000000</td>\n",
" <td>2.00000000</td>\n",
" <td>1993.00000000</td>\n",
" <td>2</td>\n",
" <td>6.00000000</td>\n",
" <td>339</td>\n",
" <td>2069.00000000</td>\n",
" <td>1769.90002441</td>\n",
" <td>1214.00000000</td>\n",
" <td>...</td>\n",
" <td>1.00000000</td>\n",
" <td>7.63482046</td>\n",
" <td>9.00000000</td>\n",
" <td>1.00000000</td>\n",
" <td>101.00000000</td>\n",
" <td>1.00000000</td>\n",
" <td>0.00000000</td>\n",
" <td>0.00000000</td>\n",
" <td>0.00000000</td>\n",
" <td>1.00000000</td>\n",
" </tr>\n",
" <tr>\n",
" <th>19451</th>\n",
" <td>56.00000000</td>\n",
" <td>15.00000000</td>\n",
" <td>2.00000000</td>\n",
" <td>1995.00000000</td>\n",
" <td>2</td>\n",
" <td>2.00000000</td>\n",
" <td>379</td>\n",
" <td>598.79998779</td>\n",
" <td>527.70001221</td>\n",
" <td>1214.00000000</td>\n",
" <td>...</td>\n",
" <td>1.00000000</td>\n",
" <td>6.39492750</td>\n",
" <td>11.00000000</td>\n",
" <td>1.00000000</td>\n",
" <td>101.00000000</td>\n",
" <td>1.00000000</td>\n",
" <td>0.00000000</td>\n",
" <td>0.00000000</td>\n",
" <td>0.00000000</td>\n",
" <td>1.00000000</td>\n",
" </tr>\n",
" <tr>\n",
" <th>19477</th>\n",
" <td>56.00000000</td>\n",
" <td>15.00000000</td>\n",
" <td>2.00000000</td>\n",
" <td>1996.00000000</td>\n",
" <td>2</td>\n",
" <td>2.00000000</td>\n",
" <td>388</td>\n",
" <td>584.79998779</td>\n",
" <td>515.50000000</td>\n",
" <td>1214.00000000</td>\n",
" <td>...</td>\n",
" <td>1.00000000</td>\n",
" <td>6.37126970</td>\n",
" <td>12.00000000</td>\n",
" <td>1.00000000</td>\n",
" <td>101.00000000</td>\n",
" <td>1.00000000</td>\n",
" <td>0.00000000</td>\n",
" <td>0.00000000</td>\n",
" <td>0.00000000</td>\n",
" <td>1.00000000</td>\n",
" </tr>\n",
" <tr>\n",
" <th>19576</th>\n",
" <td>56.00000000</td>\n",
" <td>15.00000000</td>\n",
" <td>2.00000000</td>\n",
" <td>2000.00000000</td>\n",
" <td>2</td>\n",
" <td>2.00000000</td>\n",
" <td>420</td>\n",
" <td>558.70001221</td>\n",
" <td>476.20001221</td>\n",
" <td>1214.00000000</td>\n",
" <td>...</td>\n",
" <td>1.00000000</td>\n",
" <td>6.32561255</td>\n",
" <td>16.00000000</td>\n",
" <td>1.00000000</td>\n",
" <td>101.00000000</td>\n",
" <td>1.00000000</td>\n",
" <td>0.00000000</td>\n",
" <td>0.00000000</td>\n",
" <td>0.00000000</td>\n",
" <td>1.00000000</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"<p>737 rows × 39 columns</p>\n",
"</div>"
],
"text/plain": [
" fip age race year sex totcase \\\n",
"19 1.00000000 15.00000000 2.00000000 1985.00000000 2 5683.00000000 \n",
"39 1.00000000 15.00000000 2.00000000 1986.00000000 2 5344.00000000 \n",
"71 1.00000000 15.00000000 2.00000000 1987.00000000 2 4983.00000000 \n",
"89 1.00000000 15.00000000 2.00000000 1988.00000000 2 5276.00000000 \n",
"106 1.00000000 15.00000000 2.00000000 1989.00000000 2 5692.00000000 \n",
"... ... ... ... ... ... ... \n",
"19380 56.00000000 15.00000000 2.00000000 1992.00000000 2 5.00000000 \n",
"19401 56.00000000 15.00000000 2.00000000 1993.00000000 2 6.00000000 \n",
"19451 56.00000000 15.00000000 2.00000000 1995.00000000 2 2.00000000 \n",
"19477 56.00000000 15.00000000 2.00000000 1996.00000000 2 2.00000000 \n",
"19576 56.00000000 15.00000000 2.00000000 2000.00000000 2 2.00000000 \n",
"\n",
" totpop rate totrate id ... female \\\n",
"19 106187 6527.50000000 5351.89990234 14.00000000 ... 1.00000000 \n",
"39 106831 6351.20019531 5002.29980469 14.00000000 ... 1.00000000 \n",
"71 106496 5759.10009766 4679.00000000 14.00000000 ... 1.00000000 \n",
"89 105238 6139.60009766 5013.39990234 14.00000000 ... 1.00000000 \n",
"106 102956 5951.50000000 5528.60009766 14.00000000 ... 1.00000000 \n",
"... ... ... ... ... ... ... \n",
"19380 375 1242.19995117 1333.30004883 1214.00000000 ... 1.00000000 \n",
"19401 339 2069.00000000 1769.90002441 1214.00000000 ... 1.00000000 \n",
"19451 379 598.79998779 527.70001221 1214.00000000 ... 1.00000000 \n",
"19477 388 584.79998779 515.50000000 1214.00000000 ... 1.00000000 \n",
"19576 420 558.70001221 476.20001221 1214.00000000 ... 1.00000000 \n",
"\n",
" lnr t younger fa pi wm15 \\\n",
"19 8.78377914 1.00000000 1.00000000 1.00000000 0.00000000 0.00000000 \n",
"39 8.75639915 2.00000000 1.00000000 1.00000000 0.00000000 0.00000000 \n",
"71 8.65853691 3.00000000 1.00000000 1.00000000 1.00000000 0.00000000 \n",
"89 8.72251511 4.00000000 1.00000000 1.00000000 1.00000000 0.00000000 \n",
"106 8.69139862 5.00000000 1.00000000 1.00000000 1.00000000 0.00000000 \n",
"... ... ... ... ... ... ... \n",
"19380 7.12463903 8.00000000 1.00000000 101.00000000 1.00000000 0.00000000 \n",
"19401 7.63482046 9.00000000 1.00000000 101.00000000 1.00000000 0.00000000 \n",
"19451 6.39492750 11.00000000 1.00000000 101.00000000 1.00000000 0.00000000 \n",
"19477 6.37126970 12.00000000 1.00000000 101.00000000 1.00000000 0.00000000 \n",
"19576 6.32561255 16.00000000 1.00000000 101.00000000 1.00000000 0.00000000 \n",
"\n",
" wf15 bm15 bf15 \n",
"19 0.00000000 0.00000000 1.00000000 \n",
"39 0.00000000 0.00000000 1.00000000 \n",
"71 0.00000000 0.00000000 1.00000000 \n",
"89 0.00000000 0.00000000 1.00000000 \n",
"106 0.00000000 0.00000000 1.00000000 \n",
"... ... ... ... \n",
"19380 0.00000000 0.00000000 1.00000000 \n",
"19401 0.00000000 0.00000000 1.00000000 \n",
"19451 0.00000000 0.00000000 1.00000000 \n",
"19477 0.00000000 0.00000000 1.00000000 \n",
"19576 0.00000000 0.00000000 1.00000000 \n",
"\n",
"[737 rows x 39 columns]"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"abortion_bf15"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Unweighted Model"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This _mostly_ matches the unweighted model from `lm_robust`. Notice the intercept, the year coefficients, etc. The estimates differ, however, for the parameters that aren't well identified because of the multicollinearity issues. The default, using SVD (via the Moore-Penrose psuedo-inverse) more or less spreads the effect across collinear columns."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Intercept 7.30380009\n",
"C(repeal)[T.1.0] -1.26707970\n",
"C(year)[T.1986.0] -0.02291998\n",
"C(year)[T.1987.0] -0.16124192\n",
"C(year)[T.1988.0] -0.06897324\n",
"C(year)[T.1989.0] -0.10139403\n",
"C(year)[T.1990.0] -0.21528376\n",
"C(year)[T.1991.0] -0.26288372\n",
"C(year)[T.1992.0] -0.46602701\n",
"C(year)[T.1993.0] -0.57354959\n",
"C(year)[T.1994.0] -0.64230542\n",
"C(year)[T.1995.0] -0.80247626\n",
"C(year)[T.1996.0] -1.11839159\n",
"C(year)[T.1997.0] -1.13284577\n",
"C(year)[T.1998.0] -1.09917115\n",
"C(year)[T.1999.0] -1.13235332\n",
"C(year)[T.2000.0] -1.33885860\n",
"C(fip)[T.2.0] -0.69910141\n",
"C(fip)[T.4.0] -0.55579753\n",
"C(fip)[T.5.0] 0.10271236\n",
"C(fip)[T.6.0] 0.15957167\n",
"C(fip)[T.8.0] -0.68531986\n",
"C(fip)[T.9.0] -1.08677769\n",
"C(fip)[T.10.0] -0.40742094\n",
"C(fip)[T.11.0] -1.70341156\n",
"C(fip)[T.12.0] -0.80864615\n",
"C(fip)[T.13.0] -0.48397227\n",
"C(fip)[T.15.0] -0.55944406\n",
"C(fip)[T.16.0] -1.52890281\n",
"C(fip)[T.17.0] -0.51842334\n",
"C(fip)[T.18.0] 0.03096018\n",
"C(fip)[T.19.0] -0.03016982\n",
"C(fip)[T.20.0] 0.29883407\n",
"C(fip)[T.21.0] 0.07057992\n",
"C(fip)[T.22.0] -0.64551321\n",
"C(fip)[T.23.0] -2.18688140\n",
"C(fip)[T.24.0] -1.04120293\n",
"C(fip)[T.25.0] -1.31587892\n",
"C(fip)[T.26.0] -0.60239822\n",
"C(fip)[T.27.0] -0.14719801\n",
"C(fip)[T.28.0] -0.08205431\n",
"C(fip)[T.29.0] 0.06320215\n",
"C(fip)[T.30.0] -1.36163932\n",
"C(fip)[T.31.0] -0.00580614\n",
"C(fip)[T.32.0] -1.60278891\n",
"C(fip)[T.33.0] -3.51971156\n",
"C(fip)[T.34.0] -1.78144294\n",
"C(fip)[T.35.0] -0.86451223\n",
"C(fip)[T.36.0] -0.76510790\n",
"C(fip)[T.37.0] 0.04399084\n",
"C(fip)[T.38.0] -2.01218528\n",
"C(fip)[T.39.0] -0.19859718\n",
"C(fip)[T.40.0] 0.41343931\n",
"C(fip)[T.41.0] -0.67007500\n",
"C(fip)[T.42.0] -0.14780248\n",
"C(fip)[T.44.0] -0.54828921\n",
"C(fip)[T.45.0] -1.01635253\n",
"C(fip)[T.46.0] -1.36518847\n",
"C(fip)[T.47.0] 0.18565636\n",
"C(fip)[T.48.0] -0.43811643\n",
"C(fip)[T.49.0] -0.58979804\n",
"C(fip)[T.50.0] -1.94278109\n",
"C(fip)[T.51.0] -0.65685829\n",
"C(fip)[T.53.0] 0.59700200\n",
"C(fip)[T.54.0] -0.45096998\n",
"C(fip)[T.55.0] -0.16934096\n",
"C(fip)[T.56.0] -1.86812527\n",
"C(repeal)[T.1.0]:C(year)[T.1986.0] -0.00510910\n",
"C(repeal)[T.1.0]:C(year)[T.1987.0] -0.02788920\n",
"C(repeal)[T.1.0]:C(year)[T.1988.0] -0.27242669\n",
"C(repeal)[T.1.0]:C(year)[T.1989.0] -0.23836481\n",
"C(repeal)[T.1.0]:C(year)[T.1990.0] -0.15710470\n",
"C(repeal)[T.1.0]:C(year)[T.1991.0] -0.08983659\n",
"C(repeal)[T.1.0]:C(year)[T.1992.0] 0.02352184\n",
"C(repeal)[T.1.0]:C(year)[T.1993.0] 0.19688597\n",
"C(repeal)[T.1.0]:C(year)[T.1994.0] 0.27794833\n",
"C(repeal)[T.1.0]:C(year)[T.1995.0] 0.15218360\n",
"C(repeal)[T.1.0]:C(year)[T.1996.0] -0.14035806\n",
"C(repeal)[T.1.0]:C(year)[T.1997.0] 0.06981847\n",
"C(repeal)[T.1.0]:C(year)[T.1998.0] -0.32847326\n",
"C(repeal)[T.1.0]:C(year)[T.1999.0] -0.24171052\n",
"C(repeal)[T.1.0]:C(year)[T.2000.0] -0.08237028\n",
"acc 0.00163723\n",
"ir -0.00001754\n",
"pi -0.01573799\n",
"alcohol 0.45832689\n",
"crack -0.06620465\n",
"poverty -0.00485309\n",
"income 0.00005270\n",
"ur 0.00038949\n",
"Length: 90, dtype: float64"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"formula = \"lnr ~ C(repeal)*C(year) + C(fip) + acc + ir + pi + alcohol + crack + poverty + income + ur\"\n",
"\n",
"reg = smf.ols(formula, data=abortion_bf15).fit(method='pinv')\n",
"\n",
"reg.params"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Aside: Effects of collinearity with SVD"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"rng = np.random.default_rng(seed=8675309)\n",
"\n",
"X = rng.uniform(0, 100, size=(100, 1))\n",
"beta = 1.5\n",
"\n",
"y = X.dot(beta) + rng.normal(scale=3.5, size=(100, 1))\n",
"\n",
"# make sure X is 2 dimensional\n",
"#X = "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Estimate beta using the SVD"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([[1.49117109]])"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"np.linalg.pinv(X).dot(y)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Double the X"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"XX = np.column_stack([X, X])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Aha"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([[0.74558555],\n",
" [0.74558555]])"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"np.linalg.pinv(XX).dot(y)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We'll come back to this."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Aside: Effects of collinearity with QR\n",
"\n",
"Not so clean. The column-pivoting rank-revealing QR decomposition used by estimatr can get around this by removing the column(s) causing the problem."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([ 1.18027578e+13, -1.18027578e+13])"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"sm.OLS(y, XX).fit(method='qr').params"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You can really see the collinearity issues if you try to fit the model using the QR decomposition. These are the coefficients that are sharing the weight of this effect in the default SVD (pinv) model."
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Intercept 7.29673318\n",
"C(repeal)[T.1.0] 3133303217854.80468750\n",
"C(year)[T.1986.0] -0.02604028\n",
"C(year)[T.1987.0] -0.16100534\n",
"C(year)[T.1988.0] -0.06935126\n",
"C(year)[T.1989.0] -0.10164296\n",
"C(year)[T.1990.0] -0.21590577\n",
"C(year)[T.1991.0] -0.26324314\n",
"C(year)[T.1992.0] -0.46611059\n",
"C(year)[T.1993.0] -0.57413943\n",
"C(year)[T.1994.0] -0.64322522\n",
"C(year)[T.1995.0] -0.80421520\n",
"C(year)[T.1996.0] -1.11997730\n",
"C(year)[T.1997.0] -1.13440763\n",
"C(year)[T.1998.0] -1.10079150\n",
"C(year)[T.1999.0] -1.13450157\n",
"C(year)[T.2000.0] -1.34122408\n",
"C(fip)[T.2.0] -3133303217856.74511719\n",
"C(fip)[T.4.0] -0.54606608\n",
"C(fip)[T.5.0] 0.11277881\n",
"C(fip)[T.6.0] -3133303217855.88867188\n",
"C(fip)[T.8.0] -0.67683213\n",
"C(fip)[T.9.0] -1.07977003\n",
"C(fip)[T.10.0] -0.40025487\n",
"C(fip)[T.11.0] -1.69513262\n",
"C(fip)[T.12.0] -0.79952382\n",
"C(fip)[T.13.0] -0.47528417\n",
"C(fip)[T.15.0] -3133303217856.61035156\n",
"C(fip)[T.16.0] -1.51872750\n",
"C(fip)[T.17.0] -0.50905375\n",
"C(fip)[T.18.0] 0.03854467\n",
"C(fip)[T.19.0] -0.02154076\n",
"C(fip)[T.20.0] 0.30723087\n",
"C(fip)[T.21.0] 0.08021348\n",
"C(fip)[T.22.0] -0.63548223\n",
"C(fip)[T.23.0] -2.17938423\n",
"C(fip)[T.24.0] -1.03385196\n",
"C(fip)[T.25.0] -1.30897829\n",
"C(fip)[T.26.0] -0.59375295\n",
"C(fip)[T.27.0] -0.13966967\n",
"C(fip)[T.28.0] -0.07078731\n",
"C(fip)[T.29.0] 0.07109987\n",
"C(fip)[T.30.0] -1.35232676\n",
"C(fip)[T.31.0] 0.00169566\n",
"C(fip)[T.32.0] -1.59525562\n",
"C(fip)[T.33.0] -3.51315049\n",
"C(fip)[T.34.0] -1.77336362\n",
"C(fip)[T.35.0] -0.85296674\n",
"C(fip)[T.36.0] -3133303217856.81494141\n",
"C(fip)[T.37.0] 0.05300960\n",
"C(fip)[T.38.0] -2.00472302\n",
"C(fip)[T.39.0] -0.19054343\n",
"C(fip)[T.40.0] 0.42397366\n",
"C(fip)[T.41.0] -0.66081709\n",
"C(fip)[T.42.0] -0.13920995\n",
"C(fip)[T.44.0] -0.54108121\n",
"C(fip)[T.45.0] -1.00759513\n",
"C(fip)[T.46.0] -1.35551917\n",
"C(fip)[T.47.0] 0.19596176\n",
"C(fip)[T.48.0] -0.42797090\n",
"C(fip)[T.49.0] -0.58223248\n",
"C(fip)[T.50.0] -1.93446186\n",
"C(fip)[T.51.0] -0.64849003\n",
"C(fip)[T.53.0] -3133303217855.45312500\n",
"C(fip)[T.54.0] -0.44065852\n",
"C(fip)[T.55.0] -0.16253513\n",
"C(fip)[T.56.0] -1.85977103\n",
"C(repeal)[T.1.0]:C(year)[T.1986.0] 0.00032629\n",
"C(repeal)[T.1.0]:C(year)[T.1987.0] -0.02826032\n",
"C(repeal)[T.1.0]:C(year)[T.1988.0] -0.27277086\n",
"C(repeal)[T.1.0]:C(year)[T.1989.0] -0.23979109\n",
"C(repeal)[T.1.0]:C(year)[T.1990.0] -0.15720672\n",
"C(repeal)[T.1.0]:C(year)[T.1991.0] -0.08992207\n",
"C(repeal)[T.1.0]:C(year)[T.1992.0] 0.02363056\n",
"C(repeal)[T.1.0]:C(year)[T.1993.0] 0.19710696\n",
"C(repeal)[T.1.0]:C(year)[T.1994.0] 0.27839444\n",
"C(repeal)[T.1.0]:C(year)[T.1995.0] 0.15275266\n",
"C(repeal)[T.1.0]:C(year)[T.1996.0] -0.14013719\n",
"C(repeal)[T.1.0]:C(year)[T.1997.0] 0.07051265\n",
"C(repeal)[T.1.0]:C(year)[T.1998.0] -0.32775278\n",
"C(repeal)[T.1.0]:C(year)[T.1999.0] -0.24110242\n",
"C(repeal)[T.1.0]:C(year)[T.2000.0] -0.08113890\n",
"acc 0.00164493\n",
"ir -0.00001790\n",
"pi -0.01423322\n",
"alcohol 0.45899571\n",
"crack -0.06621404\n",
"poverty -0.00510021\n",
"income 0.00005268\n",
"ur 0.00024773\n",
"Length: 90, dtype: float64"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"formula = \"lnr ~ C(repeal)*C(year) + C(fip) + acc + ir + pi + alcohol + crack + poverty + income + ur\"\n",
"\n",
"reg = smf.ols(formula, data=abortion_bf15).fit(method='qr')\n",
"\n",
"reg.params"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Weighted Model\n",
"\n",
"Two differences here. The first difference is that OLS in Python really is ordinary. It doesn't take weights. We have to use the weighted least-squares class WLS. Notice that we're still using the default `pinv` method to fit the model. The collinear columns are not the same as `estimatr`, but everything else should look good.\n",
"\n",
"The second difference is that we're using SVD vs the PQR method used by estimatr when it detects rank-deficiency."
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Intercept 7.47534301\n",
"C(repeal)[T.1.0] -1.35959607\n",
"C(year)[T.1986.0] -0.02082737\n",
"C(year)[T.1987.0] -0.28646648\n",
"C(year)[T.1988.0] -0.34946485\n",
"C(year)[T.1989.0] -0.40722001\n",
"C(year)[T.1990.0] -0.48641174\n",
"C(year)[T.1991.0] -0.53530339\n",
"C(year)[T.1992.0] -0.78662355\n",
"C(year)[T.1993.0] -0.97671663\n",
"C(year)[T.1994.0] -1.06436659\n",
"C(year)[T.1995.0] -1.35624511\n",
"C(year)[T.1996.0] -1.52006222\n",
"C(year)[T.1997.0] -1.57184605\n",
"C(year)[T.1998.0] -1.54512931\n",
"C(year)[T.1999.0] -1.58741354\n",
"C(year)[T.2000.0] -1.67257156\n",
"C(fip)[T.2.0] -0.47612414\n",
"C(fip)[T.4.0] -0.95893856\n",
"C(fip)[T.5.0] 0.08603103\n",
"C(fip)[T.6.0] 0.09818574\n",
"C(fip)[T.8.0] -0.82460495\n",
"C(fip)[T.9.0] -1.68089178\n",
"C(fip)[T.10.0] -0.74575439\n",
"C(fip)[T.11.0] -2.11235472\n",
"C(fip)[T.12.0] -1.09474761\n",
"C(fip)[T.13.0] -0.72157625\n",
"C(fip)[T.15.0] -0.38965365\n",
"C(fip)[T.16.0] -1.79371939\n",
"C(fip)[T.17.0] -0.72282047\n",
"C(fip)[T.18.0] -0.14315369\n",
"C(fip)[T.19.0] -0.14979696\n",
"C(fip)[T.20.0] 0.08004447\n",
"C(fip)[T.21.0] -0.11037085\n",
"C(fip)[T.22.0] -0.74395266\n",
"C(fip)[T.23.0] -2.02122539\n",
"C(fip)[T.24.0] -1.47456155\n",
"C(fip)[T.25.0] -1.89432424\n",
"C(fip)[T.26.0] -0.76799387\n",
"C(fip)[T.27.0] -0.32031988\n",
"C(fip)[T.28.0] -0.02758883\n",
"C(fip)[T.29.0] -0.16692162\n",
"C(fip)[T.30.0] -1.26616853\n",
"C(fip)[T.31.0] -0.11756566\n",
"C(fip)[T.32.0] -1.96480850\n",
"C(fip)[T.33.0] -3.59097336\n",
"C(fip)[T.34.0] -2.27555404\n",
"C(fip)[T.35.0] -1.18140339\n",
"C(fip)[T.36.0] -0.98777640\n",
"C(fip)[T.37.0] -0.08849572\n",
"C(fip)[T.38.0] -1.82649080\n",
"C(fip)[T.39.0] -0.48922418\n",
"C(fip)[T.40.0] 0.30486166\n",
"C(fip)[T.41.0] -1.12281228\n",
"C(fip)[T.42.0] -0.63106777\n",
"C(fip)[T.44.0] -1.20962341\n",
"C(fip)[T.45.0] -1.11089996\n",
"C(fip)[T.46.0] -2.02895172\n",
"C(fip)[T.47.0] 0.05924823\n",
"C(fip)[T.48.0] -0.60548288\n",
"C(fip)[T.49.0] -1.26672374\n",
"C(fip)[T.50.0] -1.99735806\n",
"C(fip)[T.51.0] -0.94936056\n",
"C(fip)[T.53.0] 0.39577237\n",
"C(fip)[T.54.0] -0.52137261\n",
"C(fip)[T.55.0] -0.48367762\n",
"C(fip)[T.56.0] -2.19210632\n",
"C(repeal)[T.1.0]:C(year)[T.1986.0] -0.25905382\n",
"C(repeal)[T.1.0]:C(year)[T.1987.0] -0.34173798\n",
"C(repeal)[T.1.0]:C(year)[T.1988.0] -0.61055324\n",
"C(repeal)[T.1.0]:C(year)[T.1989.0] -0.78500028\n",
"C(repeal)[T.1.0]:C(year)[T.1990.0] -0.63163148\n",
"C(repeal)[T.1.0]:C(year)[T.1991.0] -0.55260264\n",
"C(repeal)[T.1.0]:C(year)[T.1992.0] -0.44241661\n",
"C(repeal)[T.1.0]:C(year)[T.1993.0] -0.30606258\n",
"C(repeal)[T.1.0]:C(year)[T.1994.0] -0.11788972\n",
"C(repeal)[T.1.0]:C(year)[T.1995.0] 0.02102913\n",
"C(repeal)[T.1.0]:C(year)[T.1996.0] -0.12353414\n",
"C(repeal)[T.1.0]:C(year)[T.1997.0] 0.02085286\n",
"C(repeal)[T.1.0]:C(year)[T.1998.0] -0.03605240\n",
"C(repeal)[T.1.0]:C(year)[T.1999.0] 0.01470912\n",
"C(repeal)[T.1.0]:C(year)[T.2000.0] 0.04145081\n",
"acc 0.00287774\n",
"ir 0.00044389\n",
"pi -0.03902207\n",
"alcohol 0.44684213\n",
"crack 0.05282875\n",
"poverty -0.00276098\n",
"income 0.00005847\n",
"ur -0.02782371\n",
"Length: 90, dtype: float64"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"formula = (\n",
" \"lnr ~ C(repeal)*C(year) + C(fip)\"\n",
" \" + acc + ir + pi + alcohol + crack + poverty + income + ur\"\n",
")\n",
"\n",
"reg = (\n",
" smf\n",
" .wls(formula, data=abortion_bf15, weights=abortion_bf15.totpop.values)\n",
" .fit(\n",
" cov_type='cluster', \n",
" cov_kwds={'groups': abortion_bf15.fip.values}, \n",
" method='pinv')\n",
")\n",
"\n",
"reg.params"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Numerical Issues"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The covariance matrix estimation, however, has problems, and you'll notice that there's also a weight of 0 included so the weighted log-likelihood isn't defined. "
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"reg.model.weights.min()"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/Users/skipper/.miniconda3/lib/python3.7/site-packages/statsmodels/base/model.py:1834: ValueWarning: covariance of constraints does not have full rank. The number of constraints is 89, but rank is 27\n",
"/Users/skipper/.miniconda3/lib/python3.7/site-packages/statsmodels/regression/linear_model.py:764: RuntimeWarning: divide by zero encountered in log\n"
]
},
{
"data": {
"text/html": [
"<table class=\"simpletable\">\n",
"<caption>WLS Regression Results</caption>\n",
"<tr>\n",
" <th>Dep. Variable:</th> <td>lnr</td> <th> R-squared: </th> <td> 0.857</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Model:</th> <td>WLS</td> <th> Adj. R-squared: </th> <td> 0.837</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Method:</th> <td>Least Squares</td> <th> F-statistic: </th> <td> 1898.</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Date:</th> <td>Fri, 12 Mar 2021</td> <th> Prob (F-statistic):</th> <td>1.66e-66</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Time:</th> <td>09:22:28</td> <th> Log-Likelihood: </th> <td> -inf</td>\n",
"</tr>\n",
"<tr>\n",
" <th>No. Observations:</th> <td> 737</td> <th> AIC: </th> <td> inf</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Df Residuals:</th> <td> 648</td> <th> BIC: </th> <td> inf</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Df Model:</th> <td> 88</td> <th> </th> <td> </td> \n",
"</tr>\n",
"<tr>\n",
" <th>Covariance Type:</th> <td>cluster</td> <th> </th> <td> </td> \n",
"</tr>\n",
"</table>\n",
"<table class=\"simpletable\">\n",
"<tr>\n",
" <td></td> <th>coef</th> <th>std err</th> <th>z</th> <th>P>|z|</th> <th>[0.025</th> <th>0.975]</th> \n",
"</tr>\n",
"<tr>\n",
" <th>Intercept</th> <td> 7.4753</td> <td> 1.176</td> <td> 6.356</td> <td> 0.000</td> <td> 5.170</td> <td> 9.781</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(repeal)[T.1.0]</th> <td> -1.3596</td> <td> 0.412</td> <td> -3.296</td> <td> 0.001</td> <td> -2.168</td> <td> -0.551</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(year)[T.1986.0]</th> <td> -0.0208</td> <td> 0.058</td> <td> -0.361</td> <td> 0.718</td> <td> -0.134</td> <td> 0.092</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(year)[T.1987.0]</th> <td> -0.2865</td> <td> 0.111</td> <td> -2.579</td> <td> 0.010</td> <td> -0.504</td> <td> -0.069</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(year)[T.1988.0]</th> <td> -0.3495</td> <td> 0.131</td> <td> -2.667</td> <td> 0.008</td> <td> -0.606</td> <td> -0.093</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(year)[T.1989.0]</th> <td> -0.4072</td> <td> 0.173</td> <td> -2.349</td> <td> 0.019</td> <td> -0.747</td> <td> -0.067</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(year)[T.1990.0]</th> <td> -0.4864</td> <td> 0.217</td> <td> -2.244</td> <td> 0.025</td> <td> -0.911</td> <td> -0.062</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(year)[T.1991.0]</th> <td> -0.5353</td> <td> 0.223</td> <td> -2.398</td> <td> 0.016</td> <td> -0.973</td> <td> -0.098</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(year)[T.1992.0]</th> <td> -0.7866</td> <td> 0.269</td> <td> -2.927</td> <td> 0.003</td> <td> -1.313</td> <td> -0.260</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(year)[T.1993.0]</th> <td> -0.9767</td> <td> 0.287</td> <td> -3.403</td> <td> 0.001</td> <td> -1.539</td> <td> -0.414</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(year)[T.1994.0]</th> <td> -1.0644</td> <td> 0.323</td> <td> -3.297</td> <td> 0.001</td> <td> -1.697</td> <td> -0.432</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(year)[T.1995.0]</th> <td> -1.3562</td> <td> 0.375</td> <td> -3.612</td> <td> 0.000</td> <td> -2.092</td> <td> -0.620</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(year)[T.1996.0]</th> <td> -1.5201</td> <td> 0.415</td> <td> -3.667</td> <td> 0.000</td> <td> -2.333</td> <td> -0.707</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(year)[T.1997.0]</th> <td> -1.5718</td> <td> 0.468</td> <td> -3.359</td> <td> 0.001</td> <td> -2.489</td> <td> -0.655</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(year)[T.1998.0]</th> <td> -1.5451</td> <td> 0.536</td> <td> -2.884</td> <td> 0.004</td> <td> -2.595</td> <td> -0.495</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(year)[T.1999.0]</th> <td> -1.5874</td> <td> 0.579</td> <td> -2.743</td> <td> 0.006</td> <td> -2.722</td> <td> -0.453</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(year)[T.2000.0]</th> <td> -1.6726</td> <td> 0.658</td> <td> -2.543</td> <td> 0.011</td> <td> -2.962</td> <td> -0.383</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(fip)[T.2.0]</th> <td> -0.4761</td> <td> 0.267</td> <td> -1.784</td> <td> 0.074</td> <td> -0.999</td> <td> 0.047</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(fip)[T.4.0]</th> <td> -0.9589</td> <td> 0.524</td> <td> -1.829</td> <td> 0.067</td> <td> -1.986</td> <td> 0.068</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(fip)[T.5.0]</th> <td> 0.0860</td> <td> 0.065</td> <td> 1.319</td> <td> 0.187</td> <td> -0.042</td> <td> 0.214</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(fip)[T.6.0]</th> <td> 0.0982</td> <td> 0.192</td> <td> 0.512</td> <td> 0.608</td> <td> -0.277</td> <td> 0.474</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(fip)[T.8.0]</th> <td> -0.8246</td> <td> 0.418</td> <td> -1.971</td> <td> 0.049</td> <td> -1.645</td> <td> -0.005</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(fip)[T.9.0]</th> <td> -1.6809</td> <td> 0.736</td> <td> -2.285</td> <td> 0.022</td> <td> -3.122</td> <td> -0.239</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(fip)[T.10.0]</th> <td> -0.7458</td> <td> 0.484</td> <td> -1.540</td> <td> 0.123</td> <td> -1.695</td> <td> 0.203</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(fip)[T.11.0]</th> <td> -2.1124</td> <td> 1.246</td> <td> -1.695</td> <td> 0.090</td> <td> -4.555</td> <td> 0.331</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(fip)[T.12.0]</th> <td> -1.0947</td> <td> 0.514</td> <td> -2.132</td> <td> 0.033</td> <td> -2.101</td> <td> -0.088</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(fip)[T.13.0]</th> <td> -0.7216</td> <td> 0.246</td> <td> -2.930</td> <td> 0.003</td> <td> -1.204</td> <td> -0.239</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(fip)[T.15.0]</th> <td> -0.3897</td> <td> 0.104</td> <td> -3.740</td> <td> 0.000</td> <td> -0.594</td> <td> -0.185</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(fip)[T.16.0]</th> <td> -1.7937</td> <td> 0.199</td> <td> -8.994</td> <td> 0.000</td> <td> -2.185</td> <td> -1.403</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(fip)[T.17.0]</th> <td> -0.7228</td> <td> 0.459</td> <td> -1.575</td> <td> 0.115</td> <td> -1.622</td> <td> 0.177</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(fip)[T.18.0]</th> <td> -0.1432</td> <td> 0.187</td> <td> -0.765</td> <td> 0.444</td> <td> -0.510</td> <td> 0.224</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(fip)[T.19.0]</th> <td> -0.1498</td> <td> 0.196</td> <td> -0.764</td> <td> 0.445</td> <td> -0.534</td> <td> 0.235</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(fip)[T.20.0]</th> <td> 0.0800</td> <td> 0.218</td> <td> 0.368</td> <td> 0.713</td> <td> -0.346</td> <td> 0.506</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(fip)[T.21.0]</th> <td> -0.1104</td> <td> 0.059</td> <td> -1.866</td> <td> 0.062</td> <td> -0.226</td> <td> 0.006</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(fip)[T.22.0]</th> <td> -0.7440</td> <td> 0.283</td> <td> -2.629</td> <td> 0.009</td> <td> -1.299</td> <td> -0.189</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(fip)[T.23.0]</th> <td> -2.0212</td> <td> 0.148</td> <td> -13.617</td> <td> 0.000</td> <td> -2.312</td> <td> -1.730</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(fip)[T.24.0]</th> <td> -1.4746</td> <td> 0.501</td> <td> -2.943</td> <td> 0.003</td> <td> -2.457</td> <td> -0.492</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(fip)[T.25.0]</th> <td> -1.8943</td> <td> 0.538</td> <td> -3.519</td> <td> 0.000</td> <td> -2.949</td> <td> -0.839</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(fip)[T.26.0]</th> <td> -0.7680</td> <td> 0.361</td> <td> -2.129</td> <td> 0.033</td> <td> -1.475</td> <td> -0.061</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(fip)[T.27.0]</th> <td> -0.3203</td> <td> 0.403</td> <td> -0.794</td> <td> 0.427</td> <td> -1.111</td> <td> 0.470</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(fip)[T.28.0]</th> <td> -0.0276</td> <td> 0.146</td> <td> -0.189</td> <td> 0.850</td> <td> -0.313</td> <td> 0.258</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(fip)[T.29.0]</th> <td> -0.1669</td> <td> 0.286</td> <td> -0.583</td> <td> 0.560</td> <td> -0.728</td> <td> 0.394</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(fip)[T.30.0]</th> <td> -1.2662</td> <td> 0.188</td> <td> -6.719</td> <td> 0.000</td> <td> -1.636</td> <td> -0.897</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(fip)[T.31.0]</th> <td> -0.1176</td> <td> 0.308</td> <td> -0.381</td> <td> 0.703</td> <td> -0.722</td> <td> 0.487</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(fip)[T.32.0]</th> <td> -1.9648</td> <td> 1.045</td> <td> -1.880</td> <td> 0.060</td> <td> -4.013</td> <td> 0.084</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(fip)[T.33.0]</th> <td> -3.5910</td> <td> 0.874</td> <td> -4.108</td> <td> 0.000</td> <td> -5.304</td> <td> -1.878</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(fip)[T.34.0]</th> <td> -2.2756</td> <td> 0.646</td> <td> -3.525</td> <td> 0.000</td> <td> -3.541</td> <td> -1.010</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(fip)[T.35.0]</th> <td> -1.1814</td> <td> 0.346</td> <td> -3.419</td> <td> 0.001</td> <td> -1.859</td> <td> -0.504</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(fip)[T.36.0]</th> <td> -0.9878</td> <td> 0.139</td> <td> -7.105</td> <td> 0.000</td> <td> -1.260</td> <td> -0.715</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(fip)[T.37.0]</th> <td> -0.0885</td> <td> 0.128</td> <td> -0.689</td> <td> 0.491</td> <td> -0.340</td> <td> 0.163</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(fip)[T.38.0]</th> <td> -1.8265</td> <td> 0.135</td> <td> -13.543</td> <td> 0.000</td> <td> -2.091</td> <td> -1.562</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(fip)[T.39.0]</th> <td> -0.4892</td> <td> 0.247</td> <td> -1.980</td> <td> 0.048</td> <td> -0.973</td> <td> -0.005</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(fip)[T.40.0]</th> <td> 0.3049</td> <td> 0.073</td> <td> 4.182</td> <td> 0.000</td> <td> 0.162</td> <td> 0.448</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(fip)[T.41.0]</th> <td> -1.1228</td> <td> 0.345</td> <td> -3.252</td> <td> 0.001</td> <td> -1.800</td> <td> -0.446</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(fip)[T.42.0]</th> <td> -0.6311</td> <td> 0.330</td> <td> -1.914</td> <td> 0.056</td> <td> -1.277</td> <td> 0.015</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(fip)[T.44.0]</th> <td> -1.2096</td> <td> 0.522</td> <td> -2.317</td> <td> 0.021</td> <td> -2.233</td> <td> -0.186</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(fip)[T.45.0]</th> <td> -1.1109</td> <td> 0.189</td> <td> -5.876</td> <td> 0.000</td> <td> -1.481</td> <td> -0.740</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(fip)[T.46.0]</th> <td> -2.0290</td> <td> 0.702</td> <td> -2.889</td> <td> 0.004</td> <td> -3.405</td> <td> -0.653</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(fip)[T.47.0]</th> <td> 0.0592</td> <td> 0.103</td> <td> 0.574</td> <td> 0.566</td> <td> -0.143</td> <td> 0.262</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(fip)[T.48.0]</th> <td> -0.6055</td> <td> 0.327</td> <td> -1.849</td> <td> 0.064</td> <td> -1.247</td> <td> 0.036</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(fip)[T.49.0]</th> <td> -1.2667</td> <td> 0.207</td> <td> -6.134</td> <td> 0.000</td> <td> -1.671</td> <td> -0.862</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(fip)[T.50.0]</th> <td> -1.9974</td> <td> 0.220</td> <td> -9.067</td> <td> 0.000</td> <td> -2.429</td> <td> -1.566</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(fip)[T.51.0]</th> <td> -0.9494</td> <td> 0.325</td> <td> -2.923</td> <td> 0.003</td> <td> -1.586</td> <td> -0.313</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(fip)[T.53.0]</th> <td> 0.3958</td> <td> 0.061</td> <td> 6.534</td> <td> 0.000</td> <td> 0.277</td> <td> 0.514</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(fip)[T.54.0]</th> <td> -0.5214</td> <td> 0.127</td> <td> -4.094</td> <td> 0.000</td> <td> -0.771</td> <td> -0.272</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(fip)[T.55.0]</th> <td> -0.4837</td> <td> 0.552</td> <td> -0.876</td> <td> 0.381</td> <td> -1.566</td> <td> 0.598</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(fip)[T.56.0]</th> <td> -2.1921</td> <td> 0.583</td> <td> -3.759</td> <td> 0.000</td> <td> -3.335</td> <td> -1.049</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(repeal)[T.1.0]:C(year)[T.1986.0]</th> <td> -0.2591</td> <td> 0.063</td> <td> -4.089</td> <td> 0.000</td> <td> -0.383</td> <td> -0.135</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(repeal)[T.1.0]:C(year)[T.1987.0]</th> <td> -0.3417</td> <td> 0.168</td> <td> -2.029</td> <td> 0.042</td> <td> -0.672</td> <td> -0.012</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(repeal)[T.1.0]:C(year)[T.1988.0]</th> <td> -0.6106</td> <td> 0.203</td> <td> -3.011</td> <td> 0.003</td> <td> -1.008</td> <td> -0.213</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(repeal)[T.1.0]:C(year)[T.1989.0]</th> <td> -0.7850</td> <td> 0.237</td> <td> -3.307</td> <td> 0.001</td> <td> -1.250</td> <td> -0.320</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(repeal)[T.1.0]:C(year)[T.1990.0]</th> <td> -0.6316</td> <td> 0.176</td> <td> -3.594</td> <td> 0.000</td> <td> -0.976</td> <td> -0.287</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(repeal)[T.1.0]:C(year)[T.1991.0]</th> <td> -0.5526</td> <td> 0.139</td> <td> -3.963</td> <td> 0.000</td> <td> -0.826</td> <td> -0.279</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(repeal)[T.1.0]:C(year)[T.1992.0]</th> <td> -0.4424</td> <td> 0.160</td> <td> -2.757</td> <td> 0.006</td> <td> -0.757</td> <td> -0.128</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(repeal)[T.1.0]:C(year)[T.1993.0]</th> <td> -0.3061</td> <td> 0.181</td> <td> -1.692</td> <td> 0.091</td> <td> -0.661</td> <td> 0.048</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(repeal)[T.1.0]:C(year)[T.1994.0]</th> <td> -0.1179</td> <td> 0.207</td> <td> -0.570</td> <td> 0.569</td> <td> -0.523</td> <td> 0.287</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(repeal)[T.1.0]:C(year)[T.1995.0]</th> <td> 0.0210</td> <td> 0.223</td> <td> 0.094</td> <td> 0.925</td> <td> -0.416</td> <td> 0.458</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(repeal)[T.1.0]:C(year)[T.1996.0]</th> <td> -0.1235</td> <td> 0.207</td> <td> -0.598</td> <td> 0.550</td> <td> -0.529</td> <td> 0.282</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(repeal)[T.1.0]:C(year)[T.1997.0]</th> <td> 0.0209</td> <td> 0.264</td> <td> 0.079</td> <td> 0.937</td> <td> -0.497</td> <td> 0.539</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(repeal)[T.1.0]:C(year)[T.1998.0]</th> <td> -0.0361</td> <td> 0.348</td> <td> -0.104</td> <td> 0.917</td> <td> -0.717</td> <td> 0.645</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(repeal)[T.1.0]:C(year)[T.1999.0]</th> <td> 0.0147</td> <td> 0.358</td> <td> 0.041</td> <td> 0.967</td> <td> -0.686</td> <td> 0.716</td>\n",
"</tr>\n",
"<tr>\n",
" <th>C(repeal)[T.1.0]:C(year)[T.2000.0]</th> <td> 0.0415</td> <td> 0.388</td> <td> 0.107</td> <td> 0.915</td> <td> -0.718</td> <td> 0.801</td>\n",
"</tr>\n",
"<tr>\n",
" <th>acc</th> <td> 0.0029</td> <td> 0.001</td> <td> 2.247</td> <td> 0.025</td> <td> 0.000</td> <td> 0.005</td>\n",
"</tr>\n",
"<tr>\n",
" <th>ir</th> <td> 0.0004</td> <td> 0.000</td> <td> 1.056</td> <td> 0.291</td> <td> -0.000</td> <td> 0.001</td>\n",
"</tr>\n",
"<tr>\n",
" <th>pi</th> <td> -0.0390</td> <td> 0.066</td> <td> -0.589</td> <td> 0.556</td> <td> -0.169</td> <td> 0.091</td>\n",
"</tr>\n",
"<tr>\n",
" <th>alcohol</th> <td> 0.4468</td> <td> 0.328</td> <td> 1.362</td> <td> 0.173</td> <td> -0.196</td> <td> 1.090</td>\n",
"</tr>\n",
"<tr>\n",
" <th>crack</th> <td> 0.0528</td> <td> 0.035</td> <td> 1.525</td> <td> 0.127</td> <td> -0.015</td> <td> 0.121</td>\n",
"</tr>\n",
"<tr>\n",
" <th>poverty</th> <td> -0.0028</td> <td> 0.014</td> <td> -0.197</td> <td> 0.844</td> <td> -0.030</td> <td> 0.025</td>\n",
"</tr>\n",
"<tr>\n",
" <th>income</th> <td> 5.847e-05</td> <td> 4.43e-05</td> <td> 1.320</td> <td> 0.187</td> <td>-2.84e-05</td> <td> 0.000</td>\n",
"</tr>\n",
"<tr>\n",
" <th>ur</th> <td> -0.0278</td> <td> 0.037</td> <td> -0.762</td> <td> 0.446</td> <td> -0.099</td> <td> 0.044</td>\n",
"</tr>\n",
"</table>\n",
"<table class=\"simpletable\">\n",
"<tr>\n",
" <th>Omnibus:</th> <td>368.324</td> <th> Durbin-Watson: </th> <td> 0.949</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Prob(Omnibus):</th> <td> 0.000</td> <th> Jarque-Bera (JB): </th> <td>8710.830</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Skew:</th> <td>-1.715</td> <th> Prob(JB): </th> <td> 0.00</td>\n",
"</tr>\n",
"<tr>\n",
" <th>Kurtosis:</th> <td>19.489</td> <th> Cond. No. </th> <td>1.03e+16</td>\n",
"</tr>\n",
"</table><br/><br/>Notes:<br/>[1] Standard Errors are robust to cluster correlation (cluster)<br/>[2] The smallest eigenvalue is 1.95e-16. This might indicate that there are<br/>strong multicollinearity problems or that the design matrix is singular."
],
"text/plain": [
"<class 'statsmodels.iolib.summary.Summary'>\n",
"\"\"\"\n",
" WLS Regression Results \n",
"==============================================================================\n",
"Dep. Variable: lnr R-squared: 0.857\n",
"Model: WLS Adj. R-squared: 0.837\n",
"Method: Least Squares F-statistic: 1898.\n",
"Date: Fri, 12 Mar 2021 Prob (F-statistic): 1.66e-66\n",
"Time: 09:22:28 Log-Likelihood: -inf\n",
"No. Observations: 737 AIC: inf\n",
"Df Residuals: 648 BIC: inf\n",
"Df Model: 88 \n",
"Covariance Type: cluster \n",
"======================================================================================================\n",
" coef std err z P>|z| [0.025 0.975]\n",
"------------------------------------------------------------------------------------------------------\n",
"Intercept 7.4753 1.176 6.356 0.000 5.170 9.781\n",
"C(repeal)[T.1.0] -1.3596 0.412 -3.296 0.001 -2.168 -0.551\n",
"C(year)[T.1986.0] -0.0208 0.058 -0.361 0.718 -0.134 0.092\n",
"C(year)[T.1987.0] -0.2865 0.111 -2.579 0.010 -0.504 -0.069\n",
"C(year)[T.1988.0] -0.3495 0.131 -2.667 0.008 -0.606 -0.093\n",
"C(year)[T.1989.0] -0.4072 0.173 -2.349 0.019 -0.747 -0.067\n",
"C(year)[T.1990.0] -0.4864 0.217 -2.244 0.025 -0.911 -0.062\n",
"C(year)[T.1991.0] -0.5353 0.223 -2.398 0.016 -0.973 -0.098\n",
"C(year)[T.1992.0] -0.7866 0.269 -2.927 0.003 -1.313 -0.260\n",
"C(year)[T.1993.0] -0.9767 0.287 -3.403 0.001 -1.539 -0.414\n",
"C(year)[T.1994.0] -1.0644 0.323 -3.297 0.001 -1.697 -0.432\n",
"C(year)[T.1995.0] -1.3562 0.375 -3.612 0.000 -2.092 -0.620\n",
"C(year)[T.1996.0] -1.5201 0.415 -3.667 0.000 -2.333 -0.707\n",
"C(year)[T.1997.0] -1.5718 0.468 -3.359 0.001 -2.489 -0.655\n",
"C(year)[T.1998.0] -1.5451 0.536 -2.884 0.004 -2.595 -0.495\n",
"C(year)[T.1999.0] -1.5874 0.579 -2.743 0.006 -2.722 -0.453\n",
"C(year)[T.2000.0] -1.6726 0.658 -2.543 0.011 -2.962 -0.383\n",
"C(fip)[T.2.0] -0.4761 0.267 -1.784 0.074 -0.999 0.047\n",
"C(fip)[T.4.0] -0.9589 0.524 -1.829 0.067 -1.986 0.068\n",
"C(fip)[T.5.0] 0.0860 0.065 1.319 0.187 -0.042 0.214\n",
"C(fip)[T.6.0] 0.0982 0.192 0.512 0.608 -0.277 0.474\n",
"C(fip)[T.8.0] -0.8246 0.418 -1.971 0.049 -1.645 -0.005\n",
"C(fip)[T.9.0] -1.6809 0.736 -2.285 0.022 -3.122 -0.239\n",
"C(fip)[T.10.0] -0.7458 0.484 -1.540 0.123 -1.695 0.203\n",
"C(fip)[T.11.0] -2.1124 1.246 -1.695 0.090 -4.555 0.331\n",
"C(fip)[T.12.0] -1.0947 0.514 -2.132 0.033 -2.101 -0.088\n",
"C(fip)[T.13.0] -0.7216 0.246 -2.930 0.003 -1.204 -0.239\n",
"C(fip)[T.15.0] -0.3897 0.104 -3.740 0.000 -0.594 -0.185\n",
"C(fip)[T.16.0] -1.7937 0.199 -8.994 0.000 -2.185 -1.403\n",
"C(fip)[T.17.0] -0.7228 0.459 -1.575 0.115 -1.622 0.177\n",
"C(fip)[T.18.0] -0.1432 0.187 -0.765 0.444 -0.510 0.224\n",
"C(fip)[T.19.0] -0.1498 0.196 -0.764 0.445 -0.534 0.235\n",
"C(fip)[T.20.0] 0.0800 0.218 0.368 0.713 -0.346 0.506\n",
"C(fip)[T.21.0] -0.1104 0.059 -1.866 0.062 -0.226 0.006\n",
"C(fip)[T.22.0] -0.7440 0.283 -2.629 0.009 -1.299 -0.189\n",
"C(fip)[T.23.0] -2.0212 0.148 -13.617 0.000 -2.312 -1.730\n",
"C(fip)[T.24.0] -1.4746 0.501 -2.943 0.003 -2.457 -0.492\n",
"C(fip)[T.25.0] -1.8943 0.538 -3.519 0.000 -2.949 -0.839\n",
"C(fip)[T.26.0] -0.7680 0.361 -2.129 0.033 -1.475 -0.061\n",
"C(fip)[T.27.0] -0.3203 0.403 -0.794 0.427 -1.111 0.470\n",
"C(fip)[T.28.0] -0.0276 0.146 -0.189 0.850 -0.313 0.258\n",
"C(fip)[T.29.0] -0.1669 0.286 -0.583 0.560 -0.728 0.394\n",
"C(fip)[T.30.0] -1.2662 0.188 -6.719 0.000 -1.636 -0.897\n",
"C(fip)[T.31.0] -0.1176 0.308 -0.381 0.703 -0.722 0.487\n",
"C(fip)[T.32.0] -1.9648 1.045 -1.880 0.060 -4.013 0.084\n",
"C(fip)[T.33.0] -3.5910 0.874 -4.108 0.000 -5.304 -1.878\n",
"C(fip)[T.34.0] -2.2756 0.646 -3.525 0.000 -3.541 -1.010\n",
"C(fip)[T.35.0] -1.1814 0.346 -3.419 0.001 -1.859 -0.504\n",
"C(fip)[T.36.0] -0.9878 0.139 -7.105 0.000 -1.260 -0.715\n",
"C(fip)[T.37.0] -0.0885 0.128 -0.689 0.491 -0.340 0.163\n",
"C(fip)[T.38.0] -1.8265 0.135 -13.543 0.000 -2.091 -1.562\n",
"C(fip)[T.39.0] -0.4892 0.247 -1.980 0.048 -0.973 -0.005\n",
"C(fip)[T.40.0] 0.3049 0.073 4.182 0.000 0.162 0.448\n",
"C(fip)[T.41.0] -1.1228 0.345 -3.252 0.001 -1.800 -0.446\n",
"C(fip)[T.42.0] -0.6311 0.330 -1.914 0.056 -1.277 0.015\n",
"C(fip)[T.44.0] -1.2096 0.522 -2.317 0.021 -2.233 -0.186\n",
"C(fip)[T.45.0] -1.1109 0.189 -5.876 0.000 -1.481 -0.740\n",
"C(fip)[T.46.0] -2.0290 0.702 -2.889 0.004 -3.405 -0.653\n",
"C(fip)[T.47.0] 0.0592 0.103 0.574 0.566 -0.143 0.262\n",
"C(fip)[T.48.0] -0.6055 0.327 -1.849 0.064 -1.247 0.036\n",
"C(fip)[T.49.0] -1.2667 0.207 -6.134 0.000 -1.671 -0.862\n",
"C(fip)[T.50.0] -1.9974 0.220 -9.067 0.000 -2.429 -1.566\n",
"C(fip)[T.51.0] -0.9494 0.325 -2.923 0.003 -1.586 -0.313\n",
"C(fip)[T.53.0] 0.3958 0.061 6.534 0.000 0.277 0.514\n",
"C(fip)[T.54.0] -0.5214 0.127 -4.094 0.000 -0.771 -0.272\n",
"C(fip)[T.55.0] -0.4837 0.552 -0.876 0.381 -1.566 0.598\n",
"C(fip)[T.56.0] -2.1921 0.583 -3.759 0.000 -3.335 -1.049\n",
"C(repeal)[T.1.0]:C(year)[T.1986.0] -0.2591 0.063 -4.089 0.000 -0.383 -0.135\n",
"C(repeal)[T.1.0]:C(year)[T.1987.0] -0.3417 0.168 -2.029 0.042 -0.672 -0.012\n",
"C(repeal)[T.1.0]:C(year)[T.1988.0] -0.6106 0.203 -3.011 0.003 -1.008 -0.213\n",
"C(repeal)[T.1.0]:C(year)[T.1989.0] -0.7850 0.237 -3.307 0.001 -1.250 -0.320\n",
"C(repeal)[T.1.0]:C(year)[T.1990.0] -0.6316 0.176 -3.594 0.000 -0.976 -0.287\n",
"C(repeal)[T.1.0]:C(year)[T.1991.0] -0.5526 0.139 -3.963 0.000 -0.826 -0.279\n",
"C(repeal)[T.1.0]:C(year)[T.1992.0] -0.4424 0.160 -2.757 0.006 -0.757 -0.128\n",
"C(repeal)[T.1.0]:C(year)[T.1993.0] -0.3061 0.181 -1.692 0.091 -0.661 0.048\n",
"C(repeal)[T.1.0]:C(year)[T.1994.0] -0.1179 0.207 -0.570 0.569 -0.523 0.287\n",
"C(repeal)[T.1.0]:C(year)[T.1995.0] 0.0210 0.223 0.094 0.925 -0.416 0.458\n",
"C(repeal)[T.1.0]:C(year)[T.1996.0] -0.1235 0.207 -0.598 0.550 -0.529 0.282\n",
"C(repeal)[T.1.0]:C(year)[T.1997.0] 0.0209 0.264 0.079 0.937 -0.497 0.539\n",
"C(repeal)[T.1.0]:C(year)[T.1998.0] -0.0361 0.348 -0.104 0.917 -0.717 0.645\n",
"C(repeal)[T.1.0]:C(year)[T.1999.0] 0.0147 0.358 0.041 0.967 -0.686 0.716\n",
"C(repeal)[T.1.0]:C(year)[T.2000.0] 0.0415 0.388 0.107 0.915 -0.718 0.801\n",
"acc 0.0029 0.001 2.247 0.025 0.000 0.005\n",
"ir 0.0004 0.000 1.056 0.291 -0.000 0.001\n",
"pi -0.0390 0.066 -0.589 0.556 -0.169 0.091\n",
"alcohol 0.4468 0.328 1.362 0.173 -0.196 1.090\n",
"crack 0.0528 0.035 1.525 0.127 -0.015 0.121\n",
"poverty -0.0028 0.014 -0.197 0.844 -0.030 0.025\n",
"income 5.847e-05 4.43e-05 1.320 0.187 -2.84e-05 0.000\n",
"ur -0.0278 0.037 -0.762 0.446 -0.099 0.044\n",
"==============================================================================\n",
"Omnibus: 368.324 Durbin-Watson: 0.949\n",
"Prob(Omnibus): 0.000 Jarque-Bera (JB): 8710.830\n",
"Skew: -1.715 Prob(JB): 0.00\n",
"Kurtosis: 19.489 Cond. No. 1.03e+16\n",
"==============================================================================\n",
"\n",
"Notes:\n",
"[1] Standard Errors are robust to cluster correlation (cluster)\n",
"[2] The smallest eigenvalue is 1.95e-16. This might indicate that there are\n",
"strong multicollinearity problems or that the design matrix is singular.\n",
"\"\"\""
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"reg.summary()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Recreate the Plot"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"At this point, though, you can re-create the _point estimates_ of the plot. The small thinko here was that you were using the main effects for the years and not those interacted with `repeal1`! See how I changed the indexing below? I'm not a ggplot person, so I didn't mess with the graph below."
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAj0AAAHICAYAAAClJls2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB5qUlEQVR4nO3dd1QU198G8GfpTbCggCAIdqMoltgVFCyYKPbeRTQxDSXWqESDURNLbBg1Yi8xtl9iRcWSqIm9Y2IDFRRQWJC+zPsH7LwsC4KwsMA8n3M4urN3Z+6dvTv73Tu3yARBEEBERERUzuloOwNEREREJYFBDxEREUkCgx4iIiKSBAY9REREJAkMeoiIiEgSGPQQERGRJDDoISIiIklg0ENERESSwKCHiIiIJIFBTwmQyWQYPXq0trNRIlxdXVGzZs1i2ffZs2fRpk0bVKhQATKZDEFBQQCAmJgYjB49GtWrV4dMJoOrq2uxHL+kJCUl4auvvoKDgwN0dXU1ej6L8/0pipo1a2r9fXvy5AlkMhnmzZun1XzkJbdzVBrez3nz5kEmk+HJkydazQcVr9L++SgoPW1noKSFhITAzc3tnWkeP378XheSJ0+eICgoCF5eXmjatGnRMliMli9fjooVK5bJAOzNmzfw8vKCra0tfvzxR5iYmKBt27YAgClTpmDnzp2YNWsWnJycYGVlVWz5CAkJQUhICL788ktUrFixWI6xaNEiLF++HFOnToWzszMqVKhQLMcpSbGxsVi+fDlcXV21HtzQ+wsKCkJsbCy+/PJLbWeFqEgkF/Qo9e/fH7179871uapVq77Xvp48eQJ/f3/UrFkz16AnKSkJurq6hcmmRi1fvhw1a9Ysk0HPP//8gzdv3mDjxo3o06ePynPHjx9H9+7dMWfOnGLPR0hICPz9/TF69OhiC3qOHz8OZ2dnLFmypFj2rw2xsbHw9/cHgFyDntDQUMhkshLOVdl3/PhxlMTyiUFBQXjy5EmuQc/s2bMxffp0GBoaFns+iIpKskFPkyZNMHz48BI5lpGRUYkcpzyLjIwEAFSqVCnX53LbXlZFRkbC3t5e29nQiISEBJiZmeWbjl+YhWNgYKDtLEBPTw96epL9KikX0tPToVAocv0cCoKAxMRELeSqeLBPTz4uXryIjz/+GNWrV4ehoSGqVauGTp064dChQwAy72crb5eNGTMGMplMrV9Jbn16lNvOnDmD9u3bw9TUFFZWVpg2bRoUCgVSUlIwffp01KhRA4aGhmjRogUuXLigso+MjAwEBATA1dUVNjY2MDAwgK2tLUaNGoWwsDAxnfJe7NOnT3HmzBkxjznvwz98+FDsG2NgYAA7Ozt88skniI6OVjsv//33H/r06QMLCwtUqFABXbt2xY0bN977/J4+fRo9evRApUqVYGhoiAYNGmDRokVQKBQq52rUqFEAADc3NzHvo0ePhkwmgyAI2Lx5s7hd2denoPtXevz4Mby9veHg4ABDQ0NYWVmha9euOHHiBIDMFgpla4Wjo6N4vILc437+/DnGjx8PW1tb8dxOmDABERERYhpl34jHjx+rvE/57X/37t3w8vKCg4MDjIyMULlyZXTv3h3nz5/P8zWPHz9G3759UalSJZiZmaFr1664fv26WrqMjAz89NNPaNKkCYyNjWFubo7OnTuL5yQ7ZZ+TmzdvomfPnqhUqRIqVKiAoKAgODo6AgD8/f3FcmW/hZxXn56jR4/Czc0N5ubmMDY2RtOmTbF69Wq11g1lXYiLi8OkSZNgZWUFQ0NDNG/eHMePH3/n+SuI3377DZ06dRLz4eLigg0bNuSaduvWrXB2doahoSFsbW3h6+uLu3fvqr2XBf38vkvOPj1BQUEqn++cf9mPv3btWnTr1g12dnYwMDBAtWrV0K9fP9y+fVvlGDKZDGfOnMHTp09V9hUSEgIg7z49BanzQGbrqfJzGxQUhMaNG8PIyAi2traYNWtWrp/VvJR0fZHL5Zg8eTKsra1hbGyMZs2a4ddff83znNy/fx+DBw8Wj+fk5ISpU6dCLperpFO+j6dOncIPP/yAOnXqiOmXLVumkbLHxMRgwoQJsLGxgaGhIS5cuCAeNzg4GAsXLkTdunVhaGio1ur8+++/48MPP4SxsTGqVauGiRMn5hoYvXz5Ep999hlq1qwJAwMDWFlZYfjw4WrnJT4+Ht988w1at26NqlWrwsDAADVr1sTkyZPx+vXrAr4bBSPZ8DwxMTHXL3M9PT3xtsWDBw/g7u6OqlWr4tNPP4WNjQ2ioqJw+fJlXLx4Eb169ULfvn2RlpaGgIAATJgwAR06dACAAvUruXbtGg4dOoTx48dj+PDhOHr0KBYvXgxdXV3cunULcrkcU6dOxdu3b/Hjjz/io48+wpMnT8Q+HqmpqVi8eDH69++Pjz/+GBUqVMDNmzfxyy+/4OTJk7h58yYqV66MqlWrYuvWrfjqq69gaWmJWbNmiXlQ3sq7fv06XF1dYWJigrFjx8LBwQH//vsv1q5di5MnT+Lvv/+GhYUFACA8PBxt27ZFbGwsJk6ciAYNGuDChQtwdXVFlSpVCvwe/PLLLxg/fjxcXFwwffp0VKxYEX/++SdmzJiBa9euYdeuXQAyv0TOnTuHn3/+GTNnzkSDBg0AALVq1YK7uztGjBiBDh06YMKECQAg9vUp6P4B4OrVq+jSpQsSExMxduxYNGnSBHK5HBcvXkRwcDA8PDwwa9YsVK5cGfv378eyZctgaWkJAHB2dn5nOZ8/f46WLVvi1atX8Pb2hrOzM27cuIH169fj6NGj+Oeff2BlZYW+ffuidu3aau9TfvtftWoVLC0tMWHCBFhZWSE8PBwbN26Em5sbzpw5I54Ppbdv38LV1RXNmjXDggULEB4ejtWrV6NDhw64cOECGjVqJKYdPXo0tm7divbt22PhwoWIj4/Hhg0b0K1bN2zZskWttTQ8PBydOnWCl5cXFi5ciMjISHTs2BHLli3DV199hT59+qBv374AkG8L0MaNG+Ht7Y2aNWvCz88PZmZm2Lt3LyZPnowbN27g559/VntNt27dUKlSJcyaNQuJiYlYvnw5Pv74Y/z777+Fbj2bN28e/P394ebmhrlz58LY2BjHjh2Dt7c3/vvvP3z//fdi2tWrV2Py5MmoX78+/P39oa+vj507d+LcuXNq+y3o5/d9dOzYEVu3blXbvnbtWvz111+wtrYWty1ZsgRt2rTBZ599hipVquDBgwfYsGEDTpw4gWvXrqFWrVoAMj9/3333HaKjo1W+cJWfw9wUtM5nFxgYiBcvXmD8+PGoWrUq9u3bh4CAAFSoUAHTp0/Pt+wlXV/S09PRvXt3XLhwAQMGDICrqyvCw8Mxbtw41KtXTy399evX0bFjRygUCkyaNAlOTk44f/48fvzxRwQHB+Ovv/6CiYmJymtmzpyJ+Ph4jB07FmZmZtiyZQt8fX1hY2ODwYMHF6ns7u7uqFy5MqZPn46MjAxYW1uLwYifnx/evn2LkSNHolq1aqhRo4b4uiNHjmDVqlXw8fHB2LFjcfLkSaxbtw5A5nuopPyeSEhIwLhx41C3bl08f/4ca9euxfHjx3H58mXxHD9//hzr169H//79MWTIEBgaGuLSpUtYt24dzp8/j3/++Qf6+vrvfD8KTJCY06dPCwDy/KtXr56YdsWKFQIA4dKlSwXa56ZNm3J9HoAwatQotW0ymUz4888/Vba7uLgIMplM6Nmzp5CRkSFu379/vwBACAwMFLdlZGQIiYmJasc7ceKEAEBYvHixynYHBwehU6dOueaxadOmgqOjoxATE6Oy/eLFi4Kurq4wb948cduIESMEAML+/ftV0i5cuFAAIDg4OOR6jOwiIiIEIyMjwcvLS6WcgiAIP/zwgwBACAkJEbdt2rRJACCcPn1abV+5nd/32X9GRobQqFEjQU9PT/jnn3/U9q9QKMT/z507VwAgPH78ON8yKinP1+7du1W2b968WQAgjBs3TmX7u96n3CQkJKhti4iIEKpUqSJ4enqqbO/UqZMAQPj0009Vtl+6dEmQyWRCly5dxG0nT54UAAgfffSRkJ6eLm5/9eqVUK1aNaFixYpCfHy8Sr4BCGvXrlXLz+PHjwUAwty5c3MtQ84yx8bGCmZmZoKtra1KnUxLSxM8PDwEAMK5c+fE7aNGjRIACBMmTFDZ74ULFwQAwowZM3I9bn55vHr1qiCTyYTPP/9cLf3kyZMFHR0d4eHDh4IgCMKbN28EU1NTwcnJSZDL5WK65ORkoWXLlmr71sTnt1OnTvl+3pT1bNCgQSqfhdzqze3btwV9fX3hk08+KfBxcvtMvE+dV14/ra2thdevX4vbFQqF0KBBA8HGxuad5RME7dSXn3/+WQAg+Pn5qWy/fPmyIJPJ1M5Jhw4dBJlMJly4cEElvb+/vwBAmD9/vrhNeb1zdnYWkpOTxe0JCQlClSpVhDZt2hS57IMHD1a7NiqPW6tWLZXPtiD8/+fD2NhYrPNK3bp1E/T19VXqVO/evYVKlSqppX38+LFgZmYmjB49WtyWkpIipKamCjmtX79eACDs2bNH7bnCkuztrdGjR+PEiRNqf5s2bRLTKPuJ7N+/H0lJSRrPQ5s2bdR+hXfo0AGCIOCLL75Q6djZqVMnAMC///4rbpPJZDA2NgaQ2VQeGxuL6OhoNG3aFBYWFrh06VKB8nH79m1cv34dgwcPRkZGBqKjo8W/WrVqoXbt2jh27Jh4nAMHDqBhw4bw8vJS2c+XX35ZoP4bALB3714kJydj/PjxiImJUTnmRx99BADiMQvjffZ/48YN3L59G8OHD0eLFi3U9qWjU/iPifJ81a9fHwMHDlR5bsSIEahVqxb27dtXpM6opqam4v/j4+MRExMDPT09tGrVKs86MHPmTJXHH374Idzd3XHq1Cm8efMGQOYtHSCzo2r2jvjKls/Y2FicPHlSZT+VK1eGt7d3ocuidPz4cSQkJOCzzz5Tae3Q09PD7NmzVfKX3dSpU1Uet27dGmZmZnjw4EGh8rF9+3YIgoBx48ap1KHo6Gj06tULGRkZCA4OFvP89u1bfPLJJyoj7gwNDeHr66u2b019ft/lxIkTGD9+PDp27CjeAlZS1htBECCXyxEdHQ0rKyvUq1evSMcubJ0fO3asSt88HR0ddOnSBREREUhISHjnMbVRX/bv3w8AmDZtmsr25s2bo2vXrirboqKicO7cOXTr1g2tW7dWy4OpqWmu+Zs8ebJKPxtTU1O0adNGJX+FLfu0adPyHDwwefLkPK/lffr0gZOTk8o2Dw8PpKWl4fHjxwCAuLg4/O9//4OnpyfMzc1VPjdmZmZo3bq1yvXdwMBAbMlJT08XPwudO3cGAI18FpQke3tLeWvkXQYPHozdu3fj+++/x7Jly/Dhhx+iY8eOGDx4sMotgMLKWXGA/w+0cj6n3B4TE6Oyff/+/ViyZAmuXLmC1NRUlecKei/03r17AICFCxdi4cKF78zrq1evEB8fj4YNG6qlMTIyQq1atRAbG1vgYyoDkNy8fPky3/1oYv/KC4iLi0uhj5eXqKgoxMfH44MPPlB7TiaT4YMPPsChQ4fw5s2b976VoXTz5k3MmTMHp06dQnx8vNoxcqpYsSKqV6+utr1hw4Y4ceIEHj16hObNm+PRo0cAkGveGzduDCCzH1h2tWrV0shIxcIcG8j9M1WlShW1z01BKetRkyZN8kyjrEfKPNevX18tTW7bAM18fvNy48YN9OvXD7Vr18aBAwfUOqmePXsW3377Lf766y+1H3XKPliFUdg6n9d7B2Re9971g0ob9eXRo0ewtLTM9ZZ+/fr1Vb7U35U/ExMT1KpVq9D5K2zZ69atq7atIM/l9z4BmdfUjIwMbN++Hdu3b891Pzl/TK5fvx5r1qzB7du3kZ6ervKcJvv1SDboKQh9fX38/vvvuHr1Ko4dO4bz589j2bJlCAgIwJIlSzBlypQi7f9dXw55PZf919GBAwfQt29ftGjRAkuXLoW9vb34y1HZalMQynSfffYZevXqlWsa5X7zU9AWC+UxN2zYAAcHh1zT5PbFXFDvs/+itLLkR7nv4hqOHR4ejvbt28PMzAwzZsxA/fr1YWpqCh0dHSxcuBCnTp16730q81qY85KzT0Jhveu8vetcFuRz8z6U9ej333/Pc4RZzi+Bgr7Xmvr85iYsLAyenp4wNTXFkSNH1EY3Xr58GV26dIGTkxO+++47ODk5wcTEBDKZDF988QXevn1b6GMXts6/63qY3/unrfpS0DIWxzkpyL7fdbx3fVbf9VxB3idl3R04cGCBWn6XL1+Or776Cl26dMGaNWvEgUMKhQLdu3cv0mchJwY9BdCsWTM0a9YMQGbE2bp1a8yaNQtffPEF9PT0tDa/yJYtW2BkZIQzZ86oVNK3b9+Ktyiyyyuf2aP6/Fq/qlWrhgoVKuDu3btqzyUnJ+PRo0cF6sysPGalSpXyPWZhvM/+lZ0Or127lu9+3/e9Vp6vnCNigMwLxJ07d1CpUqVCD7nfv38/4uPjceDAAbEpWCl7h/XsYmNj8eLFC7WgUjnCSPkrX3lb886dO2jVqpVKWmV5lJ1d8/O+56127dricXr27Kny3K1bt97r2EVRt25dHD16FDY2NuI1IC/K4OfevXvw9PRUeU7ZYpTd+35+Cyo2NhY9evSAXC7HmTNncg36d+zYgfT0dBw5ckQtaIuJiVGbZuN93r/irvO50UZ9cXJyQmhoKGJiYtSueTnfb+WxczsnSUlJePTokViG91VaPis586Sjo4OkpKQCXd+3bNmCmjVr4vjx4yotQLl9bopKsn16CiK30V2VK1eGo6MjUlJSxF9DymZXTQ+ty4+uri5kMplaFDx//vxcI2MzM7Nc89i0aVM4Oztj48aNuVYyQRAQFRUFILNJsnfv3rh79y4OHDigkm758uX53ntXGjhwIIyMjDBv3rxcX5OUlKR2q+Z9vM/+mzRpgkaNGmHbtm24fPmyWtrs5/J932sdHR14eXnh/v372Lt3r8pz27dvx8OHD9G3b99CB87KX105f5kePXoUf//9d56vCwgIUHn8999/Izg4GG5ubuKXkXKUVUBAgMo5iI6OxurVq1GxYkV06dKlQPl83/Pm4eEBMzMzrF69WiUAUCgU+O677wAA/fr1K9C+imLkyJEAgBkzZiAtLU3t+bi4OKSkpAAAunbtChMTE6xZs0al7qakpOQ6zPh9P78FkZKSgt69e+PBgwfYs2dPnoFaXvUmMDAw19vKZmZmePPmTYFaQIq7zudGG/VF2adx0aJFKtuvXLmiNqVD1apV0aFDBxw7dkztc/njjz8iISGh0PkrLZ+V7KpUqYKePXvijz/+wOnTp3NNk72eKetj9novCAK+/fZbjedNsi09N27cwLZt23J9rkuXLrCxscGCBQtw9OhRfPTRR3BycoKOjg5CQkJw/PhxcY4aILMvRIUKFbBmzRqYmJigYsWKqFatmtovb00bMGAA9u7di06dOmH06NEQBAFHjx7FvXv3xOHU2bVu3RobNmzAN998gwYNGkBHRwcff/wxTE1NsXXrVnTu3BnNmjXDmDFj0KhRI7Fj2oEDBzB69Ghxjg/leRk0aBAmTpyI+vXr4+LFizh06BBq1aqldj82N7a2tli3bh3Gjh2LevXqYdSoUXBycsLr169x79497N+/HwcOHCj0kgXvs3/lHCGdO3dGu3btxCHrb9++xaVLl+Do6Che2JSdEKdNm4Zhw4bByMgIjRo1emcfr4CAAAQHB2PIkCE4ffo0GjduLA7frVGjhnhhKowePXrA1NQUI0aMwKeffgpLS0tcuXIFO3bsQOPGjcVfetlZWlrif//7H54/fw4PDw+EhYVh9erVMDExwdKlS8V0nTt3xogRI7B161a4ubmhT58+SEhIwIYNG/Dq1Sts2bKlwB3Xq1Spgtq1a2PXrl2oVasWrKysYGpqio8//jjX9BYWFli+fDm8vb3RokULjB07Fqampti7dy/+/PNPeHt7o3379oU7ae+hefPmWLBgAWbPno1GjRphyJAhsLOzw6tXr3Dz5k0cOnQId+/eRc2aNVGxYkUsXLgQX3zxBT788EOMHj0a+vr62LFjh3hRz/5F/76f34KYO3cuzp49C09PT8TExKhd45ydneHs7Iy+ffti6dKl6NGjByZMmAATExOcO3cOx48fz/Uz3Lp1a/z++++YPHky2rZtC11dXXTu3BnVqlXLNR/FWedzo436MmbMGGzcuBFLlizBkydPxCHra9asQfPmzXH58mWV9/unn35Cx44d0blzZ5Uh6zt27ECTJk1y7exeEKXls5LT2rVr0b59e3h4eGDYsGFo0aIFdHR08PTpU/zxxx9o2bKlOKfagAEDMG3aNHTr1g39+/dHYmIi9u3bp9bPTSM0Ng6sjMhvyDoA4ciRI2LaQYMGCTVr1hSMjY0Fc3NzoUmTJsLixYuFpKQklf3+8ccfgouLi2BoaCgAUBlaijyGrOfcJgjvHhKd22s2btwoNGrUSDAyMhKqVq0qDB06VAgPD891eOvLly+Fvn37CpUqVcp1SGV4eLjw6aefCk5OToKBgYFQsWJFoXHjxsIXX3wh3LlzR2VfoaGhQu/evYUKFSoIZmZmgoeHh3Dt2rUCDaHN7uLFi0L//v0FKysrQV9fX7CyshLatGkjzJ8/X2X45fsOWX/f/QuCIPz777/CqFGjBBsbGzFtt27dhODgYJV0ixYtEhwdHQU9Pb13DsPOLjw8XBg3bpxgY2Mj6OnpCdWrVxe8vb2FFy9eqKV93yHr58+fFzp27CiYm5sLFSpUEDp37iycP39eHJqanfL9efTokeDl5SVYWFgIJiYmgru7u3DlyhW1fSsUCmH58uVC48aNBUNDQ8HMzExwc3MTjh079t75vnTpktC2bVvBxMREbWqDvF57+PBhoVOnToKZmZlgaGgoODs7CytXrlQbaptbWQuaL6V3Das/evSo4OnpKVSpUkXQ19cXqlevLri5uQk//vij2rUgKChI+OCDDwQDAwOhevXqwtSpU4VLly4JAIRFixappH2fz29Bhqwrz0Nef9nLdujQIaFFixaCiYmJUKlSJeHjjz8W7ty5k+tnOCEhQRg7dqxQrVo1QUdHR+WzmNc1q6B1/l1TfrzvFBElWV8EIXOagkmTJgnVqlUTDA0NhWbNmgn79u0TfH19BQDCy5cvVdLfvXtXGDhwoGBpaSno6+sLDg4Ogq+vrxAbG6uS7l3Xu7zyromyv+u47/p85PW6169fC9OnTxfq168vGBoaChUqVBDq168veHt7CxcvXhTTKRQKYdGiRUKdOnUEQ0NDoXr16sKkSZOE169fv/P6XhgyQSiBhVuIiCTs119/xcCBA7Fr1y4MGjRI29mhYtazZ0+cOXMGcrm8SFNekObx3SAi0pDk5GS1fi8pKSlYsmQJ9PX1xSVrqHzIbemFy5cv4+jRo3B3d2fAUwpJtk8PEZGmnT9/Hp9//jn69+8Pe3t7vHjxAjt37sT9+/cxZ86cPPvAUNk0adIkyOVytGvXThyxtmHDBhgbG2P+/Pnazh7lgre3iIg05NGjR/Dz88Pff/+NqKgo6Onp4YMPPsDEiRMxZswYbWePNGz79u1YvXo1QkNDIZfLUblyZXTs2BFz587VyAS2pHkMeoiIiEgSeMORiIiIJIFBDxEREUkCgx4iIiKSBAY9REREJAkMeoiIiEgSGPQQERGRJHBywhxyW1m9qGQyGe7p6SE9Pb1AqxRLna6uLhQKhbazUarJZDLoSaxONS7Ca1mnCkaK9aqwWKcKpiTrVEEW6mVLTwnR19fXdhbKDOWK1PRurFMFxzpVcKxXBcM6VXClqU4x6CEiIiJJYNBDREREksCgh4iIiCSBQQ8RERFJAoMeIiIikgQGPURERCQJDHqIiIhIEhj0EBERkSQw6CEiIiJJYNBDREREksCgh4iIiCSBQQ8RERFJAoMeIiIikgQGPURERCQJetrOABEREZUvcrkccrkcMpkMRkZGSE5OhiAIMDc3h7m5udbyxZYeIiIi0qjAwEC4uLigadOmqF+/Ppo2bQoXFxcEBgZqNV9s6SEiIiKNmjhxIoYOHYrIyEj06NEDR44cgbW1tVZbeQAGPURERKRhyttYMpkMAGBjYwNbW1st54q3t4iIiEgiGPQQERGRJDDoISIiIklg0ENERESSwKCHiIiIJIFBDxEREUkCh6znoKurC11dXY3uUyaTQVdXFwYGBhAEQaP7Lo/09Fgt8yPFOmVQhNeyThWMFOtVYbFOFYy+vr74r4FBUT7FmsF3LQeFQgGFQqHRfcpkMih0dZGamsoLSQGlpqZqOwulmvLLSUp1qqg1gnUqf1KsV0XBOpW/tLQ08d/ScL54e4uIiIgkgUEPERERSQJvbxEREZVTytXOc9L2aufawpYeIiKickq52nnOP22vdq4tbOkhIiIqp5SrnUdERMDT0xOHDx+GjY2NJFt5AAY9RERE5VbO21g2Njaws7PTYo60i7e3iIiISBIY9BAREZEkMOghIiIiSWDQQ0RERJLAoIeIiIgkgUEPERERSQKDHiIiIpIEBj1EREQkCZyckIiIqJhxDazSgS09RERExYxrYJUObOkhIiIqZlwDq3Rg0ENERFTMuAZW6cDbW0RERCQJDHqIiIhIEhj0EBERkSQw6CEiIiJJYNBDREREksCgh4iIiCSBQQ8RERFJAoMeIiIikgQGPURERCQJDHqIiIhIEhj0EBERkSQw6CEiIiJJYNBDREREksBV1omISFLkcjnkcrnKtpyroFP5xJYeIiKSlMDAQLi4uKj8BQYGajtbVALY0kNERJIyceJEuLm5wdPTE4cPH4aNjQ1beSSiTAU9CQkJWL16Na5evQpjY2MMHDgQnp6euabt1asXDA0NIZPJAAANGzbEvHnzSjC3RERUGpmbm8PGxgYAYGNjAzs7Oy3niEpKmQp61q1bB4VCgU2bNiEiIgJz5syBnZ0dnJ2dc02/bNkyVmYiIiICUIb69CQnJ+PPP//E8OHDYWJiglq1aqFz584IDg7WdtaIiIioDCgzLT3Pnz8HANjb24vbnJyccODAgTxfM3v2bCgUCtSpUwejR49WeS0RERFJS5kJepKTk2FsbKyyzdTUFElJSbmmDwgIQL169ZCWloZ9+/Zhzpw5WLNmDUxMTFTSRUREICIiQnxsaGgIa2trjeZd2a9I+S+9m0wm47nKhxTrVFFKyjpVMFKqV9nLWpjyFuV1RXl9YUntuHkpM0GPkZGRWoDz9u1btUBIqVGjRgAAfX19DB8+HKdPn8a9e/fQvHlzlXTr1q2Dv7+/+HjatGmYM2eOhnOfmQ8qGF1dXejq6mo7G6We1OpU7p/0gmGdKjip1CsjIyPx37y+R96lsHWqqMctLG0d19DQUPy3JI+blzIT9Nja2gIAwsPDUaNGDQDA48eP4eDgUKDX5xVh+vj4oFevXuJjQ0PDPFuPCksmkyENma1VgiBodN/lkYGBAVJTU7WdjVJNWZ+lVKeK8qlknSoYKdWr5ORk8d/CXPMLW6eKetzC0tZxU1JSxH+L+7g57+TkpswEPUZGRmjXrh22b9+Ozz//HC9fvsTJkyfx9ddfq6UNCwtDWloaatasifT0dPz2229ITU1FvXr11NLa2NiIQxcBIDo6utg+7IIglPsLiSbwPBWclM5VUUoppfOkCVI4X8ryFbasRXldUV5fWFI7bl7KTNADZLbKrFq1CqNHj4aJiQmGDRuGJk2aAAAGDhyIuXPn4oMPPkBsbCzWrl2L6OhoGBgYoHbt2vD394eZmZmWS0BERETaUqaCHjMzM0yfPj3X5/bs2SP+39nZGWvXri2pbBEREVEZUKaCHiIi0jzlApwymQxGRkZinx4uwknlTZmZnJCIiIqHcgHOpk2bon79+mjatCkX4aRyiS09REQSN3HiRAwdOhSRkZHo0aMHjhw5Amtra7byULnDoIeISOKUt7GUQ9ZtbGzEaUKIyhPe3iIiIiJJYNBDREREksCgh4iIiCSBQQ8RERFJAoMeIiIikgQGPURERCQJDHqIiIhIEhj0EBERkSQw6CEiIiJJYNBDREREksCgh4iIiCSBQQ8RERFJAoMeIiIikgQGPURERCQJDHqIiIhIEhj0EBERkSQw6CEiIiJJYNBDREREksCgh4iIiCSBQQ8RERFJAoMeIiIikgQGPURERCQJDHqIiIhIEhj0EBERkSToaTsDRESUSS6XQy6Xq203NzeHubm5FnJEVL4w6CEiKiUCAwOxZMkSte1+fn74+uuvtZCj4sUgj0oagx4iolJi4sSJGDp0KCIiIuDp6YnDhw/Dxsam3AYAUgvySPsY9BARlRI5WzhsbGxgZ2enxRwVL6kFeaR9DHqIiEgrpBbkkfZx9BYRERFJAoMeIiIikgTe3spBV1cXurq6Gt2nTCaDrq4uDAwMIAiCRvddHunpsVrmR4p1yqAIry1rdUpfX1/818CgKCXncYvrmIWtU1I6x9o8bl7K1pWgBCgUCigUCo3uUyaTQaGri9TUVMl8QRVVamqqtrNQqimDHinVqaLWiLJUp9LS0sR/SzLfUjquJo5ZmNdJ6Rxr87h54e0tIiIikgQGPURERCQJDHqIiIhIEhj0EBERkSQw6CEiIiJJYNBDREREksCgh4iIiCSBQQ8RERFJAoMeIiIikgQGPURERCQJDHqIiIhIEhj0EBERkSQw6CEiIiJJYNBDREREksCgh4iIiCSBQQ8RERFJAoMeIiIikgQGPURERCQJDHqIiIhIEhj0EBERkSQw6CEiIiJJYNBDREREksCgh4iIiCSBQQ8RERFJAoMeIiIikgQGPURERCQJDHqIiIjKgYMHD2LkyJFo0qQJHBwc0KlTJ2zatAkZGRkq6U6cOAE3NzfY2dmhZcuW+OWXX3Ld3+rVq9GsWTPY2dnBw8MDf/75p1qahIQETJkyBXXr1oWDgwOGDx+O8PDwYimfJjDoISIiKgfWrl0LAwMDzJs3D9u3b0ePHj0wc+ZMfPvtt2KaGzduYOTIkXB2dsauXbswePBgzJgxA1u3blXZ1+rVq/Hdd99h3Lhx2LlzJxwdHTF48GDcvXtXJd2ECRNw7NgxfP/991i/fj0iIiLQr18/JCUllUiZ35eetjNARERERbdt2zZYWlqKj9u3b4+3b99i48aNGDlyJADg559/hrOzM1asWCGmef78ORYtWoRhw4ZBR0cHKSkpWLp0KXx8fPDpp58CANq2bYuOHTti2bJlWL9+PQDgypUrOHHiBHbs2AEPDw8AQMOGDdGyZUvs3r0bo0ePLsHSFwxbeoiItODRo0eYOnUqXF1dYW1tjQ4dOqilEQQBK1euRPPmzWFra4sOHTpg//79aunkcjl8fX1Rr1492Nvbo3fv3rh165ZKmrCwMFStWlXtr3v37sVWxuxKuryLFy/OtbxVq1bF1KlTi62c2pQ94FFq3LgxkpOTERcXBwD4559/0KdPH5U0/fr1w8uXL8Vz+M8//0Aul6uk09XVhZeXF4KDgyEIAgAgODgYFhYWcHd3F9PZ2dmhVatWOHHiBKKjo/H7778DAJ49e6bZwhYSW3qIiHKQy+WQy+Vq283NzWFubq6RY4SGhuLEiRNo1qwZBEFQ63cBAJs3b8bq1avh6+uLli1b4siRI/Dx8YGJiQm6desmpvPx8cH169cxZ84cVK1aFevWrUPfvn0REhICW1tblX3OmjUL7du3Fx+bmZlppDz5KenyDh8+HJ07d1bZ/4ULF/Dtt9+iS5cuxVtYZAZ5a9asweXLl3H//n3UqVMH586dU0kjCAJWrVqFoKAgREZGwsnJCb6+vmpBiVwux7x58/DHH38gKSkJLi4uWLBgARo3bqyS7sGDB5gzZw4uXLgAfX19dO3aFQBQqVIlVK5cGQCQlpaGOnXqqLyuXr164uubNGmCBw8eAIBaurp16yIhIQERERGoXr06Hjx4gNq1a0Mmk6mlO378ONq2bYuEhAQAQN++fbFt2za4ubm997nUJAY9REQ5BAYGYsmSJWrb/fz88PXXX2vkGN26dUOPHj0AAJMnT8aNGzfU0qxfvx7e3t7w8/MDALi6uiI8PBwLFy4Ug4DLly8jODgY27ZtE7e1b98eLVq0wOrVqxEQEKCyTycnJ7Ro0UIjZXgfJV3e6tWro3r16ir7DwoKQsWKFdGlSxe8evWq2MoKFCzIW7VqFQICAjQS5MXHx6NPnz6wsbFBYGAgkpKSMGfOHLx8+RLTpk2Drq6uuD8LCwuVfFSsWBEA8ObNGwBAbGwsDA0NYWxsnGe66tWrIy4uLtcfARUrVkRkZCQAQKFQAABSU1Mxbtw43L9/HwYGBu9zKjWKt7eIiHKYOHEirl27hsOHDwMADh8+jGvXrmHixIkaO4aOTv6X38TERLVfxp07d8adO3fE2wW3bt2CTCaDq6urmMbExAStWrXC8ePH3ytPyl/lxUHb5U1OTsbhw4fx8ccfl8iXbrdu3XDjxg1s2rRJrUUGyGxxWbp0qRjkubq6YtGiRXB3d8fChQvFdMogb/ny5Rg2bBi6du2KrVu3Ql9fH6tXrxbT/fLLL4iPj8fOnTvRo0cPtGvXTrwNlbPFJmfLTG7bc0uj3F9B0ikUCjHgUYqPj0dERESuxy4pDHqIiHIwNzeHnZ0dbGxsAAA2Njaws7PT2K2t95HzC9rQ0BAAxFsQKSkp0NHRUfklr0wXFhamNorm66+/hpWVFRo0aICvvvoKb968wb///os2bdqIt72WLVuWa8tESdB0eZWOHz+O+Ph49OvXD8eOHRODJjc3t/cODgsivyAvPDwcCQkJGgvybt26hUaNGqFq1aqQy+UYPHgwLCwsULFiRZw8eVIMWIDMlpzslI+VLTkVK1ZEcnIykpOTVdIp+wUp01lYWIjbcqbL+f4AmQFSlSpV8j4pJYBBDxFRKSWTyXD16lWVbZcvXwbw/19UTk5OUCgUuHnzppgmIyMD165dgyAI4peSgYEBxowZg2XLlmHfvn345JNPcPDgQXh5eaFPnz54/Pix+Ppt27aptCKUFE2WN6d9+/bBxsYGJiYmGDVqlJguNjYWI0eOVNlfSUhJSQGguSAvJSUF+vr6SE5OxvDhwxEVFYXdu3fDyMgIISEhaN26tfjaS5cuqewrNDQUQGZfnOz/KvOg9ODBA5iZmYk/BurWrYv//vtPJaBSpmvYsKFKK5COjg6mTJlSYn3I8lKug56EhAQsWrQIgwYNwujRo8WmaiKisqBnz55YuXIlgoODERsbi927d4ujmZRfKG5ubnBycsLUqVNx9+5dREVFYe7cuXj69CmA/29xsLa2xuLFi8VbH5999hkCAwNx9+5dvHz5UuVWhEKhwK+//lrCpdVsebOTy+UIDg5Gnz59cPjwYbU0urq6Jf79YG9vDx0dHY0FebVq1cLdu3cxZswY3LlzB7t374ZMJsPLly8RERGB1NRU8fWrV69WaRHbv38/rKysxNtwLVu2hLm5OQ4cOCCmUSgUOHjwINzd3cX3wt3dHXFxcTh16pSY7vnz57h06RJGjhyJnTt3ikPZ586dq7H+cEVRroOedevWQaFQYNOmTfjmm2+wffv2Eo/miYgKa+rUqWjatCmGDBmCOnXqYO7cuZg+fToAoFq1agAAfX19bNiwAYmJiejUqRMaNmyIs2fPwsfHB/r6+qhUqVKe+/fw8ICRkVGJlKUgiqu8hw4dQkpKCvr165fnbbuSvp1namqKAQMGaCzIGzFihBjcjRs3Dg8fPsSoUaPUWmGAzE7F48ePx59//omlS5di69atmDZtmrgvQ0ND+Pr6IjAwEGvWrMH58+fx6aef4unTp/jqq6/E/TRv3hweHh748ssvsX//fpw4cQKjRo1CjRo1MGjQIHTp0kUcENC7d+88+xKVpHI7eis5ORl//vknli9fDhMTE9SqVQudO3dGcHAwnJ2dtZ09IqJ8WVhYYM+ePYiMjMSbN2/g5OSEo0ePwsDAQOU61rhxY1y4cAGPHj0CkNk6MG3aNDRp0gT6+vrvPIaOjg7MzMyQlJQktvbo6upiyJAhxVewPBRXefft24c6derA2dkZaWlpWLVqlcrz6enp4siykjR//ny8evVKPNdVqlTB9OnTMXfuXLUgz9vbG506dQKQOQGgj48P1q9fLwZ5tWrVQsWKFfH69WssW7Ys32P/+++/GDhwIKpXr46AgACMGDFC5flPPvkEgiBg/fr1iIqKQoMGDbBz5040bNhQJd26devEVpzU1FR06NABmzZtUhv5VVqU26Dn+fPnADKbEJWcnJxUmusAICIiQqU3eWpqqsbfLJlMhkhDQ6SkpOQadZMqfX19pKWlaTsbpZpMJoOhxOrU8yK8trB1SjnsNjIyslh/pSYmJiItLU28buV2XHNzc7x8+RKBgYHw8PDIdS4hZavN7du3sW/fPnz55ZfiPnMTEhKCxMREfP3119i3bx/+++8/AED//v3x0UcfvfO1RVGS5Y2KisL58+cxceJEPH/+HNbW1vj+++/xzTffICUlBUZGRliwYAGqVav2XuV9nzqVvbzZy2pjY4Nly5bh1atXkMvlsLe3R0hICPT19WFpaSnmp3Llyti7dy/CwsIgCAIcHBywcOFC1K9fX2Xo/alTp5CamoqwsDCYm5ujWrVqcHNzE4eiA5nXDktLS2zfvh0mJibi9tzK3qdPH7U5g3JLN2XKFEyZMiXXdNnLW9xym5xRjVBO3b59Wxg6dKjKtkuXLgne3t4q2+bOnSsAEP/atWun8ph//OMf//jHP/6V/r+CKLctPUZGRmpDF9++favWiuPj44NevXqJj4urpeexBH6VJ8bHIzGXeT5MzMxgUqFCgffDlp78SbGlp34RXlvYOhUREYEePXrgyJEj4ogVTXn+/Dl69uyZ63Pff/89pk+fjqlTp2Lv3r148eIFjI2N0b59e3z++eewsrJSSf/jjz/i2LFjePPmDSwtLdGzZ094e3uLI4GAzM6qu3fvRnh4OJKTk8VWgE8++UQcURMZGYnu3bvj6NGjsLa2LtPlVfr4448RERGBTp06ISwsDBkZGfjtt99U3tujR49i1apV8Pb2RpMmTRASEoLdu3dj+fLlKkPFJ0+ejDt37oh5CgoKwv3797Fnzx7xfEVHR6N///6oXr06WrZsCX19ffz666+Ii4vDZ599hgoVKiAgIABHjhzB1atXkZKSgho1aiAqKgp79+7F8+fPsW3bNlStWlU87vr161GjRg1UqVIFT548wcaNG+Ho6IjVq1eL/XASExMRGBiIZs2awdDQEDdv3sQvv/yC8ePHw9vbG2FhYejVqxcOHz6sNmljcSrOOlUYMqGcXjGTk5MxdOhQrFixAjVq1ACQOXlTbGwsfH1983xddHS0xvMik8nwn7ExkpKSyvUXVNDixdiSyyy2I/38MPo9eu0bGBiojDQgdTKZDMYSqFPZqU/vVnCFrVPPnj2Di4sLrl27Bjs7uyLkoGwc9/nz52jatCmuX7+utnxFcSrO8mZkZIiBgXIm6HPnzonH/Pvvv9G5c2eMGDFCZTXyoUOH4sWLFwgJCQGQOaqqR48e4kzQBgYGiI2NRYsWLeDl5SXOBP3rr7/ik08+waFDh1R+UOd07do1XLhwAcuWLUNYWBhMTU3h7u6OWbNmqQUlc+fOxYEDBxAdHQ0rKyv0798fvr6+Kp3Qk5KSMGrUKNy4cQNv375F7dq14ePjI/YXkkKdKsjtrXLd0tOuXTts374dn3/+OV6+fImTJ0+WiiFz5VX/iRPRY+hQREdE4DNPT6w8fBiWNjYw1cKEbkREQNEmCZwxYwaePXsGOzu7fCcJVAY9yhbF+vXrIyoqSkzbqFEjcYj377//jp9++glz587FgAED8i2Dv78//P3935nG2NgYe/bsyXdfUleuh6z7+PgAAEaPHg1/f38MGzYMTZo00XKuyi8zc3NY2dnBMus2gKWNDazs7GDGoIeISilNTxLYo0cPVK1aFbNnz0ZkZCRev36NxYsXIz4+Hr/99ps4H9DWrVvh5eXFVu0SVm5beoDM1YOVczwQERHllH2SwHbt2onb3zVJYLNmzQCoTxJobGyMSpUq4X//+x+GDRsmTvZnbm6ONm3aICQkRJwPKD09Hbdv30ZwcDA8PT1LsMTSVq5beoiIiN5F05MERkVFiRP07dixA7/++is8PDxw5swZtQU49fT0VG6BUfEr1y09RERE+dHkJIGrVq1CbGwsgoODxY7GnTp1wrlz5xAVFaUy8CAlJQWNGjUqyaJKHlt6iIhI0ipVqoQ9e/bg1q1bOHv2LG7cuAFbW9s8Z4K+ePEirly5gpCQECQlJanMBB0aGoo6deqojKySyWTo1KkTjIyMVDpWT506Fc2bNy+5ghKDHiIiIiBzUdYGDRpAT08PQUFB8PLyQoUcc4zJZDLUqlULtWvXRkxMDA4ePIjhw4eLz9eoUQMPHjxQmScuIyMDd+/eRdu2bcUlIrZt28bRxFrA21tERFRuJSYmIjg4GEDmXDXx8fE4dOgQYmJixDR79+5FUlISnJycEBkZic2bNyMsLAxr165V2dfSpUvh6OiIatWq4fHjx/jxxx/h7Oyssk7ZiBEjsG3bNgwfPhze3t7Q09PDjh07cOfOHXzzzTeoV68eAOCDDz4ogdJTTgx6iIio3IqOjsa4ceNUtuV8LAgC1q5dqzJJYGBgoNoMwnFxcZg3bx6io6NhbW2NAQMGwNfXV+WWlbOzM/bu3YslS5bgiy++QHp6OurWrYstW7agS5cuePbsWfEVlvLFoIeIiMote3v7XEdIKWcoBoABAwa89ySB75rlu127dirD36n0YJ8eIpK8R48eYerUqXB1dYW1tTU6dOiglkYQBKxcuRLNmzeHra0tOnToIA5rzk4ul8PX1xf16tWDvb09evfujVu3br3z+CNGjEDVqlWxevVqjZWJiNQx6CEiyQsNDcWJEyfg6Ogo9rnIadWqVQgICMDgwYOxfft2tG3bFj4+Pjh27JhKOh8fHxw5cgRz5szBhg0boKenh759++L58+e57jc4OBhXrlzReJnyoq0A7+LFi+jTpw8cHBzg5OSEnj174uHDhxovH9G7MOghIsnr1q0bbty4gU2bNomz6GaXlpaGpUuXwtvbG35+fnB1dcWiRYvg7u6OhQsXiukuX76M4OBgLF++HMOGDUPXrl2xdetW6Ovr59qKk5KSgpkzZ2L27NnFWr7sChLgBQUFaTTACwkJQb9+/cS+LT///DM6duyoMsKJqCQw6CEiySvKopR37twRO6fmtyhlTqtXr4aFhYXK6J/Xr1/j/PnzAFAsQUF+AR4ArF+/XmMBXnp6Or766itMmjQJixYtQqdOneDu7o5p06ZxYj4qcUUOel69eoWwsDC1PyKi8kLTi1ICmR1pV6xYgYCAAHGpg/DwcLRs2RJffPEFgMwOtpoe7ZNfgAdkDvPWVIAXEhKCZ8+eYfz48bke69q1azhw4AAAqC3TQKRphQp6YmJiMHToUBgZGcHGxgaOjo7iX82aNeHo6KjpfBIRaU32RSmze9eilEo5F6VUmj17Nj766CO0bNlS3LZz507Ex8eLi1JGRkZi0qRJxVWsd9JUgHflyhVUrlwZV69eRatWrWBtbY22bdti//79WLRoEbp16ya2IE2YMEEMMImKQ6GGrI8fPx4hISHw8/NDw4YN1T4cRETlSfZFKRs0aIAWLVrg2LFj71yUctWqVahatSp++ukntUUpT58+jZCQEFy8eFHlOImJiSqPFQoFrl27VtzFUyOTyTS26virV6+QmJiIL7/8EjNmzECtWrWwa9cuTJgwQdy3cuj3zZs3sWrVKkyZMqWESkpSU6ig5/Tp0/jpp58wcuRITeeHiEiUIJfjrVwOAKiUbbu5uTnMzc1LNC+aXJRy5syZ8Pb2hrGxsUrrT27MzMyKsVS5++ijjzQW4GVkZCA5ORn+/v4YM2YMAKBDhw44f/48Xr58KbZqAZn9f3K2phFpUqFub1WsWBGWlpaazgsRkYq9gYEY4uKCIS4ucMn2FxgYWOJ5KcyilBcvXsx1Ucr//vsPy5cvR+3atcW/3MhkMkybNq1Eypfd1KlT0bRpUwwZMgR16tTB3LlzMX36dABQC/ASExPRqVMnNGzYEGfPnoWPjw/09fXFAE/5b/ah8TKZDB988IFKwANkBko5Z0Em0qRCtfT4+flh5cqV6Nq1K/T0OKlzaZP913F2pubmMCvhX8dERdF/4kT0GDoU0RER+MzTE4cPH4aNjU2Jt/JkZ21tDWtraygUinwXpQQyl0E4ePAg5syZIz6v7LibnZeXF0aOHIn09HRcvnwZDx48wLfffiu2jpQkCwsL7NmzB5GRkXjz5g2cnJxw9OjRPAO8R48eAci85TVt2jSVAK9OnTq5HsPe3l5Mk5aWBiCzP9Cnn35anEUjiStUxHL//n3cvXsXtWrVQqdOnVCxYkWV52UyGVasWKGJ/FEh7A0MxJYlS9S2j/Tzw2iu6ktliFmOQN3GxgZ2dnYaP05xLUqpbNHJuShlXksUODk54dNPPxWXSPjoo480Xtb3oYkAz83NDXp6ejhz5owYAAmCgIsXL8LV1RX16tXDxYsXcfnyZezcuRNOTk4lV0CSnEIFPb///rt4v/bcuXNqzzPo0a6cv45XHj4MSxsbmLKVhyhXxbUopZWVVa6LUmpTXgEeAHHk7R9//AFjY2ONBHjW1tYYO3YsFixYAEEQxI7M9+/fx08//YQmTZqIQZ6Dg0MJnQWSqkIFPY8fP9Z0PkiDcv46trSxgVUx/DomKi+Ka1HK95Hb8YvDuwK89evXA9B8gOfv7w8zMzP89NNPeP36NerVq4ft27ejSZMmxVhSInXskENEJCF5BXgAxOUjPvroI/j4+OS7r4IGeHp6epgxYwZmzJjxfpkl0rAiBT3//fcfHjx4gOTkZLXn+vbtW5RdExEREWlUoYIeuVyOvn374vTp0wAym0KB/5+/AeB04kRERFS6FKpn3bRp0xAREYFz585BEATs378fISEhGDduHBwdHdVmGSUiIiLStkIFPUePHsWsWbPQqlUrAED16tXRsWNH/Pzzz/Dy8sKPP/6o0UwSERERFVWhgp5Xr16hRo0a0NXVhampqcpcFj169MDRo0c1lkEiIiIiTShU0FOjRg1ER0cDyJxtUznHAwD89ddfMDIy0kzuiIiIiDSkUB2ZPTw8EBwcjD59+uCrr77CqFGjcOnSJRgYGODvv//mCrlERERU6hQq6Fm0aBESExMBACNGjICZmZk4RfuqVasKNL8DkaZwrTEiIiqIQgU9JiYmMDExER/36dMHffr00VimiN4H1xojIqKCKNJiMPfu3cPWrVsREBCAyMhIAJkTFsbHx2skc2XNo0ePMHXqVLi6usLa2hodOnRQSyMIAnatXImhzZujm60txnbogNP796uli4uJwbKpUzHYxQU97O0xpl07/LZunTgnUna3Ll6Eb58+8HRwwMdOTvi8Z0+EP3xYLGUsjfpPnIid165h5eHDAICVhw9j57Vr6D9xopZzRkREpUmhWnoSExMxfvx47NmzB0DmF3n37t1hbW2NGTNmwNHREYsXL9ZoRsuC0NBQnDhxAs2aNYMgCMjIyFBLs3vVKmwMCMBwX1980LIl/jxyBAt8fGBoYoK23bqJ6eaMHo3njx5h3KxZsKpRA1dCQrB69mxkZGRgwKRJYrrLISGYNWwYPIcPx7Avv0R6WhruXbmC1KQkGBgalki5tY1rjRERUUEUKuiZOnUqTp06hd9//x0dOnRAhQoVxOc8PT2xbNkySQY93bp1Q48ePQAAkydPxo0bN1SeT0tNxbalS9HX2xuj/PwAAC1cXfEyPBy/LFwoBj1RL17g1sWL8FuxAj2GDgUANOvQAQ/v3MHpAwfEoEeRno4fvvoKAyZNwvjZs8XjtHJ3BwC8fPaseAtMRERUhhTq9tbevXuxaNEidO/eXW14es2aNfHkyRNN5K3MybmycE4vnjxBYkICWri5qWxv2bkzHt25IwYp6WlpADI74mZnZmGhcnvrckgIXj17Bq/x49WOpVAo8CLrfchIT3/vshAREZU3hQp6EhISYGNjk+tzb9++LVKGyrPUrIVZ9Q0MVLbrZ92GCnvwAABg4+CA5q6u2L5sGR7fu4fEhASc++MPnD98GH3GjRNfd+/KFZhXroz7V69iZKtWcLe2xui2bXEoKAgT3NwwJatzud+AAWz1ISIiySvU7S1nZ2f89ttv6Nq1q9pzf/zxB1q0aFHkjGmLrq4udHV1NbIfmUwGAwMDyGQy6Orqoma9etDR0cG/N27gw2ytPaFXrwIAEhMSYJAVEH2/YwfmjhmDcR07AshczHXSt9/ioxEjxNfFRkcjOTERP3z5Jbxnz0aN2rVxZMcOLPfzUylDZHg45o4ejV/Oni1yuQpCX19f/NcgR4BXEHp6haqWRT5uWaKsUwYGBrl2bi9vpFaneNzSf0zWqdJ93LwU6l375ptv0Lt3byQmJmLAgAGQyWT4+++/sXPnTvzyyy84nDWKpixSKBQaWSFeoVBAEASkpqZCJpNBoasLfUNDeAwYgO3Ll8O+bl00bNECF44dQ/BvvwEAMhQKpKamQhAELJg4EWH//otZgYGwtLbG9b/+wvr582FiZgbP4cMBAOnp6UhNTsYkf3/0HDkSAFC/eXMc271bpQwZCgUe3LiB2NevYWJmVuSy5Sct6/ZcWloaUlNTC7WPwrxOE8ctK5RBj7K+lHdSq1M8btk4JutU6T1uXgoV9PTs2RO7du2Cn58ftm/fDgD45JNPYGdnh+3bt6NLly4azWR5Mmn+fLx+9QozhgwBAFhUqYIx06cjcO5cVK5WDQBw8cQJnDl0COtDQlDrgw8AAE3atUNCXBwC581D96FDoaOjA/NKlQAALtmGxuvq6UEmk6l/Ecpk0MuKuImIiKSo0PP09O/fH48fP8b9+/dx/vx53L17F2FhYejfv78m81fumFeqhEV79mDPrVvYePYsdt+4gWq2ttA3MEBtZ2cAwNPQUOjo6sKpYUOV19Zu1AgJcXGQv34NALCvU0dt/7q6unCsX191m54eug0aJJkh7FR2PX/0CMumToW3qyvcra0xNo+5rlauXInmzZvD1tYWHTp0wP5c5rqSy+Xw9fVFvXr1YG1tjd69e+PWrVsqacLDwzFs2DA4OzvDzs4OjRo1wtixY/FQQvNcEUlJ4W5KIvNiceDAAYSHhyM5q4Oukkwmw4oVK4qcufLM0toaltbWUCgUOBQUBFcvL5hmDf23srNDhkKB/27dQp2sQAgAQq9fh7GpKSyqVAEAtHRzg66eHq6cOSMGQIIgQABQ3dERb+PiEPf6NdwHDMBXucxYTFRQJbXUx5PQUFw8cQIN3jHX1apVqxAQEABfX1+0bNkSR44cgY+PD0xMTNAt21xXPj4+uH79OubMmYPq1atj1apV6Nu3L0JCQmBrawsgc+CFlZUV+vbtC2tra7x8+RIrVqyAl5cXQkJCNFYuIiodChX07NmzByNGjEBGRgaqVaum1jlJqkFPYmIigoODAQDPnj1DfHw8Dh06BJlMBpvOnWFkZobgvXuRkpQEWycnxERG4n+bNyMyLAyz1q4V99PawwPW9vbwHzcOI/38YGltjWvnzuFQUBAGfvIJZDIZAKCKtTV6jx2LDQsWQBAE1KhVC0d37cLT0FCsOX4cFlWqYIiLC0Z//TVbecoJZfAhk8lgZGSE5ORkCIJQ7OuMldRSH226dUO7rLmuFk2ejNCcc12lpWHp0qXw9vaGX9ZcV66urggPD8fChQvFoOfy5csIDg7Gtm3b0K1bNxgYGKBVq1Zo0aIFVq9ejYCAAABA/fr1sXTpUpVjNGnSBK1bt0ZISAhatWqlsbIRkfYVKuiZOXMmvLy88PPPP8PCwkLTeSqzoqOjMS7bkHIA4uPVR46gYcuWEAQBv65di4iwMBibmqKVuztmBQaiirW1+BpjMzP88Ntv2BgQgA0LFiA+NhbW9vYYP3s2+udYzHWSvz9MzMyw86efIH/9Gg716uG77dtRt0kTDlMvh7S1zlj/iRPRY+hQREdE4DNPT6w8fBiWNjZqc0kVVX5zXYWHhyMhIQFuOea66ty5M2bMmIFnz57Bzs4Ot27dgkwmg6urq5jGxMQErVq1wvHjx8WgJzeVK1cGALx48QJbt24FAPz777+w4yzfRGVeoYKeqKgoTJgwgQFPDvb29oiKilLbLpPJ8J+xMZKSkuAxYAA8BgzId1/Va9bENz//nG86XT09jJ0xA2NnzChUnqlsUQYfMZGRmNyjB1YdOYIq1tYaDz5yKi1LfaSkpACAWuuyYVZL5oMHD2BnZ4eUlBTo6OioTT9haGiIsLAwJCUlwdjYWNyekZEBhUKBiIgIfPfdd7C0tMSiRYvEAQGDBw/G5s2b0b179+IsHhEVs0J1ZO7RowcuXryo6bwQUT7MzM1hZWcHy6zJQZXBR3He2ipN7O3toaOjg6tZc1spXb58GQAQGxsLAHBycoJCocDNmzfFNBkZGbh27RoEQUBcXJzK6z/99FNUr14dzZs3x9WrV6Gnp4fU1FRxiG1GRgYmTpyokeksiEh7CtXSs3btWgwePBiJiYno0qULKlasqJamWbNmRc0bFdDzR4+wZ80a3L18GY/v34d9nTr45dw5lTSCIGD3qlU4FBSEmMhI2Do5YYSvL9yyZm1WiouJwS8LF+LSyZOIi4mBdY0a+GjkSPSdMEHsSxS0eHGut1gA4ONRozD0yy+LpZxEpqamGDBgAFauXIkGDRqgRYsWOHbsmDh6S1lH3dzc4OTkhKlTp2LVqlWwtbXFDz/8gKdPnwJQv402ffp0TJgwAc+ePUNgYCD+/vtvtWO/ffsWUVFRsM52K5qIypZCBT1yuRwJCQlYuHAhvv/+e5XnBEHInIyPv4hKTEFGvGhydfeew4fjw86dVfZ/88IF/Pztt/iQczRRMZs/fz5evXqFIVlzXVWpUgXTp0/H3LlzUS1rrit9fX1s2LAB3t7e6NSpEwCgYcOG8PHxwfr161Epa44rJQcHBzg4OMDFxQWdO3eGo6Oj2lxX+vr6qJI1cpKIyqZCBT0jRoxAeHg4Vq5cibp165aKqaWlLL8RL+lpaRpd3b1q9eqoWr26yjH+FxSEChUr4sMuXfD61atiLW9pUVLDuElVpUqVsGfPHkRGRuLNmzdwcnLC0aNHYWBgAOdsUzw0btwYFy5cwKNHj6Cvr48aNWpg2rRpaNKkiTg1fm5MTU1Rt25dhIaGQkdHR/wRsXDhwne+johKv0IFPZcvX8aOHTvg5eWl4exQYeQ34uVleHieq7uvnDEDL589g5Wd3TtXd4/LmhAxN6nJyTh/+DBcvbxw7o8/EHLwIADgv1u3tNLZtaRoayQVZbK2toZ11lxXQUFB8PLyQoWsua6UZDIZatWqBQMDA7x48QIHDx7EnDlz3rnfuLg4REVFwdPTE2ZmZtizZw9WrVqFQYMGFWdxiKgEFCroqV27Nm9flSFpWSNe3rW6u5Wdncrq7jVq1cq8vXXmDM4fPgzfH37Ic/8Xjh/H2/h4pKWmImDiRPGX8ZxRo/D97t1omSPYKi9Kahi3lCQnJuJS1lxXL589Q2J8PM4cOoS4mBgxzd69e5GUlAQnJydERkZi8+bNCAsLw9psc10BwNKlS+Ho6Ihq1arh8ePH+PHHH+Hs7CzeFgOAxYsXQy6X48MPP0SVKlUQHh6OdevWIS0tDXPmzIGhoSH27NmDdu3alcwJIKJiVaig54cffsD06dPRuHFj1K1bV9N5Ig2zyhrxcv/qVTTNdvG+lzXiRZ414gUAvg0Kwnxvb5XV3SfMnYuu7/iVe3LfPlSxssLx3btVtguCgJUzZmBLOR3pV1qGcZcnsdHR8M8x11XOx4IgYO3atQgLC4OpqSnc3d0RGBio1sE4Li4O8+bNQ3R0NKytrTFgwAD4+vqqtIw6Oztj7dq1+PXXX/H27VvY2NigdevW2LRpE2rWrIlnnOuKqFwpVNAzZcoUREREoGHDhqhevbra6C2ZTIYbOfqVkPYYm5rCY8AA7Fq5Eo4NGoiru5/KGvGikzXiRRAELP78czx7+FBldfdfAgJQwcJCXN09uwS5HJeCg+HWpw+O79ql9nxMZGTxFo7KFWt7e5zKZa6rl8+eYYiLCwBgwIABGFCAua78/f3h7+8PIHNen9xWeO7evTvn3iGSkEIFPc2bNxeHhlLZoMnV3bM7e+gQ0lJS0HP4cJzetw9p2b5YZDo6sGdLIBERlRKFCnqCgoI0nA0qbsrV3aMjIxH/5g1snZzw19GjBV7d/bes1d0rWlqqPHdy3z7Y16mDxq1aYdqqVQiYOBG6+vpIS0mBsakpvv7ppxIrIxER0bsUakZmKrssra3h2KABdPX03rm6e3Y5V3dXiomMxI0//0SXfv0AAJ379MG6kycxePJkAMAPv/0Gx/r1S6BURERE+StUSw+VLgUZ8aLJ1d2VTh04gIyMDHTu21fcVqtRI5hVrIitP/6ISlWrFnPJiYiICo5BTzlQ0BEvmlzdHQBO/fYb6jdrBltHx+IpGBERkQYx6CkHco54Uc4UrJw/JjoiAs5t2qBNt275zhRc0NXdAWDtiRNFyjcREVFJYp+ecmhvYCCGuLjgM09PAMBnnp4Y4uKCvYGBWs4ZERGR9rClpxxSzhScE2cKJiIiKWPQUw7lnCmYiIiIeHuLiHLx/NEjLJs6Fd6urnC3tsbYDh3U0giCgF0rV2Jo8+boZmuLsR064HTWLN/ZJcjl+NHXF1716qGHvT2+6t1bbVqE+9euYckXX2DEhx+ih709hrdsibVz5uBtfHyxlZGIpIctPUSk5kloKC6eOIEGzZpBEARxEdnsdq9ahY0BARju64sPWrbEn0eOYIGPDwxNTNC2Wzcx3Xc+Pgi9fh0+c+agUtWq2LtuHab07Yv1ISGoZmsLADh94ADC//sPAz/9FDVq18azhw+x6fvvcfvvv7Hy8OESKzcRlW8MeohITZtu3dCuRw8AwKLJkxGaYy299LQ0bFu6FH29vTHKzw8A0MLVFS/Dw/HLwoVi0HP38mVcCg7Ggm3bxG1N27fHsBYtsGf1akwOCAAADPnsM5XZvpu2a4dK1arhmxEjcOvCBVg7OBR7maVMLpdDLpcjMmutvIiICAiCAHNzc5jzVjmVI7y9RURqcq6xltPL8HAkJiSghZubyvaWnTvj0Z07eJm1Ovm/t25BJpOhhaurmMbIxASNW7XChePHxW05lzcBgDqNGwMAol++RGx0NAAgLS2tUOWhdwsMDISLiwt6ZAW6PXr0gIuLCwI54pPKGbb0ENF7S0tJAQDoGxiobNc3NAQAhD14ACs7O6SmpECmowNdXV21dJFhYUhJSoKhsXGux7h18SIA4Mj27bh69iwAoGfPnti9ezcaZwVEpBkTJ07E0KFDIZPJYGRkhOTkZLGlh6g8YUsPEb03K3t76Ojo4P7Vqyrb712+DACQx8YCAGo4OSFDocCDmzfFNBkZGQi9dg2CICAhLi7X/SfExWFjQACqVq+OmxcuiNujo6MxcOBAvH37VsMlkjZzc3PY2dnBzs4ONWrUEP/PoIfKGwY9RPTejE1N4TFgAHatXIlLwcGIj43F8d27cSpr9JZO1jptLdzcYOfkhGVTp+LR3bt4ExWFwLlz8eLpUwCALJfbaAqFAt9NnIjkxEQIgoD0bLe0BEFAdHQ07t+/XwKlJKLyhkEPERXKpPnzUa9pU8wYMgS969TB2rlzMWb6dABA5WrVAAB6+vqYs2EDkhMTMb5TJ/Rr2BBXz55FPx8f6Onrw7xSJbX9LvX1xY2//kLAjh0wyuPWl0GO22pERAXBPj1EVCjmlSph0Z49iI6MRPybN7B1csJfR49C38AAtZ2dxXS1GzfG5gsX8PzRIwgA7Jyc8NO0aajbpAn09PVV9rnO3x/H9+zB/C1bUN/FBV7jx2PNN98gQ6EAAOjp6aF+/fpo0KBBSRaViMoJBj1EVCSW1tawtLaGQqHAoaAguHp5wbRCBZU0MpkMdrVqAQBio6Nx+uBB+MyZo5Jm18qV2LN6Nb7+6Se09vAAAPQZPx7paWnYsWwZ5LGxaN26NdavXw89PV66iOj98cpBRGqSExNxKTgYAPDy2TMkxsfjzKFDiIuJEdME792LlKQk2Do5ISYyEv/bvBmRYWGYtXatyr62LV0KW0dHVKpWDeH//Yfty5ejrrMzug0ZIqY5+dtv+Pnbb+HWpw9q1K6Nu1kdogHAzcsLnXr1whAXF6xcuRKWuQxvJyIqCAY9RKQmNjoa/uPGqWzL+VgQBPy6di0iwsJgbGqKVu7umBUYiCrW1irp4uPiEDhvHmKjo1HZygoeAwZghK+vylxAl0+fBgCc3r9fbSmLkX5+uS6gS0T0vspM0LNt2zYcOXIEGRkZ6NChAyZMmJBnE/fMmTMRGhqqMjfInj17SiqrRGWetb09TkVFqW1/+ewZhri4AAA8BgyAx4AB+e5rkr8/Jvn7vzPNtFWrMG3VqjyfV052SERUFGUi6Dl+/DjOnj2LpUuXwsjICPPnz8eePXsw9B2//saPHy/OLkpERERUJoasBwcHw8vLC1ZWVrCwsMDAgQMRnNXfgIiIiKggykRLT1hYGGrWrCk+dnR0RHR0NN6+fQtTU9NcX7N9+3Zs27YN1tbWGDRoED788MNc00VERCAiIkJ8bGhoCOscfRKKSpY1UZvy3/Iue3kLU+aivK4ory8sbRxXSmXVxHFZp97/uCWppMubfXHVyMhIyGSy915clXWqdB83L1oPehRZ82/kRVdXF8nJySrBjfL/SUlJuQY9o0aNQo0aNaCvr49//vkHS5YswXfffYe6deuqpV23bh38s/U3mDZtGubkGEqrCfo55iMpz4yMjMR/jfOYXO5ddHV11dZqKonjFpY2jmuYtcaVoaFhuS+rJo7LOlVw2rhWlXR5ly5dioCAAAAQu0HMnDkTs2bNKvA+WKcKRlvXqrxoPej55ptvcPv27Vyfq1ixIrZs2QIjIyMkJiaK25X/z+sE1qtXT/x/27ZtcenSJVy4cCHXoMfHxwe9evUSHxsaGiIpKalQZcmLTCZDGiAu4lfeJScni/8W5lwaGBggNTW1xI9bWNo4bkrWgp8pKSnlvqyaOC7rVMEof4mX9LWqpMs7fvx4DBw4UGWbubn5ex2bdapgSvJaZWJikm8arQc9ymj7Xezt7fH48WNxFtbHjx/D0tIyz1tbOeno6OT5AbaxsYGNjY34ODo6utg+7IIgSCLoUZaxsOUtyuuK8vrCKs7jPn/0CHvWrMHdy5fx+P592Nepg1/OnVM5ZkZGBnavWoVDQUGIiYyErZMTRvj6wq1PH5V9JcjlWDdvHs798QdSkpJQ38UFny5YgNrZVixPS03FLwsX4u7ly/j35k0kJyZi//37sKhSpdjL+i6sUyWrvJe3QoUKqJBjAs3s+SgI1qnSfdy8lImOzF26dMGhQ4fw6tUryOVy7N69G+7u7rmmTUhIwJUrV5CSkgKFQoFLly7h/PnzefbpISrNnoSG4uKJE7B1dETNbC2Y2e1etQobAwLQbfBgfLd9O5q0bYsFPj7469gxlXTf+fjgzyNH4DNnDuZu2ABdPT1M6dsXr54/F9OkJCXhj61bYWBkhMatWhVr2UidXC7Hs2fPxH6GERERePbsGeRyuZZzRlQ+aL2lpyC6du2KqKgofPXVV1AoFOjYsaNK0+S8efPQsGFDDBw4EAqFAjt27MCzZ88gk8lgY2ODr776Cg0bNtRiCYgKp023bmiX1edg0eTJCL1xQ+X59NRUbFu6FH29vTHKzw8A0MLVFS/Dw/HLwoVo260bAODu5cu4FByMBdu2iduatm+PYS1aYM/q1Zic1eJqZmGBg//+C5lMhqM7d+KfrEkDqWQEBgZiyZIl4mNPT08AgJ+fH77++mttZavYyOVyyOVylSAPwHt3KiYqqDIR9MhkMgwfPhzDhw/P9fl58+aJ/7ewsMCPP/5YQjkjKl7ZZy3Ozctnz5CYkIAWbm4q21t27oyVM2bg5bNnsLKzw7+3bkEmk6GFq6uYxsjEBI1btcKF48fFoAd498id1Kx+ARn5DECgwpk4cWKu84+V1wBAakEeaV+ZCHqIKHepWZ0E9Q0MVLbrZ42YCHvwAFZ2dkhNSYFMR0dttIm+oSEiw8KQkpQEw3eMrBAEAUGLF2Nr1g+KiV26YP7WrXBu00aTxZE8qbVwSC3II+0rE316iCh31jVqQEdHB/evXlXZfi9rwU55bCwAoIaTEzIUCjy4eVNMk5GRgdBr1yAIAhLi4t55nMPbtmH7smVAVkfE+Lg4TB88GFHZ5rgiel/m5uaws7NT+2PQQ8WFQQ9RGWZsZgaPAQOwa+VKXAoORnxsLI7v3o1TWYt26mTdqmrh5gY7JycsmzoVj+7exZuoKATOnYsXT58CAGT53EY787//qd3SUqSn4+ZffxVDqYiIigeDHqIybtL8+ajXtClmDBmC3nXqYO3cuRgzfToAoHK1agAAPX19zNmwAcmJiRjfqRP6NWyIq2fPop+PD/T09WFeqdI7j5Hb4r6CIBRqcjYiIm1hnx6iMs68UiUs2rMH0ZGRiH/zBrZOTvjr6FHoGxigtrOzmK5248bYfOECnj96BAGAnZMTfpo2DXWbNIFePrPwdh86FJdOnoSQkQEgs2XIpEIFuHTsWJxFIyLSKLb0EJUTltbWcGzQALp6ejgUFARXLy+Y5piATSaTwa5WLdSoVQtxMTE4ffAgeuYxKjK7jh99BN8ffoBpVl8LO0dHLDt4EBaVKxdLWYiIigNbeqhQ8popODtBELBr5cp8ZwqOi4nBLwsX4tLJk4iLiYGNvT16jhiBvhMmqAyffnjnDjZ+9x3uX7uGtJQU1KxfHyOmTMGHXbqUSJm1ITkxEZeCgwFkDU+Pj8eZQ4cQ9/q1mCZ4716kJCXB1skJMZGR+N/mzYgMC8OstWtV9rVt6VLYOjqiUrVqCP/vP2xfvhx1nZ3RbcgQlXSXgoORnJiI0OvXAQB/HTsGEzMzfPDhh9hw5gyGuLhg8d69sLKzK97CExFpGIMeKhTlTMENmjUTl0LI6ffNm7F79WoM9/XFBy1b4s8jR7DAxweGJibiBHkAMGf0aDx/9AjjZs2CVY0auH7uHFbPno2MjAwMmDQJAPD61StM7dcPNg4OmLJ0KQwMDXFw0ybMGj4cP/3xBxo0a1ZiZS9JsdHR8B83TmVbzseCIODXtWsRERYGY1NTtHJ3x6zAQFSxtlZJFx8Xh8B58xAbHY3KVlbwGDAAI3x91eYCWv7113gZHi4+XvLFFwCAkX5+6JHL8GIiorKCQQ8VSn4zBQPAvvXr850pOOrFC9y6eBF+K1aIX6itu3TBg5s3cfrAATHouXLmDOJiYrDm2DHYODgAyJxRuF/Dhjj3++/lNuixtrfHqagote2vnj/H4KZNAQAeAwbAY8CAfPc1yd8fk/z98023M8fw9+xePnuW7+uJiEor9umhQslvpmAg89ZMbjMFP7pzR/zyTE9LAwCxr4iSmYWFyuJ0uaXTNzCAoZERMjIycHjbNvz87bcAoLKWlKY8f/QIy6ZOhberK9ytrTG2Qwe1NMrbeUObN0c3W1uM7dABp7OGjmeXIJfjR19feNWrhx729viqd2/8d+uWWrrEhAQsnTIFXnXrwtPBAbOGD0dkthYYIiJ6Pwx6qFi9a6ZgALBxcEBzV1dsX7YMj+/dQ2JCAs787384f/gw+mS7jdOuRw9UqloVa2bPRnRkJOJev0bQ4sVITEjAi0ePsMzPTwwwvh4wAGH//qvRchRk4c/fN2/W2MKfALBgwgRcOHYMn33/PeasX4/oiAj49euHlKQkjZaNiEgqeHuLio1MJsP9q1fRtF07cVvOmYIB4NugIMz39sa4rOHPMpkME+bORddBg8Q05pUqYcX//oeZw4ZhYOPGADJbfSb6+2PZ1Kkqx01NScE6f398t22bxsqiqdt5BV34896VK7h44gQCduxAaw8PAIBjw4YY3rIlju3eLW4jIqKCY0sPFZsOPXvmO1OwIAhY/PnnePbwIWYFBmLZgQMYM306fgkIwOFsQcubqCh8M2oUrGvUQMCOHVj8669o7eGBtXPmQCfHBHlCRgYiw8I0WhZN3c7Lb+FPpUvBwTCzsEArd3dxm5WdHRq3aoWLJ07gbXw8ACApIaEoxSIikhS29FCxGTF1KpLevsWMrCHRFlWqYMz06QicO1ecKfjiiRM4c+gQ1oeEoNYHHwAAWrq5Ie71awTOm4fuQ4dCR0cHu1atQkJsLNYFB8PAyAgA0LxTJ4xp3168Vaakq6eHmvXrl2BJ/5+mFv58+uABatSurbbiuUPdujh98CDGZ7WKjevUCVOWLuWoKiKiAmBLDxUbMwsLLNqzB3tu3cLGs2ex+8YNVLO1VZkp+GloKHR0deHUsKHKa2s3aoSEuDjIs+ajeRoaihp16ogBD5B5G6xBs2awqFIFMh0dMeCoUKkSJmV1ai5Jytt52RV24c+EuDiY5bLoovzNG8hfvxY7eWcoFPjhyy9x559/iqNIRETlCoMeKnbvminYys4OGQqF2uil0OvXYWxqCosqVTLT1aiBsAcPVDrxZmRk4N9bt1DfxQVLfv0Vfb29AQA/7N0Lyxxz1JSEgtzOe5+FP3O28gDAiydP1LbpGRjgypkzxVAiIqLyhbe3qFDymikYAKrXrAkAOH/4MIyMjd85U3BrDw9Y29vDf9w4jPTzg6W1NW789RcOBQVh4CefiF/8H40YgcPbtmHW8OHo4+0NPT09HNmxA4/u3MGEb75Bs44dYevkhN2rV6sNfy8pBbmdp1z481tvb4zv1AkA4NSwIfr5+GD/+vXiwp9mFha5Dr1X5FjpHAAgCGq31YiISB2DHiqUd80U/M369QAKNlOwsZkZfvjtN2wMCMCGBQsQHxsLGwcHjJ89G/19fMR0dZydsWTvXmxZsgRLvvgCivR0ONSti/lbtpSaZSiUt/M0sfCnQ926uHLmDARBUGnx0dXVBWQyQDmHkUwGHV1ddOrVqySLSkRUJjHooULJa6Zg4P9n7e3Qs6dK4JKX6jVr4puffxYfGxgYIDU1VS1d03btVIa/l1aW1tawtLaGQqHId+FPIDOAPH3wIHzmzBGfb+Xuji0//IB/Tp0Sg7pXz5/jv9u30T9rnp+Ip09Ro1YtzFizRmxdIyKivDHoISoATd3OAwq28GeD5s3R2sMDS778EpO+/RYmZmYIWrQI1jVqYNzMmeg/cSIGN22KJXv3opqtbcmcBCKiMo5BD1EBaOp2HlDwhT9nrVuHwLlzseLrr5GWmgqXDh0wb9MmGBobF2NJCQDkcjnkcjkiIiIAQPzX3Nwc5lrqM0ZERceghzQmQS7HW7kc0VlfEMp/Tc3Ncx1+XZZo8nZeQRf+NK1QAVOWLsWUpUvfL7NUZIGBgViyZIn42NPTEwDg5+eHr7/+WlvZIqIiYtBDGrM3MBBbsn1RfJb1RTHSzw+j+UVBZcjEiRMxNJcJH9nKQ1S2Meghjek/cWKuMwNrawg5UWHxNhZR+cSghzTGrBzcxiIiKk/YP00VZ2QmIiIqpwIDA+Hi4iL2S/P09ISLiwsCAwO1nDPtYEsPERFROcX+aaoY9BAREZVTUr2NlRfe3iIiIiJJYEsPERFRMZNah2JleSMjIwFkllcQBK2Xly09RERExUxqHYqV5e3RowcAoEePHqWivGzpISIiKmZS61CsLK9MJoORkRGSk5PFlh5tYtBDRERUzLR9W6ekKcsrk8lgbGyMpKQkCIKg7Wzx9hYRERFJA4MeKlOeP3qEZVOnwtvVFe7W1hjboYNaGkEQsGvlSgxt3hzdbG0xtkMHnN6/Xy1dXEwMlk2disEuLuhhb48x7drht3Xr1H6NPH3wANMHD4angwN61a6NgE8+QVxMTLGVkYiIigdvb+Wgq6sLXV1dje5TJpNBV1cXBgYGpaJ5r7TT08u7Wj57+BCXgoPRsHlzAEBGRgYMDAwAAPr6+gCAw1u3YufKlRjl54fGH36Ic4cPY4GPD0zNzdE+q1MdAMwbOxbPHj7EhDlzYF2jBv45fRqrZ8+GTCbD4MmTAQBv5XJM7dsXljY2mLthA5KTkrDO3x8zhw3DuuBg6OjoiMfV19cX81LctHFMbRw3IS4OCXI5YqOjAQDR0dHQ19eHubk5LCwsCryfd9Up+n+8VhUc61TBlLY6xXctB4VCAYVCodF9ymQyKHR1kZqaWire9LIgNTU11+0tu3TBruvXAQCLJk9G6I0bYtq0tDQAwN5169DX2xvDfX0BAE3at8eLJ0/w8/z5+LBLFwBA1IsXuPHXX/BbsQIeAwcCABq3aYMHN28i+Lff0HfCBADAr+vWIUEux8+nT6NS1aoAAGt7e0zy8MDpgwfRoWdP8bhpaWl55lvTtHFMbRx3x8qV2LJkifjYw8MDAODn54evv/76vfZVkueprFJ+QfFaVTCsU/krbXWKQQ+VKTo6+d+RTU5MRAs3N5VtLTt3xsoZM/Dy2TNY2dkhPevLO+cK8GYWFoh7/Vp8/N+tW6jVqJEY8ABAvaZNYV65Mi4cO4YOPXsiJSmpKEWid+g/cSJ6ZI14qZ9tu5Q6hBKR5jDooXJJP8etF31DQwBA2IMHsLKzg42DA5q7umL7smWoUasWrGrUwJUzZ3D+8GH4/vCD+LrUlBTxlk7O/d+7ehUDGjVCzMuXAICQAwcwKOu2GGmGmbk5zLICHDst54WIyj52ZKZyRyaT4f7Vqyrb7l2+DACQx8aK274NCkIVKyuM69gRHzk6Yt6YMRgzfTq6DhokprGrVQuP791Tac15+ewZXr98ibAHD/D61Stx+zp/f/x98mQxlYqIiIqKQQ+VOx169sSulStxKTgY8bGxOL57N05ljd7SkckAZI7wWvz553j28CFmBQZi2YEDGOnnh18CAnB42zZxXx+NGIHEhAQsnToVUREReP74MRZ99hkgk0HI2k92pw8cKKliEhHRe+LtLSp3RkydiqS3bzFjyBAAgEWVKhgzfToC585F5WrVAAAXT5zAmUOHsD4kBLU++AAA0KRdOyTExSFw3jx0HzoUOjo6sKtVC1//9BNWzpiBE3v2AMgMqmrWq4cn9+8jZ7e80tBRj4iIcsegh8odMwsLLNqzB9GRkYh/8wa2Tk746+hR6BsYoLazMwDgaWgodHR14dSwocprazdqhN/i4iB//RoVLS0BAO79+6NTr1549vAhzCpWRFUbG4z48EPIdHQgZGQA2QId1169Sq6gRET0Xnh7i8otS2trODZoAF09PRwKCoKrlxdMK1QAAFjZ2SFDocB/t26pvCb0+nUYm5rCokoVle36BgZwbNAAVW1scPXcOUSGh8NvxQpUzJZu/OzZaN21a/EXjIiICoUtPVSmJCcm4lJwMIDMDsWJ8fE4c+gQAKB6zZoAgPOHD8PI2Bi2Tk6IiYzE/zZvRmRYGGatXSvup7WHB6zt7eE/bhxG+vnB0toa186dw6GgIAz85BPIsvr+JL19i81LlsC5TRsYGBri3pUr2LFiBUb5+aHboEHoOnAgnoSGYlyHDujSr1/JngwiInovDHqoTImNjob/uHEq25SPv1m/HkBmv5pf165FRFgYjE1N0crdHbMCA1HF2lp8jbGZGX747TdsDAjAhgULEB8bC2t7e4yfPRv9fXzEdDo6Onh89y6O7tyJpLdvYV+7Nr5YtAjds/oLyWQymJiZFXexiYhIAxj0UJlibW+PU1FRuT738tkzAJkdjbMHLnmpXrMmvvn553emMTQ2xqKsDsxERFS2MeghonwlyOV4K5cjOiICAMR/TbNNHkhEVNqxIzMR5WtvYCCGuLjgM09PAMBnnp4Y4uKCvYGBWs4ZEVHBsaWHiPKVfQ2s7HKuXUZEVJox6CGifJnxNhYRlQO8vUVERESSwKCHiIiIJIFBDxEREUkCgx4iIiKSBAY9REREJAkMeoiIiEgSGPQQERGRJDDoISIiIklg0ENERESSwKCHiIiIJIFBDxEREUkCgx4iIiKSBAY9REREJAkMeoiIiEgSGPQQERGRJDDoISIiIklg0ENERESSwKCHiIiIJIFBDxEREUkCgx4iIiKSBD1tZ4CICi5BLsdbuRwxkZEAgOiICAiCAFNzc5iZm2s5d0REpVu5DHpu3ryJ3bt34+HDhzAwMMCWLVu0nSUijdgbGIgtS5aIjyf36AEAGOnnh9Fff62tbBERlQnlMugxMjKCu7s7OnXqhG3btmk7O0Qa03/iRPQYOhQymQxGRkZITk4WW3qIiOjdymXQU7duXdStWxe3bt3SdlaINMos6zaWTCaDsbExkpKSIAiCtrNFRFQmlMug531EREQgIiJCfGxoaAhra2uNHkMmk6n8S+8mk8kKda6yn+eSPNfaOK4U61RRSlrSdaKskmK9KizWqYIpbXWqzAU9CoXinc/r6uq+1/7WrVsHf39/8fG0adMwZ86cQuXtXfT19TW+z/JKV1f3vd9HIPO2pvJfY2NjTWer1B1XanWqKGe2sHVKiqRWrwqLdargSlOdKnNBzzfffIPbt2/n+lzFihXfu9Oyj48PevXqJT42NDREUlJSkfKYk0wmQxog9r+gdzMwMEBqamqB0ytHNEVntdg9e/wYycnJJTaiKTk5WfxX03UnL8pfTVKqU0U5s+9bp6RKivWqsFinCqYk65SJiUm+acpc0BMQEKDR/dnY2MDGxkZ8HB0dXWxvjCAIvJAUwPuep1/XrtXqiCZlXrXx/kqpThWllFI6T5rA85U/nqP3U1rOV5kLegoiIyMD6enpSE9PBwCkpqZCJpOVqiY20hzliKacOKKJiIiyK5dBz507dzBr1izxcf/+/VGtWjVs2LBBi7mi4mLGifmIiKgAymXQ07hxYxw6dEjb2SAiIqJShGtvERERkSQw6CEiIiJJYNBDREREksCgh4iIiCSBQQ8RERFJAoMeIiIikgQGPURERCQJDHqIiIhIEhj0EBERkSQw6CEiIiJJYNBDREREksCgh4iIiCSBQQ8RERFJAoMeIiIikgQGPURERCQJDHqIiIhIEhj0EBERkSQw6CEiIiJJYNBDREREksCgh4iIiCSBQQ8RERFJAoMeIiIikgQGPURERCQJDHqIiIhIEhj0EBERkSQw6CEiIiJJYNBDREREksCgh4iIiCSBQQ8RERFJAoMeIiIikgQGPURERCQJDHqIiIhIEhj0EBERkSQw6CEiIiJJ0NN2BojKqgS5HG/lckRHRACA+K+puTnMzM21mTUiIsoFg54cdHV1oaurq9F9ymQy6OrqwsDAAIIgaHTf5ZGeXtmolvvXr8em778XH3/m6QkAGDN9OsbNmFGsx5ZinTIowmvLSp3SNinWq8JinSqY0lanZEJpyEUpEh0drfF9ymQy/GdsjKSkpFLxppd2BgYGSE1N1XY28qVs6cmpJFp6ZDIZjCVWpxoX4bVlpU5pmxTrVWGxThVMSdYpS0vLfNMwVCUqJDPexiIiKlPYkZmIiIgkgUEPERERSQKDHiIiIpIEBj1EREQkCQx6iIiISBIY9BAREZEkMOghIiIiSWDQQ0RERJLAoIeIiIgkgUEPERERSQKDHiIiIpIEBj1EREQkCQx6iIiISBIY9BAREZEkyARBELSdifIuIiIC69atg4+PD2xsbLSdHSoHWKeoOLBekaaVtjrFlp4SEBERAX9/f0RERGg7K1ROsE5RcWC9Ik0rbXWKQQ8RERFJAoMeIiIikgQGPSXAxsYGc+fOLRX3M6l8YJ2i4sB6RZpW2uoUOzITERGRJLClh4iIiCSBQQ8RERFJAoMeIiIikgQ9bWegrPn9999x6tQpPHnyBG3atIGfn5/43JUrV7B582ZERkbCxsYGkyZNQv369VVee/DgQcjlclhZWWHUqFFo3rw5AODkyZNYuXIlDAwMxPSffPIJXF1dS6xspB1FqVP79+/HH3/8gfj4eDRq1AiTJ09GpUqVxOe3bduGI0eOICMjAx06dMCECROgp8ePvRQUV73itUqa0tLSEBgYiBs3biA+Ph6WlpYYMGCA+L4/ffoUK1euxJMnT2BtbY1Jkybhgw8+EF//559/IigoCLGxsWjQoAG++OILVKlSRXy+xK5VAr2XP//8U7hw4YKwdu1aYfHixeL2Fy9eCAMHDhSuX78upKenC8eOHROGDh0qxMfHC4IgCPfu3RP69+8vhIaGChkZGcK5c+eE/v37C3K5XBAEQQgODhamTJmilTKRdhW2Tp0+fVrw9vYWIiIihJSUFGHFihXCzJkzxdcfO3ZM8Pb2FiIjI4XY2FhhypQpwvbt20u8fKQdxVWveK2SpqSkJGHbtm1CRESEoFAohDt37giDBg0S7t27J6SlpQnjx48X9u7dK6SmpgqnTp0ShgwZItap8PBwYeDAgcK1a9eE5ORkYe3atcL06dPFfZfktYq3t95T27Zt0bp1a5ibm6tsv3r1KurVq4cmTZpAV1cXXbt2hbGxMS5evAgAiIyMhL29PerWrQuZTIb27dtDR0cHkZGR2igGlSKFrVMXL16Eu7s7rK2tYWBggMGDB+PWrVtinQoODoaXlxesrKxgYWGBgQMHIjg4uMTLR9pRXPWKpMnIyAjDhg2DtbU1dHR00LBhQzRo0AD37t3DrVu3kJKSgj59+kBfXx9ubm6wsrLCX3/9BQA4ffo0mjVrhqZNm8LQ0BDDhg3D/fv3xVmaS/JaxaBHQwRBgJDL6P8nT54AAFq0aIH09HTcu3cPCoUCZ86cgZmZGezt7VXSDh8+HBMmTEBQUBBSUlJKKvtUCuVXp3I+p3ysfD4sLAw1a9YUn3d0dER0dDTevn1bLPmlsqGo9Ur5f16rpC05ORn//fcfHBwcEBYWBgcHB+jo/H9I4ejoiLCwMACZt74cHR3F5ypUqICqVavi6dOnAEr2WsWgR0OaNm2Ke/fu4erVq0hPT8fRo0cRFRUlXgxMTEzQrl07zJo1C/369cOaNWswefJkGBoaAgA++OADrFq1Clu2bMGcOXNw+/ZtBAUFabFEpG351anmzZvjxIkTePHiBVJSUrBz507IZDLx+eTkZJiamor7U/4/KSmp5AtDpUZR6xWvVSQIAlasWIE6derAxcUFSUlJKtcaIPN6o7zWJCcnw8TE5J3Pl9S1ij0aNcTOzg6+vr7YtGkTYmJi0LJlSzRp0kTsqHX8+HEcP34cy5cvh52dHUJDQxEQEAB/f384OTnB2tpaZV8jR47EDz/8AB8fH20VibQsvzrl7u6OmJgYzJkzB6mpqfDy8oKxsTEsLS0BZDZHJyYmivtT/t/Y2LjkC0OlRlHrFa9V0iYIAtasWYOYmBh8++23kMlkMDY2VrnWAJnXG+W1Jue1CADevn2b5/PFea1i0KNBbdu2Rdu2bQEACoUC3t7e6NOnD4DM5uAWLVqIt7MaNGiAOnXq4MaNG3ByclLbl46OTq5N0CQt76pTOjo6GDJkCIYMGQIAePbsGXbt2gUHBwcAgL29PR4/fowGDRoAAB4/fgxLS0u1X2QkPUWpVznxWiUdgiAgMDAQjx49wvz582FkZAQg81qzb98+ZGRkiLe4Hj9+jO7duwMAHBwcVG6PJiQkIDo6WivXKt7eek8KhQKpqanIyMhARkYGUlNTkZ6eDgD4999/oVAokJCQgPXr16NatWpwcXEBANSrVw9XrlzB8+fPAQAPHjzAvXv3xPuYV65cwevXrwFkdnrevHkz2rRpU/IFpBJX2DqVkJCAFy9eQBAEREZGYtWqVejduzfMzMwAAF26dMGhQ4fw6tUryOVy7N69G+7u7lorJ5Ws4qpXvFZJ17p16xAaGgp/f3+V21WNGzeGvr4+Dhw4gLS0NJw5cwaRkZFivXB1dcWVK1dw48YNpKSkYPv27ahXr564HldJXqu49tZ72rFjB3bt2qWyrXPnzvjyyy8xc+ZMPHz4EDo6OmjZsiXGjx8vjpwQBAG7du3CyZMnER8fDwsLC/Ts2RO9e/cGAGzatAmnT59GUlISzM3N0bZtWwwbNkyMpKn8KmydioiIwLfffouoqCiYmZmha9euGDx4sPhLSxAEbN++HUeOHIFCoUDHjh05T4+EFFe94rVKml69eoXx48dDX18furq64vb+/ftj4MCBePLkCVatWoUnT57AysoKkyZNQqNGjcR058+fx+bNm/HmzRs0bNhQZZ6ekrxWMeghIiIiSeDtLSIiIpIEBj1EREQkCQx6iIiISBIY9BAREZEkMOghIiIiSWDQQ0RERJLAoIeIiIgkgUEPERERSQKDHiIiIpIEBj1EREQkCQx6iIiQuUBnWlqatrNBRMWIQQ8RlQqHDh2CTCbDv//+q7I9Li4OJiYm+OmnnwAAFy5cQOfOnWFqagoLCwsMHToUr169UnnN9OnT0bhxY5iZmcHW1hZDhgxBRESEShpXV1d89NFH2Lx5M+rVqwdDQ0Ncv369WMtIRNrFoIeISoWePXvC1tYWv/zyi8r2nTt3IiMjA8OHD8eFCxfg6uoKCwsL7N69Gz///DP++ecf9OrVS+U1r169wsyZM/HHH39gxYoVePLkCTp16oT09HSVdJcvX8aPP/6I+fPn4/Dhw6hRo0axl5OItIerrBNRqfHNN9/gl19+QVhYGHR1dQEAH374IZycnLBr1y4xcDl//jxkMhkA4M6dO2jcuDF+//13eHp6qu1ToVAgMjISdnZ2OHbsGLp27Qogs6XnwoULePjwIezs7EqukESkNWzpIaJSY9y4cYiIiMDRo0cBALdv38Y///yDcePGITExEX/++ScGDBgAhUKB9PR0pKeno169erCxscE///wj7ufIkSNo27YtLCwsoKenJwY1Dx48UDmes7MzAx4iCWHQQ0SlRs2aNeHh4YGNGzcCADZu3AgHBwd06dIFb968gUKhwFdffQV9fX2VvxcvXiA8PBwAxNtd1atXx9atW3HhwgVcvHgRAJCcnKxyvGrVqpVsAYlIq/S0nQEiouy8vb0xdOhQPH/+HNu3b8enn34KHR0dVKxYETKZDDNnzoSXl5fa6ywtLQEA+/fvh4WFBfbs2QMdnczfdU+fPs31WMpbZEQkDQx6iKhU6d27NypVqoShQ4ciJiYGY8aMAQCYmpqiTZs2uHfvHhYsWJDn65OSkqCvr68S0Gzfvr3Y801EpR9vbxFRqaKvr49Ro0bh7NmzcHd3h729vfjckiVL8Mcff2DQoEHYv38/QkJCsG3bNowaNQohISEAAA8PD0RGRuKzzz7DyZMnsWDBAmzevFlLpSGi0oRBDxGVOn369AGQ2bE5u7Zt2+L8+fNISEjAmDFj4OnpiW+//RYmJiaoXbs2AMDT0xOLFi3CwYMH0atXL5w9exa///57iZeBiEofDlknolJnzpw5WLNmDZ4/fw5DQ0NtZ4eIygn26SGiUiM0NBShoaFYuXIlPv30UwY8RKRRbOkholLD1dUVFy9eRPfu3bF9+3aYmppqO0tEVI4w6CEiIiJJYEdmIiIikgQGPURERCQJDHqIiIhIEhj0EBERkSQw6CEiIiJJYNBDREREksCgh4iIiCSBQQ8RERFJwv8B0NB3JiYFfioAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 640x480 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"<ggplot: (8778923054321)>"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"abortion_plot = pd.DataFrame(\n",
" {\n",
" 'sd': reg.bse['C(repeal)[T.1.0]:C(year)[T.1986.0]':'C(repeal)[T.1.0]:C(year)[T.2000.0]'],\n",
" 'mean': reg.params['C(repeal)[T.1.0]:C(year)[T.1986.0]':'C(repeal)[T.1.0]:C(year)[T.2000.0]'],\n",
" 'year': np.arange(1986, 2001)\n",
" })\n",
"abortion_plot['lb'] = abortion_plot['mean'] - abortion_plot['sd']*1.96\n",
"abortion_plot['ub'] = abortion_plot['mean'] + abortion_plot['sd']*1.96\n",
"\n",
"(\n",
" p.ggplot(abortion_plot, p.aes(x = 'year', y = 'mean')) + \n",
" p.geom_rect(p.aes(xmin=1985, xmax=1992, ymin=-np.inf, ymax=np.inf), fill=\"cyan\", alpha = 0.01) +\n",
" p.geom_point() +\n",
" p.geom_text(p.aes(label = 'year'), ha='right') +\n",
" p.geom_hline(yintercept = 0) +\n",
" p.geom_errorbar(p.aes(ymin = 'lb', ymax = 'ub'), width = 0.2,\n",
" position = p.position_dodge(0.05)) +\n",
" p.labs(title= \"Estimated effect of abortion legalization on gonorrhea\")\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Remaining Differences\n",
"\n",
"The variance-covariance matrix isn't the same, because we're not using the PQR decomposition. I'll work that out but haven't had a chance yet. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Rank-Revealing Column-Pivoted QR (PQR) Decomposition\n",
"\n",
"Statsmodels doesn't yet have this option, though I'll make a PR.\n",
"\n",
"References: \n",
"[[1](https://johnwlambert.github.io/least-squares/)],\n",
"[[2](https://cran.r-project.org/web/packages/RcppEigen/vignettes/RcppEigen-Introduction.pdf)]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Unweighted Model"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This matches the no weights version of `estimatr` exactly"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Intercept 7.30380009\n",
"C(repeal)[T.1.0] -1.96618112\n",
"C(year)[T.1986.0] -0.02291998\n",
"C(year)[T.1987.0] -0.16124192\n",
"C(year)[T.1988.0] -0.06897324\n",
"C(year)[T.1989.0] -0.10139403\n",
"C(year)[T.1990.0] -0.21528376\n",
"C(year)[T.1991.0] -0.26288372\n",
"C(year)[T.1992.0] -0.46602701\n",
"C(year)[T.1993.0] -0.57354959\n",
"C(year)[T.1994.0] -0.64230542\n",
"C(year)[T.1995.0] -0.80247626\n",
"C(year)[T.1996.0] -1.11839159\n",
"C(year)[T.1997.0] -1.13284577\n",
"C(year)[T.1998.0] -1.09917115\n",
"C(year)[T.1999.0] -1.13235332\n",
"C(year)[T.2000.0] -1.33885860\n",
"C(fip)[T.2.0] NaN\n",
"C(fip)[T.4.0] -0.55579753\n",
"C(fip)[T.5.0] 0.10271236\n",
"C(fip)[T.6.0] 0.85867309\n",
"C(fip)[T.8.0] -0.68531986\n",
"C(fip)[T.9.0] -1.08677769\n",
"C(fip)[T.10.0] -0.40742094\n",
"C(fip)[T.11.0] -1.70341156\n",
"C(fip)[T.12.0] -0.80864615\n",
"C(fip)[T.13.0] -0.48397227\n",
"C(fip)[T.15.0] 0.13965735\n",
"C(fip)[T.16.0] -1.52890281\n",
"C(fip)[T.17.0] -0.51842334\n",
"C(fip)[T.18.0] 0.03096018\n",
"C(fip)[T.19.0] -0.03016982\n",
"C(fip)[T.20.0] 0.29883407\n",
"C(fip)[T.21.0] 0.07057992\n",
"C(fip)[T.22.0] -0.64551321\n",
"C(fip)[T.23.0] -2.18688140\n",
"C(fip)[T.24.0] -1.04120293\n",
"C(fip)[T.25.0] -1.31587892\n",
"C(fip)[T.26.0] -0.60239822\n",
"C(fip)[T.27.0] -0.14719801\n",
"C(fip)[T.28.0] -0.08205431\n",
"C(fip)[T.29.0] 0.06320215\n",
"C(fip)[T.30.0] -1.36163932\n",
"C(fip)[T.31.0] -0.00580614\n",
"C(fip)[T.32.0] -1.60278891\n",
"C(fip)[T.33.0] -3.51971156\n",
"C(fip)[T.34.0] -1.78144294\n",
"C(fip)[T.35.0] -0.86451223\n",
"C(fip)[T.36.0] -0.06600649\n",
"C(fip)[T.37.0] 0.04399084\n",
"C(fip)[T.38.0] -2.01218528\n",
"C(fip)[T.39.0] -0.19859718\n",
"C(fip)[T.40.0] 0.41343931\n",
"C(fip)[T.41.0] -0.67007500\n",
"C(fip)[T.42.0] -0.14780248\n",
"C(fip)[T.44.0] -0.54828921\n",
"C(fip)[T.45.0] -1.01635253\n",
"C(fip)[T.46.0] -1.36518847\n",
"C(fip)[T.47.0] 0.18565636\n",
"C(fip)[T.48.0] -0.43811643\n",
"C(fip)[T.49.0] -0.58979804\n",
"C(fip)[T.50.0] -1.94278109\n",
"C(fip)[T.51.0] -0.65685829\n",
"C(fip)[T.53.0] 1.29610341\n",
"C(fip)[T.54.0] -0.45096998\n",
"C(fip)[T.55.0] -0.16934096\n",
"C(fip)[T.56.0] -1.86812527\n",
"C(repeal)[T.1.0]:C(year)[T.1986.0] -0.00510910\n",
"C(repeal)[T.1.0]:C(year)[T.1987.0] -0.02788920\n",
"C(repeal)[T.1.0]:C(year)[T.1988.0] -0.27242669\n",
"C(repeal)[T.1.0]:C(year)[T.1989.0] -0.23836481\n",
"C(repeal)[T.1.0]:C(year)[T.1990.0] -0.15710470\n",
"C(repeal)[T.1.0]:C(year)[T.1991.0] -0.08983659\n",
"C(repeal)[T.1.0]:C(year)[T.1992.0] 0.02352184\n",
"C(repeal)[T.1.0]:C(year)[T.1993.0] 0.19688597\n",
"C(repeal)[T.1.0]:C(year)[T.1994.0] 0.27794833\n",
"C(repeal)[T.1.0]:C(year)[T.1995.0] 0.15218360\n",
"C(repeal)[T.1.0]:C(year)[T.1996.0] -0.14035806\n",
"C(repeal)[T.1.0]:C(year)[T.1997.0] 0.06981847\n",
"C(repeal)[T.1.0]:C(year)[T.1998.0] -0.32847326\n",
"C(repeal)[T.1.0]:C(year)[T.1999.0] -0.24171052\n",
"C(repeal)[T.1.0]:C(year)[T.2000.0] -0.08237028\n",
"acc 0.00163723\n",
"ir -0.00001754\n",
"pi -0.01573799\n",
"alcohol 0.45832689\n",
"crack -0.06620465\n",
"poverty -0.00485309\n",
"income 0.00005270\n",
"ur 0.00038949\n",
"Length: 90, dtype: float64"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from scipy import linalg\n",
"\n",
"X = reg.model.exog\n",
"y = reg.model.endog\n",
"\n",
"Q, R, P = linalg.qr(X, pivoting=True)\n",
"rank = np.linalg.matrix_rank(R)\n",
"\n",
"R_inv = np.linalg.solve(R[:rank, :rank], np.eye(rank))\n",
"\n",
"effects = np.dot(Q.T, y)\n",
"\n",
"beta = np.zeros(X.shape[1])\n",
"beta.fill(np.nan)\n",
"beta[:rank] = R_inv.dot(effects[:rank])\n",
"beta = beta[np.argsort(P)]\n",
"\n",
"pd.Series(\n",
" beta,\n",
" index=reg.model.exog_names\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Weighted Model"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Intercept 7.47534301\n",
"C(repeal)[T.1.0] -1.74924971\n",
"C(year)[T.1986.0] -0.02082737\n",
"C(year)[T.1987.0] -0.28646648\n",
"C(year)[T.1988.0] -0.34946485\n",
"C(year)[T.1989.0] -0.40722001\n",
"C(year)[T.1990.0] -0.48641174\n",
"C(year)[T.1991.0] -0.53530339\n",
"C(year)[T.1992.0] -0.78662355\n",
"C(year)[T.1993.0] -0.97671663\n",
"C(year)[T.1994.0] -1.06436659\n",
"C(year)[T.1995.0] -1.35624511\n",
"C(year)[T.1996.0] -1.52006222\n",
"C(year)[T.1997.0] -1.57184605\n",
"C(year)[T.1998.0] -1.54512931\n",
"C(year)[T.1999.0] -1.58741354\n",
"C(year)[T.2000.0] -1.67257156\n",
"C(fip)[T.2.0] -0.08647049\n",
"C(fip)[T.4.0] -0.95893856\n",
"C(fip)[T.5.0] 0.08603103\n",
"C(fip)[T.6.0] 0.48783939\n",
"C(fip)[T.8.0] -0.82460495\n",
"C(fip)[T.9.0] -1.68089178\n",
"C(fip)[T.10.0] -0.74575439\n",
"C(fip)[T.11.0] -2.11235472\n",
"C(fip)[T.12.0] -1.09474761\n",
"C(fip)[T.13.0] -0.72157625\n",
"C(fip)[T.15.0] NaN\n",
"C(fip)[T.16.0] -1.79371939\n",
"C(fip)[T.17.0] -0.72282047\n",
"C(fip)[T.18.0] -0.14315369\n",
"C(fip)[T.19.0] -0.14979696\n",
"C(fip)[T.20.0] 0.08004447\n",
"C(fip)[T.21.0] -0.11037085\n",
"C(fip)[T.22.0] -0.74395266\n",
"C(fip)[T.23.0] -2.02122539\n",
"C(fip)[T.24.0] -1.47456155\n",
"C(fip)[T.25.0] -1.89432424\n",
"C(fip)[T.26.0] -0.76799387\n",
"C(fip)[T.27.0] -0.32031988\n",
"C(fip)[T.28.0] -0.02758883\n",
"C(fip)[T.29.0] -0.16692162\n",
"C(fip)[T.30.0] -1.26616853\n",
"C(fip)[T.31.0] -0.11756566\n",
"C(fip)[T.32.0] -1.96480850\n",
"C(fip)[T.33.0] -3.59097336\n",
"C(fip)[T.34.0] -2.27555404\n",
"C(fip)[T.35.0] -1.18140339\n",
"C(fip)[T.36.0] -0.59812275\n",
"C(fip)[T.37.0] -0.08849572\n",
"C(fip)[T.38.0] -1.82649080\n",
"C(fip)[T.39.0] -0.48922418\n",
"C(fip)[T.40.0] 0.30486166\n",
"C(fip)[T.41.0] -1.12281228\n",
"C(fip)[T.42.0] -0.63106777\n",
"C(fip)[T.44.0] -1.20962341\n",
"C(fip)[T.45.0] -1.11089996\n",
"C(fip)[T.46.0] -2.02895172\n",
"C(fip)[T.47.0] 0.05924823\n",
"C(fip)[T.48.0] -0.60548288\n",
"C(fip)[T.49.0] -1.26672374\n",
"C(fip)[T.50.0] -1.99735806\n",
"C(fip)[T.51.0] -0.94936056\n",
"C(fip)[T.53.0] 0.78542602\n",
"C(fip)[T.54.0] -0.52137261\n",
"C(fip)[T.55.0] -0.48367762\n",
"C(fip)[T.56.0] -2.19210632\n",
"C(repeal)[T.1.0]:C(year)[T.1986.0] -0.25905382\n",
"C(repeal)[T.1.0]:C(year)[T.1987.0] -0.34173798\n",
"C(repeal)[T.1.0]:C(year)[T.1988.0] -0.61055324\n",
"C(repeal)[T.1.0]:C(year)[T.1989.0] -0.78500028\n",
"C(repeal)[T.1.0]:C(year)[T.1990.0] -0.63163148\n",
"C(repeal)[T.1.0]:C(year)[T.1991.0] -0.55260264\n",
"C(repeal)[T.1.0]:C(year)[T.1992.0] -0.44241661\n",
"C(repeal)[T.1.0]:C(year)[T.1993.0] -0.30606258\n",
"C(repeal)[T.1.0]:C(year)[T.1994.0] -0.11788972\n",
"C(repeal)[T.1.0]:C(year)[T.1995.0] 0.02102913\n",
"C(repeal)[T.1.0]:C(year)[T.1996.0] -0.12353414\n",
"C(repeal)[T.1.0]:C(year)[T.1997.0] 0.02085286\n",
"C(repeal)[T.1.0]:C(year)[T.1998.0] -0.03605240\n",
"C(repeal)[T.1.0]:C(year)[T.1999.0] 0.01470912\n",
"C(repeal)[T.1.0]:C(year)[T.2000.0] 0.04145081\n",
"acc 0.00287774\n",
"ir 0.00044389\n",
"pi -0.03902207\n",
"alcohol 0.44684213\n",
"crack 0.05282875\n",
"poverty -0.00276098\n",
"income 0.00005847\n",
"ur -0.02782371\n",
"Length: 90, dtype: float64"
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from scipy import linalg\n",
"\n",
"X = reg.model.wexog\n",
"y = reg.model.wendog\n",
"\n",
"Q, R, P = linalg.qr(X, pivoting=True)\n",
"rank = np.linalg.matrix_rank(R)\n",
"\n",
"R_inv = np.linalg.solve(R[:rank, :rank], np.eye(rank))\n",
"\n",
"effects = np.dot(Q.T, y)\n",
"\n",
"beta = np.fiull(X.shape[1], np.nan)\n",
"beta[:rank] = R_inv.dot(effects[:rank])\n",
"beta = beta[np.argsort(P)]\n",
"\n",
"pd.Series(\n",
" beta,\n",
" index=reg.model.exog_names\n",
")"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.6"
}
},
"nbformat": 4,
"nbformat_minor": 4
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment