Merge pull request #372 from MarioniLab/devel

Resolve matrix commute error in model using random effect and GRM at the same time.

Author: MikeDMorgan

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"

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)

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,

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;

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;
     }