use COD to create sufficient statistic

Author: bgoodri

inst/include/CODOLS.hpp

@@ -0,0 +1,23 @@
+#ifndef RSTANARM__CODOLS_HPP
+#define RSTANARM__CODOLS_HPP
+
+/*
+ * Compute ordinary least squares coefficients,
+ * even in the situation where X is rank deficient
+ * See https://eigen.tuxfamily.org/dox/classEigen_1_1CompleteOrthogonalDecomposition.html
+ */
+
+template <typename T2__, typename T3__>
+inline
+Eigen::Matrix<typename boost::math::tools::promote_args<T2__, T3__>::type,
+              Eigen::Dynamic, 1>
+CODOLS(const Eigen::Matrix<T2__, Eigen::Dynamic, Eigen::Dynamic>& X,
+       const Eigen::Matrix<T3__, Eigen::Dynamic, 1>& y,
+       std::ostream* pstream__) {
+  typename boost::math::tools::promote_args<T2__, T3__>::type T1__;
+  using namespace Eigen;
+  CompleteOrthogonalDecomposition<MatrixXd> cod(X);
+  return cod.solve(y);
+}
+#endif
+

inst/include/meta_header.hpp

@@ -2,4 +2,5 @@
 #define RSTANARM__META_HEADER_HPP
 
 #include "csr_matrix_times_vector2.hpp"
+#include "CODOLS.hpp"
 #endif

src/stan_files/continuous.stan

@@ -24,6 +24,8 @@ functions {
             - N * (log(sigma) + log(sqrt(2 * pi())));
   }
 
+  vector CODOLS(matrix X, vector y); // implemented in C++
+
   /** 
   * test function for csr_matrix_times_vector
   *
@@ -94,13 +96,9 @@ transformed data {
     matrix[N, has_intercept + K + K_smooth ] X_ = has_intercept ? append_col(rep_vector(1.0, N), 
                                                   (K_smooth > 0 ? append_col(X[1], S) : X[1])) : 
                                                   (K_smooth > 0 ? append_col(X[1], S) : X[1]);
-    matrix[cols(X_), cols(X_)] R = qr_thin_R(X_);
-    if (tail(diagonal(R), 1)[1] > 1e-16) {
-      matrix[N, cols(R)] Q = qr_thin_Q(X_);
-      XtX = crossprod(X_);
-      OLS = mdivide_right_tri_low(y' * Q, R')';
-      SSR = dot_self(y - X_ * OLS);
-    } else can_do_OLS = 0;
+    XtX = crossprod(X_);
+    OLS = CODOLS(X_, y);
+    SSR = dot_self(y - X_ * OLS);
   }
   if (can_do_normalidglm) {
     XS = K_smooth > 0 ? append_col(X[1], S) : X[1];

tests/testthat/test_stan_glm.R

@@ -466,3 +466,20 @@ test_that("contrasts attribute isn't dropped", {
                  chains = 1, refresh = 0)
   expect_equal(fit$contrasts, contrasts)
 })
+
+test_that("returns something with collinear predictors", {
+  N <- 100
+  y <- rnorm(N)
+  z <- sample(c(0,1), N, replace=TRUE)
+  x1 <- rnorm(N)
+  x2 <- 2*x1
+
+  fit_1 <- stan_glm(
+    y ~ z * (x1 + x2),
+    data = data.frame(y, z, x1, x2),
+    prior = normal(location = 0, scale = 0.1),
+    prior_intercept = normal(location = 0, scale = 0.1),
+    chains = CHAINS, iter = ITER, refresh = REFRESH
+  )
+  expect_stanreg(fit_1)  
+})