Skip to content

Commit

Permalink
Switch any call to pchisq to .pchisq_filter_extreme
Browse files Browse the repository at this point in the history
  • Loading branch information
@amstilp
amstilp committed Mar 20, 2024
1 parent 503ba35 commit c4ad859
Show file tree
Hide file tree
Showing 6 changed files with 18 additions and 18 deletions.
4 changes: 2 additions & 2 deletions R/admixMap.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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))
Expand Down Expand Up @@ -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
Expand Down
4 changes: 2 additions & 2 deletions R/jointScoreTest.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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)
Expand Down
4 changes: 2 additions & 2 deletions R/prepareOutputForNullModel.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
10 changes: 5 additions & 5 deletions R/testGeno.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand All @@ -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)
# }

Expand Down Expand Up @@ -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))
Expand Down Expand Up @@ -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))
Expand Down
2 changes: 1 addition & 1 deletion R/testSingleVarianceComponent.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
}
12 changes: 6 additions & 6 deletions R/variantSetTests.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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)
Expand All @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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)
}

0 comments on commit c4ad859

Please sign in to comment.