diff --git a/src/paramEst.cpp b/src/paramEst.cpp
index 786a2ec..6b2cd07 100644
--- a/src/paramEst.cpp
+++ b/src/paramEst.cpp
@@ -358,7 +358,8 @@ arma::vec estHasemanElstonConstrained(const arma::mat& Z, const arma::mat& PREML
     arma::vec _he_update(c+1);
 
     // use RcppML NNLS (needs casting Eigen <-> arma)
-    _he_update = solveQP(vecZ, Ybig, _he_update);
+    // _he_update = solveQP(vecZ, Ybig, _he_update);
+    _he_update = nnlsSolve(vecZ, Ybig, _he_update, Iters);
 
     return _he_update;
 }
@@ -388,7 +389,8 @@ 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 = solveQP(vecZ, Ybig, _he_update);
+    _he_update = nnlsSolve(vecZ, Ybig, _he_update, Iters);
 
     return _he_update;
 }
@@ -421,7 +423,8 @@ 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 = solveQP(vecZ, Ybig, _he_update);
+    _he_update = nnlsSolve(vecZ, Ybig, _he_update, Iters);
 
     return _he_update;
 }
@@ -451,7 +454,8 @@ 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 = solveQP(vecZ, Ybig, _he_update);
+    _he_update = nnlsSolve(vecZ, Ybig, _he_update, Iters);
 
     return _he_update;
 }
diff --git a/vignettes/milo_glmm.Rmd b/vignettes/milo_glmm.Rmd
index 9fabd6a..c6d31d7 100644
--- a/vignettes/milo_glmm.Rmd
+++ b/vignettes/milo_glmm.Rmd
@@ -167,7 +167,7 @@ set.seed(42)
 rownames(pbmc.design) <- pbmc.design$ObsID
 
 da_results <- testNhoods(pbmc.milo, design = ~ adjmfc.time + (1|tenx_lane), design.df = pbmc.design, fdr.weighting="graph-overlap",
-                         glmm.solver="Fisher", REML=TRUE, norm.method="TMM", BPPARAM = bpparam, fail.on.error=FALSE, intercept.type="fixed")
+                         glmm.solver="Fisher", REML=TRUE, norm.method="TMM", BPPARAM = bpparam, fail.on.error=FALSE)
 table(da_results$SpatialFDR < 0.1)
 ```