.Rhistory

} else {
# check whether the covariates are positive for the box-cox--------
if(any(covariate <= 0)){
stop("covariate statistics need to be positive")
}
bc <- boxcox(covariate ~ test)
lambda <- bc$x[which.max(bc$y)]
if(lambda == 0){
model <- lm(log(covariate + .0001) ~ test)
} else {
model <- lm(covariate**lambda ~ test)
}
mean_covariateEffect <- model$coef[[1]] + model$coef[[2]]*mean_testEffect
}
message("computing rank probabilities")
# compute the probability of the rank of the covariate given the mean effect
ranksProb <- sapply(1:m, prob_rank_givenEffect, et = mean_covariateEffect,
ey = mean_covariateEffect, nrep = nrep, m0 = m0, m1 = m1)
message("finished computing the rank probabilities")
}
# compute the weights (always right-tailed)------------
message("computing weights")
if(effectType == "continuous"){
wgt = weight_continuous_nwt(alpha = alpha, et = mean_testEffect,
ranksProb = ranksProb)$w
# wgt = weight_continuous(alpha = alpha, et = mean_testEffect, m = m,
#                         tail = 1, delInterval = delInterval,
#                         ranksProb = ranksProb)
} else {
wgt = weight_binary_nwt(alpha = alpha, m1 = m1,
et = mean_testEffect, ranksProb = ranksProb)
# wgt = weight_binary(alpha = alpha, et = mean_testEffect, m = m,
#                     m1 = m1, tail = 1, delInterval = delInterval,
#                     ranksProb = ranksProb)
}
message("finished computing the weights")
}
message("comparing pvalues with thresholds")
if(method == "BH"){
padj <- p.adjust(Ordered.pvalue/wgt, method = "BH")
OD <- add_column(OD, padj)
rejections_list = OD[which((padj <= alpha) == TRUE), ]
} else {
rejections_list = OD[which((Ordered.pvalue <= alpha*wgt/m) == TRUE), ]
}
# outputs--------------
n_rejections = dim(rejections_list)[1]
return(list(totalTests = m, nullProp = nullProp,
ranksProb = ranksProb, weight = wgt,
rejections = n_rejections, rejections_list = rejections_list))
}
m = 1000
set.seed(3)
covariates = runif(m, min = 0, max = 2.5)          # covariate statistics
H = rbinom(m, size = 1, prob = 0.1)             # hypothesis true or false
tests = rnorm(m, mean = H * covariates)            # Z-score
pvals = 1 - pnorm(tests)                        # pvalue
# general use
results <- opw(pvalue = pvals, covariate = covariates, effectType = "continuous",
method = "BH")
library(qvalue)
results <- opw(pvalue = pvals, covariate = covariates, effectType = "continuous",
method = "BH")
library(tibble)
results <- opw(pvalue = pvals, covariate = covariates, effectType = "continuous",
method = "BH")
library(MASS)
results <- opw(pvalue = pvals, covariate = covariates, effectType = "continuous",
method = "BH")
results
mod <- lm(log(covariates) ~ tests)
et = mean(tests)
ey = mod$coef[[1]] + mod$coef[[2]]*et
results2 <- opw(pvalue = pvals, covariate = covariates,
mean_covariateEffect = ey, mean_testEffect = et, tail = 2,
effectType = "continuous", method = "BH")
m = length(pvalue)
pvalue = pvals
covariate = covariates
mean_covariateEffect = ey
mean_testEffect = et
tail = 2
effectType = "continuous"
method = "BH"
m = length(pvalue)
nullProp = qvalue(p = pvalue, pi0.method = "bootstrap")$pi0
m0 = ceiling(nullProp*m)
m1 = m - m0
Data = tibble(pvalue, covariate)
OD <- Data[order(Data$covariate, decreasing=TRUE), ]
Ordered.pvalue <- OD$pvalue
test <- qnorm(pvalue/tail, lower.tail = FALSE)
test[which(!is.finite(test))] <- NA
test_effect_vec <-  sort(test, decreasing = TRUE)[1:m1]
mean_testEffect <- mean(test_effect_vec, na.rm = TRUE)
mean_covariateEffect <- mean_covariateEffect
bc <- boxcox(covariate ~ test)
lambda <- bc$x[which.max(bc$y)]
if(lambda == 0){
model <- lm(log(covariate + .0001) ~ test)
} else {
model <- lm(covariate**lambda ~ test)
}
mean_covariateEffect <- model$coef[[1]] + model$coef[[2]]*mean_testEffect
ranksProb <- sapply(1:m, prob_rank_givenEffect, et = mean_covariateEffect,
ey = mean_covariateEffect, nrep = nrep, m0 = m0, m1 = m1)
nrep=10000
ranksProb <- sapply(1:m, prob_rank_givenEffect, et = mean_covariateEffect,
ey = mean_covariateEffect, nrep = nrep, m0 = m0, m1 = m1)
ranksProb
alpha = alpha
alpha=.05
et = mean_testEffect
ranksProb = ranksProb
m = length(ranksProb)
prob <- ranksProb/sum(ranksProb, na.rm = TRUE)
f <- function(c)
{
sum(pnorm(et/2 + c/et - log(alpha*prob)/et, lower.tail = FALSE)) - alpha
}
df <- function(c)
{
sum(-dnorm(et/2 + c/et - log(alpha*prob)/et)/et)
}
nmax <- 1000
n = 1
if(is.null(x0)) {x0 = 0} else {x0 = x0}
while(f(x0) < 0 && (n <= nmax))
{
if(f(0) > f(.5)){x0 = x0 - .5} else {x0 = x0 + .5}
n = n + 1
}
x0 = NULL
if(is.null(x0)) {x0 = 0} else {x0 = x0}
while(f(x0) < 0 && (n <= nmax))
{
if(f(0) > f(.5)){x0 = x0 - .5} else {x0 = x0 + .5}
n = n + 1
}
if(f(x0) < 0 && (n >= nmax)) {
cat("f(x0) = ", f(x0))
stop("f(x0) is negative, modify initial guess x0 so that f(x0) >= 0")
}
x <- rep(0, nmax)
ex <- rep(0, nmax)
x[1] <- x0 - f(x0) / df(x0)
ex[1] <- abs((x[1] - x0)/x[1])
k <- 2
while ((ex[k - 1] >= .0001) && (k <= nmax))
{
x[k] <- x[k - 1] - (f(x[k - 1]) / df(x[k - 1]))
ex[k] <- abs((x[k] - x[k - 1])/x[k])
k <- k + 1
}
c <- x[k - 1] # could return this
w <- (m/alpha)*pnorm(et/2 + c/et - log(alpha*prob)/et, lower.tail = FALSE)
lambda <- exp(c)
library(qvalue)
ranks <- 1:m
nullProp = qvalue(p = pvals, pi0.method = "bootstrap")$pi0
m0 = ceiling(nullProp*m)
m1 = m - m0
probs <- sapply(ranks, prob_rank_givenEffect, et = ey, ey = ey,
nrep = 10000, m0 = m0, m1 = m1)
results3 <- opw(pvalue = pvals, covariate = covariates, ranksProb = probs,
effectType = "continuous", tail = 2, method = "BH")
# supply weight externally
wgt <- weight_continuous(alpha = .05, et = et, m = m, ranksProb = probs)
results4 <- opw(pvalue = pvals, covariate = covariates, weight = wgt,
effectType = "continuous", alpha = .05, method = "BH")
results = weight_continuous_nwt(alpha = .05, et = et, ranksProb = prob)
results
results2 <- opw(pvalue = pvals, covariate = covariates,
mean_covariateEffect = ey, mean_testEffect = et, tail = 2,
effectType = "continuous", method = "BH")
m = 1000
set.seed(3)
covariates = runif(m, min = 0, max = 2.5)          # covariate statistics
H = rbinom(m, size = 1, prob = 0.1)             # hypothesis true or false
tests = rnorm(m, mean = H * covariates)            # Z-score
pvals = 1 - pnorm(tests)                        # pvalue
mod <- lm(log(covariates) ~ tests)
et = mean(tests)
ey = mod$coef[[1]] + mod$coef[[2]]*et
results2 <- opw(pvalue = pvals, covariate = covariates,
mean_covariateEffect = ey, mean_testEffect = et, tail = 2,
effectType = "continuous", method = "BH")
rm(list=ls())
m = 1000
set.seed(3)
covariates = runif(m, min = 0, max = 2.5)          # covariate statistics
H = rbinom(m, size = 1, prob = 0.1)             # hypothesis true or false
tests = rnorm(m, mean = H * covariates)            # Z-score
pvals = 1 - pnorm(tests)                        # pvalue
results <- opw(pvalue = pvals, covariate = covariates, effectType = "continuous",
method = "BH")
results
mod <- lm(log(covariates) ~ tests)
et = mean(tests)
ey = mod$coef[[1]] + mod$coef[[2]]*et
results2 <- opw(pvalue = pvals, covariate = covariates,
mean_covariateEffect = ey, mean_testEffect = et, tail = 2,
effectType = "continuous", method = "BH")
ey
et
m = 1000
set.seed(3)
covariates = runif(m, min = 1, max = 3)          # covariate statistics
H = rbinom(m, size = 1, prob = 0.1)             # hypothesis true or false
tests = rnorm(m, mean = H * covariates)            # Z-score
pvals = 1 - pnorm(tests)                        # pvalue
# general use
results <- opw(pvalue = pvals, covariate = covariates, effectType = "continuous",
method = "BH")
mod <- lm(log(covariates) ~ tests)
et = mean(tests)
ey = mod$coef[[1]] + mod$coef[[2]]*et
results2 <- opw(pvalue = pvals, covariate = covariates,
mean_covariateEffect = ey, mean_testEffect = et, tail = 2,
effectType = "continuous", method = "BH")
et
ey
library(qvalue)
ranks <- 1:m
nullProp = qvalue(p = pvals, pi0.method = "bootstrap")$pi0
m0 = ceiling(nullProp*m)
m1 = m - m0
mod <- lm(log(covariates) ~ tests)
et = mean(sort(tests, decreasing = TRUE)[1:m1], na.rm = TRUE)
et
ey = mod$coef[[1]] + mod$coef[[2]]*et
ey
results2 <- opw(pvalue = pvals, covariate = covariates,
mean_covariateEffect = ey, mean_testEffect = et, tail = 2,
effectType = "continuous", method = "BH")
probs <- sapply(ranks, prob_rank_givenEffect, et = ey, ey = ey,
nrep = 10000, m0 = m0, m1 = m1)
results3 <- opw(pvalue = pvals, covariate = covariates, ranksProb = probs,
effectType = "continuous", tail = 2, method = "BH")
# supply weight externally
wgt <- weight_continuous(alpha = .05, et = et, m = m, ranksProb = probs)
results4 <- opw(pvalue = pvals, covariate = covariates, weight = wgt,
effectType = "continuous", alpha = .05, method = "BH")
prob_cont
knitr::opts_chunk$set(tidy = FALSE, cache = FALSE, autodep = TRUE)
# compute weihgts
wgt = weight_continuous_nwt(alpha = .1, et = mean_testEffect, ranksProb = prob_cont)
set.seed(123)
# compute ranks probability
prob_cont <- sapply(1:m, prob_rank_givenEffect, et = mean_covariateEffect,
ey = mean_covariateEffect, nrep = 1000, m0 = m0, m1 = m1)
# initial stage--------
pvals = de_res_air$pvalue
# install OPWeight package
source("https://bioconductor.org/biocLite.R")
biocLite("OPWeight", suppressUpdates=TRUE)
# load packages
library("OPWeight")
library("ggplot2")
library("airway")
library("DESeq2")
library(cowplot)
library(tibble)
library(qvalue)
library(MASS)
data("airway", package = "airway")
dds <- DESeqDataSet(se = airway, design = ~ cell + dex)
dds <- DESeq(dds)
de_res_air <- results(dds)
colnames(de_res_air)
dim(de_res_air)
covariates = de_res_air$baseMean + .0001  # add a small constant to make all values positive
set.seed(123)
# one should use more nrep to obtain more acurate results
opw_results <- opw(pvalue = de_res_air$pvalue, covariate = covariates, nrep = 1000,
alpha = .1, tail = 2, effectType = "continuous", method = "BH")
# initial stage--------
pvals = de_res_air$pvalue
tests = qnorm(pvals/2, lower.tail = FALSE)
covariates = de_res_air$baseMean + .0001    # to ensure covariates are postive for the box-cox
# formulate a data set-------------
Data = tibble(pvals, covariates)
OD <- Data[order(Data$covariates, decreasing=TRUE), ]
Ordered.pvalue <- OD$pvals
# estimate the true null and alternative test sizes------
m = length(Data$pvals); m
nullProp = qvalue(Data$pvals, pi0.method="bootstrap")$pi0; nullProp
m0 = ceiling(nullProp*m); m0
m1 = m - m0; m1
# fit box-cox regression
#--------------------------------
bc <- boxcox(covariates ~ tests)
lambda <- bc$x[which.max(bc$y)]; lambda
model <- lm(covariates^lambda ~ tests)
# etimated test efects of the true altrnatives------------
test_effect = if(m1 == 0) {0
} else {sort(tests, decreasing = TRUE)[1:m1]}		# two-tailed test
# for the continuous effects etimated mean effects
mean_testEffect = mean(test_effect, na.rm = TRUE)
mean_testEffect
mean_covariateEffect = model$coef[[1]] + model$coef[[2]]*mean_testEffect
mean_covariateEffect
set.seed(123)
# compute ranks probability
prob_cont <- sapply(1:m, prob_rank_givenEffect, et = mean_covariateEffect,
ey = mean_covariateEffect, nrep = 1000, m0 = m0, m1 = m1)
# compute weihgts
wgt = weight_continuous_nwt(alpha = .1, et = mean_testEffect, ranksProb = prob_cont)
alpha = .1
padj <- p.adjust(Ordered.pvalue/wgt, method = "BH")
wgt=wgt
alpha = .1
padj <- p.adjust(Ordered.pvalue/wgt, method = "BH")
# compute weihgts
wgt = weight_continuous_nwt(alpha = .1, et = mean_testEffect, ranksProb = prob_cont)
alpha = .1
padj <- p.adjust(Ordered.pvalue/wgt, method = "BH")
wgt
# compute weihgts
wgt = weight_continuous_nwt(alpha = .1, et = mean_testEffect, ranksProb = prob_cont)$w
alpha = .1
padj <- p.adjust(Ordered.pvalue/wgt, method = "BH")
rejections_list = OD[which((padj <= alpha) == TRUE), ]
n_rejections = dim(rejections_list)[1]
n_rejections
opw_results2 <- opw(pvalue = pvals, covariate = covariates, weight = wgt,
effectType = "continuous", alpha = .1, method = "BH")
opw_results2$rejections
library(OPWeight)
OPWeight::opw
opw
knitr::opts_chunk$set(tidy = FALSE, cache = FALSE, autodep = TRUE)
opw
# install OPWeight package
# if (!requireNamespace("BiocManager", quietly = TRUE))
#     install.packages("BiocManager")
# BiocManager::install("OPWeight")
# source("https://bioconductor.org/biocLite.R")
# biocLite("OPWeight", suppressUpdates=TRUE)
# load packages
library(tidyverse)
library(OPWeight)
library(ggplot2)
library(airway)
library(DESeq2)
library(dplyr)
library(rlang)
library(cowplot)
library(tibble)
library(qvalue)
library(MASS)
opw
data("airway", package = "airway")
dds <- DESeqDataSet(se = airway, design = ~ cell + dex)
dds <- DESeq(dds)
de_res_air <- as.data.frame(results(dds))
colnames(de_res_air)
dim(de_res_air)
# add a small constant to make all values positive
covariates = de_res_air$baseMean + .0001
set.seed(123)
# one should use more nrep to obtain more acurate results
opw_results <- opw(pvalue = de_res_air$pvalue, covariate = covariates, nrep = 1000,
alpha = .1, tail = 2, effectType = "continuous", method = "BH")
names(opw_results)
opw_results$rejections
m = opw_results$totalTests
testRanks = 1:m
probs = opw_results$ranksProb
weights = opw_results$weight
dat = tibble(testRanks, probs, weights)
names(opw_results)
opw
m = opw_results$totalTests
testRanks = 1:m
probs = opw_results$dataOut$ranksProb
weights = opw_results$dataOut$weight
dat = tibble(testRanks, probs, weights)
p_prob = ggplot(dat, aes(x = testRanks, y = probs)) + geom_point(size = .5, col="firebrick4") +
labs(x = "Test rank" , y = "p(rank | effect)")
p_weight = ggplot(dat, aes(x = testRanks, y = weights)) + geom_point(size = .5, col="firebrick4")+
labs(x = "Test rank" , y = "Weight")
plot_grid(p_prob, p_weight, labels = c("a", "b"))
Data <- tibble(pval=de_res_air$pvalue, covariate=de_res_air$baseMean)
hist_pval <- ggplot(Data, aes(x = Data$pval)) +
geom_histogram(colour = "#1F3552", fill = "#4281AE")+
labs(x = "P-values")
hist_covariate <- ggplot(Data, aes(x = Data$covariate)) +
geom_histogram( colour = "#1F3552", fill = "#4274AE")+
labs(x = "Covariate statistics")
pval_covariate <- ggplot(Data, aes(x = rank(-Data$covariate), y = -log10(pval))) +
geom_point()+
labs(x = "Ranks of covariates", y = "-log(pvalue)")
p_ecdf <- ggplot(Data, aes(x = pval)) +
stat_ecdf(geom = "step")+
labs(x = "P-values", title="Empirical cumulative distribution")+
theme(plot.title = element_text(size = rel(.7)))
p <- plot_grid(hist_pval, hist_covariate, pval_covariate, p_ecdf,
labels = c("a", "b", "c", "d"), ncol = 2)
# now add the title
title <- ggdraw() + draw_label("Airway: data summary")
plot_grid(title, p, ncol = 1, rel_heights=c(0.1, 1))
# initial stage--------
pvals = de_res_air$pvalue
tests = qnorm(pvals/2, lower.tail = FALSE)
covariates = de_res_air$baseMean + .0001    # to ensure covariates are postive for the box-cox
# formulate a data set-------------
Data = tibble(pvals, covariates)
OD <- Data[order(Data$covariates, decreasing=TRUE), ]
Ordered.pvalue <- OD$pvals
# estimate the true null and alternative test sizes------
m = length(Data$pvals); m
nullProp = qvalue(Data$pvals, pi0.method="bootstrap")$pi0; nullProp
m0 = ceiling(nullProp*m); m0
m1 = m - m0; m1
# fit box-cox regression
#--------------------------------
bc <- boxcox(covariates ~ tests)
lambda <- bc$x[which.max(bc$y)]; lambda
model <- lm(covariates^lambda ~ tests)
# etimated test efects of the true altrnatives------------
test_effect = if(m1 == 0) {0
} else {sort(tests, decreasing = TRUE)[1:m1]}		# two-tailed test
# for the continuous effects etimated mean effects
mean_testEffect = mean(test_effect, na.rm = TRUE)
mean_testEffect
mean_covariateEffect = model$coef[[1]] + model$coef[[2]]*mean_testEffect
mean_covariateEffect
# initial stage--------
pvals = de_res_air$pvalue
tests = qnorm(pvals/2, lower.tail = FALSE)
covariates = de_res_air$baseMean + .0001    # to ensure covariates are postive for the box-cox
# formulate a data set-------------
Data = tibble(pvals, covariates)
OD <- Data[order(Data$covariates, decreasing=TRUE), ]
Ordered.pvalue <- OD$pvals
# estimate the true null and alternative test sizes------
m = length(Data$pvals); m
nullProp = qvalue(Data$pvals, pi0.method="bootstrap")$pi0; nullProp
m0 = ceiling(nullProp*m); m0
m1 = m - m0; m1
# etimated test efects of the true altrnatives------------
# keep top m1 tests for regression ---------
test[which(!is.finite(test))] <- NA
# initial stage--------
pvals = de_res_air$pvalue
tests = qnorm(pvals/2, lower.tail = FALSE)
covariates = de_res_air$baseMean + .0001    # to ensure covariates are postive for the box-cox
# formulate a data set-------------
Data = tibble(pvals, covariates)
OD <- Data[order(Data$covariates, decreasing=TRUE), ]
Ordered.pvalue <- OD$pvals
# estimate the true null and alternative test sizes------
m = length(Data$pvals); m
nullProp = qvalue(Data$pvals, pi0.method="bootstrap")$pi0; nullProp
m0 = ceiling(nullProp*m); m0
m1 = m - m0; m1
# etimated test efects of the true altrnatives------------
# keep top m1 tests for regression ---------
test[which(!is.finite(tests))] <- NA
# initial stage--------
pvals = de_res_air$pvalue
tests = qnorm(pvals/2, lower.tail = FALSE)
covariates = de_res_air$baseMean + .0001    # to ensure covariates are postive for the box-cox
# formulate a data set-------------
Data = tibble(pvals, covariates)
OD <- Data[order(Data$covariates, decreasing=TRUE), ]
Ordered.pvalue <- OD$pvals
# estimate the true null and alternative test sizes------
m = length(Data$pvals); m
nullProp = qvalue(Data$pvals, pi0.method="bootstrap")$pi0; nullProp
m0 = ceiling(nullProp*m); m0
m1 = m - m0; m1
# etimated test efects of the true altrnatives------------
# keep top m1 tests for regression ---------
tests[which(!is.finite(tests))] <- NA
Data2 = add_column(Data, tests)
OD2 <- Data2[order(Data2$tests, decreasing=TRUE), ][1:m1, ]
test_effect_vec <-  OD2$tests
# fit box-cox regression
#--------------------------------
bc <- MASS::boxcox(covariates ~ tests, data = OD2)
lambda <- bc$x[which.max(bc$y)]
model <- lm(covariates**lambda ~ tests, data = OD2)
# for the continuous effects etimated mean effects
mean_testEffect = mean(test_effect_vec, na.rm = TRUE)
mean_testEffect
mean_covariateEffect = model$coef[[1]] + model$coef[[2]]*mean_testEffect
mean_covariateEffect
set.seed(123)
# compute ranks probability
prob_cont <- sapply(1:m, prob_rank_givenEffect, et = mean_covariateEffect,
ey = mean_covariateEffect, nrep = 1000, m0 = m0, m1 = m1)
# compute weihgts
wgt = weight_continuous_nwt(alpha = .1, et = mean_testEffect, ranksProb = prob_cont)$w
alpha = .1
padj <- p.adjust(Ordered.pvalue/wgt, method = "BH")
rejections_list = OD[which((padj <= alpha) == TRUE), ]
n_rejections = dim(rejections_list)[1]
n_rejections
opw_results2 <- opw(pvalue = pvals, covariate = covariates, weight = wgt,
effectType = "continuous", alpha = .1, method = "BH")
opw_results2$rejections
install.packages("tidyverse")
getwd()
purl("OPWeight.Rmd")
library(knitr)
purl("OPWeight.Rmd")
getwd()
setwd(C:/Users/MHASAN/Dropbox/My R packages/OPWeight/OPWeight/vignettes)
setwd("C:/Users/MHASAN/Dropbox/My R packages/OPWeight/OPWeight/vignettes")
purl("OPWeight.Rmd")