test/rescreen.R

In mclements/prostata: Prostate Cancer Microsimulation for Sweden

## Shuang's parameters
## in shell, "ssh -X marcle@vector" and (ess-remote) on that buffer

## A: Discount funtion
## B: Events summary
## BASE CASE (S2): Simulation by HSV set PORPUS-U
## Sensitivity (S0): Simulation by HSV set EQ-5D
## Sensitivity (S1): Simulation by HSV set Heijnsdijk 2012
## Sensitivity (S3): Simulation by chaning RP price according to 2018 price of Karolinska Sjukhuset
## Backup sensitivity (S4): Simulation by MRI test characteristics updates (NOT USED, NOT RUN YET)

## A. discount function
plot(0:100,(1/1.03)^(0:100), type="l")
discounted <- function(object,discountRate) {
  LE <- local({
    pt <- xtabs(pt~age,data=object$summary$pt)
    ages <- as.numeric(names(pt))+0.5 # mid-point approximation
    sum((1/(1+discountRate))^ages*pt)/object$n
  })
  societal.costs <- local({
    costs <- xtabs(costs~age,data=object$societal.costs)
    ages <- as.numeric(names(costs))+0.5 # mid-point approximation
    undiscounted.costs <- costs*(1+object$simulation.parameters$discountRate.costs)^ages
    sum((1/(1+discountRate))^ages*undiscounted.costs)/object$n
  })
  healthsector.costs <- local({
    costs <- xtabs(costs~age,data=object$healthsector.costs)
    ages <- as.numeric(names(costs))+0.5 # mid-point approximation
    undiscounted.costs <- costs*(1+object$simulation.parameters$discountRate.costs)^ages
    sum((1/(1+discountRate))^ages*undiscounted.costs)/object$n
  })
  utilities <- local({
    utilities <- xtabs(ut~age,data=object$summary$ut)
    ages <- as.numeric(names(utilities))+0.5 # mid-point approximation
    undiscounted.utilities <- 
      utilities*(1+object$simulation.parameters$discountRate.effectiveness)^ages
    sum((1/(1+discountRate))^ages*undiscounted.utilities)/object$n
  })
  list(utilities=utilities,
       societal.costs=societal.costs,
       healthsector.costs=healthsector.costs,
       LE=LE)
}

#Example
# discounted(fit1,0)
# discounted(fit1,0.03)
# discounted(fit1,0.05)
# discounted.costs(fit1,0)
# discounted.costs(fit1,0.05)


## B.Events summary 
## check the event names
table(fit1$summary$events$event)

## check that we got the names right
xtabs(n~event,
      data=subset(fit2$summary$events, 
                  event %in% c("toClinicalDiagnosticBiopsy","toScreenInitiatedBiopsy")))

myreport <- function(object,names,per=1e4,subset=TRUE) {
  sum(subset(object$summary$events, event %in% names & subset)$n)/object$n*per
}

#Example
myreport(fit1,c("toClinicalDiagnosticBiopsy","toScreenInitiatedBiopsy"))
myreport(fit2,c("toClinicalDiagnosticBiopsy","toScreenInitiatedBiopsy"))
myreport(fit3,c("toClinicalDiagnosticBiopsy","toScreenInitiatedBiopsy"))

myreport(fit1,"toMRI")
myreport(fit2,"toMRI")
myreport(fit3,"toMRI")

myreport(fit1,c("toClinicalDiagnosis","toScreenDiagnosis"))
myreport(fit2,c("toClinicalDiagnosis","toScreenDiagnosis"))
myreport(fit3,c("toClinicalDiagnosis","toScreenDiagnosis"))

myreport(fit1,c("toClinicalDiagnosis","toScreenDiagnosis"),
         subset=(fit1$summary$event$age>=55 & fit1$summary$event$age<70))
myreport(fit2,c("toClinicalDiagnosis","toScreenDiagnosis"),
         subset=(fit2$summary$event$age>=55 & fit2$summary$event$age<70))
myreport(fit3,c("toClinicalDiagnosis","toScreenDiagnosis"),
         subset=(fit3$summary$event$age>=55 & fit3$summary$event$age<70))


# BASE CASE(S2): HSV sets by using PORPUS-U
ShuangParametersPorpusU <- prostata:::ShuangParameters

library(prostata)
n.sim <- 1e7
mc.cores <- 30
summary(fit1.S2 <- callFhcrc(n.sim, screen="noScreening", mc.cores=mc.cores,
                             flatPop=TRUE,pop=1995-55,
                             parms=c(ShuangParametersPorpusU,
                                     formal_compliance=1)))
summary(fit2.S2 <- callFhcrc(n.sim, screen="regular_screen",
                             flatPop=TRUE,pop=1995-55,
                             parms=c(ShuangParametersPorpusU,
                                     formal_compliance=1,
                                     start_screening=55,
                                     screening_interval=4),
                             mc.cores=mc.cores))
summary(fit3.S2 <- callFhcrc(n.sim, screen="regular_screen",
                             flatPop=TRUE,pop=1995-55,
                             parms=c(ShuangParametersPorpusU,
                                     MRI_screen=TRUE,
                                     start_screening=55,
                                     screening_interval=4,
                                     formal_compliance=1),
                             mc.cores=mc.cores))

ICER(fit2.S2,fit1.S2)
ICER(fit3.S2,fit2.S2)
ICER(fit3.S2,fit1.S2)

d.S2 <- d2.S2 <- rbind("Un-screened"=unlist(ICER(fit1.S2,fit2.S2)[2:3]),
                       "PSA+SBx"=c(0,0),
                       "PSA+MRI+TBx/SBx"=unlist(ICER(fit3.S2,fit2.S2)[2:3]))
d2.S2[,"delta.QALE"] <- d.S2[,"delta.QALE"]*10000
plot(d2.S2, xlim=1.05*range(d2.S2[,1]), ylim=1.05*range(d2.S2[,2]),
     xlab="Incremental QALYs (years) per 10000 men", ylab="Incremental costs (€)")
text(d2.S2,labels=rownames(d2.S2),pos=c(4,2,2))
abline(v=0,lty=2)
abline(h=0, lty=2)

d.S2 <- d2.S2 <- rbind("Un-screened"=unlist(ICER(fit1.S2,fit2.S2)[2:3]),
                       "PSA+SBx"=c(0,0),
                       "PSA+MRI+TBx/SBx"=unlist(ICER(fit3.S2,fit2.S2)[2:3]))
d2.S2[,"delta.QALE"] <- d.S2[,"delta.QALE"]*10000
plot(d2.S2, xlim=c(-1.6,0.4), ylim=c(-120,60),
     xlab="Incremental QALYs (years) per 10000 men", ylab="Incremental costs (€)",
     main="Cost-effectiveness plane")
text(d2.S2,labels=rownames(d2.S2),pos=c(4,2,4))
abline(v=0,lty=2)
abline(h=0, lty=2)

save(fit1.S2, file="fit1.S2_PORPUS-U.RData")
save(fit2.S2, file="fit2.S2_PORPUS-U.RData")
save(fit3.S2, file="fit3.S2_PORPUS-U.RData")

## BASE CASE(S2): Save the results (currently the base case)
sink("output_PORPUSU.txt")
print(summary(fit1.S2))
print(summary(fit2.S2))
print(summary(fit3.S2))
print(ICER(fit2.S2,fit1.S2))
print(ICER(fit3.S2,fit2.S2))
print(ICER(fit3.S2,fit1.S2))
sink()

## BASE CASE(S2): Discounted LY and costs for S2 - PORPUS-U (currenetly the base case)
discounted(fit1.S2,0)
discounted(fit1.S2,0.03)
discounted(fit1.S2,0.05)

discounted(fit2.S2,0)
discounted(fit2.S2,0.03)
discounted(fit2.S2,0.05)

discounted(fit3.S2,0)
discounted(fit3.S2,0.03)
discounted(fit3.S2,0.05)

sink("Output_Discounted LY and costs_S2_PORPUS-U_20191219.txt")
discounted.S2 <- matrix(c(discounted(fit1.S2,0)$LE, discounted(fit2.S2,0)$LE,discounted(fit3.S2,0)$LE,
                          discounted(fit1.S2,0.03)$LE, discounted(fit2.S2,0.03)$LE, discounted(fit3.S2,0.03)$LE,
                          discounted(fit1.S2,0.05)$LE, discounted(fit2.S2,0.05)$LE, discounted(fit3.S2,0.05)$LE,
                          discounted(fit1.S2,0)$utilities, discounted(fit2.S2,0)$utilities,discounted(fit3.S2,0)$utilities,
                          discounted(fit1.S2,0.03)$utilities, discounted(fit2.S2,0.03)$utilities, discounted(fit3.S2,0.03)$utilities,
                          discounted(fit1.S2,0.05)$utilities, discounted(fit2.S2,0.05)$utilities, discounted(fit3.S2,0.05)$utilities,
                          discounted(fit1.S2,0)$healthsector.costs, discounted(fit2.S2,0)$healthsector.costs,discounted(fit3.S2,0)$healthsector.costs,
                          discounted(fit1.S2,0.03)$healthsector.costs, discounted(fit2.S2,0.03)$healthsector.costs, discounted(fit3.S2,0.03)$healthsector.costs,
                          discounted(fit1.S2,0.05)$healthsector.costs, discounted(fit2.S2,0.05)$healthsector.costs, discounted(fit3.S2,0.05)$healthsector.costs,
                          discounted(fit1.S2,0)$societal.costs, discounted(fit2.S2,0)$societal.costs,discounted(fit3.S2,0)$societal.costs,
                          discounted(fit1.S2,0.03)$societal.costs, discounted(fit2.S2,0.03)$societal.costs, discounted(fit3.S2,0.03)$societal.costs,
                          discounted(fit1.S2,0.05)$societal.costs, discounted(fit2.S2,0.05)$societal.costs, discounted(fit3.S2,0.05)$societal.costs),
                        ncol=3, byrow=TRUE)
colnames(discounted.S2) <- c("No screening", "SBx", "MRI")
rownames(discounted.S2) <- c("LY, undiscounted", "LY, 3% discounted", "LY, 5% discounted",
                             "QALYs, undiscounted","QALYs, 3% discounted", "QALYs, 5% discounted",
                             "Health sector, undiscounted","Health sector, 3% discounted", "Health sector, 5% discounted",
                             "societal, undiscounted","societal, 3% discounted", "societal, 5% discounted")
discounted.S2 <- as.table(discounted.S2)
print(discounted.S2)
sink()

## BASE CASE(S2): Number of events including No.Biopsy, No.MRI and No.Diagnosis life time (currently the base case)
sink("Output_myreport_S2_PORPUS-U_20191219.txt")
myreport.S2 <- matrix(c(myreport(fit1.S2,"toMRI"),
                        myreport(fit2.S2,"toMRI"),
                        myreport(fit3.S2,"toMRI"),
                        myreport(fit1.S2,c("toClinicalDiagnosticBiopsy","toScreenInitiatedBiopsy")),
                        myreport(fit2.S2,c("toClinicalDiagnosticBiopsy","toScreenInitiatedBiopsy")),
                        myreport(fit3.S2,c("toClinicalDiagnosticBiopsy","toScreenInitiatedBiopsy")),
                        myreport(fit1.S2,c("toClinicalDiagnosis","toScreenDiagnosis")),
                        myreport(fit2.S2,c("toClinicalDiagnosis","toScreenDiagnosis")),
                        myreport(fit3.S2,c("toClinicalDiagnosis","toScreenDiagnosis")),
                        myreport(fit1.S2,c("toClinicalDiagnosis","toScreenDiagnosis"),
                                 subset=(fit1.S2$summary$event$age>=55 & fit1.S2$summary$event$age<70)),
                        myreport(fit2.S2,c("toClinicalDiagnosis","toScreenDiagnosis"),
                                 subset=(fit2.S2$summary$event$age>=55 & fit2.S2$summary$event$age<70)),
                        myreport(fit3.S2,c("toClinicalDiagnosis","toScreenDiagnosis"),
                                 subset=(fit3.S2$summary$event$age>=55 & fit3.S2$summary$event$age<70))),
                      ncol=3, byrow=TRUE)
colnames(myreport.S2) <- c("No screening", "SBx", "MRI")
rownames(myreport.S2) <- c("MRI", "Biopsy", "Diagnosis","Diagnosis-screening")
myreport.S2 <- as.table(myreport.S2)
print(myreport.S2)
sink()


## SENSITIVITY(S0): Simulation by HSV set EQ-5D
library(prostata)
ShuangParameters_S0 <- prostata:::ShuangParameters
ShuangParameters_S0$utility_estimates <-  c("Invitation" = 1,                    # Heijnsdijk 2012
                                            "Formal PSA" = 0.99,                 # Heijnsdijk 2012
                                            "Formal panel" = 0.99,               # Heijnsdijk 2012
                                            "Opportunistic PSA" = 0.99,          # Heijnsdijk 2012
                                            "Opportunistic panel" = 0.99,        # Heijnsdijk 2012
                                            "Biopsy" = 0.90,                     # Heijnsdijk 2012
                                            "Combined biopsy" = 0.90,            # Heijnsdijk 2012
                                            "Cancer diagnosis" = 0.80,           # Heijnsdijk 2012
                                            "Prostatectomy part 1" = 0.829,      # Hall 2015
                                            "Prostatectomy part 2" = 0.893,      # Extended review by SH (Glazener 2011, Korfarge 2005, Hall 2015)
                                            "Radiation therapy part 1" = 0.818,  # Hall 2015
                                            "Radiation therapy part 2" = 0.828,  # Extended review by SH (Korfarge 2005, Hall 2015)
                                            "Active surveillance" = 0.9,         # Loeb 2018
                                            "Postrecovery period" = 0.861,       # Extended review by SH (Torvinen 2013, Waston 2016) 
                                            "ADT+chemo" = 0.727,                 # Extended review by SH (Hall 2019, Diels 2015, Loriot 2015, Skaltsa 2014, Wu 2007, Chi 2018)
                                            "Palliative therapy" = 0.62,         # Magnus 2019
                                            "Terminal illness" = 0.40,           # Heijnsdijk 2012
                                            "Death" = 0.00)

n.sim <- 1e7
mc.cores <- 30
summary(fit1.S0 <- callFhcrc(n.sim, screen="noScreening", mc.cores=mc.cores,
                          flatPop=TRUE,pop=1995-55,
                          parms=c(ShuangParameters_S0,
                                  formal_compliance=1)))
summary(fit2.S0 <- callFhcrc(n.sim, screen="regular_screen",
                          flatPop=TRUE,pop=1995-55,
                          parms=c(ShuangParameters_S0,
                                  formal_compliance=1,
                                  start_screening=55,
                                  screening_interval=4),
                          mc.cores=mc.cores))
summary(fit3.S0 <- callFhcrc(n.sim, screen="regular_screen",
                          flatPop=TRUE,pop=1995-55,
                          parms=c(ShuangParameters_S0,
                                  MRI_screen=TRUE,
                                  start_screening=55,
                                  screening_interval=4,
                                  formal_compliance=1),
                          mc.cores=mc.cores))

# xtabs(I(costs/1e7)~item+type,fit1$societal.costs)
# xtabs(I(costs/1e7)~item+type,fit2$societal.costs)
# xtabs(I(costs/1e7)~item+type,fit3$societal.costs)

# prostata:::FhcrcParameters$cost_parameters
# prostata:::ShuangParameters$cost_parameters
# prostata:::FhcrcParameters$production
# prostata:::ShuangParameters$production
# prostata:::FhcrcParameters$lost_production_years
# prostata:::ShuangParameters$lost_production_years
# prostata:::FhcrcParameters$utility_estimates
# prostata:::ShuangParameters$utility_estimates
# prostata:::FhcrcParameters$utility_duration
# prostata:::ShuangParameters$utility_duration
# prostata:::FhcrcParameters$currency_rate
# prostata:::ShuangParameters$currency_rate

ICER(fit2.S0,fit1.S0)
ICER(fit3.S0,fit2.S0)
ICER(fit3.S0,fit1.S0)

d.S0 <- d2.S0 <- rbind("Un-screened"=unlist(ICER(fit1.S0,fit2.S0)[2:3]),
           "PSA+SBx"=c(0,0),
           "PSA+MRI+TBx/SBx"=unlist(ICER(fit3.S0,fit2.S0)[2:3]))
d2.S0[,"delta.QALE"] <- d.S0[,"delta.QALE"]*10000
plot(d2.S0, xlim=1.05*range(d2.S0[,1]), ylim=1.05*range(d2.S0[,2]),
     xlab="Incremental QALYs (years) per 10000 men", ylab="Incremental costs (€)")
text(d2.S0,labels=rownames(d2.S0),pos=c(4,2,2))
abline(v=0,lty=2)
abline(h=0, lty=2)

save(fit1.S0, file="fit1.S0_EQ-5D.RData")
save(fit2.S0, file="fit2.S0_EQ-5D.RData")
save(fit3.S0, file="fit3.S0_EQ-5D.RData")

sink("Output_EQ-5D.txt")
print(summary(fit1.S0))
print(summary(fit2.S0))
print(summary(fit3.S0))
print(ICER(fit2.S0,fit1.S0))
print(ICER(fit3.S0,fit2.S0))
print(ICER(fit3.S0,fit1.S0))
sink()

# plot(c(0,10),c(0,10),type="n")
# segments(1:2,1:2,3:4,5:6)
# points(c(1:2,3:4),c(1:2,5:6))

## SENSITIVITY(S0): Discounted LY and costs for S0 - EQ-5D
discounted(fit1.S0,0)
discounted(fit1.S0,0.03)
discounted(fit1.S0,0.05)

discounted(fit2.S0,0)
discounted(fit2.S0,0.03)
discounted(fit2.S0,0.05)

discounted(fit3.S0,0)
discounted(fit3.S0,0.03)
discounted(fit3.S0,0.05)

sink("Output_Discounted LY and costs_S0_EQ-5D.txt")
discounted.S0 <- matrix(c(discounted(fit1.S0,0)$LE, discounted(fit2.S0,0)$LE,discounted(fit3.S0,0)$LE,
                          discounted(fit1.S0,0.03)$LE, discounted(fit2.S0,0.03)$LE, discounted(fit3.S0,0.03)$LE,
                          discounted(fit1.S0,0.05)$LE, discounted(fit2.S0,0.05)$LE, discounted(fit3.S0,0.05)$LE,
                          discounted(fit1.S0,0)$utilities, discounted(fit2.S0,0)$utilities,discounted(fit3.S0,0)$utilities,
                          discounted(fit1.S0,0.03)$utilities, discounted(fit2.S0,0.03)$utilities, discounted(fit3.S0,0.03)$utilities,
                          discounted(fit1.S0,0.05)$utilities, discounted(fit2.S0,0.05)$utilities, discounted(fit3.S0,0.05)$utilities,
                          discounted(fit1.S0,0)$healthsector.costs, discounted(fit2.S0,0)$healthsector.costs,discounted(fit3.S0,0)$healthsector.costs,
                          discounted(fit1.S0,0.03)$healthsector.costs, discounted(fit2.S0,0.03)$healthsector.costs, discounted(fit3.S0,0.03)$healthsector.costs,
                          discounted(fit1.S0,0.05)$healthsector.costs, discounted(fit2.S0,0.05)$healthsector.costs, discounted(fit3.S0,0.05)$healthsector.costs,
                          discounted(fit1.S0,0)$societal.costs, discounted(fit2.S0,0)$societal.costs,discounted(fit3.S0,0)$societal.costs,
                          discounted(fit1.S0,0.03)$societal.costs, discounted(fit2.S0,0.03)$societal.costs, discounted(fit3.S0,0.03)$societal.costs,
                          discounted(fit1.S0,0.05)$societal.costs, discounted(fit2.S0,0.05)$societal.costs, discounted(fit3.S0,0.05)$societal.costs),
                        ncol=3, byrow=TRUE)
colnames(discounted.S0) <- c("No screening", "SBx", "MRI")
rownames(discounted.S0) <- c("LY, undiscounted", "LY, 3% discounted", "LY, 5% discounted",
                             "QALYs, undiscounted","QALYs, 3% discounted", "QALYs, 5% discounted",
                             "Health sector, undiscounted","Health sector, 3% discounted", "Health sector, 5% discounted",
                             "societal, undiscounted","societal, 3% discounted", "societal, 5% discounted")
discounted.S0 <- as.table(discounted.S0)
print(discounted.S0)
sink()

## SENSITIVITY(S0): Number of events including No.Biopsy, No.MRI and No.Diagnosis life time - EQ-5D
sink("Output_myreport_S0_EQ-5D.txt")
myreport.S0 <- matrix(c(myreport(fit1.S0,"toMRI"),
                        myreport(fit2.S0,"toMRI"),
                        myreport(fit3.S0,"toMRI"),
                        myreport(fit1.S0,c("toClinicalDiagnosticBiopsy","toScreenInitiatedBiopsy")),
                        myreport(fit2.S0,c("toClinicalDiagnosticBiopsy","toScreenInitiatedBiopsy")),
                        myreport(fit3.S0,c("toClinicalDiagnosticBiopsy","toScreenInitiatedBiopsy")),
                        myreport(fit1.S0,c("toClinicalDiagnosis","toScreenDiagnosis")),
                        myreport(fit2.S0,c("toClinicalDiagnosis","toScreenDiagnosis")),
                        myreport(fit3.S0,c("toClinicalDiagnosis","toScreenDiagnosis")),
                        myreport(fit1.S0,c("toClinicalDiagnosis","toScreenDiagnosis"),
                                 subset=(fit1.S0$summary$event$age>=55 & fit1.S0$summary$event$age<70)),
                        myreport(fit2.S0,c("toClinicalDiagnosis","toScreenDiagnosis"),
                                 subset=(fit2.S0$summary$event$age>=55 & fit2.S0$summary$event$age<70)),
                        myreport(fit3.S0,c("toClinicalDiagnosis","toScreenDiagnosis"),
                                 subset=(fit3.S0$summary$event$age>=55 & fit3.S0$summary$event$age<70))),
                      ncol=3, byrow=TRUE)
colnames(myreport.S0) <- c("No screening", "SBx", "MRI")
rownames(myreport.S0) <- c("MRI", "Biopsy", "Diagnosis","Diagnosis-screening")
myreport.S0 <- as.table(myreport.S0)
print(myreport.S0)
sink()


## SENSITIVITY(S1). Changing HSV set by using Heijnsdijk 2012 
ShuangParameters_S1 <- prostata:::ShuangParameters
ShuangParameters_S1$utility_estimates <-  c("Invitation" = 1,                    # Heijnsdijk 2012
                                            "Formal PSA" = 0.99,                 # Heijnsdijk 2012
                                            "Formal panel" = 0.99,               # Heijnsdijk 2012
                                            "Opportunistic PSA" = 0.99,          # Heijnsdijk 2012
                                            "Opportunistic panel" = 0.99,        # Heijnsdijk 2012
                                            "Biopsy" = 0.90,                     # Heijnsdijk 2012
                                            "Cancer diagnosis" = 0.80,           # Heijnsdijk 2012
                                            "Prostatectomy part 1" = 0.67,       # Heijnsdijk 2012
                                            "Prostatectomy part 2" = 0.77,       # Heijnsdijk 2012
                                            "Radiation therapy part 1" = 0.73,   # Heijnsdijk 2012
                                            "Radiation therapy part 2" = 0.78,   # Heijnsdijk 2012
                                            "Active surveillance" = 0.97,        # Heijnsdijk 2012
                                            "Postrecovery period" = 0.95,        # Heijnsdijk 2012 
                                            "ADT+chemo" = 0.727,                 # Review by Shuang
                                            "Palliative therapy" = 0.60,         # Heijnsdijk 2012
                                            "Terminal illness" = 0.40,           # Heijnsdijk 2012
                                            "Death" = 0.00)

n.sim <- 1e7
mc.cores <- 30
summary(fit1.S1 <- callFhcrc(n.sim, screen="noScreening", mc.cores=mc.cores,
                          flatPop=TRUE,pop=1995-55,
                          parms=c(ShuangParameters_S1,
                                  formal_compliance=1)))
summary(fit2.S1 <- callFhcrc(n.sim, screen="regular_screen",
                          flatPop=TRUE,pop=1995-55,
                          parms=c(ShuangParameters_S1,
                                  formal_compliance=1,
                                  start_screening=55,
                                  screening_interval=4),
                          mc.cores=mc.cores))
summary(fit3.S1 <- callFhcrc(n.sim, screen="regular_screen",
                          flatPop=TRUE,pop=1995-55,
                          parms=c(ShuangParameters_S1,
                                  MRI_screen=TRUE,
                                  start_screening=55,
                                  screening_interval=4,
                                  formal_compliance=1),
                          mc.cores=mc.cores))

ICER(fit2.S1,fit1.S1)
ICER(fit3.S1,fit2.S1)
ICER(fit3.S1,fit1.S1)

d.S1 <- d2.S1 <- rbind("Un-screened"=unlist(ICER(fit1.S1,fit2.S1)[2:3]),
                 "PSA+SBx"=c(0,0),
                 "PSA+MRI+TBx/SBx"=unlist(ICER(fit3.S1,fit2.S1)[2:3]))
d2.S1[,"delta.QALE"] <- d.S1[,"delta.QALE"]*10000
plot(d2.S1, xlim=1.05*range(d2.S1[,1]), ylim=1.05*range(d2.S1[,2]),
     xlab="Incremental QALYs (years) per 10000 men", ylab="Incremental costs (€)")
text(d2.S1,labels=rownames(d2.S1),pos=c(4,2,2))
abline(v=0,lty=2)
abline(h=0, lty=2)

# d.S1 <- d2.S1 <- rbind("PSA+SBx"=unlist(ICER(fit2.S1,fit1.S1)[2:3]),
#                        "Un-screened"=c(0,0), #unscreened as reference
#                        "PSA+MRI+TBx/SBx"=unlist(ICER(fit3.S1,fit1.S1)[2:3]))
# d2.S1[,"delta.QALE"] <- d.S1[,"delta.QALE"]*10000
# plot(d2.S1, xlim=1.05*range(d2.S1[,1]), ylim=1.05*range(d2.S1[,2]),
#      xlab="Incremental QALYs (years) per 10000 men", ylab="Incremental costs (€)")
# text(d2.S1T,labels=rownames(d2.S1T),pos=c(4,2,2))
# abline(v=0,lty=2)
# abline(h=0, lty=2)

save(fit1.S1, file="fit1.S1_Heijnsdijk 2012.RData")
save(fit2.S1, file="fit2.S1_Heijnsdijk 2012.RData")
save(fit3.S1, file="fit3.S1_Heijnsdijk 2012.RData")

## SENSITIVITY(S1): Save the results - Heijnsdijk 2012
sink("Output_Heijnsdijk 2012 HSV.txt")
print(summary(fit1.S1))
print(summary(fit2.S1))
print(summary(fit3.S1))
print(ICER(fit2.S1,fit1.S1))
print(ICER(fit3.S1,fit2.S1))
print(ICER(fit3.S1,fit1.S1))
sink()

## SENSITIVITY(S1): Discounted LY and costs for S1 - Heijnsdijk 2012
discounted(fit1.S1,0)
discounted(fit1.S1,0.03)
discounted(fit1.S1,0.05)

discounted(fit2.S1,0)
discounted(fit2.S1,0.03)
discounted(fit2.S1,0.05)

discounted(fit3.S1,0)
discounted(fit3.S1,0.03)
discounted(fit3.S1,0.05)

sink("Output_Discounted LY and costs_S1_Heijnsdijk 2012 HSV.txt")
discounted.S1 <- matrix(c(discounted(fit1.S1,0)$LE, discounted(fit2.S1,0)$LE,discounted(fit3.S1,0)$LE,
                          discounted(fit1.S1,0.03)$LE, discounted(fit2.S1,0.03)$LE, discounted(fit3.S1,0.03)$LE,
                          discounted(fit1.S1,0.05)$LE, discounted(fit2.S1,0.05)$LE, discounted(fit3.S1,0.05)$LE,
                          discounted(fit1.S1,0)$utilities, discounted(fit2.S1,0)$utilities,discounted(fit3.S1,0)$utilities,
                          discounted(fit1.S1,0.03)$utilities, discounted(fit2.S1,0.03)$utilities, discounted(fit3.S1,0.03)$utilities,
                          discounted(fit1.S1,0.05)$utilities, discounted(fit2.S1,0.05)$utilities, discounted(fit3.S1,0.05)$utilities,
                          discounted(fit1.S1,0)$healthsector.costs, discounted(fit2.S1,0)$healthsector.costs,discounted(fit3.S1,0)$healthsector.costs,
                          discounted(fit1.S1,0.03)$healthsector.costs, discounted(fit2.S1,0.03)$healthsector.costs, discounted(fit3.S1,0.03)$healthsector.costs,
                          discounted(fit1.S1,0.05)$healthsector.costs, discounted(fit2.S1,0.05)$healthsector.costs, discounted(fit3.S1,0.05)$healthsector.costs,
                          discounted(fit1.S1,0)$societal.costs, discounted(fit2.S1,0)$societal.costs,discounted(fit3.S1,0)$societal.costs,
                          discounted(fit1.S1,0.03)$societal.costs, discounted(fit2.S1,0.03)$societal.costs, discounted(fit3.S1,0.03)$societal.costs,
                          discounted(fit1.S1,0.05)$societal.costs, discounted(fit2.S1,0.05)$societal.costs, discounted(fit3.S1,0.05)$societal.costs),
                        ncol=3, byrow=TRUE)
colnames(discounted.S1) <- c("No screening", "SBx", "MRI")
rownames(discounted.S1) <- c("LY, undiscounted", "LY, 3% discounted", "LY, 5% discounted",
                             "QALYs, undiscounted","QALYs, 3% discounted", "QALYs, 5% discounted",
                             "Health sector, undiscounted","Health sector, 3% discounted", "Health sector, 5% discounted",
                             "societal, undiscounted","societal, 3% discounted", "societal, 5% discounted")
discounted.S1 <- as.table(discounted.S1)
print(discounted.S1)
sink()

## SENSITIVITY(S1): Number of events including No.Biopsy, No.MRI and No.Diagnosis life time - Heijnsdijk 2012
sink("Output_myreport_S1_Heijnsdijk 2012 HSV_20191219.txt")
myreport.S1 <- matrix(c(myreport(fit1.S1,"toMRI"),
                        myreport(fit2.S1,"toMRI"),
                        myreport(fit3.S1,"toMRI"),
                        myreport(fit1.S1,c("toClinicalDiagnosticBiopsy","toScreenInitiatedBiopsy")),
                        myreport(fit2.S1,c("toClinicalDiagnosticBiopsy","toScreenInitiatedBiopsy")),
                        myreport(fit3.S1,c("toClinicalDiagnosticBiopsy","toScreenInitiatedBiopsy")),
                        myreport(fit1.S1,c("toClinicalDiagnosis","toScreenDiagnosis")),
                        myreport(fit2.S1,c("toClinicalDiagnosis","toScreenDiagnosis")),
                        myreport(fit3.S1,c("toClinicalDiagnosis","toScreenDiagnosis")),
                        myreport(fit1.S1,c("toClinicalDiagnosis","toScreenDiagnosis"),
                                 subset=(fit1.S1$summary$event$age>=55 & fit1.S1$summary$event$age<70)),
                        myreport(fit2.S1,c("toClinicalDiagnosis","toScreenDiagnosis"),
                                 subset=(fit2.S1$summary$event$age>=55 & fit2.S1$summary$event$age<70)),
                        myreport(fit3.S1,c("toClinicalDiagnosis","toScreenDiagnosis"),
                                 subset=(fit3.S1$summary$event$age>=55 & fit3.S1$summary$event$age<70))),
                      ncol=3, byrow=TRUE)
colnames(myreport.S1) <- c("No screening", "SBx", "MRI")
rownames(myreport.S1) <- c("MRI", "Biopsy", "Diagnosis","Diagnosis-screening")
myreport.S1 <- as.table(myreport.S1)
print(myreport.S1)
sink()


## SENSITIVITY(S3). Changing the price of prostatectomy to 96,438 SEK 
library(prostata)
ShuangParameters_S3 <- prostata:::ShuangParameters
ShuangParameters_S3$cost_parameters["Prostatectomy"] <-  96438+6302.19*20*0.25+1460*1

n.sim <- 1e7
mc.cores <- 30
summary(fit1.S3 <- callFhcrc(n.sim, screen="noScreening", mc.cores=mc.cores,
                             flatPop=TRUE,pop=1995-55,
                             parms=c(ShuangParameters_S3,
                                     formal_compliance=1)))
summary(fit2.S3 <- callFhcrc(n.sim, screen="regular_screen",
                             flatPop=TRUE,pop=1995-55,
                             parms=c(ShuangParameters_S3,
                                     formal_compliance=1,
                                     start_screening=55,
                                     screening_interval=4),
                             mc.cores=mc.cores))
summary(fit3.S3 <- callFhcrc(n.sim, screen="regular_screen",
                             flatPop=TRUE,pop=1995-55,
                             parms=c(ShuangParameters_S3,
                                     MRI_screen=TRUE,
                                     start_screening=55,
                                     screening_interval=4,
                                     formal_compliance=1),
                             mc.cores=mc.cores))

ICER(fit2.S3,fit1.S3)
ICER(fit3.S3,fit2.S3)
ICER(fit3.S3,fit1.S3)

d.S3 <- d2.S3 <- rbind("Un-screened"=unlist(ICER(fit1.S3,fit2.S3)[2:3]),
                       "PSA+SBx"=c(0,0),
                       "PSA+MRI+TBx/SBx"=unlist(ICER(fit3.S3,fit2.S3)[2:3]))
d2.S3[,"delta.QALE"] <- d.S3[,"delta.QALE"]*10000
plot(d2.S3, xlim=1.05*range(d2.S3[,1]), ylim=1.05*range(d2.S3[,2]),
     xlab="Incremental QALYs (years) per 10000 men", ylab="Incremental costs (€)")
text(d2.S3,labels=rownames(d2.S3),pos=c(4,2,2))
abline(v=0,lty=2)
abline(h=0, lty=2)

save(fit1.S3, file="fit1.S3_RP price update.RData")
save(fit2.S3, file="fit2.S3_RP price update.RData")
save(fit3.S3, file="fit3.S3_RP price update.RData")

## SENSITIVITY(S3): Save the results - RP price updated
sink("Output_RP price update.txt")
print(summary(fit1.S3))
print(summary(fit2.S3))
print(summary(fit3.S3))
print(ICER(fit2.S3,fit1.S3))
print(ICER(fit3.S3,fit2.S3))
print(ICER(fit3.S3,fit1.S3))
sink()

## SENSITIVITY(S3): Discounted LY and costs for S3 - RP price updated
discounted(fit1.S3,0)
discounted(fit1.S3,0.03)
discounted(fit1.S3,0.05)

discounted(fit2.S3,0)
discounted(fit2.S3,0.03)
discounted(fit2.S3,0.05)

discounted(fit3.S3,0)
discounted(fit3.S3,0.03)
discounted(fit3.S3,0.05)

sink("Output_Discounted LY and costs_S3_RP price update.txt")
discounted.S3 <- matrix(c(discounted(fit1.S3,0)$LE, discounted(fit2.S3,0)$LE,discounted(fit3.S3,0)$LE,
                          discounted(fit1.S3,0.03)$LE, discounted(fit2.S3,0.03)$LE, discounted(fit3.S3,0.03)$LE,
                          discounted(fit1.S3,0.05)$LE, discounted(fit2.S3,0.05)$LE, discounted(fit3.S3,0.05)$LE,
                          discounted(fit1.S3,0)$utilities, discounted(fit2.S3,0)$utilities,discounted(fit3.S3,0)$utilities,
                          discounted(fit1.S3,0.03)$utilities, discounted(fit2.S3,0.03)$utilities, discounted(fit3.S3,0.03)$utilities,
                          discounted(fit1.S3,0.05)$utilities, discounted(fit2.S3,0.05)$utilities, discounted(fit3.S3,0.05)$utilities,
                          discounted(fit1.S3,0)$healthsector.costs, discounted(fit2.S3,0)$healthsector.costs,discounted(fit3.S3,0)$healthsector.costs,
                          discounted(fit1.S3,0.03)$healthsector.costs, discounted(fit2.S3,0.03)$healthsector.costs, discounted(fit3.S3,0.03)$healthsector.costs,
                          discounted(fit1.S3,0.05)$healthsector.costs, discounted(fit2.S3,0.05)$healthsector.costs, discounted(fit3.S3,0.05)$healthsector.costs,
                          discounted(fit1.S3,0)$societal.costs, discounted(fit2.S3,0)$societal.costs,discounted(fit3.S3,0)$societal.costs,
                          discounted(fit1.S3,0.03)$societal.costs, discounted(fit2.S3,0.03)$societal.costs, discounted(fit3.S3,0.03)$societal.costs,
                          discounted(fit1.S3,0.05)$societal.costs, discounted(fit2.S3,0.05)$societal.costs, discounted(fit3.S3,0.05)$societal.costs),
                        ncol=3, byrow=TRUE)
colnames(discounted.S3) <- c("No screening", "SBx", "MRI")
rownames(discounted.S3) <- c("LY, undiscounted", "LY, 3% discounted", "LY, 5% discounted",
                             "QALYs, undiscounted","QALYs, 3% discounted", "QALYs, 5% discounted",
                             "Health sector, undiscounted","Health sector, 3% discounted", "Health sector, 5% discounted",
                             "societal, undiscounted","societal, 3% discounted", "societal, 5% discounted")
discounted.S3 <- as.table(discounted.S3)
print(discounted.S3)
sink()

## SENSITIVITY(S3): Number of events including No.Biopsy, No.MRI and No.Diagnosis life time - RP price update
sink("Output_myreport_S3_RP price update.txt")
myreport.S3 <- matrix(c(myreport(fit1.S3,"toMRI"),
                        myreport(fit2.S3,"toMRI"),
                        myreport(fit3.S3,"toMRI"),
                        myreport(fit1.S3,c("toClinicalDiagnosticBiopsy","toScreenInitiatedBiopsy")),
                        myreport(fit2.S3,c("toClinicalDiagnosticBiopsy","toScreenInitiatedBiopsy")),
                        myreport(fit3.S3,c("toClinicalDiagnosticBiopsy","toScreenInitiatedBiopsy")),
                        myreport(fit1.S3,c("toClinicalDiagnosis","toScreenDiagnosis")),
                        myreport(fit2.S3,c("toClinicalDiagnosis","toScreenDiagnosis")),
                        myreport(fit3.S3,c("toClinicalDiagnosis","toScreenDiagnosis")),
                        myreport(fit1.S3,c("toClinicalDiagnosis","toScreenDiagnosis"),
                                 subset=(fit1.S3$summary$event$age>=55 & fit1.S3$summary$event$age<70)),
                        myreport(fit2.S3,c("toClinicalDiagnosis","toScreenDiagnosis"),
                                 subset=(fit2.S3$summary$event$age>=55 & fit2.S3$summary$event$age<70)),
                        myreport(fit3.S3,c("toClinicalDiagnosis","toScreenDiagnosis"),
                                 subset=(fit3.S3$summary$event$age>=55 & fit3.S3$summary$event$age<70))),
                      ncol=3, byrow=TRUE)
colnames(myreport.S3) <- c("No screening", "SBx", "MRI")
rownames(myreport.S3) <- c("MRI", "Biopsy", "Diagnosis","Diagnosis-screening")
myreport.S3 <- as.table(myreport.S3)
print(myreport.S3)
sink()



# BACKUP SENSITIVITY (S4). Changing MRI related test characteristics
ShuangParameters_S4 <- prostata:::ShuangParameters
ShuangParameters_S4$pMRIposG0=0.457               # Pr(MRI+ | ISUP 0 || undetectable)
ShuangParameters_S4$pMRIposG1=0.669               # Pr(MRI+ | ISUP 1 && detectable)
ShuangParameters_S4$pMRIposG2=0.939               # Pr(MRI+ | ISUP 2+ && detectable)
ShuangParameters_S4$pSBxG0ifG1=0.108              # Pr(SBx gives ISUP 0 | ISUP 1)
ShuangParameters_S4$pSBxG0ifG2=0.096              # Pr(SBx gives ISUP 0 | ISUP 2)

n.sim <- 1e7
mc.cores <- 1
summary(fit1.S4 <- callFhcrc(n.sim, screen="noScreening", mc.cores=mc.cores,
                             flatPop=TRUE,pop=1995-55,
                             parms=c(ShuangParameters_S4,
                                     formal_compliance=1)))
summary(fit2.S4 <- callFhcrc(n.sim, screen="regular_screen",
                             flatPop=TRUE,pop=1995-55,
                             parms=c(ShuangParameters_S4,
                                     formal_compliance=1,
                                     start_screening=55,
                                     screening_interval=4),
                             mc.cores=mc.cores))
summary(fit3.S4 <- callFhcrc(n.sim, screen="regular_screen",
                             flatPop=TRUE,pop=1995-55,
                             parms=c(ShuangParameters_S4,
                                     MRI_screen=TRUE,
                                     start_screening=55,
                                     screening_interval=4,
                                     formal_compliance=1),
                             mc.cores=mc.cores))

ICER(fit2.S4,fit1.S4)
ICER(fit3.S4,fit2.S4)
ICER(fit3.S4,fit1.S4)

d.S4 <- d2.S4 <- rbind("Un-screened"=unlist(ICER(fit1.S4,fit2.S4)[2:3]),
                       "PSA+SBx"=c(0,0),
                       "PSA+MRI+TBx/SBx"=unlist(ICER(fit3.S3,fit2.S4)[2:3]))
d2.S4[,"delta.QALE"] <- d.S4[,"delta.QALE"]*1000
plot(d2.S4, xlim=1.05*range(d2.S4[,1]), ylim=1.05*range(d2.S4[,2]),
     xlab="Incremental QALYs (years) per 1000 men", ylab="Incremental costs (€)")
text(d2.S4,labels=rownames(d2.S4),pos=c(4,2,2))
abline(v=0,lty=2)
abline(h=0, lty=2)

save(fit1.S4, file="fit1.S4_MRI test characteristics") #NOT RUN YET, NOT SAVED
save(fit2.S4, file="fit2.S4_MRI test characteristics") #NOT RUN YET, NOT SAVED
save(fit3.S4, file="fit3.S4_MRI test characteristics") #NOT RUN YET, NOT SAVED

sink("output_MRI test characteristics.txt")
print(summary(fit1.S4))
print(summary(fit2.S4))
print(summary(fit3.S4))
print(ICER(fit1.S4,fit2.S4))
print(ICER(fit3.S4,fit2.S4))
print(ICER(fit3.S4,fit1.S4))
sink()


# 
library(prostata)
n.sim <- 1e5
mc.cores <- 1
summary(fit1 <- callFhcrc(n.sim, screen="noScreening", mc.cores=mc.cores,
                          flatPop=TRUE,pop=1995-55,
                          parms=c(Parameters,
                                  formal_compliance=1)))
summary(fit2 <- callFhcrc(n.sim, screen="regular_screen",
                          flatPop=TRUE,pop=1995-55,
                          parms=c(Parameters,
                                  formal_compliance=1,
                                  start_screening=55,
                                  screening_interval=4),
                          mc.cores=mc.cores))
ICER(fit2,fit1)

## Monte Carlo uncertainty in the ICER
mc.uncertainty <- function(fit1,fit2,R=1000,do.boot=FALSE) {
    out <- local({
        X <- fit2$indiv_costs-fit1$indiv_costs
        Y <- fit2$indiv_utilities-fit1$indiv_utilities
        muX <- mean(X)
        muY <- mean(Y)
        n <- length(X)
        list(icer=muX/muY,
             sd_log_ratio=sqrt((sd(X)/muX/sqrt(n-1))^2+(sd(Y)/muY/sqrt(n-1))^2))
    })
    out$icer.lower = with(out, icer*exp(-1.96*sd_log_ratio))
    out$icer.upper = with(out, icer*exp(1.96*sd_log_ratio))
    if (do.boot) {
        stopifnot(require(boot))
        d <- data.frame(x=fit2$indiv_costs-fit1$indiv_costs,
                        y=fit2$indiv_utilities-fit1$indiv_utilities)
        boot1 <- boot(d,
                      function(d, w) log(mean(d$x[w])/mean(d$y[w])),
                      R=R)
        out <- c(list(boot=boot1),out)
    }
    out
}
mc.uncertainty(fit1.S2,fit2.S2)
mc.uncertainty(fit1.S2,fit3.S2)
plot(density(fit1$indiv_utilities-fit2$indiv_utilities))
plot(density(fit1$indiv_costs-fit2$indiv_costs))
t.test(fit2$indiv_costs-fit1$indiv_costs)
t.test(fit2$indiv_utilities-fit1$indiv_utilities)

xtabs(pt~age,data=fit1$summary$pt)

## Check using the old code (log(50000)=10.8)
library(prostata)
mc.cores <- 10
n <- 1e5
fit1.old <- callFhcrc(n, mc.cores=mc.cores,flatPop=TRUE,pop=1995-55,
                 parms=list(formal_compliance=1))
fit2.old <- callFhcrc(n, mc.cores=mc.cores,flatPop=TRUE,pop=1995-55,
                  parms=list(formal_compliance=1,
                             start_screening=55,
                             screening_interval=4),
                  screen="regular_screen")
summary(fit1.old)
summary(fit2.old)
ICER(fit2.old,fit1.old)
xtabs(I(costs/1e5)~item+type,fit1.old$societal.costs)
xtabs(I(costs/1e5)~item+type,fit2.old$societal.costs)
library(boot)
boot(data.frame(x=fit2$indiv_costs-fit1$indiv_costs,
                y=fit2$indiv_utilities-fit1$indiv_utilities),
     function(d, w) log(mean(d$x[w])/mean(d$y[w])),
     R=1000)


xtabs(I(costs/1e5)~item+type,fit1.old$societal.costs)
xtabs(I(costs/1e5)~item+type,fit1$societal.costs)
##
xtabs(I(costs/1e5)~item+type,fit2.old$societal.costs)
xtabs(I(costs/1e5)~item+type,fit2$societal.costs)

xtabs(I(ut/1e5)~dx,fit1.old$summary$ut)
xtabs(I(ut/1e5)~dx,fit2.old$summary$ut)


## proportion of deaths that are prostate cancer
other=sum(subset(fit1$summary$events,event=="toOtherDeath")$n)
pca=sum(subset(fit1$summary$events,event=="toCancerDeath")$n)
pca/(pca+other)
other=sum(subset(fit2$summary$events,event=="toOtherDeath")$n)
pca=sum(subset(fit2$summary$events,event=="toCancerDeath")$n)
pca/(pca+other)

refresh
require(foreign)
require(rstpm2)
require(dplyr)
require(ggplot2)
require(lattice)
require(maxLik)
##
temp <- read.dta("~/Downloads/rescreen-20141024.dta")
temp2 <- temp %>%
    filter((is.na(diadate) | sample_date<diadate) & time>0) %>% # restrict to non-cancers to PSA before cancer
    mutate(total_cat=cut(total,c(0,1,3,10,Inf))) # generate PSA categories
temp4 <- temp2 %>%
    mutate(age_start = as.numeric((sample_date-dob)/365),
           age_finish = as.numeric((pmin(next_sample_date,finish,na.rm=TRUE)-dob)/365)) %>%
    as.data.frame()
temp4 <- filter(temp4, !is.na(total_cat)) %>%
    mutate(time = time/365.25)
temp2b <- read.dta("/home/marcle/Downloads/psa_rates_nopca_20141017.dta") # (#test | no pca)
## xyplot(psa_rate ~ year | factor(age5), data=temp2b, subset=age5>=40 & year<=2012, type="l")
## Results are moderately stable for 2008-2011
expit <- rstpm2:::expit
Call <- function(args,fun) do.call(deparse(substitute(fun)),as.list(args))
##

cureModel <- function(formula,data=NULL,formula.shape=~1,formula.cure=~1,par=NULL,maxLik=FALSE,...) {
    formula.scale <- formula
    rstpm2:::lhs(formula.scale) <- NULL
    mf <- model.frame(formula,data)
    y <- model.extract(mf,"response")
    if (ncol(y)==2) { # right censored
        delayed <- FALSE
        time <- y[,1]
        status <- y[,2]
    }
    if (ncol(y)==3) { # left truncated and right censored
        entry <- y[,1]
        time <- y[,2]
        status <- y[,3]
        if (delayed <- any(entry>0))
            i <- which(entry>0)
    }
    formulae <- list(shape=formula.shape,
                     scale=formula.scale,
                     cure=formula.cure)
    X <- lapply(formulae,model.matrix,data)
    stopifnot(nrow(X$scale) == nrow(X$shape))
    n <- lapply(X,ncol)
    cumn <- cumsum(n)
    index <- mapply(function(a,b) (a+1):b, c(0,cumn[-length(cumn)]),cumn,SIMPLIFY=FALSE)
    names(index) <- names(formulae)
    ll <-function(beta) {
        lp <- mapply(function(x,idx) as.vector(x %*% beta[idx]), X, index, SIMPLIFY=FALSE)
        shape <- exp(lp$shape)
        scale <- exp(lp$scale)
        cure <- rstpm2:::expit(lp$cure)
        S <- cure+(1-cure)*pweibull(time,shape,scale,lower.tail=FALSE)
        f <- (1-cure)*dweibull(time,shape,scale)
        value <- sum(ifelse(status==1, log(f), log(S)))
        if (delayed) {
            S0 <- cure[i]+(1-cure[i])*pweibull(entry[i],shape[i],scale[i],lower.tail=FALSE)
            value <- value - sum(log(S0))
        }
        value
    }
    ## logLik <-function(beta) {
    ##     lp <- mapply(function(x,idx) as.vector(x %*% beta[idx]), X, index, SIMPLIFY=FALSE)
    ##     shape <- exp(lp$shape)
    ##     scale <- exp(lp$scale)
    ##     cure <- rstpm2:::expit(lp$cure)
    ##     S <- cure+(1-cure)*pweibull(time,shape,scale,lower.tail=FALSE)
    ##     f <- (1-cure)*dweibull(time,shape,scale)
    ##     value <- ifelse(status==1, log(f), log(S))
    ##     if (delayed) {
    ##         S0 <- cure[i]+(1-cure[i])*pweibull(entry[i],shape[i],scale[i],lower.tail=FALSE)
    ##         value[i] <- value[i] + log(S0)
    ##     }
    ##     value
    ## }
    if (is.null(par)) {
        aft <- survreg(formula,data)
        ## AFT -> Weibull
        sigma <- aft$scale # AFT "scale"
        shape <- 1/sigma # Weibull "shape"
        logscale <- coef(aft)[1]
        logbeta <- (coef(aft)[-1])
        par <- c(shape=c(log(shape),rep(0,n$shape-1)),
                 scale=c(logscale,logbeta),
                 cure=c(-7,rep(0,n$cure-1)))
        print(par)
        ## par <- c(shape=rep(0,n.shape), # shape
        ##                scale=c(log(sum(status)/sum(time)),rep(0,n.scale-1)), # scale
        ##                cure=c(logit((1-mean(status))/2), rep(0,n.cure-1)))
    }
    if (maxLik) obj <- maxLik(ll, start = par, ...) else {
        fit <- .Call("fitCureModel",time,status,X$shape,X$scale,X$cure,par,package="rstpm2") # fast
        coef <- fit$coef
        obj <- list(maximum=-fit$Fmin,
                    estimate=coef,
                    gradient=rstpm2:::grad(ll,coef),
                    gradientObs=matrix(NA,length(time),length(coef)), # hack
                    hessian=fit$hessian,
                    code=if (fit$fail==0) 1 else 4,
                    message = "",
                    last.step=NULL,
                    fixed = rep(FALSE,length(coef)),
                    iterations=0,
                    type="",
                    constraints=NULL,
                    varcovar=solve(fit$hessian))
        class(obj) <- c("maxLik","maxim","list")
        ## ## alternative approach using bbmle::mle2
        ## negll <- function(beta) -ll(beta)
        ## rownames(fit$hessian) <- colnames(fit$hessian) <- parnames(negll) <- names(coef)
        ## obj <- mle2(negll, start = par, eval.only = TRUE, vec.par = TRUE) # SLOW
        ## obj@details$convergence <- fit$fail
        ## obj@vcov <- solve(fit$hessian)
    }
    obj
}
temp50 <- filter(temp4, age5==50)
## debug(cureModel)

system.time(fit <- cureModel(Surv(time,event)~1, data=temp50, par=c(shape=0.127633659438804, scale=0.730222998678452, cure=-2.0680040438727)))
system.time(fit <- cureModel(Surv(time,event)~1, data=temp50))

summary(fit2 <- cureModel(Surv(time,event)~total_cat, data=temp50))
## Log-Likelihood: -113312
summary(fit2 <- cureModel(Surv(time,event)~total_cat, data=temp50, par=c(shape=0.147551370026859, scale=c(0.884863166386332, -0.115910931108631, -1.09037223824094, -1.52389713434596), cure=-2.18162256296465)))

## Simple models:
temp30plus <- filter(temp4, age5>=30) %>% group_by(age5, total_cat)
models <- do(temp30plus,cureModel(Surv(time,event)~1, data=.) %>% coef %>% Call(data.frame))
models2 <- mutate(models, cure=expit(cure), shape=exp(shape), scale=exp(scale..Intercept.),
                  scale..Intercept.=NULL, psa=switch(unclass(total_cat),"1"=0,"2"=1,"3"=3,"4"=10)) %>%
    ungroup() %>% mutate(total_cat=psa,psa=NULL)
dput(as.data.frame(models2))

## Weibull regression formulated in terms of the *weibull functions (cf. survreg)
WeibullRegression <- function(formula,data=NULL,formula.shape=~1,par=NULL,...) {
    formula.scale <- formula
    rstpm2:::lhs(formula.scale) <- NULL
    mf <- model.frame(formula,data)
    y <- model.extract(mf,"response")
    if (ncol(y)==2) { # right censored
        delayed <- FALSE
        time <- y[,1]
        status <- y[,2]
    }
    if (ncol(y)==3) { # left truncated and right censored
        entry <- y[,1]
        time <- y[,2]
        status <- y[,3]
        if (delayed <- any(entry>0))
            i <- which(entry>0)
    }
    formulae <- list(shape=formula.shape,
                     scale=formula.scale)
    X <- lapply(formulae,model.matrix,data)
    stopifnot(nrow(X$scale) == nrow(X$shape))
    cumn <- cumsum(sapply(X,ncol))
    index <- mapply(function(a,b) (a+1):b, c(0,cumn[-length(cumn)]),cumn,SIMPLIFY=FALSE)
    names(index) <- names(formulae)
    ll <-function(beta) {
        lp <- mapply(function(x,idx) as.vector(x %*% beta[idx]), X, index, SIMPLIFY=FALSE)
        shape <- exp(lp$shape)
        scale <- exp(lp$scale)
        S <- pweibull(time,shape,scale,lower.tail=FALSE)
        f <- dweibull(time,shape,scale)
        value <- sum(ifelse(status==1, log(f), log(S)))
        if (delayed) {
            S0 <- pweibull(entry[i],shape[i],scale[i],lower.tail=FALSE)
            value <- value - sum(log(S0))
        }
        value
    }
    if (is.null(par))
        par <- c(shape=rep(0,n.shape), # shape
                       scale=c(log(sum(status)/sum(time)),rep(0,n.scale-1)))
    maxLik(ll, start = par, ...)
}
system.time(fit2 <- WeibullRegression(Surv(time,event)~total_cat, data=temp50, par=c(shape=0.147551370026859, scale=c(0.884863166386332, -0.115910931108631, -1.09037223824094, -1.52389713434596))))
summary(fit2)
summary(fit.init2 <- survreg(Surv(time,event)~total_cat,data=temp50))

## AFT -> Weibull
sigma <- fit.init2$scale # AFT "scale"
shape <- 1/sigma # Weibull "shape"
logscale <- coef(fit.init2)[1]
logbeta <- (coef(fit.init2)[-1])
c(log(shape),logscale,logbeta)
summary(fit2)

## maxLik output
solve(vcov(fit))
exp(coef(fit)[1:2])
expit(coef(fit)[3])

## > solve(vcov(fit$maxLik))
##           shape     scale      cure
## shape 81839.971  4831.236 -4540.198
## scale  4831.236 56432.327  7144.990
## cure  -4540.198  7144.990  4205.504


## plot of the distributions
x <- seq(0,20,length=101)
rescreening2 <- lapply(1:nrow(rescreening),
                       function(i) data.frame(age5=rescreening$age5[i],
                                              total_cat=rescreening$total_cat[i],
                                              x=x,
                                              p=(1-rescreening$cure[i])*pweibull(x,rescreening$shape[i],rescreening$scale[i])))
rescreening2 <- do.call("rbind",rescreening2)
xyplot(p ~ x | factor(total_cat)+factor(age5),data=rescreening2,type="l",subset=age5>=55)
xyplot(p ~ x | factor(age5),data=rescreening2,group=factor(total_cat),type="l",subset=age5>=55,auto.key=TRUE)


fit2 <- cureModel(Surv(time,event)~total_cat, data=temp50, par=c(shape=0.147551370026859, scale=c(0.884863166386332, -0.115910931108631, -1.09037223824094, -1.52389713434596), cure=-2.18162256296465))
fit3 <- cureModel(Surv(time,event)~total_cat,
                  formula.shape=~total_cat,
                  data=temp50,
                  par=c(shape=c(0.27,-0.1,-0.4,-0.54),
                      scale=c(0.85,-0.1,-1.2,-1.8),
                      cure=-2.05)) # best fit?
fit4 <- cureModel(Surv(time,event)~total_cat,
                  formula.cure=~total_cat,
                  data=temp50,
                  par=c(0.154696663520933, 0.865091437590591, -0.111246732543863, -1.05240102943909,
-1.48716715943421, -2.04615245632196, -0.0207990798623829, -0.639229272169032,
-1.17752470254301))
fit5 <- cureModel(Surv(time,event)~total_cat,
                  formula.shape=~total_cat,
                  formula.cure=~total_cat,
                  data=temp50,
                  par=c(shape=c(0.27,-0.1,-0.4,-0.54),
                      scale=c(0.85,-0.1,-1.2,-1.8),
                      cure=c(-2.05,0,0,0)))

sapply(list(fit,fit2,fit3,fit4,fit5),AIC)


## check the screening results from the simulation
require(microsimulation)
require(sqldf)
require(mgcv)
set.seed(12345)
sim <- callFhcrc(1e6, screen="screenUptake")
events <- sim$summary$events
pt <- sim$summary$pt
rates <- sqldf('with pop as (select age, year, sum(pt) as pt from pt where year between 1980 and 2020 and age between 40 and 89 group by age, year), screen as (select age, year, count(*) as n from events where event in ("toScreen","toBiopsyFollowUpScreen") group by age, year) select pop.age, pop.year, pop.pt, coalesce(n,0) as n, 1.0*coalesce(n,0)/pt as rate from pop left natural join screen')
rates$pred <- predict(gam(n~s(age,year)+offset(log(pt)),data=rates,family="poisson"),newdata=transform(rates,pt=1),type="response")
ratestab <- xtabs(rate~age+year,rates)
## ratestab[ratestab==0] <- NA
filled.contour(as.numeric(dimnames(ratestab)$age),
               as.numeric(dimnames(ratestab)$year),
               ratestab)
ratestab <- xtabs(n~age+year,rates)
filled.contour(as.numeric(dimnames(ratestab)$age),
               as.numeric(dimnames(ratestab)$year),
               ratestab)

## rllogis
require(Rcpp)
require(inline)
src <- '
#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
  double rllogis(double shape, double scale) {
    double u = R::runif(0.0,1.0);
    return scale*exp(-log(1.0/u-1.0)/shape);
  }'
sourceCpp(code=src)
tmp <- sapply(1:1e4,function(i) rllogis(3.3,1.0))
plot(density(tmp,from=0))
require(eha)
tmp2 <- eha::rllogis(1e5,3.3,1.0)
lines(density(tmp2,from=0),lty=2)



## rllogis_trunc
require(Rcpp)
require(inline)
src <- '
#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
  double rllogis_trunc(double shape, double scale, double left) {
    double S0 = 1.0/(1.0+exp(log(left/scale)*shape));
    double u = R::runif(0.0,1.0);
    return scale*exp(log(1.0/(u*S0)-1.0)/shape);
  }'
sourceCpp(code=src)
tmp <- sapply(1:1e5,function(i) rllogis_trunc(3.3,1.0,3))
plot(density(tmp,from=3))
require(eha)
tmp2 <- eha::rllogis(1e6,3.3,1.0)
tmp2 <- tmp2[tmp2>=3]
lines(density(tmp2,from=3),lty=2)



## Screening uptake
require(Rcpp)
require(inline)
src <- '
#include <Rcpp.h>
using namespace Rcpp;
  double rllogis(double shape, double scale) {
    double u = R::runif(0.0,1.0);
    return scale*exp(-log(1.0/u-1.0)/shape);
  }
  double rllogis_trunc(double shape, double scale, double left) {
    double S0 = 1.0/(1.0+exp(log(left/scale)*shape));
    double u = R::runif(0.0,1.0);
    return scale*exp(log(1.0/(u*S0)-1.0)/shape);
  }
// [[Rcpp::export]]
NumericVector testRcpp(NumericVector cohort) {
      //double pscreening = cohort[0]>=1932.0 ? 0.9 : 0.9-(1932.0 - cohort[0])*0.01;
      double pscreening = 0.9;
      double shapeA = 3.8;
      double scaleA = 15.0;
      double shapeT = 2.16;
      double scaleT = 11.7;
      double uscreening = R::runif(0.0,1.0);
      double first_screen;
      if (cohort[0] > 1960.0) {
	first_screen = 35.0 + rllogis(shapeA,scaleA); // (i) age
      } else if (cohort[0] < 1945.0) {
	first_screen = (1995.0 - cohort[0]) + rllogis(shapeT,scaleT); // (ii) period
      } else {
	double age0 = 1995.0 - cohort[0];
	double u = R::runif(0.0,1.0);
	if ((age0 - 35.0)/15.0 < u) // (iii) mixture
	  first_screen = age0 + rllogis_trunc(shapeA,scaleA,age0-35.0);
	else first_screen = age0 + rllogis(shapeT,scaleT);
      }
      if (uscreening<pscreening)
	return wrap(first_screen);
      else return wrap(-1.0);
}
'
sourceCpp(code=src)
tmp <- sapply(rep(1950,1000),testRcpp)
mean(tmp==-1) # ~10%
i <- tmp != -1
plot(density(tmp[i]),xlim=c(30,120))

## Re-screening
require(Rcpp)
require(inline)
src <- '
#include <Rcpp.h>
#include "/home/marcle/src/R/microsimulation/src/rcpp_table.h"
using namespace Rcpp;
  template<class T>
  T bounds(T x, T a, T b) {
    return (x<a)?a:((x>b)?b:x);
  }
// [[Rcpp::export]]
DataFrame testRcpp(DataFrame rescreening, double time, double psa) {
  typedef Table<pair<double,double>,double> TableDDD; // as per TableBiopsyCompliance
  TableDDD rescreen_shape, rescreen_scale, rescreen_cure;
  rescreen_shape = TableDDD(rescreening, "age5", "total", "shape");
  rescreen_scale = TableDDD(rescreening, "age5", "total", "scale");
  rescreen_cure  = TableDDD(rescreening, "age5", "total", "cure");
  TableDDD::key_type key = TableDDD::key_type(bounds<double>(time,30.0,90.0),psa);
  double prescreened = 1.0 - rescreen_cure(key);
  double shape = rescreen_shape(key);
  double scale = rescreen_scale(key);
  double u = R::runif(0.0,1.0);
  double t = time + R::rweibull(shape,scale);
  bool brescreened = u<prescreened;
  return wrap(Rcpp::DataFrame::create(_["u"]=u,_["prescreened"]=prescreened,_["t"]=t,_["shape"]=shape,_["scale"]=scale,_["brescreened"]=brescreened));
}
'
sourceCpp(code=src)
rescreening$total <- rescreening$total_cat
tmp <- do.call("rbind",lapply(1:1e3,function(i) testRcpp(rescreening,71,0.1)))
head(tmp)
mean(tmp$brescreened)
plot(density(tmp$t[tmp$brescreened],from=71))

subset(rescreening,age5==60)
subset(rescreening,age5==30)

require(Rcpp)
require(RcppArmadillo)
require(inline)
src <- '
#include <RcppArmadillo.h>
using namespace Rcpp;
using namespace arma;
// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::export]]
NumericVector testRcpp(NumericVector inX) {
vec X = as<vec>(inX);
//X.resize(X.size()-1);
//return wrap(X==math::inf() || X==-math::inf());
return wrap(X);
}
'
sourceCpp(code=src)
testRcpp(as.double(c(-Inf,1,Inf,NA)))

require(Rcpp)
require(RcppArmadillo)
require(inline)
src <- '
#include <RcppArmadillo.h>
using namespace Rcpp;
using namespace arma;
// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::export]]
NumericVector testRcpp(NumericVector inX) {
rowvec X = as<rowvec>(inX);
return wrap(exp(X));
}
'
sourceCpp(code=src)
testRcpp(as.double(c(-Inf,1,Inf,NA)))


require(Rcpp)
require(RcppArmadillo)
require(inline)
src <- '
#include <RcppArmadillo.h>
using namespace Rcpp;
using namespace arma;
// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::export]]
NumericMatrix testRcpp(NumericMatrix inX) {
mat X = as<mat>(inX);
return wrap(expmat(X));
}
'
sourceCpp(code=src)
testRcpp(rbind(c(-0.1,0.1),c(0.2,-0.2)))



require(Rcpp)
require(inline)
src <- '
#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
List testRcpp() {
List list;
list["a"].push_back(1); // FAILS
return wrap(list);
}
'
sourceCpp(code=src)
testRcpp()

require(Rcpp)
require(inline)
src <- '
#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
SEXP myeval(SEXP sexp) {
Environment env = Environment::global_env();
switch(TYPEOF(sexp)) {
 case REALSXP:
 case INTSXP:
 case LGLSXP:
  return sexp;
  break;
 case SYMSXP:
  return Rf_eval(sexp,env);
  break;
 case LANGSXP:
  //
  break;
};
return wrap(0);
}
'
sourceCpp(code=src)
x=3:4
myeval(as.name("x"))
myeval(1:10)
myeval(T)

## We can readily evaluate values and expressions using R
require(Rcpp)
require(inline)
src <- '#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
SEXP myeval(SEXP sexp) {
  Environment env = Environment::global_env();
  return Rf_eval(sexp,env);
}'
sourceCpp(code=src)
x=3:4
myeval(as.name("x"))
myeval(1:10)
myeval(T)
myeval(exp(x))



require(Rcpp)
require(inline)
src <- '
#include <Rcpp.h>
// [[Rcpp::export]]
  double rnormPos1(double mean, double sd) {
    double x;
    while ((x=R::rnorm(mean,sd))<0.0) { }
    return x;
  }
'
sourceCpp(code=src)
y <- sapply(1:1e4,function(i) rnormPos1(1,2))
mean(y)
sd(y)
rnormPos <- function(n,mean=0,sd=1,lbound=0) {
    if (length(mean)<n) mean <- rep(mean,length=n)
    if (length(sd)<n) sd <- rep(sd,length=n)
    x <- rnorm(n,mean,sd)
    while(any(i <- which(x<lbound)))
        x[i] <- rnorm(length(i),mean[i],sd[i])
    x
}
mean(rnormPos(1e5,1,2))
sd(rnormPos(1e5,1,2))
trunc_mean <- function(mu,sigma) {
    alpha <- -mu/sigma
    lambda <- dnorm(alpha)/pnorm(alpha,lower.tail=FALSE)
    mu + sigma*lambda
}
trunc_var <- function(mu,sigma) {
    alpha <- -mu/sigma
    lambda <- dnorm(alpha)/pnorm(alpha,lower.tail=FALSE)
    delta <- lambda*(lambda-alpha)
    sigma^2*(1-delta)
}
trunc_mean(1,2)
sqrt(trunc_var(1,2))

require(Rcpp)
src <- "#include <Rcpp.h>
using namespace Rcpp ;
// [[Rcpp::export]]
SEXP fn_impl( Language call, Environment env){
    for (int i=0; i<1e5; ++i)
       Rf_eval( call, env ) ;
    return Rf_eval( call, env ) ;
}"
sourceCpp(code=src)
fn <- function(call, list){
    fn_impl( call, list2env(list) )
}
fn(substitute(a*2), list(a = 1)) # fast
fn(substitute(exp(x)), list(x = 1)) # fast
## system.time(fn(substitute(lm(y~x)), list(x = 1:10, y=1:10))) # slow
system.time(for (i in 1:1e5) { x = 1:10; y=1:10; x+y })
system.time(fn(substitute(x+y), list(x = 1:10, y=1:10))) # 3-4 times faster

##
dots <- function(...) {
    eval(substitute(alist(...)))
}
##
let <- function(...) {
    args <- as.list(sys.call())[-1]
    n <- length(args)
    eval(args[[n]], args[-n])
}
let(a=1,b=2,a+b)
## let(a=1,b=2,let(c=3,a+b+c))
##
let <- function(..., expr) {
    expr <- substitute(expr)
    dots <- list(...)
    eval(expr, dots)
}
let(a=1,b=2,expr=a+b)
## let(a=1,b=2,expr=let(c=3,expr=a+b+c))
##
let <- function(args, body) {
    f <- function() {}
    body(f) <- substitute(body)
    formals(f) <- args
    f()
}
let(list(a=1,b=2),a+b)
let(list(a=1,b=2),let(list(c=1),a+b+c)) # FAILS

let <- function(args, body) {
    bquote(function(a=1,b=3) {.(body)})
}
let(list(a=1,b=2),a+b)
let(list(a=1,b=2),let(list(c=1),a+b+c))()()


require(Rcpp)
src <- "#include <Rcpp.h>
using namespace Rcpp ;
// [[Rcpp::export]]
SEXP fn_impl( List args){
    Language call = args[0] ;
    List data = args[1] ;
    // now construct the call to eval:
    Language eval_call( \"eval\", call, data ) ;
    // and evaluate it
    return Rf_eval( eval_call, R_GlobalEnv ) ;
}"
sourceCpp(code=src)
fn_impl(list(substitute(a*2),list(a=1))) # FAILS

## stsplit/splitLexis using Rcpp
require(Rcpp)
require(inline)
src <- '
#include <Rcpp.h>
#include <vector>
#include <algorithm>
using namespace Rcpp;
using namespace std;
inline double min(double x, double y) { return (x<y ? x : y); }
inline double max(double x, double y) { return (x<y ? y : x); }
// [[Rcpp::export]]
DataFrame splitLexisCpp(NumericVector splits, NumericVector starts, NumericVector finishes) {
  vector<double> split(splits.begin(), splits.end());
  vector<int> _row, _event, _lower;
  vector<double> _t0, _t;
  int j;
  for (int i=0; i<starts.size(); ++i) {
    j = lower_bound(split.begin(), split.end(), starts[i]) - split.begin();
    if (starts[i] < split[j]) --j;
    for(; split[j] < finishes[i]; ++j) {
      _row.push_back(i+1);
      _lower.push_back(split[j]);
      _t0.push_back(max(split[j],starts[i]));
      _t.push_back(min(split[j+1],finishes[i]));
      _event.push_back(split[j]<=finishes[i] && finishes[i] < split[j+1] ? 1 : 0); // indicator for end of segment
    }
  }
  return DataFrame::create(_(".row")=wrap(_row),
                           _(".lower")=wrap(_lower),
                           _(".t0")=wrap(_t0),
                           _(".t")=wrap(_t),
                           _(".event")=wrap(_event));
}
'
sourceCpp(code=src)
splitLexis2 <- function(df,start,finish,splits) {
    start <- deparse(substitute(start))
    finish <- deparse(substitute(finish))
    stopifnot(all(df[[start]]<=df[[finish]]))
    stopifnot(max(df[[finish]])<=max(splits))
    stopifnot(min(splits)<=min(df[[start]]))
    index <- splitLexisCpp(splits,df[[start]],df[[finish]])
    cbind(df[index$.row,],index[,-1])
}
## system.time(out <- temp4 %>%
##     splitLexis2(age_start,age_finish,c(0,seq(40,90,by=5),120)) %>%
##     group_by(total_cat,.lower) %>%
##     summarise(pt = sum(.t-.t0), event=sum(event*.event)) %>%
##     mutate(rate = event/pt, age = .lower))
## NOT NEEDED: split by age group with time since last test as the primary time scale
## temp5 <- temp4 %>%
##     splitLexis2(age_start,age_finish,c(0,seq(40,90,by=5),120)) %>%
##     mutate(event = event * .event,
##            age = .lower,
##            .t0 = .t0 - age_start,
##            .t = .t - age_start)

## old versions
cureModel <- function(formula,data=NULL,formula.shape=~1,formula.cure=~1,par=NULL,...) {
    formula.scale <- formula
    rstpm2:::lhs(formula.scale) <- NULL
    y <- model.extract(model.frame(formula,data),"response")
    if (ncol(y)==2) { # right censored
        delayed <- FALSE
        time <- y[,1]
        status <- y[,2]
    }
    if (ncol(y)==3) { # left truncated and right censored
        entry <- y[,1]
        time <- y[,2]
        status <- y[,3]
        if (delayed <- any(entry>0))
            i <- which(entry>0)
    }
    X.shape <- model.matrix(formula.shape,data)
    X.scale <- model.matrix(formula.scale,data)
    X.cure <- model.matrix(formula.cure,data)
    n.shape <- ncol(X.shape)
    n.scale <- ncol(X.scale)
    n.cure <- ncol(X.cure)
    negll <-function(beta) {
        shape <- exp(as.vector(X.shape %*% beta[1:n.shape]))
        scale <- exp(as.vector(X.scale %*% beta[(1:n.scale)+n.shape]))
        cure <- rstpm2:::expit(as.vector(X.cure %*% beta[(1:n.cure)+n.scale+n.shape]))
        S <- cure+(1-cure)*pweibull(time,shape,scale,lower.tail=FALSE)
        f <- (1-cure)*dweibull(time,shape,scale)
        nll <- - sum(ifelse(status==1, log(f), log(S)))
        if (delayed) {
            S0 <- cure[i]+(1-cure[i])*pweibull(entry[i],shape[i],scale[i],lower.tail=FALSE)
            nll <- nll - sum(-log(S0))
        }
        nll
    }
    if (is.null(par))
        par <- c(shape=rep(0,n.shape), # shape
                       scale=c(log(sum(status)/sum(time)),rep(0,n.scale-1)), # scale
                       cure=c(logit((1-mean(status))/2), rep(0,n.cure-1)))
    fit <- optim(par, negll, hessian=TRUE, ...)
    fit$se.fit <- sqrt(diag(solve(fit$hessian)))
    fit
}
cureModel <- function(formula,data=NULL,formula.shape=~1,formula.cure=~1,par=NULL,...) {
    formula.scale <- formula
    rstpm2:::lhs(formula.scale) <- NULL
    mf <- model.frame(formula,data)
    y <- model.extract(mf,"response")
    if (ncol(y)==2) { # right censored
        delayed <- FALSE
        time <- y[,1]
        status <- y[,2]
    }
    if (ncol(y)==3) { # left truncated and right censored
        entry <- y[,1]
        time <- y[,2]
        status <- y[,3]
        if (delayed <- any(entry>0))
            i <- which(entry>0)
    }
    formulae <- list(shape=formula.shape,
                     scale=formula.scale,
                     cure=formula.cure)
    X <- lapply(formulae,model.matrix,data)
    stopifnot(nrow(X$scale) == nrow(X$shape))
    cumn <- cumsum(sapply(X,ncol))
    index <- mapply(function(a,b) (a+1):b, c(0,cumn[-length(cumn)]),cumn,SIMPLIFY=FALSE)
    names(index) <- names(formulae)
    ll <-function(beta) {
        lp <- mapply(function(x,idx) as.vector(x %*% beta[idx]), X, index, SIMPLIFY=FALSE)
        shape <- exp(lp$shape)
        scale <- exp(lp$scale)
        cure <- rstpm2:::expit(lp$cure)
        S <- cure+(1-cure)*pweibull(time,shape,scale,lower.tail=FALSE)
        f <- (1-cure)*dweibull(time,shape,scale)
        value <- sum(ifelse(status==1, log(f), log(S)))
        if (delayed) {
            S0 <- cure[i]+(1-cure[i])*pweibull(entry[i],shape[i],scale[i],lower.tail=FALSE)
            value <- value - sum(log(S0))
        }
        value
    }
    if (is.null(par))
        par <- c(shape=rep(0,n.shape), # shape
                       scale=c(log(sum(status)/sum(time)),rep(0,n.scale-1)), # scale
                       cure=c(logit((1-mean(status))/2), rep(0,n.cure-1)))
    fit <- maxLik(ll, start = par, ...)
    fit
}

require(Rcpp)
require(inline)
src <- '
#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
using namespace Rcpp;
using namespace arma;
// [[Rcpp::export]]
double logLikeCureModel(NumericVector nvtime, NumericVector nvstatus, NumericMatrix nmXshape,
NumericMatrix nmXscale, NumericMatrix nmXcure, NumericVector nvbeta) {
  double ll = 0.0;
  mat Xshape = as<mat>(wrap(nmXshape));
  mat Xscale = as<mat>(wrap(nmXscale));
  mat Xcure = as<mat>(wrap(nmXcure));
  vec time = as<vec>(wrap(nvtime));
  vec status = as<vec>(wrap(nvstatus));
  vec beta = as<vec>(wrap(nvbeta));
  int n0=Xshape.n_cols;
  int n1=n0+Xscale.n_cols;
  int n2=n1+Xcure.n_cols;
  vec shape = exp(Xshape * beta(span(0,n0-1)));
  vec scale = exp(Xscale * beta(span(n0,n1-1)));
  vec cure = 1.0/(1.0+exp(-Xcure * beta(span(n1,n2-1))));
  for (int i=0; i<nvtime.size(); ++i) {
    ll += status(i)==1.0 ?
      log(1.0-cure(i)) + R::dweibull(time(i),shape(i),scale(i),1) :
      log(cure(i)+(1.0-cure(i)) * R::pweibull(time(i),shape(i),scale(i),0,0));
  }
  return ll;
}
'
sourceCpp(code=src)
cureModel <- function(formula,data=NULL,formula.shape=~1,formula.cure=~1,par=NULL,...) {
    formula.scale <- formula
    rstpm2:::lhs(formula.scale) <- NULL
    mf <- model.frame(formula,data)
    y <- model.extract(mf,"response")
    if (ncol(y)==2) { # right censored
        delayed <- FALSE
        time <- y[,1]
        status <- y[,2]
    }
    if (ncol(y)==3) { # left truncated and right censored
        entry <- y[,1]
        time <- y[,2]
        status <- y[,3]
        if (delayed <- any(entry>0))
            i <- which(entry>0)
    }
    formulae <- list(shape=formula.shape,
                     scale=formula.scale,
                     cure=formula.cure)
    X <- lapply(formulae,model.matrix,data)
    stopifnot(nrow(X$scale) == nrow(X$shape))
    cumn <- cumsum(sapply(X,ncol))
    index <- mapply(function(a,b) (a+1):b, c(0,cumn[-length(cumn)]),cumn,SIMPLIFY=FALSE)
    names(index) <- names(formulae)
    ll <- llNew <-function(beta) {
        logLikeCureModel(time, status, X$shape,
                         X$scale, X$cure, beta)
    }
    llOld <-function(beta) {
        lp <- mapply(function(x,idx) as.vector(x %*% beta[idx]), X, index, SIMPLIFY=FALSE)
        shape <- exp(lp$shape)
        scale <- exp(lp$scale)
        cure <- rstpm2:::expit(lp$cure)
        S <- cure+(1-cure)*pweibull(time,shape,scale,lower.tail=FALSE)
        f <- (1-cure)*dweibull(time,shape,scale)
        value <- sum(ifelse(status==1, log(f), log(S)))
        if (delayed) {
            S0 <- cure[i]+(1-cure[i])*pweibull(entry[i],shape[i],scale[i],lower.tail=FALSE)
            value <- value - sum(log(S0))
        }
        value
    }
    if (is.null(par))
        par <- c(shape=rep(0,n.shape), # shape
                       scale=c(log(sum(status)/sum(time)),rep(0,n.scale-1)), # scale
                       cure=c(logit((1-mean(status))/2), rep(0,n.cure-1)))
    fit <- maxLik(ll, start = par, ...)
    fit
}


## "screening prevalence" ~= (overall screening rate)/(rescreening rate)
## = (number of screens)/pop / ((number of re-screens) / (person-time re-screened))
out <- temp4 %>%
    filter(format(start,"%Y")>="2008") %>%
    splitLexis2(age_start,age_finish,c(0,seq(40,90,by=5),120)) %>%
    group_by(.lower) %>%
    summarise(pt = sum(.t-.t0), event=sum(event*.event)) %>%
    mutate(overall_rate = event/pt, age5 = .lower)
temp2c <- temp2b %>%
    filter(year>=2008 & year<2012) %>%
    group_by(age5) %>%
    summarise(pt=sum(pt), n_tests=sum(n_tests)) %>%
    mutate(rescreen_rate=n_tests/pt)
temp2c %>%
    select(c(age5, rescreen_rate)) %>%
    inner_join(out %>% select(c(age5,overall_rate))) %>%
    mutate(prev = overall_rate / rescreen_rate)
## The rescreening rate is biased, as there is limited follow-up and declining rates for increasing duration since the previous test.

## cf. SQL
agesplits <- mutate(data.frame(start=c(0,seq(40,90,by=5))),finish=c(start[-1],120))
require(sqldf)
system.time(agg1 <- sqldf("select t1.total_cat, t2.start, t2.finish, sum(min(t1.age_finish,t2.finish) - max(t1.age_start,t2.start)) as pt, sum(event*(t1.age_finish<=t2.finish)*(t1.age_start<=t2.finish)) as event from temp4 as t1, agesplits as t2 where min(t1.age_finish,t2.finish) > max(t1.age_start,t2.start) group by t1.total_cat, t2.start, t2.finish"))

## 2D splits
require(Rcpp)
require(inline)
src <- '
#include <Rcpp.h>
#include <vector>
#include <algorithm>
using namespace Rcpp;
using namespace std;
inline double min(double x, double y) { return (x<y ? x : y); }
inline double max(double x, double y) { return (x<y ? y : x); }
// [[Rcpp::export]]
DataFrame splitLexis2DCpp(NumericVector splits0, NumericVector starts0, NumericVector finishes0,
                          NumericVector splits1, NumericVector starts1) {
  vector<double> split0(splits0.begin(), splits0.end());
  vector<double> split1(splits1.begin(), splits1.end());
  vector<int> _row, _event, _lower0, _lower1;
  vector<double> _t0, _t;
  int j, k;
  double u0, u1, offset;
  for (int i=0; i<starts0.size(); ++i) {
    offset = starts1[i] - starts0[i];
    j = lower_bound(split0.begin(), split0.end(), starts0[i]) - split0.begin();
    if (starts0[i] < split0[j]) --j;
    for(; split0[j] < finishes0[i]; ++j) {
      u0 = max(split0[j],starts0[i]) + offset;
      u1 = min(split0[j+1],finishes0[i]) + offset;
      k = lower_bound(split1.begin(), split1.end(), u0) - split1.begin();
      if (u0 < split1[k]) --k;
      for(; split1[k] < u1; ++k) {
        _row.push_back(i+1);
        _lower0.push_back(split0[j]);
        _lower1.push_back(split1[k]);
        _t0.push_back(max(split1[k],u0) - offset);
        _t.push_back(min(split1[k+1],u1) - offset);
        _event.push_back(split1[k]<=finishes0[i]+offset && finishes0[i]+offset < split1[k+1] ? 1 : 0);
      }
    }
  }
  return DataFrame::create(_(".row")=wrap(_row),
                           _(".lower1")=wrap(_lower0),
                           _(".lower2")=wrap(_lower1),
                           _(".t0")=wrap(_t0),
                           _(".t")=wrap(_t),
                           _(".event")=wrap(_event));
}
'
sourceCpp(code=src)
splitLexis2D <- function(df,start1,finish1,splits1,start2,splits2) {
    start1 <- deparse(substitute(start1))
    start2 <- deparse(substitute(start2))
    finish1 <- deparse(substitute(finish1))
    ##
    stopifnot(all(df[[start1]]<=df[[finish1]]))
    stopifnot(max(df[[finish1]])<=max(splits1))
    stopifnot(min(splits1)<=min(df[[start1]]))
    stopifnot(min(splits2)<=min(df[[start2]]))
    ##
    index <- splitLexis2DCpp(splits1,df[[start1]],df[[finish1]],splits2,df[[start2]])
    cbind(df[index$.row,],index[,-1])
}
splitLexis2D(data.frame(age0=0,age1=15,year0=1990.5),age0,age1,0:20,year0,seq(1900,2100,by=5))


require(lattice)
xyplot(rate ~ age, data=filter(out,age>=40), group=total_cat, type="l",
       auto.key=TRUE)

temp3 <- temp2 %>% filter(!is.na(total) & total>0 & age5>=50 & age5<54)
fit <- coxph(Surv(time,event)~ns(log(total),df=3),data=temp3, model=TRUE)
fit %>% summary
##
totals <- exp(with(filter(temp2,total>0),seq(min(log(total)), max(log(total)), length=301)))
X <- model.matrix(fit,data=data.frame(total=totals))
Sigma <- vcov(fit)
fitted <- as.vector(X %*% coef(fit))
se.fit <- sqrt(colSums(t(X) * (Sigma %*% t(X))))
matplot(totals,fitted+cbind(0,-1.96*se.fit,1.96*se.fit),lty=c(1,2,2),col=1,type="l",xlim=c(0,10))

require(lattice)
rmNA <- function(x) x[!is.na(x)]
d <- expand.grid(y=0:1,x=c(0,max(temp2$time/365)),group=rmNA(unique(temp2$total_cat)))
xyplot(y~x|group,
       data=d,
       panel=function(x,y,subscript) {
           plot(survfit(Surv(time/365,event)~total_cat,
                        data=temp2[subscript,],
                        subset=!is.na(total) & total>0 & age5>=70 & age5<74),
                xlab="Time from previous PSA test",
                ylab="Survival",
                col=1:4)
           })

library(survival)
library(ggplot2)
library(scales)
fortify.survfit <- function(survfit.data)
    data.frame(time = survfit.data$time,
               n.risk = survfit.data$n.risk,
               n.event = survfit.data$n.event,
               n.censor = survfit.data$n.censor,
               surv = survfit.data$surv,
               std.err = survfit.data$std.err,
               upper = survfit.data$upper,
               lower = survfit.data$lower,
               strata = rep(names(survfit.data$strata), survfit.data$strata))

d.survfit <- survfit(Surv(time/365,event)~total_cat,
                     data=temp2,
                     subset=!is.na(total) & total>0 & age5>=70 & age5<74)
ggplot(data = d.survfit) +
    geom_line(aes_string(x = 'time', y = 'surv', colour = 'strata')) +
    geom_ribbon(aes_string(x = 'time', ymin = 'lower', ymax = 'upper', fill = 'strata'), alpha = 0.5) +
    scale_y_continuous(labels = scales::percent)
## fortified <- ggplot2::fortify(d.survfit)

do.callR <- function(args,what,...)
    do.call(what,args,...)

fortified <-
    do.call("rbind",
            lapply(seq(45,85,by=5), function(age)
                   cbind(age5=sprintf("%i-%i",age,age+4),
                         ggplot2::fortify(survfit(Surv(time/365,event)~total_cat,
                                                  data=temp2,
                                                  subset=!is.na(total) & total>0 & age5==age)))))
fortified <- seq(45,85,by=5) %>%
    lapply(function(age)
           cbind(age5=sprintf("%i-%i",age,age+4),
                 ggplot2::fortify(survfit(Surv(time/365,event)~total_cat,
                                          data=temp2,
                                          subset=!is.na(total) & total>0 & age5==age)))) %>%
    do.callR("rbind")

fortified <-
    filter(temp2, !is.na(total) & total>0 & age5>=40) %>%
    group_by(age5) %>%
    do(cbind(age5=sprintf("%i-%i",unique(.$age5),unique(.$age5)+4),
                  ggplot2::fortify(survfit(Surv(time/365,event)~total_cat,
                                           data=.)))) %>%
    mutate(strata = factor(strata, levels=c("total_cat=(0,1]", "total_cat=(1,3]","total_cat=(3,10]", "total_cat=(10,Inf]")))

ggplot(data = fortified) +
    facet_wrap( ~ age5) +
    geom_line(aes_string(x = 'time', y = 'surv', colour = 'strata')) +
    geom_ribbon(aes_string(x = 'time', ymin = 'lower', ymax = 'upper', fill = 'strata'), alpha = 0.5) +
    scale_y_continuous(labels = scales::percent)

ggplot(data = fortified) +
    facet_wrap( ~ strata) +
    geom_line(aes_string(x = 'time', y = 'surv', colour = 'age5')) +
    geom_ribbon(aes_string(x = 'time', ymin = 'lower', ymax = 'upper', fill = 'age5'), alpha = 0.5) +
    scale_y_continuous(labels = scales::percent)





## re-screening rates - irrespective of time from last screen
## http://stackoverflow.com/questions/23026145/dplyr-how-to-number-label-data-table-by-group-number-from-group-by
group_number = (function(){i = 0; function() i <<- i+1 })()
temp4 <- temp2 %>% mutate(age_start = as.numeric((start-dob)/365), age_finish = as.numeric((finish-dob)/365)) %>% group_by(total_cat) %>% mutate(label = group_number())
temp4 <- as.data.frame(temp4)