From c1a01b9c1f97d55e54a6d1384532e23c15442cdc Mon Sep 17 00:00:00 2001 From: MikeDMorgan Date: Thu, 28 Aug 2025 18:46:58 +0100 Subject: [PATCH] fix bug in GLMM with multiple random effects including with a GRM --- R/glmm.R | 1 - R/testNhoods.R | 19 ++++++++++++++++--- src/fitGeneticPLGlmm.cpp | 2 ++ src/paramEst.cpp | 1 + src/pseudovarPartial.cpp | 9 +++++---- 5 files changed, 24 insertions(+), 8 deletions(-) diff --git a/R/glmm.R b/R/glmm.R index ac686fa..9313e3c 100644 --- a/R/glmm.R +++ b/R/glmm.R @@ -179,7 +179,6 @@ fitGLMM <- function(X, Z, y, offsets, init.theta=NULL, Kin=NULL, geno.I <- diag(nrow(full.Z)) colnames(geno.I) <- paste0("CovarMat", seq_len(ncol(geno.I))) full.Z <- do.call(cbind, list(full.Z, geno.I)) - # add a genetic variance component sigma_g <- Matrix(runif(1, 0, 1), ncol=1, nrow=1, sparse=TRUE) rownames(sigma_g) <- "CovarMat" diff --git a/R/testNhoods.R b/R/testNhoods.R index c943fcd..354f023 100644 --- a/R/testNhoods.R +++ b/R/testNhoods.R @@ -572,12 +572,25 @@ testNhoods <- function(x, design, design.df, kinship=NULL, "SE"= unlist(lapply(lapply(fit, `[[`, "SE"), function(x) x[ret.beta])), "tvalue" = unlist(lapply(lapply(fit, `[[`, "t"), function(x) x[ret.beta])), "PValue" = unlist(lapply(lapply(fit, `[[`, "PVALS"), function(x) x[ret.beta])), - matrix(unlist(lapply(fit, `[[`, "Sigma")), ncol=length(rand.levels), byrow=TRUE), "Converged"=unlist(lapply(fit, `[[`, "converged")), "Dispersion" = unlist(lapply(fit, `[[`, "Dispersion")), "Logliklihood"=unlist(lapply(fit, `[[`, "LOGLIHOOD"))) - rownames(res) <- 1:length(fit) - colnames(res)[6:(6+length(rand.levels)-1)] <- paste(names(rand.levels), "variance", sep="_") + # need to know how many variance components there are to get proper data frame dimensions + n.sigmas <- unique(unlist(lapply(lapply(fit, `[[`, "Sigma"), length))) + varcomps <- as.data.frame(matrix(unlist(lapply(fit, `[[`, "Sigma")), ncol=n.sigmas, byrow=TRUE)) + + if(n.sigmas == length(rand.levels)){ + colnames(varcomps) <- paste(names(rand.levels), "variance", sep="_") + } else{ + colnames(varcomps) <- paste(c(names(rand.levels), "CovarMat"), "variance", sep="_") + } + + res <- do.call(cbind.data.frame, list(res[, c("logFC", "logCPM", "SE", "tvalue", "PValue")], + varcomps, + res[, c("Converged", "Logliklihood")])) + + rownames(res) <- c(1:n.nhoods) + # colnames(res)[6:(6+length(rand.levels)-1)] <- paste(names(rand.levels), "variance", sep="_") } else { # need to use legacy=TRUE to maintain original edgeR behaviour fit <- glmQLFit(dge, x.model, robust=robust, legacy=TRUE) diff --git a/src/fitGeneticPLGlmm.cpp b/src/fitGeneticPLGlmm.cpp index 96c95a7..9dc67b8 100644 --- a/src/fitGeneticPLGlmm.cpp +++ b/src/fitGeneticPLGlmm.cpp @@ -220,6 +220,7 @@ List fitGeneticPLGlmm(const arma::mat& Z, const arma::mat& X, const arma::mat& K Vmu = computeVmu(muvec, curr_disp, vardist); W = computeW(curr_disp, Dinv, vardist); Winv = W.i(); + // pre-compute matrics: X^T * W^-1, Z^T * W^-1 arma::mat xTwinv = X.t() * Winv; arma::mat zTwin = Z.t() * Winv; @@ -272,6 +273,7 @@ List fitGeneticPLGlmm(const arma::mat& Z, const arma::mat& X, const arma::mat& K if(REML){ arma::mat VstarZ = V_star_inv * Z; VP_partial = precomp_list["PZZt"]; + VS_partial = pseudovarPartial_VG(u_indices, Z, VstarZ, K); score_sigma = sigmaScoreREML_arma(VP_partial, y_star, P, diff --git a/src/paramEst.cpp b/src/paramEst.cpp index 6b2cd07..064e41f 100644 --- a/src/paramEst.cpp +++ b/src/paramEst.cpp @@ -30,6 +30,7 @@ arma::vec sigmaScoreREML_arma (const Rcpp::List& pvstar_i, const arma::vec& ysta for(int i=0; i < c; i++){ const arma::mat& P_pvi = pvstar_i(i); // this is Vstar_inv * partial derivative const arma::mat& Pdifi = remldiffV(i); + double lhs = -0.5 * arma::trace(Pdifi); arma::mat mid1(1, 1); mid1 = arma::trans(ystarminx) * P_pvi * Vstarinv * ystarminx; diff --git a/src/pseudovarPartial.cpp b/src/pseudovarPartial.cpp index c177d34..40aa506 100644 --- a/src/pseudovarPartial.cpp +++ b/src/pseudovarPartial.cpp @@ -156,14 +156,15 @@ List pseudovarPartial_VG(const List& u_indices, const arma::mat& Z, const arma:: arma::uvec u_idx = u_indices[i]; arma::mat omat(n , n); + arma::mat VsZcols = VstarZ.cols(u_idx-1); + arma::mat Zcols = Z.cols(u_idx-1).t(); + if(i == c - 1){ - omat = VstarZ.cols(u_idx-1) * K * Z.cols(u_idx-1); + omat = VstarZ.cols(u_idx-1) * K * Zcols; } else{ - omat = VstarZ.cols(u_idx-1) * Z.cols(u_idx-1); + omat = VstarZ.cols(u_idx-1) * Zcols; } - - outlist[i] = omat; }