From ca168b2c7f54f4785b6b49be531cacbd5ca53be8 Mon Sep 17 00:00:00 2001 From: Lucy D'Agostino McGowan Date: Sun, 29 Oct 2023 19:50:55 -0400 Subject: [PATCH] update chapter --- .../chapter-11/execute-results/html.json | 4 +- .../figure-html/unnamed-chunk-8-1.png | Bin 0 -> 46446 bytes chapters/chapter-11.qmd | 203 +++++++++++++++++- 3 files changed, 204 insertions(+), 3 deletions(-) create mode 100644 _freeze/chapters/chapter-11/figure-html/unnamed-chunk-8-1.png diff --git a/_freeze/chapters/chapter-11/execute-results/html.json b/_freeze/chapters/chapter-11/execute-results/html.json index f014312..a477e45 100644 --- a/_freeze/chapters/chapter-11/execute-results/html.json +++ b/_freeze/chapters/chapter-11/execute-results/html.json @@ -1,7 +1,7 @@ { - "hash": "52cdc5c043504620e3c65dd4015038d2", + "hash": "35aeb2819ee49bb26e30753bbfdeca1e", "result": { - "markdown": "# Fitting the outcome model {#sec-outcome-model}\n\n\n\n\n\n## Using matched data sets\n\n\n::: {.cell}\n\n```{.r .cell-code}\nrnorm(5)\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n[1] -0.5790 -0.5911 0.3102 0.3006 1.4746\n```\n\n\n:::\n:::\n\n\n## Using weights in outcome models\n\n## Estimating uncertainty\n", + "markdown": "# Fitting the outcome model {#sec-outcome-model}\n\n\n\n\n\n## Using matched data sets\n\nWhen fitting an outcome model on matched data sets, we can simply subset the original data to only those who were matched and then fit a model on these data as we would otherwise. For example, re-performing the matching as we did in @sec-using-ps, we can extract the matched observations in a dataset called `matched_data` as follows.\n\n\n::: {.cell}\n\n```{.r .cell-code}\nlibrary(broom)\nlibrary(touringplans)\nlibrary(MatchIt)\n\nseven_dwarfs_9 <- seven_dwarfs_train_2018 |> \n filter(wait_hour == 9) \n\nm <- matchit(\n park_extra_magic_morning ~ park_ticket_season + park_close + park_temperature_high,\n data = seven_dwarfs_9\n)\nmatched_data <- get_matches(m)\n```\n:::\n\n\nWe can then fit the outcome model on this data. For this analysis, we are interested in the impact of extra magic morning hours on the average posted wait time between 9 and 10am. The linear model below will estimate this in the matched cohort.\n\n\n::: {.cell}\n\n```{.r .cell-code}\nlm(wait_minutes_posted_avg ~ park_extra_magic_morning, data = matched_data) |>\n tidy(conf.int = TRUE)\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n# A tibble: 2 × 7\n term estimate std.error statistic p.value conf.low\n \n1 (Inte… 67.0 2.37 28.3 1.94e-54 62.3 \n2 park_… 7.87 3.35 2.35 2.04e- 2 1.24\n# ℹ 1 more variable: conf.high \n```\n\n\n:::\n:::\n\n\nRecall that by default `{MatchIt}` estimates an average treatment effect among the treated. This means among days that have extra magic hours, the expected impact of having extra magic hours on the average posted wait time between 9 and 10am is 7.9 minutes (95% CI: 1.2-14.5). \n\n## Using weights in outcome models\n\nNow let's use propensity score weights to estimate this same estimand. We will use the ATT weights so the analysis matches that for matching above.\n\n\n::: {.cell}\n\n```{.r .cell-code}\nlibrary(propensity)\n\nseven_dwarfs_9_with_ps <-\n glm(\n park_extra_magic_morning ~ park_ticket_season + park_close + park_temperature_high,\n data = seven_dwarfs_9,\n family = binomial()\n ) |>\n augment(type.predict = \"response\", data = seven_dwarfs_9)\nseven_dwarfs_9_with_wt <- seven_dwarfs_9_with_ps |>\n mutate(w_att = wt_att(.fitted, park_extra_magic_morning))\n```\n:::\n\n\nWe can fit a *weighted* outcome model by using the `weights` argument. \n\n\n::: {.cell}\n\n```{.r .cell-code}\nlm(wait_minutes_posted_avg ~ park_extra_magic_morning, \n data = seven_dwarfs_9_with_wt,\n weights = w_att) |>\n tidy()\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n# A tibble: 2 × 5\n term estimate std.error statistic p.value\n \n1 (Intercept) 68.7 1.45 47.3 1.69e-154\n2 park_extra_ma… 6.23 2.05 3.03 2.62e- 3\n```\n\n\n:::\n:::\n\n\nUsing weighting, we estimate that among days that have extra magic hours, the expected impact of having extra magic hours on the average posted wait time between 9 and 10am is 6.2 minutes. \nWhile this approach will get us the desired estimate for the point estimate, the default output using the `lm` function for the uncertainty (the standard errors and confidence intervals) are not correct.\n\n\n## Estimating uncertainty\n\nThere are three ways to estimate the uncertainty:\n\n1. A bootstrap\n2. A sandwich estimator that only takes into account the outcome model\n3. A sandwich estimator that takes into account the propensity score model\n\nThe first option can be computationally intensive, but should get you the correct estimates.\nThe second option is computationally the easiest, but tends to overestimate the variability. There are not many current solutions in R (aside from coding it up yourself) for the third; however, the `{PSW}` package will do this. \n\n### The bootstrap\n\n1. Create a function to run your analysis once on a sample of your data\n\n\n::: {.cell}\n\n```{.r .cell-code}\nfit_ipw <- function(split, ...) {\n .df <- analysis(split)\n \n # fit propensity score model\n propensity_model <- glm(\n park_extra_magic_morning ~ park_ticket_season + park_close + park_temperature_high,\n data = seven_dwarfs_9,\n family = binomial()\n ) \n \n # calculate inverse probability weights\n .df <- propensity_model |>\n augment(type.predict = \"response\", data = .df) |>\n mutate(wts = wt_att(.fitted, \n park_extra_magic_morning,\n exposure_type = \"binary\"))\n \n # fit correctly bootstrapped ipw model\n lm(wait_minutes_posted_avg ~ park_extra_magic_morning, \n data = .df, weights = wts) |>\n tidy()\n}\n```\n:::\n\n\n\n2. Use {rsample} to bootstrap our causal effect\n\n\n::: {.cell}\n\n```{.r .cell-code}\nlibrary(rsample)\n\n# fit ipw model to bootstrapped samples\nbootstrapped_seven_dwarfs <- bootstraps(\n seven_dwarfs_9, \n times = 1000, \n apparent = TRUE\n)\n\nipw_results <- bootstrapped_seven_dwarfs |> \n mutate(boot_fits = map(splits, fit_ipw)) \n\nipw_results\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n# Bootstrap sampling with apparent sample \n# A tibble: 1,001 × 3\n splits id boot_fits \n \n 1 Bootstrap0001 \n 2 Bootstrap0002 \n 3 Bootstrap0003 \n 4 Bootstrap0004 \n 5 Bootstrap0005 \n 6 Bootstrap0006 \n 7 Bootstrap0007 \n 8 Bootstrap0008 \n 9 Bootstrap0009 \n10 Bootstrap0010 \n# ℹ 991 more rows\n```\n\n\n:::\n:::\n\n\n\nLet's look at the results.\n\n\n::: {.cell}\n\n```{.r .cell-code}\nipw_results |>\n mutate(\n estimate = map_dbl(\n boot_fits,\n \\(.fit) .fit |>\n filter(term == \"park_extra_magic_morning\") |>\n pull(estimate)\n )\n ) |>\n ggplot(aes(estimate)) +\n geom_histogram(bins = 30, fill = \"#D55E00FF\", color = \"white\", alpha = 0.8) + \n theme_minimal()\n```\n\n::: {.cell-output-display}\n![](chapter-11_files/figure-html/unnamed-chunk-8-1.png){width=672}\n:::\n:::\n\n\n\n3. Pull out the causal effect\n\n\n::: {.cell}\n\n```{.r .cell-code}\n# get t-based CIs\nboot_estimate <- int_t(ipw_results, boot_fits) |> \n filter(term == \"park_extra_magic_morning\")\nboot_estimate\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n# A tibble: 1 × 6\n term .lower .estimate .upper .alpha .method\n \n1 park_extra_m… -0.0637 6.95 11.6 0.05 studen…\n```\n\n\n:::\n:::\n\n\nWe estimate that among days that have extra magic hours, the expected impact of having extra magic hours on the average posted wait time between 9 and 10am is 7 minutes, 95% CI (-0.1, 11.6). \n\n### The outcome model sandwich\n\nThere are two ways to get the sandwich estimator. The first is to use the same weighted outcome model as above along with the `{sandwich}` package. Using the `sandwich` function, we can get the robust estimate for the variance for the parameter of interest, as shown below.\n\n\n::: {.cell}\n\n```{.r .cell-code}\nlibrary(sandwich)\nweighted_mod <- lm(wait_minutes_posted_avg ~ park_extra_magic_morning, \n data = seven_dwarfs_9_with_wt,\n weights = w_att) \n\nsandwich(weighted_mod)\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n (Intercept)\n(Intercept) 1.488\npark_extra_magic_morning -1.488\n park_extra_magic_morning\n(Intercept) -1.488\npark_extra_magic_morning 8.727\n```\n\n\n:::\n:::\n\n\nHere, our robust variance estimate is 8.727. We can then use this to construct a robust confidence interval.\n\n\n::: {.cell}\n\n```{.r .cell-code}\nrobust_var <- sandwich(weighted_mod)[2, 2]\npoint_est <- coef(weighted_mod)[2]\nlb <- point_est - 1.96 * sqrt(robust_var)\nub <- point_est + 1.96 * sqrt(robust_var)\nlb\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\npark_extra_magic_morning \n 0.4383 \n```\n\n\n:::\n\n```{.r .cell-code}\nub\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\npark_extra_magic_morning \n 12.02 \n```\n\n\n:::\n:::\n\nWe estimate that among days that have extra magic hours, the expected impact of having extra magic hours on the average posted wait time between 9 and 10am is 6.2 minutes, 95% CI (0.4, 12). \n\nAlternatively, we could fit the model using the `{survey}` package. To do this, we need to create a design object, like we did when fitting weighted tables.\n\n\n::: {.cell}\n\n```{.r .cell-code}\nlibrary(survey)\n```\n\n::: {.cell-output .cell-output-stderr}\n\n```\nLoading required package: grid\n```\n\n\n:::\n\n::: {.cell-output .cell-output-stderr}\n\n```\nLoading required package: Matrix\n```\n\n\n:::\n\n::: {.cell-output .cell-output-stderr}\n\n```\n\nAttaching package: 'Matrix'\n```\n\n\n:::\n\n::: {.cell-output .cell-output-stderr}\n\n```\nThe following objects are masked from 'package:tidyr':\n\n expand, pack, unpack\n```\n\n\n:::\n\n::: {.cell-output .cell-output-stderr}\n\n```\nLoading required package: survival\n```\n\n\n:::\n\n::: {.cell-output .cell-output-stderr}\n\n```\n\nAttaching package: 'survey'\n```\n\n\n:::\n\n::: {.cell-output .cell-output-stderr}\n\n```\nThe following object is masked from 'package:graphics':\n\n dotchart\n```\n\n\n:::\n\n```{.r .cell-code}\ndes <- svydesign(\n ids = ~1,\n weights = ~w_att,\n data = seven_dwarfs_9_with_wt\n)\n```\n:::\n\n\nThen we can use `svyglm` to fit the outcome model.\n\n\n::: {.cell}\n\n```{.r .cell-code}\nsvyglm(wait_minutes_posted_avg ~ park_extra_magic_morning, des) |>\n tidy(conf.int = TRUE)\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n# A tibble: 2 × 7\n term estimate std.error statistic p.value conf.low\n \n1 (Int… 68.7 1.22 56.2 6.14e-178 66.3 \n2 park… 6.23 2.96 2.11 3.60e- 2 0.410\n# ℹ 1 more variable: conf.high \n```\n\n\n:::\n:::\n\n### Sandwich estimator that takes into account the propensity score model\n\nThe correct sandwich estimator will also take into account the propensity score model. The `{PSW}` will allow us to do this. This package has some quirks, for example it doesn't work well with categorical variables, so we need to create dummy variables for `park_ticket_season` to pass into the model. *Actually, the code below isn't working because it seems there is a bug in the package. Stay tuned!*\n\n\n::: {.cell}\n\n```{.r .cell-code}\nlibrary(PSW)\n\nseven_dwarfs_9 <- seven_dwarfs_9 |>\n mutate(park_ticket_season_regular = ifelse(park_ticket_season == \"regular\", 1, 0),\n park_ticket_season_value = ifelse(park_ticket_season == \"value\", 1, 0)\n )\npsw(data = seven_dwarfs_9, \n form.ps = \"park_extra_magic_morning ~ park_ticket_season_regular + park_ticket_season_value + park_close + park_temperature_high\",\n weight = \"ATT\",\n wt = TRUE,\n out.var = \"wait_minutes_posted_avg\")\n```\n:::\n", "supporting": [ "chapter-11_files" ], diff --git a/_freeze/chapters/chapter-11/figure-html/unnamed-chunk-8-1.png b/_freeze/chapters/chapter-11/figure-html/unnamed-chunk-8-1.png new file mode 100644 index 0000000000000000000000000000000000000000..2e4c475966d886e8c3b5245bbdc232ae1bc12135 GIT binary patch literal 46446 zcmeFaXH-;6w>H`Yl?)1qk`)yYNlH%Lf`Wr zD_%K<^K>8Ie|z=~Z}5@(n04~ZBmTG9>TgN!@~ky&D7|TBBA+TamsRvE-=?IYcYQQ! z4msBEJljwCIRA{n%l?zZXPDFYM_PumYDcT8q+CAwHc|UF@$-egIQ`-j5zBBzS@+Bw zXH9>{O5Q5oVk%8ti72S=erJkA{qZ4l*bZZLrlCLS`OP}=j+<-bRnoRmoI{5yb6`40 zVrb_xv9F(;4KluXh9@a_r>ZCDz!=wxKIsAH`JPZ2(f6sm#~Utw5=)LyT6~%Y~=_WM^U&UuSZdh(z3(84V~uC@LJ_! z&rq>C9n=t_V<}d<5VV#*diXZ6?<~KCl~hNQWZdnC=g-VO8eTZ4ogi1GlNb{&zirH) z$53-pHoz!~%B4=ob!&OUO(*&_+0au-CQ9bulY-hgIQcq1ogD6WbjGWR>0Xyy_1XB^LFZpc+Jz6$w-_;}BAno-Dv{=#EUS}Pv|lh_0b=5fn6<8hX6*4n;J z4^;53@4)gjt#|6C#=i>5$K$!?o02N#EVa0bsgRzQ$*xW0x9@PZ3SSny{PgY`>+xsb zg^nJJ<$TKhN#j=Vt*=s7XhdsY9(TUASn##glZfj&TKrM*qjC51H9Ho!UVY2PR$oDj z({*Ayu{QNDI_!<4eC>4}6)9??i`tzhTuxXpTl2frb+<-4Veu8lV;vdBW8(T$I|q;? zep5Z#d7@JypS3A`3mpuqlon-GZ^^~+YuTD+#lF=O)!rcHXn^?^codOY>^Io@zWNFN z)aKLbQ;t)9!~Ngo1{5oSLAhXW;xyssVv)_j;e9S*~;v@>0|68BF6x)xO4aD zz@yS@)3{criOdF?py6^s;I5rIEu@CSntV~BTu*@wW8>|T?U6K>`nth@$;Az+s;oY(Xqn(m{_ zzMps(oqucZq< z-$+kJB&JdL88#iiLY8)%lk-esUGn)WoP)=7$-bYj+mzesQo4P1byUux@pYH`sIzYP z=P2hDxw7r9vBv1l66f39UU`Ba<^5qq#H3_&oQORas8bA<4i&;j5B=vY@XJ{uqWAPD zsK??N>HNiJTA7bPJx8oQ2Js9Lu?in2%&t38ROmmZ0Y)LkG(wF<8ljDVb!ThXnEZR{ zg#RW(M3X@ME|TgG%eG2rFaOtaf4qUtU+MtV@0R*RM3yb8si=d$k?8MPal%|_eka}E zzc*zQ6F0?2>yvQpCF|e!5EyzQ)D_h1?Ci9?l9JMhbsqWm(@kw{TC?d{EtTp=D?@ts z$;6MSCRGxV%EVLq-_$xrj1x}H%M)Mg(j!Y(9XP}jZGp4jipHd-hESvANxi(hHbhEo}CGu_iNj}8Q3 zTZ)C=3QfM8^Qa7YKa7Dp1G&ACVbu~p4LfJiNHts=COB*2rK(?ka{jVf z*bo$H-(N()^dH6KLhaQFeomN*mluQ59=7~RKYs*2psxN>_DKlha$14yYHo>@fZ`_W zgLMAwZxi@;2inVzmnq2*&tllCp|)B55n!CdvUcPUdwoeC0gG5#;tAb5*v}QA^Z^6s zj9V!~EutAL;$4E0Bh>0!k8{FOFRB?pE#m26I{#C$vJW9y>o8b7zgPDm2wr`h?GF=~ zn7GCYwfY<&G&WybUP3M6-9A#xpDp}9V+#=$&24SEy+Eid?=u-48U6g4UZCLn)AQav zndZv*>!XpCf`01_8MnVg2>I=7W!xF8xHI&CJ0LVPl&D9m<#u^!$_WUe!3d$^sS)1Z z-j(LzhKFN3h8{F^Nv;WXq#3d_!tq=4LyrV+*PXllnasV{QE6qloP_pX^&_T0+;Hvc zYk%wRrHOJjC{CR`-t;`Y+@U+0cWZk9@4L2+^^-d@G!kwcb04k5>Ubf4Cby>|HLst! z6u&(em{+zvGuZDvN3CP*yR7P~*>XEUH3fo?5tuuucuK^u`83C7XKlsPlgjemn@yG8 z3!k_v6=-D{PN_y6Ny{I4C}h{2P31FJvPxr&-{{C{7zh_IWNuG4HjW9h?tacvbRKF) z4fh48l1=GfnK`f-hF=St@4_O*wy^juwV?Fzt==&?@mmS2z0PX6+vBleH1wygj%PU( zGwjS)($X5_*;5OX^7FT4^X4oILdn=Nqh%7zWqGpAG`(bng2~Y3+FHF%=^?isyS!tf z;rs?K?ZQka5Xw}D;f1x~uyApmmiXRIYDmJK7FtUnY0C-rp>f;7)ds@1Ti|R93Sp--sfk{ZHEfDj6q_s|N+mp*b_1R1x8SmAt$0_&n#mf(O(n46o?GtD{&a|N`%5a-v za!i@iXXxy|9-3PO$llZFn~6 zN+{J|B0fNcJG(a6KPLCF`GrUT=@7QD*`Tg!%<4thQ8UXUSJ+61aaS7F zHrI65x+M2YypRY&6i=5tpQ#!kXMT}Gbv>pnmLW*ztX5IBEsXi&5Q&|SA9jsexaqn+ z+;idk!YDyV!T{Gi{raAp87GSz`XI{_d?#~S`x!!JAj(GcptVdeV&X%qrx?}qFQVj{ zrSN`T>5kk~&y%*E>#rkf&`QIh56X%rPp!1ojh*)lj;oV<8-k9$NHHZ%p(z(%&7|16 z&~WM+u>yIObAws8x_s`)G4%~!&J{+EmYQzo92u%5oHXiX38@@S*;07?^`Y1mKihtdkJjFuC$=N~%c`^*-M1t<#?oVejR@?!Ag0S>Cg zF)Xc`4DN!PbgAU*7#zPY=Rj(t5?R>z{w}1Ka@cZ0$6xXOYUygpU{Z1#h4biVa`ypW z_x#Zao6PxhrDvCtjkm2OzCTu|*NtQOZxj{1<)n^(4Q2tjek|IHO$xNs;(aa}Zl5+E z36>riI>1(q$S(1n_sNeCt2})30F;JuVziuuxny4_xPmwv5oI7!EoG1*_1?&2K3qb3 zAmJsHBSd#e9=KOkB0PUVc_zn!dv)>si{Bh5)E@)<6$NBvWkY0&&GLc} zGT3Mv{g(G~*qq$l+$W;hOsH!V&l^RBg9Kn?Ws1p=7xrbdeqc>LjpaA7v0PUpL=unH zx;M4c@$B2e35v%VR~VRmSJZPrQL~8M3JgYt-A7y_##vZQw2_8vtjk7x(2}zXPA3h4 zudKq{mu z)R#!{Kz@TW4ZqN~ za3NDM}^c1=#*Jf*;dUts$XiS$E7A<5dALM75lO10kD(-Ln1;-)^6~L&R0r~r(Pog3a zaX_j7xVL<+d=0|W7Wn{@st>q3njk#%>_Tv_4GSBFi13=KK;raN%nb#I1ZdMFxTg&# z2}b{CNPwUR8Qw?_C+ZN0d&%IQ!8Vl;5@Hd57Wiie|1Y>dbSx2NK$(-Vt|q(9>a#3>Lga`d}{)I|6zf2-{dP7MD~CWfI70&}Qz z)m#Z`8Q@8nOUjp5UJy2&?g%jIE>)UtP*oHF^RXE9J|_vIIybQSF`(lO?C+i7qz3g~=0yy}KO1FcWHS7)wz7sXOw@;UIKI)= zgL-$*n-AV0_mOgbbLgM{{^wmzpsYMfl_U*yuNloPaGAFE&y} ziaW5Y7XTT!;ob)-D&qZQxVESl4Io>eRrKgEKSeE@`;C$Kg_-rmv2X|Is-PJu$n*Q+ zHwWPFH$Yw=eII>*!{^T0v<1A!qU<`1xe}|PDjA0rbN}`FJr|dh$^$v6g_y$RHZ0c~ z-An1+y3E5fJa9uo_PzOX-qp(AhBHk*{v|in2T5LjK;E+T@ySUFSc*BN_1 zuVlt+gmVL!{R`d*XyM0vAx`(u!VV-Kl}L{mP)};zf2=g z`$CT7(-E?-#TG~Mi7SGGC9J`!#`wIoG^6^be z<}!h@t%Z3|j5vM6O*LuC=Z#tU=4Ya8o}E&m&`41Ivi}fA;4p_9ONonJdTg<683p(K zh;C`O{3^b?*WV}FF13efHs*QARWF_vz*#3Uqzn}p>VsM^2ME7(3#SH5Nl$aMXs7Bk zc6WCNSec$fO_C0o*=$`I`H)h0&GJPdV#KU_vMpIQK2FZGvQ-=uR*13+O(V=B_hA^X zl&w_@%5NRxK$l0~sUT-GaO*u&g17K_>kgvYfcun8w^@cs+lWHnOMa`^eY60_6q1R_}3rR_aL@#m8Q z0AOHQHOQ8Ak<`d{xovvJZ0OFWpU2HfQkwvsr9JG+kh ztGSs)qPFa8xSWYe9gQ87ox}(crNjDES>TLLK)R1EY^^T_R+TO`N>@^g_tKSa_RZ1c z4{Wavpi|*%-|imhwcqtqe$$`J^mP9AdC1Ac1o-bKyurIl_I>B;4 ziTJMgt_=n%_}sU=(I`0OZcHQu9$@$2lmW%h z!C?-(hiidgGhRMgHnfnNSHkZzdJavl#P0cB-4NPHL=y&kGt4JB+;J(K zoc`3U@O-puyWSvxn$Um^{!-Dxs}UcZO}PnKHa!3|EOlLUbaa<@0L^OY8&Oty4Ja#U zO*en+lIS<9E5FhCZ`*_D28oC4F(p?BQ$*Iv3D8_bYLJ6A1G5`13r~VFpP2oGnR!I) zx_NuW@^LOebIq(@V}l~DC;`n?-trvEwaW%*E-P%<dv+JtoBc8>P%h{t z15hzu*?JEN$YTba@@EU782HZ?{*gim5&b_*3Y+v441P&HzaYVdZt9BA z=q>?pc*Jxo&(Apf2qke5KpN`uXzU3T^tHyIsiXc=Ng@_0J@#2HVm}=Gp?Og+ z4~@@E(AfAndiT}zwzO!-Y2%d#xy#J19f3QiFSrrhCAxjNr8`cY^GlVNO|HU+zBi59 z$vk~Ud6D2FBOHNt9I{=$`nTALATW7Fe8pS%Qr$c*ho0P2pOiHu=0$^?g)uKImDDj3 zhPrCfvZ&U+>}buiVig#jAq!)Jw+iI`P=qL9!EcAwgyHx1?<(zI%LW(aCg4>cZ>2z_ z%@HP_fkLg_VqpHy1%qGig4q1Pvkx&pZ>q;4##;NVr_Vd@;R$O1?{7c-UiPQ9Hlah| z&aJZFdepYZUbHF$6C!9;3+YMB5=Xnj-;(UwGkx`MIcTJyOR9J3AQyUPGvOY&_w2yD z6vIDD!x(`Dd0m$N8PGH}6Weg;Gv{utbOo?n^{O=`ZVJHuMf}T_P4WWbcv(vyhB7Fw zgz3NN`S|O<_guVB00t$dJMXi6XZ>v2*5ft+88O9<$nO?z9& z;aIk~v(aN3uD2QN;U~TO)=A)Yq~_%h|3lxp^B9}Hq!YwVO?}aomM*(6ULI9&(xAF) z@-+l|qxJ*VXOSVIQG?pv-^KISCj1cHyI||7V_<<$0)>$P+v}L?b!zV_e=;8CBiN<= zYW|h&ceI}`n3gfXP> z!}f=t0`qbruzLtb{1F}qs*e8J!^#)-aj$zGA@?1Ol&^UINY9Uu-VEumzU@L}H2Z)0QFq@<)WpJ>T-rH-j@a=0sz(#xy4s~tLa6`zsh z_hotMG+S~)!tto?(3kQmCELVV4A!STU+jd4E~jT@r(VO?Eb{G3t~c=*t_d*(XI4Uv z^K#@7awY)Y*tM-0RbOpac5o~ek8Q}Y>Cm1X5(8&yHqh?(SuZ4z=2C)CgoU*xAF(CL z>=&b*rP+}b!DV&sCJKZZ&+$?+LyD3!NBq`;S;PeZrdTwA-BrQF0;QF#ZO%8$M{&C` zs$h|(ov!XO=Z-M%&QnR>jio5BR+W0<7wOlNRvugv9ZoK`j>S9rt$MW%EHhVj?eru| z8l~EGce@ow@I}O_Yw+eb*5tDBgLKBGY=P>V_??e8QxfwQP9!gXf1kT;kX{`6MME+H zs)h@E6sHSWRUjlNtjqr-|~Lz~w-g9AF7Z#eV=qoReLF^Ted1|-uU0#QGImyc{6_k_B0}wYc1pB z-CdCoqEAHp3CR94^znVW8fA~x=Ts-02-BCG+56kBG+#0TUc?py?_-EjQPw~&wAEWd z215%#^uH%uLhVs>zm$q5A5dr7s*qXMr%d2r5BwqpWKMu^3N6Q}`yVyid=jwE+I$I+@$8>1{9nF>2#Xx*oL!c&y|I45W3aI7 zGV#riBXAnK8Me!7~?gKtl$ZCP69DQvO z5F09W`iyP{27hW5%8GKT>!5l`n$9&!#HNwlfeav5^<_@%_5vWII`!BartgUXM<7MF zR{ikqo%xPtp!S#4EAnl|y?B_-bZ2Xw#6IgOgik_59Lvaa0v7USY+$=@;6BIBd>|Kv zobp9Mqhg60UK4R~MLT=*v*qvi)nL0VTJ1!Bn)L`Bz)CIv6*>4^LzP?&zU{J@boA&R zpW+jRi)!R+TRV!Mbet3fll20{%omkd56R6y_|~A1Lw~W5MPpP`d%M_oGFpjm{>Soq zpxgu)>MPlIcK?1q{B{?PR(7fB&RDd2z|nKJgX7sZ_j+oa;;6ROKGt^Gg_hopg`2WFp_^SPCXUjokM-Vy5sBDgQwJDqU;`~!*YJ(kzHzIo)8!8$2fA9;Wx|B zA!7w&F%3J*a413cdh8m|n*A7wZtD#V*z(KBu@P94n9T$n2?3N0G;L0t& z#@<_=>WGC3Ed3D*q)yv(?*Z+Lm%87nI;v(%(8$gwlXQL$Qv-_Ef7`Zi;(bTgn2 zocyXgU!xPo?mhiMVb+Xm&rbfMN)%6nJ=`D^gbB4SK@PvRB(1?v=~A-fM3|PJk1>|aYgXR;lKYwpbCbkKD5ed0#5 z(f|&Z+qCVD>Y|C6X>$b`(xF1qC|Bk5yu8Lj(^mst~Le};%!o)!Dha6R0 zR`LbPB$LzmJ^RfJy`%m$I{7hvJL`khOvaxhzPu9V>UFGT>*^q(m3e$WA9PmWNK%H} zl=pWLvqa-#a;pXNZ)kljo0^M?iYi$giVte1-8hO|isdlJ13uX?T{->P5P%_?+%DP$ zg-_R-e|>v#5U?uCyb!J}`bEP(C$!fpm`dqDka*Y#D3Ck2YNkN%xf1akpgcQRK`EGG zJ0u3hRn&BKvM|}5wP7Wh;0A%BGf{5cR$V6-jo0gSDztPAMTc@dXLKnLW^11Hya`l2 zGr-^-W%n9?Bjh+#b*Gw9=N_SiXE(GNj@+Iy#@1xHOf=7bQbV2v0k&)B&h&18_G^>D z1<0{dZ!@f|J;bgd-{9_sfl|A}Jy3#8&VNGF>|9JLCjg3(M{b3cCeHz>wXOcD#9>GS zSN5(_9;AhD&1bkywbLo6Q$npKgwfIx#vI5#z_++P=jRtcpM~B|Nn-}78PmPsix0yi z>o}6o{RRMm@zU;Lco}VW6--1R48F65h)jhtJ(y790%Kxhm6XKIv-t*u-7exOn-K z&m0Kfqle}=UyQv*qhH9|kMXNZ4MbUZHcX~p$DB2IkQeRnLdtB;?ei4CQzV_1c>v3H z(Ob`V98O}skP~uc4!;h#=N>o*F7P=*ly6h$-B^sbKvnvhsZ!ZSTYAMsCPPiRDsWak zvnLLdRXmu=5=ZNFdZxo~^U^G#n?MRUYk#8XaHg6MqO;i3`Izt2b3`0si)d*wIbA!= zDH3svxG~E|Z7_8+J&t@lfd#} z6W4t6t=p37wYnZDY^2-?L=X|jBry6XH4yD z1;n;v?NIwORVl0Yw+kGkr&6-yw^}zq!xnRnDorUw%1NY@uMeXom3-N1$(%WnrBAi; z*3#-4#g~u7jm;Zt{qb!AxXcRuZXNI5PJz;s2v$EIV@*obThgHe_g}nJjGl@c-pni5 zH)N)NW77BjdFSMSXaVWOm9x)}M^EZ$jEuZhU=z1zvG#t#?z=JT-lh)7p8f{EEw|h; zp^9hT7)aMR>Pw^HnO04$r_zH>e)n=_$kHw*f)n@$pvY)YX=k+-Yu>n{J&rwuc zvW`X0XP~!VM?|DPUORE38XnP>1&*+@kLBVW^ORync*DXxxAVfq+i{8pJ)zt?xB?Xl zQCH2Kvrp?J+l+9v-!vN+ramQ)C>sonk){kJ(t^Ke=IGJjz^wzXO?w9s&S6w;`n6Ni68BIF52RdKDfy8 zk`pan_*swry_gRiK<^HA+oV(L+;F^j_wza-@kcbWG-RV{O*k#pG%L5c;EAq533KjV~zmyOkT?e zBf#mlhFd1F&E81D34MdSo7E4eOU`eKKb#sq+sdcBndb;`as!6CN$Nz4P$n`J9vSVi z*_w~zM*Ezt!Z_oGs-<0#06$_eMI*QNZ6S9p<_&!PlUgRN{YIUV+;t-khpN>V`9RTR zy2jCJPG;SehBErCS{VEA@-ompW$;5SJEfZJvS!fSCv7jYnZnTAfsdg;i`c4A@uPiL zy5tNvkZXf<#>be$of{<+iusYpeSC7pe;)mBm>LdYvClu|eQ|!S|8;dD8Fp`NeuF0I^=Q7oE+(9vbutvDNw#Au8} zehmLFo1aj4?z7a<(OHWuTaToJe;WsZWgqWWFjp{Izyo-(=5w*@BZ7VfX?F$K&LH2K z%D7R%hq+PJn3V{dv`+$N0~ToCOglp4p_XT3bl2$mGM2HB^@Pydiq+>Qd3KTq2 zl+~&3*oZG=kce?4z7GIYkBql&9%ipGt=5reCn&*Wp5eGb(1e5v3&;C-EPYj~91erC zUoo|!qKB`12x~LES#giN>*TCOw5ztq`bcG$Mlq!R`1cCg77~+Ak`Uy8m;3ByLZIDW`q!e4{WsOW$OGz<`Hu&};pQbv9~C^QRpC&ao}G(GBeDb$zX#}+bM9*;E!?LnGtPI|`>y^!>)#{h zkHZQ7zXijr8AUjDWtdh5*#vV0GlT=H5t1Lm1{v}SMg4fDY8AMvLeirEs0Q!-(9kmn zaADIcbKCuMFFFc#i*|cmI2cNQpA$90OqkFbNV_F2+=y8_MOWC?mmtyg97vb(H-R)o zWo(cW6}nD{v%&sKfbBjB90UPY-=(jE8&m1DIG|qyr#jMU!wuc^*1jEpAMx_`R`6LK zJ8BX2)ldAth}(0@=lV0UCP$^kfkeqxhLMV4%@c`XBDzdP0ELY~XB>^g79#|RiU{I19qe+EMBfKCu{~tojHj>LIfOATQUaR8H zq;vTaAv!o1$kAB3(NO@}kIbh!`10j%m6QcyIA zs~f>+1_$T`$dWqY>E7*J1TlX>anjkXt!hE!H6(4*0u6M=iBkzxLEHp-q6qpBIDF2l zqiFw`RnS}sTM24V6-w}puNhv;lSHc@pwk;UWDJM=jr_Kk2SIxz+PpqD0}lX0%u++9 zpUF%};QoPNZi!d(5r1_; z!+}axO>X`Byc&K#jtT$`2^&l!-g`#@y$uPuBy_j~%mI>J2yJYa1)2>+2mBw4u#E7$l51dKuo4r4!+_>OPXqvXzR?XF1g*Ff*%`K*2Kg-~UPlekC1Mu+kOD$@&1}0l{?+5`=7ZU6t z^MAL*0~nAqG71XMOgAo*WA*(09bYn2(R{odip-2jX10TQSc4;lF? zK-^eAg&z;w;gQ)Gj~xW&{!F^_T}YaHk}<42ZCEqHa7W#d51PXj095&zH6MUCuXuO) znY_;n%PlXgcRxMBbCIceLkl^j=@aU@Fy=6Z0ZNh5igQRkNWH4=fRYzLrLS_bLhn8a z0XB5~17vwIo(_OPszsPdp`z;lJ05hQ$!`7^@U_VjMcv@X=v9;_S{N*ZTy|5NLaXfOr)$NeScFbI$PWX+imJP4;BXA2Nk4=1GX@8)r%-;2Do&?a2hbn)kxYhko z0M*u%Wj>5|m{To3lUQFli0;2V@cO8tSm2X`bRUV1U%hac=wKzcqK`+=`3qbZ{=(CC zyvj!f9nw5*y><28I{e}BLx_;$H>rq;&tEiUxOHvC@Fl9NBW5;T&8n}oL{%_vGsa*! zdF=H>VNR46b}qHlcXLKVYuN9bkoDKHbGinpM+J@`6;$mO}e8x%Z_(~xLilZ!Lg*}X*C^;Z;mF#uSE zD1w5(fvFN(78V@bTYLKT@>m9S7~4rq-`=(DdOnwQ5JVz*kWnXY*56%U>^Pz^_f|h$ zrDogw^q6LwRScu<^*Uo$DqI*>TY9bRwQW`*{0(ik4v>PD*U->_8-&fgeL&BxO+z=F z7@C-QGs@zvnp`#Ul{5XX6}cZriY<&~eMkKSwik@2$RM8&{on0fBJ*^A!tFSpzm;U&`4oy^&77RU6%>S50B+q#^7 zo}Swl!kCff@|%#3=Z4tsO*^3EDQ%HEiFfdlwSmCZtv@;MC}OFhfj#(k7?rVE0H@)9 zmJm|=hV@6X?8S{?KDcdO~PV57j0)L8EO@oz`o zfapkx%t9P$O3lD%tH70xL*epG?g2~5buWu$@7s1AjR|03La%GZ=8n(bZodcu&-P-+1S#Bc)KkxhUlxl(u@yw&R`Irp2^9 zV;*0g#Rudb?OE2dCnCIP3B)BcHaE$4U;4Lnj_Z8;)X+$2Gq*8->p{45 z{)kEt2Z`BEl0rqE<+-FX`%0B90m^1@p+S(AkPatg{3}62gvzv{S!Pv1r$it~m($Z&1;5t5OjqS z=b`FH$;J^7jw}L0LmVgIZ{*sh)m4Wjz(cnB2DT0F!w;VaCF|CC@JWgQ33(BFcR(YZ z_OX3Ux?3_JdCk}GW!G)pWlX~UK7NTvFOk}UkNDF z=m?5I5fE{cIF6lYBTmql*BjOhDmNP{4U1Sfnm^rWpueA;06T8#?q)vBum3#NxWvoV zZNi6D()wpIW$z3K?SZZ{ASILn77&o81Kx29GA9HCw;OlINIkN6kBE-DDc zK+x5G&irEQ=;)|z*JbW#hES==_w_4Wue(6?4ZwgFFS)TT%K>E5F?P5dL;VGllzl)+ zMcbH#+*nu{-qed?UYF(H@%e-=KU`)pwCPfk(AXFw{Ekk;;h9%;>uP@VizcbU-i@^M z6!x;$X^T8)0Q-?DRG6gT4_Pec)05P=RLxs+O}2V5G3#{Wr0cMxF`G^>@BPK%h|?om zAr>(KPASt-NfnRMBO@XLGN)eIJ*$^&TW;5U{Zuj`sb6Mz?mf=EQ!(OZ@z=~J@C4DP z?s{gaJEcKb?@xYko@AxiRp}ig20u9SC(TLXaT62c?Ma;{k!v7xK>l_%N+jb z9Bm;$8w$D4_R?_#EP?NiSl|h*JML`GH^vRJRH==6Iz}nrbZvu;b?+%U^9Pm0fI`^ayb<*RT=i(@@fBRWMLa{aoe^YCni(~!LD^LoLTi@D21H2nHhO-Qm zI{=L2F{9F5K3|KNdV~i>CM+r)k>4CeZjC9bL8NCrT_xYLJzlS#<8i} z2>=UTwvM0}9tb4;_~$6etlE)9T=V1MP-y{{eaw_^0G8rY5GBUW&d#(u+e$W6+DH&& zR+#kl^}QH;Yo|IxHP}&Gf)2RN+mo${BNO>Cp+yMf<+ksip4D!5=17_Ns5K_6$)0Og zmRk_B#^p!_c$nT=p3=`Yek4|NE`IW$#9g+$0?CjqS3|7+F;!gD!;jRXy5z`U3$NOy z&E+E;r30hM*qr?H!HHSUx581&-u}#Iu8u?{Nqxyu&etF$EmV#@S-oo{878cawrCSx z?RM!ksyQN6wA1E0iJYQRF~q$yX=z<{K8Zec^_6qi9p(mUo#D?>U~ zLU`$vl*Yn8)((&-!41T3ht;z!TlRZ1BJe)5RpR}i^NirucWyX2F0l5sMa)qe0N0kY z&}Ll+NJ{S&y%<3W&;%yy2?oHp)mH{hPked+xGJ4p!^b(OHgsF$;1NDlF^stk9%W+| zm(s^WH_$Vz`F!yV68mocO=>|w5f6*6Mi4=^LM*PfRfDuf$i+~fhxR2MV^cMYHfPjF zG@{gw#g8Vvx``5d@XBi&ckKSFh2jxov#Yr4lU{0VW16SuG}?sYoMLL<>McK?mTYr1 zH$YONUSX{pCf90hkgVy$E0U_K1{SxU%J1AnE+lTelw8kLTw|uxxVRqGPS$T8 zEZm=VT+m+?;0KadpMO1%Rsc3siV|S?=;b!Ph6Ri?ga-nE`mOs-07LA3VXPJarqMG@ zb8n6+X?UIEI2arKn{xX*^U-<%v$U@X;K2gH82~PDc!9(;A=-@s4eGa*pcpK}Im@VX zkBX z3+kv_wU}z`k!8@NiKTQo)!{zCL*x&r5}?Ul3neSv$pg!WfAh|MT4>Pq_NQI>QJ@Bm zf~e7(y`eos?%i#+=|@x)|M-{29s{P(OfbHX`Og zOE4+<)Jy;-p>*(SiQNBvz%H+^BGbGfbD7YFDu3jqLFVTD((M@sS^|(_We%WqXgiGk zd~<)a4-4U8%Q(*mK~lOSOsM9=h+&QH;mWWDUb&;T7+$}dx;Hq`D@=3r%z}5zz(ymwku#6vFmHNU$IKtupu=0rfq`~dm9 zz4)jCe56OQHRi_#bIPI^JB$yv!<}^cNFod74w}~R(Ik7b67R*bAYv(==YntB-2)A> z0}S1|PHyTvTFL(5(@v13_>Je8E<5bF*gW%66m|TVMu7pnAiyOqD@o3rrAxgpksgrq1P2q~!PME8bsyKhBAfFZ3$FKYfi zZiIwaH@c96j?kD3AM211o0oq1r-R*$@w@Qq;w-Z4>ImfVRuBr zJ4NZ#L6@$#r=SRSNdr`tz7z*N-21@KX``6Qc>LR8*XW^^@eb6vB5oU9fUK}nlQJM| z{NIy92AvF+D(B8z`=PW@Un#G?5;9#1vHbn3tIYaS{3JIZJ=|LfoCqq2}0@V}{xvIC)`{j<9L|Gq)kTahHN(XqExd!cUHg3WY^ydVa}T9|7VA)UMCd_~NNMY_6f@tOGjR3z=jeXiL{n?)^9yJG z`qvAG2{q6kJ!cS!ghs$(DLY?W<0i8jf49-7R?-{9Xnn;h$g{{_2t1^icnpe%uG53W z1wYpiltYZU4b~zW+e!x|3W?+(aq;pX?%!)9P*eRJ5Wa|O)jk0Zcr ztdAa~g<8KTAYAIvxB#eSTv7q5ASXB|5NiDqfN<@>*A<|a@ea%;;?Fkz*~VSC!X9<- z|A#apY`*9>tA8OCv+HJPPct;gS_x=r0N?up)1!3j+)L?`^cvxy4DM`tX6T`0P>Xs`V_w*1)oR6?S*ARWWZs%5YKc)A z{081J9lR1Fit8O9&9a#vZ90uzEXFEqcXxI9s)Cr!Ayo>iaC|QMi`LvgA~;!Arc~`y zyO^5Qu{x}dy{m7p!`9Ir`oR;Ffn#_X+*I+6_HiQfy z&R8}cWE;f5CsOb0YKaq~;)EEl(>78&BcMx=q?zlh!aNqg{M_95Rk6scSGC$O9s(|* zD6Ji)fj27!8zk4ZrsmaGSsI?DZlt#o2+J>vonnac8Z(hG>jpTHk@=Fao@L>p%(Q}n zlMkGe5)-|@E1VwlR9?0WzR)inncV;E`g;GAEVE#YLv&w)?#K<6x!04cjV~`WE{-Ol zg%iK07g$$a`9W6Ue$S~n58?R9YZ%;^vnz75$08-X`c*_yDM#&w2T88+a>>hjZ;OZM zgZjyGTm;{w3A^FwG7exN`}=N}`{Qw$o5ON=6x^}EZPfHK+> z>HLs|SO&`&5!X8YDtXNk%4r3f_5fhPi-CTwOi=$<0EF1y7s0K1)QN7nc)_vmeD`A* zs)EeLm_+N9zyukTDkt4JVM6-;-0zMP8QuHD@_<@AUCcHP&zAA!N6Y`ONA>lURxXSY zvCLXjhjJrjvaiOpW0@b3(eXS2bJ07SivC@#n?#Z3rg_Sx6bxICu=C4x@*h4wtl&}W zH7PTQuGygL=0N*62t9KyKB|Oke)nYyQn(BZh!ceZYOmeEDj|_c@-@u`v?OpXX6tI znO}VW_l3F(px#>xT|C(d9>Vw6NMjEj`uD$F0Lqd{IC%pyKsW`oI6Vr6GNqJRjHRrs zshwk4L9Wxx81yW5HFcj|$a~Y3j&y#tlNqOHz~|?xqFb_6c%^7L;&uMwtm$S4`g(e$ zw7|dK*raUo{Y9V=wdZ95&$YunEOb``jYIM)L7&CM>7Mp>?bjT$pR>ca&CSDA)ZmBO zw`Sw(r@X$nWV^4X8VO7M$9Vi8*0!f~{-rDSTrd!V7!&d~HYjCO2ozNj|7R#Gp)Ldn z_t!~H$WKd*1L)UGt`L;+Lwp3hMfp-4y77boCda#-m57*^qSvqAK%sV&+Ld;udTzt# zRsQy)DfU3vZ??w7cHhogpWwZ=vLY*6+OzlC69n`Rx%o`h24=*mkWmo~LVj2dXws1} zvd*ngTeIvM78~oU+-FAjdyEuZg4pb4s(3jTb7YgmZkWtqbT&VPKjoB=2-tbrFQ2sc z&}-M+kcnV`O(8j@QU`UD3It;=ma4f>p&^Ed5XuW#$wP+nmOw8kIA=m;f~fzW_RcgM z>c0Q`Wl5zhkqV*BF3J*Oh$1CZsU%yHsAQ>R?9Ew`#&%XB+mIsJmr&NMXCy*dn+e&s z#?B1$|9qXNuCA*&{txcszMtI3)gvi>^ZTx!_4O8ySt>KUU zy+z}n+t*eTCTo>`crb%Af@?E)oc;BYb__#7%qS0jzTxh#0g&aVm(;^l`u1LwVzhMl zkHflJ`AGC4qkr@Pms3?i%!JXxs=o~*us%#-+kaYPXq;h2SnFd|S{cT>R>6sOt>pO{ z1}EAId@Ojt9RFui8RKjl!QhE8d`!o9L_a@!TAI7(T*8wlPXI3TL`6|iaUW`!ozdmN?&h+`*eRl(o4Gjr>YWm8#`kH*-7xh3q zQ_JLJnPu4y7tZXrdjOZOeB^6pV|OZ=6w4$kq=( zSI6}ble6|{!vJ>&Zf$m&*KReZ_hL|(!7OxjUW$WV28vbY&;A-k#XnoEws40oYv(nc z@z$lUVr=nx18e>3nZ>;{%KlC3Y7(J&@#d{tm#F1bywbj)MXmlXt?fo@Vb}G$#t2wUo9@5KM(oTWSZx{8#eCxX@Z43%t zp@$eYQlRDyKsj`xc0wS4oXweuRJlce%4oqs^4mZaJKF zroPs@iF+qN$8K<1{kzh9%rv~t4mou|gk^U1;n{>XV-_hvzQB9%5sSuc_W-cY`h6sv z)oep(5ibU6M*?Gd#)mTp5`9Y_5X~PR^xdfW=iJ-9q?{8bFp?WuOU+-F052AaOrT;^-hvuXc-H# zb}-)zYoVAEoxrI+Q|#VbGr_iVD*uxbpFL2Jpqj=~z$Bn{FvTSI()u*>!neI_w)XZM z(3V@%7L{EiWVt9nCx;7LDXh;iAV%ryACtKjlKNqPZ+3la(Tnl7=JMC%bCLwzrVn^; zjURoJJ?@qVO$t#VX5r!Ck`rZ7PYWE!s&_X|(`xF~T-|z{3oo_x7QDRP@5j!!>xMzk zf!F%s;rjZuJa2VB)p-I0EihncT<+?$qDHfiTlBG()cxM_7IQ|H4Fn7Gtg||8^^Xsw zy&P7x&Rux;mqks}C=t)oI~sOSf6K`UA_*HFZc?cjPg5i8nZz`Qi_r)A%uRE@c6tl6 zHK+P#IY)`w7pGsRK@rs!(Nl$~$(#G0>g_MgzFscDZh=XiOA^d$);GL5Q`hxn%Z*3( zC8Dn$tNkiW?%PTo$G)+a$WYv;YIU^^aTQ&dfmj0Yo;GF zH0U3>10XbTHC*A%!>Ex!H3>U23wnG1v9MgjUGYz;#k=)8GN^qw&`;|d=Na|fH{_Gq za%*vXLfYXuU!Td1=CMQt zq12>=$}dJcMuqRM-&k<|x69R!y?iejkr~S41dVYMy|wHqDPO(U2^mQ3#tZwL(QUS9 zv7*QuG$zbwG)N0~<^jFHGTV-!980Cv`HR+#8TX&BZ1WFwtV$@Zs%d1CbnLAW)9){` zv&u7zPjISEd$Do=KqH+D<8Se`#nRQcXSJzduOgW5R~eunS6y3b0#JABNI6W9Q!TMqRLT zHBF4feumDRQ(e_jC;cSokxw=5KQjFAp~UdY9N|wTEq#z9_A*HZQJcH7`|EL73Dwhm zvLR!IgJE4Bl1}Rty(#K4q{Wn(x{ocP=djuViC?O=CKvKrJ)rbDy||7!UB%nyD_|s8 z{sro_TxthZpl_&^W=pS$y0E`2;Bt1AQ1#-}D~a{Haz5F~hZFVI#*{q{g3+M5bKKOp z=P_sdK$Rkpd4~MBY}ouu8g7_&XUcIVR~;X_`g?In#f@tU-qM?^GwclQ@dYiTYEm<< za$gE9-dwkSaTbw4GBA)1E5Cf57{9|;^J&i&n1a~zMO>3o=*rRrOAwF4i!BN1%*cGLu$awT? z3%^oRxfWL2&OpIp1)hKH=9Q{UAGA4CcA7fAW3_zu2kR0I17Mk!XMpIU4IC9(dL4%9 zA9FhYr?LIZ7ss!z+@Yx}{uZ$W5#1pF`vQXFP&yB1h_rPi$E#_b%x7R%zC-~jalKq9 zz~rJo{q4$w>+3MK#Ia)X`u6R=8dJ*;zOoZiHaiY8u2_r)Pryxuo14*$DQ-q|Q3G3>aA2&m_l{ggcEsW*+^#_(e_F-@V6h4)&CNNsbpHR4$bex^hg4_vCN4V{U zyNrD3l72Ykbsb)^GJ1^-s8Ls8Pbs4fd47^9CC)($J#ce)2^9AUq~>)!_<;aWV8Hmm zH}(qx^ve<`nc7!NP>|yA!Sjfd`;4%gUm+Zn97TtLx^~Cj3JLKIN;pn^Hu81w=oaYY zA2F%FD^1mbwieTaaLeAU4Rrmm5ToM;@%2Q_0jB1=RQ4-D~ijh5yL6 zYN@n!7PPKSSz!r)-$XXNJs?LMfo$cA>M@6N$zJO}3%}bSa@-yd<#QV4Gq2yp13?PB z?^PE8qUyXc4=D!D^$~?3Usnc{l2!h{eu$z3q-TmKI_Dj~F|8v-g5g4(*OG?R*}T*_Hjl%FZ9B{cNXZ%p?;6 zSYAq}FH&uGv{IP?O6BBDKwgCASJo%87SW!fgo?U%qQuTsjlc`RpJMTukLK&H)hLT{ zs@#H@DxPoUKpq2ABQ5R&1F>iPC;DUG-%XK`h>_$OA6j&b{f@5~TFJQ=Cf{AJm9gkLPb{!4lJ?GbKb`KXQf^ zx!3xL^27&u)UMD*_c1%_=_3$a32MC?xHo@vKetYW&YV-{;J4E;vJ7^Mjic|5<4Sn< z(KliO3vXchpf0{87CVFL6S-vN&%VBL8j;Va$+>9nCax(zsw4b&wf#BQ`zC_NI?CN65x|*)4C^`ix;PUv^5vqI?e(k-vekRJPV7@L^NO7h}nt&;s9zkW^vno zmjZ=aG@SVURrn3;YJjaZ9w~kX$GH|HhB4{f=95ozJ0!9N4~$mpjCE&L zL_t0==UCCTz3Ngm zgHm!$T{ixE>88E%rk*!?`&m}6la`B|$!{+RY#a(wY*O8L}wl51YgCWa18%2Hr4L=PP_Ikl|)q zCs*u10>M?b4z#*fKi}^7Z#C_O0o{k|f30cvW3rXh;{ke18Q+gB_AM80!O>)iKRaB2 zHZ`1Bp$imMFZ-H=+2jXUuklBykLAPeM?stv&(BHbA!*aEML-l(_V9M+JDRVUpkTOG zi@>>hi|vCMy1(J#0r!{@X^0DHu4Qui4(97nLy5bg^gG+MQ|#gcQEO# zf%%@;ZsYx~r{WWiuEqm1%Wo!9EPPsUE}E>J_KCXFi;1OL?3=vl-1XEdmg?Pd-itJU z5llzJDK(U0kWFwqC?bKU`5yKtb@#Y(9e`Isvq=1Rsquc`J+PifGd{ATVH8I=;)_^M z4Ina=mEL2g-kpW0l1)t)Gig;m-HEtW1};=hEhbFauE*_&ry)ZFw8=tQ_SO!JG@1Q{nUNy>m5(x5G-^~h4o-9d2K7wv-f za$y}&^a*w&+O!YyA3D1TvZMQ&A{>fTOp2P#-KX+NTJS#<{Q+P1% zw@I7iNYmPj7ff&V(0u1-V1mp+ieXmc~~0wa5W~Rhpe-&(>!HP zS*#VNmrf(FzRA9Ni_V78r~Pn>#Fbn%p9oKGhxm{g?i-%XjSIC%Z42ZIXhghzp`ilw z`f+P!*EQR2;gMa*#Jon@ICCe0`aVgpbz3`S>L%KDOJm> z1Cm|$^j-vV9X9$_ugzI2wE23S`pjJr)819yb%2KF(LGT?O)3}Oi+ZEpZskHmkXI~S z!gdJJmX7Fcz+P2f(!MTUx&%**IH@j7HkQGrF@GIujD{x~4Ns^%H;+htaS70|N{wxz zHj-eHY>Is=HO>h;H=T%m@nTm8t~;QUQYFj@(CDi;GL_Wq&HO~EK-17@DnF~!Jr-id zJ7v=$asbZMrn^MwhiEO-v?`(a)<%#fO6iSr~A4XHA@IX^gg6p3$bdQ!3n2H2K zG`ixT41I391w=1`lQn5O3wf_J03-)Y8olgK- z-#@}XBvtT`R|&luyM!)rb|;(Pw*@aSw`GXItkgZ;YD zDK@G%HUYSvqah2%qfBp(fL!0|PYblR`vSrD!?*2Nx zN21n4CME)7-D6@Hil5{AU7>e@aKh1|`>&RZReZJu*xSrA6~AOxJXsJ8UdgVd-sTPK z&(+X_s=QUklS6-g=U-nOkvPv{E76$R_|LcS89e>N!|RZ-I1A^E0Nk>v>17=I?;_?k zvbEvzbB$mNq>{km>V#nXU%Sn`TWieux{85Wujn3p!kKxA5)RBfNrI3!jbca6JC7PD zA&M0n7q^2p-(37L$1h_T&fAa?VZeYQE%ZIm+x55N+c-{d0fe-Jr#ca|xGc$h>J`RG z1B|p*r3`PxjZE{JN-8^QKK}4XK65Ryg*xFDs$Hz;-ffyECk<*ELu6pQx(`o-6pBv^ z5tpodsJVYb0xLV;l06iWw&;*O{rTF{<;)1LPRx+@LqC{ zO7nHwKg)uc_Mg&=9#oP={^*IeIW!Vc84DgDdxH7gr#$C12JqgAT1g9CD~rqJ%aIDm z#OmC><@rG2XLQ2f{98VG!5-J%PU3^o!O(yvWNB|f9h~*2vp?wUuKd&A0aqqa>Wi6- z{qeNUBfe_oKKqWa6e{9F%8W@qyM%SNK%_-rtMbDg9OC&EUb}Q8mmD1s=`T=uBHD?p z+@@dv{7Wm`CxZFo+Nw`P?WiuNK+Jrn?>IOR;-wTT~11@UQ@(fHZILR{#Dz zR|Pw_B%s!90mbmbm-DhDD}_M0fJS!N0=ePEfG);aEoS7Ivc3gU}56s7rB@-4lx(l4AHCQvj8R^vvmZ?87*UV`BOd zRJSNE-hxgMrf29{$~&l;R}!&iaeF3Og)Ao6#G(D#u(jvMTs7hZFX7ae`gcu8<4-)C zh-q_+P8E<4pe7639*`MkyV8Dr`6xI6hdkk--oo;&+N39f(es3KS?ebY$^)Z|ahtmj zS1vAA{Ba)4h{8U@DNlKNufw8fKK4>`&#-Z9_Rz^>yR6QRMtKJpWsB}dgVYG{*G?O& z)juL77(|_GjURfrQ0*Z&*EaPd-tN@U&wM{JGbNd z(djYNt9v&6{DE73b8EyRxdPS3|MJqk%+Z5uZ;O9huTq~w*+5&2PZevP$<3?E97-$3 z-25t_DDd75S$$x@6<67Q5)$B~9FiIK#fe8OA2a-C;8M3V6mwi2HGU=}t8t&n59^y2 zZJ^TgI33PSA6YL7q4QNzYIU@~P_v8rcX1CSAF?O|ei!g$$!qoxr)y88Dn2EGHq^1U zfv<;~?*x^?%uXXql9Z$+sN(k5Hs6agTLWecwvL{XqIL>3PcJVe9WWlQJ-72iAqDik z;FuDF{UfLR1_j8oIp4hUaL<+m z6x4|wN*r8b07G8Oa~gZ?fd2$E?{iQ{F-{K_2%+U;a7?!r8paaKFRh%+F4Z6>4AB zfW+%?!+D_d{;~4%0pYh+8m>hBa0Lm%C}?kRK(Pa}jvG>Nfm3-OXs~qlMmT)TA+Q!N z!4PgUbF6l7Z*QOfI%6&CZan+J24iSEQ)981x!1yzUQ5I4g@iaq`_0h&4LvqDfF=XU zuBe#a_171VK~Ov@fS}p6?CH5Zp38NCyvDFD{%ZZaOzkL9UR{1Pk311Qb5VD`*1?Ns z5zf0pa$G8K>LOixpKc@ZOo(aB>B0S~o;lLh__g~*7cwO>#hz6u%}2aqzn$jxIY~0& zdqs2JdC1m;^hkWLBNeEyTOla-Kmu`e-}71C*|VEJCJcfkg0(CCf(an2`BOZdO?sEd zc?LWGB36UtY)e}u@Els6&*?w{SBH^^k`5q#WBqaL&@@;gC# z%knBU)t*M=&Tk}Gk&_7dvIRf1iu9(Pr-xfZn&A$+ zl{XL{k!iB>JG<^xVs1CoXD+Sxn3oVrK2)3RwY&o>PoUBWqtN~Q)GXE{| znJ1z%>bDSx8IGG0P2k#0G_}C9Sq8ld(-x~LL-<@ zJ{w*`ar$nBso;K2(HnGvtHifdfb_%1akkUf7o0FD@^gC&4XB2=r*fnBVpR|}GYFd& z6grO_WzLDdI!@PZIIwVHUr%9k`{;1%H}4?f(i&-dX8}f`F*Jna$f}j%UZ*KsT=r;= zAIFbk16fC6?3dqbwSQ1a5**cB1Q#gAzXFM|8KwZvL)EhFBygF_;HKbRoY0`lmp`)} z8vx|{2c3hJ)N_^Ge`3`r2a#&tt$|svV7HGSB?G?INm zpB49vjRw_tsZCzH%B$k&f|ZjJMp`oPRzTbZ@S|OPK*sJ;yD41%>NDYT0V4+AH>I}a#}_m8Vk-!%Jdqp_TnnMIW>m8~|y7nBrX3~)u z2tIUz(9wFJ^L@~{Q`t3eJRX-g2UvctF6;5np)H)&ud`W~~*q9n$G~=b5~`Lx24rF$f$d vpzNx(Y#0CfkLMl18T{`L{(mqyc)D + filter(wait_hour == 9) + +m <- matchit( + park_extra_magic_morning ~ park_ticket_season + park_close + park_temperature_high, + data = seven_dwarfs_9 +) +matched_data <- get_matches(m) ``` +We can then fit the outcome model on this data. For this analysis, we are interested in the impact of extra magic morning hours on the average posted wait time between 9 and 10am. The linear model below will estimate this in the matched cohort. + +```{r} +lm(wait_minutes_posted_avg ~ park_extra_magic_morning, data = matched_data) |> + tidy(conf.int = TRUE) +``` + +Recall that by default `{MatchIt}` estimates an average treatment effect among the treated. This means among days that have extra magic hours, the expected impact of having extra magic hours on the average posted wait time between 9 and 10am is 7.9 minutes (95% CI: 1.2-14.5). + ## Using weights in outcome models +Now let's use propensity score weights to estimate this same estimand. We will use the ATT weights so the analysis matches that for matching above. + +```{r} +library(propensity) + +seven_dwarfs_9_with_ps <- + glm( + park_extra_magic_morning ~ park_ticket_season + park_close + park_temperature_high, + data = seven_dwarfs_9, + family = binomial() + ) |> + augment(type.predict = "response", data = seven_dwarfs_9) +seven_dwarfs_9_with_wt <- seven_dwarfs_9_with_ps |> + mutate(w_att = wt_att(.fitted, park_extra_magic_morning)) +``` + +We can fit a *weighted* outcome model by using the `weights` argument. + +```{r} +lm(wait_minutes_posted_avg ~ park_extra_magic_morning, + data = seven_dwarfs_9_with_wt, + weights = w_att) |> + tidy() +``` + +Using weighting, we estimate that among days that have extra magic hours, the expected impact of having extra magic hours on the average posted wait time between 9 and 10am is 6.2 minutes. +While this approach will get us the desired estimate for the point estimate, the default output using the `lm` function for the uncertainty (the standard errors and confidence intervals) are not correct. + + ## Estimating uncertainty + +There are three ways to estimate the uncertainty: + +1. A bootstrap +2. A sandwich estimator that only takes into account the outcome model +3. A sandwich estimator that takes into account the propensity score model + +The first option can be computationally intensive, but should get you the correct estimates. +The second option is computationally the easiest, but tends to overestimate the variability. There are not many current solutions in R (aside from coding it up yourself) for the third; however, the `{PSW}` package will do this. + +### The bootstrap + +1. Create a function to run your analysis once on a sample of your data + +```{r} +fit_ipw <- function(split, ...) { + .df <- analysis(split) + + # fit propensity score model + propensity_model <- glm( + park_extra_magic_morning ~ park_ticket_season + park_close + park_temperature_high, + data = seven_dwarfs_9, + family = binomial() + ) + + # calculate inverse probability weights + .df <- propensity_model |> + augment(type.predict = "response", data = .df) |> + mutate(wts = wt_att(.fitted, + park_extra_magic_morning, + exposure_type = "binary")) + + # fit correctly bootstrapped ipw model + lm(wait_minutes_posted_avg ~ park_extra_magic_morning, + data = .df, weights = wts) |> + tidy() +} +``` + + +2. Use {rsample} to bootstrap our causal effect + +```{r} +library(rsample) + +# fit ipw model to bootstrapped samples +bootstrapped_seven_dwarfs <- bootstraps( + seven_dwarfs_9, + times = 1000, + apparent = TRUE +) + +ipw_results <- bootstrapped_seven_dwarfs |> + mutate(boot_fits = map(splits, fit_ipw)) + +ipw_results +``` + + +Let's look at the results. + +```{r} +ipw_results |> + mutate( + estimate = map_dbl( + boot_fits, + \(.fit) .fit |> + filter(term == "park_extra_magic_morning") |> + pull(estimate) + ) + ) |> + ggplot(aes(estimate)) + + geom_histogram(bins = 30, fill = "#D55E00FF", color = "white", alpha = 0.8) + + theme_minimal() +``` + + +3. Pull out the causal effect + +```{r} +# get t-based CIs +boot_estimate <- int_t(ipw_results, boot_fits) |> + filter(term == "park_extra_magic_morning") +boot_estimate +``` + +We estimate that among days that have extra magic hours, the expected impact of having extra magic hours on the average posted wait time between 9 and 10am is `r round(boot_estimate$.estimate, 1)` minutes, 95% CI (`r round(boot_estimate$.lower, 1)`, `r round(boot_estimate$.upper, 1)`). + +### The outcome model sandwich + +There are two ways to get the sandwich estimator. The first is to use the same weighted outcome model as above along with the `{sandwich}` package. Using the `sandwich` function, we can get the robust estimate for the variance for the parameter of interest, as shown below. + +```{r} +library(sandwich) +weighted_mod <- lm(wait_minutes_posted_avg ~ park_extra_magic_morning, + data = seven_dwarfs_9_with_wt, + weights = w_att) + +sandwich(weighted_mod) +``` + +Here, our robust variance estimate is `r round(sandwich(weighted_mod)[2,2], 3)`. We can then use this to construct a robust confidence interval. + +```{r} +robust_var <- sandwich(weighted_mod)[2, 2] +point_est <- coef(weighted_mod)[2] +lb <- point_est - 1.96 * sqrt(robust_var) +ub <- point_est + 1.96 * sqrt(robust_var) +lb +ub +``` +We estimate that among days that have extra magic hours, the expected impact of having extra magic hours on the average posted wait time between 9 and 10am is `r round(point_est, 1)` minutes, 95% CI (`r round(lb, 1)`, `r round(ub, 1)`). + +Alternatively, we could fit the model using the `{survey}` package. To do this, we need to create a design object, like we did when fitting weighted tables. + +```{r} +library(survey) + +des <- svydesign( + ids = ~1, + weights = ~w_att, + data = seven_dwarfs_9_with_wt +) +``` + +Then we can use `svyglm` to fit the outcome model. + +```{r} +svyglm(wait_minutes_posted_avg ~ park_extra_magic_morning, des) |> + tidy(conf.int = TRUE) +``` +### Sandwich estimator that takes into account the propensity score model + +The correct sandwich estimator will also take into account the propensity score model. The `{PSW}` will allow us to do this. This package has some quirks, for example it doesn't work well with categorical variables, so we need to create dummy variables for `park_ticket_season` to pass into the model. *Actually, the code below isn't working because it seems there is a bug in the package. Stay tuned!* + +```{r} +#| eval: false +library(PSW) + +seven_dwarfs_9 <- seven_dwarfs_9 |> + mutate(park_ticket_season_regular = ifelse(park_ticket_season == "regular", 1, 0), + park_ticket_season_value = ifelse(park_ticket_season == "value", 1, 0) + ) +psw(data = seven_dwarfs_9, + form.ps = "park_extra_magic_morning ~ park_ticket_season_regular + park_ticket_season_value + park_close + park_temperature_high", + weight = "ATT", + wt = TRUE, + out.var = "wait_minutes_posted_avg") +``` +