From f75313ef7c6be18613304d535300d48c30885740 Mon Sep 17 00:00:00 2001 From: Hannah Frick Date: Tue, 9 Dec 2025 13:47:06 +0000 Subject: [PATCH 1/9] add Royston metric --- DESCRIPTION | 1 + NAMESPACE | 3 ++ R/surv-royston.R | 108 ++++++++++++++++++++++++++++++++++++++++ man/royston_survival.Rd | 55 ++++++++++++++++++++ 4 files changed, 167 insertions(+) create mode 100644 R/surv-royston.R create mode 100644 man/royston_survival.Rd diff --git a/DESCRIPTION b/DESCRIPTION index 419fb8e6..84663fb6 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -122,6 +122,7 @@ Collate: 'surv-concordance_survival.R' 'surv-roc_auc_survival.R' 'surv-roc_curve_survival.R' + 'surv-royston.R' 'template.R' 'validation.R' 'yardstick-package.R' diff --git a/NAMESPACE b/NAMESPACE index e998640c..f24d3531 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -108,6 +108,7 @@ S3method(roc_aunp,data.frame) S3method(roc_aunu,data.frame) S3method(roc_curve,data.frame) S3method(roc_curve_survival,data.frame) +S3method(royston_survival,data.frame) S3method(rpd,data.frame) S3method(rpiq,data.frame) S3method(rsq,data.frame) @@ -244,6 +245,8 @@ export(roc_aunu) export(roc_aunu_vec) export(roc_curve) export(roc_curve_survival) +export(royston_survival) +export(royston_survival_vec) export(rpd) export(rpd_vec) export(rpiq) diff --git a/R/surv-royston.R b/R/surv-royston.R new file mode 100644 index 00000000..196464a6 --- /dev/null +++ b/R/surv-royston.R @@ -0,0 +1,108 @@ +#' Royston-Sauerbei D statistic +#' +#' Compute the Royston-Sauerbei D statistic +#' +#' @family linear pred survival metrics +#' @templateVar fn royston_survival +#' @template return +#' @details +#' +#' Royston and Sauerbrei proposed $R^2_D$ as a measure of explained variation +#' on the log relative hazard scale based on the authors’ D statistic. +#' D measures prognostic separation of survival curves, and is closely related +#' to the standard deviation of the prognostic index. +#' +#' Larger values of the score are associated with better model performance. +#' +#' @inheritParams brier_survival +#' +#' @param estimate The column identifier for the predicted linear predictor, +#' this should be a numeric variable. This should be an unquoted column name +#' although this argument is passed by expression and supports +#' [quasiquotation][rlang::quasiquotation] (you can unquote column names). For +#' `_vec()` functions, a numeric vector. +#' +#' @param ... Currently not used. +#' +#' @author Hannah Frick +#' +#' @references +#' +#' Royston, P., Sauerbrei, W., "A new measure of prognostic separation in +#' survival data", Statistics in Medicine, 23, 723-748, 2004. +#' +#' @examples +#' royston_survival( +#' data = lung_surv, +#' truth = surv_obj, +#' estimate = .pred_linear_pred +#' ) +#' @export +royston_survival <- function(data, ...) { + UseMethod("royston_survival") +} +royston_survival <- new_linear_pred_survival_metric( + royston_survival, + direction = "maximize" +) + +#' @export +royston_survival.data.frame <- function( + data, + truth, + estimate, + na_rm = TRUE, + case_weights = NULL, + ... +) { + linear_pred_survival_metric_summarizer( + name = "royston_survival", + fn = royston_survival_vec, + data = data, + truth = !!enquo(truth), + estimate = !!enquo(estimate), + na_rm = na_rm, + case_weights = !!enquo(case_weights) + ) +} + +#' @export +royston_survival_vec <- function( + truth, + estimate, + na_rm = TRUE, + case_weights = NULL, + ... +) { + check_linear_pred_survival_metric(truth, estimate, case_weights) + + if (na_rm) { + result <- yardstick_remove_missing(truth, estimate, case_weights) + + truth <- result$truth + estimate <- result$estimate + case_weights <- result$case_weights + } else if (yardstick_any_missing(truth, estimate, case_weights)) { + return(NA_real_) + } + + royston_survival_impl(truth, estimate, case_weights) +} + +royston_survival_impl <- function(truth, estimate, case_weights) { + bns <- normal_score_blom(estimate) + + fit <- survival::coxph(truth ~ bns) + est <- unname(coef(fit)) + est^2 / (est^2 + pi^2 / 6) +} + +normal_score_blom <- function(x) { + tibble::tibble(x = x) %>% + dplyr::mutate( + x_first = rank(.data$x, ties.method = "first"), + z = qnorm((.data$x_first - 3 / 8) / (dplyr::n() + 0.25)) + ) %>% + dplyr::mutate(s = mean(.data$z), .by = "x") %>% + dplyr::pull("s") +} diff --git a/man/royston_survival.Rd b/man/royston_survival.Rd new file mode 100644 index 00000000..1b703abf --- /dev/null +++ b/man/royston_survival.Rd @@ -0,0 +1,55 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/surv-royston.R +\name{royston_survival} +\alias{royston_survival} +\title{Royston-Sauerbei D statistic} +\usage{ +royston_survival(data, ...) +} +\arguments{ +\item{data}{A \code{data.frame} containing the columns specified by \code{truth} and +\code{...}.} + +\item{...}{Currently not used.} + +\item{estimate}{The column identifier for the predicted linear predictor, +this should be a numeric variable. This should be an unquoted column name +although this argument is passed by expression and supports +\link[rlang:topic-inject]{quasiquotation} (you can unquote column names). For +\verb{_vec()} functions, a numeric vector.} +} +\value{ +A \code{tibble} with columns \code{.metric}, \code{.estimator}, +and \code{.estimate} and 1 row of values. + +For grouped data frames, the number of rows returned will be the same as +the number of groups. + +For \code{royston_survival_vec()}, a single \code{numeric} value (or \code{NA}). +} +\description{ +Compute the Royston-Sauerbei D statistic +} +\details{ +Royston and Sauerbrei proposed $R^2_D$ as a measure of explained variation +on the log relative hazard scale based on the authors’ D statistic. +D measures prognostic separation of survival curves, and is closely related +to the standard deviation of the prognostic index. + +Larger values of the score are associated with better model performance. +} +\examples{ +royston_survival( + data = lung_surv, + truth = surv_obj, + estimate = .pred_linear_pred +) +} +\references{ +Royston, P., Sauerbrei, W., "A new measure of prognostic separation in +survival data", Statistics in Medicine, 23, 723-748, 2004. +} +\author{ +Hannah Frick +} +\concept{linear pred survival metrics} From 43d8d3ea3f8a1c8785d3896724da75912dfe710d Mon Sep 17 00:00:00 2001 From: Hannah Frick Date: Mon, 15 Dec 2025 11:15:53 +0000 Subject: [PATCH 2/9] add linear pred and test --- data-raw/lung_surv.R | 10 +++++++--- data/lung_surv.rda | Bin 11443 -> 12800 bytes tests/testthat/test-surv-royston.R | 25 +++++++++++++++++++++++++ 3 files changed, 32 insertions(+), 3 deletions(-) create mode 100644 tests/testthat/test-surv-royston.R diff --git a/data-raw/lung_surv.R b/data-raw/lung_surv.R index 98c957ca..1f23ffcf 100644 --- a/data-raw/lung_surv.R +++ b/data-raw/lung_surv.R @@ -11,7 +11,7 @@ options(pillar.advice = FALSE, pillar.min_title_chars = Inf) lung_data <- survival::lung |> - select(time, status, age, sex, ph.ecog) + dplyr::select(time, status, age, sex, ph.ecog) model_fit <- survival_reg() |> @@ -30,7 +30,11 @@ lung_surv <- predict(model_fit, lung_data, type = "time"), # We'll need the surv object lung_data |> transmute(surv_obj = Surv(time, status)) - ) |> - .censoring_weights_graf(model_fit, .) + ) +lung_surv <- .censoring_weights_graf(model_fit, lung_surv) +lung_surv <- lung_surv |> + bind_cols( + predict(model_fit, lung_data, type = "linear_pred") + ) usethis::use_data(lung_surv, overwrite = TRUE) diff --git a/data/lung_surv.rda b/data/lung_surv.rda index 2aa569a1b772c33f1d5a66cd8857e47e4e8d7c10..b8c40a888e5cbe50bd2e41ad132d3a9b598bb178 100644 GIT binary patch literal 12800 zcmV+bGXKp&T4*^jL0KkKS%2$7mH>aT|NsC0|NsC0|NsC0|NsC0|NsC0|NsC0|NsC0 z|NsC0|Nr1EUw4#Ewewa)H^uG%+^;-%&B@<)B>cmWh{ zLG;mv*4yu~&dGY*%`;-y0(Wxpo!N?|?t?PE$H*Ip;RIT^k>*c-O^mmJ`>h|3m zD{Zd8HOXhQ)7!b#;4qS;gCR7;(W-bd0~Jr&Xr7HJ?E?q|^a^^X`bVj!rbNPE(<*O5 zU=z@dqMj2P0*^8$sp!y|Hb#@kX`q`)v>F(MW{H^q(@j+UsfH%0?Nky~5tK1g@?_M) zYJQ_-Jfqb-r?pSW$rIBCfiX|?jHdLB6!w~o^qQvd6xwQO>UwGE7?WzA$*O*-quNaL zB6(BOGNztW^+)MXZB*DbFvy;%?J2z}`U#4Dq<)HedTCG7RP{Yi8mFo14JZK&hDm}n zX{M)wDeVCo85=}o)b%!_PbBn4)byUCPf?NSlzIeYG->D{Xf$Zak*BHX5rGW{sh~9V z8a+%+8&gd*$)-$AG#Y3a)X+#t5F%umG}H{544988wHY#GGHphK6UurPNoHv_k*IX#k%4rHwRk);1DjS*`5(CSMBpVMuD(I>BrsVD1v180Ko@rd`>uX z4>%GKKmff*4=>Eeq98Ew6Z zIP80pa)?2E=?R9ZwWz9LvYey^b7?OJ9BzsTsYdK7*zGc&?Mc}nItwtpwBCT#? zOv|KFR;%rt&FLp%;54#}tkLX8_s(!oSBTrPLA=4w>}cL{kb?O%rsm~&OxIprXRfHC zxitob)wn@B0OvSRqBVTuGFMEbw2NzjT*;6GJ z&9hEO(=nkX#bAh2?k+FNfwKgj2|v6u4Hjt`D~Jx6&d&eL6u z|7)(n#OSWIulF#){4SKaN*&VfKI+anbj`v;Uh;Q0Wod18EIOw#?$=dp`-WZ%@6MS{ zThiSxD%uOPxjz-gHe%11D=nG_eb-jav&{3=_IY{pF6{0-^3u$5OuXa6;cEZsf;9NyO8 zOcsB7)wyz8TrF8#alWUF6QDS%xJIaqzW?OJNI^zOqAlvKZk|YCD^03$Ttr84^qwWi zygH!xJA(1ECxgDJaHg&=x0K7_flj4;l3<7g77#=hylSqzr^{ z1VSW23_%dYA}Z310t>mk?o-L(6Tjo(4r#oou|16yQ$e7KyHWnL(h&v7zyubx09%QX ziva-<63GETTqNnWj|#(I3ljjNIgLhV+6+6>GSfh2gR7q|#>~X;9)@TD8UPfq2T@k7 zJ`GWhVnCB{3fJE@CRQ~S!;rM4RIam3Y?li+-zT3sYkTVxpVMe8poJSt<%IKuv-w3& zNLdE@7$*Sc4^v`JX6w}v78=!WlG!GL)Lr<7%?=n|OQBv6P+KxzM`Nw+zKgeCjIxpJ z;FxM!r$GN%N_aAiO&Dj^VGg z1)TG{cJsf>8>3%|D_d@+)-TtjLXD!`8j)5F2Vd>NEWWip)piKTw8x%Lht0f)8A+7j zMes~S0DvOK+m`H770}lKp3b~sam(DnxzCV}Fjsb^!ht?a6(M5g>PxYE9Y(L$4v1#O z1=ctNPi}yGDez882cYR+_-4!TLS*T28BTey?9Ka%P!_VN53maaQ4%Rra?FEAkJ(d5 z0R^eQrLmK-mahhuS~Iu_JtYAkqDBPHt_i)P$% z8|5H+=8||c?jM@Z6_v zv}L~jo@G6+pqQ434$!M48oy81eAa68(2LQOWS8FtY{(&uH^q{-0#&Ed8}d2F9&ER zE=;+1Hh~fRoZ&1Mppao0Iye{R4gi5=I#bw;?}!%&9K+`W?+D&?P$&rLLLQZSh8LVR z>@9W6>6ISqrUfaVBtI|QoAgIq(ifB+tfLM8@zY~mPJJFeD7JqU(^sT~6`Y$}HXEWC zo%&0a_lY$~5uuI#hTlvC`E1wU^%)C4c_`=@pX45{p^y|gl~&pvN7hZ?Vsky+=qtHK z(;J>>C_4E*`@Qp{{4M*!LPZk%JPP7R|Dv^Oa8Gsc_p7mK!jr1NDxg?s`zFd_IJ^__ zZY3{jYM?M;V-4jQ9>Np!B)e-ijTtrK@m7)^U1a^^>Kr2gC|Rh}P)IBM7&Dy8Mk@cv zdA%Uw-b7Qfv`!x1UnB)cdAwZVAsTxlU3G@S-M3Mu@DkmuXe2X;^ALS&kKX8>qwRUpKve_EKwk;X zbSZPIr^rl`98;A88NF7jm0SIuHj1z)j=}5%S1)d~BSVZ&@<##;)XN~Y7GpVn6`)rG zm!n05wuWpg5z2Wi?DYyCSAs`+qz(zS)?7OmpmtyB^^TydG&p~fkS*9Sq7AwH&Y32) z@oi?|37^Z#$CP37F+57pvnR>%k-I%`xQtjL{p9_xea*3*3&5w_Nt-J*uyfNaO+&aV zbZMQfHT>{=;DgeG#%$VDq%pPj^T`X&82So|hpN+p> zF3kO&($tz)sUNUwdben0X|30>=3d?K_I{P=MLr?NH7Y4#uMX(N*;Ste7g%_xg7A2} zPiwUV7aahBEFY)8#?qN&n zl~BJ{iP{gAzMkpJW+M?VbdzOyDdQk*J?(!=3xT-#vkkz|=kUU85Gj+6?Q{r@Ew8MC z#^gb-4;#yc@ODCoO?7>GA4!7#@B6UQ9DFBGxeKTLZAQJe74~b4yOE`U*br^S&0!O{ zs0j@ycy>aNKHreS?E@^F^by1dgw0M7I2#1v!P|UlT;!EE39#D1hKX=FaMh*W`bRJC ziCJkiSem6ahChsZHM3zEW6GQlvjSx>F5yKciY`Mrx3A_4WnQjCa5rB$Mp)?^fu7p^ zDsMDBhf8{WDjE`((Sr>Bm2@g)Rl=mhnlhk^=E?z38DdvkvA!Ot5OzK}b8%hmui2Qs zI^*IRpdF;mzr$2VW*${VyVD<7GrdD@3^$Bly%;dD)=;)XaFe=v-AkqC6>DnJYtmUfgu7I>*zVTA>vJNP3^&6_J0MWic zQPLG9qplv;YeELR`Cxfbw&swR4jG9c@1rJBfN)o3sTUwbCypWS;*0JOfLw8~4kNUf z^dmPcs10V=0|GT1=Z-}tm=l73GzH??{$Nxd2{)eb8!_B)a&+R3Y(%4U%9KY;g*OCp zD758Wavu+A%7A9d32tQz0n_=`%$6i%_3ycuMWmPsi zc(7{zIIb@(>j~8~oO^Yuz97`AC-3(rmVJnW*4#>NO4&I|lBWDGGw>OneJoZGk|heM zeoEF2i0x5N?GnC1M+jF@>9I@CnbE(F>-|zl87#F%30)3NdY0vQg17!I*nUAPrv!%i zBhtLPruO!{5k;vL?bBsJ)26tDS%W82;?N1_n&#jq+xVO4+8XH&;PNK2>h}ejnpe)a zfd!1cT zV<+UM|MlD(0xSpt?teUwV%?6SuuBpg^gf#0PUN?K9jHr$`zZM5;wW#@ym`4iZ+MHm zbO#W=%zhA+_idjBGe#Y8Oh$N_lg(G5ay56AZ5BRv9>dGU9M1crJ%dPMFDqg2~=BX8w}xRBqxBSOIoNcM#P+b?1&Kr>=n??FDIyp72tZr_qJm-j5Q_jn zEQ2)4Li-G)hFO;Fk4+L64lrt1v#5k3&j6_a!n@s35CIU25Qu;z;%&dB5=k}y5s|hX zCyan4W;V60D8@7ZI&mwVI%|$}%C+-80fi(3i%|(6;u!VD4T2;uh{baqQx=9Uz`_C( zqbF<-ZpizH;0!PjS;7s#!+j#7^Ms6$ke4z6Qv2v#;@l7%LJ#5u0rdn2Y(d3f3c(T6 z>OBi)%z_XCFpb^jc6XlZm4^s1F_H-(5mG=9fEqOBW0|I^Hb$n?cCw$z>t_~jy6b>I z5oTTUtvTi#b1%HNzM8lYyg_ul&H0`bl%ElXIL1!RxvkUCzWQ7(op?xHq->|MRm)n- zhsdpC&v)aJ<6~;5)Am}LyVR_9=1Q|bpifw;7 z>18mG{-+P;WI~j2FXE{)xZ*3ki^F4e$%JUZpRnU>RLfRG%%?p*{cl!=47TfdpTjO0 z%Rs|9ooD*Sk3l?Xy1Dupn;I675qQ=_lk6QGHdWkh>oa;-X;VtnSS3FzX)e05NL8h$ zcR2~&O-k_k5EYV$*{LJ~qqNMA)!(2*e`SP&&)$T5wa69Y%bcyDKJfIB%4 zw5;f&5xsw_Xn()7PbH$-A^S85vX{6#I#(|}gnoU)c|QmJn!RCQzTfw-NdSZ)n%_r! zh(^HBZ+W-WF`z^QL~|%vXkcwAWM?pyhiP9-uTnbWKwb|Zjuoa=dq4JRT&DgBQXd87Ykb3aY^vpS!+G9oWj6*f7j8=~z zh8zC8Yf|9UY-_B}kdK~pn?^v+l*?18=w8`|4bIE?Dq6$Vemic>Ga?@Y=R0;O9m0JXL)fGXtRotv#qYFN~n8q_ksVlvs=M&$euua zQi^r+UmW_3PNZ{_dV{YHz#D&KP_n2X;1!bw;|XL5D8ivhDm`JwNnG%wm&@z`#Hu1& z5zaP5zJ$1D2a_N{EXt7w@3yMDZeSOZxo_BIy!s6m76lb6^c&%yZ1kElMNkQM4E}W7 zJr7RiwOR{Y$<4R6ivpls=NjLSq`mO}+Gc`9T&F!j)Ac&1fx3e;-2` z>pz?@-Pk~lp?EAIK)n^r2wgFH8|17ORK1l`CQ>9VCuA!#Aen#t{-U=ai zdlt@$hZL|5Y?O&>G$CklG(|{g(CqR~6v`Tz3u{4Xt%RqysiQD)U(P?=sHmZMMgeLd z11^4Nxtfq1PXJ#&#c^BY^8N0xeTkyA0)vH?8st>P60;*M;b--cVg#P-&A6hHpJ2cl&}Po`S9<4zN34Z4J)fPY_sr070M&IuvFn3H!mCYy?G%xEb$3vM<+ZT@YamQ z98g%X=tO~lf`0gsxirKQcm=zm?q2W%Gj#6uR`Ux&-FrflFya!K{mVah-qohPKhE z2H_aP>OS1ef{Da5CTXV$&KPsTb(l%de*e**=2O+*we`^f`0U46%0{Q8V~ycqGfhbKfF7-cmr zK~PF&PM`DY_15&kczIEbPPB5tjLYb*9R*hBJlSSbq&$YNn90j`C(lGRj|&rRL>H!; z3>NWP@zdFf&_RCY0D{h<3t4DqZAgOEJD%7Bd84o=gmQSKX+wI^t6{tw1WA`5nIoFf z;a_v8gb7)n>SZCZMA{~*i68Hfh?i(==9@Tc!+CZL)ccH;XYjY z5P}QI*%U!a8AllKjzQD$ca`&jLLv*xtEWU3w1_Ql@!eP2J`ifpAFBC@nvj|o!2Ouh z(CNa>)GCOT!)x|SU3@=p7RXGYeVhbAe`;_LTv!M#QPF}6sC?Howbhh1m2-Wcmp1U} zN}NU{Imywm9^RLS@`U1I5e2kJcKOZQO3f2+@(L6|ctq6K8tDnAu4JCm-`rzQM$lwq zvm2e1X~ni@P(gLF2TwX)PQOPW=w=Yo1JndzpyzWtgwccYYNm5B&8((hNH?!@`e2-@Vs>x7-8Jw!~GPMDd+oHiY@NRhin1~H7A)jeQJ1x}I+GC*a*06*{FCMC-zBKcg!V}LmxHptSQ>RS;QXh%wIuBPSIR;zR|c z>WGrs*CG$<1Thf7Iexj1)#rSUYDN0eM0f{d5L9ecSdrZGeiamqKEp;42H!tYFmF~# z4iRI?TZsmC0lH==Z>WONhlB_%(10MYA`52DKZc=DT@IFp@mlGdw9x=fWr&bi$*&!( zi);;{mjf#h!D8-sNG6LD$z-Jn_b`K>`$557>*RA0R@`x&)1Q- z?J9fKR-B9A$lA#Kh&RqYR5?}3T)TM;<6_B8KURFJhRUh}3t&VR+32DRHi#{H!h#Dl zSi+`SrX_x%DU}X4h!H$cA-g%mLEHyV2tb9I;z$e+wGez=^E(~vhsiXFG%>%fnVTa+ zqbr}ifKXfZIm3TYmjfs>1V@Vm*j)AziIIak8NwjiamS-0U)}}q-!T}ftj0hPUdRw# zv|xh!GQ%N;{ea$JnV+ov4dHz*fhj0p2YV0tdeyjCX+5_R%NI}(UBpb+s6{bMQE%i2 zylAIS5=VE3CgXVc9DNjxSJSNInFdSKFl^);(dKHKk)LFS9EbeEYwl7JP%?*%K8`T? z{&w5NP{@&hx)^IldoawmvA8*fQARN`*_Pq@j(Z7)9`bUS43SPRAGrlUx#Po;-_XD& zq3^)SYA;NegCmdh*DK5d2q~j<_IdVL+X^?@y8YRf1ep;0ZS?T6 zT?j6&7O24mH)u>xWsBT5Ddz?n>dDTIj|q)A?eI?bIm!~oRZ#vWc-^7Td(Nzy`dz{o>j!6wpeHjbN`i)7LpO|}h< zrqgYu^Mt{JZG&uSv~y^TSaqHV>evh*cKv5w4Mo`LIbBm}F-1lUTn`y{Z>+_agD)nB zh0M-d4rc^}3v;ft0c(gfJyg*iL#WR58@``YEpi}9r6O}MI}>LV^_;!NY~-ky=<4oSC9y56%w_20p_I-bX3Y3@UGOn-)Z3-zh=R z_Zu$3B6l-g0;AQ|J1PkTfubLhZakVQtY(I0r=a*JRGXD>20*g^(|?!B77L)`%^m*Q z!uJduS2y#;2)4<|vq&3bQD(58u{>)Dw&}C`11{?RFR|4JIZ?PQ`8WYk&rNr)2SY13 zw5R?CRgFLqughKcBh0%a%8I!GmIirvk*#Na9G5>Il8R=gI&t#)7y_%C+#Jb7how2= z7sSFdI}uUo5n(=rgSK@p|CKIs9A!xtEV)so$TL2rjg36Rs%RCzs+!Z)Ovz)MT?K8F zHZ9Pn(g$N-ZhlTP&PTP#Yid}t0%%;VN8{`mo2(PGRJ)XKKF? zby|+qX9iMz`WS|25;XBnL>Rb4R=w(3I~>&=vg`ZpLx(@!?pa+WzmBQnRy{+hnFpyM z_G?M|3kr4P%I`YAsgqPR5dg&CE3ez+geOP^o3#ysn%k zGdjSPH9~)z*#grHy`{3Kvz#m9!%8Wm-A7ivKo=zwwkK##rVbhsX z)Zfcty)KP5Zff#h`Zo8~a41JDPdI!=5%_IPv`<%sDxGmw%Te}jpEyz-LL)Y7l5<(Z1fZ(7Py|Yl}pL1tK9d!XVUvHKI zo6t3U`FnR3P;0yS%a;1nVt!nY6Qmy>;;D&@3RL{QS)y#1BG!MtSApsmtYk4SNUv}j zE(a)5`mKx0ro`9*Pk?sajw0JXZ8zZpP)zzGHyZDOz#Ptou^8>;EAVVf$H0|{E zl2<3p@nuYE=8a@^DVmi_0h%}nT&*>5ZHBdCOcfZ+j>im^r2iFi5{+x? zSHg69VX&9yc{y#ifws@G(|+b4vnOeT-=I5%UNdv5p;kQuxIVY7V-SzX-}UVQa6?&} z2Jnl+Y}i}isGGcZqrC&p3JF}4-x}5W7a65#fI7Lk#(AxtFFmiw7r@h!CltSBv&%0K zu5cy9At5wTpnN|G+H#3>4PI}raj0GxgHGeP1^x2v5@iMpbZtYPY^VB54l$dtx0anN z_H#Gd{zfXNYB-Sv25Q7}4I@|2mAEQXzql22o7YIL-SmM6y~$jkKMHJqC)~@}w}qAR z|9uRKLPZh4_tzNJl~fBkDUNyg4$g`7b-!iEU;UbfQ9gUKBH41>6(ez@k`zq27_D|qKww&)T#WK<7L^j zlS#8Ru+XnbC0x1t0hkOb@-QY3kl%C$E@ZDb<+#WQ{hGkn1besnfk_Ba*S%XM*MK+% zFSehmN-QFFkqB)w+u;KHxhvO|x6Bh#9ccO01m_h<52c80jaZ^=XgZM3k7q2C30gub z!hEfT^w4Nen5*DBbTTO<@o{Gv>;t;c{Fr7%brlcc-){^q#6JM_FBrNfsV%lg%L)9! zF*rFhOUOJSfhMuMplp0Xzi#^>Oo97^j&wk7wck0&eUj9xtVv)p_qs*Uei@zM%nxv4 z!EG{0#~*j0U*Pa)n?lQsIi$JhRe}4q(?+Ecl(PMe#8_7=R17H2NT0~cmYRED@7nAZ zU%gKc<6<#m%!MUIu<;V6g7n-#9R;Lcj_2Z}6;B>O)~LtJqN(JhBzz(l^w0#ShpjK; z9YS&ZFk$nrX1Ee-_f>=&^nY1TK0oqIlZ}&3fnC4(GhR+Tza`y-1F$S%mRe323MR#c9oryglp9%={0`zUSy1 zxN5kfhAS4r(8^-eU{0(;mkMsL4em}LlwOP78&$^_FgO%J)cWI!y#eNI+5^DBZTpUO z&D~QSNPnsl;$}GM0Eh_3SMzJ;u>Mg6iYhEzYxJ8hkIV5La%bE~rUjS?Gnu(sd$MOU zH8^l=%+_kHi$=*J*$@ZcOWxHcAJzO~)Eb#kBv!f;ULGu@{jbS=t zbN)br&?trXE=B14vvOMSS^>{Pi@|8=(==$I2! zd#qx49+~Qz(&rs#_Uta5pUWem{a3nQ6`}R^W3z}!2c=N>_XDaezfqU>DeI@Ns+;Wj z>W`^7o1uuc?3G<%_-rBP7hzdfR+7&&bMf~V0o(_SY=-Nt8fJd#WAZ%_Rsja7 zq@HKO=T+7~#NW3E@(f-^8}Yu187e>@mOdhjyR@2xme)Uk82;W3E8U?e@bTaSJ* z?k##RO^`p{9VuZZjcp-Vj-VxrnFAZ`#rSSClF|n3*Qg$Sqn&(%?U}xu66`4Gnp7(z z^H$}{8*-s-Z;J0XlDioGc$?Q!v2Ha`=<^^X53-&el!ZWUbcYBYnFrXm+ z18rG4HbXlM2SdM|hE|N*Pn?rkK-&8@UA9cBG*0NZ85?L%>y;6X%P}rXM#ly$A2+1#5mtsd5g^9(qKISn0jY0rPD4&3nE9+J~chNCSO>_MKHO;3F{7k zOl&^FL5s*l)OC`OjG=jtm0i+nhVVgV`Z8 zVGk@oA0ZE`^(_cns$mpIlFHS2N`Z1(WU#X%K9Ndz!P)`=-2D?x8zf zVIcDuHo>CJc=r;B>&~it+{Yd_SIIbj_E+w-^2aY*I5n*sY&c}y_o+ob?92@9i)AJD z2sx7SWl4e6YB=oAhBI(JDWtZu4>W)*>=$zegq**=NPl4$Kle&Z?;*j%-=A@RY#s01 zN-l7H6UYkb-FKhbbGoJKLJUl7hDj)v4BjUCnK0wE*g5Psh-ZdPP8#BV> z-ww5lk=WcGFP088!@4KU4^TDR175A3gQ$x<9nlvVk{uw?;C)nEc+aZzj4rZUjHwz} zUI^HWd)1HZ`Hr?uA4)_V<_4NnN?CNgznzf!jdkg-3=fn#u|>4^f{4;SfX6xT)=qzh zhn8u3qal z21JB-P%Ul#K!Et%!3D!&IeI+b`+?3J7Pheqtch!oJ)I2Gwbg89hzOH@aassf`SN|cqANx6CocznIYJ_*#?z|)^k0i!SmQfOXPioKcWWt9LWzvLSTXonIG99 zH)IVlZIO=1G96TPFc!iTEs^3Pp=t_poy=_+)x(=vT;R&X>)mAHcU*azK)#X<{hQ>N z^wY$hBm-5D<4g`GdPsLL-=QCV6fablB4_Ud{a^n{G-^n5&@H@$&zhn_UNcwt>ImS@ zCN7FVQ9{~-uq1~X129K`LL5e)E_KK21?WO``Lm3?as$DD`oj$_l!DU&6j1eTJPkU` zf(^VH;-6;CGH5rpm&#P9ca8$hYeFZz948%inyQ6#0l63g-a;G`5iFj1O)21ZA#9+W zqMV4^Lt1EJ4SVjEx4g%JU$>WrTsSQuvn}!gV&NU#ZaUXwU^fObvYp zD&5JUaK1d5fuJBMaFRtr+;I25M)zx%e?t{Qd~- z`vh+g0BtssO)&^*q#=xKV8+JA#>U3R#=(qiU;u8Js2cn{25&fYCfA@aiZycThbLA+ zn7pGA;Xu_W_CUxk)f|H2yC9Y|X^Q?59@2EZKxS*cp`k6|H7Q!SX*mW?-zC7EPlv#e zqMo62pySh)h@#N>hNPG-fEmg2l5(sPi>c!=0X>k%kc=+o`Z6#JZD4T%$ZM)`Ra~Xi zd2j~TmcXZ z1OP+|q#&@3Hims=u}80Wo?+F|Qs*WG16PDD8%BVUv#Ze<=fR9wq9MZ6X8dHN!i9uS z;DIe-qc!-YF@2DzV#H&Zk#n^SP9m>DOMDa!+~5}G_Z`G3!9q|?Z)jTNoVG`pd?B+g~IcDHh) z8Hdn*ph0-Zg6)9?uZk+tgncxLdK~x{fTUGv1H3SnuT#o+1S(*ZP(%ukd z0{cM-A}h`0kuoF@TZo8-AO;D5(-c7Tdyy?wR<%F-B)T_r(_JO^+EraImr(9f;?Yp7 zsaRG`SWyKE%H-9`!ns~6)ynlrrlPzkgDx3yL@p?xgei0l3N(? z93hTiv)E-zWlUw+Jmw6PV##SKkcAXcK}Cy-weJhW@eN~kc^x&4;CBqO$S}w;ye}!F z0hkC9)T6yM{VOy~6AA)@Q{ex+2oN^2SJ9#Y1FHF+g`i<%2-pAvlj~`M51tW_1p_sS SVZT9t{9VZu;X*9x2LVm+p6mB?cK+%I{Vw+zE>Md=m02FKxxz5+P$sp z=Ye4*QW_=@!ZgqVp3)Jh{HohjV8oQ36oRPG%-!7>Sm@0 z!zP*yHBWkhH1yJ*$QcyW`bieuCa zspUOAGA3$cYMVlO22*76CQ5ruRQ8)nk0UZ<#Ka~OYEK$vHjzCO(OW0SO-(216ZA|PrZJ&FX-y3gpbbn20g|3iCaL)iQ!zbHN;LG*w3(#H zYI=_nY?DW%8flPd$)TWnfDH{aG}A*ygH5C~G}@k{O#{>hnl#!+s5Bm?hJ$Ko4^Si| zh-gH}nm}Y}>S**6)X+Uaf*Lg1Ks01A4K&k7q|g8jG7N@>Y6gL|8Udiuhz3oh(8z5- zG%{ppJws{@Gynh*RDx*K)jg3sB*2(c$>lK9MD!+UCxr5+lV}QhO`>|9nre6ik4dl! zev=xfv?+}=(@jrLMAH=TnHi#F3FSPDnlfs7dMT&1GHH;|pP-Wy(X^ROlnA7djUY;T zZ8X&ICRF_zYJQa3Wih5kihD+s{70&J7$ym-c+ovDQ^Hf(15AcVup>;F8mFhJ=*Vf3 zYI&xK`ks>%^)`vJXa<^-@?>c|qtg`iJR*5AWl!<|JXjz>fK?zsbyb@!#<%Nbuw@@2 zNAc0yrMa_Acgr0^<#Z=>+IQN=&g9KnIBBepPjbd*r)>;ZXVX9YL{XL;xZjz{){>he*8T z7RU9u1PM>)N$H2$3D7{39X$)gz8}~G32E>;$+_&6coxD4ECTdO3NQ%-3vK++#4`wp z-KZ*}-QuE%1XMv{#1t0pC@t0~00f!9g3bgMU|G7L1we)Z1qG6Vj3A2ukb(>0^1Mh; z7V-lA@PSA|0mFz$1^_dG0}v4l0TMia38QWfIZTbm2%~Uy)`tv%bbfT)lJg&*VYusq z0R6?+%Gg};VwpI95YMvd!5uFYM{^7HTQRx=gYF;)IseLc{2nZ70huoj*>cJqHKFI=c+F6sR8UOI8; z>_$Qt5;QfCmWMujd5s(ocPPDHF3F&B-Tc=^vZ$=#2%<0DO2w0YW?^L1JjFL;2*0ETsjk5L#@q2rS?bSpo}6sDe^LfCh>IQBcL8C@cVi zDpg=6O*8G}X2eqA|WsBs3IaB7)bT zLeE+y7@8+{$2sDKvh)z%aFx59I zw;~Y*)EG^fi3Kh9osf7PaM1)8HwX#~gN0q;gBY7Y#x{o18ygTYj9}Q<3}b1uHZ}m* z4WSL7O^v3|inft}5M1TLp)+>P;B^ll$*(c#GjN&fcUy}Q1QfLqaDZB|0EmHTJ*iyX zEUT;rAmTZ%L@v#Vq38iN2tLE8@OH%H04M+(U(9QsGmO}ki*pEq(E-El4-CI29ITEY z{~2!Y5@g*M4$D;PxFdPyirRvYrW4zxgExeYl%VKL7~;x;t^5_YOUq23MXE!F7SrM> z5EK%1L>RdqOitgu4#U?1m@JyFP1Fn_Xi1d;&`;UJcFE<%7ql`D#YHP9g5}i_+#c~T zCXx{a;!(8!;8+d5Wqdve!;kq_EA6aYzEJo`@tTFMF@l~OOFf%zb^nFsoO2A3;g-G! z+*E>GTYbn7B?FdP^MV+%F@?s>WG5vv$as+JT5jp}dI#bA2t(SX_d(oGVnCHTY-e$A z?16if3gp(6BXUF6jITWoGlWnUwx|s_EC`|`D6vf35jNO;uknZi3s`1C_*Kv`wsK$8 zG?a7a$IoX>2ztjv>jwEZ7&F1R68Hm$^>K6-swRMxiA^b@#EJ{v7>K`up14i{*WbH#8pUOV|$!rrwI%Q$%tbAos3< zzsjtZ0GR&gL2wTY0K0i3z?O-C%`jg7fw4~(1#Y8+ZO(wrt(?~_VY{5 zo)cvD$GzJ>m%1PMr|8E@akBymuPOFCNgd|uf)v#%i1%e1-7$MG3-*xI)ssq?APl*r zF_QQ>^Xoy4K0@MaB48=A9)1=?*`Cqy zDrKyl9!_AG9YZ(_<=YjjDY-;lGlRpR#eV$2NT`@Fgty}bm^#4@Q1?hhI&f`%lOR@Z z2J0qMH8kS9OActIDZ9I=c7T;4+u@hex3^q#`PBviN#AL1XcN}*Clo%yEV)<`p^jz~ zVWD1{h%p2S#g#Hvz+DG**ds}>{jxU!m)e!=A<02TKmtV)+XrjNa7E`GRpCIiaSRzl zI_Y;}xG(U1vqmYoAqS#Yi4Vk(F6M59WR-we!V6SP0SHZDA~CjKJUFCDX9bZGIwnMNh)wbE%LbphH-F~SgHqgFis2dqy~hop{;js z^h_i)2Lm7gV)A0X6(y)=I}!y0DY4L@Xn2(^|DDsc7cx# zr%5hDrkz5y^Pdq-NC#x)OkIoFT6B}Rc;)~*Y?|q+(rt!8!AC4hXE&bBmv`TX{rZ70-OU$1+E3=Taxs4zt z(d=nLUBF@)FreVA>8h4+8oU1+7^7yeqUsLn`sqxZQIgDY z;}a^$0t>P;=z{br>iHPZ5Lnh4jldwcP(fJ%|7$tq;Me@Zzz;^IWAzHj+CnkV`qKfK zjH5@Dh{s3tUSQ`pzM18&s$A-O@^>7rGPyC2mx`ZSca9|}>9i$|FMOJL9s@M=b8;G3 zWPU=Nsrw!26is!nM+)9^u(ulzAK}gKJB@7@1j9!OBVrpsLD>}d8IZ)sT&NIPtfm2X zCjK`d+KQuT%i>j!&r+eiH}{#+?kUecv}lpKUbma&PGKhGyHVoNlw4yAAQ?5RjN%ur z$(hEoA6vHaRrYe8iFZ#emuwmM-<17dz>P9)!qHyag1hGCdZ> zd4ZZg8&clW7WV>HYM{$a={TT1t()Wk`Vi{Mjl!2wYDsKYBq)X^P}`BgT+44SHwv@48n6Vg=&$|ry?m5KbaqBOx#=38}xp`dU-&?ux_Ik4br+G`@KTEz#< z#B1h?jXr+sT|?-3$f$I~wuWBE7lfO1O;5M8F%N3i&Mp=2%YBIz3HrCv)MQROsK{h4 zt<(cJ2shx&hoq>CZtIc2Rg{h?$?SxReHA~5DQ`uJFIHp7jmYse4jH30_hn14G-GDY zcOY{WI1cH69?O&n>GysOUl|o#?C_$K#<4e;lO^}uD9mZ6k9vGqQ+ri3FD+T6tTz+a zD_pE`p7>05*rxC+4|Z4T_#$S(Ai{+n^`3t`FQIC|?-U#tu}U{;6t>=BzkecV^vJ>v z*Yo2XWL7*KOY168+-J%mW#H~-A`f3$AQTZ$7wB0)vu50#0M{A-6#}1%`Z%1y1-2#Q z0kmnY{02-U3nz&{2E0PWk@||R7=i-NtDPNCC8(>*v(j%rTpijoAM)CJ#VY=^(QU$Dr0e z9ToQ4;f}VF?QU0-f2Vkr{Xlu#rt%mNvFZUzYspppuXw)aVWt}6c5I)roIB2Ss4s_H zidJRg6MTz=0B5Wu`2Cl&LnY%Zg5JukhqDaoz15W{e;JO^mvnbOBEXu-nqpt#cm@G1 za~1Tos?bSeORo)7QpHl5Ep7h#iVB)0|Z?RNU7t5 zt*(a&2@LrP%*`oQPFNSPqpKdCYFW9gPR#7p?@65J?${Fd@us%gCoG>3d)KUJPBN9V z(RQ_fb}0|cOY6(E<;QX$A?sAC9;3)+kxn;oi% zLMLLq(TEv#=?vN2Q)AR$0x{DB74kGzf8!vP025P2gF(2N!A zByU87^T-8Hf7H0eI6!>_9H1Zypn%8_ehdPzh~xj(!qt)?2m)ZG*7BcAfa3cNs5KQ9 zE-yF~2;B_XyhmS`9((x7k2e#clWT1cKcS_kr=|!J*?m{Dgz`sf-OZNv8vXOl0*WtG z6ctbiP4}N+}4WB9VX~pg@%hhNL2GHkj0E4LKABpLzX$Qp(Le$71g?gd~KJkiLqek|IJx zupmmqkYgJ@iwJX_T=0?82;0!%=XHk%-OcIQCA`%t&u@I5$&giy=W393KV=du>Q}A( zONz=r#mhkBAB^%;Oef0>1gUWj4^Cie-PT`&EZ`hQLh+5Hm+AwyFb18fL?SF<;Jc(a z2NHd5H!@S3qvQH)uXnR_3SU%HqTC@9yuU~#2*AmhfJz|(S#}!jm9^K#-sjWvRc6~8 z*GabR^;o#tt=C_^WYfJ`2^j>K+g9q9&%q5m%}WvhgFpxZDNi!6;?rC^tEs|tF`&t* z-bl=-D3ro;YM8A-w5O9O<)j{(=h!kT;>2K!*$02QV8U;K2|-?fiNnq4SQ)Kigd4js zc!gRL^(GR~;XyMg47(a6@>}b+wysQi>nWsakEKQIqXE;`7^;O~0h|o|C1W^d&Z(gN zXL%HF?#hS%ec#>oJ1aJ?;+zn_9DYlDgo=>QD*N;nXGnjPVj< zi{kaGpSETbXhf-hy@&ORA6^jt$@r$UJ4d%Z;A~jt&Q?OZX|)@Wlik}$X5n7bK-hkYF6u`vv9S+25<&2C z$wF7joYX7AkyHwSni4A=YDRr)$;o*8rkmjiyTs7VJb!wMOi5YSFgNi6}YDK+k`1+IHN9v;Y-EzN!3 zZCTkvaENu7!xMq3Dx){zOyPyPqdjkh3YpeknL0K_u!DBS*n7fpjGsRGP~eWfoDg(` z&1LbM7+Nz=7LLi!GX&vLES#<36(F{!bnLjv%*5jr{S8-oAe@1>H!BFb6YyMXpzH10 zCZxfOP-Jve`cnV%jeq6!we<3kI2KD7AkRK4a0oAKUj7-ViDEVzC-Uo!QNYze2rpZFR)X3b zCmPoN@_vG^_H1x(#VUK1q!3!-y3|2ysDkKzz8#6pa)a zc=qDm=^-VVK?UjgU?9E-2rQGYbwGma-+L#8D%J&jJSUku40R3v0u~wLRa!rH19lVB z0EjJqL$~Lj^3*W|j)*J?w>I#E^gjj?Q@I;tG}mLHS&k=48ex1C5MQccsprm+xOx-+ zkaA`6x5v(wblhF8xbEa_6Z8dHJ#3FR4CzooWJwKs))~ea;^cz&4L#{tK#G}n3^#aS zR$h)n+#)8P3M(K8E&l115e0+K)fr$xZMH|aJ#NJ01!OH9G=o30-E7SxJTUsmA}eR_ zx&-{K`ayU#`yl3##9R0vTVdf}4k4#?+(KDGtolU{@JR~-(|68{$iF55&y^m{x^?!L_8unasxAdshWU)9jYCU=FjpFp#*GWL$$vEg4kn7 zaQY-7l&GWF1|JCQNXD*EE(Ke_SbdaEmNv}L{Oh)qL=%C|F)5I%ealSKM=9Dh zxA6!PacMCjf_|`=vZ^DwkR4eG%)n(pi~LkL-le7P0=!VxzC6F@oFrk+!aOREGH=KzAj2rYXk*sn^0<}=lM6-(tk z9N_>@KFUZqRs*w`Bn+4t976%Vj+}b2C<(S@9W)I|-}>wgAZev{z#vOI@?yM50R@kT zD*8Upcljz6T2dFOiJBqjpk*7!LnsqDpO`4yGr!OC0!`y~LAb1dg4huSo7ExI^tQ-CJXU!EEwB(;nSlk;P#Mg8 ztPR2wx#61<O60m*Ou1}QCy)fNH4P`xbNzlk+jl# zGM=nkfdF3Sf_nG2o{pAXJhNpnQlzA-w^FjgSeUZIUz}2v3 zqcxr98nj7Px=q3vq67gJMfq}IAO^1KkVc$nDfJrZgWCS`Bx51(P8^Gw6}S$n*^}8$ z*&*6nW`CYqT$mdL*Xrs2oiRYj!FnLPN7d$%udshDIhJNhH!@ooEx)!!)$(4@n5e$H zGW>k~&BV~oY2&$X(xF0x^i3>5ReRboiE8)asShgA81jGjjW42^flrJWR*cNdCR?40 z`fO*;go&ejDV_n(gx-3|!;(!{_u~)QIt44;LC92F74@q?Hx06!3GtrhpqZvOufW@rs?E&mn?a+? z&!8Uf=y1f`Uiw&wHyLBpCgftgT$a(FjH*Sjv0_1nn(3Yn0(+Ecw~~fBuTXM|b+rNh z#6KIA+Po_t7HUtdbUpYoenp_5_G~9E?yVsT+<(z<&#O~#*{36hvbxIVzf~ik^zpL! zzEa_AJkwGSyRoSObgW}4)v==tP@};t>T-7l3VG_sW;Op}_5}EawgDov~0{T=g7( z6rsG^c^@d*Oap7c5;rnFKg7K`qeN`=0ix!v3YWb_P*S@FdIq^>75P5%A{7eU1**;t zRJ{KOrbHx!FCzo0Y(+6NpiLg{kvJbxhn72B1$?!T`X`!yFGAlEs=mFs%h$TCQ9;Y- zL>@P15{ZetZ60>x`a{%gFZRHcm}4|uH>Qja-S8V@opEIytvN9xEPEkrvq1JvETxLE zv2XEnT;ushrEv1xc(lWID5Su;yDd`~o$p8p){#J*B=JF1tt>i1akyNti@d{L=M z)63m8!LeyZG4AM4?QhQTcj|72q@fS|7Da6}BfSZ{B4I_kT2>*cq2iq4%0GEZ_SH`p zR4e~SDq`s6*RbL(`^xSiBo&#>l7_sCGOwplsn6$L7Y_V?sC3-= zf374B#F;z&6n@e@^b#!vcGt=B9BA2ES_AMqt)W-wPm8JRce_v_lQB5VYN^T_ z?Djr|9V2Nn&fM7C-KD;0j^HF#-Wk-e2F6#5dIXPfJ*Qya_M<|ktf#`-9@fV@4uA3U zifUO~5%`w^r0}*o=P&$ePzFW?23+>bt}&uHwM$wCfe8V-0qP64zC2{%`fZtl_7&RqRhI3+(Mf(K4Ri_#<(@!C`Dr0c<0&#?4@U3)tt-W z(ldqyJoXO`w6-r{r8)_|H{J}e4#<}6GYmszxeBYzG*67ZK02wt_h-r$JHgM()JSHd z{KOD#ij(ulM`xFz*sFh<3O*Wi2t9P59$p1wvGa6Q>C5P8FBe_O%e|-4TwGFGs0G*H zKLV!fUhJDAhZbJMx5QqDSPVC%rkEr^Uf1yc9PHXzIwYrL8-dq=U+z5Z~d&5 z&zMYEp_JOC6nzevF<pvR7+&>hK;?U9zTDZS z*KdVV$xHM5M~QQ*hCot5Gnu0Uv&I_I7#U|A7%2{gg;5*rSHbBuHPcc{eN4}uW#pLi zSN8VhbmGT&rH4Wqq&D7!@4!Mf4naNS?jgCEAi$5zK zl1m^yhAuyz=7f+p#q`WWAgc*?&Bu8=mwkX-7KDk6Cdb3WKwEQh3h2l_aztKj1E0+v zg5S!$2j913Ji7gb(HN~m5r2{V%9K7AL*fDrP%SOFEp5kE4YvIPH?mgae)CMkF^3$} z9woo^1pXR7F5@Q%6X-WugN`pZq7N6$4=nz5K7Np=<~jP4P-?CB6VY<6cTH&c2CHqI zaaL)jt8${9d=nT2=w!`IrXh$%?8Es`lHise)TW(npH%2{Ip0ZFLviA{shpVqexpuc zJdJNzx9m}I)i)5?v|QJT1{Xrvu<9NNHw~Tc-4IoyFKYCHQ{dLLrq{r_+>Kll06>lSvtKH+9}N62s`k1aAM{BC ziCjk1p@1_0&Mo&oPcMh~!%R_8V&kOVKUc2hc&(ym%1EXKma~$Mq)eN+o4zvG_3F7@=R+* zX)bgU+XjOMU+p?GBJE}aW(>aB!!6mTWx|<`ooiYd#sd&+y~MRGnX3q_?AS(%c?LWcz;md+z`RaVo^jd83C z_OgvOjm&W{I!1c7siK;lWQxLzG}JBpsPgQaSPWgf)oZffDivF#n#iNamRg7utpqkz^nZte8#VjVa2j_OudCp}e)N z(-EyVe}TOhXFJywG_95%q%+=afZ&G4gJvr3n$tu-OU3-|8})%uV{Lc&b1uLTh>+p; z(4j<&JvY+{YJr1w(>@~yNZCma4yF6&&l?qutrDOh+uzU^IW{7fNK+C7BPe>aGektm z-@*~DfXL{-7uqs4j8{cU6m`DkMCY*%fYi)4+`fsJbTMKmnIhRK*nFr<}IV@OD-<2e%}p0S42-!!d* z@b3@Kbg#`J1fZVRDf_NPyjCenA8Hh)}olLe6Zk1_GiR|1)& zZt-&cOpjzbi3Rst^GvkmaMo8!+dZ|#;_P3#(?|Wv4VwhZ51?_T#OIE5C~4IZjcztv zrpn3qa4e0PLiYr`%Zdb&Jb;O(k@PhU+{fzOQ+JJp*TIGhq#6fH!12`~fXoCcDOy z%}3%86#SOs=55LOdto3uW;}8b-VkDP29JR``=b{L{!KX`&oF!84}& z|52=mu#WQOP-jP7JavvSvy*MT8V?$EYhNYQ5OYPm#ghfkK`tx4mv;l1Gra3CmF5e1 z;-M$Gx(3LYaRhW+`7Z7n5bjH>xh6U&av{tf*5fZ@ymz$W44IkC@_|Sg`{l)gSn##j zxO{ivc1)YLT-yB+axaNqp=92;;_c->Rxs*8uM*Bno{QKGxhD@zBT7P6@Dh4L4?+^v zIkPiWM@T`>Hj(4LbLdb0*-}BoV9QXF7C0-XBP%oJzpP9qQy9#{7`J8QsCLA~K5{Bdg`-JE$obPvM|-)qkQ^2zMtdE$_M zBoCRpkn3HW92}7$_~09)2fGR#EtjAs z%z?tXrpFhv4hup0_m?A>a1;EJ11SOi z0xsH;81LS0PXt=CF z0+4>rIW#1MDXrDaVM-*B5wDconM#5xkGtf-Nhukr%82XW5OssY&9ejj7zNZ88<<`K zaB69B4@gc6Pf+rv4Cw)Rgj6R(uWQq=)~^8^ol;$7FjpT4He$XA+WKrB$7V*o!bvSI zZ3-!X18KC9X^2BjAq-<<1~xV}Ha0dkHVk8900X&bplhY@dOXSBuNQ(_*E(?TQ^$}< z=K9Jd=$;e9ADewzD0VD9p2l4?59#y^_RhLTZp|e08~)(}$$Eqns}zKeDU@&Q9j)O^ zG$6nx^aphBy|X$qK*$T9W=J6ZmnncPT)9^Q)?~MrO-rw$1}yaS#d!NsiDni+k=0Kv zMRXK~=)dJS6&YXg>kHn3e}NA<1Q%tf*}a|}?-~Gv0ZD$w&{&mP{>-(Ji2Ec|E>mqt zS?kqd?1RJ_N45JLw=gIdr~)7s2mpu`NI_sma)$6n54#tom6UN>v7PKsouw>bQj?iL zVgniRnm^R`w18+OL|x z!I@3E=M*L_fHmC>#_#M3w_vrkj^O$kO~5?el?C;qwNnGqnsQxno=Gf=Z~8mz1uxPF zFX6#}5Lgrijwq{X2s!%X<&Jy{n0izSu3$9fR=+gu$gBck;SfErAOTL4K6xMusR%(4 zTLu*qJV6DVh=@`EfM5tejZ6z`gqKYUev2j+t<9QxItR0?LOuTcBqeD@E&zk_gvtN_ diff --git a/tests/testthat/test-surv-royston.R b/tests/testthat/test-surv-royston.R new file mode 100644 index 00000000..2114718e --- /dev/null +++ b/tests/testthat/test-surv-royston.R @@ -0,0 +1,25 @@ +test_that("comparison test with survival::royston", { + lung_data <- survival::lung |> + dplyr::select(time, status, age, sex, ph.ecog) + reference_fit <- survival::coxph( + survival::Surv(time, status) ~ age + sex + ph.ecog, + data = lung_data + ) + royston_ref <- survival::royston(reference_fit) + + lung_surv <- data_lung_surv() + + res <- royston_survival( + data = lung_surv, + truth = surv_obj, + estimate = .pred_linear_pred + ) + + expect_equal( + res[[".estimate"]], + royston_ref["R.D"], + ignore_attr = TRUE, + tolerance = 1e-4 + ) +}) + From b8ab312c2ef232b2ff6158257e5236b3feae7c4b Mon Sep 17 00:00:00 2001 From: Hannah Frick Date: Mon, 15 Dec 2025 11:16:14 +0000 Subject: [PATCH 3/9] add case weights for Royston metric --- R/surv-royston.R | 26 ++++++++--- tests/testthat/test-surv-royston.R | 72 ++++++++++++++++++++++++++++++ 2 files changed, 92 insertions(+), 6 deletions(-) diff --git a/R/surv-royston.R b/R/surv-royston.R index 196464a6..7d924ab9 100644 --- a/R/surv-royston.R +++ b/R/surv-royston.R @@ -90,19 +90,33 @@ royston_survival_vec <- function( } royston_survival_impl <- function(truth, estimate, case_weights) { - bns <- normal_score_blom(estimate) + if (is.null(case_weights)) { + case_weights <- rep(1, length(estimate)) + } else { + case_weights <- vec_cast(case_weights, to = double()) + } + + bns <- normal_score_blom(estimate, case_weights) - fit <- survival::coxph(truth ~ bns) + fit <- survival::coxph(truth ~ bns, weights = case_weights) est <- unname(coef(fit)) est^2 / (est^2 + pi^2 / 6) } -normal_score_blom <- function(x) { - tibble::tibble(x = x) %>% +normal_score_blom <- function(x, case_weights) { + tibble::tibble( + .row = rep(seq_along(x), times = case_weights), + x = rep(x, times = case_weights), + ) |> dplyr::mutate( x_first = rank(.data$x, ties.method = "first"), + # does not need kappa (from Royston & Sauerbrei 2004) because it'll + # "disappear" into the baseline hazard of the Cox model in + # `royston_survival_impl()` z = qnorm((.data$x_first - 3 / 8) / (dplyr::n() + 0.25)) - ) %>% - dplyr::mutate(s = mean(.data$z), .by = "x") %>% + ) |> + # average over ties + dplyr::mutate(s = mean(.data$z), .by = "x") |> + dplyr::slice(1, .by = .row) |> dplyr::pull("s") } diff --git a/tests/testthat/test-surv-royston.R b/tests/testthat/test-surv-royston.R index 2114718e..dc68bab8 100644 --- a/tests/testthat/test-surv-royston.R +++ b/tests/testthat/test-surv-royston.R @@ -23,3 +23,75 @@ test_that("comparison test with survival::royston", { ) }) +test_that("`normal_score_blom()` works with case weights", { + # weights without ties + x <- 1:10 + 0.5 + case_weights <- rep(2, 10) + nsb <- normal_score_blom(x, case_weights) + expect_length(nsb, 10) + expect_equal( + nsb[1], + mean(qnorm((1:2 - 3 / 8) / (sum(case_weights) + 0.25))) + ) + + # weights and ties + x <- c(x, x[1:5], x[1:3]) + case_weights <- c(case_weights, rep(1, 8)) + nsb <- normal_score_blom(x, case_weights) + expect_length(nsb, 18) + expect_equal( + nsb[1], + mean(qnorm((1:4 - 3 / 8) / (sum(case_weights) + 0.25))) + ) +}) + +test_that("case weights works with equal weights", { + lung_surv <- data_lung_surv() + lung_surv$wts <- rep(1, nrow(lung_surv)) + + res <- royston_survival( + data = lung_surv, + truth = surv_obj, + estimate = .pred_linear_pred + ) + + res_wts <- royston_survival( + data = lung_surv, + truth = surv_obj, + estimate = .pred_time, + case_weights = wts + ) + + expect_equal( + res[[".estimate"]], + res_wts[[".estimate"]] + ) +}) + +test_that("works with hardhat case weights", { + lung_surv <- data_lung_surv() + lung_surv$case_wts <- rep(2, nrow(lung_surv)) + + df <- lung_surv + + df$imp_wgt <- hardhat::importance_weights(lung_surv$case_wts) + df$freq_wgt <- hardhat::frequency_weights(lung_surv$case_wts) + + expect_no_error( + royston_survival( + df, + truth = surv_obj, + .pred_time, + case_weights = imp_wgt + ) + ) + + expect_no_error( + royston_survival( + df, + truth = surv_obj, + .pred_time, + case_weights = freq_wgt + ) + ) +}) From adfecb5f854b7e83039908da644fa87053cfe6cc Mon Sep 17 00:00:00 2001 From: Hannah Frick Date: Mon, 15 Dec 2025 11:22:30 +0000 Subject: [PATCH 4/9] polish --- R/surv-royston.R | 6 ++++-- man/royston_survival.Rd | 20 ++++++++++++++++++++ 2 files changed, 24 insertions(+), 2 deletions(-) diff --git a/R/surv-royston.R b/R/surv-royston.R index 7d924ab9..651327e8 100644 --- a/R/surv-royston.R +++ b/R/surv-royston.R @@ -46,6 +46,7 @@ royston_survival <- new_linear_pred_survival_metric( direction = "maximize" ) +#' @rdname royston_survival #' @export royston_survival.data.frame <- function( data, @@ -66,6 +67,7 @@ royston_survival.data.frame <- function( ) } +#' @rdname royston_survival #' @export royston_survival_vec <- function( truth, @@ -99,7 +101,7 @@ royston_survival_impl <- function(truth, estimate, case_weights) { bns <- normal_score_blom(estimate, case_weights) fit <- survival::coxph(truth ~ bns, weights = case_weights) - est <- unname(coef(fit)) + est <- unname(stats::coef(fit)) est^2 / (est^2 + pi^2 / 6) } @@ -113,7 +115,7 @@ normal_score_blom <- function(x, case_weights) { # does not need kappa (from Royston & Sauerbrei 2004) because it'll # "disappear" into the baseline hazard of the Cox model in # `royston_survival_impl()` - z = qnorm((.data$x_first - 3 / 8) / (dplyr::n() + 0.25)) + z = stats::qnorm((.data$x_first - 3 / 8) / (dplyr::n() + 0.25)) ) |> # average over ties dplyr::mutate(s = mean(.data$z), .by = "x") |> diff --git a/man/royston_survival.Rd b/man/royston_survival.Rd index 1b703abf..f5b8b1ba 100644 --- a/man/royston_survival.Rd +++ b/man/royston_survival.Rd @@ -2,9 +2,15 @@ % Please edit documentation in R/surv-royston.R \name{royston_survival} \alias{royston_survival} +\alias{royston_survival.data.frame} +\alias{royston_survival_vec} \title{Royston-Sauerbei D statistic} \usage{ royston_survival(data, ...) + +\method{royston_survival}{data.frame}(data, truth, estimate, na_rm = TRUE, case_weights = NULL, ...) + +royston_survival_vec(truth, estimate, na_rm = TRUE, case_weights = NULL, ...) } \arguments{ \item{data}{A \code{data.frame} containing the columns specified by \code{truth} and @@ -12,11 +18,25 @@ royston_survival(data, ...) \item{...}{Currently not used.} +\item{truth}{The column identifier for the true survival result (that +is created using \code{\link[survival:Surv]{survival::Surv()}}.). This should be an unquoted column name +although this argument is passed by expression and supports +\link[rlang:topic-inject]{quasiquotation} (you can unquote column names). For +\verb{_vec()} functions, an \code{\link[survival:Surv]{survival::Surv()}} object.} + \item{estimate}{The column identifier for the predicted linear predictor, this should be a numeric variable. This should be an unquoted column name although this argument is passed by expression and supports \link[rlang:topic-inject]{quasiquotation} (you can unquote column names). For \verb{_vec()} functions, a numeric vector.} + +\item{na_rm}{A \code{logical} value indicating whether \code{NA} +values should be stripped before the computation proceeds.} + +\item{case_weights}{The optional column identifier for case weights. +This should be an unquoted column name that evaluates to a numeric column +in \code{data}. For \verb{_vec()} functions, a numeric vector, +\code{\link[hardhat:importance_weights]{hardhat::importance_weights()}}, or \code{\link[hardhat:frequency_weights]{hardhat::frequency_weights()}}.} } \value{ A \code{tibble} with columns \code{.metric}, \code{.estimator}, From f1d90ea71f7021dcce580a11f1ec7d7373eecea1 Mon Sep 17 00:00:00 2001 From: Hannah Frick Date: Mon, 15 Dec 2025 14:01:38 +0000 Subject: [PATCH 5/9] add `royston_survival()` to pkgdown index --- _pkgdown.yml | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/_pkgdown.yml b/_pkgdown.yml index 2d123afa..c323e164 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -90,6 +90,10 @@ reference: contents: - concordance_survival +- title: Linear Predictor Survival Metrics + contents: + - royston_survival + - title: Curve Survival Functions contents: - roc_curve_survival From 1f9f64db03f97ae75964c86d7408f021379de5f8 Mon Sep 17 00:00:00 2001 From: Hannah Frick Date: Wed, 17 Dec 2025 14:08:21 +0000 Subject: [PATCH 6/9] Update R/surv-royston.R Co-authored-by: Emil Hvitfeldt --- R/surv-royston.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/surv-royston.R b/R/surv-royston.R index 651327e8..50b8b914 100644 --- a/R/surv-royston.R +++ b/R/surv-royston.R @@ -115,7 +115,7 @@ normal_score_blom <- function(x, case_weights) { # does not need kappa (from Royston & Sauerbrei 2004) because it'll # "disappear" into the baseline hazard of the Cox model in # `royston_survival_impl()` - z = stats::qnorm((.data$x_first - 3 / 8) / (dplyr::n() + 0.25)) + z = stats::qnorm((.data$x_first - 3 / 8) / (dplyr::n() + 1 / 4)) ) |> # average over ties dplyr::mutate(s = mean(.data$z), .by = "x") |> From 6b2520762a03bfb844dd13c2c74658754db6524a Mon Sep 17 00:00:00 2001 From: Hannah Frick Date: Wed, 17 Dec 2025 16:12:01 +0000 Subject: [PATCH 7/9] more details on scaling --- R/surv-royston.R | 3 +++ 1 file changed, 3 insertions(+) diff --git a/R/surv-royston.R b/R/surv-royston.R index 50b8b914..80de7782 100644 --- a/R/surv-royston.R +++ b/R/surv-royston.R @@ -102,6 +102,9 @@ royston_survival_impl <- function(truth, estimate, case_weights) { fit <- survival::coxph(truth ~ bns, weights = case_weights) est <- unname(stats::coef(fit)) + # the regression coefficient is sigma* in Royston & Sauerbrei 2004 + # because `normal_score_blom()` returns their u_i (not z_i) + # we thus don't need to scale the coefficient with kappa here est^2 / (est^2 + pi^2 / 6) } From 8f32b67d0775b0b759cbf5fce6547cad5f821aaa Mon Sep 17 00:00:00 2001 From: Hannah Frick Date: Wed, 17 Dec 2025 16:17:13 +0000 Subject: [PATCH 8/9] the test is a little flaky --- tests/testthat/test-surv-royston.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/testthat/test-surv-royston.R b/tests/testthat/test-surv-royston.R index dc68bab8..72b97cec 100644 --- a/tests/testthat/test-surv-royston.R +++ b/tests/testthat/test-surv-royston.R @@ -19,7 +19,7 @@ test_that("comparison test with survival::royston", { res[[".estimate"]], royston_ref["R.D"], ignore_attr = TRUE, - tolerance = 1e-4 + tolerance = 1e-3 ) }) From 77c1d6354505e2497575107b7a5461ee8ed78e6c Mon Sep 17 00:00:00 2001 From: Hannah Frick Date: Fri, 19 Dec 2025 15:38:40 +0000 Subject: [PATCH 9/9] allow case weights of zero --- R/surv-royston.R | 9 +++++++-- tests/testthat/test-surv-royston.R | 9 ++++++++- 2 files changed, 15 insertions(+), 3 deletions(-) diff --git a/R/surv-royston.R b/R/surv-royston.R index 80de7782..0499866c 100644 --- a/R/surv-royston.R +++ b/R/surv-royston.R @@ -109,7 +109,10 @@ royston_survival_impl <- function(truth, estimate, case_weights) { } normal_score_blom <- function(x, case_weights) { - tibble::tibble( + # includes observations with zero weights + x_0 <- tibble::tibble(.row = seq_along(x), x = x) + + rankits <- tibble::tibble( .row = rep(seq_along(x), times = case_weights), x = rep(x, times = case_weights), ) |> @@ -122,6 +125,8 @@ normal_score_blom <- function(x, case_weights) { ) |> # average over ties dplyr::mutate(s = mean(.data$z), .by = "x") |> - dplyr::slice(1, .by = .row) |> + dplyr::slice(1, .by = .row) + + dplyr::left_join(x_0, rankits, by = c(".row")) |> dplyr::pull("s") } diff --git a/tests/testthat/test-surv-royston.R b/tests/testthat/test-surv-royston.R index 72b97cec..d460e37e 100644 --- a/tests/testthat/test-surv-royston.R +++ b/tests/testthat/test-surv-royston.R @@ -43,6 +43,13 @@ test_that("`normal_score_blom()` works with case weights", { nsb[1], mean(qnorm((1:4 - 3 / 8) / (sum(case_weights) + 0.25))) ) + + # weights of zero + x <- 1:10 + 0.5 + case_weights <- c(0, 0, rep(2, 8)) + nsb <- normal_score_blom(x, case_weights) + expect_length(nsb, 10) + expect_true(all(is.na(nsb[1:2]))) }) test_that("case weights works with equal weights", { @@ -70,7 +77,7 @@ test_that("case weights works with equal weights", { test_that("works with hardhat case weights", { lung_surv <- data_lung_surv() - lung_surv$case_wts <- rep(2, nrow(lung_surv)) + lung_surv$case_wts <- c(rep(0, 10), rep(2, nrow(lung_surv) - 10)) df <- lung_surv