Nothing
R/Performance.CondPoisson.R
In Sequential: Exact Sequential Analysis for Poisson and Binomial Data
Defines functions
Performance.CondPoisson
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 Oct. 27, 2023, 1:07 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions Performance.CondPoisson
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.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.