Merge branch 'devel'

Author: MikeDMorgan

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: miloR
 Type: Package
 Title: Differential neighbourhood abundance testing on a graph
-Version: 2.5.1
+Version: 2.5.2
 Authors@R:
     c(person("Mike", "Morgan", role=c("aut", "cre"), email="michael.morgan@abdn.ac.uk",
      comment=c(ORCID="0000-0003-0757-0711")),
@@ -15,9 +15,7 @@ URL: https://marionilab.github.io/miloR
 BugReports: https://github.com/MarioniLab/miloR/issues
 biocViews: SingleCell, MultipleComparison, FunctionalGenomics, Software
 LinkingTo: Rcpp,
-           RcppArmadillo,
-           RcppEigen,
-           RcppML
+           RcppArmadillo
 Depends: R (>= 4.0.0),
          edgeR
 Imports: BiocNeighbors,

NAMESPACE

@@ -97,6 +97,7 @@ importFrom(BiocParallel,bpstopOnError)
 importFrom(BiocParallel,bptry)
 importFrom(BiocSingular,bsparam)
 importFrom(Matrix,Matrix)
+importFrom(Matrix,chol)
 importFrom(Matrix,colSums)
 importFrom(Matrix,crossprod)
 importFrom(Matrix,diag)

R/glmm.R

@@ -76,7 +76,7 @@
 #' @name fitGLMM
 #'
 #' @importMethodsFrom Matrix %*%
-#' @importFrom Matrix Matrix solve crossprod kronecker
+#' @importFrom Matrix Matrix solve crossprod kronecker chol
 #' @importFrom stats runif var
 #' @importFrom BiocParallel bpstopOnError
 #' @export
@@ -131,7 +131,10 @@ fitGLMM <- function(X, Z, y, offsets, init.theta=NULL, Kin=NULL,
         }
 
         # create full Z with expanded random effect levels
+        # add a residual variance parameter
         full.Z <- initializeFullZ(Z=Z, cluster_levels=random.levels)
+        # random.levels <- c(random.levels, list("Residual"=rownames(Z)))
+
         if(is.null(glmm.control[["init.u"]])){
             curr_u <- matrix(runif(ncol(full.Z), 0, 1), ncol=1)
         } else{
@@ -172,13 +175,18 @@ fitGLMM <- function(X, Z, y, offsets, init.theta=NULL, Kin=NULL,
         } else{
             curr_sigma <- Matrix(glmm.control[["init.sigma"]], ncol=1, sparse=TRUE)
         }
-        rownames(curr_sigma) <- colnames(Z)
+
+        rownames(curr_sigma) <- names(random.levels)
 
         ## add the genetic components
         ## augment Z with I
-        geno.I <- diag(nrow(full.Z))
+        ## Use K=LD^1/2(LD^1/2)^T instead of ZZ^T for GRM
+        ## chol return L^T
+        geno.I <- t(chol(Kin))
+        # 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"
@@ -201,7 +209,10 @@ fitGLMM <- function(X, Z, y, offsets, init.theta=NULL, Kin=NULL,
 
         # if we only have a GRM then Z _is_ full.Z?
         # full.Z <- initializeFullZ(Z, cluster_levels=random.levels)
-        full.Z <- Z
+        ## Use K=LD^1/2(LD^1/2)^T instead of ZZ^T for GRM
+        ## chol return L^T
+        full.Z <- t(chol(Kin))
+        # full.Z <- Z
         # should this be the matrix R?
         colnames(full.Z) <- paste0(names(random.levels), seq_len(ncol(full.Z)))
 
@@ -255,6 +266,7 @@ fitGLMM <- function(X, Z, y, offsets, init.theta=NULL, Kin=NULL,
     } else if(is.null(Kin)){
         # create full Z with expanded random effect levels
         full.Z <- initializeFullZ(Z=Z, cluster_levels=random.levels)
+        # random.levels <- c(random.levels, list("Residual"=rownames(Z)))
 
         # random value initiation from runif
         if(is.null(glmm.control[["init.u"]])){
@@ -456,7 +468,7 @@ initialiseG <- function(cluster_levels, sigmas, Kin=NULL){
             diag(G[c(i:(i+x.q-1)), c(i:(i+x.q-1)), drop=FALSE]) <- sigmas[x, ] # is this sufficient to transform the sigma to the model scale?
         } else{
             if(rownames(sigmas[x, , drop=FALSE]) %in% c("Genetic")){
-                diag(G[c(i:(i+x.q-1)), c(i:(i+x.q-1)), drop=FALSE]) <- sigmas[x, ] * Kin
+                diag(G[c(i:(i+x.q-1)), c(i:(i+x.q-1)), drop=FALSE]) <- sigmas[x, ] #* Kin
             }else{
                 diag(G[c(i:(i+x.q-1)), c(i:(i+x.q-1)), drop=FALSE]) <- sigmas[x, ] # is this sufficient to transform the sigma to the model scale?
             }
@@ -565,7 +577,12 @@ initializeFullZ <- function(Z, cluster_levels, stand.cols=FALSE){
 
         i.z.list[[colnames(Z)[i]]] <- i.z
     }
+    # # Add the residual variance
+    # resid.eye <- diag(nrow(Z))
+    # colnames(resid.eye) <- rownames(Z)
+    # i.z.list[["Residual"]] <- resid.eye
     full.Z <- do.call(cbind, i.z.list)
+
     return(full.Z)
 }
 

R/testNhoods.R

@@ -539,7 +539,7 @@ testNhoods <- function(x, design, design.df, kinship=NULL,
             if(isTRUE(geno.only)){
                 message("Running genetic model with ", nrow(kinship), " individuals")
             } else{
-                message("Running genetic model with ", nrow(z.model), " observations")
+                message("Running genetic+ model with ", nrow(z.model), " observations")
             }
         }
 

src/paramEst.cpp

@@ -2,7 +2,6 @@
 #include "paramEst.h"
 #include "computeMatrices.h"
 #include "utils.h"
-#include "solveQP.h"
 #include<RcppArmadillo.h>
 #ifdef _OPENMP
 #include <omp.h>
@@ -358,8 +357,7 @@ arma::vec estHasemanElstonConstrained(const arma::mat& Z, const arma::mat& PREML
     // solve by linear least squares
     arma::vec _he_update(c+1);
 
-    // use RcppML NNLS (needs casting Eigen <-> arma)
-    // _he_update = solveQP(vecZ, Ybig, _he_update);
+    // Tried RcppML NNLS but arma <> eigen cast is too expensive
     _he_update = nnlsSolve(vecZ, Ybig, _he_update, Iters);
 
     return _he_update;
@@ -390,7 +388,6 @@ arma::vec estHasemanElstonConstrainedML(const arma::mat& Z, const Rcpp::List& u_
 
     // solve by linear least squares
     arma::vec _he_update(c+1);
-    // _he_update = solveQP(vecZ, Ybig, _he_update);
     _he_update = nnlsSolve(vecZ, Ybig, _he_update, Iters);
 
     return _he_update;
@@ -424,7 +421,6 @@ arma::vec estHasemanElstonConstrainedGenetic(const arma::mat& Z, const arma::mat
     arma::mat vecZ = vectoriseZGenetic(Z, u_indices, PREML, PZ, Kin); // projection already applied
 
     arma::vec _he_update(c+1);
-    // _he_update = solveQP(vecZ, Ybig, _he_update);
     _he_update = nnlsSolve(vecZ, Ybig, _he_update, Iters);
 
     return _he_update;
@@ -455,11 +451,10 @@ arma::vec estHasemanElstonConstrainedGeneticML(const arma::mat& Z,
     arma::mat vecZ = vectoriseZGeneticML(Z, u_indices, Kin); // projection already applied
 
     arma::vec _he_update(c+1);
-    // _he_update = solveQP(vecZ, Ybig, _he_update);
     _he_update = nnlsSolve(vecZ, Ybig, _he_update, Iters);
 
-    return _he_update;
-}
+    return _he_update;}
+
 
 
 arma::vec nnlsSolve(const arma::mat& vecZ, const arma::vec& Y, arma::vec nnls_update, const int& Iters){
@@ -700,7 +695,12 @@ arma::mat vectoriseZGenetic(const arma::mat& Z, const Rcpp::List& u_indices,
 
         // always set the last component to the genetic variance if there is a kinship matrix
         if(i == c-1){
-            arma::vec _vecZ = Kin(lower_indices);
+            // arma::vec _vecZ = Kin(lower_indices);
+            // Cholesky LL^T of kinship == ZZ^T
+            arma::mat _ZZT = PZ.cols(u_idx - 1) * Z.cols(u_idx - 1).t() * P.t(); // REML projection
+
+            // vectorise
+            arma::vec _vecZ = _ZZT(lower_indices);
             vecMat.col(i+1) = _vecZ;
             // vecMat.col(i) = _vecZ;
         } else{
@@ -741,7 +741,12 @@ arma::mat vectoriseZGeneticML(const arma::mat& Z, const Rcpp::List& u_indices,
 
         // always set the last component to the genetic variance if there is a kinship matrix
         if(i == c-1){
-            arma::vec _vecZ = Kin(lower_indices);
+            // arma::vec _vecZ = Kin(lower_indices);
+            // Cholesky LL^T of kinship == ZZ^T
+            arma::mat _ZZT = Z.cols(u_idx - 1) * Z.cols(u_idx - 1).t();
+
+            // vectorise
+            arma::vec _vecZ = _ZZT(lower_indices);
             vecMat.col(i+1) = _vecZ;
             // vecMat.col(i) = _vecZ;
         } else{

src/pseudovarPartial.cpp

@@ -75,7 +75,8 @@ List computePZList_G(const List& u_indices, const arma::mat& PZ, const arma::mat
         arma::mat _pzz(n, n);
 
         if(i == c - 1){
-            _pzz = PZ.cols(u_idx-1) * K * Z.cols(u_idx-1).t(); // convert 1-based to 0-based
+            // _pzz = PZ.cols(u_idx-1) * K * Z.cols(u_idx-1).t(); // convert 1-based to 0-based
+            _pzz = PZ.cols(u_idx-1) * Z.cols(u_idx-1).t(); // convert 1-based to 0-based
         } else{
             _pzz = PZ.cols(u_idx-1) * Z.cols(u_idx-1).t(); // convert 1-based to 0-based
         }
@@ -160,7 +161,8 @@ List pseudovarPartial_VG(const List& u_indices, const arma::mat& Z, const arma::
         arma::mat Zcols = Z.cols(u_idx-1).t();
 
         if(i == c - 1){
-            omat = VstarZ.cols(u_idx-1) * K * Zcols;
+            //omat = VstarZ.cols(u_idx-1) * K * Zcols;
+            omat = VstarZ.cols(u_idx-1) * Zcols;
         } else{
             omat = VstarZ.cols(u_idx-1) * Zcols;
         }
@@ -181,7 +183,8 @@ List pseudovarPartial_G(arma::mat Z, const arma::mat& K, List u_indices){
     for(unsigned int i = 0; i < items; i++){
         if(i == items - 1){
             arma::uvec icols = u_indices[i];
-            arma::mat _omat(Z.cols(icols - 1) * K * Z.cols(icols - 1).t());
+            // arma::mat _omat(Z.cols(icols - 1) * K * Z.cols(icols - 1).t());
+            arma::mat _omat(K);
             outlist[i] = _omat; // K is equivalent to ZZ^T
         } else{
             arma::uvec icols = u_indices[i];