Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Inv fix jk #108

Open
wants to merge 2 commits into
base: master
Choose a base branch
from
Open

Inv fix jk #108

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
ee659d3
use pseudoinversion when regular inversion fails
juhis Aug 26, 2019
2ed9116
an optional argument to remove binary covariates with a small count
juhis Aug 26, 2019
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
12 changes: 9 additions & 3 deletions R/SAIGE_fitGLMM_fast.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -602,6 +602,7 @@ solveSpMatrixUsingArma = function(sparseGRMtest){
#' @param isDiagofKinSetAsOne logical. Whether to set the diagnal elements in GRM to be 1. By default, FALSE
#' @param useSparseSigmaforInitTau logical. Whether to use sparse GRM to estimate the initial values for fitting the null GLMM. By default, FALSE
#' @param useSparseSigmaConditionerforPCG logical. Whether to use sparse GRM to construct a conditoner for PCG. By default, FALSE. Current this option is deactivated.
#' @param minCovariateCount integer. If binary covariates have a count less than this, they will be excluded from the model to avoid convergence issues. By default, -1 (no covariates will be excluded)
#' @return a file ended with .rda that contains the glmm model information, a file ended with .varianceRatio.txt that contains the variance ratio values, and a file ended with #markers.SPAOut.txt that contains the SPAGMMAT tests results for the markers used for estimating the variance ratio.
#' @export
fitNULLGLMM = function(plinkFile = "",
Expand Down Expand Up @@ -639,7 +640,8 @@ fitNULLGLMM = function(plinkFile = "",
isCovariateTransform = TRUE,
isDiagofKinSetAsOne = FALSE,
useSparseSigmaConditionerforPCG = FALSE,
useSparseSigmaforInitTau = FALSE){
useSparseSigmaforInitTau = FALSE,
minCovariateCount = -1){

useSparseSigmaConditionerforPCG = FALSE
if(useSparseSigmaConditionerforPCG){
Expand Down Expand Up @@ -832,7 +834,7 @@ fitNULLGLMM = function(plinkFile = "",

#check for perfect separation
if(traitType == "binary"){
out_checksep = checkPerfectSep(formula.null, data=dataMerge_sort)
out_checksep = checkPerfectSep(formula.null, data=dataMerge_sort, minCovariateCount)
covarColList <- covarColList[!(covarColList %in% out_checksep)]
formula = paste0(phenoCol,"~", paste0(covarColList,collapse="+"))
formula.null = as.formula(formula)
Expand Down Expand Up @@ -2322,7 +2324,7 @@ getsubGRM = function(sparseGRMFile=NULL,
}


checkPerfectSep<-function(formula, data){
checkPerfectSep<-function(formula, data, minCovariateCount){
X1<-model.matrix(formula,data=data)
X_name = colnames(X1)
X1 = as.matrix(X1[,-1])
Expand All @@ -2339,6 +2341,10 @@ checkPerfectSep<-function(formula, data){
colnamesDelete = c(colnamesDelete, X_name[i])
cat("perfect seperation is detected! ", X_name[i], " will be excluded in the model\n")
}
if(sum(sumTable < minCovariateCount) > 0){
colnamesDelete = c(colnamesDelete, X_name[i])
cat("less than ", minCovariateCount, " samples in a covariate detected! ", X_name[i], " will be excluded in the model\n")
}
}
}

Expand Down
8 changes: 7 additions & 1 deletion src/SAIGE_fitGLMM_fast.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2230,7 +2230,13 @@ Rcpp::List getCoefficients(arma::fvec& Yvec, arma::fmat& Xmat, arma::fvec& wVec,
}

arma::fmat Xmatt = Xmat.t();
arma::fmat cov = inv_sympd(Xmatt * Sigma_iX);
arma::fmat cov;
try {
cov = arma::inv_sympd(Xmatt * Sigma_iX);
} catch (const std::exception& e) {
cov = arma::pinv(Xmatt * Sigma_iX);
cout << "inv_sympd failed, inverted with pinv" << endl;
}
arma::fmat Sigma_iXt = Sigma_iX.t();
arma::fvec SigmaiXtY = Sigma_iXt * Yvec;
arma::fvec alpha = cov * SigmaiXtY;
Expand Down