################################################################################ # Calculate PSe and PSp for any number of stages and K=2 # # The calculations for PSe and PSp are performed by breaking up the large # summation for each of them into parts. These parts then are simply put into # long vectors and elementwise multiplication is performed. Elements in # the resulting vector are then summed to produce the desired results. # # Notes on understanding the code # 1) The individual "i" index of the paper is denoted by an "a" here instead. # 2) The object S.all is "S" in the paper. The object S is "L" in the paper. # # Author: Chris Bilder ################################################################################ # Additional functions needed for PSeAllStages() # Functions help to find the correct indices. For K >2, it may be easier to use # matrix multiplication like with MRCV research. # Could find K within functions here but more efficient to find only once in the function which calls these. # Row number in joint.p is for an infection k positive and all other infections are negatives # This is the same function as row.for.k() in Exp_tests.R - If changes are made to it here # changes may need to be made for Exp_tests.R row.for.k <- function(k, K, ordering) { row.numbers <- 1:2^K # Row numbers in ordering row.numbers[ordering[,k] == 1 & rowSums(ordering) == 1] } # Row in joint.p where disease kbar is 0's # This is the same function as all.zero.kbar() in Exp_tests.R - If changes are made to it here # changes may need to be made for Exp_tests.R all.zero.kbar <- function(k, K, ordering) { row.numbers <- 1:2^K # Row numbers in ordering row.numbers[ordering[,-k] == 0] } # Row in joint.p matrix disease k is 1 # This is the same function as all.one.k() in Exp_tests.R - If changes are made to it here # changes may need to be made for Exp_tests.R all.one.k <- function(k, K, ordering) { row.numbers <- 1:2^K # Row numbers in ordering row.numbers[ordering[,k] == 1] } # Used in testing error part # Calculates 1 - P(G_{Lj+}^{(s')} = 0 | G~_{Ljkbar}^{(s')} = g~_{Ljkbar}^{(s')}, Y~_ak = 1) # g~_{Ljkbar}^{(s') = 0 part comes before 1 part test.error <- function(Se, Sp, s, k, kbar) { 1 - (1 - Se[k,s]) * c(Sp[kbar,s], (1 - Se[kbar,s])) } # Used in transition probability part # For k = 1, prod(p_b+0) over b <> a in stage s for group of interest lambda.calc <- function(joint.p, group.mem, s, a, k, infections, ordering) { K <- length(infections) I.11 <- ncol(group.mem) # I.11 is I because this function only called for S > 2 # Included "& !is.na(group.mem[s,])" to deal with NAs for algorithms where it is not possible for everyone # to be tested in all stages index.p.part.no.a <- group.mem[s,] == group.mem[s,a] & !is.na(group.mem[s,]) & (1:I.11 != a) # All people in group except individual a index.pb.plus.AllZeros <- all.zero.kbar(k = k, K = K, ordering = ordering) # prod.pb.plus.AllZeros.SubGroup # Added as.matrix() for the case when subgroup size was only 2 leading to # joint.p being only one individual (a vector) and colSums() will not work for it prod(colSums(as.matrix(joint.p[index.pb.plus.AllZeros, index.p.part.no.a]))) } ################################################################################ # Called by PSeAllStages() to calculate PSe # Calculate PSe for K = 2 infections; There are some coding conventions used here to help generalize to K > 2 in future PSeAllStages <- function(joint.p, group.mem, a, k, infections = 1:2, Se = Se, Sp = Sp, ordering) { K <- length(infections) # Number of infections - this is here for future generalizations kbar <- infections[-k] # Infections other than k S.all <- nrow(group.mem) # Number of total stages for algorithm S <- sum(!is.na(group.mem[,a])) # Number of stages for individual a where assume group.mem is coded correctly - NAs only appear at end ############################################################# #Testing error calculations excluding the leading S_e:Ljk test.error.part <- test.error(Se = Se, Sp = Sp, s = 1, k = k, kbar = kbar) # Stage 1 if (S > 2) { for(s in 2:(S-1)) { test.error.part.newstage <- test.error(Se = Se, Sp = Sp, s = s, k = k, kbar = kbar) test.error.part <- kronecker(test.error.part, test.error.part.newstage) } } # Order for four stages: (stage 1, stage2, stage3) = (0,0,0), (0,0,1), (0,1,0), (0,1,1), # (1,0,0), (1,0,1), (1,1,0), (1,1,1) # Could also use expand.grid() and take prod() of rows ######################################################## # Calculations needed for P(G~_{LjKBAR}^{(1)} = g~_{LjKBAR}^{(1)} | Y~_ak = 1) # = P(G~_11k = g~_11k | Y~_ak = 1) # Note by including the ordering object here, this is how I can look at disease 2 when needed. # Thus, disease 2 essentially becomes disease 1. index.pa.1andAllZeros <- row.for.k(k = k, K = K, ordering = ordering) # For k = 1, p_a10 pa.1andAllZeros <- joint.p[index.pa.1andAllZeros, a] #The "AllZeros" terminology is in anticipation of generalizing for K > 2. For now, it is just one 0. index.pa.1andPlus <- all.one.k(k = k, K = K, ordering = ordering) # For k = 1, p_a1+ pa.1andPlus <- sum(joint.p[index.pa.1andPlus, a]) index.pb.plus.AllZeros <- all.zero.kbar(k = k, K = K, ordering = ordering) # For k = 1, prod(p_b+0) over b <> a prod.pb.plus.AllZeros <- prod(colSums(joint.p[index.pb.plus.AllZeros,-a])) # P(G~_11kbar = 0 | Y~_ak = 1) end.prob <- pa.1andAllZeros * prod.pb.plus.AllZeros / pa.1andPlus end.prob.vec <- rep(x = c(end.prob, (1 - end.prob)), each = 2^(S-2)) # Order for four stages: (stage 1, stage2, stage3) = (0,0,0), (0,0,1), (0,1,0), (0,1,1), # (1,0,0), (1,0,1), (1,1,0), (1,1,1) ######################################################## # Transition probability calculations # P(G~_{Ljkbar}^{(s'+1)} = g~_{LjKBAR}^{(s'+1)} | G~_{Ljkbar}^{(s')} = g~_{LjKBAR}^{(s')}, Y~_ak = 1) nu <- 1 # Need this here for S = 2 so that calculations work when execute # save.stage <- test.error.part*nu*end.prob.vec nu.trans.prob <- list() # Use a list storage format so that it can be dynamically extended to a new size depending on S # nu.trans.label <- list() # For diagnostic purposes only, not printed or saved outside of function if (S > 2) { for (s in 1:(S-2)) { # For k = 1, prod(p_b+0) over b <> a in stage s' + 1 for group of interest lambda.minus.low <- lambda.calc(joint.p = joint.p, group.mem = group.mem, s = s + 1, a = a, k = k, infections = infections, ordering = ordering) # For k = 1, prod(p_b+0) over b <> a in stage s' for group of interest lambda.minus.high <- lambda.calc(joint.p = joint.p, group.mem = group.mem, s = s, a = a, k = k, infections = infections, ordering = ordering) # nu notation: nu.(stage s').(stage s'+1) # P(G~_{Ljkbar}^{(s'+1)} = 0 | G~_{Ljkbar}^{(s')} = 0, Y~_ak = 1) nu.0.0 <- 1 # P(G~_{Ljkbar}^{(s'+1)} = 0 | G~_{Ljkbar}^{(s')} = 1, Y~_ak = 1) nu.1.0 <- pa.1andAllZeros * (lambda.minus.low - lambda.minus.high) / (pa.1andPlus - pa.1andAllZeros * lambda.minus.high) # P(G~_{Ljkbar}^{(s'+1)} = 1 | G~_{Ljkbar}^{(s')} = 0, Y~_ak = 1) nu.0.1 <- 0 # P(G~_{Ljkbar}^{(s'+1)} = 1 | G~_{Ljkbar}^{(s')} = 1, Y~_ak = 1) nu.1.1 <- (pa.1andPlus - pa.1andAllZeros * lambda.minus.low) / (pa.1andPlus - pa.1andAllZeros * lambda.minus.high) nu.trans.prob[[s]] <- c(nu.0.0, nu.0.1, nu.1.0, nu.1.1) # nu.trans.label[[s]] <- paste("Transition from", s , "to", s+1) # Diagnostic purposes } nu <- nu.trans.prob[[1]] if (S > 3) { for (s in 2:(S-2)) { nu <- rep(x = nu.trans.prob[[s]], times = 2) * rep(x = nu, each = 2) } } } # Example of nu ordering for S = 4 # g1 g2 g3 nu.g2->g3 nu.g1->g2 # 0 0 0 0->0 0->0 # 0 0 1 0->1 0->0 # 0 1 0 1->0 0->1 # 0 1 1 1->1 0->1 # 1 0 0 0->0 1->0 # 1 0 1 0->1 1->0 # 1 1 0 1->0 1->1 # 1 1 1 1->1 1->1 # In the PSe expression, the first summation is over g1, the second summation is over g2, ..., # so this controls the ordering of terms in the nu vector save.stage <- test.error.part * nu * end.prob.vec # All vectors are in the correct order for elementise multiplication Se[k,S] * sum(save.stage) } ################################################################################ # Additional functions needed for PSePSpAllStages() # Testing error calculations: 1 - P(G_{sj} = (0,0) | G~_{sj} = g~_{sj}) # Order for diseases 1 and 2 in the conditioning part # 1 2 # 0 0 # 0 1 # 1 0 # 1 1 test.error2 <- function(Se, Sp, s, k, kbar) { # The Product k = 1, ... , K is done because there are both the k and kbar diseases # in kronecker() 1 - kronecker(c(Sp[k,s], (1 - Se[k,s])), c(Sp[kbar,s], (1 - Se[kbar,s]))) } # Row in joint.p where disease k is 0's all.zero.k <- function(k, K, ordering) { row.numbers <- 1:2^K # Row numbers in ordering row.numbers[ordering[,k] == 0] } ################################################################################ # Calculate PSp and call PSeAllStages() to calculate PSe; all calculations are for K = 2 # The computational approach here is to find all possible terms in the large summations for PSp and PSe and # then simply sum these terms up. # For example, consider the case of finding P(Y_ik = 1) and S = 3. # We need to sum over g~_{3j}^{(1)}, g~_{3j}^{(2)}, and g~_3jk. There are a total of # 4 * 4 * 2 = 32 separate terms in the summation. At intermediate steps, these terms are put # in the order of # # g~_{3jk}^(1) g~_{3jkbar}^(1) g~_{3jk}^(2) g~_{3jkbar}^(2) g~_3jk # 0 0 0 0 0 # 0 0 0 0 1 # 0 0 0 1 0 # 0 0 0 1 1 # 0 0 1 0 0 # 0 0 1 0 1 # 0 0 1 1 0 # 0 0 1 1 1 # # 0 1 0 0 0 # 0 1 0 0 1 # #... # # 1 1 1 1 0 # 1 1 1 1 1 # # Once all of the terms are found, they are summed to find P(Y_ik = 1). PSePSpAllStages <- function(joint.p, group.mem, a, k, infections = 1:2, Se = Se, Sp = Sp, ordering) { K <- length(infections) # Number of infections - this is here for future generalizations kbar <- infections[-k] # Infections other than k S.all <- nrow(group.mem) # Number of total stages for algorithm S <- sum(!is.na(group.mem[,a])) # Number of stages for individual a where assume group.mem is coded correctly - NAs only appear at end I.11 <- ncol(group.mem) # I=I.11 for S > 2 and I for S = 2 ################################ # Testing error calculations for P(Y_ik = 1). This will compute all needed 1 - Prod((1 - Se)^g~ * Sp^(1-g~)) values # and put into the correct order test.error.part <- test.error2(Se = Se, Sp = Sp, s = 1, k = k, kbar = kbar) # Stage 1 # Stages 2 to S - 1 if (S > 2) { for(s in 2:(S-1)) { test.error.part.newstage <- test.error2(Se = Se, Sp = Sp, s = s, k = k, kbar = kbar) test.error.part <- kronecker(test.error.part, test.error.part.newstage) } } # Last stage test.error.part <- kronecker(test.error.part, c((1-Sp[k,S]), Se[k,S])) #Example ordering when S = 2 for G11 and g2ak # g1k g1kbar g2ak # 0 0 0 # 0 0 1 # 0 1 0 # 0 1 1 # 1 0 0 # 1 0 1 # 1 1 0 # 1 1 1 ################################### # theta_g~_11 = P(G~_{Lj}^{(1)} = g~_{Lj}^{(1)}) = P(G~_11 = g~_11) - # Note by including the ordering object here, this is how I can look at disease 2 when needed. # Thus, disease 2 essentially becomes disease 1. # k is the first infection and kbar is the second infection given in my code labels below index.all.zeros <- rowSums(ordering) == 0 theta00 <- prod(as.matrix(joint.p[index.all.zeros,])) index.all.zero.k <- all.zero.k(k = k, K = K, ordering = ordering) theta01 <- prod(colSums(as.matrix(joint.p[index.all.zero.k,]))) - theta00 index.all.zero.kbar <- all.zero.kbar(k = k, K = K, ordering = ordering) theta10 <- prod(colSums(as.matrix(joint.p[index.all.zero.kbar,]))) - theta00 theta11 <- 1 - theta01 - theta10 - theta00 theta <- rep(x = c(theta00, theta01, theta10, theta11), each = 2 * 4^(S-2)) # Reason for "each = 2 * 4^(S-2)": These four items need to be repeated each 2 * 4^(S-2) times. # In the case of S = 3, this means each is repeated 8 times. The comment before the main function shows # that each set of g1k and g1kbar is repeated 8 times. ########################################################### # Transition from stage S-1 to stage S, i.e., P(G~_Sjk = g~_Sjk | G~_{Sj}^{(S-1)} = g~_{Sj}^{(S-1)}) # Example notation for S = 2: prob01to0 means for g1k = 0, g1kbar = 1, g2ak = 0 # i.e., prob k kbar to k prob00to0 <- 1 prob00to1 <- 0 # All people in group containing individual a except individual a index.p.part.no.a <- group.mem[S-1,] == group.mem[S-1,a] & !is.na(group.mem[S-1,]) & (1:I.11 != a) # All people in group containing individual a index.p.part.a <- group.mem[S-1,] == group.mem[S-1,a] & !is.na(group.mem[S-1,]) #If k = 1, this would be Prod(p_b+0, b <> a). prod.p.bplus0.no.a <- prod(colSums(as.matrix(joint.p[index.all.zero.kbar, index.p.part.no.a]))) #If k = 1, Prod(p_b+0) prod.p.bplus0 <- prod(colSums(as.matrix(joint.p[index.all.zero.kbar, index.p.part.a]))) # Prod(p_b00) prod.pb00 <- prod(as.matrix(joint.p[index.all.zeros, index.p.part.a])) #p_a00 pa00 <- joint.p[index.all.zeros, a] # For k = 1, p_a0+ pa.0andPlus <- sum(joint.p[index.all.zero.k, a]) #If k = 1, Prod(p_b0+) prod.p.b0plus <- prod(colSums(as.matrix(joint.p[index.all.zero.k, index.p.part.a]))) # P(G~_Sjk = 0 | G~_{Sj}^{(S-1)} = (0,1)) with k = 1, kbar = 2 # If want k = 2, the code takes this into account with how the ordering object is used # to re-index what is considered k and kbar (i.e., "2" becomes "1") prob01to0 <- 1 # P(G~_Sjk = 1 | G~_{Sj}^{(S-1)} = (0,1)) with k = 1, kbar = 2 prob01to1 <- 0 # P(G~_Sjk = 0 | G~_{Sj}^{(S-1)} = (1,0)) with k = 1, kbar = 2 # k=1, S=2: P(G~_2a1 = 0| G~_11=(1,0)) prob10to0 <- (pa00 * prod.p.bplus0.no.a - prod.pb00)/(prod.p.bplus0 - prod.pb00) # P(G~_Sjk = 1 | G~_{Sj}^{(S-1)} = (1,0)) with k = 1, kbar = 2 # Could do compliment of above prob10to1 <- (prod.p.bplus0 - pa00 * prod.p.bplus0.no.a)/(prod.p.bplus0 - prod.pb00) # P(G~_Sjk = 0 | G~_{Sj}^{(S-1)} = (1,1)) with k = 1, kbar = 2 prob11to0 <- (pa.0andPlus - pa00 * prod.p.bplus0.no.a - prod.p.b0plus + prod.pb00) / (1 - prod.p.bplus0 - prod.p.b0plus + prod.pb00) # P(G~_Sjk = 1 | G~_{Sj}^{(S-1)} = (1,1)) with k = 1, kbar = 2 prob11to1 <- (1 - pa.0andPlus - prod.p.bplus0 + pa00 * prod.p.bplus0.no.a) / (1 - prod.p.bplus0 - prod.p.b0plus + prod.pb00) #Example ordering for G11 and g2ak # g1k g1kbar g2ak # 0 0 0 # 0 0 1 # 0 1 0 # 0 1 1 # 1 0 0 # 1 0 1 # 1 1 0 # 1 1 1 prob.Sminus1.k.kbar.to.Sk <- rep(x = c(prob00to0, prob00to1, prob01to0, prob01to1, prob10to0, prob10to1, prob11to0, prob11to1), times = 4^(S-2)) # "times = 4^(S-2)" corresponds to: 1) There are 2^K = 2^2 = 4 possible values for each G~_sj vector. # Because these probabilties correspond to stages S and S-1, the probabilities need to be repeated # over the first S - 2 stages. Thus, these probabilities correspond to the last 3 columns of the example # ordering given in the comments for this function. # #Example ordering for Gstage1, Gstage2, and g3ak # g1k g1kbar g2k g2kbar g3ak # 0 0 0 0 0 # 0 0 0 0 1 # 0 0 0 1 0 # 0 0 0 1 1 # 0 0 1 0 0 # 0 0 1 0 1 # 0 0 1 1 0 # 0 0 1 1 1 # 0 1 0 0 0 # 0 1 0 0 1 # ... ################################### #Transition from stage s-1 to stage s # P(G~_{Sj}^{(s'+1)} = g~_{Sj}^{(s'+1)} | G~_{Sj}^{(s')} = g~_{Sj}^{(s')}) # For example, with S = 3, we need to find P(G~_{3j}^(2) = g~_{3j}^(2) | G~_11 = g_11) # Notation for S = 3: prob00to00 means for g1k = 0, g1kbar = 0 to g2ak = 0, g2akbar = 0 # prob k kbar at stage s' to k kbar at stage s'+1 trans.prob <- 1 trans.prob.list <- list() if (S > 2) { for (s in 1:(S-2)) { # For K > 2, it may be better to put all probabilties in a matrix or K-dimensional array # P(G~_{Sj}^{(s'+1)} = (0,0) | G~_{Sj}^{(s')} = (0,0)) prob00to00 <- 1 # P(G~_{Sj}^{(s'+1)} = (0,1) | G~_{Sj}^{(s')} = (0,0)) prob00to01 <- 0 prob00to10 <- 0 prob00to11 <- 0 # All people in group at stage s for individual a index.stage.s <- group.mem[s,] == group.mem[s,a] & !is.na(group.mem[s,]) # All people in group at stage s+1 for individual a index.stage.splus1 <- group.mem[s+1,] == group.mem[s+1,a] & !is.na(group.mem[s+1,]) # S=3, prod(p_b00, b in B_{Sj}^{(s')} where j is group containing individual a prod.pb00.s <- prod(as.matrix(joint.p[index.all.zeros, index.stage.s])) # S=3, prod(p_b00, b in B_{Sj}^{(s'+1)}) prod.pb00.splus1 <- prod(as.matrix(joint.p[index.all.zeros, index.stage.splus1])) # S=3, prod(p_b0+, b in B_{Sj}^{(s')} prod.pb0plus.s <- prod(colSums(as.matrix(joint.p[index.all.zero.k, index.stage.s]))) # S=3, prod(p_b0+, b in B_{Sj}^{(s'+1)} prod.pb0plus.splus1 <- prod(colSums(as.matrix(joint.p[index.all.zero.k, index.stage.splus1]))) # S=3, prod(p_b+0, b in B_{Sj}^{(s')} prod.pbplus0.s <- prod(colSums(as.matrix(joint.p[index.all.zero.kbar, index.stage.s]))) # S=3, prod(p_b+0, b in B_{Sj}^{(s'+1)} prod.pbplus0.splus1 <- prod(colSums(as.matrix(joint.p[index.all.zero.kbar, index.stage.splus1]))) #Who was in group at stage s that did not make it into stage s + 1 group? group.numb <- group.mem[s+1,a] # Group number of individual a stage s + 1 group.numb.prev.stage <- group.mem[s,a] bar.group.mem <- group.mem[s,] == group.numb.prev.stage & !(group.numb == group.mem[s+1,]) & !is.na(group.mem[s+1,]) # prod(p_b0+, b in Bbar_{Sj}^{(s'+1)} prod.pb0plus.bar.splus1 <- prod(colSums(as.matrix(joint.p[index.all.zero.k, bar.group.mem]))) # prod(p_b+0, b in Bbar_{Sj}^{(s'+1)} prod.pbplus0.bar.splus1 <- prod(colSums(as.matrix(joint.p[index.all.zero.kbar, bar.group.mem]))) # P(G~_{Sj}^{(s'+1)} = (0,0) | G~_{Sj}^{(s')} = (0,1)) prob01to00 <- (prod.pb00.splus1 * prod.pb0plus.bar.splus1 - prod.pb00.s) / (prod.pb0plus.s - prod.pb00.s) # P(G~_{Sj}^{(s'+1)} = (0,1) | G~_{Sj}^{(s')} = (0,1)) prob01to01 <- (prod.pb0plus.s - prod.pb00.splus1 * prod.pb0plus.bar.splus1) / (prod.pb0plus.s - prod.pb00.s) # prob01to00 + prob01to01 = 1 as it should - could use 1 - prob01to00 then prob01to10 <- 0 prob01to11 <- 0 prob10to00 <- (prod.pb00.splus1 * prod.pbplus0.bar.splus1 - prod.pb00.s) / (prod.pbplus0.s - prod.pb00.s) prob10to01 <- 0 prob10to10 <- (prod.pbplus0.s - prod.pb00.splus1 * prod.pbplus0.bar.splus1) / (prod.pbplus0.s - prod.pb00.s) prob10to11 <- 0 denominator <- 1 - prod.pbplus0.s - prod.pb0plus.s + prod.pb00.s prob11to00 <- (prod.pb00.splus1 - prod.pb00.splus1 * prod.pbplus0.bar.splus1 - prod.pb00.splus1 * prod.pb0plus.bar.splus1 + prod.pb00.s) / denominator prob11to01 <- (prod.pb0plus.splus1 - prod.pb00.splus1 - prod.pb0plus.s + prod.pb00.splus1 * prod.pb0plus.bar.splus1) / denominator prob11to10 <- (prod.pbplus0.splus1 - prod.pb00.splus1 - prod.pbplus0.s + prod.pb00.splus1 * prod.pbplus0.bar.splus1 ) / denominator prob11to11 <- (1 - prod.pbplus0.splus1 - prod.pb0plus.splus1 + prod.pb00.splus1) / denominator # prob11to00 + prob11to01 + prob11to10 + prob11to11 # = 1 as it should trans.prob.list[[s]] <- c(prob00to00, prob00to01, prob00to10, prob00to11, prob01to00, prob01to01, prob01to10, prob01to11, prob10to00, prob10to01, prob10to10, prob10to11, prob11to00, prob11to01, prob11to10, prob11to11) } trans.prob <- rep(x = trans.prob.list[[1]], each = 2 * 4^(S-3)) # This is for stage 1 to stage 2. There are S - 3 transitions remaining that each have 4 possible values. # At stage S, there are 2 possible values for g_Sak. if (S > 3) { for (s in 2:(S-2)) { trans.prob <- rep(x = rep(x = trans.prob.list[[s]], each = 2 * 4^(S-s-2)), times = 4^(s-1)) * trans.prob # For rep(x = trans.prob.list[[s]], each = 2 * 4^(S-s-2)): There are (S-s-2) stages remaining that each # have 4 possible values. # For rep(x = rep(... ), times = 4^(s-1)): There are 4^(s-1) stages prior to this. # Example with S = 4: The "each = 2 * 4^(S-s-2)" = 2 due to two possible values of g4ak. # Also, the "times = 4^(s-1)" is due to the four possible values of g~_1 = (g~_11, g~_12). } } } # Example ordering for S = 3 # g1k g1kbar g2k g2kbar g3ak # 0 0 0 0 0 # 0 0 0 0 1 each = 2 - notice duplicates for each value of g3ak # 0 0 0 1 0 # 0 0 0 1 1 # 0 0 1 0 0 # 0 0 1 0 1 # 0 0 1 1 0 # 0 0 1 1 1 # 0 1 0 0 0 # 0 1 0 0 1 # ... # Need to include new transition probabilities in the sum probYak <- sum(test.error.part * prob.Sminus1.k.kbar.to.Sk * trans.prob * theta) ################################## # Calculate PSp PSe <- PSeAllStages(joint.p = joint.p, group.mem = group.mem, a = a, k = k, infections = infections, Se = Se, Sp = Sp, ordering = ordering) # From PSeAllStages() - could have this function return a list instead and grab these values out of it index.pa.1andPlus <- all.one.k(k = k, K = K, ordering = ordering) pa.1andPlus <- sum(joint.p[index.pa.1andPlus, a]) # For k = 1, p_a1+ PSp <- 1 - (probYak - PSe*pa.1andPlus)/(1-pa.1andPlus) PPPV <- pa.1andPlus * PSe / (pa.1andPlus * PSe + (1 - pa.1andPlus) * (1 - PSp)) PNPV <- (1 - pa.1andPlus) * PSp / ((1 - pa.1andPlus) * PSp + pa.1andPlus * (1 - PSe)) list(PSe = PSe, PSp = PSp, PPPV = PPPV, PNPV = PNPV) } ################################################################################ #