# NOT RUN {
##############################################################
# Simulated data example.
# - Data generated using "perfect" simulation.
# - one binary predictor.
# - 250 observations in each
# combination of predictor and
# status.
##############################################################
# Functions required for simulation
findq <- function(true.cdf,target,low,upp,
epsilon=0.0000001)
{
plow <- true.cdf(low)
pupp <- true.cdf(upp)
pcenter <- true.cdf((upp+low)/2)
err <- abs(pcenter-target)
i <- 0
while(err > epsilon)
{
i <- i + 1
if(target< pcenter)
{
upp <- (upp+low)/2
pupp <- pcenter
pcenter <- true.cdf((upp+low)/2)
err <- abs(pcenter-target)
}
if(target>= pcenter)
{
low <- (upp+low)/2
plow <- pcenter
pcenter <- true.cdf((upp+low)/2)
err <- abs(pcenter-target)
}
}
return((upp+low)/2)
}
true.cdf.nond1 <- function(x)
{
pnorm(x,2.1,sqrt(0.0324))
}
true.cdf.nond2 <- function(x)
{
0.5*pnorm(x,1.85,sqrt(0.005))+
0.5*pnorm(x,2.25,sqrt(0.005))
}
true.cdf.d1 <- function(x)
{
0.5*pnorm(x,1.95,sqrt(0.005))+
0.5*pnorm(x,2.35,sqrt(0.005))
}
true.cdf.d2 <- function(x)
{
pnorm(x,2.5,sqrt(0.0324))
}
# Simulating the data
nsim <- 250
qq <- seq(1,nsim)/(nsim+1)
y.nond1 <- rep(0,nsim)
for(i in 1:nsim)
{
aa <- findq(true.cdf.nond1,qq[i],
low=-6,upp=6,epsilon=0.0000001)
y.nond1[i] <- aa
}
y.nond2 <- rep(0,nsim)
for(i in 1:nsim)
{
aa <- findq(true.cdf.nond2,qq[i],
low=-6,upp=6,epsilon=0.0000001)
y.nond2[i] <- aa
}
y.nond <- c(y.nond1,y.nond2)
trt.nond <- c(rep(0,nsim),rep(1,nsim))
y.d1 <- rep(0,nsim)
for(i in 1:nsim)
{
aa <- findq(true.cdf.d1,qq[i],
low=-6,upp=6,epsilon=0.0000001)
y.d1[i] <- aa
}
y.d2 <- rep(0,nsim)
for(i in 1:nsim)
{
aa <- findq(true.cdf.d2,qq[i],
low=-6,upp=6,epsilon=0.0000001)
y.d2[i] <- aa
}
y.d <- c(y.d1,y.d2)
trt.d <- c(rep(0,nsim),rep(1,nsim))
# Design matrices
z.d <- cbind(rep(1,2*nsim),trt.d)
colnames(z.d) <- c("(Intercept)","trt")
z.nond <- cbind(rep(1,2*nsim),trt.nond)
colnames(z.nond) <- c("(Intercept)","trt")
# design matrix for posterior predictive inference
zpred <- rbind(c(1,0),c(1,1))
# Prior information
prior <- list(a0=10,
b0=1,
nu=4,
m0=rep(0,2),
S0=diag(100,2),
psiinv=diag(1,2),
tau1=6.01,
taus1=6.01,
taus2=2.01)
# Initial state
state <- NULL
# MCMC parameters
nburn <- 5000
nsave <- 5000
nskip <- 4
ndisplay <- 500
mcmc <- list(nburn=nburn,
nsave=nsave,
nskip=nskip,
ndisplay=ndisplay)
# Fitting the model
fit1 <- LDDProc(y.d=y.d,z.d=z.d,
y.nond=y.nond,z.nond=z.nond,
zpred.d=zpred,
prior.d=prior,
prior.nond=prior,
mcmc=mcmc,
state=state,
status=TRUE,
compute.band=TRUE)
fit1
summary(fit1)
plot(fit1)
# Ploting the conditional
# ROC curve for x=c(1,0),
# along with the true curve
par(cex=1.7,mar=c(4.1, 4.1, 1, 1))
plot(fit1$rocgrid,fit1$rocp.h[1,],type="l",
lty=2,lwd=2,ylim=c(0,1),xlim=c(0,1),
xlab="False positive rate",
ylab="True positive rate")
lines(fit1$rocgrid,fit1$rocp.l[1,],lty=2,lwd=2)
lines(fit1$rocgrid,fit1$rocp.m[1,],lty=1,lwd=2)
nn <- length(fit1$rocgrid)
tt <- rep(0,nn)
for(i in 1:nn)
{
tt[i] <- findq(true.cdf.nond1,
1-fit1$rocgrid[i],
low=-6,upp=6,
epsilon=0.0000001)
}
true.roc1 <- 1.0 - true.cdf.d1(tt)
lines(fit1$rocgrid,true.roc1,
lty=1,lwd=3,col="red")
# Ploting the conditional
# ROC curve for x=c(1,1),
# along with the true curve
par(cex=1.7,mar=c(4.1, 4.1, 1, 1))
plot(fit1$rocgrid,fit1$rocp.h[2,],type="l",
lty=2,lwd=2,ylim=c(0,1),xlim=c(0,1),
xlab="False positive rate",
ylab="True positive rate")
lines(fit1$rocgrid,fit1$rocp.l[2,],lty=2,lwd=2)
lines(fit1$rocgrid,fit1$rocp.m[2,],lty=1,lwd=2)
nn <- length(fit1$rocgrid)
tt <- rep(0,nn)
for(i in 1:nn)
{
tt[i] <- findq(true.cdf.nond2,
1-fit1$rocgrid[i],
low=-6,upp=6,
epsilon=0.0000001)
}
true.roc2 <- 1.0 - true.cdf.d2(tt)
lines(fit1$rocgrid,true.roc2,lty=1,lwd=3,col="red")
# }
Run the code above in your browser using DataLab