Skip to content

Commit

Permalink
minor speedups for GEE and GLM modes
Browse files Browse the repository at this point in the history
  • Loading branch information
@jr-leary7
jr-leary7 committed Oct 10, 2024
1 parent 75c7d93 commit cff92ec
Show file tree
Hide file tree
Showing 6 changed files with 41 additions and 38 deletions.
1 change: 1 addition & 0 deletions NAMESPACE
View file
Original file line number Diff line number Diff line change
Expand Up @@ -44,6 +44,7 @@ importFrom(Matrix,rowMeans)
importFrom(Matrix,solve)
importFrom(Matrix,t)
importFrom(Rcpp,sourceCpp)
importFrom(RcppEigen,fastLmPure)
importFrom(bigstatsr,as_FBM)
importFrom(broom.mixed,tidy)
importFrom(doSNOW,registerDoSNOW)
Expand Down
8 changes: 3 additions & 5 deletions R/backward_sel_WIC.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@
#' @importFrom gamlss gamlss
#' @importFrom geeM geem
#' @importFrom MASS negative.binomial
#' @importFrom withr with_output_sink
#' @importFrom stats vcov
#' @param Y The response variable. Defaults to NULL.
#' @param B_new The model matrix. Defaults to NULL.
#' @param is.gee Is the model a GEE? Defaults to FALSE.
Expand Down Expand Up @@ -43,10 +43,8 @@ backward_sel_WIC <- function(Y = NULL,
fit <- gamlss::gamlss(Y ~ B_new - 1,
family = "NBI",
trace = FALSE)
withr::with_output_sink(tempfile(), {
fit_sum_mat <- as.matrix(summary(fit))
})
wald_stat <- fit_sum_mat[, 3][-c(1, nrow(fit_sum_mat))]^2
vcov_mat <- stats::vcov(fit, type = "all")
wald_stat <- unname((vcov_mat$coef / vcov_mat$se)[-c(1, length(vcov_mat$coef))]^2)
}
return(wald_stat)
}
2 changes: 1 addition & 1 deletion R/marge2.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -610,7 +610,7 @@ marge2 <- function(X_pred = NULL,

if (trunc.type == 2) {
B2a <- matrix(rep(B2[, best.var], 2), ncol = 2)
B_temp <- cbind(B, B2a*cbind(b1_new, b2_new))
B_temp <- cbind(B, B2a * cbind(b1_new, b2_new))
B_new <- B2a * cbind(b1_new, b2_new)
var_name3 <- var_name2[best.var]
colnames(B_new) <- rep(var_name3, 2)
Expand Down
20 changes: 10 additions & 10 deletions R/score_fun_gee.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,8 +4,8 @@
#' @author Jakub Stoklosa
#' @author David I. Warton
#' @author Jack Leary
#' @importFrom stats lm.fit
#' @importFrom Matrix diag t solve
#' @importFrom RcppEigen fastLmPure
#' @importFrom Matrix solve
#' @importFrom MASS ginv
#' @description Calculate the score statistic for a GEE model.
#' @param Y The response variable. Defaults to NULL.
Expand Down Expand Up @@ -42,7 +42,7 @@ score_fun_gee <- function(Y = NULL,
# check inputs
if (is.null(Y) || is.null(N) || is.null(n_vec) || is.null(VS.est_list) || is.null(AWA.est_list) || is.null(J2_list) || is.null(Sigma2_list) || is.null(J11.inv) || is.null(JSigma11) || is.null(mu.est) || is.null(V.est) || is.null(B1) || is.null(XA)) { stop("Some inputs to score_fun_gee() are missing.") }
# generate score statistic
reg <- try({ stats::lm.fit(B1, Y) }, silent = TRUE) # This is not the model fit!! It just checks whether any issues occur for a simple linear regression model.
reg <- try({ RcppEigen::fastLmPure(B1, Y) }, silent = TRUE) # This is not the model fit!! It just checks whether any issues occur for a simple linear regression model.
if (inherits(reg, "try-error") || any(is.na(reg$coefficients))) {
score <- NA_real_
} else {
Expand All @@ -58,15 +58,15 @@ score_fun_gee <- function(Y = NULL,
AWA.est_i <- AWA.est_list[[i]]
J2_i <- J2_list[[i]]
Sigma2_i <- Sigma2_list[[i]]
D.est_i <- eigenMapMatMult(A = Matrix::diag((mu.est[(sum(n_vec1[seq(i)]) + 1):k]), nrow = n_vec[i], ncol = n_vec[i]),
D.est_i <- eigenMapMatMult(A = diag(mu.est[(sum(n_vec1[seq(i)]) + 1):k], nrow = n_vec[i], ncol = n_vec[i]),
B = XA[(sum(n_vec1[seq(i)]) + 1):k, ],
n_cores = 1)
D_est_i_transpose <- Matrix::t(D.est_i)
D_est_i_transpose <- t(D.est_i)
J21 <- J21 + eigenMapMatMult(A = D_est_i_transpose,
B = Matrix::t(J2_i),
B = t(J2_i),
n_cores = 1)
Sigma21 <- Sigma21 + eigenMapMatMult(A = D_est_i_transpose,
B = Matrix::t(Sigma2_i),
B = t(Sigma2_i),
n_cores = 1)
B.est <- B.est + eigenMapMatMult(A = D_est_i_transpose,
B = VS.est_i,
Expand All @@ -82,12 +82,12 @@ score_fun_gee <- function(Y = NULL,
B = J11.inv,
n_cores = 1)
temp_prod_1 <- eigenMapMatMult(A = temp_prod_1,
B = Matrix::t(Sigma21),
B = t(Sigma21),
n_cores = 1)
temp_prod_2 <- eigenMapMatMult(A = Sigma21,
B = J11.inv,
n_cores = 1)
J21_transpose <- Matrix::t(J21)
J21_transpose <- t(J21)
temp_prod_2 <- eigenMapMatMult(A = temp_prod_2,
B = J21_transpose,
n_cores = 1)
Expand All @@ -102,7 +102,7 @@ score_fun_gee <- function(Y = NULL,
if (inherits(Sigma_inv, "try-error")) {
Sigma_inv <- MASS::ginv(Sigma)
}
temp_prod <- eigenMapMatMult(A = Matrix::t(B.est),
temp_prod <- eigenMapMatMult(A = t(B.est),
B = Sigma_inv,
n_cores = 1)
score <- eigenMapMatMult(A = temp_prod,
Expand Down
16 changes: 8 additions & 8 deletions R/score_fun_glm.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,8 +4,8 @@
#' @author Jakub Stoklosa
#' @author David I. Warton
#' @author Jack Leary
#' @importFrom stats lm.fit
#' @importFrom Matrix t
#' @importFrom RcppEigen fastLmPure
#' @importFrom Matrix solve
#' @importFrom MASS ginv
#' @description Calculate the score statistic for a GLM model.
#' @param Y The response variable. Defaults to NULL.
Expand All @@ -31,7 +31,7 @@ score_fun_glm <- function(Y = NULL,
# check inputs
if (is.null(Y) || is.null(VS.est_list) || is.null(A_list) || is.null(B1_list) || is.null(mu.est) || is.null(V.est) || is.null(B1) || is.null(XA)) { stop("Some inputs to score_fun_glm() are missing.") }
# generate score statistic
reg <- try({ stats::lm.fit(B1, Y) }, silent = TRUE) # This is not the model fit!! It just checks whether any issues occur for a simple linear regression model.
reg <- try({ RcppEigen::fastLmPure(B1, Y) }, silent = TRUE) # This is not the model fit!! It just checks whether any issues occur for a simple linear regression model.
if (inherits(reg, "try-error") || any(is.na(reg$coefficients))) {
score <- NA_real_
} else {
Expand All @@ -41,28 +41,28 @@ score_fun_glm <- function(Y = NULL,
B_list_i <- eigenMapMatMult(A = B1_list_i,
B = XA,
n_cores = 1)
D_list_i <- eigenMapMatMult(A = Matrix::t(XA),
D_list_i <- eigenMapMatMult(A = t(XA),
B = (XA * c(mu.est^2 / V.est)),
n_cores = 1)
temp_prod <- eigenMapMatMult(A = Matrix::t(B_list_i),
temp_prod <- eigenMapMatMult(A = t(B_list_i),
B = A_list_i,
n_cores = 1)
temp_prod <- eigenMapMatMult(A = temp_prod,
B = B_list_i,
n_cores = 1)
inv.XVX_22 <- D_list_i - temp_prod
B.est <- eigenMapMatMult(A = Matrix::t(mu.est * VS.est_i),
B.est <- eigenMapMatMult(A = t(mu.est * VS.est_i),
B = XA,
n_cores = 1)
XVX_22 <- try({ solve(inv.XVX_22) }, silent = TRUE)
XVX_22 <- try({ Matrix::solve(inv.XVX_22) }, silent = TRUE)
if (inherits(XVX_22, "try-error")) {
XVX_22 <- MASS::ginv(inv.XVX_22)
}
temp_prod <- eigenMapMatMult(A = B.est,
B = XVX_22,
n_cores = 1)
score <- eigenMapMatMult(A = temp_prod,
B = Matrix::t(B.est),
B = t(B.est),
n_cores = 1)
}
res <- list(score = score)
Expand Down
32 changes: 18 additions & 14 deletions R/stat_out_score_gee_null.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@
#' @importFrom geeM geem
#' @importFrom MASS negative.binomial ginv
#' @importFrom stats fitted.values
#' @importFrom Matrix diag t solve
#' @importFrom Matrix solve
#' @description A function that calculates parts of the score statistic for GEEs only (it is used for the full path for forward selection).
#' @param Y The response variable Defaults to NULL.
#' @param B_null The design matrix matrix under the null model Defaults to NULL.
Expand Down Expand Up @@ -38,7 +38,8 @@ stat_out_score_gee_null <- function(Y = NULL,
corstr = cor.structure,
family = MASS::negative.binomial(theta.hat),
scale.fix = FALSE,
sandwich = FALSE)
sandwich = FALSE,
maxit = 10)
alpha_est <- ests$alpha
sigma_est <- ests$phi
mu_est <- as.matrix(stats::fitted.values(ests))
Expand All @@ -51,10 +52,10 @@ stat_out_score_gee_null <- function(Y = NULL,
k <- sum(n_vec[seq(i)])
# set up working correlation matrix structure
if (cor.structure == "independence") {
R_alpha <- Matrix::diag(1, n_vec[i], n_vec[i])
R_alpha <- diag(1, n_vec[i], n_vec[i])
} else if (cor.structure == "exchangeable") {
R_alpha <- matrix(alpha_est, nrow = n_vec[i], ncol = n_vec[i])
R_alpha <- R_alpha + Matrix::diag(c(1 - alpha_est), nrow = n_vec[i], ncol = n_vec[i])
R_alpha <- R_alpha + diag(1 - alpha_est, nrow = n_vec[i], ncol = n_vec[i])
} else if (cor.structure == "ar1") {
R_alpha <- alpha_est^outer(seq(n_vec[i]), seq(n_vec[i]), \(x, y) abs(x - y))
} else {
Expand All @@ -63,32 +64,35 @@ stat_out_score_gee_null <- function(Y = NULL,
# compute iteration values for each statistic
temp_seq_n <- (sum(n_vec1[seq(i)]) + 1):k
mu_est_sub <- mu_est[temp_seq_n]
diag_sqrt_V_est <- Matrix::diag(sqrt(V_est[temp_seq_n]),
nrow = n_vec[i],
ncol = n_vec[i])
diag_sqrt_V_est <- diag(sqrt(V_est[temp_seq_n]),
nrow = n_vec[i],
ncol = n_vec[i])
temp_prod <- eigenMapMatMult(A = diag_sqrt_V_est,
B = R_alpha,
n_cores = 1)
V_est_i <- eigenMapMatMult(A = temp_prod,
B = diag_sqrt_V_est,
n_cores = 1)
V_est_i_inv <- eigenMapMatrixInvert(A = V_est_i, n_cores = 1)
S_est_i <- Matrix::t(Y)[temp_seq_n] - mu_est_sub
V_est_i_inv <- try({ Matrix::solve(V_est_i) }, silent = TRUE)
if (inherits(V_est_i_inv, "try-error")) {
V_est_i_inv <- MASS::ginv(V_est_i)
}
S_est_i <- t(Y)[temp_seq_n] - mu_est_sub
temp_prod <- eigenMapMatMult(A = S_est_i,
B = Matrix::t(S_est_i),
B = t(S_est_i),
n_cores = 1)
temp_prod <- eigenMapMatMult(A = V_est_i_inv,
B = temp_prod,
n_cores = 1)
AWA_est_i <- eigenMapMatMult(A = temp_prod,
B = V_est_i_inv,
n_cores = 1)
D_est_i <- eigenMapMatMult(A = Matrix::diag(mu_est_sub,
nrow = n_vec[i],
ncol = n_vec[i]),
D_est_i <- eigenMapMatMult(A = diag(mu_est_sub,
nrow = n_vec[i],
ncol = n_vec[i]),
B = B_null[temp_seq_n, ],
n_cores = 1)
D_est_i_transpose <- Matrix::t(D_est_i)
D_est_i_transpose <- t(D_est_i)
temp_prod <- eigenMapMatMult(A = D_est_i_transpose,
B = V_est_i_inv,
n_cores = 1)
Expand Down

0 comments on commit cff92ec

Please sign in to comment.