pairedbinCUSUM     package:surveillance     R Documentation(latin1)

_P_a_i_r_e_d _b_i_n_a_r_y _C_U_S_U_M _a_n_d _i_t_s _r_u_n-_l_e_n_g_t_h _c_o_m_p_u_t_a_t_i_o_n

_D_e_s_c_r_i_p_t_i_o_n:

     CUSUM for paired binary data as described in Steiner et al.
     (1999).

_U_s_a_g_e:

     pairedbinCUSUM(stsObj, control = list(range=NULL,theta0,theta1,
                                           h1,h2,h11,h22))
     pairedbinCUSUM.runlength(p,w1,w2,h1,h2,h11,h22, sparse=FALSE)

_A_r_g_u_m_e_n_t_s:

  stsObj: Object of class 'sts' containing the paired responses for
          each of the, say n, patients. The observed slot of 'stsObj'
          is thus a n times 2 matrix.

 control: Control object as a list containing several parameters.

        '_r_a_n_g_e' Vector of indices in the observed slot to monitor.

        '_t_h_e_t_a_0' In-control parameters of the paired binary CUSUM.

        '_t_h_e_t_a_1' Out-of-control parameters of the paired binary CUSUM.

        '_h_1' Primary control limit (=threshold) of 1st CUSUM.

        '_h_2' Primary control limit (=threshold) of 2nd CUSUM.

        '_h_1_1' Secondary limit for 1st CUSUM.

        '_h_2_2' Secondary limit for 2nd CUSUM.

       p: Vector giving the probability of the four different possibile
          states, i.e. c((death=0,near-miss=0),(death=1,near-miss=0),
          (death=0,near-miss=1),(death=1,near-miss=1)).

      w1: The parameters 'w1' and 'w2' are the sample weights vectors
          for the two CUSUMs, see eqn. (2) in the paper. We have that
          'w1' is equal to deaths 

      w2: As for 'w1'

      h1: decision barrier for 1st individual cusums

      h2: decision barrier for 2nd cusums

     h11: together with 'h22' this makes up the joing decision barriers

     h22: together with 'h11' this makes up the joing decision barriers

  sparse: Boolean indicating whether to use sparse matrix computations
          from the 'Matrix' library (usually much faster!). Default:
          'FALSE'.

_D_e_t_a_i_l_s:

     For details about the method see the Steiner et al. (1999)
     reference listed below. Basically, two individual CUSUMs are run
     based on a logistic regression model. The combined CUSUM not only
     signals if one of its two individual CUSUMs signals, but also if
     the two CUSUMs simultaneously cross the secondary limits.

_V_a_l_u_e:

     An 'sts' object with 'observed', 'alarm', etc. slots trimmed to
     the 'control$range' indices.

_A_u_t_h_o_r(_s):

     S. Steiner and M. Hoehle

_R_e_f_e_r_e_n_c_e_s:

     Steiner, S. H., Cook, R. J., and Farewell, V. T. (1999),
     Monitoring paired binary surgical outcomes using cumulative sum
     charts, Statistics in Medicine, 18, pp. 69-86.

_S_e_e _A_l_s_o:

     'categoricalCUSUM'

_E_x_a_m_p_l_e_s:

     #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))

     #Compute ARL using non-sparse matrix operations
     ## Not run: 
     pairedbinCUSUM.runlength(p,w1=c(-1,37,-9,29),w2=c(-1,7),h1=70,h2=32,h11=38,h22=17)
     ## End(Not run)

     #Sparse computations don't work on all machines (e.g. the next line
     #might lead to an error. If it works this call can be considerably (!) faster
     #than the non-sparse call.
     ## Not run: 
     pairedbinCUSUM.runlength(p,w1=c(-1,37,-9,29),w2=c(-1,7),h1=70,h2=32,
                              h11=38,h22=17,sparse=TRUE)
     ## End(Not run)

     #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.years=FALSE)


     ######################################################################
     #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
     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)

     ######################################################################
     # 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 caued 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.years=FALSE)

