Nothing
###############################################
# --------------------------------------------#
# Functions for the doubly t-distribution #
# --------------------------------------------#
###############################################
#' Approximation of the distribution function of the doubly noncentral t-distribution
#'
#' Computes the value of the distribution function of the doubly noncentral t-distribution at \code{x}.
#'
#' @inheritParams overview
#' @keywords internal
#' @return Distribution value of the doubly noncentral t-distribution at \code{x}.
doublyT <- function(x, df, lambda1, lambda2, lb = 0, ub) {
values <- lb:ub
suppressWarnings(sum(dpois(values, lambda2/2)*pt(x*(1+2*values/df)^(0.5), df+2*values, lambda1)))
}
#' Calculation of the NCPs of each randomization sequence for the doubly noncentral t-distribution
#'
#' Computes the noncentrality parameters delta and lambda for the doubly noncentral t-distribution of each randomization sequence.
#'
#' @param randSeq object of the class randSeq.
#' @param bias object of the class bias.
#' @param endp object of the class endpoint.
#' @keywords internal
#' @return matrix containing the noncentrality parameters delta and lambda of all randomization sequences.
genNcps <- function(randSeq, bias, endp) {
stopifnot(is(randSeq, "randSeq"), randSeq@K == 2,
is(bias, "issue"), is(endp, "endpoint"), sum(duplicated(endp@sigma)) == 1)
# all randomization sequences
allRandSeq <- randSeq@M
# all expectation values in form of a matrix
allExp <- getExpectation(randSeq, bias, endp)
# number of randomization sequences
r <- nrow(allRandSeq)
P <- lapply(1:r, function(i) {
randSeq <- allRandSeq[i, ]
exp <- allExp[i, ]
splitGroups <- split(exp, randSeq)
# average expectations in both groups
avExp <- lapply(splitGroups, mean)
# variance of the expectations in both groups
varExp <- lapply(splitGroups, function(x) sum((x-mean(x))^2))
# number of assigend patients to both of the groups
numAssPat <- lapply(splitGroups, length)
# defining P
P <- matrix(0, nrow = 1, ncol = 5)
# updating the matrix P
P[1, 3] <- do.call("sum", numAssPat)
P[1, 4] <- ifelse(!is.null(numAssPat$"0"), numAssPat$"0", 0)
P[1, 5] <- ifelse(!is.null(numAssPat$"1"), numAssPat$"1", 0)
if (!is.null(numAssPat$"0") && !is.null(numAssPat$"1")) {
P[1, 1] <- 1/endp@sigma[1] * sqrt((P[1, 4]*P[1, 5]) / (P[1, 4]+P[1, 5])) * (avExp$"0" - avExp$"1")
P[1, 2] <- 1/endp@sigma[1]^2 * (varExp$"0" + varExp$"1")
return(P)
} else {
P[1, 1] <- 0
P[1, 2] <- sum((exp - mean(exp))^2)/endp@sigma[1]^2
return(P)
}
}
)
P <- do.call("rbind", P)
colnames(P) <- c("lambda1", "lambda2", "N", "nA", "nB")
P
}
#' Calculation of the biased type-one-error probability (resp. power) of Student`s t-test
#'
#' Computes the biased type-one-error probability (resp. power) of Student`s t-test due to shifts in the expectation vectors
#' in both treatment groups.
#'
#' @param randSeq object of the class randSeq.
#' @param bias object of the class bias.
#' @param endp object of the class endpoint.
#' @keywords internal
#' @return the biased type-one-error probability (resp. power) of all randomization sequences.
doublyTValues <- function(randSeq, bias, endp) {
stopifnot(is(randSeq, "randSeq"), randSeq@K == 2,
is(bias, "issue"), is(endp, "endpoint"), sum(duplicated(endp@sigma)) == 1)
# calculation of of the matrix containing the ncp
ncps <- genNcps(randSeq, bias, endp)
# variable for the defined alpha quantile
alpha <- bias@alpha
# function for calculating the p values of the singular randomization sequences
p.value <- sapply(1:nrow(ncps), function(i) {
x <- ncps[i, ]
# return zero if in one treatment group was no observation
if( x[4] == 0 || x[5] == 0)
return(0)
# lower boundary
lb <- max(floor(x[2]/2 - qpois(.995, x[2]/2)), 0)
# upper boundary
ub <- as.vector(ceiling(x[2]/2 + qpois(.995, x[2]/2)))
# degrees of freedom
df <- as.vector(x[3]) - 2
# t quantiles
tQuantLow <- qt(alpha/2, df)
tQuantUpper <- -tQuantLow
p.value.less <- doublyT(tQuantLow , df, as.vector(x[1]), as.vector(x[2]), lb, ub)
p.value.greater <- 1 - doublyT(tQuantUpper, df, as.vector(x[1]), as.vector(x[2]), lb, ub)
p.value <- p.value.less + p.value.greater
return(p.value)
}
)
p.value
}
#' Calculation of the NCPs of each randomization sequence for the doubly noncentral t-distribution
#'
#' Computes the noncentrality parameters delta and lambda for the doubly noncentral t-distribution of each randomization sequence.
#'
#' @param randSeq a list of randSeq(rCrSeq or others) with possible varying N's
#' @param bias a list of biases - corresponding to the different randSeq's
#' @param endp object of the class endpoint.
#' @param weight if set to TRUE the weight will be set to 1, according to the paper
#'
#' @return a list containing the noncentrality parameters delta and lambda of all randomization sequences.
#' @export
genNcps_new <- function(randSeq, bias, endp, weight = FALSE){
# we need to iterate over all randSeq and bias objects
iter = length(randSeq)
# number of patients allocated to treatment 0
n_C <- lapply(1:iter, function(i){
ass_matrix <- randSeq[[i]]@M
sapply(1:dim(ass_matrix)[1], function(x) sum(ass_matrix[x,]))
})
# number of patients allocated to treatment 1
n_E <- lapply(1:iter, function(i){
temp_N <- rep(randSeq[[i]]@N, dim(randSeq[[1]]@M)[1])
temp_N - n_C[[i]]
})
# get expectatoins based on the bias chosen
expectations <- list()
for(i in 1:iter){
expectations[[i]] <- matrix(0, dim(randSeq[[i]]@M)[1], randSeq[[i]]@N)
for(x in 1:dim(randSeq[[i]]@M)[1]){
temp_randSeq <- genSeq(crPar(randSeq[[i]]@N))
temp_randSeq@M <- t(matrix(randSeq[[i]]@M[x, ]))
expectations[[i]][x, ] <- getExpectation(temp_randSeq, bias[[i]], endp);
}
}
#
avgExp <- list()
varExp <- list()
for(i in 1:iter){
avgExp[[i]] <- matrix(0, dim(randSeq[[i]]@M)[1], 2)
varExp[[i]] <- matrix(0, dim(randSeq[[i]]@M)[1], 2)
for(x in 1:dim(randSeq[[i]]@M)[1]){
splitGroups <- split(expectations[[i]][x,], randSeq[[i]]@M[x,])
# average expectations in both groups
avgExp[[i]][x,1] <- mean(splitGroups$`0`)
avgExp[[i]][x,2] <- mean(splitGroups$`1`)
varExp[[i]][x,1] <- sum((splitGroups$`0`-mean(splitGroups$`0`))^2)
varExp[[i]][x,2] <- sum((splitGroups$`1`-mean(splitGroups$`1`))^2)
}
avgExp[[i]][which(is.na(avgExp[[i]]))] <- 0
}
temp_delta <- list()
norm_factor <- list()
for(i in 1:iter){
if(weight){
temp_delta[[i]] <- sapply(1:dim(randSeq[[i]]@M)[1], function(x){(avgExp[[i]][x,1]-avgExp[[i]][x,2])})
}else{
temp_delta[[i]] <- sapply(1:dim(randSeq[[i]]@M)[1], function(x){(n_E[[i]][x]*n_C[[i]][x]/(n_E[[i]][x]+n_C[[i]][x]))*(avgExp[[i]][x,1]-avgExp[[i]][x,2])})
}
temp_delta[[i]][which(is.na(temp_delta[[i]]))] <- 0
norm_factor[[i]] <- sapply(1:dim(randSeq[[i]]@M)[1], function(x){n_E[[i]][x]*n_C[[i]][x]/(n_E[[i]][x]+n_C[[i]][x])})
}
weight_no_star <- Reduce(`+`, norm_factor)
if(weight)
weight_no_star <- rep(1, length(norm_factor[[1]]))
# delta - non centrality parameter
#delta = do.call(sum, temp_delta)/(endp@sigma[1]*sqrt(do.call(sum, norm_factor)))
delta = Reduce(`+`, temp_delta)/(endp@sigma[1]*sqrt(weight_no_star^2/Reduce(`+`, norm_factor)))
delta[which(is.na(delta))] <- 0
temp_expSum <- list()
for(i in 1:iter){
temp_expSum[[i]] <- rowSums(expectations[[i]]^2) - n_E[[i]]*(avgExp[[i]][,1]^2) - n_C[[i]]*(avgExp[[i]][,2]^2)
}
# lambda - non centrality parameter
lambda = (1/endp@sigma[1]^2) * Reduce(`+`, temp_expSum)
return(list(delta, lambda))
}
#' Calculating p values
#'
#' Computes the p values based on the noncentrality parameters delta and lambda for the doubly noncentral t-distribution
#'
#' @param delta The first noncentrality parameter
#' @param lambda The second noncentrality parameter
#' @param N the amount of patients in the trial
#' @param alpha significance level
#' @param df degrees of freedom
#'
#' @return a p value
#' @export
get_p_values_new <- function(delta, lambda, N, alpha = 0.05, df = sum(N-2)){
# lower boundary
lb <- max(floor(lambda/2 - qpois(.995, lambda/2)), 0)
# upper boundary
ub <- as.vector(ceiling(lambda/2 + qpois(.995, lambda/2)))
# degrees of freedom
#df <- sum(N-2)
# t quantiles
tQuantLow <- qt(alpha/2, df)
tQuantUpper <- -tQuantLow
p.value.less <- doublyT(tQuantLow , df, delta, lambda, lb, ub)
p.value.greater <- 1 - doublyT(tQuantUpper, df, delta, lambda, lb, ub)
p.value <- p.value.less + p.value.greater
return(p.value)
}
# Function to generate the noncentrality parameters and p-values for a stratified sequence
# Input:
# sequences - the stratified sequence
# bias - the bias for the analysis
# weight - binary value that specifies whether the analysis should be weighted or not
doublyTValues_new<-function(sequences,bias,endp,weight){
tempN <- 0
for(i in 1:length(sequences)){
tempN <- sequences[[i]]@N + tempN
}
N <- tempN
ncps <- genNcps_new(sequences,bias,endp,weight = weight)
delta <- ncps[[1]]
lambda <- ncps[[2]]
pvals <- sapply(1:length(delta), function(i) get_p_values_new(delta = delta[i], lambda = lambda[i], N = N, alpha = 0.05))
pvals
}
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.