From 8d7e38e9fe72804ba9c4e3b810aeb0f795bfd047 Mon Sep 17 00:00:00 2001 From: Dr Nathan Green Date: Thu, 29 Aug 2024 15:27:02 +0100 Subject: [PATCH] drafted adding in mim code. need to finish --- R/IPD_stats.R | 15 +++++++ R/mim.R | 71 ++++++++++++++++++++++++++++++++++ R/{ModStanR.R => outstandR.R} | 0 R/strategy_.R | 16 ++++++++ mime.png | Bin 24645 -> 0 bytes 5 files changed, 102 insertions(+) create mode 100644 R/mim.R rename R/{ModStanR.R => outstandR.R} (100%) delete mode 100644 mime.png diff --git a/R/IPD_stats.R b/R/IPD_stats.R index c01c2ec..4701401 100644 --- a/R/IPD_stats.R +++ b/R/IPD_stats.R @@ -120,3 +120,18 @@ IPD_stats.gcomp_stan <- function(strategy, var = var(hat.delta.AC)) } + +#' @rdname IPD_stats +#' @section Multiple imputation marginalisation: +#' +#' @export +#' +IPD_stats.mim <- function(strategy, + ipd, ald) { + + hat.delta.AC <- mim(ipd, ald) + + list(mean = mean(hat.delta.AC), + var = var(hat.delta.AC)) +} + diff --git a/R/mim.R b/R/mim.R new file mode 100644 index 0000000..f4d3fe4 --- /dev/null +++ b/R/mim.R @@ -0,0 +1,71 @@ + +#' Multiple imputation marginalization (MIM) +#' +#' @param ipd.index,ipd.target Index RCT and target covariate datasets +#' @param M Number of syntheses used in analysis stage (high for low Monte Carlo error) +#' +mim <- function(ipd.index, + ipd.target, + M = 1, + n.chains = 2, + warmup = 10, + iters = 1000) { + + ## SYNTHESIS STAGE ## + + # first-stage logistic regression model fitted to index RCT using MCMC (Stan) + outcome.model <- stan_glm( + y ~ (trt * x1) + (trt * x2), + data = ipd.index, + family = binomial, + algorithm = "sampling", + iter = iters, + warmup = warmup, + chains = n.chains, + # thin to use M independent samples in analysis stage + thin = (n.chains * (iters - warmup)) / M + ) + + # create augmented target dataset + target.t1 <- target.t0 <- ipd.target + target.t1$trt <- 1 + target.t0$trt <- 0 + + aug.target <- rbind(target.t0, target.t1) + + # complete syntheses by drawing binary outcomes from their posterior predictive distribution + y.star <- posterior_predict(outcome.model, newdata = aug.target) + + ## ANALYSIS STAGE ## + + # fit second-stage regression to each synthesis using maximum-likelihood estimation + reg2.fits <- lapply(1:M, function(m) + glm(y.star[m, ] ~ trt, data = aug.target, family = binomial)) + + # treatment effect point estimates in each synthesis + hats.delta <- unlist(lapply(reg2.fits, function(fit) + coef(fit)["trt"][[1]])) + + # point estimates for the variance in each synthesis + hats.v <- unlist(lapply(reg2.fits, function(fit) + vcov(fit)["trt", "trt"])) + + # quantities originally defined by Rubin (1987) for multiple imputation + bar.delta <- mean(hats.delta) # average of treatment effect point estimates + bar.v <- mean(hats.v) # "within" variance (average of variance point estimates) + b <- var(hats.delta) # "between" variance (sample variance of point estimates) + + # pooling: average of point estimates is marginal log odds ratio + hat.Delta <- bar.delta + + # pooling: use combining rules to estimate the variance + hat.var.Delta <- (1 + (1 / M)) * b - bar.v + + # Wald-type interval estimates constructed using t-distribution with nu degrees of freedom + nu <- (M - 1) * (1 + bar.v / ((1 + 1 / M) * b)) ^ 2 + + lci.Delta <- hat.Delta + qt(0.025, df = nu) * sqrt(hat.var.Delta) + uci.Delta <- hat.Delta + qt(0.975, df = nu) * sqrt(hat.var.Delta) + + hat.Delta +} diff --git a/R/ModStanR.R b/R/outstandR.R similarity index 100% rename from R/ModStanR.R rename to R/outstandR.R diff --git a/R/strategy_.R b/R/strategy_.R index 83b19cb..ff20f24 100644 --- a/R/strategy_.R +++ b/R/strategy_.R @@ -153,6 +153,22 @@ strategy_gcomp_stan <- function(formula = NULL) { do.call(new_strategy, c(strategy = "gcomp_stan", args)) } +#' @rdname strategy +#' +#' @section Multiple imputation marginalization (MIM): +#' +#' @return `mim` class object +#' @export +# +strategy_mim <- function(formula = NULL) { + check_formula(formula) + + default_args <- formals() + args <- c(formula = formula, as.list(match.call())[-c(1,2)]) + args <- modifyList(default_args, args) + do.call(new_strategy, c(strategy = "mim", args)) +} + #' @name strategy #' @title New strategy objects #' diff --git a/mime.png b/mime.png deleted file mode 100644 index 45413f2fd766e1e98984eaa99fd995ffaaa0c3c8..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 24645 zcmXtg1z1$y^Zwmkq#HzOL{gCMR$5X~x;LcK7B`q4ggpHb)`qoy)w4u;2EYbuCV(P{U=#fea>R7ocW562{V7EbG^@z z{}d5H-5&A9g8Eh^aDKaz2Ne8XP56i$8ib=3tf02r>Vlr)#$bB21iF&udE~_Y z=Y!Y5$+oK}w11kozNYPXML-QYn~X;62~AlWTVumfqfm_m&U$aP}meg;*%P zz~{CIMkIxAo66Gs-&o+m74d^;Ug*v(CTQISL*$P~+z|3Ub1O>$0YZ#QJEjmlGn6Bw z#(^Rf@s*D!PNy-kHJ;1kp|+E*c7 zZG;O#zTZo}pt(T2;PnjtJ32qqc)@4kQ9CTYzW_&pP4tOd_TJFs;aKm4!-y^y_2-I) zpzlkR5l^0Q^FAnSsEK<2*&2h6)c7&P&%!UCJ$vbQKqQx@r)+Xx$;2dF*~Bqg$-((e z&ezq^E|@9uN15p9v90Msm-)yYq8B;Gll8!_Gqyt5@&KDP8x*x_2%54Z13R znz7B`&OO07nxU%WAHj!<(?YG^+Njx&VAFuNfc>}+6GlZwb8gIG$?{DdIzJe3JM%v^YGkh9a23+7c^Gho9Zm>)Pl-so+u=wflGB+ zAAn+S5Wl|f7t{T1B!UzizB<3xtUp?*ar;)VazN-ux1~z1Uo?aE!?W{RStu^H0|P{ zN-|j`R-UjT*e=Iwpik2g5l#wAnY`-;GW5;BWSG-&=CN zmCtowiWgN{9Y@%bJdHH5{I_jO!jh(HfY+z;!8QAcpIysfFZA&}spp;rueHdU{v73N z)bCK!q?|G)GxSWaYNTHDJU9MCmpJVI(3Q1N5a8zBX_rn^!N+zG@Y} z{KnM!b!Ry8r?c#jUqex?$5R|5BmFJn7B7zzl>7b$2=BxX|4JNgv`wyQ!jJG0z&~Be zD(92y9pNijW&E-;mHsWu|I5E&jon@0{-#w8FHG&7#7ZaE+dIi)T6NON{owWydU_ZB z2tC%Iqv*Kro|TfV^KMAaV{+ss0XJEy0$gw%L2Lz<#RnUG;X^`eK$)S`wTj zfp%o+JJ$&5eWWUKF#vCbhLNXvq!cdSbo@`_W&KWISvY?mMQm^C)wD*5@f&|>Y&()D zFzF*#($S63KAeTyvw3hBHz#MxzpU1cc{m{If3EkOa#EMQn`Xc_1!LbDxNImJb1JqF zDkh!Op9>CBzO_dI;T<;7HB`d1U47jUMuJKvq^GRI)hKm&X!vb+jts(jOEQm_7qX1@ z^GiwD-B->Tkzwp;Rt0hF(uk1H3A1MGHEd&~H1uRyBmIa`CQ4VE@3634^UvTCihR$& zQf_BSAwG8F8(W$+`|pV%;K@&(zf%Zu*$cA5s-~VHUAAaZaq;SRhdyg<=bs0JzE=m@ z1nQeGJ~;aN__O;2N>f2i%ALV;V&}6nMtiO0oHc%p!$cNK zcUo9F$acjIiPIE_pw^e?;o|)~fx36Z^iF52B0-roo`#NY)A3um%?Z!Y%qIf1 zh2*w?0Tug%>m)#p<1a~x_8K7B*!UdtS1bEo=f<#0PTi*Z$k_NcD^28Zh3}o}A~~I! zv7cL4uPoVgUGjw0U7i4JR6%aVS?_A++$KI{9Z6f`fFEr9pTGLoaR;8`<62;Ax+$%( zy!FkdQinXuqIx(!zn-BTMlSfS`*uBb$xOrPpS2!bo=70_$rD-;YvN_j?7lcb)zQi- zo<#=#=WMla-@bKlie8W)Rc0@sPWf2BK6 z!e<~bPkeBfyXdzrHLFXlH}|PC^Z)qnuTuMcmEZ6^ffuWZ z3vQMM+oY7?*ynhrdut+fWhORx7#9e9>aW^JA|;?sR6qMx!>cc5t(vj7k$WZiQa0xu zPx8TFaj3b#$Ha8ia>$3j@GiBc@gvPAEv$7SqT<7|m)XMp*2DUCt$Vx!Bt4s_|BxLG zs`wf-VDRhvLjdn~W3RYfPiNFzW^7>UreXEkXjwUi>RfJTxGnUyqqFt?d&v@u#-t`Q zX$b!M%J14jPH-7MV{&p=R`1ko$}4s8hbz6GYlizRDP=ch&Uj;1*VWmg?%cX2!2Hqd zdkiikZ490>Fo<(vCdCZPf;dCuB%b5lF~_u<9b}w}lA02EoY1MUGI7#G##`jjx8;|g zhe#Ql$Q#jl4J?rv3GJGVmgS)YEYloS69T|!t5pP%2V$iAF|~jN-VPc|CI=tDk`e|l zi3cqR<@^o^W8JVvKMeGhjOpC^%3t9-Lo>R=j3;M59iGq34p=|H$sqD$mS|$!rfj&A z$qC3It}roA*eJz=7yM@ZuUM;6VA8MQ3v6=U&t3KW8LQ2`dXA11@ZhQnKfgIX5dZ-~ z`S*t2sQDn1h>6JUqJV#gDe#@C<$vTE(E&`#q~9tl1CxIsQgF<#hVgha;m?VVe8cX} z1&OQsiC6aBjQ2=F;EzhtropT2cr}xiwKXczD-63dHKFk7mk&EpPQ?rhDx=zLt>A!8 zZfF@S8HRwIPi?h0Pmre1h1=H5w!ns<9^bwdezRu=XJ;o^cT{r)efa}7q%F{Q8A}qE z6d@V2%y6BqgoC`i{J^Xx>U#7vh*UuFH|Z_Z9WAN2n;>{%5ZL_Dg=Ft_4oBLD+aJGO zmZ_Pnt}?b{DZ3d-&RviRFfJAUZHp5h8^V2}x_5+_>Sg z>FPjI4n{Q}U1*$t(=roK#n8Y6|A_c`%pj2Zlqv>e;}Bot`GpohCL9nxXr?cndWa)u zGZL7O#|}*83GY0|3e5(~;mzwL*?OOyF?etFF}R<52f4q!#DBzA8tncSpLtCdF%+n% z8QR!Y10!ft9Aa1Jz+a}>Yh+|DIzBicx-0Te4L6of<54buA}K5FeNN8Mu7GvuRGN&NC^67A?5fD)c-v3AeiL098$+|& zsN*{$yKTBi-2f1EvQk7=y8q%RR}2aaMqQXo$1RK!E7%c2(9nPo9lfjIZgyfu;TTA# zDszobT9{YjXvl>9nE!d@R%;_iX9kJPn~^}AbV_9Z=&c53OlSJrn*mI+$sPA&*;I8d zSk;GmU_@QNj`uHCqhrs!BH76gh?{%T2VDspM@%3j4-rw zC$g=Rp@V9d%TJ+km7t7}U8sC-D_g!&=+81@dg=JAE4(xsMhXhD@3{g5 z>kHBlfG{G0dvh|)XSc|iDYW^6not}i@(9ke*ZG`e91QQOQJHFBB-9hl4`Wa;q1xq} z^G#&~dY5)+B7L_N-fIvjoD(k6oHT+M@ZHyqj=jpnFFRwLw>!gsSOs-jE$w%rQ=|d`pZ}6vxX9zhBuph2V%zySLn8)1ZcYF57DsWF zldyxP&Ce!6qS^D7b-9#cuRq)s^JE1v3Wp6oMdhl1G;lJOtMO4CDlq+8kwA(C zJ}&a}h|t@&hhe=k2FKfY{fEikM>?H$e3i_rnxP=m{k3Dyg9p+Pnm6XS`Lq# zFqn&*zsR|0-3QG6ZWd=@lVsMhxLVt++3j*FyYNf?AN>^8rO1L;*U{^ieKq0d(+S~UWN3NvVibzYv_hfqOK>hpRF zHbpkd?F@)O`f;&{BX37O4okxmr1Tq)Kky%}JvwrF;)rQ-@d2S5(<-yd|f_sRspV zpAjhFK7*3;HlhqPrMS3x?2q!Pww3O)2Cg{aDi9$35-5?Dl4PCl=~bD;(^~xU;Qp+X zuI6^YbgPqP`#mi>4FPTlescLIk+Ofo=}cF+(-F!_&<8vBAZj)4!wUg2*{p0G1vo4c z-;jMAj8R!0?+bL5*uF5Z(8SiZn@QJ`?0_5NMgfL(SeC`e7zlZa%Eo9qDdJfodaWgx z3J+yl#ZF`O6MHB%L8~o86MI0{bKoJRAgD0xpXu}xp%kR8o64`mNdy)B0%#vhHz>X| zX!6)}76i@k&9LNZ6bqEb(J@pz3)x86Op?{_o`H<6SL)bX@DWg>)83SLPY8itb^OMtSTCG1MD z0v`+=+O$vt@mp3^K!tPBu-LWt|J?YrJai(2{pG6IGI>1mt*eXyfU``CZ2bQZf4=Sn z-80=SOiiKtvmirk^Tlr!X@G&;esV)e9g)m4ct>@)5z7Vb1- zEh4hlrjAEx`UJRqu~KWL?VL3rp~QyLqygwS03t}9A`TlEMNWUQK=w?f)s#lUZIXd4 zF%7`(X05;!f|b&;(u=rv|KwVp#;L#x!-76MIBh=a1oGFVGxAOzX5`!t8wM(9?j z;`J4)SaEPmR&Yc1420wMMg*857ujS)Z^w&RS7Fh6S6E!jT2P%gpgOCT0z{T?t;pk# zEUT~xOkF95;^u~YCpa>D9S2JHN`XqRsGw!1LjL$Ly=I&AzCB@0I*mQ>9kg9J z$>snJCQJ=QJBCG9eo&He7HY8rgfb`rZ0X)~g>%YNZygE(31&2Q0?Pc!gQcSg2Zg{T zKxyBY5p~*t#XdwioGZN7AwXjFeViKnt7HDs^gZ8K$_VPvX70R&^ucD|0!`zRHW6~% z3?KfJ6O3d{({EKMR$7jDwYyG>GqnSUJJ!EX z+DmeASnk?+JKc>`nW}izIAEyMW&iZU26ipP$Cyrb;ts48J46S1q?MGjk)@9QLnB*$ zN{anV^)s}YN%TCH6PX6eIM)T;?KQH|oZiK_I$=6<=3&Tv)_@4QmT}>XYb?2LU_37< zMj}hZvBj&6W;?zxR4bXhU^QgROw9lsPOlNL9XgGc-g+M5ui5hd1m7P`ZSkl&*SFrOmVWQLG^@~jFN6I#(fdvPS@*_hUhJN%T*mTQ$vq17z@X*d42}o00hp80!fO8~ zK$UI`8uH`O7}Og@P{`zQ#^?ntWCV75=txhtB;GlhSgazhR>}#d#tF6X_zWE## z=6sn{0gvxTY$K2TZDb3{(9LyDR-Hdza+uBHd~sC+EfL3LjyDKxn)=H-GTZ;HoB7#u z)5jpOCe}n9FBZW`DI)e~c-C5MyM3bP1&?MOFb);ZF#wFz%YhL^;XOyo&{9&HVg32E z6R2LcUyyP2h9ZR|aj%<0blo$wZtu=CV~xu8?G2$>tURpq)BS|D#{qaO+i&^yh?BAN z-xme0LSyZ~JYwYzO~B22?>Hw zLm3sJ00qi$oO0QH_>$ldCw>4~L#MS`g31jnGeI(guq`-_9BhinrJT7v*nkD-D3FOZSkAIPeuRbqBd@j4|yiHo|h#ohj= zl*i}%miJ(K6Y=EN4<0TSaq(E!n0Bfs;a(YDrkiP9WnR1yiTLx1Rj{r3y0LtskJy2? zpJudNf>HMyq;dT25C-=RxVt~;b$jlr#x|+_-M@g#=n06Me?hJRzwTWHB{nA|SBv}8 z+$t=N>hFLN3x-p=J>hm3MSP7Q{oG}B7MNcw0BF26K;)DmG7Ja32`%Gx2sJa>7E z9Z1{pUcjF}+700-bl{-5)Q_Q|A`=Xg{s&3|!d|!bmM0PfVkKryl#hig9Eaq^AsR>o42v?m z*UuXnf*DZ;Fj4Nke`RMP82~p%$C$ zO7Oj`z4_}O`zmDn-6$zb3TIJNH%b*{Y8k5Dn_vX;mU25s*I5SL=H25GPL2p!d>L=e zMw?8-rS=qOyA}0!$Eji@2kqdjj`xQ4ett?%>HTp(wvMVF5fnI#sz-98#Ad6Zd|m$J zalCmrEX^KJXIy=VJ90CSua9Bxhud=ERI3_GaOaLOK$X|4T|w2~aT{0%Vqx@w=v*6A zbp$1HOZV0+v;-%B_(**W<%8v)dN3gXxho`7@bXF;HUTg%4H8CatS&5(fzDJ)wPym> z0`C()IP-owRls9b{7p|XU(IsmuDMxg%v~qWxD$fqLBf^~>G%wbs^keEq!7}(+dvFo z-K{-63eeHPCJf95Bj=BI#N(C{o@or%nBa?37$TeX)mDRVGUk3^N_g`S8d(;AE3^`3 z6+p_r#4}2Zbd~5w-D^XH9Jo#p`s4pqh$RL;ta}i$9H@3RI6p5#qaFo>vMkY$Dh2{R z$BdTq;f%6@MS*yMtl2gc8aE8VlUd9=M9)JgCH~dPV)Va8&b%u){7~EuC}+cM!k_&q)IJd^lFoXwY+aeh`WCK)=VKLOQFm* zkVXUh_HTSno24(F?2@rza95#=-Yj=?bPAQ%@UaE@)kG_Y-;|X*-X;~|@s-ezl@T`p zb{UHfRJ*u!&Z$BA*LHybYTj%zie+8vg|4ry4u2h)m}_x5=2;+Y{}6Wqs%PE<9Rh;; zn7$P1;rYP#uRR~#29#%88w;EC-yP0+Dq$W-eayGD(ZxR{a`kzDOC{~W;ClWAnuU{~ zPxnSwOEe{0CfvOz#@aSJ%k(V0{Qy-{FLV52IQPTqU!{@Q4L$&QK*}gv-+CUD*XXx& zcdc};D+zPw5CccsUzHgj=8*D|!95QgW*_8Jx`M}!BGXQ|nXD6>fi(DZE5FhdKFG}= zYzyT2&l`!U7|UWz^s`26@RXo^9uzF8{>gMEc`HGkQJbTgZvIV&dG0!P(@P&OcWUds z9fyS&dS%+Kl9tSidth4QzITej^ET81+B>r)KuzejU#aSTxISp(5DU{06jNr+y5Ki} zvOIPeO6l{j`?FB}r^14gnKW49r2|5nb~^jkrr3VsdXezz6+7-Q<8O<*tS|n5FM#U; z2KBK{nx+}VoajEoc{||}1Nimcn1*D)WUj`s80`R%NTfw7`iH1uwtA>6 zT7y#VJX|-y`US3b~oQcC~@7h1@o;Ec2-c;<)R_WYqNTWlQpv6 zlb>k;3cbQ|h^7UEVbD>P*m6+I!2Os&)x`T{Wo)x+EL#c^#ENg^Vdv9-zvwv<@L%qw z&w8IJudJ@re2pF^RtTm9C8ylL_{2oP9TJB(p!J@rA3ZuI~$vYiQ|U zmmJ8Mh!QsRsHplEKk(i>Ntf{N-}9_}{f7>qa1SBce)5ah98zP2YiVI0d$+&F#XmkP zv&2839f#p{v2Z2Gg1KHX5Fz9heLZZ}w$b)K9?_gLMvMbtslu81{N;OX1>?W8v#dnT zybotA>~vogth-?-pgT3g@smY3EUDT!gCg7ch<65s*QYAqq;RcI=;TEBhl_(wcWt2^yoPBNPV1YzVguLIa-<{v%J;B6TXsew&QK=M zMwe4$K{YiZKRT|%o_i;i69CY1+lDc;zh`RFt0BwUT}&^a#L4{DuKh0B68*3V#Mgm~ zOEPw6iFc`kVotLt_aT|S%f{-dMt|wd2RyWr=W2tJSZEWbB7ot8VU6WHgZU|Dw8XG3 zX|0gsTl1J}KLMv(&(!);(`QR?udJ1zt4s(`AGF6Thv2h8oU10*K9X!$4RoobOuDNF z@zt?hMJ?K*Kz-c_?^k!0 zn*n%lM~ygB^$E9#qH3sP?%t45D!bS1UO(NQzpF6t>pkkbH)jX59m-qS_iUEeZur`_ z%O+wQ&o$Iz)TTeTmBnB&Q-_mnPOsD3um<(mAcSnP`TFSHZD_?2p=EhYx3cNGS;JB1 z!D@=gRM3PBSq3|+L7lUEUq_vK5?+)8YeSg>@=QKCW=>fFOM-R-UiT#->giXDI~G-* z$KO0_8;=X=9i5+K=l#=?)2OotfHGM(nX^Wf%_zB-JT!vf!FKn+71QS<?E5_4|2#aW$0ZT#@ls)HxHzX2^|-#Cz3ziW4jOeQpUBs zwfI4+NrG_CJXfZGC4tyNC35R10Y@8 z>FAF$=*MhTTu%uaK=S@+!J;fgj3X8V>q&ke#*$3aDiF`dX5z63cNKUNk?c~+} zJ22^vZP|&yZc@W2I!iZzZ#kS zQVH&(1^FQY1@*vTxmiV_CRea%NDySpx!ir0rNfWMRd+koD-R@`vgC1s&%8kyL9Rb8 zVngfzrT6!yvb(ISU&pP7GFK{sLX2g_HLWS=5`7!~)8WKk(Atyaz`RerL>Icy?o&=e zZ<1gzp!pgc+EGHT`&Lv*M!p>onT$|ud3vJ+^ZTYkyL30d_kF>N$NS2>gswv7E+^!l zaZvUk&9;7~SYbCq`oHI$zlNJwtu3lJm&bR$h&t|gpXhDf2yrr}cbA*1{qlLCFv=O~VGLG2HfKM0Bt2Yz(M3XuwB?0Hu~W`s#$^)u^kJ_-=AHv| zzFv=VlGo{i`}Ovh}4~(__|beJ_Pz zam3gF7eRXXM4HEbQ>kHB6j)I?O;KN-!5+7abq0b!BY%___1o%xOD+9QY3Af`Iac97 zmMMJucz2d;*zk|L$*58QB4*k3BC@9%qcR~?2@?CREz#1VYkmJpXBu7p5rZ*%kU`2( z1m?O^B%eOwX|3vLM9Z1&qtZOHy0YSW;+?JiSgDUf^xKJ}IK=B!hiZ$xzTcZWdkF=u zT{i_DS7fi2)M*Dg+txQT^H6kTgYoOv$TnA6{|hn)gLEcYU$>2q2JGCdj>Yut*!y7X z#q`+Pl16zB937w@@uO`#45R1n6cCl3_cUPd==rE~R@VZ)#|=oSkg*=|QuT+4%f7Mq z95`x<*M1F_R}LC5i@n@@+>WD|JJd99hZJL<>Kdl5{8({SiK6QqsGpl(+hvyVxGUdw zWuPAQLoSgPNDoE^oLp{}%(jizq$4K7l)>#y&bw>Xps+S?FuxG&O%x3i%!H^F;Q+Jr zeTn><9uTTB%HC6s?_2b{LbMGGUGIvAn#DN4!^z2e~cRvFU#*U zs(RNDhgfXG_WQ!MX7jd6#KEB~CjM}PKYp*_6E?5bWc0wgxjcEzuFh%h@ynDl>f3iN z7-MQW%gm%5onL*{n0SD$e>JYk@EI2~coBfr*tr~r%hbx+DY%$lt$obujbW!W-78Y{ z?km~iM{(rvTh?2K-vzPk9vf0xOWNDFkMJ68hVg0(gss%MyN$Jde!6Rt^8RzhmoPBR z`n#4LBIUWIq=jK>Ltl1W#@CB3U1g=c)?Fc7m{u*CY&e6>20LE|f~RTfHExT`bYEii z+CevZe*LT>=hAeSK!FzqqB(j^4eJ^@UT8liwM#mtGg_zR2TVBn!S!}6-Dzob%zrqL z8qc|6co<#S+FbyjTi&=b3ve!5ceA01qztzeP80>q6HbGcU9W70ZUD}S%uYZXkry^L zE+6!FQ@ucLF)5#;0l35l9=xiy&VN<<@DOZCIW4Wsop+yk@5uoiKcwJmTNm8%c8^&X z#0V3e7J)$h27h`o?mT+bGk@+tOVxHrHE+~ZHQ+=@YzG({f0C`g9glVRf-#57mot_F za`hDem5F9I$^tY<_C_b-j(uiqn;Q$J&lS04!@;5$$m=kf8n*9!=2j ztaa$d^5%|(+X|tdu^9$WXJ)_a z9o}(6_j|WAziJKsLCQsF8rIQQm>>lvt-YoTPi8wsN(@{x+TT8$D7}sjA{~z!uliUQZa|cV$)>Q1jm4PI(Udr=vN^3D{_YY}h1N z7HO4?q00lpqTp{N%s?;QGUZ_OEQFJocnwq7zpBJ*30C2z-Im|#69k_5r@4S z@MbEBCj2Nnuy{BCTLZd`{}n$uEMNi7f_DR1^N&m_JV}C?0~F>@w$(Qf<%Y_5#v!By z1y}|A&ENeV)Dd;hjymwD%M7&}B)x;bW@IpYo_pm$CmCV5KA3%49klpR5Ioz^%QU>B z1M(!?5H%NZ7&8)mlgfr$9or^}tLxmPl{sjK!Z2@D|1 zmjb(_x*H??F6Dq{xhzxK(B)K#=2d59H+SNt&^K!?YNEj*4c1pQzm~{I{f+G@MHjPH z!~uA2_PqZ~E>L(m6Z~J<9oPGDHF^r9E3nQ9ApG;fn|Al-3|xEfj}9sO#3NVpg1Ifc zn$obQJkEe;ffS(r@TU2S^pyQK*{RGP1GiVU`1)LZ^t5WTXRxb4fM(uY(;4G?29=Vg z>0NRl?AmQ| zvxj|m-;76zM@>t6pDFoZDygCZ4PUZXk}_VADw5tGDdN`-$U>!4nDLq ziCAk_a;qn{AUguIA=kM1-MKL9*ar~1`fU_5#_#qLAA<~q1dRmrud@`&_i*j^$gE~= z$DcQ&e+9wC*+tLqxTf}0s$K}!)o6UqK^Dg6ai=BLq_10bOT8;>+XXOF(|N&%_A@d# z4m^919WAWZ>@9!pvIMqAU=TTZ1(}~btTmt&z!&E}ElANeF)>Ly+?iqU{jQz`i%dy= z2LB$Oq?4_}HAb}m{T>~MCfjtkwz`da+gh$_oSunBZt0VxZE}8|yq&H@2hj}EOzB@= z3JqlUURX|Lm0P}vjzbOe0Fp^QmV7&s=wox~G*PgLNa6F?+!xtp-p|Q{G!Bto%D3Ke zW+H8oyx9ed5F%q&=3 zNv;7{Pbgm&gcu6I#b)ZA3aT^=cmkD*zlG_)Fxt=%N>5DD(ZF8+Zb$s=iRu|31lq59 ziaEp)$PE>U+W2oCz>#O@@D%J|-O#hBqm5Lvq}2!3jt$HO9L-~WEn|OJfsK-NpGdsQ zG^D_s{NYX2BL-I3&0~@M>c}V!^2N5GCj6%4!-ac{Dr!5GDWz}{PVT4CO3efsJcb{O zfepNr-#rq063oDQg0ELrZ+UU?%jrKUZ+ciMB?@=3*@$zuQSzfQt!Xf)tm)O`h|cyu z$OiK%f_(kXKFpndN~Ea5PcrhjBQWuuwvnuue7$_9NJdekdibARJSw7Rpb#r>C=e>j zlady+R|6JCz$-ubty%li6}fF$DHf^L!ZiC16N8MrelBvm1KGb$oJ_7rIq#=xW83#? z_KBye3mMfs!DjyLRY?fowc>)Aicp;Hb7Lr!x#Ay`a->xBefRR9>%Qc1Y9+&PbA%oT zWA`!w5vTOsmCS_ZwW#40b8N5oPM@k)2)d8+mscc1^3RN$Ue{w}QMFs=qJX=cv2tuy z>y`?8=$}86KW^Ba^(^m$I5hmSugDd|?ZcHeMbyg4wAOB4BE7O9sV7U-eplE`7s~G+ zP^(~@!pSf(vc>>Z(2U$|5@)5c8|e-5!GdJY?aU{J9&Mq&w28L)Qlv%@2cFLxo&Uq~A=ne4X2-`?u+R zRr^&vjO#E}6I)-Wa@q?+2$C#9z#^P9c+aZ*I;5Eim6+zzHu&D#oj2q@UX{nbYSs*& zutWO8rI=>mcd8`*stt)I^Mv2RJ}_ARLzwf;Bz=vFw| zk?#Ab!o?l0ur2V~gQf^Xg(vlTe-Bn-@!I3^E5rLfty#gw=9dz9)f?EfQ46cV8JABX zZYf|t?~Sd#qG31wu^}9Bfy{*+0l%($DRo=w~0@P;N(p>v-*MyW)fw-0_GV`;s_n0AzmIn zK0ybQH-ZWMzv)6#dyBBiYaBc0eHbmpNs`OH8Y&81D9gzuz1KwkhnpqV!(N6KxZ>pK zv|ne!Abd{N8g12m8?ZCVRK3TDVgw$b@0~NJ>uV|;T2X8S2em4u3EL%dba&) zG|#|gjn`ORB4aKws8ZsOa!cz)w(Q|MK6yk{bi&W<&Wh>6AAz=*B9nv-6_mr6UVcGv zarq?_WV@3e|CB4%=XJphFk6BIdw;e~O|kHmvms#wsz3$n4T0IiOIMZ4`x6|s>6v@uA!yzCoNE_J2cs{`8d_ClKvGkvy@mu!@kch!Bu*mJs<|6ohS%{xD1 z)~E27Uo1fA8N(SOz)LT|(GAVesH+_V_`zkh&eY!rJVB4=;pb~BwgzeBcn)D+>F0>c zj`qEk`4h<5Ny<@MgiYRk_&*fhr)eK3&WwZakNzrRlshr4P#yT#DOsoyo)u znq02W4m;1x$bXpJx6;5LbMYMBND-6O;0F`zdUS z9piCIG96GG0|9;gw*KfGIK!X2j@U=?V8rXBZ?jvD%Wl0UlzCMw9A5|C0ScIm_CuP&R|`eCykEIs>p4bDmUW*JsWqRwA$>i1EOKHZVE{e7Y5dhdcTA4Oq#EPe~aHvWck%r@AZn ziWuv43i<5BOPG|MJtRven_O8Ae0Hn<Nn%UCLmj#cNY=(tf z3(%azXi{{BT(mA%fBv#|#ju&m-#?8Bxt%I|;69>|?|8fOs*$wJ>Wt2aDmGAea}Pe@sRUceLePv9T(UehH_Z_cKXvdQ51Er3QdKJS>fn(}Dp~dEfmz4HrY0#FP~KY7_PQ zq`K%O^4Gm!5i8bx2423*OlJ77-++t(##=M>HO#>07Rw_cqQ4|jm1B#a$!HF)FB>!5 z28+PXrTVebe1cv#qN|$qHT3MCAIAoM$=o`|bDx;;_;Bndx1M(efmh@l4g$FirL(DrC#8Pd3GPM`h`KRdWH}UbGcd_K=E?6GXf|ZN zZfXN>mwzE@xmOrbHb&Hplk}l^p<jDPA)bFs z9^ZE1gMZqzt*!iKm<}Llf}qkBYPVy0wNkrU$>g&tZ|_zO%(Kdei_KrccV{~=WvvMt z{5t9?5;~AO5UFp^pc>@>+uvkZVjy^hHS6j`lS4-)A;m)NCyBx{u4H7B=j^5ZkYDDZ z4ntx|I+m&-H~`F`6u{={4NHoh>x~z6((`RjJlNg$?u~1|QilIbIyC#L2Rh8Qo=H13 zob@id1;E!In(lAtO|}Sf+wWb%8gTYp|HP@i+PIC(Pn1K6nIUkw9TJUYH&j|M5f zk*A`2LISKtFk%Lopyq2=epMYA4kN`dRRTOZR;QV@`d4R^&P~@E{Rp_7kmC8RJu?T4 zq#|qE2r)804|joLfhz|ggge=2BYorhb#JD-8VIPkA3V_R5g5@fb%Up)7ilhPYeWuv zu@~Iv#&e}sU^%DOljWm4Xnu@xrAUbO@ky3m_m>0a1}AU3%{w0!VKHqEzWk1nEtX8k&eAN=KxVNG}4?3H2`D z_r1UO{{Fl7d6Fl4a`v3v`ONI@?Ci{7*Vz&^z8=b*hoJhZKwwieS0**oBa$5w86y1j zmYz8C(wmWixheDoRel_E=1UG?V+ZC8ME||Rk&W^C2&JtH&Y+W#nV(~$g@}Tyf8m%| zEz)a=d;mh42#~BnEa{znYMy%3$iUsyvUWE(O`Nn0oxF9HPiOX~?6~H|8l7^}VM9A6 zy_*d&Ihj#feP=FZKRclmrubA5MEBGPyV}6-7P$`s-#gdQv8rWZdZYV+_l8rs3Z%P3 zmeS<8j?ql3G|{-iMeUuOhi3_33~?kpw9C$3ZL|syD{@y)&;om%^hyF(Hdp&ng33R` z$;v{gL}j7YL|7Pp(etB^gN)wj{jK(@cI239j{jaDD{z4?H>#5wHB)r8i4TnCJ9&f; zhSdZ+u3^{-`B>^p86b}@222hV`%ii!-lQaBsqVk+I(Tx}eUQExS6{!g^pPHWNE3H+ zBf=(H;n^vn6NfAB-yF8En1VzA?6t69f{^V%Jas*-u7#2Q1x^IkTYhh6U5v0{wK?o(35-aAHgVeGB;o9d+h97Y@4;(}M9JMm56HcI8I$fm znV#i-tp7yj>bE5U^G00XC3(MUapo&YO+ZdYU1f1l+@l9i=Pd!MM(}0QqP}ZcYT*DR$F$L zJZICl83hIV0weJ^lg?Q>24E($Vo^sR8-XPXWJ5S0{yFarwEX-y#%{2l+f0W)dd28L z3sj0z@-h{B2XEn^O{t~+^gc;po0~Sjm_&Y#VOf5i8{@T+3+zi$QniBHcuiFbmjn1fWYF1 z*Pmy9cZh zIL3tsfk+fmh0S&gsuxx5ojk|8iZux;Lm&`N@g%227cSPv_UN%0g5w{w2^SB0?tx=Z z0epc2NkTg@mI;3MtKfJ=L7^rd99dB%2mn)qf&v$C%crEGx{Q}rNGvXtm*B>>eac8ODkTibyvG?3DJFi)UwuBDJ!#xIo-SRDrs)c&{7$ z#C5tuX8QT8B~nCWdo3}!<%6fFux`24_L37rCMEFtBVGTM{0RQsZ3h#Wb+ed)PLV)J zSpd-kHacBO*~&3ZE5*MJJ%74FNKw(EHN?S9!Ko`%rYzU=&43jEZ(}%7qT?BKVW3Q4{8vRIY6{qiaOABB_8@BI2ziz&M=!3Pp(sri*k$BJkH&7}m0~BT9bH9eLK4 zV*BdF!|>xP-iP-Bxb6$&9CaNoEQRy(yB00;za-}27hu1klo_sYJ00|lR5vYE+058E zGZBxfhm{YP3&>ssU_w;XKUr`P2>pk^B_rxG?1KIB7VryCd=IVQFt6~m5D)WX|NO2q z;}IjxqqMa3MHdj=6Cbl>G8e<}C=VUwG1$VWJILS8qFuRYE-g^*!EGIcm?-J?1BlP< zM%&_f#1-b>(^ABD-_XQ8c%<;`k%C2FB&{I_;qqc)Qm&f#jkpW+Y=F`ujU@->3uH!m z4o8t0VGQC?Ct;!jC3ktSzQ(P$of&kW`My^oj6au_RsLyxVKmU(X$*g%twF1cOY;(T zPr;Lx2tkOTwn;O)?~RtVSFchES$&Q|S5_{2^fSYX0{10HMzq!1T}BN)WkqG$&pWi& zDAK7?>f!m53sjx)m7Fy!==f2aXhQJ;$?!e+3mhek6hc?IVws5NErEsutT0%=U9_U++2yQHa zD&7kSi=5hSlVg|TxgMz=D=bW1RpQnury?4iab#TIPadKkZyM{-=`GQFBUUQwaG4Tg z{g%{+DGg5w=J;VO!}2yBgl%HR;9@*Yx%xjTCep@iVGCGEBZ*#omSg3zhu`xei`FCC zUekM?@#H;(8_Kb)890u;5fX5c#uj(t=h+qUTzqRLp+maE5o6NvLtK*x5e$M=izyo^ z8wqXFG)TuaQm-G3_JgXeiplb>22@g-NlH1>q zdcsAm+L~+Zt1BtRMu^W&$wf&h7e9IN>s(spvC73C%LZnQR2GUW2XcnX)J`Vo<-yx8 zL?GubEbUwD8SPYkEuvMD6qhg@2Un`DxAVqSqftX!kshHq%~kkxq=(_SSsmvi{w+4{ z0-;#E6qs&aQqZ&VqXWMwL7!i##}n_BuiDsSJTQ))!i6q8(9vfXS7)N-VeBW_I0AUFA7-LV0-84sz#IEZUuTv? zV=aurNkwS&!i5>x^`F4l{T4x(bxx=Lppt+FRE#>*DgQ7NiXQ#UK<0B*5K5_|g=#bj zFacS9$6so&xm6!$LYI3l{F67XblqarBlN7DiApO1a*jV&PGancqH=3FSWVUHt})ea z7kO&?oyW*cMxEISMRiNL0cM>#5^U##A-)JbCj0ifnN45d2EUaqw!fg>Y5F4r8PlEF zqOFS3Qp>aEs-($*gB_!m!$n5&^C`YUzUGP<^O>{iL|xJ;(~ zl=0YoM=RngHj)Cu#HO|H`RYHIJ=cSJyYyu-CMkYY@_>Twx(*hfg)2GE0TheDm74^- zDFQxWiC@^b>NXOTa=(3>SDgS&^Pqde-SVTiteh@Pc5~yDX`zE5%`HK@irr2nj!twS!-x182LuRs$k9n1G-T{qMr>C+i&sFlL+*jcBE&g32r=Y=4XQWomoe*)*NPB7I^K{wyi^48R@xz^WpY>y8tj>OIe6%mOVhHR> zsc|oC-3}a{io9?3{K)p1z1rlaH)QdZt>dqSs`;k!z=xHYhNpD6n@}6Pu;^W5K}sd; zZ()cyiRabau=UA&da+qw*<~Z-Yt?gE8uCaZqeFF;7`N`-T~g#5Bi5{j;OyicS)4|X zG1dk8pTgq_hoVlip||FREQc#JS7&S0_>ff4dB|=S)%I09Xv>a5D)P6%8RH%M&c!6> zxeTU@C#W;YG=Xewd{WjW`xAm?LDPNjOfs(|>H1B99sWx>EMDTiDcb27q^$ zsDG-`es5Ejrv%KWJRmA>$YNlzoTQ5V2BWBcGm~p>PiW#ZVrbS4;rgM(TqUy((`-<5 z`>*SKE^f_ibwW24Ay*JqQGnIv6m6rfoxd0YHcj-W^Spdz4u+!)GHX8aELf$1Hu%xm z!W6IlE>n_K_h?`=0T1LO&xGO7{4(ji(o5!jO5rqRXVq-^(RSjcd1^zZ;MsHO1}VvI zb){Lai>!w92FIo{#tli&!hmj}7k773ts?}}Z?0au>v`+AJof~e*;<2vKDx3Eb}8hA zov?f_$lNM4gJ;>nsF_7!xpGseyEyubgzuV2+mwcN zG5;1%Me#y`A8XA)p`r^j`6tQD#93ygX$3EEF-=3Vb={vs9+?d)$L(~tPfn0g>LF93PD+J!6?afIO*A%55^>5JQUo0w-xwQ;A- z-^E{cbIG3UiK{_@+lcj(G#|J`6XSte#7KWr@Y@jh)}ELwW=$#gdcvX_>FJ0^mHWxA zhtU`jE>6~+O+>))MYJjV!<>*qUmj-b>G!%LHQG#YBR8V{jgZ6;TIlxJzBtJRmmqWW zO|v;#MiG0*Vs31cj??5V&+v3@^)263vdW_|4}B(_&5V7xmE?nt;b~I7SoCk?5#2(81XynEAtLRHFOreuWDYw z?lRgwBGrVxJsE>nUi;&B$o36mmQPy8>6b|V+Wi18m8tBoK^OhI5yfh5;2ZHFqh3pA z%ylUs&VB01;pcNyHnKm0ada~5A=iCz%O!L4-0T^H>yc}N)%^}OvX7iEcp=N2!j5!j zVleqw{AsGZ!4`j)N2W1%sAGRvDR(GnM@vT7^?i$=LmT>gyuqf;*U&Ft?GM)cc;J0i z7P($KyM2N&{(CFJw&Uab%^VgXNc>5vFo5kqQ)AYN?0)&~)9tPs6EVVp-nS}9h02=2 zGC-cbqg&XH4Os?h6W3a~_p&_cw*s-rA2MENip+bb8d{Yb-yw#dA~WU)E4c5OTxq_=lO%DmA-%QSko`^wkJ=#SGo*#jqezHH;nR6 z3qEq3da>_tfj*I8jjJuI8P+?)f6Z8HF{LBM-H%%do31$L`ntlJ(d#boH+19La2|Rv zfqAfQ+$)=#WcID)Ho+#-PH^Znvoo`FC=n+)XNJYnNH*^WpDY67fJk1o)-1X)p`W8L zqnK4O&w};44!T0i8>u#oD?Aq2xKcJKA~#!N_S!6G6s*TklE19+U|Zj_aYt<$hz-)J zV=YAVgLquq8ON8T&AFY`f#aQp8(!}98~6)K893jVeLHT|XU6w71g;~6a?b_4+x{}@ zzGCdvkKE6zHaq)xP@Qj1FMLH{sZzhH)_wlaZl=L!QbSHyu3dyd4#w{B3L5sKrB^d( zuwojMI=~r|-lkUCzZ;Kc8Oqg=rrRkc7`y(#u^6AEO zvBx?$ZC2@f&nzVjyi=Gqjr~__f!KtCX^HSR>Femmr}brt6UvX%S_^Ppg_!?n{mfET z6kGv6$Uepq_6#MF>-UYh?zgrW)~D)Ld)eR_AX94V%CqV+8$Rr<@B7yv;Z~?#BSoR7 z>o>EeaLRz=mp2%hd2Dar_jNU9Cu{rfjf!1<>kX(TDp>n(ZLv(fKFZ$2`9{j=D8`OR zd?4>6#txq7+&7;>Aix?K52Ja~2cD)?`9rXs7{ISLJWuD<}x8SfgO=n&5H}2}dBV{ILO-BX^+#Q8naeGo{#f z4itX*w!6l@=>sW(b@f3F_$)Hv#sS9@&hzxz+B!SaBWizVzPdOy%srm}sv>UJ%dkLi zd2S_?8&LDc*W>A7<8Q6?!4l#8SZT93w5B$EJG?2NQol7+qBiFda`&vy3(=zv8dt1g z%JDc?d?j@(p49rY=T~2*Gz;`1PW+{EWtJ`KwTUQpY0B60Sbo}jAltvzxncvkR-nIo z#nAoB=qt$!^R>~#ube`x+2`>^aJkgVnWo<EV}rE8 z0gu#RAMH+{B=^Cx4 zC$*cT$za4>^cSpmt>vsQ?%PPbhQ`NkPLLV!Lr3>P2DC1X?%Q+47W(JDL}~**;prEF ziK1n?6V7D8?{q380{2R?{Ezis`ms8qSSH}tRU|~<=> z4Bo>HV|=oQGXj6;>#G44jNVu$gCxNrLOjR;$i@$4>-^~!sr$UoH&H(FRkH!Bqo^A} zE!#Tvj^{aeyJ*Xy&O`kgTMp$^fkMx2a?X<;5K-&rrbfOw5(ij)%XHB6*_FE6PYs6o z>E<_HGm5&=K#AuN*CsHzr2b@Lm(6f#rr55&rjqH_OO!9M6h0(-3Z^S`_0EsQ@@H6X zCb9D66P3*IqYeD;-yfHhR9f5UGeAS9{_gG`Jp(u1R#$WPk&xJFj3*utcoqb7tA2DI zxieoS9UqYNQuh|dc@L;-2s{Xd5FfxPuAOcvz$fogdikubuvpsxh+`L$KaqoU7Ip+r z{e9ECl%hmZWqRWM9zVh^KZF^}3m&9}6uL%4q!)~cAtWyp^>DIaPOCyvu0^>j%c%L4 ziNf$Tx-lrq8{CRIQs8^AE{cb_th$=biC9g_n1Y_FYYpFnbqKY3C~yUc;IPZ)LBMA7 zNE7-i1m3>aPerBnuQC-C72ELB4kR~AcKg(kVD>Z7TNBFJ$+~Y0X4(zBr*B8J&-!L= z+~Jbf($SuD{*)|HYqPPR7TY2M^O|2{`X}z3JP7s|Dy4q!etkpLeC=~`Bwo@HYX3Lx4SfOmJg0V(xqnk*G^a4*=i;~ z&r`^Dgo-Uc+xYaFMolu(EI+nk=7nVxo=gAYb|hQ%%!A?Q2NL1CLMS8{?V3i4`dHQc zlh%v=%<<|?hkQ|PY0%>1=eX>Yr0VMG1C}?uu8=F0Z=x%JyYjlMd=GO|Q&YWxD>rlO zpZ&w8=H{1azt8Z)7#|++wYnqih%Y;&4grJfHV2>>Fam+l;Q=6Aasw0qg-7w1R9r&y z?+YWMDKb!${m&namlj)oAQB1+{ty2H!zH8uhJRBJY{D`H0Fy#c?GoO9q5A(?#=C4V z+(``M@rWS!<&Jw!;9l-EzmnQf{(W4B548(FKgcE{Y(TikG@1`+%;|k z|5O#ov%PMWKX2`0ZXtn=)l544neAA#t>1Z T*8c)2|Di{UstT3&O+)_+F(gNC