simulation.Rmd

---
title: Validation simulation
author: Terry Therneau
output: pdf_document
---

This is a simulation based off of the results in censoring.Rmd. 
I realized that there will always be the concern that I "selected a good seed"
so as to get the results I like.  I frankly worry about it myself.
This variant starts with the same data, models, and prediction, but then
the entire data creation and fitting process is set up as a function that can
be called multiple times.   This is run separately from the validation.Rnw
file, which uses some of the figures generated by this function.

First the intitial fits
```{r, initial}
# Setup
library(survival)
library(splines)
library(survey)
library(parallel)
options(mc.cores=5)

# Initial fits
rott2 <- rotterdam
# recurrence free survival = rfs = earlier of recurrence or death
rott2$ryear <- with(rotterdam, pmin(rtime, dtime))/365.25
rott2$rfs   <- with(rotterdam, pmax(recur, death))

# These are both pretty good models
rfit1 <- coxph(Surv(ryear, rfs) ~ meno + size + grade + pmin(nodes,10), rott2)
rfit2 <- survreg(Surv(ryear, rfs) ~ meno + size + grade + pmin(nodes,10),
                 data=rott2, dist="loglogistic")

cdata <- subset(gbsg, select=c(size, grade, nodes, meno)) # covariates
cdata$size <- cut(cdata$size, c(-1, 20, 50, 1000),
                  c("<=20", "20-50", ">50"))

# Originally I made it larger, but keep it the size of GBSG
# cdata <- cdata[rep(1:nrow(cdata), ncopy),]
ndata <- nrow(cdata)
tau <- 4

# For speed, get a single coxph baseline survival
#  A plot of the baseline is very close to a line, from 0-6 years
baseline <- survfit(rfit1, se.fit=FALSE, censor=FALSE)
basefit <- with(baseline, lm(cumhaz ~ time, subset= (time < 6)))$coefficients
plot(baseline, cumhaz=TRUE, lwd=2, xlab="Years", ylab="Cum hazard")
abline(basefit, col=3, lwd=2)
```

There will be a fixed set of $\eta = X\beta$ values for each of the 
two models.
The predicted values are the same for all iterations, so create two
separate starter data sets,
new1 for the coxph model and new2 for the aft model.
Now, for simulation purposes we can just use exponential survival for the
Cox portion, all that matters is that the data be proportional hazards, and
it would be a touch faster.  This will be new1.  The results of rfit1 will be
new0, which I'll fill in later. (Looks better to use Cox, but harder to get the
simulated times exactly right.)

```{r, basedata}
new0 <- cdata
new0$eta <- predict(rfit1, newdata=cdata, type="lp")
chaz <- summary(baseline, times=tau)$cumhaz
temp <- chaz * exp(new0$eta)  # predicted cumulative hazard
new0$phat <- 1- exp(-temp)    # P(death) by time tau

new1 <- new0
chaz <- basefit[2]*tau
temp <- chaz * exp(new1$eta)  # predicted cumulative hazard
new1$phat <- 1- exp(-temp)    # P(death) by time tau

new2 <- cdata
new2$eta <- predict(rfit2, newdata=cdata, type='lp')
new2$phat <- psurvreg(tau, new2$eta, rfit2$scale, rfit2$dist)
```

Functions to create the actual data sets.  They create a random time
for each subject, add censoring from one of five distributions, create
the time and status variables, and add the cumulative hazard at the censored
time. Censoring is
  * Fixed at year 15, which is as far as the coxph baseline reaches.
  * Administrative censoring from year 3 to 7
  * Higher censoring for those with a higher P(death) by $\tau$.
  * Lower censoring for those with high P(death) by $\tau$.
  * Administrative + 10\% malicious.

Here is the distribution of $\hat p$ for the two samples, along with the
probability of censoring by 2 years for each pattern.

```{r, patterns}
d1 <- density(new1$phat); d2 <- density(new2$phat)
matplot(cbind(d1$x, d2$x), cbind(d1$y, d2$y), type='l', lwd=2, col=1:2, lty=1,
        xlab="P(death by 4)", ylab="Density")
legend(.7, 3, c("Cox model", "AFT model"), lwd=2, lty=1, col=1:2)

phat <- seq(0, 1, length=51)
matplot(phat, cbind(0, .25, 1- exp(-2*(4^phat)/12),
                            1- exp(-2* (4^(.9-phat))/15)),
        type='l', lwd=2, col=1:5, lty=1,
        xlab="P(death by 4)", ylab="P(censor before 2)")
```

To simulate data from the Cox model, we use the cumulative hazard function.
A random exponential times the invervse CH will give data from the 
distribution.  The CH for a subject with risk score $\eta$ is $\exp(\eta)$
times the baseline CH; a random exponential will intersect this curve at the
same point as and exponential with parameter $\exp(-eta)$ intersects the original
curve.  The findInterval function can do the inverse prob very fast. 

```{r, newdata}
# two functions to create data: first a PH outcome
ndata0 <- function(ndata, type=0, threshold, base=baseline, tau=4) {
    n <- nrow(ndata)
    if (missing(threshold)) temp <- rexp(n, exp(-ndata$eta))
    else                    temp <- rexp(n, exp(-pmin(ndata$eta, threshold)))
    yy <- baseline$time[findInterval(temp, baseline$cumhaz, left.open= TRUE)]
                 
    if (type == 0) censor <- rep(100, n)
    else if (type==1) censor <- runif(n, 3, 7)
    else if (type==2) censor <- .3 + rexp(n, (4^ndata$phat)/12)
    else if (type==3) censor <- .3 + rexp(n, (4^(.9- ndata$phat))/15)
    else censor <- ifelse(runif(n) > .9 & yy > .5, 
               yy- runif(n, .2, .5), .3 + runif(n, 3.5, 7))

    ndata$ryear <- pmin(yy, censor)
    ndata$rfs   <- 1* (yy <= censor)
    
    # Add the cumulative hazard at ryear, for the corresponding hazard
    index <- findInterval(ryear, baseline$time)
    ndata$expect <- c(0, baseline$cumhaz)[index]
    ndata
}

# This version does simple exponential: faster and more likely correct
ndata1 <- function(ndata, type=0, threshold, base=basefit[2], tau=4) {
    n <- nrow(ndata)
    if (missing(threshold)) yy <- rexp(n, base*exp(ndata$eta))
    else                    yy <- rexp(n, base*exp(pmin(ndata$eta, threshold)))
                 
    if (type == 0) censor <- rep(100, n)
    else if (type==1) censor <- runif(n, 3, 7)
    else if (type==2) censor <- .3 + rexp(n, (4^ndata$phat)/12)
    else if (type==3) censor <- .3 + rexp(n, (4^(.9- ndata$phat))/15)
    else censor <- ifelse(runif(n) > .9 & yy > .5, 
               yy- runif(n, .2, .5), .3 + runif(n, 3.5, 7))

    ndata$ryear <- pmin(yy, censor)
    ndata$rfs   <- 1* (yy <= censor)
    
    # Add the cumulative hazard at ryear, which is simply rate*time
    ndata$expect <- ndata$ryear * base*exp(ndata$eta)
    ndata$expect4 <- 4*base*exp(ndata$eta)  # expect at 4 years
    ndata
}


# second an AFT outcome
ndata2 <- function(ndata, type=0, threshold, fit=rfit2) {
    n <- nrow(ndata)
    if (missing(threshold)) yy <- rsurvreg(n, ndata$eta, fit$scale, 
                                           distribution = fit$dist)
    else yy <-  rsurvreg(n, pmax(ndata$eta, threshold),  fit$scale,
                                           distribution = fit$dist)
    if (type == 0) censor <- rep(100, n)
    else if (type==1) censor <- runif(n, 4, 7)
    else if (type==2) censor <- .3 + rexp(n, (4^ndata$phat)/14)
    else if (type==3) censor <- .3 + rexp(n, (4^(.9- ndata$phat))/17)
    else censor <- ifelse(runif(n) > .9 & yy > .5, 
               yy- runif(n, .2, .5), .4 + runif(n, 4, 7))

    ndata$ryear <- pmin(yy, censor)
    ndata$rfs   <- 1* (yy <= censor)
    ndata$expect <- -log(1-psurvreg(ndata$ryear, ndata$eta, fit$scale, fit$dist))
    ndata$expect4 <- -log(1-psurvreg(4, ndata$eta, fit$scale, fit$dist))
    ndata
}
```

The intention of the censoring patterns is to have about 50% censoring in the
data.  Verify this with a small number of calls.  The result is close enough.

```{r, censorcheck}
cat("ndata1\n")
test1 <- rbind(
    sapply(1:4, function(i) {temp <- ndata1(new1, i); 1-mean(temp$rfs)}),
    sapply(1:4, function(i) {temp <- ndata1(new1, i); 1-mean(temp$rfs)}),
    sapply(1:4, function(i) {temp <- ndata1(new1, i); 1-mean(temp$rfs)}),
    sapply(1:4, function(i) {temp <- ndata1(new1, i); 1-mean(temp$rfs)}),
    sapply(1:4, function(i) {temp <- ndata1(new1, i); 1-mean(temp$rfs)}))
dimnames(test1) <- list(paste("PH sample", 1:5), letters[2:5])
round(test1, 3)

cat("\nndata2\n")
test2 <-    rbind(
        sapply(1:4, function(i) {temp <- ndata2(new2, i); 1-mean(temp$rfs)}),
        sapply(1:4, function(i) {temp <- ndata2(new2, i); 1-mean(temp$rfs)}),
        sapply(1:4, function(i) {temp <- ndata2(new2, i); 1-mean(temp$rfs)}),
        sapply(1:4, function(i) {temp <- ndata2(new2, i); 1-mean(temp$rfs)}),
        sapply(1:4, function(i) {temp <- ndata2(new2, i); 1-mean(temp$rfs)}))

dimnames(test2) <- list(paste("AFT sample", 1:5), letters[2:5])
round(test2, 3)
```

Now for the four simulation functions.  Each creates a matrix of results,
one column per simulations.

```{r, sim}
nsim <- 1000
set.seed(1953)

rttrfun <- function(data, tau=4, debug=FALSE, long=FALSE) {
    data$rwt <- rttright(Surv(ryear, rfs) ~ 1, data, times=tau)
    data$y4  <- ifelse(data$ryear <=4 & data$rfs==1, 1, 0)
    data2 <- subset(data, rwt>0)
    data2$id <- 1:nrow(data2)
    if (cor(data2$eta, data2$phat) <0) data2$eta <- - data2$eta

    observed <- sum(data2$rwt * data2$y4)
    expected <- sum(data2$rwt * data2$phat)
    cfit <-concordance(y4 ~ eta, weights=rwt, data= data2) 
    data2$llp <- log(-log(1- data2$phat))
    sdata <- svydesign(ids= ~ id, weights=~rwt, variables = ~ y4 + eta + phat +
                      eta + llp, data=data2)
    fit0 <- svyglm(y4 ~ offset(llp), design=sdata, 
                   family= quasibinomial(link='cloglog'))
    fit1 <- svyglm(y4 ~ llp, design=sdata, 
                       family= quasibinomial(link='cloglog'))
    fit2 <- svyglm(y4 ~ ns(llp,3), design=sdata, 
                       family= quasibinomial(link='cloglog'))

    # Austin measures
    newphat <- predict(fit2, type='response')
    delta <- as.vector(newphat)- data2$phat
    wmean <- function(x, wt=data2$rwt)  sum(x*wt)/sum(wt)
    wtq <- function(x, p, wt=data2$rwt) {
        ii <- order(x)
        xx <- x[ii]; ww <- wt[ii]/sum(wt)
        keep <- !duplicated(xx, fromLast=TRUE)  # avoid duplicates
        approx(cumsum(ww)[keep], xx[keep], p)$y
    }        
    austin <- c(wmean(abs(delta)), wtq(abs(delta), c(.5, .9)),
                wmean(delta))
    
    # Brier score
    p0 <- mean(data2$rwt * data2$y4)   # this will equal the KM at 4
    brier1 <- with(data2, sum(rwt * (y4 - phat)^2))
    brier2 <- with(data2, sum(rwt * (y4 - p0)^2))
    if (debug) browser()
    if (long) list(O=observed, E=expected, C= cfit, fit0= fit0, fit1=fit1,
                   fit2 = fit2, brier.n= brier1, brier.d= brier2, austin=austin)
    else c(O=observed, E=expected, C= cfit$concordance, C.std = sqrt(cfit$var),
         fit0 = summary(fit0)$coefficients[1:2],
         fit1 = summary(fit1)$coefficients[2, 1:2],
         nonlin = anova(fit1, fit2)$p,
         brier.n= brier1, brier.d= brier2, austin=austin)
}

tdat <- ndata1(new1, type=1)
tdat2 <- ndata1(new1, type=1, threshold=1)
rtest <- rttrfun(tdat)

temp1 <- mclapply(1:5, function(i) sapply(1:nsim, function(x) 
                      rttrfun(ndata1(new1, type= i-1))))
temp2 <- mclapply(1:5, function(i) sapply(1:nsim, function(x) 
                      rttrfun(ndata2(new2, type= i-1))))
rsimc <- array(unlist(temp1), dim=c(length(rtest), nsim, 5))
rsima <- array(unlist(temp2), dim=c(length(rtest), nsim, 5))

martfun <- function(data, tau=4) {
    newtime <- pmin(data$ryear, tau)
    newstat <- ifelse(data$ryear <=tau, data$rfs, 0)

    exp4 <- ifelse(newtime >= tau, data$expect4, data$expect)

    expected <- sum(exp4)
    observed <- sum(newstat)
    cfit <- concordance(newtime ~ eta, reverse=TRUE, data=data, ymax=tau)

    fit0 <- glm(newstat ~ offset(log(exp4)), family=poisson, data,
                subset= (exp4 > 0))
    fit1 <- glm(newstat ~ eta + offset(log(exp4)), family=poisson, data= data,
                subset= (exp4 > 0))
    fit2 <- glm(newstat ~ nsk(eta ,3) + offset(log(exp4)),
                subset = (exp4 > 0), family=poisson, data= data)

    c(O=observed, E= sum(expected), C= cfit$concordance, C.std = sqrt(cfit$var),
         fit0= summary(fit0)$coefficients[1:2],
         fit1= summary(fit1)$coefficients[2, 1:2],
         nonlin= anova(fit1, fit2, test="Chisq")[[5]][2])
}   

mtest <- martfun(tdat)
temp <- mclapply(1:5, function(i) sapply(1:nsim, function(x)
                                         martfun(ndata1(new1, type= i-1)))) 
msimc <- array(unlist(temp), dim=c(length(mtest), nsim, 5))
temp <- mclapply(1:5, function(i) sapply(1:nsim, function(x)
                                         martfun(ndata1(new1, type= i-1)))) 
msima <- array(unlist(temp), dim=c(length(mtest), nsim, 5))

modelfun <- function(data, tau=4, long=FALSE) {
    if (cor(data$eta, data$phat)< 0) data$eta <- - data$eta
    # fit0 doesn't result in a coefficient
    # fit0 <- coxph(Surv(ryear, rfs) ~ offset(log(-log(1-phat))), data)
    llp <- log(-log(1- data$phat))
    fit1 <- coxph(Surv(ryear, rfs) ~ llp, data)
    fit2 <- coxph(Surv(ryear, rfs) ~ ns(llp,3), data)
    fit3 <- coxph(Surv(ryear, rfs) ~ eta, data)
    # phat1 <- summary(survfit(fit1, newdata=data), time=tau)$surv # slow
    temp1 <- summary(survfit(fit1), time=tau)
    phat1 <- 1 - exp(-temp1$cumhaz* exp(fit1$linear))  # much faster
    temp2 <- summary(survfit(fit2), time=tau)
    phat2 <- 1 - exp(-temp2$cumhaz * exp(fit2$linear))
    obs = sum(phat2)  # this is a fake observed
    # austin measures
    delta <- phat2 - data$phat
    austin <- c(mean(abs(delta)), quantile(abs(delta), c(.5, .9)), mean(delta))
    cfit <- concordance(fit1)

    if (long)list(O= obs, E=sum(data$phat), C= cfit, eta=data$eta,
         fit1= fit1, fit2=fit2, fit3=fit3, 
         phat= data$phat, phat1=phat1, phat2=phat2, austin= austin)
    else c(O= obs, E=sum(data$phat), C= cfit$concordance, 
           C.std = sqrt(cfit$var),
           fit1=summary(fit1)$coefficients[c(1,3)], 
           fit3=summary(fit3)$coefficients[c(1,3)],
           nonlin= anova(fit1, fit2)[[4]][2], austin = austin)
}

htest <- modelfun(tdat)
temp <- mclapply(1:5, function(i) sapply(1:nsim, function(x)
                                         modelfun(ndata1(new1, type= i-1)))) 
hsimc <- array(unlist(temp), dim=c(length(htest), nsim, 5))
temp <- mclapply(1:5, function(i) sapply(1:nsim, function(x)
                                         modelfun(ndata1(new1, type= i-1)))) 
hsima <- array(unlist(temp), dim=c(length(htest), nsim, 5))


impfun <- function(data, tau=4) {
    km <- survfit(Surv(ryear, rfs) ~ 1, data)
    p4 <- summary(km, time=tau)$surv
    censor <- (data$ryear < tau & data$rfs==0)
    ncens <- sum(censor)

    temp <- summary(km, times= data$ryear)$surv
    iphat <- double(ncens)
    iphat[order(data$ryear)] <- temp
    pp <- p4/iphat  # probability of alive at tau, given alive at censor
    data$y4 <- ifelse(censor, 1-pp, ifelse(data$ryear <= tau & data$rfs==1, 1,0))

    expected <- sum(data$phat)
    observed <- sum(data$y4)  # a psuedo observed

    if (cor(data$eta, data$phat) < 0) data$eta <- -data$eta

    fit0 <- glm(y4 ~ offset(log(-log(1-phat))), data=data,
                   family= quasibinomial(link= 'cloglog'))
    fit1 <- glm(y4 ~ eta, data=data,
                   family= quasibinomial(link= 'cloglog'))
    fit2 <- glm(y4 ~ ns(eta,3), data= data,  
                   family= quasibinomial(link= 'cloglog'))
    cfit <- concordance(fit1)

    c(O= observed,  E= expected, C= cfit$concordance,
      C.std = sqrt(cfit$var),
      fit0 = summary(fit0)$coefficients[1:2],
      fit1= summary(fit1)$coefficients[2, 1:2], 
      nonlin= anova(fit1, fit2, test="Chisq")[[5]][2])
}

itest <- impfun(tdat)
temp <- mclapply(1:5, function(i) sapply(1:nsim, function(x)
                                         impfun(ndata1(new1, type= i-1)))) 
isimc <- array(unlist(temp), dim=c(length(itest), nsim, 5))
temp <- mclapply(1:5, function(i) sapply(1:nsim, function(x)
                                         impfun(ndata1(new1, type= i-1)))) 
isima <- array(unlist(temp), dim=c(length(itest), nsim, 5))

# pseudo-values
psfun <- function(data, tau=4) {
    zz <- data   # avoid a stupid issue with formulas
    fit <-survfit(Surv(ryear, rfs) ~ 1, data=zz)
    ps <- pseudo(fit, times=tau, type="pstate")
    llp <- log(-log(1- data$phat))

    observed <- sum(ps)
    expected <- sum(data$p4)

    fit0 <- glm(ps ~ offset(llp), family= gaussian(link= bcloglog()))
    fit1 <- glm(ps ~ llp, family= gaussan(link=bcloglog()))
    fit2 <- glm(ps ~ ns(llp, 3), family= gaussian(link= bcloglog()))
}
ptest <- psfun(tdat)

                

save(rsima, rsimc, isima, isimc, hsima, hsimc, msima, msimc, 
     file="censoring2.rda")
```

Look at the non-linear part of the equations.

```{r, nonlin}
rtest1 <- rttrfun(tdat,  long=TRUE)
rtest2 <- rttrfun(tdat2, long=TRUE)

# termplot(rtest2$fit2)   # doesn't work with svyglm objects
svgplot <- function(x, type="response") {
    fit <- x$fit2
    resid <- predict(fit, type='terms', se.fit=TRUE)
    pp <- predict(fit)   # this includes the intercept
    eta <- fit$data[["llp"]]

    ii <- order(eta)
    yy <- drop(resid$fit) + outer(drop(resid$se.fit), c(0, -1.96, 1.96), '*')
    yy <- yy + mean(pp)  # add the intercept back in
    p2 <- 1- exp(-exp(pp))
    if (type=='response')
        matplot(eta[ii], yy[ii,], type='l', lty=c(1,2,2), col=1, 
                xlab="Linear predictor", ylab="Fitted")
    else plot(fit$data[["phat"]], p2, ylim=range(c(fit$data[["phat"]], p2)),
              xlab="P(death), reference model", 
              ylab= "P(death), validation, RTTR")
}

pdf("../figures/validatea1.pdf")
oldpar <- par(mfrow=c(2,2), mar=c(5,5,1,1))
svgplot(rtest1); abline(0, 1, lty=3)
text(0, -1, paste("p=", round(anova(rtest1$fit1, rtest1$fit2)$p, 2)))
svgplot(rtest2); abline(0,1, lty=3)
text(0, -1, paste("p=", round(anova(rtest2$fit1, rtest2$fit2)$p, 2)))

svgplot(rtest1, type='phat'); abline(0,1, lty=3)
text(.7, .4, paste0("error= ", round(rtest1$austin[2], 2), ", ", 
                             round(rtest1$austin[3], 2)))
svgplot(rtest2, type='phat'); abline(0,1, lty=3)
text(.7, .4, paste0("error= ", round(rtest2$austin[2], 2), ", ", 
                             round(rtest2$austin[3], 2)))
par(oldpar)
dev.off()
```

mtest1 <- modelfun(tdat, long=TRUE)
mtest2 <- modelfun(tdat2, long=TRUE)

svgplot2 <- function(x, type="response") {
	# this version for a Cox model
	fit <- x$fit2
	temp <- termplot(fit, term=1, se=TRUE, plot=FALSE)[[1]]
	yy <- temp$y + outer(temp$se, c(0, -1.96, 1.96), '*')
    if (type=='response')
        matplot(temp$x, yy, type='l', lty=c(1,2,2), col=1, 
                xlab="Linear predictor", ylab="Fitted")
    else plot(x$phat, x$phat2, ylim=range(c(x$phat, x$phat2)),
              xlab="P(death), reference model", 
              ylab= "P(death), validation, Cox")
}

oldpar <- par(mfrow=c(2,2), mar=c(5,5,1,1))
svgplot2(mtest1)
text(0, -0.5, paste("p=", round(anova(mtest1$fit1, mtest1$fit2)[2,4], 2)))
svgplot2(mtest2); abline(0,1, lty=3)
text(0, -1, paste("p=", round(anova(mtest2$fit1, mtest2$fit2)$p, 2)))

svgplot2(mtest1, type='phat'); abline(0, 1, lty=3)
abline(0,1, lty=3)
text(.7, .4, paste0("error= ", round(mtest1$austin[2], 2), ", ", 
                             round(mtest1$austin[3], 2)))
svgplot2(mtest2, type='phat'); abline(0,1, lty=3)
text(.7, .4, paste0("error= ", round(mtest2$austin[2], 2), ", ", 
                             round(mtest2$austin[3], 2)))
par(oldpar)


This next section creates a pair of draft figures for the paper.

```{r, draft}
bplot <- function(sim, parm=0, zname='x') {
    if (parm==0) {
        z <- sim[1,,]/sim[2,,]
        if (missing(zname)) zname <- "Observed/Expected"
    }
    else {
        z <- sim[parm,,]
    }
    grp <- factor(col(z),1:5, c("None", "Independent", "High","Low", "Biased"))

    boxplot(c(z) ~ grp, ylab= zname, xlab="Censoring")
}

dplot <- function(sim, parm=0, yname='x') {
     if (parm==0) {
        y = sim[1,,]/sim[2,,]
        if (missing(yname)) yname = "Observed/Expected"
     }
     else {
         y <- sim[parm,,]
     }

     dd <- apply(y, 2, density) 
     matplot(sapply(dd, function(x) x$x), sapply(dd, function(x) x$y),
             type='l', lwd=2, col=1:5, xlab=yname)
}


pdf("draft1a.pdf")
oldpar <- par(mfrow=c(2,2), mar=c(5,5,1,1))
bplot(rsimc, 0)
bplot(rsimc, 7, 'Slope')
bplot(rsima, 0)
bplot(rsima, 7, 'Slope')
par(oldpar)
dev.off()

qf1 <- function(sim) {
    o.e <- sim[1,,]/sim[2,,]
    t(apply(o.e, 2, quantile, c(.1, .25, .5, .75, .9)))
}

oe.c <- rbind(qf1(rsimc), qf1(isimc), qf1(msimc), qf1(hsimc))
oe.a <- rbind(qf1(rsima), qf1(isima), qf1(msima), qf1(hsima))

oeplot <- function(oe, xlab= "Excess deaths") {
    yy <- c(23:19, 17:13, 11:7, 5:1)
    oldpar <- par(mar=c(5, 7, 1, 1))
    matplot(oe, yy, type='n', xlab=xlab, ylab="", yaxt='n')
#    matpoints(oe[,c(1,3,5)], yy, pch=c(1,19,1), col=1)
    points(oe[,3], yy, pch=19)
    ylab<- c("RTTR:   none", "random", "a", "b", "c",
             "Impute:  none", "random", "a", "b", "c",
             "Counting:  none", "random", "a", "b", "c",
             "Model:   none", "random", "a", "b", "c")

    segments(oe[,1], yy, oe[,5], yy)
    axis(2, yy, ylab, las=2, cex.axis=.8)

    abline(v=1, lty=3)
    par(oldpar)
}
oeplot(oe.c)
oeplot(oe.a)

qf2 <- function(sim) {
    t(apply(sim, 2, quantile, c(.1, .25, .5, .75, .9)))
}
slope.c <- rbind(qf2(rsimc[7,,]), qf2(isimc[7,,]), qf2(1+ msimc[7,,]),
                 qf2(hsimc[5,,]))
slope.a <- rbind(qf2(rsima[7,,]), qf2(isima[7,,]), qf2(1+ msima[7,,]),
                 qf2(hsima[5,,]))

oeplot(slope.c)
oeplot(slope.a)

qf3 <- function(sim) {  # brier, only avail for rttr
    R2 <- 1- sim[10,,]/sim[11,,]
    t(apply(R2, 2, quantile, c(.1, .25, .5, .75, .9)))
}
round(qf3(rsimc),3)
round(qf3(rsima), 3)



# figures for the paper
pdf("../figures/validatez1.pdf")
oeplot(oe.c)
dev.off()
pdf("../figures/validatez2.pdf")
oeplot(slope.c, "Calibration slope")
dev.off()

pdf("../figures/validatez3.pdf")
oeplot(oe.a)
dev.off()

pdf("../figures/validatez4.pdf")
oeplot(slope.a, "Calibration slope")
dev.off()
```