From 5f776a6997ab8f8d84b107b835a675b6f50a0e25 Mon Sep 17 00:00:00 2001
From: Ben Dilday <ben.dilday.phd@gmail.com>
Date: Fri, 9 Mar 2018 19:08:51 -0500
Subject: [PATCH] model results as data frame

---
 NAMESPACE       |  4 ++++
 R/RcppExports.R |  4 ++++
 R/main.R        |  3 +++
 R/mnre_fit.R    | 54 ++++++++++++++++++++++++++++++++++++++++++++++++-
 src/mnre.cpp    |  4 ++++
 5 files changed, 68 insertions(+), 1 deletion(-)

diff --git a/NAMESPACE b/NAMESPACE
index 6773b62..496f149 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -4,9 +4,13 @@ export(mnre_dim_and_class_to_index)
 export(mnre_expand_matrix)
 export(mnre_fit)
 export(mnre_fit_sparse)
+export(mnre_left_covar_factor)
+export(mnre_lk_glm)
+export(mnre_lk_penalty)
 export(mnre_make_covar)
 export(mnre_simulate_ev_data)
 export(mnre_simulate_multinomial_data_factors)
+export(mnre_step_sparse)
 export(nd_min_fun)
 import(Matrix)
 importFrom(Rcpp,sourceCpp)
diff --git a/R/RcppExports.R b/R/RcppExports.R
index 77fcc53..b9d0ca3 100644
--- a/R/RcppExports.R
+++ b/R/RcppExports.R
@@ -44,6 +44,7 @@ mnre_expand_matrix <- function(x1, k_class, direction) {
     .Call('_mnre_mnre_expand_matrix', PACKAGE = 'mnre', x1, k_class, direction)
 }
 
+#' @export
 mnre_left_covar_factor <- function(x1) {
     .Call('_mnre_mnre_left_covar_factor', PACKAGE = 'mnre', x1)
 }
@@ -60,10 +61,12 @@ mnre_fit_sparse <- function(fixed_effects, random_effects, y, theta_mat, Lind, b
     .Call('_mnre_mnre_fit_sparse', PACKAGE = 'mnre', fixed_effects, random_effects, y, theta_mat, Lind, beta_fixed, beta_random, verbose)
 }
 
+#' @export
 mnre_lk_penalty <- function(beta_random, theta_norm, Lind) {
     .Call('_mnre_mnre_lk_penalty', PACKAGE = 'mnre', beta_random, theta_norm, Lind)
 }
 
+#' @export
 mnre_lk_glm <- function(fixed_effects, random_effects, beta_fixed, beta_random, y, Lind) {
     .Call('_mnre_mnre_lk_glm', PACKAGE = 'mnre', fixed_effects, random_effects, beta_fixed, beta_random, y, Lind)
 }
@@ -80,6 +83,7 @@ mnre_mu_x <- function(fe_x, re_x, beta_fixed, beta_random) {
     .Call('_mnre_mnre_mu_x', PACKAGE = 'mnre', fe_x, re_x, beta_fixed, beta_random)
 }
 
+#' @export
 mnre_step_sparse <- function(fixed_effects, random_effects, y, beta_fixed, beta_random, lambda_norm, Lind) {
     .Call('_mnre_mnre_step_sparse', PACKAGE = 'mnre', fixed_effects, random_effects, y, beta_fixed, beta_random, lambda_norm, Lind)
 }
diff --git a/R/main.R b/R/main.R
index 4c22a8c..d074e63 100644
--- a/R/main.R
+++ b/R/main.R
@@ -203,6 +203,9 @@ nd_min_fun <- function(ev) {
     ev$verbose <- 1
   }
   
+  # make sure the data is a data frame, not a tibble
+  ev$fr <- as.data.frame(ev$fr)
+  
   function(mval) {
     glf <- lme4::glFormula(ev$frm,
                            data=ev$fr, family='binomial')
diff --git a/R/mnre_fit.R b/R/mnre_fit.R
index 875f88d..032eaee 100644
--- a/R/mnre_fit.R
+++ b/R/mnre_fit.R
@@ -13,4 +13,56 @@ mnre_fit <- function(frm, data, verbose=0, off_diagonal=0.0) {
   
   ans = optim(mval, nf, method = "L-BFGS", lower=1e-8)
 
-}
\ No newline at end of file
+  mnre_fit_to_df(frm, data, ans$par, verbose = verbose, off_diagonal = off_diagonal)
+  
+}
+
+#' 
+mnre_fit_to_df <- function(frm, data, mval, verbose=0, off_diagonal=0.0) {
+  data <- as.data.frame(data)
+  glf <- lme4::glFormula(frm,
+                         data, family='binomial')
+  fe <- fixed_effects <- (glf$X)
+  re <- random_effects <- Matrix::t(glf$reTrms$Zt)
+  
+  y <- matrix(data[,all.vars(frm)[[1]]], ncol=1)
+  k_class <- max(y)
+  k <- max(y)
+  Lind = matrix(glf$reTrms$Lind, ncol=1)
+  
+  theta_mat <- matrix(mval, ncol=k_class)
+  covar_mat = mnre_make_covar(theta_mat, Lind, off_diagonal = 0.0)
+  left_factor <- mnre_left_covar_factor(covar_mat)
+  
+  fe_sp <- Matrix::Matrix(fe, sparse = TRUE)
+
+  beta_re <- matrix(rnorm(ncol(re) * k_class), ncol=k_class)
+  beta_fe <- matrix(rnorm(ncol(fe) * k_class), ncol=k_class)
+  
+  zz <- mnre_fit_sparse(fe_sp, re, y, theta_mat, Lind, beta_fe, beta_re, verbose = verbose)
+  
+  lk1 <- zz$loglk + zz$loglk_det
+  bpar = matrix(left_factor %*% matrix(zz$beta_random,ncol=1), ncol=k_class)
+  
+  lvs <- unlist(sapply(glf$reTrms$flist, levels))
+  cc1 <- as.data.frame(cbind(bpar, Lind=Lind))
+
+  ranef_labels <- names(glf$reTrms$cnms)
+  df_names <- sapply(1:k_class, function(i) {sprintf("class%02d", i)})
+  df_names <- c(df_names, "Lind")
+  df_names <- c(df_names, "ranef_label")
+  df_names <- c(df_names, "ranef_level")
+  
+  cc1$ranef <- ranef_labels[cc1[,ncol(cc1)]]
+  cc1$lv <- matrix(lvs, ncol=1)
+  
+  mvalX <- t(sapply(1:max(Lind), function(i) {
+    idx = which(cc1[,k_class+1] == i)
+    tmp <- matrix(bpar[idx,], ncol=k_class)
+    apply(tmp, 2, sd)
+  }))
+  
+  names(cc1) <- df_names
+  list(ranef=cc1, fixef=zz$beta_fixed, theta=mval)
+  
+}
diff --git a/src/mnre.cpp b/src/mnre.cpp
index 787276e..fe00b7d 100644
--- a/src/mnre.cpp
+++ b/src/mnre.cpp
@@ -114,6 +114,7 @@ arma::sp_mat mnre_expand_matrix(const arma::sp_mat& x1, int k_class, int directi
   return expanded_mat;
 }
 
+//' @export
 // [[Rcpp::export]]
 arma::sp_mat mnre_left_covar_factor(arma::sp_mat& x1) {
   arma::sp_mat left_factor;
@@ -304,6 +305,7 @@ Rcpp::List mnre_fit_sparse(const arma::sp_mat& fixed_effects,
 } // end function
 
 
+//' @export
 // [[Rcpp::export]]
 double mnre_lk_penalty(const arma::mat& beta_random,
                            const arma::mat& theta_norm,
@@ -328,6 +330,7 @@ double mnre_lk_penalty(const arma::mat& beta_random,
 }
 
 
+//' @export
 // [[Rcpp::export]]
 double mnre_lk_glm(const arma::sp_mat& fixed_effects,
                        const arma::sp_mat& random_effects,
@@ -448,6 +451,7 @@ arma::mat mnre_mu_x(const arma::sp_mat &fe_x,
   return mu;
 }
 
+//' @export
 // [[Rcpp::export]]
 Rcpp::List mnre_step_sparse(const arma::sp_mat &fixed_effects,
                                 const arma::sp_mat &random_effects,