#################################################################################### #################################################################################### # # # R functions corresponding to the manuscript entitled "Group testing for multiple # # infections with applications to the Infertility and Prevention Project" # # # #################################################################################### #################################################################################### ############################################ # Functions used to find E(T_k) # # Input: # p00,p10,p01,p11 = Cell probabilities # Se1,Se2 = Sensitivities # Sp1,Sp2 = Specificities # c = Pool size # # Output: # res = E(T_k) # ET<-function(p00,p01,p10,Se1,Se2,Sp1,Sp2,c){ theta00<-p00^c theta01<-(p01+p00)^c-p00^c theta10<-(p10+p00)^c-p00^c theta11<-1-theta00-theta01-theta10 res<-1+c*(1-Sp1*Sp2*theta00-Sp1*(1-Se2)*theta01-(1-Se1)*Sp2*theta10-(1-Se1)*(1-Se2)*theta11) return(res) } ############################################## # Functions used to find the optimal pool size # # Input: # P=c(p00,p10,p01,p11) = Cell probabilities # Se=c(Se1,Se2) = Sensitivities # Sp=c(Sp1,Sp2) = Specificities # # Output: # psz = Optimal pool size # Opt.pool.size<-function(P,Se,Sp){ p00<-P[1] p10<-P[2] p01<-P[3] p11<-P[4] Se1<-Se[1] Se2<-Se[2] Sp1<-Sp[1] Sp2<-Sp[2] Et0<-10 Et1<-1 c<-1 while(Et1 <= Et0){ Et0<-Et1 c<-c+1 Et1<-ET(p00,p01,p10,Se1,Se2,Sp1,Sp2,c)/c } psz<-c-1 return(psz) } ############################################ # Functions used to find PS_{e:1} # # Input: # p00,p10,p01,p11 = Cell probabilities # Se1,Se2 = Sensitivities # Sp1,Sp2 = Specificities # c = Pool size # # Output: # res = PS_{e:1} # PSe1<-function(p00,p10,p11,Se1,Se2,Sp2,c){ res<-Se1^2 +(1-Se1)*Se1*(Se2+(1-Se2-Sp2)*((p00+p10)^(c-1))*p10/(p10+p11)) return(res) } ############################################ # Functions used to find PS_{e:2} # # Input: # p00,p10,p01,p11 = Cell probabilities # Se1,Se2 = Sensitivities # Sp1,Sp2 = Specificities # c = Pool size # # Output: # res = PS_{e:2} # PSe2<-function(p00,p01,p11,Se1,Se2,Sp1,c){ res<-PSe1(p00,p01,p11,Se2,Se1,Sp1,c) return(res) } ############################################ # Functions used to find PS_{p:1} # # Input: # p00,p10,p01,p11 = Cell probabilities # Se1,Se2 = Sensitivities # Sp1,Sp2 = Specificities # c = Pool size # # Output: # res = PS_{p:1} # PSp1<-function(p00,p01,p10,Se1,Se2,Sp1,Sp2,c){ theta00<-p00^c theta01<-(p01+p00)^c-p00^c theta10<-(p10+p00)^c-p00^c theta11<-1-theta00-theta01-theta10 a<-(1-Sp1)*(1-Sp1)*(p00+p01)^(c-1) b<-Se1*(1-Sp1)*(1-(p00+p01)^(c-1)) c<-Sp1*(1-Sp2)*(1-Sp1)*theta00/(p01+p00) d<-(1-Se1)*(1-Sp2)*(1-Sp1)*(theta10-((p10+p00)^(c-1))*p10)/(p00+p01) e<-Sp1*Se2*(1-Sp1)*theta01/(p00+p01) f<-(1-Se1)*Se2*(1-Sp1)*(theta11-(1-(p00+p01)-(p00+p10)^(c-1)*p10))/(p00+p01) res<-1-a-b-c-d-e-f return(res) } ############################################ # Functions used to find PS_{p:2} # # Input: # p00,p10,p01,p11 = Cell probabilities # Se1,Se2 = Sensitivities # Sp1,Sp2 = Specificities # c = Pool size # # Output: # res = PS_{p:2} # PSp2<-function(p00,p01,p10,Se1,Se2,Sp1,Sp2,c){ res<-PSp1(p00,p10,p01,Se2,Se1,Sp2,Sp1,c) return(res) } ################################################################################# # EM algorithm for group testing data under the Iowa IPP pooling algorithm # # Note: This function makes use of two external fortran 77 dlls, specifically "multigibbs.dll" and "louisvar.dll." # # Input: # data = Matrix(\widetilde{Y}_{i1k},\widetilde{Y}_{i2k},Z_{1k},Z_{2k},Y_{i1k},Y_{i2k},c_k, c(Vector that identifies the other individuals in the kth pool)) # p0 = initial guess at parameter # Se = Vector of test sensitivities # Sp = Vector of test specificities # GI = Number of Gibbs iterates # a = Burn in # eps = Converegence criteria # # Output: # p = Estimate of the parameter Note: p=(p00,p10,p01,p11) # var = Covariance matrix of the estimates Note: var=Var((p00,p10,p01)') Dual.GT.DOrfman.EM<-function(data,p0,Se,Sp,GI,a,eps){ dyn.load("multigibbs.dll") dyn.load("louisvar.dll") p1<-p0 p0<-p0+2*eps N<-dim(data)[1] ########################### # Start iterating while(max(abs(p1-p0))>eps){ p0<-p1 W<-matrix(0,nrow=N,ncol=4) u<-matrix(runif(GI*N),nrow=N,ncol=GI) res<-.C("multigibbs_",as.integer(N),as.integer(data),as.integer(dim(data)[2]),as.integer(GI),as.integer(a),Se,Sp,p0,as.integer(W),u) Wi<-matrix(res[[9]],nrow=N,ncol=4) Wi<-Wi/(GI-a) p1[1]<-mean(Wi[,1]) p1[2]<-mean(Wi[,2]) p1[3]<-mean(Wi[,3]) p1[4]<-mean(Wi[,4]) } dyn.unload("multigibbs.dll") ####################################### # First step in finding the varaince of # the parameter estimates samp<-matrix(0,nrow=GI,ncol=3) u<-matrix(runif(GI*N),nrow=N,ncol=GI) res<-.C("louisvar_",as.integer(N),as.integer(data),as.integer(dim(data)[2]),as.integer(GI),as.integer(a),Se,Sp,p1,samp,u) samp<-matrix(res[[9]],nrow=GI,ncol=3) Q.mat<-matrix(0,3,3) for(i in 1:N){ Q.mat<-Q.mat+diag(as.vector(-Wi[i,1:3]/(p1[1:3]^2)))-(Wi[i,4]/(p1[4]^2)) } dyn.unload("louisvar.dll") var<--(cov(samp)+Q.mat) return(list("p"=p1,"var"=var)) }