source("http://myweb.uiowa.edu/pbreheny/7210/f17/notes/fun.R")

# OK
set.seed(2)
n <- 500
tt <- rexp(n)
a <- runif(n)
t <- pmin(tt, 1-a)
d <- 1*(tt < 1-a)
N <- 299
x <- seq(0.1, 2, len=N)
l <- numeric(N)
for (i in 1:N) {
  l[i] <- sum(dexp(t[d==1], rate=x[i], log=TRUE)) + sum(pexp(t[d==0], rate=x[i], low=FALSE, log=TRUE))
}
plotL(x, exp(l-max(l)), xlim=c(0, 2))
abline(v=1, lwd=3, lty=2, col='gray75')

# Bias
set.seed(2)
n <- 500
tt <- rexp(n)
a <- runif(n, 0, 1*(tt < 1) + 0.5*(tt >= 1))
t <- pmin(tt, 1-a)
d <- 1*(tt < 1-a)
N <- 299
x <- seq(0.1, 2, len=N)
l <- numeric(N)
for (i in 1:N) {
  l[i] <- sum(dexp(t[d==1], rate=x[i], log=TRUE)) + sum(pexp(t[d==0], rate=x[i], low=FALSE, log=TRUE))
}
plotL(x,exp(l-max(l)), xlim=c(0, 2))
abline(v=1, lwd=3, lty=2, col='gray75')
dev.off()

# Conditional independence
set.seed(1)
n <- 500
z <- rbinom(n, 1, 0.5)
tt <- rexp(n, 2*(z==0) + 2/5*(z==1))
a <- runif(n, 0, 1*(z==1) + 0.5*(z==0))
t <- pmin(tt, 1-a)
d <- 1*(tt < 1-a)
N <- 299
x <- seq(0.1, 2, len=N)
l <- numeric(N)
for (i in 1:N) {
  l[i] <- sum(dexp(t[d==1 & z==1], rate=2/5*x[i], log=TRUE)) +
    sum(pexp(t[d==0 & z==1], rate=2/5*x[i], low=FALSE, log=TRUE)) +
    sum(dexp(t[d==1 & z==0], rate=2*x[i], log=TRUE)) +
    sum(pexp(t[d==0 & z==0], rate=2*x[i], low=FALSE, log=TRUE))
}
pdf('cond-ind.pdf', 6, 5)
par(mar=c(5,5,1,1))
plotL(x, exp(l-max(l)), xlim=c(0,2))
abline(v=1, lwd=3, lty=2, col='gray75')
dev.off()

# Verifying that lam=1 is "correct" when the individual rates are 2 and 0.4
n <- 1e5
z <- rbinom(n, 1, 0.5)
tt <- rexp(n, 2*(z==0) + 2/5*(z==1))
z <- rbinom(n, 1, 0.5)
a <- runif(n, 0, 0.5*(z==1) + 1*(z==0))
t <- pmin(tt, 1-a)
d <- 1*(tt < 1-a)
N <- 299
x <- seq(0.1, 2, len=N)
l <- numeric(N)
for (i in 1:N) {
  l[i] <- sum(dexp(t[d==1], rate=x[i], log=TRUE)) + sum(pexp(t[d==0], rate=x[i], low=FALSE, log=TRUE))
}
plotL(x, exp(l-max(l)), xlim=c(0,2))
abline(v=1, lwd=3, lty=2, col='gray75')

# Conditional independence, ignoring the covariate
set.seed(6)
n <- 500
z <- rbinom(n, 1, 0.5)
tt <- rexp(n, 2*(z==0) + 2/5*(z==1))
a <- runif(n, 0, 0.5*(z==1) + 1*(z==0))
t <- pmin(tt, 1-a)
d <- 1*(tt < 1-a)
N <- 299
x <- seq(0.1, 2, len=N)
l <- numeric(N)
for (i in 1:N) {
  l[i] <- sum(dexp(t[d==1], rate=x[i], log=TRUE)) + sum(pexp(t[d==0], rate=x[i], low=FALSE, log=TRUE))
}
pdf('ignore.pdf', 6, 5)
par(mar=c(5,5,1,1))
plotL(x, exp(l-max(l)), xlim=c(0,2))
abline(v=1, lwd=3, lty=2, col='gray75')
dev.off()

# Informative censoring
set.seed(4)
n <- 100
tt <- rexp(n)
cc <- rexp(n)
t <- pmin(tt, cc)
d <- 1*(tt < cc)
N <- 299
x <- seq(0.1, 2, len=N)
l1 <- l2 <- numeric(N)
for (i in 1:N) {
  l1[i] <- sum(dexp(t[d==1], rate=x[i], log=TRUE)) + sum(pexp(t[d==0], rate=x[i], low=FALSE, log=TRUE))
  l2[i] <- sum(dexp(t, rate=2*x[i], log=TRUE))
}
L <- cbind(exp(l1-max(l1)), exp(l2-max(l2)))
colnames(L) <- c("Uninformative", expression("Using cens information"))
pdf('uninf.pdf', 6, 5)
par(mar=c(5,5,1,1))
plotL(x, L, xlim=c(0, 2))
abline(v=1, lwd=3, lty=2, col='gray75')
dev.off()