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

# Slide 9: Likelihood is correct
set.seed(2)
n <- 500
tt <- rexp(n)
a <- runif(n)
t <- pmin(tt, 1-a)
d <- 1*(tt < 1-a)
N <- 299
theta <- seq(0.1, 2, len=N)
l <- numeric(N)
for (j in 1:N) {
  l[j] <- sum(dexp(t[d==1], rate=theta[j], log=TRUE)) + sum(pexp(t[d==0], rate=theta[j], low=FALSE, log=TRUE))
}
plotL(theta, 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
theta <- seq(0.1, 2, len=N)
l <- numeric(N)
for (j in 1:N) {
  l[j] <- sum(dexp(t[d==1], rate=theta[j], log=TRUE)) + sum(pexp(t[d==0], rate=theta[j], low=FALSE, log=TRUE))
}
plotL(theta, exp(l-max(l)), xlim=c(0, 2))
abline(v=1, lwd=3, lty=2, col='gray75')

# 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
theta <- seq(0.1, 2, len=N)
l <- numeric(N)
for (j in 1:N) {
  l[j] <- sum(dexp(t[d==1 & z==1], rate=2/5*theta[j], log=TRUE)) +
    sum(pexp(t[d==0 & z==1], rate=2/5*theta[j], low=FALSE, log=TRUE)) +
    sum(dexp(t[d==1 & z==0], rate=2*theta[j], log=TRUE)) +
    sum(pexp(t[d==0 & z==0], rate=2*theta[j], low=FALSE, log=TRUE))
}
plotL(theta, 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
x <- rbinom(n, 1, 0.5)
tt <- rexp(n, 2*(x==0) + 2/5*(x==1))
a <- runif(n, 0, 0.5*(x==1) + 1*(x==0))
t <- pmin(tt, 1-a)
d <- 1*(tt < 1-a)
N <- 299
theta <- seq(0.1, 2, len=N)
l <- numeric(N)
for (j in 1:N) {
  l[j] <- sum(dexp(t[d==1], rate=theta[j], log=TRUE)) + sum(pexp(t[d==0], rate=theta[j], low=FALSE, log=TRUE))
}
plotL(theta, exp(l-max(l)), xlim=c(0,2))
abline(v=1, lwd=3, lty=2, col='gray75')

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

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