erlotinib_adrHCC827_diffexDT susceptibleCellLine_BetaUniformModel <- fitBetaUniformMixtureDistribution(erlotinib_adrHCC827_diffexDT$HCC827_pValue, nStarts = 20) fittedBetaShape <- susceptibleCellLine_BetaUniformModel@a uniformProportion <- susceptibleCellLine_BetaUniformModel@lambda fitBetaUniformMixtureDistribution(rbetauniform(100000,fittedBetaShape,nonUniformProportion)) plot.BetaUniformModel(susceptibleCellLine_BetaUniformModel, outputParameters = F, outputFormula = F) qplot(rbetauniform(10000, a = susceptibleCellLine_BetaUniformModel@a, lambda = susceptibleCellLine_BetaUniformModel@lambda)) nullModelTests <- list() for(i in 1:1E4){ randomSetSize <- rpois(1,10) random_draws <- map(1:1E4, ~{rbetauniform(randomSetSize, fittedBetaShape, uniformProportion)}) # random_draws <- map(1:1000, ~{runif(n = randomSetSize)}); uniformProportion <- 1 #random_draws <- map(1:1000, ~{rbeta(n = randomSetSize, shape1 = fittedBetaShape, shape2 = 1)}) ; uniformProportion <- 0 testStats <- map_dbl(random_draws, ~ sum(-log(.x))) empiricalGTE <- length(which(testStats[1] <= testStats[-1]))/length(testStats[-1]) fisherPval <- pgamma(testStats[1], shape = randomSetSize, scale = 1, lower.tail = FALSE) mu <- randomSetSize*(uniformProportion + (1-uniformProportion)*scaleParam) sigmaSq <- randomSetSize*(uniformProportion + (1-uniformProportion)*scaleParam^2) combinedPval <- pgamma(testStats[1], shape = randomSetSize*uniformProportion, scale = 1, lower.tail = FALSE) + pgamma(testStats[1], shape = randomSetSize*(1-uniformProportion), rate = fittedBetaShape, lower.tail = FALSE) # hyperexponentialSumTransitionMatrix <- constructHyperexponentialSumTransitionMatrix(randomSetSize, uniformProportion, fittedBetaShape) # # initialProb <- Matrix(0,ncol = ncol(hyperexponentialSumTransitionMatrix), nrow = 1) # initialProb[1:2] <- c(uniformProportion, 1-uniformProportion) # CDF of phase rate distribution os 1- alpha * exp(x*S)*oneVec # phaseRateP <- sum(initialProb %*% Matrix::expm(testStats[1]*hyperexponentialSumTransitionMatrix)) phaseRateP <- betaUniformPvalueSumTest(random_draws[[1]], susceptibleCellLine_BetaUniformModel) # Welch-Satterthwaite Approximation for a sum of exponentials, some with different paramenters approxPval <- pgamma(testStats[1], shape = (mu^2)/sigmaSq, scale = sigmaSq/mu, lower.tail = F) nullModelTests %<>% c(list(data.table(empiricalGTE, fisherPval, combinedPval, approxPval, phaseRateP, size = randomSetSize))) } microbenchmark::microbenchmark( combinedPval <- pgamma(testStats[1], shape = randomSetSize*uniformProportion, scale = 1, lower.tail = FALSE) + pgamma(testStats[1], shape = randomSetSize*(1-uniformProportion), rate = fittedBetaShape, lower.tail = FALSE), hyperexponentialSumTransitionMatrix <- constructHyperexponentialSumTransitionMatrix(randomSetSize, uniformProportion, fittedBetaShape), phaseRateP <- pphtype(testStats[1], prob = c(c(uniformProportion, 1-uniformProportion), rep(0,2*(randomSetSize-1))), rates = as.matrix(hyperexponentialSumTransitionMatrix), lower.tail = FALSE), # Welch-Satterthwaite Approximation for a sum of exponentials, some with different paramenters approxPval <- pgamma(testStats[1], shape = (mu^2)/sigmaSq, scale = sigmaSq/mu, lower.tail = F) ) ggplot(data = data.frame(x = c(0, 100)), aes(x)) + stat_function(fun = dgamma, args = list(rate = fittedBetaShape, shape = randomSetSize), colour = "red") + stat_function(fun = dgamma, args = list(rate = 1, shape = randomSetSize), colour = "black") + stat_function(fun = dgamma, args = list(scale = sigmaSq/mu, shape = (mu^2)/sigmaSq), colour = "pink") + geom_histogram(data = data.table(x = testStats), aes( y = ..density..), binwidth = 1, fill = "grey", color = "black") + stat_function(fun = actuar::dphtype, args = list(prob = c(c(uniformProportion, 1-uniformProportion), rep(0,2*(randomSetSize-1))), rates = as.matrix(hyperexponentialSumTransitionMatrix)), colour = "green", size = 1) + ylab("") + scale_y_continuous(breaks = NULL) + geom_vline(xintercept = testStats[1]) empiricalVsPred <- rbindlist(nullModelTests) empiricalVsPred[empiricalGTE <= 1E-2] empiricalVsPred[,qplot(x = combinedPval, y = (empiricalGTE - combinedPval))] + scale_x_log10(breaks = 10^(-seq(0,6))) + scale_y_continuous(breaks = seq(-0.13,0.1,0.01)) + theme_bw() empiricalVsPred[,qplot(x = combinedPval, y = (empiricalGTE - approxPval))] + scale_x_log10(breaks = 10^(-seq(0,6))) + scale_y_continuous(breaks = seq(-0.13,0.1,0.01)) + theme_bw() empiricalVsPred[,qplot(x = combinedPval, y = (empiricalGTE - phaseRateP))] + scale_x_log10(breaks = 10^(-seq(0,6))) + scale_y_continuous(breaks = seq(-0.13,0.1,0.01)) + theme_bw() empiricalVsPred[,qplot(x = size, y = (empiricalGTE - phaseRateP))] empiricalVsPred[empiricalGTE <= 1E-2] %>% ggplot(aes(x = empiricalGTE)) + geom_point(aes(y = fisherPval), colour = "blue") + geom_point(aes(y = combinedPval), colour = "red") + geom_point(aes(y = approxPval), colour = "green") + geom_point(aes(y = phaseRateP), colour = "pink") + geom_abline() + theme_bw() + scale_x_log10() + scale_y_log10()
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.