From c4ad85982903f8692c7b86b81b2b04f6d347a26e Mon Sep 17 00:00:00 2001 From: Adrienne Stilp Date: Wed, 20 Mar 2024 09:59:19 -0700 Subject: [PATCH] Switch any call to pchisq to .pchisq_filter_extreme --- R/admixMap.R | 4 ++-- R/jointScoreTest.R | 4 ++-- R/prepareOutputForNullModel.R | 4 ++-- R/testGeno.R | 10 +++++----- R/testSingleVarianceComponent.R | 2 +- R/variantSetTests.R | 12 ++++++------ 6 files changed, 18 insertions(+), 18 deletions(-) diff --git a/R/admixMap.R b/R/admixMap.R index 5e882665..65ccf16a 100644 --- a/R/admixMap.R +++ b/R/admixMap.R @@ -131,7 +131,7 @@ admixMap <- function(admixDataList, res[,"Est"] <- beta res[,"SE"] <- sqrt(Vbeta) res[,"Stat"] <- Stat - res[,"pval"] <- pchisq(Stat, df=1, lower.tail=FALSE) + res[,"pval"] <- .pchisq_filter_extreme(Stat, df=1, lower.tail=FALSE) }else{ Joint.Stat <- rep(NA, ncol(local)) @@ -167,7 +167,7 @@ admixMap <- function(admixDataList, res[,paste(names(admixDataList)[i],".SE", sep="")] <- SE[,i] } res[,"Joint.Stat"] <- Joint.Stat - res[,"Joint.pval"] <- pchisq(Joint.Stat, df=v, lower.tail=FALSE) + res[,"Joint.pval"] <- .pchisq_filter_extreme(Joint.Stat, df=v, lower.tail=FALSE) } # else # results data frame diff --git a/R/jointScoreTest.R b/R/jointScoreTest.R index dd910026..7ef984fb 100644 --- a/R/jointScoreTest.R +++ b/R/jointScoreTest.R @@ -51,7 +51,7 @@ jointScoreTest <- function(null.model, G) { # # Check rownames/colnames match. # Joint test statistic. Note this is the equivalent of the stat being squared, # by convention. Stat.joint <- as.numeric(crossprod(GY, beta)) - pval.joint <- pchisq(Stat.joint, lower.tail = FALSE, df = ncol(G)) + pval.joint <- .pchisq_filter_extreme(Stat.joint, lower.tail = FALSE, df = ncol(G)) # Percentage of variance explained jointly by these variants. pve.joint <- as.numeric(Stat.joint / null.model$RSS0) @@ -63,7 +63,7 @@ jointScoreTest <- function(null.model, G) { # # Check rownames/colnames match. stringsAsFactors = FALSE ) fixef$Stat <- fixef$Est / fixef$SE - fixef$pval <- pchisq(fixef$Stat^2, lower.tail = FALSE, df = 1) + fixef$pval <- .pchisq_filter_extreme(fixef$Stat^2, lower.tail = FALSE, df = 1) fixef$PVE <- fixef$Stat^2 / null.model$RSS0 rownames(fixef) <- colnames(G) diff --git a/R/prepareOutputForNullModel.R b/R/prepareOutputForNullModel.R index 62837b01..e4caad4a 100644 --- a/R/prepareOutputForNullModel.R +++ b/R/prepareOutputForNullModel.R @@ -80,7 +80,7 @@ SE <- sqrt(diag(betaCov)) Stat <- (vc.mod$beta/SE)^2 - pval <- pchisq(Stat, df = 1, lower.tail = FALSE) + pval <- .pchisq_filter_extreme(Stat, df = 1, lower.tail = FALSE) fixef <- data.frame(Est = vc.mod$beta, SE = SE, Stat = Stat, pval = pval) rownames(fixef) <- varNames @@ -181,7 +181,7 @@ SE <- sqrt(diag(betaCov)) Stat <- (vc.mod$beta/SE)^2 - pval <- pchisq(Stat, df = 1, lower.tail = FALSE) + pval <- .pchisq_filter_extreme(Stat, df = 1, lower.tail = FALSE) fixef <- data.frame(Est = vc.mod$beta, SE = SE, Stat = Stat, pval = pval) rownames(fixef) <- varNames diff --git a/R/testGeno.R b/R/testGeno.R index 9570b27a..e93a19bb 100644 --- a/R/testGeno.R +++ b/R/testGeno.R @@ -117,7 +117,7 @@ testGenoSingleVar <- function(nullmod, G, E = NULL, test = c("Score", "Score.SPA Stat <- score/score.SE res <- data.frame(Score = score, Score.SE = score.SE, Score.Stat = Stat, - Score.pval = pchisq(Stat^2, df = 1, lower.tail = FALSE), + Score.pval = .pchisq_filter_extreme(Stat^2, df = 1, lower.tail = FALSE), Est = score/GPG, Est.SE = 1/score.SE, PVE = (Stat^2)/RSS0) # RSS0 = (n-k) when gaussian; not when binary @@ -135,7 +135,7 @@ testGenoSingleVar <- function(nullmod, G, E = NULL, test = c("Score", "Score.SPA # Vbeta <- RSS/GPG # Stat <- beta/sqrt(Vbeta) # res <- data.frame(Est = beta, Est.SE = sqrt(Vbeta), Wald.Stat = Stat, -# Wald.pval = pchisq(Stat^2, df = 1, lower.tail = FALSE)) +# Wald.pval = .pchisq_filter_extreme(Stat^2, df = 1, lower.tail = FALSE)) # return(res) # } @@ -193,8 +193,8 @@ testGenoSingleVar <- function(nullmod, G, E = NULL, test = c("Score", "Score.SPA } res <- as.data.frame(res) - res$GxE.pval <- pchisq((res$GxE.Stat)^2, df = (v - 1), lower.tail = FALSE) - res$Joint.pval <- pchisq((res$Joint.Stat)^2, df = v, lower.tail = FALSE) + res$GxE.pval <- .pchisq_filter_extreme((res$GxE.Stat)^2, df = (v - 1), lower.tail = FALSE) + res$Joint.pval <- .pchisq_filter_extreme((res$Joint.Stat)^2, df = v, lower.tail = FALSE) if (GxE.return.cov.mat) { return(list(res = res, GxEcovMatList = res.Vbetas)) @@ -286,7 +286,7 @@ testGenoSingleVar <- function(nullmod, G, E = NULL, test = c("Score", "Score.SPA crossprod(GPG, betas))/RSS, error = function(e) { NA }) - res[,"Joint.pval"] <- pchisq(res[,"Joint.Stat"], df = v, lower.tail = FALSE) + res[,"Joint.pval"] <- .pchisq_filter_extreme(res[,"Joint.Stat"], df = v, lower.tail = FALSE) res <- as.data.frame(res) return(list(res = res, allelesCovMat = Vbetas)) diff --git a/R/testSingleVarianceComponent.R b/R/testSingleVarianceComponent.R index 537d4fdc..d0d28cd4 100644 --- a/R/testSingleVarianceComponent.R +++ b/R/testSingleVarianceComponent.R @@ -25,7 +25,7 @@ testSingleVarianceComponent <- function(nullmod, varCompName, covMatList, group. test.stat <- -2*(nullmod.noVarComp$logLik - nullmod$logLik) - pval <- 0.5*pchisq(test.stat, df = length(nullmod$varComp), lower.tail = FALSE) + 0.5*pchisq(test.stat, df = length(nullmod$varComp) - 1, lower.tail = F) + pval <- 0.5*.pchisq_filter_extreme(test.stat, df = length(nullmod$varComp), lower.tail = FALSE) + 0.5*.pchisq_filter_extreme(test.stat, df = length(nullmod$varComp) - 1, lower.tail = F) return(pval) } diff --git a/R/variantSetTests.R b/R/variantSetTests.R index 58f51d3a..e1a6f6ef 100644 --- a/R/variantSetTests.R +++ b/R/variantSetTests.R @@ -137,7 +137,7 @@ testVariantSet <- function( nullmod, G, weights, GG1 <- crossprod(G, G.rowSums) # WGPGW1 # O(mn) V.sum <- sum(GG1) # 1WGPGW1 burden.stat = U.sum / sqrt(V.sum) - burden.pval <- pchisq(burden.stat^2, df=1, lower.tail=FALSE) + burden.pval <- .pchisq_filter_extreme(burden.stat^2, df=1, lower.tail=FALSE) # adjust U and G for burden U <- U - GG1*U.sum/V.sum # WGPY - WGPGW1 * 1WGPY/(1WGPGW1) @@ -152,7 +152,7 @@ testVariantSet <- function( nullmod, G, weights, # # denominator for burden # V.sum <- sum(GG1) # 1WGPGW1 # # burden p-value - # burden.pval <- pchisq(U.sum^2/V.sum, df=1, lower.tail=FALSE) + # burden.pval <- .pchisq_filter_extreme(U.sum^2/V.sum, df=1, lower.tail=FALSE) # # adjust for burden # U <- U - GG1*U.sum/V.sum # V <- V - tcrossprod(GG1)/V.sum # O(m^2) @@ -163,7 +163,7 @@ testVariantSet <- function( nullmod, G, weights, err <- out$err # Fisher's method to combine p-values - smmat.pval <- tryCatch(pchisq(-2*log(burden.pval)-2*log(theta.pval), df=4, lower.tail = FALSE), error = function(e) { NA }) + smmat.pval <- tryCatch(.pchisq_filter_extreme(-2*log(burden.pval)-2*log(theta.pval), df=4, lower.tail = FALSE), error = function(e) { NA }) if(is.na(smmat.pval)) { err <- 1 smmat.pval <- burden.pval @@ -187,7 +187,7 @@ testVariantSet <- function( nullmod, G, weights, .regular <- function(Q, V, ncolG) { if(ncolG == 1){ - pv <- list(pval = pchisq(as.numeric(Q/V), df=1, lower.tail=FALSE), method = "integration") + pv <- list(pval = .pchisq_filter_extreme(as.numeric(Q/V), df=1, lower.tail=FALSE), method = "integration") }else{ lambda <- eigen(V, only.values = TRUE, symmetric=TRUE)$values @@ -381,7 +381,7 @@ testVariantSet <- function( nullmod, G, weights, # p value calculation if(length(scores) == 1){ lambdas[[i]] <- as.numeric(distMat) - pval <- pchisq(as.numeric(Q/distMat), df=1, lower.tail=FALSE) + pval <- .pchisq_filter_extreme(as.numeric(Q/distMat), df=1, lower.tail=FALSE) err <- ifelse(is.na(pval), 1, 0) }else{ lambda <- eigen(distMat, only.values = TRUE, symmetric=TRUE)$values @@ -493,7 +493,7 @@ integrateFxn <- function(x, qmin, otherParams, tau, rho){ temp.q<-(minval - mu)/sqrt(varia)*sqrt(2*degf) + degf - re<-pchisq(temp.q ,df=degf) * dchisq(x,df=1) + re<-.pchisq_filter_extreme(temp.q ,df=degf) * dchisq(x,df=1) return(re) }