# R/Performance.CondPoisson.R In Sequential: Exact Sequential Analysis for Poisson and Binomial Data

#### Documented in Performance.CondPoisson

```Performance.CondPoisson <-
function(Inference="exact",cv,StopType="Cases",T="n",K="n",cc,D=0,M=1,RR=2,Tailed="upper"){

if(Tailed!="upper"){stop("For this version of the Sequential package, Performance.CondPoisson works only for 'Tailed=upper'.",call. =FALSE)}

if(Inference!="exact"&Inference!="liberal"&Inference!="conservative"){stop("'Inference' must be settled as 'exact', 'conservative', or 'liberal'.",call. =FALSE)}

# ------------------- INPUT VARIABLE ----------------------------------------------------------
# Inference = approach for the calculations among "exact", "conservative" and "liberal".
# cv = critical value in the scale of the CMaxSPRT statistic
# StopType = "Tal" or "Cases". With StopType="Tal", the upper limit on the time of surveillance
# is given in the scale of Pk/V, and the parameter is "T". With StopType="Cases", the upper limit is given in the scale of the events.
# T = maximum length of surveillance, defined in terms of the relative person-time observed during the surveillance period
# K = maximum length of surveillance, defined in terms of the number of events in the surveillance period
# cc = number of adverse events observed in the historical period
# M = The minimum number of cases for which a signal is allowed to occur
# D = Time < T for first look at the data, defined in terms of the expected counts under H0
# Tailed="upper" (default) for H0:RR<=1, and Tailed="lower" for H0:RR>=1 or Tailed="two" for H0:RR=1
# Tests
if(length(Inference)>1|length(StopType)>1|length(T)>1|length(K)>1|length(cc)>1|length(D)>1|length(M)>1){stop("The inputs cannot be vectors.",call. =FALSE)}
if(sum(StopType==c("Tal","Cases"))==0){stop("StopType must contain a label equal to 'Tal' and 'Cases'.",call. =FALSE)}
if(sum(Inference==c("conservative","exact","liberal"))==0){stop("Inference must contain a label among 'conservative', 'exact' and 'liberal'.",call. =FALSE)}
if(StopType=="Tal"){if(D>T){stop("'D' must be a number in the '[0,T]' interval.",call. =FALSE)}
if(is.numeric(T)==FALSE){stop("T must be a positive number.",call. =FALSE)}else{if(T<=0){stop("T must be a positive number.",call. =FALSE)}}
L<- T
}

if(StopType=="Cases"){if(M>K){stop("M must be a positive integer in the '[0,K]' interval.",call. =FALSE)}
if(is.numeric(K)==FALSE){stop("K must be a positive integer.",call. =FALSE)}else{if(K<=0){stop("K must be a positive integer.",call. =FALSE)}}
}

Late<- D
MinCases<- M

####### Tests to verify the validity of the chosen parameters

teste1<- 0

if(T<=0){teste1<- 1; out<- c("SampleSize must be > 0")}
if(teste1==0 & M>100){teste1<- 1; out<- c("M must be a positive integer in the range [1,100]")}
#if(teste1==0 & T>1000){teste1<- 1; out<- c("Use T<=1000")}

if(Late>T & teste1==0){teste1<- 1; out<- c("D must be <= SampleSize") }
if(Late<0 & teste1==0){teste1<- 1; out<- c("Negative values for D does not make sense. Use 0<=D<=SampleSize.") }
if(M<1 & teste1==0){teste1<- 1; out<- c("M must be a positive integer in the range[1,100].") }

if(teste1==1){stop(out,call.=FALSE)}

####### end parameters validity tests
# ------------------------------------------------------------

####### Function to calculate cMaxSPRT
# ------------------------------------------------------------
cLLR<- function(k,cc,tal)
{
if(k/cc<=tal){return(0)}else{return(cc*log((cc*(1+tal)/(cc+k)))+k*log((k*(1+tal)/(tal*(cc+k)))))}
}

####### Function to find the thresholds in the 'tal' scale for a given cv
# ------------------------------------------------------------
cv_tal<- function(k,cc,cv)
{
t1<- 0
t2<- k/cc
cvt<- 0
while(max(abs(cvt-cv),abs(t2-t1))>0.00000000001){
tm<- (t1+t2)/2; cvt<- cLLR(k,cc,tm); if(cvt>cv){t1<- tm}else{t2<- tm}
}
return(tm)
}

if(StopType=="Cases"){ks<- matrix(seq(1,K),K,1)}

#----------------------------------------------------------------------------------------------
# Function that calculates the probability of type I error for a given set of IMPUT parameters
#----------------------------------------------------------------------------------------------
Perror_I_A1<- function(cvv){

####### Three functions for the integrated Poisson distribution with respect to the exponential distribution
# ------------------------------------------------------------
cond.ppois<- function(x,tt){if(x<0){return(0)};y<- seq(0,x); return(sum(exp(y*log(tt)+lfactorial(cc+y-1)-lfactorial(y)-lfactorial(cc-1)-(cc+y)*log(tt+1))))}
cond.dpois<- function(x,tt){if(x<0){return(0)};return(exp(x*log(tt)+lfactorial(cc+x-1)-lfactorial(x)-lfactorial(cc-1)-(cc+x)*log(tt+1)))}

cond.dpois_mult<- function(x,tt){if(x<0){return(0)}else{return(exp(x*log(tt)-lfactorial(x) ))}}

mu<- cv_tal(1,cc,cvv)
i<- 1
if(StopType=="Tal"){while(mu[i]<T){i<- i+1;mu<- cbind(mu,cv_tal(i,cc,cvv))}; imax=i ; mu[imax]<- T}else{mu = apply(ks,1,cv_tal,cc,cvv);imax=K;T<- max(mu);L<-T}

# imax is the maximum number of cases that will generate a signal.

mtemp = c(0,mu)
mmu = diff(mtemp) 		#The marginal difference of the mu[] vector

mu0<- mu
mmu0<- mmu

imin=MinCases
while (mu0[imin] < Late&imin<length(mu)) imin=imin+1
if(imin>MinCases) {
mu0[imin-1]=Late
mmu0[imin]=mu0[imin]-Late
} #END if imin>MinCases

mu<- mu*RR   # Expected number of cases under a real relative risk equal to RR
mmu<- mmu*RR

if(imin==length(mu)&mu[imin]<Late){stop(c("For this value of 'K', 'D' must be in the [0, ", round(mu[imin],2), ") interval."),call. =FALSE)}

# NOTE: If imax=1, this code will not work

if(imin<imax){

# Defining the p[][] matrix
# -------------------------

p<- matrix(0,imax,imax+1)
#p = seq(length=(imax-1)*imax, from=0, by=0)				# p[i,j] is the probability of having j-1 cases at time mu[i]
#dim(p) = c(imax-1,imax)								# i in 1:imax-1 is the rows and j in 1:imax is the column

# Calculating the first row p[imin][] in the matrix for which there is a chance to reject H0
# When MinCases=1, there is no skipping, and it is the first row in the matrix (p[1][]).
# ------------------------------------------------------------------------------------------

if(imin==MinCases){
for(s in 1:imin) p[imin,s] = cond.dpois_mult(s-1,mu[imin])			# Probability of having s-1 cases at time mu[imin], not rejectinh H0
p[imin,imin+1] = 1-cond.ppois(imin-1,mu[imin])				# Probability of having s+ cases at time mu[imin], rejectinh H0
} # end if

if(imin>MinCases){
for(s in 1:imin) p[imin-1,s]=cond.dpois_mult(s-1,mu[imin-1])		# Probability of having s-1 cases at time mu[imin-1], not rejecting H0
p[imin-1,imin+1] = 1-cond.ppois(imin-1,mu[imin-1])				# Probability of having s+ cases at time mu[imin-1], rejecting H0
for(s in 1:imin) {								            # Probability of having s-1 cases at time mu[imin], not rejectinh H0
for(k in 1:s){
p[imin,s]=p[imin,s]+p[imin-1,k]*cond.dpois_mult(s-k,mmu[imin])
for(k in 1:imin){
partial<- 0
for(xkm1 in (k-1):(imin)){partial<- partial + exp( log(cond.dpois_mult(xkm1-(k-1),mmu[imin])) + log(p[imin-1,k]) +     lfactorial(cc+xkm1-1) - (cc+xkm1)*log(1+mu[imin]) - lfactorial(cc-1) ) }
p[imin,imin+1] = p[imin,imin+1] + p[imin-1,k]*exp( lfactorial(cc+(k-1)-1) - (cc+k-1)*log(1+mu[imin-1]) - lfactorial(cc-1) )  - partial
}
}
}
} # end if

alpharef<- 0
if(imin>MinCases){alpharef=p[imin-1,imin+1]}
alpharef<- alpharef+p[imin,imin+1]

# Calculating the remaining rows in the p[][] matix
# -------------------------------------------------

if(MinCases+1<=imax-1&((imin+1)<=(imax-1))){
i<- (imin+1)
while(i<=(imax-1)){

for(j in 1:(i-1)){							# This loop calculates the p[][] matix, one column at a time, from left to right
for(k in 1:j)
p[i,j]=p[i,j]+p[i-1,k]*cond.dpois_mult(j-k,mmu[i])	# Calculates the standard p[][] cell values
}

for(k in 1:(i-1)){
p[i,i]=p[i,i]+p[i-1,k]*cond.dpois_mult(i-k,mmu[i])		# Calculates the diagonal under the absorbing states, which requires a unique formula
}

for(k in 1:(i-1)){
partial<- 0
for(xkm1 in (k-1):(i-1)){partial<- partial + exp( log(cond.dpois_mult(xkm1-(k-1),mmu[i])) + log(p[i-1,k]) +     lfactorial(cc+xkm1-1) - (cc+xkm1)*log(1+mu[i]) - lfactorial(cc-1) ) }
p[i,i+1]= p[i,i+1] + exp(log(p[i-1,k]) +lfactorial(cc+(k-1)-1) - (cc+k-1)*log(1+mu[i-1]) - lfactorial(cc-1) )  - partial   # Calculates the diagonal absorbing states where H0 is rejected
}
alpharef<- alpharef + p[i,i+1]
i<- i+1

} # end for i

}

pp=0
if(imax>imin){
for(k in 1:(imax-1)){
partial<- 0
for(xkm1 in (k-1):(imax-1)){partial<- partial + exp( log(cond.dpois_mult(xkm1-(k-1),mmu[imax])) + log(p[imax-1,k]) +     lfactorial(cc+xkm1-1) - (cc+xkm1)*log(1+mu[imax]) - lfactorial(cc-1) ) }
pp<- pp + exp(log(p[imax-1,k])+ lfactorial(cc+(k-1)-1) - (cc+k-1)*log(1+mu[imax-1]) - lfactorial(cc-1) )  - partial #Calculates the last probability to signal before time T
}
}

# Sums up the probabilities of absorbing states when a signal occurs, to get the alpha level
# ------------------------------------------------------------------------------------------

Power=0
ESignalTime<- 0
ESampleSize<- 0

if(imin>MinCases){ Power=p[imin-1,imin+1]; ESampleSize<- (imin-1)*p[imin-1,imin+1]}
for(i in imin:(imax-1)){ Power=Power+p[i,i+1] ; ESampleSize<- ESampleSize+ i*p[i,i+1]}
Power=Power+pp
ESampleSize<- ESampleSize+ (imax-1)*pp

ESignalTime<- ESampleSize/Power
ESampleSize<- ESampleSize + imax*(1-Power)

}else{alpha<- 1-cond.ppois(imax-1,mu[imax])} # end if(imin<imax)

return(list(Power,ESignalTime,ESampleSize))

} # end Perror_I_A1
#############################################################################################

#############################################################################################

Perror_I_A2<- function(cvv){

####### Three functions for the integrated Poisson distribution with respect to the exponential distribution
# ------------------------------------------------------------
cond.ppois<- function(x,tt){if(x<0){return(0)};y<- seq(0,x); rest<- apply(matrix(y,ncol=1),1,A,tt) ; return(sum(exp( y*log(tt)-lfactorial(y)+ log(rest) ) ) )}
cond.dpois<- function(x,tt){if(x<0){return(0)};return(exp( x*log(tt)-lfactorial(x)+ log( A(x,tt) ) ) )}
cond.dpois_mult<- function(x,tt){if(x<0){return(0)}else{return(exp(x*log(tt)-lfactorial(x) ))}}

# Another auxiliar function to be used inside term A
B<- function(x,ll,tt){
if(x==ll){return(log(tt+1)/tt)}else{
C<- function(gg){auxmod<- sum(gg-x-ll<0); return( ((-1)^auxmod)*((-1)^(x-ll-gg))*exp(-(x-ll+1)*log(tt)-log(abs(gg-x+ll))+lchoose(x-ll,gg) + (gg-x+ll)*log(tt+1) ) - ((-1)^auxmod)*(-1)^(x-ll-gg)*exp(-(x-ll+1)*log(tt)-log(abs(gg-x+ll))+lchoose(x-ll,gg) ) )}
g<- matrix(seq(0,x-ll-1,1),ncol=1)
return( exp( -(x-ll+1)*log(tt)+ log(log(tt+1)) )+ sum(apply(g,1,C))  )
}
}

A<- function(x,tt){

if(x<0){return(0)}else{
if(x==0){return( exp( -cc*log(tt+1)+log(B(x,0,tt)) ) )}else{  # the term A in the general expression
D<- function(ll){auxB<- B(x,ll,tt);return( (-1)^sum(auxB<0)* exp(lchoose(x,ll)+lfactorial(ll+cc-1)+lfactorial(x-ll)+log(abs(auxB)) -lfactorial(cc-1)-(ll+cc)*log(tt+1) ))}
ll<- matrix(seq(0,x,1),ncol=1)
return(sum(apply(ll,1,D)))}
}
}

mu<- cv_tal(1,cc,cvv)
i<- 1
if(StopType=="Tal"){while(mu[i]<T){i<- i+1;mu<- cbind(mu,cv_tal(i,cc,cvv))}; imax=i ; mu[imax]<- T}else{mu = apply(ks,1,cv_tal,cc,cvv);imax=K;T<- max(mu);L<-T}

# imax is the maximum number of cases that will generate a signal.

mtemp = c(0,mu)
mmu = diff(mtemp) 		#The marginal difference of the mu[] vector

mu0<- mu
mmu0<- mmu

imin=MinCases
while (mu0[imin] < Late&imin<length(mu)) imin=imin+1
if(imin>MinCases) {
mu0[imin-1]=Late
mmu0[imin]=mu0[imin]-Late
} #END if imin>MinCases

mu<- mu*RR   # Expected number of cases under a real relative risk equal to RR
mmu<- mmu*RR

if(imin==length(mu)&mu[imin]<Late){stop(c("For this value of 'K', 'D' must be in the [0, ", round(mu[imin],2), ") interval."),call. =FALSE)}

# NOTE: If imax=1, this code will not work

if(imin<imax){

# Defining the p[][] matrix
# -------------------------

p<- matrix(0,imax,imax+1)
#p = seq(length=(imax-1)*imax, from=0, by=0)				# p[i,j] is the probability of having j-1 cases at time mu[i]
#dim(p) = c(imax-1,imax)								# i in 1:imax-1 is the rows and j in 1:imax is the column

# Calculating the first row p[imin][] in the matrix for which there is a chance to reject H0
# When MinCases=1, there is no skipping, and it is the first row in the matrix (p[1][]).
# ------------------------------------------------------------------------------------------

if(imin==MinCases){
for(s in 1:imin) p[imin,s] = cond.dpois_mult(s-1,mu[imin])			# Probability of having s-1 cases at time mu[imin], not rejectinh H0
p[imin,imin+1] = 1-cond.ppois(imin-1,mu[imin])				# Probability of having s+ cases at time mu[imin], rejectinh H0
} # end if

if(imin>MinCases){
for(s in 1:imin) p[imin-1,s]=cond.dpois_mult(s-1,mu[imin-1])		# Probability of having s-1 cases at time mu[imin-1], not rejecting H0
p[imin-1,imin+1] = 1-cond.ppois(imin-1,mu[imin-1])				# Probability of having s+ cases at time mu[imin-1], rejecting H0
for(s in 1:imin) {								            # Probability of having s-1 cases at time mu[imin], not rejectinh H0
for(k in 1:s){
p[imin,s]=p[imin,s]+p[imin-1,k]*cond.dpois_mult(s-k,mmu[imin])
for(k in 1:imin){
partial<- 0
for(xkm1 in (k-1):(imin)){partial<- partial + exp( log(cond.dpois_mult(xkm1-(k-1),mmu[imin])) + log(p[imin-1,k]) +     lfactorial(cc+xkm1-1) - (cc+xkm1)*log(1+mu[imin]) - lfactorial(cc-1) ) }
p[imin,imin+1] = p[imin,imin+1] + p[imin-1,k]*exp( lfactorial(cc+(k-1)-1) - (cc+k-1)*log(1+mu[imin-1]) - lfactorial(cc-1) )  - partial
}
}
}
} # end if

alpharef<- 0
if(imin>MinCases){alpharef=p[imin-1,imin+1]}
alpharef<- alpharef+p[imin,imin+1]

# Calculating the remaining rows in the p[][] matix
# -------------------------------------------------

if(MinCases+1<=imax-1&((imin+1)<=(imax-1))){
i<- (imin+1)
while(i<=(imax-1)){

for(j in 1:(i-1)){							# This loop calculates the p[][] matix, one column at a time, from left to right
for(k in 1:j)
p[i,j]=p[i,j]+p[i-1,k]*cond.dpois_mult(j-k,mmu[i])	# Calculates the standard p[][] cell values
}

for(k in 1:(i-1)){
p[i,i]=p[i,i]+p[i-1,k]*cond.dpois_mult(i-k,mmu[i])		# Calculates the diagonal under the absorbing states, which requires a unique formula
}

for(k in 1:(i-1)){
partial<- 0
#for(xkm1 in (k-1):(i-1)){partial<- partial + exp( log(cond.dpois_mult(xkm1-(k-1),mmu[i])) + log(p[i-1,k]) +     lfactorial(cc+xkm1-1) - (cc+xkm1)*log(1+mu[i]) - lfactorial(cc-1) ) }
for(xkm1 in (k-1):(i-1)){partial<- partial + exp( log(cond.dpois_mult(xkm1-(k-1),mmu[i])) + log(p[i-1,k]) +  log(A(xkm1,mu[i])) ) }

#p[i,i+1]= p[i,i+1] + exp(log(p[i-1,k]) +lfactorial(cc+(k-1)-1) - (cc+k-1)*log(1+mu[i-1]) - lfactorial(cc-1) )  - partial   # Calculates the diagonal absorbing states where H0 is rejected
p[i,i+1]= p[i,i+1] + exp(log(p[i-1,k]) + log(A(k-1,mu[i-1])) )  - partial   # Calculates the diagonal absorbing states where H0 is rejected

}
alpharef<- alpharef + p[i,i+1]
i<- i+1

} # end for i

}

pp=0
if(imax>imin){
for(k in 1:(imax-1)){
partial<- 0
for(xkm1 in (k-1):(imax-1)){partial<- partial + exp( log(cond.dpois_mult(xkm1-(k-1),mmu[imax])) + log(p[imax-1,k]) +     lfactorial(cc+xkm1-1) - (cc+xkm1)*log(1+mu[imax]) - lfactorial(cc-1) ) }
pp<- pp + exp(log(p[imax-1,k])+ lfactorial(cc+(k-1)-1) - (cc+k-1)*log(1+mu[imax-1]) - lfactorial(cc-1) )  - partial #Calculates the last probability to signal before time T
}
}

# Sums up the probabilities of absorbing states when a signal occurs, to get the alpha level
# ------------------------------------------------------------------------------------------

Power=0
ESignalTime<- 0
ESampleSize<- 0

if(imin>MinCases){ Power=p[imin-1,imin+1]; ESampleSize<- (imin-1)*p[imin-1,imin+1]}
for(i in imin:(imax-1)){ Power=Power+p[i,i+1] ; ESampleSize<- ESampleSize+ i*p[i,i+1]}
Power=Power+pp
ESampleSize<- ESampleSize+ (imax-1)*pp

ESignalTime<- ESampleSize/Power
ESampleSize<- ESampleSize + imax*(1-Power)

}else{alpha<- 1-cond.ppois(imax-1,mu[imax])} # end if(imin<imax)

return(list(Power,ESignalTime,ESampleSize))

} # end Perror_I_A2
#############################################################################################

#############################################################################################

Perror_I_A3<- function(cvv){

####### Three functions for the integrated Poisson distribution with respect to the exponential distribution
# ------------------------------------------------------------
cond.ppois<- function(x,tt){if(x<0){return(0)};y<- seq(0,x); rest<- apply(matrix(y,ncol=1),1,A,tt) ; return(sum(exp( y*log(tt)-lfactorial(y)+ log(rest) ) ) )}
cond.dpois<- function(x,tt){if(x<0){return(0)};return(exp( x*log(tt)-lfactorial(x)+ log( A(x,tt) ) ) )}
cond.dpois_mult<- function(x,tt){if(x<0){return(0)}else{return(exp(x*log(tt)-lfactorial(x) ))}}

# Another auxiliar function to be used inside term A
B<- function(x,ll,tt){
return(1/((tt+1)^(x-ll+1)))
}

A<- function(x,tt){

if(x<0){return(0)}else{
if(x==0){return( exp( -cc*log(tt+1)+log(B(x,0,tt)) ) )}else{  # the term A in the general expression
D<- function(ll){auxB<- B(x,ll,tt);return( (-1)^sum(auxB<0)* exp(lchoose(x,ll)+lfactorial(ll+cc-1)+lfactorial(x-ll)+log(abs(auxB)) -lfactorial(cc-1)-(ll+cc)*log(tt+1) ))}
ll<- matrix(seq(0,x,1),ncol=1)
return(sum(apply(ll,1,D)))}
}
}

mu<- cv_tal(1,cc,cvv)
i<- 1
if(StopType=="Tal"){while(mu[i]<T){i<- i+1;mu<- cbind(mu,cv_tal(i,cc,cvv))}; imax=i ; mu[imax]<- T}else{mu = apply(ks,1,cv_tal,cc,cvv);imax=K;T<- max(mu);L<-T}

# imax is the maximum number of cases that will generate a signal.

mtemp = c(0,mu)
mmu = diff(mtemp) 		#The marginal difference of the mu[] vector

mu0<- mu
mmu0<- mmu

imin=MinCases
while (mu0[imin] < Late&imin<length(mu)) imin=imin+1
if(imin>MinCases) {
mu0[imin-1]=Late
mmu0[imin]=mu0[imin]-Late
} #END if imin>MinCases

mu<- mu*RR   # Expected number of cases under a real relative risk equal to RR
mmu<- mmu*RR

if(imin==length(mu)&mu[imin]<Late){stop(c("For this value of 'K', 'D' must be in the [0, ", round(mu[imin],2), ") interval."),call. =FALSE)}

# NOTE: If imax=1, this code will not work

if(imin<imax){

# Defining the p[][] matrix
# -------------------------

p<- matrix(0,imax,imax+1)
#p = seq(length=(imax-1)*imax, from=0, by=0)				# p[i,j] is the probability of having j-1 cases at time mu[i]
#dim(p) = c(imax-1,imax)								# i in 1:imax-1 is the rows and j in 1:imax is the column

# Calculating the first row p[imin][] in the matrix for which there is a chance to reject H0
# When MinCases=1, there is no skipping, and it is the first row in the matrix (p[1][]).
# ------------------------------------------------------------------------------------------

if(imin==MinCases){
for(s in 1:imin) p[imin,s] = cond.dpois_mult(s-1,mu[imin])			# Probability of having s-1 cases at time mu[imin], not rejectinh H0
p[imin,imin+1] = 1-cond.ppois(imin-1,mu[imin])				# Probability of having s+ cases at time mu[imin], rejectinh H0
} # end if

if(imin>MinCases){
for(s in 1:imin) p[imin-1,s]=cond.dpois_mult(s-1,mu[imin-1])		# Probability of having s-1 cases at time mu[imin-1], not rejecting H0
p[imin-1,imin+1] = 1-cond.ppois(imin-1,mu[imin-1])				# Probability of having s+ cases at time mu[imin-1], rejecting H0
for(s in 1:imin) {								            # Probability of having s-1 cases at time mu[imin], not rejectinh H0
for(k in 1:s){
p[imin,s]=p[imin,s]+p[imin-1,k]*cond.dpois_mult(s-k,mmu[imin])
for(k in 1:imin){
partial<- 0
for(xkm1 in (k-1):(imin)){partial<- partial + exp( log(cond.dpois_mult(xkm1-(k-1),mmu[imin])) + log(p[imin-1,k]) +     lfactorial(cc+xkm1-1) - (cc+xkm1)*log(1+mu[imin]) - lfactorial(cc-1) ) }
p[imin,imin+1] = p[imin,imin+1] + p[imin-1,k]*exp( lfactorial(cc+(k-1)-1) - (cc+k-1)*log(1+mu[imin-1]) - lfactorial(cc-1) )  - partial
}
}
}
} # end if

alpharef<- 0
if(imin>MinCases){alpharef=p[imin-1,imin+1]}
alpharef<- alpharef+p[imin,imin+1]

# Calculating the remaining rows in the p[][] matix
# -------------------------------------------------

if(MinCases+1<=imax-1&((imin+1)<=(imax-1))){
i<- (imin+1)
while(i<=(imax-1)){

for(j in 1:(i-1)){							# This loop calculates the p[][] matix, one column at a time, from left to right
for(k in 1:j)
p[i,j]=p[i,j]+p[i-1,k]*cond.dpois_mult(j-k,mmu[i])	# Calculates the standard p[][] cell values
}

for(k in 1:(i-1)){
p[i,i]=p[i,i]+p[i-1,k]*cond.dpois_mult(i-k,mmu[i])		# Calculates the diagonal under the absorbing states, which requires a unique formula
}

for(k in 1:(i-1)){
partial<- 0
#for(xkm1 in (k-1):(i-1)){partial<- partial + exp( log(cond.dpois_mult(xkm1-(k-1),mmu[i])) + log(p[i-1,k]) +     lfactorial(cc+xkm1-1) - (cc+xkm1)*log(1+mu[i]) - lfactorial(cc-1) ) }
for(xkm1 in (k-1):(i-1)){partial<- partial + exp( log(cond.dpois_mult(xkm1-(k-1),mmu[i])) + log(p[i-1,k]) +  log(A(xkm1,mu[i])) ) }

#p[i,i+1]= p[i,i+1] + exp(log(p[i-1,k]) +lfactorial(cc+(k-1)-1) - (cc+k-1)*log(1+mu[i-1]) - lfactorial(cc-1) )  - partial   # Calculates the diagonal absorbing states where H0 is rejected
p[i,i+1]= p[i,i+1] + exp(log(p[i-1,k]) + log(A(k-1,mu[i-1])) )  - partial   # Calculates the diagonal absorbing states where H0 is rejected

}
alpharef<- alpharef + p[i,i+1]
i<- i+1

} # end for i

}

pp=0
if(imax>imin){
for(k in 1:(imax-1)){
partial<- 0
for(xkm1 in (k-1):(imax-1)){partial<- partial + exp( log(cond.dpois_mult(xkm1-(k-1),mmu[imax])) + log(p[imax-1,k]) +     lfactorial(cc+xkm1-1) - (cc+xkm1)*log(1+mu[imax]) - lfactorial(cc-1) ) }
pp<- pp + exp(log(p[imax-1,k])+ lfactorial(cc+(k-1)-1) - (cc+k-1)*log(1+mu[imax-1]) - lfactorial(cc-1) )  - partial #Calculates the last probability to signal before time imax
}
}

# Sums up the probabilities of absorbing states when a signal occurs, to get the alpha level
# ------------------------------------------------------------------------------------------

Power=0
ESignalTime<- 0
ESampleSize<- 0

if(imin>MinCases){ Power=p[imin-1,imin+1]; ESampleSize<- (imin-1)*p[imin-1,imin+1]}
for(i in imin:(imax-1)){ Power=Power+p[i,i+1] ; ESampleSize<- ESampleSize+ i*p[i,i+1]}
Power=Power+pp
ESampleSize<- ESampleSize+ (imax-1)*pp

ESignalTime<- ESampleSize/Power
ESampleSize<- ESampleSize + imax*(1-Power)

}else{alpha<- 1-cond.ppois(imax-1,mu[imax])} # end if(imin<imax)

return(list(Power,ESignalTime,ESampleSize))

} # end Perror_I_A3

if(Inference=="liberal"){res<- Perror_I_A1(cvv=cv)}
if(Inference=="exact"){res<- Perror_I_A2(cvv=cv)}
if(Inference=="conservative"){res<- Perror_I_A3(cvv=cv)}

if(StopType=="Cases"&T!="n"){message("T was ignored since 'StopType=Cases'")}
if(StopType=="Tal"&K!="n"){message("K was ignored since 'StopType=Tal'")}

names(res)<- c("Power","ESignalTime","ESampleSize")
return(res)

} # End of the CV.CondPoisson function

### Example:

# Critical values from each approach

#cv1<- CV.CondPoisson(Inference="liberal",StopType="Cases",K=50,cc=50,D=0,M=1,alpha=0.05)  # This is the liberal approach (actual alpha > target alpha).

# Actual statistical power, expected time to signal, and expected length of surveillance:

#system.time(per<- Performance.CondPoisson(Inference="exact",cv=as.numeric(cv1[[2]]),StopType="Cases",K=50,cc=50,D=0,M=1,RR=2))
```

## Try the Sequential package in your browser

Any scripts or data that you put into this service are public.

Sequential documentation built on April 24, 2022, 9:05 a.m.