#Set in-control and out-of-control parameters as in paper
theta0 <- c(-2.3,-4.5,2.5)
theta1 <- c(-1.7,-2.9,2.5)
#Small helper function to compute the paired-binary likelihood
#of the length two vector yz when the true parameters are theta
dPBin <- function(yz,theta) {
exp(dbinom(yz[1],size=1,prob=plogis(theta[1]),log=TRUE) +
dbinom(yz[2],size=1,prob=plogis(theta[2]+theta[3]*yz[1]),log=TRUE))
}
#Likelihood ratio for all four possible configurations
p <- c(dPBin(c(0,0), theta=theta0), dPBin(c(0,1), theta=theta0),
dPBin(c(1,0), theta=theta0), dPBin(c(1,1), theta=theta0))
if (surveillance.options("allExamples"))
#Compute ARL using slow, non-sparse matrix operations
pairedbinCUSUM.runlength(p,w1=c(-1,37,-9,29),w2=c(-1,7),h1=70,h2=32,
h11=38,h22=17)
#Sparse computations can be considerably (!) faster
pairedbinCUSUM.runlength(p,w1=c(-1,37,-9,29),w2=c(-1,7),h1=70,h2=32,
h11=38,h22=17,sparse=TRUE)
#Use paired binary CUSUM on the De Leval et al. (1994) arterial switch
#operation data on 104 newborn babies
data("deleval")
#Switch between death and near misses
observed(deleval) <- observed(deleval)[,c(2,1)]
#Run paired-binary CUSUM without generating alarms.
pb.surv <- pairedbinCUSUM(deleval,control=list(theta0=theta0,
theta1=theta1,h1=Inf,h2=Inf,h11=Inf,h22=Inf))
plot(pb.surv, xaxis.labelFormat=NULL, ylab="CUSUM Statistic")
######################################################################
#Scale the plots so they become comparable to the plots in Steiner et
#al. (1999). To this end a small helper function is defined.
######################################################################
######################################################################
#Log LR for conditional specification of the paired model
######################################################################
LLR.pairedbin <- function(yz,theta0, theta1) {
#In control
alphay0 <- theta0[1] ; alphaz0 <- theta0[2] ; beta0 <- theta0[3]
#Out of control
alphay1 <- theta1[1] ; alphaz1 <- theta1[2] ; beta1 <- theta1[3]
#Likelihood ratios
llry <- (alphay1-alphay0)*yz[1]+log(1+exp(alphay0))-log(1+exp(alphay1))
llrz <- (alphaz1-alphaz0)*yz[2]+log(1+exp(alphaz0+beta0*yz[1]))-
log(1+exp(alphaz1+beta1*yz[1]))
return(c(llry=llry,llrz=llrz))
}
val <- expand.grid(0:1,0:1)
table <- t(apply(val,1, LLR.pairedbin, theta0=theta0, theta1=theta1))
w1 <- min(abs(table[,1]))
w2 <- min(abs(table[,2]))
S <- upperbound(pb.surv) / cbind(rep(w1,nrow(observed(pb.surv))),w2)
#Show results
opar <- par(mfcol=c(2,1))
plot(1:nrow(deleval),S[,1],type="l",main="Near Miss",xlab="Patient No.",
ylab="CUSUM Statistic")
lines(c(0,1e99), c(32,32),lty=2,col=2)
lines(c(0,1e99), c(17,17),lty=2,col=3)
plot(1:nrow(deleval),S[,2],type="l",main="Death",xlab="Patient No.",
ylab="CUSUM Statistic")
lines(c(0,1e99), c(70,70),lty=2,col=2)
lines(c(0,1e99), c(38,38),lty=2,col=3)
par(opar)
######################################################################
# Run the CUSUM with thresholds as in Steiner et al. (1999).
# After each alarm the CUSUM statistic is set to zero and
# monitoring continues from this point. Triangles indicate alarm
# in the respective CUSUM (nearmiss or death). If in both
# simultaneously then an alarm is caused by the secondary limits.
######################################################################
pb.surv2 <- pairedbinCUSUM(deleval,control=list(theta0=theta0,
theta1=theta1,h1=70*w1,h2=32*w2,h11=38*w1,h22=17*w2))
plot(pb.surv2, xaxis.labelFormat=NULL)
Run the code above in your browser using DataLab