From 2ecb9b92dcb8771d309546709cc2c5e5c706c39e Mon Sep 17 00:00:00 2001 From: Xihao Li Date: Mon, 9 Sep 2024 15:34:18 -0400 Subject: [PATCH] STAAR v0.9.7.1 --- R/fit_null_glmmkin_Binary_SPA.R | 50 +++++++++++++++------------------ 1 file changed, 22 insertions(+), 28 deletions(-) diff --git a/R/fit_null_glmmkin_Binary_SPA.R b/R/fit_null_glmmkin_Binary_SPA.R index 30c53a6..a80b77f 100644 --- a/R/fit_null_glmmkin_Binary_SPA.R +++ b/R/fit_null_glmmkin_Binary_SPA.R @@ -80,22 +80,18 @@ fit_null_glmmkin_Binary_SPA <- function(fixed, data = parent.frame(), kins, use_ obj_nullmodel$sparse_kins <- TRUE ## generate W - muhat <- obj_nullmodel$fitted.values - working <- muhat*(1-muhat) + X <- obj_nullmodel$X + muhat <- obj_nullmodel$fitted.values + working <- muhat*(1-muhat) - X <- obj_nullmodel$X - - obj_nullmodel$XW <- t(X)%*%diag(working) - obj_nullmodel$XXWX_inv <- X%*%solve(t(X)%*%diag(working)%*%X) - - ## generate Sigma_i - obj_nullmodel$XSigma_i <- crossprod(X,obj_nullmodel$Sigma_i) - obj_nullmodel$XXSigma_iX_inv <- X%*%obj_nullmodel$cov + obj_nullmodel$XW <- t(X*working) + obj_nullmodel$XXWX_inv <- X%*%solve(t(X*working)%*%X) + ## generate Sigma_i + obj_nullmodel$XSigma_i <- crossprod(X,obj_nullmodel$Sigma_i) + obj_nullmodel$XXSigma_iX_inv <- X%*%obj_nullmodel$cov }else if(!is.null(use_sparse) && use_sparse){ print(paste0("kins is a dense matrix, transforming it into a sparse matrix using cutoff ", kins_cutoff,".")) - #kins <- replace(kins, kins <= kins_cutoff, 0) - #kins_sp <- Matrix(kins, sparse = TRUE) kins_sp <- makeSparseMatrix(kins, thresh = kins_cutoff) if(inherits(kins_sp, "dsyMatrix") || kins_cutoff <= min(kins)){ stop(paste0("kins is still a dense matrix using cutoff ", kins_cutoff,". Please try a larger kins_cutoff or use_sparse = FALSE!")) @@ -110,17 +106,16 @@ fit_null_glmmkin_Binary_SPA <- function(fixed, data = parent.frame(), kins, use_ obj_nullmodel$sparse_kins <- TRUE ## generate W - muhat <- obj_nullmodel$fitted.values - working <- muhat*(1-muhat) - - X <- obj_nullmodel$X - obj_nullmodel$XW <- t(X)%*%diag(working) - obj_nullmodel$XXWX_inv <- X%*%solve(t(X)%*%diag(working)%*%X) + X <- obj_nullmodel$X + muhat <- obj_nullmodel$fitted.values + working <- muhat*(1-muhat) - ## generate Sigma_i - obj_nullmodel$XSigma_i <- crossprod(X,obj_nullmodel$Sigma_i) - obj_nullmodel$XXSigma_iX_inv <- X%*%obj_nullmodel$cov + obj_nullmodel$XW <- t(X*working) + obj_nullmodel$XXWX_inv <- X%*%solve(t(X*working)%*%X) + ## generate Sigma_i + obj_nullmodel$XSigma_i <- crossprod(X,obj_nullmodel$Sigma_i) + obj_nullmodel$XXSigma_iX_inv <- X%*%obj_nullmodel$cov }else{ print("kins is a dense matrix.") obj_nullmodel <- glmmkin(fixed = fixed, data = data, kins = kins, id = id, @@ -131,14 +126,13 @@ fit_null_glmmkin_Binary_SPA <- function(fixed, data = parent.frame(), kins, use_ tauregion = tauregion, verbose = verbose, ...) obj_nullmodel$sparse_kins <- FALSE - ## generate W - muhat <- obj_nullmodel$fitted.values - working <- muhat*(1-muhat) - - X <- obj_nullmodel$X + ## generate W + X <- obj_nullmodel$X + muhat <- obj_nullmodel$fitted.values + working <- muhat*(1-muhat) - obj_nullmodel$XW <- t(X)%*%diag(working) - obj_nullmodel$XXWX_inv <- X%*%solve(t(X)%*%diag(working)%*%X) + obj_nullmodel$XW <- t(X*working) + obj_nullmodel$XXWX_inv <- X%*%solve(t(X*working)%*%X) } obj_nullmodel$relatedness <- TRUE obj_nullmodel$use_SPA <- TRUE