library(tidyverse)
library(tableone) #create 'Table 1' to describe baseline characteristics
library(Matching) #multivariate and propensity score matching with balance optimization

Dataset : right heart catheterization

Une description complète est disponible à cette adresse http://biostat.mc.vanderbilt.edu/wiki/pub/Main/DataSets/rhc.html. Il s’agit de données sur des patients (2184 traités et 3551 contrôles) admis aux urgences dans 5 hôpitaux, la variable de traitement est swang1 (right heat catheterization vs. non) et l’outcome est death (yes or no). Nous allons considérer les variables de confusion suivantes : - cat1 : Primary disease category - age - sex`` -meanbp1`: Mean blood pressure

Matching (greedy)

Questions :

rhc = read_csv("./Data/rhc_data.csv")
rhc_small = rhc %>% mutate(treatment = swang1 ) %>% 
                dplyr::select(treatment,death, cat1, age,sex,meanbp1,-swang1)
                
glimpse(rhc_small)
unique(rhc_small$cat1)
  • Créer la “Table 1”. Attention le package tableone ne supporte pas les variables discrètes à plusieurs facteurs, il faut donc d’abord construire la table disjonctive.

Correction

rhc_small_disjunctive = as_tibble(model.matrix(~ . , data=rhc_small)[,-1])
glimpse(rhc_small_disjunctive)
names(rhc_small_disjunctive) = c("treatment","Death","CHF","Cirr",
                                 "colcan","Coma","COPD",
                                 "lungcan","MOSF","sepsis",
                                 "age","male","meanbp1")

Correction

xvars = c("CHF","Cirr","colcan","Coma","COPD",
         "lungcan","MOSF","sepsis",
         "age","male","meanbp1") # variables de confusion
table1 = CreateTableOne(vars=xvars,strata="treatment", data=rhc_small_disjunctive, test=FALSE)
print(table1,smd=TRUE)
  • Faire un matching “greedy” basé sur la distance de Mahalanobis.

Correction

greedymatch = Match(Tr=rhc_small_disjunctive$treatment,
                    M=1,
                    X=rhc_small_disjunctive[xvars],replace=FALSE)

rhc_small_disjunctive_matched = rhc_small_disjunctive %>%
                                  slice(c(greedymatch$index.control,
                                          greedymatch$index.treated)) %>%
                                  mutate( match_id = rep(c(1:2184),2))
                                

matchedtab1 = CreateTableOne(vars=xvars,strata="treatment", data=rhc_small_disjunctive_matched, test=FALSE)
print(matchedtab1, smd = TRUE)
  • Sur les données “matchées” étudiez l’outcome

Correction

y_trt = rhc_small_disjunctive_matched %>% filter(treatment==1) %>% 
                                      dplyr::select(Death) %>% pull()
y_con = rhc_small_disjunctive_matched %>% filter(treatment==0) %>%
                                      dplyr::select(Death)%>% pull()
table(y_trt,y_con)
mcnemar.test(matrix(table(y_trt,y_con),2,2))

Propensity score matching

Correction

psmodel = glm(treatment ~ . - Death ,
              family=binomial(),data=rhc_small_disjunctive)

summary(psmodel)
#create propensity scores
pscore = predict(psmodel,type="response")
  • Représenter les distributions des scores de propension chez les contrôles et les traités.

Correction

ggplot(rhc_small_disjunctive %>% mutate(propensity = pscore) ,
       aes(as.factor(treatment),pscore,fill = treatment) )+ geom_violin()
  • Calculer le logit des scores de propension, puis leur écart-type et le “caliper”

Correction

logit = function(p) {log(p)-log(1-p)}
logit_pscore = logit(pscore) # or predict(psmodel)
caliper = 0.2 *sd(logit_pscore)
  • Appliquer le greedy matching sur le logitdu score de propension avec le “capiler” calculé à la question précédente

Correction

psmatch = Match(Tr=rhc_small_disjunctive$treatment,M=1,X=logit_pscore,
             replace=FALSE,caliper=caliper)
n_matched = length(psmatch$index.control)

rhc_small_disjunctive_matched_propensity =rhc_small_disjunctive %>%
                                  slice(c(psmatch$index.control,
                                          psmatch$index.treated)) %>%
                                  mutate( match_id = rep(c(1:n_matched),2))
  • Calculer la “Table 1”

Correction

matchedtab1<-CreateTableOne(vars=xvars, strata ="treatment", 
                            data=rhc_small_disjunctive_matched_propensity, test = FALSE)
print(matchedtab1, smd = TRUE)
  • Sur les données “matchées” étudiez l’outcome

Correction

y_trt = rhc_small_disjunctive_matched_propensity %>% filter(treatment==1) %>% 
                                      dplyr::select(Death) %>% pull()
y_con = rhc_small_disjunctive_matched_propensity %>% filter(treatment==0) %>%
                                      dplyr::select(Death)%>% pull()
table(y_trt,y_con)
mcnemar.test(matrix(table(y_trt,y_con),2,2))
  • Faire varier la valeur du “caliper” et observer les changements dans la “Table 1” et le test sur l’outcome

Correction

for (cal_const in c(0.0001,0.01,0.1,0.5,1,5)){
  
  caliper = cal_const *sd(logit_pscore)
  print(paste(c("Le caliper vaut ",caliper),collapse=""))
  psmatch = Match(Tr=rhc_small_disjunctive$treatment,M=1,X=logit_pscore,
               replace=FALSE,caliper=caliper)
  n_matched = length(psmatch$index.control)
  
  rhc_small_disjunctive_matched_propensity =rhc_small_disjunctive %>%
                                    slice(c(psmatch$index.control,
                                            psmatch$index.treated)) %>%
                                    mutate( match_id = rep(c(1:n_matched),2))
  print(paste(c("Il y a ",length(psmatch$index.treated)," données matchées"),collapse=""))
  matchedtab1<-CreateTableOne(vars=xvars, strata ="treatment", 
                              data=rhc_small_disjunctive_matched_propensity, test = FALSE)
  print(matchedtab1, smd = TRUE)
  y_trt = rhc_small_disjunctive_matched_propensity %>% filter(treatment==1) %>% 
                                        dplyr::select(Death) %>% pull()
  y_con = rhc_small_disjunctive_matched_propensity %>% filter(treatment==0) %>%
                                        dplyr::select(Death)%>% pull()
  print(table(y_trt,y_con))
  print(mcnemar.test(matrix(table(y_trt,y_con),2,2)))
  }
LS0tCnRpdGxlOiAiRXhlbXBsZXMgcG91ciBsYSBjYXVzYWxpdMOpIDogYXBwYXJpZW1lbnQgZXQgc2NvcmUgZGUgcHJvcGVuc2lvbiAoY2hhcGl0cmUgMikiCm91dHB1dDogaHRtbF9ub3RlYm9vawotLS0KCmBgYHtyfQpsaWJyYXJ5KHRpZHl2ZXJzZSkKbGlicmFyeSh0YWJsZW9uZSkgI2NyZWF0ZSAnVGFibGUgMScgdG8gZGVzY3JpYmUgYmFzZWxpbmUgY2hhcmFjdGVyaXN0aWNzCmxpYnJhcnkoTWF0Y2hpbmcpICNtdWx0aXZhcmlhdGUgYW5kIHByb3BlbnNpdHkgc2NvcmUgbWF0Y2hpbmcgd2l0aCBiYWxhbmNlIG9wdGltaXphdGlvbgpgYGAKCiMgRGF0YXNldCA6IHJpZ2h0IGhlYXJ0IGNhdGhldGVyaXphdGlvbgpVbmUgZGVzY3JpcHRpb24gY29tcGzDqHRlIGVzdCBkaXNwb25pYmxlIMOgIGNldHRlIGFkcmVzc2UgaHR0cDovL2Jpb3N0YXQubWMudmFuZGVyYmlsdC5lZHUvd2lraS9wdWIvTWFpbi9EYXRhU2V0cy9yaGMuaHRtbC4gSWwgcydhZ2l0IGRlIGRvbm7DqWVzIHN1ciBkZXMgcGF0aWVudHMgKDIxODQgdHJhaXTDqXMgZXQgMzU1MSBjb250csO0bGVzKSBhZG1pcyBhdXggdXJnZW5jZXMgZGFucyA1IGjDtHBpdGF1eCwgbGEgdmFyaWFibGUgZGUgdHJhaXRlbWVudCBlc3QgYHN3YW5nMWAgKHJpZ2h0IGhlYXQgY2F0aGV0ZXJpemF0aW9uIHZzLiBub24pIGV0IGwnb3V0Y29tZSBlc3QgYGRlYXRoYCAoeWVzIG9yIG5vKS4gTm91cyBhbGxvbnMgY29uc2lkw6lyZXIgbGVzIHZhcmlhYmxlcyBkZSBjb25mdXNpb24gc3VpdmFudGVzIDoKLSBgY2F0MWAgOiBQcmltYXJ5IGRpc2Vhc2UgY2F0ZWdvcnkKLSBgYWdlYAotIGBzZXhgYAotIGBtZWFuYnAxYDogTWVhbiBibG9vZCBwcmVzc3VyZQoKIyBNYXRjaGluZyAoZ3JlZWR5KQojIyBRdWVzdGlvbnMgOiAKLSBDaGFyZ2VyIGxlcyBkb25uw6llcyBodHRwOi8vYmlvc3RhdC5tYy52YW5kZXJiaWx0LmVkdS93aWtpL01haW4vRGF0YVNldHMsIApgYGB7cn0KcmhjID0gcmVhZF9jc3YoIi4vRGF0YS9yaGNfZGF0YS5jc3YiKQpyaGNfc21hbGwgPSByaGMgJT4lIG11dGF0ZSh0cmVhdG1lbnQgPSBzd2FuZzEgKSAlPiUgCiAgICAgICAgICAgICAgICBkcGx5cjo6c2VsZWN0KHRyZWF0bWVudCxkZWF0aCwgY2F0MSwgYWdlLHNleCxtZWFuYnAxLC1zd2FuZzEpCiAgICAgICAgICAgICAgICAKZ2xpbXBzZShyaGNfc21hbGwpCnVuaXF1ZShyaGNfc21hbGwkY2F0MSkKYGBgCi0gQ3LDqWVyIGxhICJUYWJsZSAxIi4gQXR0ZW50aW9uIGxlIHBhY2thZ2UgYHRhYmxlb25lYCBuZSBzdXBwb3J0ZSBwYXMgbGVzIHZhcmlhYmxlcyBkaXNjcsOodGVzIMOgIHBsdXNpZXVycyBmYWN0ZXVycywgaWwgZmF1dCBkb25jIGQnYWJvcmQgY29uc3RydWlyZSBsYSB0YWJsZSBkaXNqb25jdGl2ZS4KYGBge3J9CgpgYGAKIyMjIyBDb3JyZWN0aW9uCmBgYHtyfQpyaGNfc21hbGxfZGlzanVuY3RpdmUgPSBhc190aWJibGUobW9kZWwubWF0cml4KH4gLiAsIGRhdGE9cmhjX3NtYWxsKVssLTFdKQpnbGltcHNlKHJoY19zbWFsbF9kaXNqdW5jdGl2ZSkKbmFtZXMocmhjX3NtYWxsX2Rpc2p1bmN0aXZlKSA9IGMoInRyZWF0bWVudCIsIkRlYXRoIiwiQ0hGIiwiQ2lyciIsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICJjb2xjYW4iLCJDb21hIiwiQ09QRCIsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICJsdW5nY2FuIiwiTU9TRiIsInNlcHNpcyIsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICJhZ2UiLCJtYWxlIiwibWVhbmJwMSIpCmBgYAoKYGBge3J9CgpgYGAKIyMjIyBDb3JyZWN0aW9uCmBgYHtyfQp4dmFycyA9IGMoIkNIRiIsIkNpcnIiLCJjb2xjYW4iLCJDb21hIiwiQ09QRCIsCiAgICAgICAgICJsdW5nY2FuIiwiTU9TRiIsInNlcHNpcyIsCiAgICAgICAgICJhZ2UiLCJtYWxlIiwibWVhbmJwMSIpICMgdmFyaWFibGVzIGRlIGNvbmZ1c2lvbgp0YWJsZTEgPSBDcmVhdGVUYWJsZU9uZSh2YXJzPXh2YXJzLHN0cmF0YT0idHJlYXRtZW50IiwgZGF0YT1yaGNfc21hbGxfZGlzanVuY3RpdmUsIHRlc3Q9RkFMU0UpCnByaW50KHRhYmxlMSxzbWQ9VFJVRSkKYGBgCgotIEZhaXJlIHVuIG1hdGNoaW5nICJncmVlZHkiIGJhc8OpIHN1ciBsYSBkaXN0YW5jZSBkZSBNYWhhbGFub2Jpcy4KYGBge3J9CgpgYGAKIyMjIyBDb3JyZWN0aW9uCmBgYHtyfQpncmVlZHltYXRjaCA9IE1hdGNoKFRyPXJoY19zbWFsbF9kaXNqdW5jdGl2ZSR0cmVhdG1lbnQsCiAgICAgICAgICAgICAgICAgICAgTT0xLAogICAgICAgICAgICAgICAgICAgIFg9cmhjX3NtYWxsX2Rpc2p1bmN0aXZlW3h2YXJzXSxyZXBsYWNlPUZBTFNFKQoKcmhjX3NtYWxsX2Rpc2p1bmN0aXZlX21hdGNoZWQgPSByaGNfc21hbGxfZGlzanVuY3RpdmUgJT4lCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBzbGljZShjKGdyZWVkeW1hdGNoJGluZGV4LmNvbnRyb2wsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIGdyZWVkeW1hdGNoJGluZGV4LnRyZWF0ZWQpKSAlPiUKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIG11dGF0ZSggbWF0Y2hfaWQgPSByZXAoYygxOjIxODQpLDIpKQogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIAoKbWF0Y2hlZHRhYjEgPSBDcmVhdGVUYWJsZU9uZSh2YXJzPXh2YXJzLHN0cmF0YT0idHJlYXRtZW50IiwgZGF0YT1yaGNfc21hbGxfZGlzanVuY3RpdmVfbWF0Y2hlZCwgdGVzdD1GQUxTRSkKcHJpbnQobWF0Y2hlZHRhYjEsIHNtZCA9IFRSVUUpCmBgYAotIFN1ciBsZXMgZG9ubsOpZXMgIm1hdGNow6llcyIgw6l0dWRpZXogbCdvdXRjb21lCmBgYHtyfQoKYGBgCiMjIyMgQ29ycmVjdGlvbgpgYGB7cn0KeV90cnQgPSByaGNfc21hbGxfZGlzanVuY3RpdmVfbWF0Y2hlZCAlPiUgZmlsdGVyKHRyZWF0bWVudD09MSkgJT4lIAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIGRwbHlyOjpzZWxlY3QoRGVhdGgpICU+JSBwdWxsKCkKeV9jb24gPSByaGNfc21hbGxfZGlzanVuY3RpdmVfbWF0Y2hlZCAlPiUgZmlsdGVyKHRyZWF0bWVudD09MCkgJT4lCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgZHBseXI6OnNlbGVjdChEZWF0aCklPiUgcHVsbCgpCnRhYmxlKHlfdHJ0LHlfY29uKQptY25lbWFyLnRlc3QobWF0cml4KHRhYmxlKHlfdHJ0LHlfY29uKSwyLDIpKQpgYGAKIyBQcm9wZW5zaXR5IHNjb3JlIG1hdGNoaW5nCgotIEVzdGltZXIgbGUgc2NvcmUgcHJvcGVuc2lvbiBncsOiY2Ugw6AgdW5lIHLDqWdyZXNzaW9uIGxvZ2lzdGlxdWUgZXQgY2FsY3VsZXIgbGVzIHNjb3JlcyBkZSBwcm9wZW5zaW9uLgpgYGB7cn0KCmBgYAojIyMjIENvcnJlY3Rpb24KYGBge3J9CnBzbW9kZWwgPSBnbG0odHJlYXRtZW50IH4gLiAtIERlYXRoICwKICAgICAgICAgICAgICBmYW1pbHk9Ymlub21pYWwoKSxkYXRhPXJoY19zbWFsbF9kaXNqdW5jdGl2ZSkKCnN1bW1hcnkocHNtb2RlbCkKI2NyZWF0ZSBwcm9wZW5zaXR5IHNjb3Jlcwpwc2NvcmUgPSBwcmVkaWN0KHBzbW9kZWwsdHlwZT0icmVzcG9uc2UiKQpgYGAKLSBSZXByw6lzZW50ZXIgbGVzIGRpc3RyaWJ1dGlvbnMgZGVzIHNjb3JlcyBkZSBwcm9wZW5zaW9uIGNoZXogbGVzIGNvbnRyw7RsZXMgZXQgbGVzIHRyYWl0w6lzLgpgYGB7cn0KCmBgYAojIyMjIENvcnJlY3Rpb24KYGBge3J9CmdncGxvdChyaGNfc21hbGxfZGlzanVuY3RpdmUgJT4lIG11dGF0ZShwcm9wZW5zaXR5ID0gcHNjb3JlKSAsCiAgICAgICBhZXMoYXMuZmFjdG9yKHRyZWF0bWVudCkscHNjb3JlLGZpbGwgPSB0cmVhdG1lbnQpICkrIGdlb21fdmlvbGluKCkKYGBgCi0gQ2FsY3VsZXIgbGUgbG9naXQgZGVzIHNjb3JlcyBkZSBwcm9wZW5zaW9uLCBwdWlzIGxldXIgw6ljYXJ0LXR5cGUgZXQgbGUgImNhbGlwZXIiCmBgYHtyfQoKYGBgCiMjIyMgQ29ycmVjdGlvbgpgYGB7cn0KbG9naXQgPSBmdW5jdGlvbihwKSB7bG9nKHApLWxvZygxLXApfQpsb2dpdF9wc2NvcmUgPSBsb2dpdChwc2NvcmUpICMgb3IgcHJlZGljdChwc21vZGVsKQpjYWxpcGVyID0gMC4yICpzZChsb2dpdF9wc2NvcmUpCmBgYAoKLSBBcHBsaXF1ZXIgbGUgZ3JlZWR5IG1hdGNoaW5nIHN1ciBsZSBgbG9naXRgZHUgc2NvcmUgZGUgcHJvcGVuc2lvbiBhdmVjIGxlICJjYXBpbGVyIiBjYWxjdWzDqSDDoCBsYSBxdWVzdGlvbiBwcsOpY8OpZGVudGUKYGBge3J9CgpgYGAKIyMjIyBDb3JyZWN0aW9uCmBgYHtyfQpwc21hdGNoID0gTWF0Y2goVHI9cmhjX3NtYWxsX2Rpc2p1bmN0aXZlJHRyZWF0bWVudCxNPTEsWD1sb2dpdF9wc2NvcmUsCiAgICAgICAgICAgICByZXBsYWNlPUZBTFNFLGNhbGlwZXI9Y2FsaXBlcikKbl9tYXRjaGVkID0gbGVuZ3RoKHBzbWF0Y2gkaW5kZXguY29udHJvbCkKCnJoY19zbWFsbF9kaXNqdW5jdGl2ZV9tYXRjaGVkX3Byb3BlbnNpdHkgPXJoY19zbWFsbF9kaXNqdW5jdGl2ZSAlPiUKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIHNsaWNlKGMocHNtYXRjaCRpbmRleC5jb250cm9sLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBwc21hdGNoJGluZGV4LnRyZWF0ZWQpKSAlPiUKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIG11dGF0ZSggbWF0Y2hfaWQgPSByZXAoYygxOm5fbWF0Y2hlZCksMikpCgoKYGBgCi0gQ2FsY3VsZXIgbGEgIlRhYmxlIDEiCmBgYHtyfQoKYGBgCiMjIyMgQ29ycmVjdGlvbgpgYGB7cn0KbWF0Y2hlZHRhYjE8LUNyZWF0ZVRhYmxlT25lKHZhcnM9eHZhcnMsIHN0cmF0YSA9InRyZWF0bWVudCIsIAogICAgICAgICAgICAgICAgICAgICAgICAgICAgZGF0YT1yaGNfc21hbGxfZGlzanVuY3RpdmVfbWF0Y2hlZF9wcm9wZW5zaXR5LCB0ZXN0ID0gRkFMU0UpCnByaW50KG1hdGNoZWR0YWIxLCBzbWQgPSBUUlVFKQpgYGAKCi0gU3VyIGxlcyBkb25uw6llcyAibWF0Y2jDqWVzIiDDqXR1ZGlleiBsJ291dGNvbWUKYGBge3J9CgpgYGAKIyMjIyBDb3JyZWN0aW9uCmBgYHtyfQp5X3RydCA9IHJoY19zbWFsbF9kaXNqdW5jdGl2ZV9tYXRjaGVkX3Byb3BlbnNpdHkgJT4lIGZpbHRlcih0cmVhdG1lbnQ9PTEpICU+JSAKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBkcGx5cjo6c2VsZWN0KERlYXRoKSAlPiUgcHVsbCgpCnlfY29uID0gcmhjX3NtYWxsX2Rpc2p1bmN0aXZlX21hdGNoZWRfcHJvcGVuc2l0eSAlPiUgZmlsdGVyKHRyZWF0bWVudD09MCkgJT4lCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgZHBseXI6OnNlbGVjdChEZWF0aCklPiUgcHVsbCgpCnRhYmxlKHlfdHJ0LHlfY29uKQptY25lbWFyLnRlc3QobWF0cml4KHRhYmxlKHlfdHJ0LHlfY29uKSwyLDIpKQpgYGAKLSBGYWlyZSB2YXJpZXIgbGEgdmFsZXVyIGR1ICJjYWxpcGVyIiBldCBvYnNlcnZlciBsZXMgY2hhbmdlbWVudHMgZGFucyBsYSAiVGFibGUgMSIgZXQgbGUgdGVzdCBzdXIgbCdvdXRjb21lCmBgYHtyfQoKYGBgCiMjIyMgQ29ycmVjdGlvbgpgYGB7cn0KZm9yIChjYWxfY29uc3QgaW4gYygwLjAwMDEsMC4wMSwwLjEsMC41LDEsNSkpewogIAogIGNhbGlwZXIgPSBjYWxfY29uc3QgKnNkKGxvZ2l0X3BzY29yZSkKICBwcmludChwYXN0ZShjKCJMZSBjYWxpcGVyIHZhdXQgIixjYWxpcGVyKSxjb2xsYXBzZT0iIikpCiAgcHNtYXRjaCA9IE1hdGNoKFRyPXJoY19zbWFsbF9kaXNqdW5jdGl2ZSR0cmVhdG1lbnQsTT0xLFg9bG9naXRfcHNjb3JlLAogICAgICAgICAgICAgICByZXBsYWNlPUZBTFNFLGNhbGlwZXI9Y2FsaXBlcikKICBuX21hdGNoZWQgPSBsZW5ndGgocHNtYXRjaCRpbmRleC5jb250cm9sKQogIAogIHJoY19zbWFsbF9kaXNqdW5jdGl2ZV9tYXRjaGVkX3Byb3BlbnNpdHkgPXJoY19zbWFsbF9kaXNqdW5jdGl2ZSAlPiUKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgc2xpY2UoYyhwc21hdGNoJGluZGV4LmNvbnRyb2wsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgcHNtYXRjaCRpbmRleC50cmVhdGVkKSkgJT4lCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIG11dGF0ZSggbWF0Y2hfaWQgPSByZXAoYygxOm5fbWF0Y2hlZCksMikpCiAgcHJpbnQocGFzdGUoYygiSWwgeSBhICIsbGVuZ3RoKHBzbWF0Y2gkaW5kZXgudHJlYXRlZCksIiBkb25uw6llcyBtYXRjaMOpZXMiKSxjb2xsYXBzZT0iIikpCiAgbWF0Y2hlZHRhYjE8LUNyZWF0ZVRhYmxlT25lKHZhcnM9eHZhcnMsIHN0cmF0YSA9InRyZWF0bWVudCIsIAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICBkYXRhPXJoY19zbWFsbF9kaXNqdW5jdGl2ZV9tYXRjaGVkX3Byb3BlbnNpdHksIHRlc3QgPSBGQUxTRSkKICBwcmludChtYXRjaGVkdGFiMSwgc21kID0gVFJVRSkKICB5X3RydCA9IHJoY19zbWFsbF9kaXNqdW5jdGl2ZV9tYXRjaGVkX3Byb3BlbnNpdHkgJT4lIGZpbHRlcih0cmVhdG1lbnQ9PTEpICU+JSAKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIGRwbHlyOjpzZWxlY3QoRGVhdGgpICU+JSBwdWxsKCkKICB5X2NvbiA9IHJoY19zbWFsbF9kaXNqdW5jdGl2ZV9tYXRjaGVkX3Byb3BlbnNpdHkgJT4lIGZpbHRlcih0cmVhdG1lbnQ9PTApICU+JQogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgZHBseXI6OnNlbGVjdChEZWF0aCklPiUgcHVsbCgpCiAgcHJpbnQodGFibGUoeV90cnQseV9jb24pKQogIHByaW50KG1jbmVtYXIudGVzdChtYXRyaXgodGFibGUoeV90cnQseV9jb24pLDIsMikpKQogIH0KYGBgCgo=