Nothing
#### Part 4: from "inside_theta_oct2020.R" ####
insideTheta <- function(n1, n2, n, r1, r2, thetaF, thetaE, cp.sm, cp.m, p0, p, q0=1-p0, q=1-p, mat, coeffs.p0, coeffs, power, alpha){
# To reduce looping, identify the rows which contain a CP < thetaF OR CP > thetaE:
theta.test <- function(SM, M)
{
current.vals <- mat[SM+1, M]
any(current.vals < thetaF | current.vals > thetaE)
}
low.cp.rows <- mapply(theta.test, SM=cp.sm, M=cp.m)
# low.cp.rows <- numeric(length(cp.sm))
# for(i in 1:length(low.cp.rows)){
# current.vals <- mat[cp.sm[i]+1, cp.m[[i]]]
# low.cp.rows[i] <- any(current.vals < thetaF | current.vals > thetaE)
# }
# Because of how CP changes, propagating right to left and upwards, we only
# need to identify the "lowest" row with CP < theta or CP > thetaE, ie the row with greatest Sm
# that contains a CP < theta or CP > thetaE. This is the row where we begin. Note: This may be all rows, or even none!
# Recall that row number = Sm-1
# mat.old <- mat
# We only need to act if there are CPs < thetaF or > thetaE -- otherwise, the matrix of CPs does not change.
if(sum(low.cp.rows)>0) # There may be some cases where there are no CPs to be changed; i.e. no stopping due to stochastic curtailment
{
begin.at.row <- max(which(low.cp.rows>0))
cp.changing.sm <- cp.sm[1:begin.at.row]
cp.changing.m <- cp.m[1:begin.at.row]
# We calculate the CP, round to zero if CP<thetaF (or to 1 if CP > thetaE), then move on:
###### OCT 2020: IMPORTANT: This clearly only works for continuous monitoring, i.e. C=1.
# If you want to use less frequent monitoring, this code (and probably more) will have to change.
for(rown in rev(cp.changing.sm+1))
{
mat.rown <- mat[rown, ]
mat.rown.plus1 <- mat[rown+1, ]
for(coln in rev(cp.changing.m[[rown]]))
{
currentCP <- q*mat.rown[coln+1] + p*mat.rown.plus1[coln+1]
if(currentCP > thetaE){
mat[rown, coln] <- 1 # If CP > thetaE, amend to equal 1
}else{
if(currentCP < thetaF){
mat[rown, coln] <- 0
}else{
mat[rown, coln] <- currentCP
}
}
mat.rown[coln] <- mat[rown, coln]
}
} # Again, this *must* be done one entry at a time -- cannot vectorise.
} # End of IF statement
###### STOP if design is pointless, i.e either failure or success is not possible:
# IF DESIGN GUARANTEES FAILURE (==0) or SUCCESS (==2) at n=1:
first.patient <- sum(mat[,1], na.rm = T)
if(first.patient==2)
{
return(c(n1, n2, n, r1, r2, 1, 1, NA, NA, thetaF, thetaE, NA))
}
if(first.patient==0)
{
return(c(n1, n2, n, r1, r2, 0, 0, NA, NA, thetaF, thetaE, NA))
}
###### At this point, the matrix "mat" contains the CPs, adjusted for stochastic curtailment.
###### The points in the path satisfying 0 < CP < 1 are the possible points for this design;
###### All other points are either terminal (i.e. points of curtailment) or impossible.
###### We now need the characteristics of this design: Type I error, expected sample size, and so on.
pascal.list <- list(1, c(1,1))
for(i in 3:(n+2)){
#column <- as.numeric(mat[!is.na(mat[,i-2]), i-2])
current.col <- mat[, i-2]
column <- as.numeric(current.col[!is.na(current.col)])
#CPzero.or.one <- which(column==0 | column==1)
CPzero.or.one <- column==0 | column==1
newnew <- pascal.list[[i-1]]
newnew[CPzero.or.one] <- 0
pascal.list[[i]] <- c(0, newnew) + c(newnew, 0)
}
pascal.list <- pascal.list[c(-1, -length(pascal.list))]
# Multiply the two triangles (A and p^b * q^c):
final.probs <- Map("*", pascal.list, coeffs)
# for finding type I error prob:
final.probs.p0 <- Map("*", pascal.list, coeffs.p0)
###### We have the probability of each path, taking into account stochastic and NS curtailment.
###### We now must tabulate these paths.
# Might be best to sort into Failure (Sm <= r) and Success (Sm = r+1):
final.probs.mat <- matrix(unlist(lapply(final.probs, '[', 1:max(sapply(final.probs, length)))), ncol = n, byrow = F)
#rownames(final.probs.mat) <- 0:n
final.probs.mat.p0 <- matrix(unlist(lapply(final.probs.p0, '[', 1:max(sapply(final.probs.p0, length)))), ncol = n, byrow = F)
#rownames(final.probs.mat.p0) <- 0:n
# Find successful probabilities first:
# m.success <- (r2+1):n
# Sm.success <- rep(r2+1, length(m.success))
# prob.success <- final.probs.mat[r2+2, m.success]
# prob.success.p0 <- final.probs.mat.p0[r2+2, m.success]
# success.deets <- cbind(Sm.success, m.success, prob.success, prob.success.p0)
# FIRST: Search for terminal points of success. These can only exist in rows where (updated) CP=1, and where Sm<=r+1:
potential.success.rows <- rowSums(mat[1:(r2+2), ]==1, na.rm = TRUE)
rows.with.cp1 <- which(as.numeric(potential.success.rows)>0)
# ^ These are the rows containing possible terminal points of success. They must have CP=1:
columns.of.rows.w.cp1 <- list()
j <- 1
for(i in rows.with.cp1)
{
columns.of.rows.w.cp1[[j]] <- which(mat[i, ]==1 & !is.na(mat[i, ]))
j <- j+1
}
# These rows and columns contain all possible terminal points of success.
# The point CP(Sm, m) is terminal if CP(Sm-1, m-1) < 1 .
# Strictly speaking, CP(Sm, m) is also terminal if CP(Sm, m-1) < 1 .
# However, CP(Sm, m-1) >= CP(Sm-1, m-1) [I think], so the case of
# CP(Sm, m) == 1 AND CP(Sm, m-1) < 1 is not possible.
# DEPRECATED -- TOO SLOW. FASTER CODE BELOW
# success <- NULL
# for(i in 1:length(rows.with.cp1))
# {
# for(j in 1:length(columns.of.rows.w.cp1[[i]]))
# {
# if(mat[rows.with.cp1[i] - 1, columns.of.rows.w.cp1[[i]][j] - 1] < 1) success <- rbind(success, c(rows.with.cp1[i]-1, columns.of.rows.w.cp1[[i]][j]))
# }
# }
index1 <- 1
#success <- NULL
# for(i in rows.with.cp1)
# {
# for(j in columns.of.rows.w.cp1[[index1]])
# {
# if(j==1) success <- rbind(success, c(i-1, j, final.probs.mat[i, j], final.probs.mat.p0[i, j]))
# else
# if(mat[i-1, j-1] < 1) success <- rbind(success, c(i-1, j, final.probs.mat[i, j], final.probs.mat.p0[i, j]))
# }
# index1 <- index1 + 1
# }
success <- vector("list", length(rows.with.cp1)*length(columns.of.rows.w.cp1))
k <- 1
for(i in rows.with.cp1)
{
current.row.probs.mat <- final.probs.mat[i, ]
current.row.probs.mat.p0 <- final.probs.mat.p0[i, ]
for(j in columns.of.rows.w.cp1[[index1]])
{
if(mat[i-1, j-1] < 1 || j==1){
success[[k]] <- c(i-1, j, current.row.probs.mat[j], current.row.probs.mat.p0[j])
k <- k+1
}
}
index1 <- index1 + 1
}
success <- do.call(rbind, success)
colnames(success) <- c("Sm", "m", "prob", "prob.p0")
#
# CAN FIND POWER AND TYPE I ERROR AT THIS POINT.
#
# IF TOO LOW OR HIGH, STOP AND MOVE ON:
#
pwr <- sum(success[,"prob"])
typeI <- sum(success[,"prob.p0"])
if(pwr < power | typeI > alpha) # | pwr > power+tol | alpha < alpha-tol)
{
return(c(n1, n2, n, r1, r2, typeI, pwr, NA, NA, thetaF, thetaE, NA))
}
# Now failure probabilities. Note that there AT MOST one failure probability in each row, and in that
# row the failure probability is the one that has the greatest m (i.e. the "furthest right" non-zero entry):
# First, in the matrix of CPs, find the earliest column containing only zeros and ones (and NAs): This is the latest point ("m") at which the trial stops:
first.zero.one.col <- which(apply(mat, 2, function(x) sum(!(x %in% c(0,1, NA))))==0)[1]
# Second, find the the row at which the final zero exists in this column:
final.zero.in.col <- max(which(mat[, first.zero.one.col]==0))
# Finally, find the rightmost non-zero probability in each row up and including the one above:
m.fail <- NULL
prob.fail <- NULL
prob.fail.p0 <- NULL
for(i in 1:final.zero.in.col)
{
m.fail[i] <- max(which(final.probs.mat[i ,]!=0))
prob.fail[i] <- final.probs.mat[i, m.fail[i]]
prob.fail.p0[i] <- final.probs.mat.p0[i, m.fail[i]]
}
Sm.fail <- 0:(final.zero.in.col-1)
# # Identify non-zero terms in each row:
# m.fail <- rep(NA, r2+1)
# prob.fail <- rep(NA, r2+1)
# prob.fail.p0 <- rep(NA, r2+1)
#
#
#
# for(i in 1:(r2+1))
# {
# m.fail[i] <- max(which(final.probs.mat[i ,]!=0))
# prob.fail[i] <- final.probs.mat[i, m.fail[i]]
# prob.fail.p0[i] <- final.probs.mat.p0[i, m.fail[i]]
# }
# Sm.fail <- 0:r2
fail.deets <- cbind(Sm.fail, m.fail, prob.fail, prob.fail.p0)
all.deets <- rbind(fail.deets, success)
all.deets <- as.data.frame(all.deets)
all.deets$success <- c(rep("Fail", length(m.fail)), rep("Success", nrow(success)))
names(all.deets) <- c("Sm", "m", "prob", "prob.p0", "success")
output <- all.deets
##################### Now find characteristics of design
#PET.failure <- sum(output$prob[output$m < n & output$Sm <= r2]) # Pr(early termination due to failure)
#PET.failure.p0 <- sum(output$prob.p0[output$m < n & output$Sm <= r2])
#PET.success <- sum(output$prob[output$m < n & output$Sm > r2]) # Pr(early termination due to success)
#PET.success.p0 <- sum(output$prob.p0[output$m < n & output$Sm > r2])
#PET <- PET.failure + PET.success # Pr(Early termination due to failure or success)
#PET.p0 <- PET.failure.p0 + PET.success.p0
#PET.df <- data.frame(PET.failure, PET.success, PET, PET.failure.p0, PET.success.p0, PET.p0)
#power <- sum(output$prob[output$success=="Success"])
#output$obsd.p <- output$Sm/output$m
#output$bias <- output$obsd.p - p
#bias.mean <- sum(output$bias*output$prob)
# wtd.mean(output$bias, weights=output$prob, normwt=TRUE) # Check
#bias.var <- wtd.var(output$bias, weights=output$prob, normwt=TRUE)
sample.size.expd <- sum(output$m*output$prob)
# wtd.mean(output$m, weights=output$prob, normwt=TRUE) # Check
#sample.size <- wtd.quantile(output$m, weights=output$prob, normwt=TRUE, probs=c(0.25, 0.5, 0.75))
sample.size.expd.p0 <- sum(output$m*output$prob.p0)
# wtd.mean(output$m, weights=output$prob.p0, normwt=TRUE) # Check
#sample.size.p0 <- wtd.quantile(output$m, weights=output$prob.p0, normwt=TRUE, probs=c(0.25, 0.5, 0.75))
#alpha <- sum(output$prob.p0[output$success=="Success"])
#print(output)
############ MARCH 2018: EFFECTIVE N. THE "EFFECTIVE N" OF A STUDY IS THE "REAL" MAXIMUM SAMPLE SIZE
######## The point is where every Sm for a given m equals zero or one is necessarily where a trial stops
cp.colsums <- apply(mat, 2, function(x) { sum(x==0, na.rm=TRUE)+sum(x==1, na.rm=TRUE)} ) # Sum the CP values that equal zero or one in each column
# Test these against m+1 -- i.e. for i=m, need sum to be m+1:
effective.n <- min(which(cp.colsums==2:(n+1)))
#output <- list(output, PET.df, mean.bias=bias.mean, var.bias=bias.var, sample.size=sample.size, expd.sample.size=sample.size.expd,
# sample.size.p0=sample.size.p0, expd.sample.size.p0=sample.size.expd.p0, alpha=alpha, power=power)
to.return <- c(n1, n2, n, r1, r2, typeI, pwr, sample.size.expd.p0, sample.size.expd, thetaF, thetaE, effective.n)
#to.return <- mat
to.return
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.