Black-market survival analysis: 2014

library(XML)
black <- readHTMLTable("http://www.gwern.net/Black-market%20survival", colClasses = c("character", "factor", "factor", "logical",
                                                               "factor", "logical", "character", "logical", "logical",
                                                               "factor", "factor", "logical", "logical",
                                                               "integer", "logical",
                                                               "character"))[[1]]
black[black$Codebase=="",]$Codebase <- NA
black$Start <- as.Date(black$Start); black$End <- as.Date(black$End)
black[is.na(black$End),]$End <- Sys.Date() # assuming the table is up to date
black$Age <- black$End - black$Start
black$Age <- as.integer(black$End - black$Start)

# for competing-risks modeling
levels(black$Closure) <- c("alive", "hacked", "raided", "scam", "voluntary") # sets "" to "alive"

sum(sort(as.Date(black$Start)) >= as.Date("2014-01-01"))
# [1] 43
sort(black[as.Date(black$Start) >= as.Date("2014-01-01"),]$Market)
#  [1] 1776                     Abraxas                  Alpaca Marketplace
#  [4] AlphaBay                 Andromeda                Area51
#  [7] Armory Vendor Market     BlackBank Market         Black Goblin Market
# [10] Breaking Bad             Cannabis Road            Cannabis Road 2
# [13] Cannabis Road 3          Cantina                  Cloud-Nine
# [16] DarkBay                  Darknet Nation           Deepzon
# [19] Diabolus                 Doge Road                drugslist
# [22] evolution                EXXTACY                  Freedom Market
# [25] GreyRoad                 Hansa                    Hydra
# [28] Middle Earth Marketplace Mr Nice Guy              Nucleus Marketplace
# [31] Onionshop                Panacea                  Pigeon Market
# [34] Red Sun Marketplace      Sanitarium Market        Silk Street
# [37] TOM                      Topix 2                  TorBazaar
# [40] TorEscrow                Tortuga 2                Underground Market
# [43] Utopia
length(black[black$End > as.Date("2014-01-01") & black$End < as.Date("2015-01-01"),]$Market)
# [1] 46
sort(black[black$End > as.Date("2014-01-01") & black$End < as.Date("2015-01-01"),]$Market)
#  [1] 1776                  Alpaca Marketplace    Andromeda             Armory Vendor Market
#  [5] BlackBox Market       Black Goblin Market   Black Services Market Blue Sky
#  [9] Breaking Bad          BuyItNow              Cannabis Road         Cannabis Road 2
# [13] Cannabis Road 3       Cantina               Cloud-Nine            DarkBay
# [17] Darknet Nation        Deepzon               Doge Road             drugslist
# [21] EXXTACY               FreeBay               Freedom Market        GreyRoad
# [25] Hansa                 Hydra                 Mr Nice Guy           Onionshop
# [29] Pandora               Pigeon Market         Pirate Market         Red Sun Marketplace
# [33] Sanitarium Market     Silk Road 2           Silk Street           The Marketplace
# [37] TOM                   Topix 2               TorBay                TorBazaar
# [41] TorEscrow             Tortuga               Tortuga 2             Underground Market
# [45] Utopia                White Rabbit

sort(black[black$Alive,]$Market)
#  [1] Abraxas                  Agora                    AlphaBay
#  [4] Area51                   BlackBank Market         Diabolus
#  [7] Dream Market             evolution                Middle Earth Marketplace
# [10] Nucleus Marketplace      Outlaw Market            Panacea
# 67 Levels: 1776 Abraxas Agora Alpaca Marketplace AlphaBay Andromeda Area51 ... White Rabbit
sort(black[black$Alive & black$Start >= as.Date("2014-01-01"),]$Market)
# [1] Abraxas                  AlphaBay                 Area51
# [4] BlackBank Market         Diabolus                 evolution
# [7] Middle Earth Marketplace Nucleus Marketplace      Panacea
# 67 Levels: 1776 Abraxas Agora Alpaca Marketplace AlphaBay Andromeda Area51 ... White Rabbit

sort(table(black[as.Date(black$End) >= as.Date("2014-01-01"),]$Closure))
#    raided    hacked     alive      scam voluntary
#         6         9        12        12        20

library(survival)
summary(survfit(Surv(Age, Alive, type="right") ~ 1, data=black))
#  time n.risk n.event survival  std.err lower 95% CI upper 95% CI
#    11     58       1  0.98276 0.017092     0.949823      1.00000
#    20     54       1  0.96456 0.024627     0.917479      1.00000
#    67     38       1  0.93918 0.034675     0.873615      1.00000
#    70     37       1  0.91379 0.042013     0.835050      0.99996
#    81     35       1  0.88768 0.048248     0.797983      0.98747
#   194     23       1  0.84909 0.059621     0.739919      0.97437
#   196     22       1  0.81049 0.068269     0.687151      0.95598
#   331     10       1  0.72945 0.098424     0.559938      0.95027
#   353      7       1  0.62524 0.128160     0.418378      0.93438
#   369      5       1  0.50019 0.151728     0.276014      0.90644
#   395      4       1  0.37514 0.157090     0.165105      0.85238
#   413      3       1  0.25010 0.146261     0.079489      0.78688
plot(survfit(Surv(Age, Alive, type="right") ~ 1, data=black))
# https://i.imgur.com/OamVzjO.png

library(cmprsk)
blackc <- with(black, cuminc(Age, Closure, cencode="alive")); blackc
# Estimates and Variances:
# $est
#                      200         400         600         800
# 1 hacked    0.1389819240 0.138981924 0.138981924 0.138981924
# 1 raided    0.0149253731 0.114655933 0.114655933 0.114655933
# 1 scam      0.2069723259 0.299188940 0.299188940 0.299188940
# 1 voluntary 0.2702539221 0.327633148 0.327633148 0.327633148
#
# $var
#                        200           400           600           800
# 1 hacked    0.001900471515 0.00190047151 0.00190047151 0.00190047151
# 1 raided    0.000223039082 0.00212491076 0.00212491076 0.00212491076
# 1 scam      0.002707980885 0.00366455525 0.00366455525 0.00366455525
# 1 voluntary 0.003240020448 0.00383384036 0.00383384036 0.00383384036

# risk of shut down by cause:
plot(blackc, lty=1, color=c(1:4), xlab="Days")
# https://i.imgur.com/xZE8O5n.png

# Generate predictions of shutdown in next 6 / 12 months for the still-running marketplaces:

# You are not expected to understand this ugly hack. I don't either. But it *seems* to work...
# for background, see http://gwern.net/Google%20shutdowns & http://www.gwern.net/Silk%20Road#the-bet-bmr-or-sheep-to-die-in-a-year-by-oct-2014
library(rms)
conditionalProbability <- function (d, followupUnits, cmodel) {
    chances <- rep(NA, nrow(d)) # stash results

    for (i in 1:nrow(d)) {

        # extract chance of particular subject surviving as long as it has:
        beginProb <- survest(cmodel, d[i,], times=(d[i,]$Age))$surv
        if (length(beginProb)==0) { beginProb <- 1 } # set to a default

        tmpFollowup <- followupUnits # reset in each for loop
        while (TRUE) {
            # extract chance of subject surviving as long as it has + an arbitrary additional time-units
            endProb <- survest(cmodel, d[i,], times=(d[i,]$Age + tmpFollowup))$surv
            # survival curve may not reach that far! 'survexp returns 'numeric(0)' if it doesn't;
            # so we shrink down 1 day and try again until 'survexp' *does* return a usable answer
            if (length(endProb)==0) { tmpFollowup <- tmpFollowup - 1} else { break }
        }

        # if 50% of all subjects survive to time t, and 20% of all survive to time t+100, say, what chance
        # does a survivor - at exactly time t - have of making it to time t+100? 40%: 0.20 / 0.50 = 0.40
        chances[i] <- endProb / beginProb
    }
    return(chances)
}
# we can't regress on '~ 1' because `survest` crashes inside `conditionalProbability` (‽),
# so let's just toss in all the covariates to make it work despite none of them being well-justified...
cpmodel <- cph(Surv(Age, Alive) ~ Multisig + I2P + Coin + Guns + Fraud, data = black, x=TRUE, y=TRUE, surv=TRUE)
# report in %s rather than raw probabilities for easier reading
black$SixMonth <- round(digits=1, conditionalProbability(black, round(365/2), cpmodel) * 100)
black$OneYear  <- round(digits=1, conditionalProbability(black,          365, cpmodel) * 100)

subset(black, Alive, select=c("Market", "SixMonth", "OneYear"))
#                      Market SixMonth OneYear
# 56             Dream Market    100.0   100.0
# 57                    Agora     65.6    65.6
# 58            Outlaw Market     72.3    72.3
# 59                evolution     17.5    17.5
# 60         BlackBank Market     13.1    13.1
# 61                   Area51     57.7    28.1
# 62 Middle Earth Marketplace     56.8    28.8
# 63                 Diabolus     34.7     0.0
# 64      Nucleus Marketplace     15.7     0.0
# 65                  Panacea     89.7    34.4
# 66                  Abraxas     85.9    51.2
# 67                 AlphaBay     92.5    48.4