diff --git a/.Rbuildignore b/.Rbuildignore new file mode 100644 index 0000000..91114bf --- /dev/null +++ b/.Rbuildignore @@ -0,0 +1,2 @@ +^.*\.Rproj$ +^\.Rproj\.user$ diff --git a/.gitignore b/.gitignore index fcff087..886732b 100644 --- a/.gitignore +++ b/.gitignore @@ -31,3 +31,4 @@ vignettes/*.pdf # Temporary files created by R markdown *.utf8.md *.knit.md +.Rproj.user diff --git a/DESCRIPTION b/DESCRIPTION new file mode 100644 index 0000000..427b345 --- /dev/null +++ b/DESCRIPTION @@ -0,0 +1,14 @@ +Package: mnre +Title: Multinomial Random Effects Models +Version: 0.1 +Authors@R: person("Ben", "Dilday", email = "ben.dilday.phd@gmail.com", role = c("aut", "cre")) +Description: Multinomial Random Effects Models +Depends: R (>= 3.4.1) +License: MIT +Encoding: UTF-8 +LazyData: true +Imports: Rcpp (>= 0.12.12) +LinkingTo: Rcpp, + RcppArmadillo, + RcppEigen +RoxygenNote: 6.0.1 diff --git a/NAMESPACE b/NAMESPACE new file mode 100644 index 0000000..4bc04ca --- /dev/null +++ b/NAMESPACE @@ -0,0 +1,13 @@ +# Generated by roxygen2: do not edit by hand + +export(mnre_dim_and_class_to_index) +export(mnre_expand_matrix) +export(mnre_fit_sparse) +export(mnre_make_covar) +export(mnre_simulate_ev_data) +export(mnre_simulate_multinomial_data_factors) +export(nd_min_fun) +import(Matrix) +importFrom(Rcpp,sourceCpp) +importFrom(magrittr,"%>%") +useDynLib(mnre) diff --git a/R/RcppExports.R b/R/RcppExports.R new file mode 100644 index 0000000..77fcc53 --- /dev/null +++ b/R/RcppExports.R @@ -0,0 +1,90 @@ +# Generated by using Rcpp::compileAttributes() -> do not edit by hand +# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 + +#' Index from dimension and class indicies +#' +#' @param i_dim integer +#' @param i_class integer +#' @param n_dim integer +#' @examples +#' mnre_dim_and_class_to_index(1, 1, 1) +#' mnre_dim_and_class_to_index(2, 3, 5) +#' @return The integer index corresponding to the dimension and class indicies +#' @export +mnre_dim_and_class_to_index <- function(i_dim, i_class, n_dim) { + .Call('_mnre_mnre_dim_and_class_to_index', PACKAGE = 'mnre', i_dim, i_class, n_dim) +} + +#' Make covariance matrix +#' +#' @param theta_mat numeric matrix +#' @param Lind unsigned integer matrix +#' @param off_diagonal numeric +#' @examples +#' \dontrun{ +#' theta_mat <- matrix(1:4, ncol=2) +#' Lind <- matrix(c(rep(1, 10), rep(2, 5)), ncol=1) +#' covar_mat <- mnre_make_covar(theta_mat, Lind) +#' } +#' @return Covariance matrix for random effects +#' @export +mnre_make_covar <- function(theta_mat, Lind, off_diagonal = 0.0) { + .Call('_mnre_mnre_make_covar', PACKAGE = 'mnre', theta_mat, Lind, off_diagonal) +} + +#' Expand matrix +#' +#' Used to duplicate the random effects coefficient matrix before applying the left factor +#' @param x1 sparse matrix +#' @param k_class integer +#' @param direction integer +#' @return Matrix expanded in the The integer index corresponding to the dimension and class indicies +#' @export +mnre_expand_matrix <- function(x1, k_class, direction) { + .Call('_mnre_mnre_expand_matrix', PACKAGE = 'mnre', x1, k_class, direction) +} + +mnre_left_covar_factor <- function(x1) { + .Call('_mnre_mnre_left_covar_factor', PACKAGE = 'mnre', x1) +} + +#' Multinomial mixed effects fit +#' +#' Used to duplicate the random effects coefficient matrix before applying the left factor +#' @param x1 sparse matrix +#' @param k_class integer +#' @param direction integer +#' @return Matrix expanded in the The integer index corresponding to the dimension and class indicies +#' @export +mnre_fit_sparse <- function(fixed_effects, random_effects, y, theta_mat, Lind, beta_fixed, beta_random, verbose = 1L) { + .Call('_mnre_mnre_fit_sparse', PACKAGE = 'mnre', fixed_effects, random_effects, y, theta_mat, Lind, beta_fixed, beta_random, verbose) +} + +mnre_lk_penalty <- function(beta_random, theta_norm, Lind) { + .Call('_mnre_mnre_lk_penalty', PACKAGE = 'mnre', beta_random, theta_norm, Lind) +} + +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) +} + +mnre_lk <- function(fixed_effects, random_effects, beta_fixed, beta_random, y, theta_norm, Lind) { + .Call('_mnre_mnre_lk', PACKAGE = 'mnre', fixed_effects, random_effects, beta_fixed, beta_random, y, theta_norm, Lind) +} + +mnre_mu <- function(fixed_effects, random_effects, beta_fixed, beta_random) { + .Call('_mnre_mnre_mu', PACKAGE = 'mnre', fixed_effects, random_effects, beta_fixed, beta_random) +} + +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) +} + +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) +} + +fill_mtwm_x <- function(x1, x2, mu) { + .Call('_mnre_fill_mtwm_x', PACKAGE = 'mnre', x1, x2, mu) +} + diff --git a/R/main.R b/R/main.R new file mode 100644 index 0000000..3fdc4a6 --- /dev/null +++ b/R/main.R @@ -0,0 +1,243 @@ + +#' @import Matrix +#' @importFrom magrittr %>% +#' @export +mnre_simulate_ev_data <- function(nlim=1000, year=2016, OBP=FALSE) { + ev_obj <- BProDRA:::generate_model_df(nlim=nlim, year=year) + dfX <- ev_obj$ev + + # the BPRO code codes it like this + # ev[cc_bip0,]$outcome <- 1 + # ev[cc_bip1,]$outcome <- 2 + # ev[cc_bip2,]$outcome <- 3 + # ev[cc_hr,]$outcome <- 4 + # ev[cc_so,]$outcome <- 5 + # ev[cc_bb,]$outcome <- 6 + + idx1 <- which(dfX$outcome == 1) + idx2 <- which(dfX$outcome == 2) + idx3 <- which(dfX$outcome == 3) + idx4 <- which(dfX$outcome == 4) + idx5 <- which(dfX$outcome == 5) + idx6 <- which(dfX$outcome == 6) + + # make batted ball for out the reference + dfX[idx1,]$outcome <- 0 + + # only 1 kind of batted ball for hit + dfX[c(idx2, idx3),]$outcome <- 1 + + # shift the other categories up + # HR + dfX[idx4,]$outcome <- 2 + + # shift the other categories up + # SO + dfX[idx5,]$outcome <- 3 + + # shift the other categories up + # BB + dfX[idx6,]$outcome <- 4 + + # OBP + if (OBP) { + dfX[c(idx1, idx5),]$outcome <- 0 + dfX[c(idx2, idx3, idx4, idx6),]$outcome <- 1 + } + + dfX$y <- dfX$outcome + + glf <- lme4::glFormula(y ~ (1|BAT_ID) + (1|PIT_ID) + (1|HOME_TEAM_ID), + data=dfX, family='binomial') + #glf <- lme4::glFormula(y ~ (1|PIT_ID) + (1|BAT_ID), + # data=dfX, family='binomial') + + dfX <- dfX %>% dplyr::mutate(pid=as.integer(as.factor(PIT_ID))) + dfX <- dfX %>% dplyr::mutate(bid=as.integer(as.factor(BAT_ID))) + + dfX <- dfX %>% dplyr::mutate(sid=as.integer(as.factor(HOME_TEAM_ID))) + dfX <- dfX %>% dplyr::mutate(bid = bid+max(pid)) + dfX <- dfX %>% dplyr::mutate(sid = sid+max(bid)) + + re <- random_effects <- Matrix::t(glf$reTrms$Zt) + fe <- fixed_effects <- (glf$X) + k_class <- max(dfX$y) + beta_re <- matrix(rep(0, k_class * dim(re)[[2]]), ncol=k_class) + beta_fe <- matrix(rep(0, k_class * dim(fe)[[2]]), ncol=k_class) + xx <- model.matrix(y ~ BAT_ID + PIT_ID + HOME_TEAM_ID, data=dfX) +# xx <- model.matrix(y ~ BAT_ID + PIT_ID, data=dfX) + y <- matrix(dfX$y, ncol=1) + theta_init <- matrix(rep(glf$reTrms$theta, k_class), ncol=k_class) + list(fr=dfX, frm=glf$formula, + ev_obj=ev_obj, xx=xx, y=y, + re=re, fe=fe, + beta_re=beta_re, beta_fe=beta_fe, + theta_init=theta_init, + Lind=matrix(glf$reTrms$Lind, ncol=1)) +} + + + +#' @export +mnre_simulate_multinomial_data_factors <- function(rseed=101, + nfct=2, + nlev=10, + K_class=3, + nobs=10000) { + + set.seed(rseed) + + k <- K_class - 1 + fcts <- lapply(1:nfct, function(i) {sprintf("%s%03d", LETTERS[i], 1:nlev)}) + sigmas <- matrix(rep(1, nfct * k), ncol=k) + + coef_null_class <- matrix(rnorm(nfct * length(fcts[[1]])), ncol=nfct) + + # the null class needs to be identically 0, otherwise there will be a + # spurious corelation between coefficients for the other classes. + coef_null_class = coef_null_class * 0 + + coefs <- lapply(1:nfct, function(i) { + tmp <- as.data.frame( + matrix(rnorm(k * length(fcts[[i]])), ncol=k) - coef_null_class[,i] + ) + tmp$fct <- as.factor(fcts[[i]]) + tmp + }) + + df_obs <- lapply(1:length(coefs), function(i) { + dplyr::data_frame(fct=sample(fcts[[i]], nobs, replace = TRUE)) + }) + + df_obs <- dplyr::bind_cols(df_obs) + + names(df_obs) <- sprintf("fct%02d", 1:nfct) + + dfX <- df_obs + for (i in 1:length(coefs)) { + dfX <- dfX %>% merge(coefs[[i]], by.x=sprintf("fct%02d", i), by.y="fct") + } + + # must be identically 0 + icpt_null_class <- 0 + icpts <- matrix(rnorm(k), ncol=k) - icpt_null_class + + oc <- sapply(1:nrow(dfX), function(i) { + + r <- dfX[i,] + seq(nfct+1,ncol(r),k) + + lams <- sapply(1:k, function(j) { + noise <- rnorm(1, 0, 0.2) + ii=seq(nfct+j,ncol(r),k) + icpts[[j]] + sum(r[ii]) + noise + }) + + probs <- c(1, exp(lams)) + + tmp <- rmultinom(1, 1, probs) + as.integer(which(tmp > 0) - 1 ) + }) + + dfX <- dfX %>% dplyr::mutate(y=oc) + fct_cut <- which(grepl('^fct|^y$', names(dfX))) + + if (length(fct_cut) > 1) { + dfY <- dfX[,fct_cut] + } else { + dfY <- dplyr::data_frame(fct01=dfX[,fct_cut]) + } + + s = ' y ~ (1|fct01) ' + + if (nfct >= 2) { + for (i in 2:nfct) { + s <- sprintf("%s + (1|fct%02d) ", s, i) + } + } + + frm <- as.formula(s) + glf <- lme4::glFormula(frm, + data=dfY, family='binomial') + re <- random_effects <- Matrix::t(glf$reTrms$Zt) + fe <- fixed_effects <- (glf$X) + k_class <- max(dfY$y) + beta_re <- matrix(rep(0, k_class * dim(re)[[2]]), ncol=k_class) + beta_fe <- matrix(rep(0, k_class * dim(fe)[[2]]), ncol=k_class) + y <- matrix(dfY$y, ncol=1) + theta_init <- matrix(rep(glf$reTrms$theta, k_class), ncol=k_class) + ans <- list(true_pars=dfX, + fr=glf$fr, y=y, + re=re, fe=fe, + frm=frm, + beta_re=beta_re, beta_fe=beta_fe, + theta_init=theta_init, + Lind=matrix(glf$reTrms$Lind, ncol=1), + off_diagonal=0.0) + + ans + +} + + +#' N-dimensional function generator +#' @param ev list +#' @return deviance function +#' @examples +#' \dontrun{ +#' ev = mnre_simulate_multinomial_data_factors(nfct=2, K_class = 2, nlev=50, nobs=20000) +#' nf <- nd_min_fun(ev) +#' nf(c(1,1)) +#' } +#' @export +nd_min_fun <- function(ev) { + + frm <- ev$frm + if ("off_diagonal" %in% names(ev)) { + ev$off_diagonal <- ev$off_diagonal + } else { + ev$off_diagonal <- 0.0 + } + + if (!"verbose" %in% names(ev)) { + ev$verbose <- 1 + } + + function(mval) { + glf <- lme4::glFormula(ev$frm, + data=ev$fr, family='binomial') + fe <- fixed_effects <- (glf$X) + re <- random_effects <- Matrix::t(glf$reTrms$Zt) + y <- matrix(ev$y, ncol=1) + k_class <- max(y) + k <- max(y) + Lind = ev$Lind + + s = 'mval ' + for (v in mval) { + s = sprintf("%s %.4e ", s, v) + } + message(s) + + theta_mat <- matrix(mval, ncol=k_class) + + fe2 <- Matrix::Matrix(fe, sparse = TRUE) + + beta_re <- ev$beta_re + beta_fe <- ev$beta_fe + + # beta_re <- matrix(rep(0, ncol(re) * k_class), ncol=k_class) + # beta_fe <- matrix(rep(0, ncol(fe) * k_class), ncol=k_class) + # + message("starting beta ", beta_fe[[1]], " ", beta_re[[1]]) + zz <- mnre_fit_sparse(fe2, re, y, theta_mat, Lind, beta_fe, beta_re, verbose = ev$verbose) + + ev$beta_fe <<- zz$beta_fixed + ev$beta_re <<- zz$beta_random + + zz$loglk + zz$loglk_det + + } +} + + diff --git a/R/mnre.R b/R/mnre.R new file mode 100644 index 0000000..adea927 --- /dev/null +++ b/R/mnre.R @@ -0,0 +1,4 @@ +#' @useDynLib mnre +#' @importFrom Rcpp sourceCpp +NULL + diff --git a/README.md b/README.md index fefabe1..91ba812 100644 --- a/README.md +++ b/README.md @@ -1,2 +1,42 @@ # mnre -multinomial random effects + +Multinomial Random Effects + +## Usage + +This code solves multinomial mixed effects models. + +A minimal example is: + +fit a binomial model + +``` +devtools::install_github('bdilday/mnre') +library(mnre) +ev = mnre_simulate_multinomial_data_factors(nfct=2, K_class = 2, nlev=50, nobs=20000) +nf = mnre::nd_min_fun(ev) +ans = optim(c(1,1), nf, method = "L-BFGS", lower=1e-8) +print(ans$par) +[1] 0.9926043 0.9537682 +``` + +Compare to `lme4` + +``` +glmer_mod <- glmer(ev$frm, data=ev$fr, family='binomial', nAGQ=0) +glmer_mod@theta +[1] 0.9925264 0.9537615 +``` + +Fit a multinomial model + +``` +ev = mnre_simulate_multinomial_data_factors(nfct=1, K_class = 3, nlev=50, nobs=20000) +ev$verbose = 0 +nf = mnre::nd_min_fun(ev) +ans = optim(c(1, 1), nf, method = "L-BFGS", lower=1e-8) +print(ans$par) +[1] 0.9351018 1.0758567 +``` + + diff --git a/man/mnre_dim_and_class_to_index.Rd b/man/mnre_dim_and_class_to_index.Rd new file mode 100644 index 0000000..9986601 --- /dev/null +++ b/man/mnre_dim_and_class_to_index.Rd @@ -0,0 +1,25 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/RcppExports.R +\name{mnre_dim_and_class_to_index} +\alias{mnre_dim_and_class_to_index} +\title{Index from dimension and class indicies} +\usage{ +mnre_dim_and_class_to_index(i_dim, i_class, n_dim) +} +\arguments{ +\item{i_dim}{integer} + +\item{i_class}{integer} + +\item{n_dim}{integer} +} +\value{ +The integer index corresponding to the dimension and class indicies +} +\description{ +Index from dimension and class indicies +} +\examples{ +mnre_dim_and_class_to_index(1, 1, 1) +mnre_dim_and_class_to_index(2, 3, 5) +} diff --git a/man/mnre_expand_matrix.Rd b/man/mnre_expand_matrix.Rd new file mode 100644 index 0000000..55ed572 --- /dev/null +++ b/man/mnre_expand_matrix.Rd @@ -0,0 +1,21 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/RcppExports.R +\name{mnre_expand_matrix} +\alias{mnre_expand_matrix} +\title{Expand matrix} +\usage{ +mnre_expand_matrix(x1, k_class, direction) +} +\arguments{ +\item{x1}{sparse matrix} + +\item{k_class}{integer} + +\item{direction}{integer} +} +\value{ +Matrix expanded in the The integer index corresponding to the dimension and class indicies +} +\description{ +Used to duplicate the random effects coefficient matrix before applying the left factor +} diff --git a/man/mnre_fit_sparse.Rd b/man/mnre_fit_sparse.Rd new file mode 100644 index 0000000..04f141f --- /dev/null +++ b/man/mnre_fit_sparse.Rd @@ -0,0 +1,22 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/RcppExports.R +\name{mnre_fit_sparse} +\alias{mnre_fit_sparse} +\title{Multinomial mixed effects fit} +\usage{ +mnre_fit_sparse(fixed_effects, random_effects, y, theta_mat, Lind, beta_fixed, + beta_random, verbose = 1L) +} +\arguments{ +\item{x1}{sparse matrix} + +\item{k_class}{integer} + +\item{direction}{integer} +} +\value{ +Matrix expanded in the The integer index corresponding to the dimension and class indicies +} +\description{ +Used to duplicate the random effects coefficient matrix before applying the left factor +} diff --git a/man/mnre_make_covar.Rd b/man/mnre_make_covar.Rd new file mode 100644 index 0000000..04571f4 --- /dev/null +++ b/man/mnre_make_covar.Rd @@ -0,0 +1,28 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/RcppExports.R +\name{mnre_make_covar} +\alias{mnre_make_covar} +\title{Make covariance matrix} +\usage{ +mnre_make_covar(theta_mat, Lind, off_diagonal = 0) +} +\arguments{ +\item{theta_mat}{numeric matrix} + +\item{Lind}{unsigned integer matrix} + +\item{off_diagonal}{numeric} +} +\value{ +Covariance matrix for random effects +} +\description{ +Make covariance matrix +} +\examples{ +\dontrun{ +theta_mat <- matrix(1:4, ncol=2) +Lind <- matrix(c(rep(1, 10), rep(2, 5)), ncol=1) +covar_mat <- mnre_make_covar(theta_mat, Lind) +} +} diff --git a/man/nd_min_fun.Rd b/man/nd_min_fun.Rd new file mode 100644 index 0000000..38fc349 --- /dev/null +++ b/man/nd_min_fun.Rd @@ -0,0 +1,24 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/main.R +\name{nd_min_fun} +\alias{nd_min_fun} +\title{N-dimensional function generator} +\usage{ +nd_min_fun(ev) +} +\arguments{ +\item{ev}{list} +} +\value{ +deviance function +} +\description{ +N-dimensional function generator +} +\examples{ +\dontrun{ +ev = mnre_simulate_multinomial_data_factors(nfct=2, K_class = 2, nlev=50, nobs=20000) +nf <- nd_min_fun(ev) +nf(c(1,1)) +} +} diff --git a/mnre.Rproj b/mnre.Rproj new file mode 100644 index 0000000..21a4da0 --- /dev/null +++ b/mnre.Rproj @@ -0,0 +1,17 @@ +Version: 1.0 + +RestoreWorkspace: Default +SaveWorkspace: Default +AlwaysSaveHistory: Default + +EnableCodeIndexing: Yes +UseSpacesForTab: Yes +NumSpacesForTab: 2 +Encoding: UTF-8 + +RnwWeave: Sweave +LaTeX: pdfLaTeX + +BuildType: Package +PackageUseDevtools: Yes +PackageInstallArgs: --no-multiarch --with-keep.source diff --git a/src/.gitignore b/src/.gitignore new file mode 100644 index 0000000..22034c4 --- /dev/null +++ b/src/.gitignore @@ -0,0 +1,3 @@ +*.o +*.so +*.dll diff --git a/src/Makevars b/src/Makevars new file mode 100644 index 0000000..e69de29 diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp new file mode 100644 index 0000000..ecb5d8d --- /dev/null +++ b/src/RcppExports.cpp @@ -0,0 +1,202 @@ +// Generated by using Rcpp::compileAttributes() -> do not edit by hand +// Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 + +#include +#include +#include + +using namespace Rcpp; + +// mnre_dim_and_class_to_index +int mnre_dim_and_class_to_index(int i_dim, int i_class, int n_dim); +RcppExport SEXP _mnre_mnre_dim_and_class_to_index(SEXP i_dimSEXP, SEXP i_classSEXP, SEXP n_dimSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< int >::type i_dim(i_dimSEXP); + Rcpp::traits::input_parameter< int >::type i_class(i_classSEXP); + Rcpp::traits::input_parameter< int >::type n_dim(n_dimSEXP); + rcpp_result_gen = Rcpp::wrap(mnre_dim_and_class_to_index(i_dim, i_class, n_dim)); + return rcpp_result_gen; +END_RCPP +} +// mnre_make_covar +arma::sp_mat mnre_make_covar(const arma::mat& theta_mat, const arma::umat& Lind, double off_diagonal); +RcppExport SEXP _mnre_mnre_make_covar(SEXP theta_matSEXP, SEXP LindSEXP, SEXP off_diagonalSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const arma::mat& >::type theta_mat(theta_matSEXP); + Rcpp::traits::input_parameter< const arma::umat& >::type Lind(LindSEXP); + Rcpp::traits::input_parameter< double >::type off_diagonal(off_diagonalSEXP); + rcpp_result_gen = Rcpp::wrap(mnre_make_covar(theta_mat, Lind, off_diagonal)); + return rcpp_result_gen; +END_RCPP +} +// mnre_expand_matrix +arma::sp_mat mnre_expand_matrix(const arma::sp_mat& x1, int k_class, int direction); +RcppExport SEXP _mnre_mnre_expand_matrix(SEXP x1SEXP, SEXP k_classSEXP, SEXP directionSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const arma::sp_mat& >::type x1(x1SEXP); + Rcpp::traits::input_parameter< int >::type k_class(k_classSEXP); + Rcpp::traits::input_parameter< int >::type direction(directionSEXP); + rcpp_result_gen = Rcpp::wrap(mnre_expand_matrix(x1, k_class, direction)); + return rcpp_result_gen; +END_RCPP +} +// mnre_left_covar_factor +arma::sp_mat mnre_left_covar_factor(arma::sp_mat& x1); +RcppExport SEXP _mnre_mnre_left_covar_factor(SEXP x1SEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< arma::sp_mat& >::type x1(x1SEXP); + rcpp_result_gen = Rcpp::wrap(mnre_left_covar_factor(x1)); + return rcpp_result_gen; +END_RCPP +} +// mnre_fit_sparse +Rcpp::List mnre_fit_sparse(const arma::sp_mat& fixed_effects, const arma::sp_mat& random_effects, const arma::vec& y, const arma::mat& theta_mat, const arma::uvec& Lind, arma::mat beta_fixed, arma::mat beta_random, int verbose); +RcppExport SEXP _mnre_mnre_fit_sparse(SEXP fixed_effectsSEXP, SEXP random_effectsSEXP, SEXP ySEXP, SEXP theta_matSEXP, SEXP LindSEXP, SEXP beta_fixedSEXP, SEXP beta_randomSEXP, SEXP verboseSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const arma::sp_mat& >::type fixed_effects(fixed_effectsSEXP); + Rcpp::traits::input_parameter< const arma::sp_mat& >::type random_effects(random_effectsSEXP); + Rcpp::traits::input_parameter< const arma::vec& >::type y(ySEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type theta_mat(theta_matSEXP); + Rcpp::traits::input_parameter< const arma::uvec& >::type Lind(LindSEXP); + Rcpp::traits::input_parameter< arma::mat >::type beta_fixed(beta_fixedSEXP); + Rcpp::traits::input_parameter< arma::mat >::type beta_random(beta_randomSEXP); + Rcpp::traits::input_parameter< int >::type verbose(verboseSEXP); + rcpp_result_gen = Rcpp::wrap(mnre_fit_sparse(fixed_effects, random_effects, y, theta_mat, Lind, beta_fixed, beta_random, verbose)); + return rcpp_result_gen; +END_RCPP +} +// mnre_lk_penalty +double mnre_lk_penalty(const arma::mat& beta_random, const arma::mat& theta_norm, const arma::umat& Lind); +RcppExport SEXP _mnre_mnre_lk_penalty(SEXP beta_randomSEXP, SEXP theta_normSEXP, SEXP LindSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const arma::mat& >::type beta_random(beta_randomSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type theta_norm(theta_normSEXP); + Rcpp::traits::input_parameter< const arma::umat& >::type Lind(LindSEXP); + rcpp_result_gen = Rcpp::wrap(mnre_lk_penalty(beta_random, theta_norm, Lind)); + return rcpp_result_gen; +END_RCPP +} +// mnre_lk_glm +double mnre_lk_glm(const arma::sp_mat& fixed_effects, const arma::sp_mat& random_effects, const arma::mat& beta_fixed, const arma::mat& beta_random, const arma::vec& y, const arma::umat& Lind); +RcppExport SEXP _mnre_mnre_lk_glm(SEXP fixed_effectsSEXP, SEXP random_effectsSEXP, SEXP beta_fixedSEXP, SEXP beta_randomSEXP, SEXP ySEXP, SEXP LindSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const arma::sp_mat& >::type fixed_effects(fixed_effectsSEXP); + Rcpp::traits::input_parameter< const arma::sp_mat& >::type random_effects(random_effectsSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type beta_fixed(beta_fixedSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type beta_random(beta_randomSEXP); + Rcpp::traits::input_parameter< const arma::vec& >::type y(ySEXP); + Rcpp::traits::input_parameter< const arma::umat& >::type Lind(LindSEXP); + rcpp_result_gen = Rcpp::wrap(mnre_lk_glm(fixed_effects, random_effects, beta_fixed, beta_random, y, Lind)); + return rcpp_result_gen; +END_RCPP +} +// mnre_lk +double mnre_lk(const arma::sp_mat& fixed_effects, const arma::sp_mat& random_effects, const arma::mat& beta_fixed, const arma::mat& beta_random, const arma::vec& y, const arma::mat& theta_norm, const arma::umat& Lind); +RcppExport SEXP _mnre_mnre_lk(SEXP fixed_effectsSEXP, SEXP random_effectsSEXP, SEXP beta_fixedSEXP, SEXP beta_randomSEXP, SEXP ySEXP, SEXP theta_normSEXP, SEXP LindSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const arma::sp_mat& >::type fixed_effects(fixed_effectsSEXP); + Rcpp::traits::input_parameter< const arma::sp_mat& >::type random_effects(random_effectsSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type beta_fixed(beta_fixedSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type beta_random(beta_randomSEXP); + Rcpp::traits::input_parameter< const arma::vec& >::type y(ySEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type theta_norm(theta_normSEXP); + Rcpp::traits::input_parameter< const arma::umat& >::type Lind(LindSEXP); + rcpp_result_gen = Rcpp::wrap(mnre_lk(fixed_effects, random_effects, beta_fixed, beta_random, y, theta_norm, Lind)); + return rcpp_result_gen; +END_RCPP +} +// mnre_mu +arma::mat mnre_mu(const arma::mat& fixed_effects, const arma::sp_mat& random_effects, const arma::mat& beta_fixed, const arma::mat& beta_random); +RcppExport SEXP _mnre_mnre_mu(SEXP fixed_effectsSEXP, SEXP random_effectsSEXP, SEXP beta_fixedSEXP, SEXP beta_randomSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const arma::mat& >::type fixed_effects(fixed_effectsSEXP); + Rcpp::traits::input_parameter< const arma::sp_mat& >::type random_effects(random_effectsSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type beta_fixed(beta_fixedSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type beta_random(beta_randomSEXP); + rcpp_result_gen = Rcpp::wrap(mnre_mu(fixed_effects, random_effects, beta_fixed, beta_random)); + return rcpp_result_gen; +END_RCPP +} +// mnre_mu_x +arma::mat mnre_mu_x(const arma::sp_mat& fe_x, const arma::sp_mat& re_x, const arma::mat& beta_fixed, const arma::mat& beta_random); +RcppExport SEXP _mnre_mnre_mu_x(SEXP fe_xSEXP, SEXP re_xSEXP, SEXP beta_fixedSEXP, SEXP beta_randomSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const arma::sp_mat& >::type fe_x(fe_xSEXP); + Rcpp::traits::input_parameter< const arma::sp_mat& >::type re_x(re_xSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type beta_fixed(beta_fixedSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type beta_random(beta_randomSEXP); + rcpp_result_gen = Rcpp::wrap(mnre_mu_x(fe_x, re_x, beta_fixed, beta_random)); + return rcpp_result_gen; +END_RCPP +} +// mnre_step_sparse +Rcpp::List mnre_step_sparse(const arma::sp_mat& fixed_effects, const arma::sp_mat& random_effects, const arma::vec& y, const arma::mat& beta_fixed, const arma::mat& beta_random, const arma::mat& lambda_norm, const arma::uvec& Lind); +RcppExport SEXP _mnre_mnre_step_sparse(SEXP fixed_effectsSEXP, SEXP random_effectsSEXP, SEXP ySEXP, SEXP beta_fixedSEXP, SEXP beta_randomSEXP, SEXP lambda_normSEXP, SEXP LindSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const arma::sp_mat& >::type fixed_effects(fixed_effectsSEXP); + Rcpp::traits::input_parameter< const arma::sp_mat& >::type random_effects(random_effectsSEXP); + Rcpp::traits::input_parameter< const arma::vec& >::type y(ySEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type beta_fixed(beta_fixedSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type beta_random(beta_randomSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type lambda_norm(lambda_normSEXP); + Rcpp::traits::input_parameter< const arma::uvec& >::type Lind(LindSEXP); + rcpp_result_gen = Rcpp::wrap(mnre_step_sparse(fixed_effects, random_effects, y, beta_fixed, beta_random, lambda_norm, Lind)); + return rcpp_result_gen; +END_RCPP +} +// fill_mtwm_x +arma::sp_mat fill_mtwm_x(const arma::sp_mat& x1, const arma::sp_mat& x2, const arma::mat& mu); +RcppExport SEXP _mnre_fill_mtwm_x(SEXP x1SEXP, SEXP x2SEXP, SEXP muSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const arma::sp_mat& >::type x1(x1SEXP); + Rcpp::traits::input_parameter< const arma::sp_mat& >::type x2(x2SEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type mu(muSEXP); + rcpp_result_gen = Rcpp::wrap(fill_mtwm_x(x1, x2, mu)); + return rcpp_result_gen; +END_RCPP +} + +static const R_CallMethodDef CallEntries[] = { + {"_mnre_mnre_dim_and_class_to_index", (DL_FUNC) &_mnre_mnre_dim_and_class_to_index, 3}, + {"_mnre_mnre_make_covar", (DL_FUNC) &_mnre_mnre_make_covar, 3}, + {"_mnre_mnre_expand_matrix", (DL_FUNC) &_mnre_mnre_expand_matrix, 3}, + {"_mnre_mnre_left_covar_factor", (DL_FUNC) &_mnre_mnre_left_covar_factor, 1}, + {"_mnre_mnre_fit_sparse", (DL_FUNC) &_mnre_mnre_fit_sparse, 8}, + {"_mnre_mnre_lk_penalty", (DL_FUNC) &_mnre_mnre_lk_penalty, 3}, + {"_mnre_mnre_lk_glm", (DL_FUNC) &_mnre_mnre_lk_glm, 6}, + {"_mnre_mnre_lk", (DL_FUNC) &_mnre_mnre_lk, 7}, + {"_mnre_mnre_mu", (DL_FUNC) &_mnre_mnre_mu, 4}, + {"_mnre_mnre_mu_x", (DL_FUNC) &_mnre_mnre_mu_x, 4}, + {"_mnre_mnre_step_sparse", (DL_FUNC) &_mnre_mnre_step_sparse, 7}, + {"_mnre_fill_mtwm_x", (DL_FUNC) &_mnre_fill_mtwm_x, 3}, + {NULL, NULL, 0} +}; + +RcppExport void R_init_mnre(DllInfo *dll) { + R_registerRoutines(dll, NULL, CallEntries, NULL, NULL); + R_useDynamicSymbols(dll, FALSE); +} diff --git a/src/mnre.cpp b/src/mnre.cpp new file mode 100644 index 0000000..787276e --- /dev/null +++ b/src/mnre.cpp @@ -0,0 +1,704 @@ + + +#include +#include +#include "mnre.h" + +// [[Rcpp::depends(RcppArmadillo)]] +// [[Rcpp::depends(RcppEigen)]] + +using namespace arma; + +#define DEBUG 0 +#define EXPAND_ROW 0 +#define EXPAND_COLUMN 1 + +//' Index from dimension and class indicies +//' +//' @param i_dim integer +//' @param i_class integer +//' @param n_dim integer +//' @examples +//' mnre_dim_and_class_to_index(1, 1, 1) +//' mnre_dim_and_class_to_index(2, 3, 5) +//' @return The integer index corresponding to the dimension and class indicies +//' @export +// [[Rcpp::export]] +int mnre_dim_and_class_to_index(int i_dim, int i_class, int n_dim) { + return i_class * n_dim + i_dim; +} + +//' Make covariance matrix +//' +//' @param theta_mat numeric matrix +//' @param Lind unsigned integer matrix +//' @param off_diagonal numeric +//' @examples +//' \dontrun{ +//' theta_mat <- matrix(1:4, ncol=2) +//' Lind <- matrix(c(rep(1, 10), rep(2, 5)), ncol=1) +//' covar_mat <- mnre_make_covar(theta_mat, Lind) +//' } +//' @return Covariance matrix for random effects +//' @export +// [[Rcpp::export]] +arma::sp_mat mnre_make_covar(const arma::mat& theta_mat, + const arma::umat& Lind, + double off_diagonal = 0.0) { + + int nr = Lind.n_rows; + int nc = theta_mat.n_cols; + arma::sp_mat covar_mat; + + int max_idx = nr * nc * nc; + arma::umat covar_mat_idx(2, max_idx); + arma::vec covar_mat_v(max_idx); + int idx1, idx2, lidx; + double theta1, theta2; + int counter = 0; + double covar_val; + + for (int ir=0; ir < nr; ir++) { + lidx = Lind(ir)-1; + for (int k1=0; k1 < nc; k1++) { + theta1 = theta_mat(lidx, k1); + idx1 = mnre_dim_and_class_to_index(ir, k1, nr); + for (int k2=0; k2 < nc; k2++) { + theta2 = theta_mat(lidx, k2); + idx2 = mnre_dim_and_class_to_index(ir, k2, nr); + + covar_val = theta1 * theta2; + if (k1 != k2) { + covar_val *= off_diagonal; + } + + covar_mat_idx(0, counter) = idx1; + covar_mat_idx(1, counter) = idx2; + covar_mat_v(counter) = covar_val; + counter += 1; + + } + } + } + + covar_mat = arma::sp_mat(covar_mat_idx, covar_mat_v); + return covar_mat; +} + +//' Expand matrix +//' +//' Used to duplicate the random effects coefficient matrix before applying the left factor +//' @param x1 sparse matrix +//' @param k_class integer +//' @param direction integer +//' @return Matrix expanded in the The integer index corresponding to the dimension and class indicies +//' @export +// [[Rcpp::export]] +arma::sp_mat mnre_expand_matrix(const arma::sp_mat& x1, int k_class, int direction) { + + arma::sp_mat expanded_mat; + + for (int k1=0; k1 < k_class; k1++) { + if (k1 == 0) { + expanded_mat = x1; + } else { + if (direction == EXPAND_ROW) { + expanded_mat = arma::join_vert(expanded_mat, x1); + } else if (direction == EXPAND_COLUMN) { + expanded_mat = arma::join_horiz(expanded_mat, x1); + } else { + throw std::invalid_argument("direction must be EXPAND_ROW or EXPAND_COLUMN"); + } + } + } + return expanded_mat; +} + +// [[Rcpp::export]] +arma::sp_mat mnre_left_covar_factor(arma::sp_mat& x1) { + arma::sp_mat left_factor; + + SimplicialLLT > solver; + + SparseMatrix rx = Rcpp::as >(wrap(x1)); + solver.compute(rx); + left_factor = Rcpp::as(wrap(solver.matrixL())); + return left_factor; +} + + +//' Multinomial mixed effects fit +//' +//' Used to duplicate the random effects coefficient matrix before applying the left factor +//' @param x1 sparse matrix +//' @param k_class integer +//' @param direction integer +//' @return Matrix expanded in the The integer index corresponding to the dimension and class indicies +//' @export +// [[Rcpp::export]] +Rcpp::List mnre_fit_sparse(const arma::sp_mat& fixed_effects, + const arma::sp_mat& random_effects, + const arma::vec& y, + const arma::mat& theta_mat, + const arma::uvec& Lind, + arma::mat beta_fixed, + arma::mat beta_random, + int verbose=1) { + arma::mat lambda_ones(theta_mat.n_rows, theta_mat.n_cols); + lambda_ones.fill(1); + + int n_data = fixed_effects.n_rows; + int n_dim_fixed = fixed_effects.n_cols; + int n_dim_random = random_effects.n_cols; + int k_class = (int) y.max(); + int n_dim = n_dim_fixed + n_dim_random; + + int D_times_K = n_dim * k_class; + + // arma::mat beta_fixed = arma::mat(n_dim_fixed, k_class); + // arma::mat beta_random = arma::mat(n_dim_random, k_class); + // beta_fixed.fill(0); + // beta_random.fill(0); + + int step_counter; + double tol = 1e-5; + double grad_tol = 1e-5; + + double loglk_det; + double loglk_det_val; + double loglk_det_sign; + + Rcpp::List ll; + int half_steps = 128; + double newton_step_factor = 1.0; + double lk_diff_stat, loglk_old, lk_new, beta_diff; + int max_step = 200; + + arma::sp_mat covar_mat = mnre_make_covar(theta_mat, Lind, 0.0); + arma::mat mu; + arma::sp_mat left_factor = mnre_left_covar_factor(covar_mat); + + arma::sp_mat re_x = mnre_expand_matrix(random_effects, k_class, EXPAND_COLUMN); + arma::sp_mat fe_x = mnre_expand_matrix(fixed_effects, k_class, EXPAND_COLUMN); + + if (DEBUG) { + Rcpp::Rcout << " left_factor " << left_factor.n_rows << " " << left_factor.n_cols << std::endl; + Rcpp::Rcout << " random_effects " << random_effects.n_rows << " " << random_effects.n_cols << std::endl; + Rcpp::Rcout << " re_x " << re_x.n_rows << " " << re_x.n_cols << " " << re_x(1, 2) << std::endl; + } + + arma::sp_mat ZLam = re_x * left_factor; + + if (DEBUG) { + Rcpp::Rcout << " ZLam " << ZLam.n_rows << " " << ZLam.n_cols << " " << ZLam(1, 2) << std::endl; + } + + arma::sp_mat rtwr; + arma::sp_mat sp_diag_mat = speye(n_dim_random * k_class, n_dim_random * k_class); + SimplicialLLT > solver; + SparseMatrix rx; + + for (int istep=0; istep < max_step; istep++) { + if (DEBUG) { + Rcpp:Rcout << " istep " << istep << std::endl; + } + + mu = mnre_mu_x(fe_x, ZLam, beta_fixed, beta_random); + + if (DEBUG) { + Rcpp::Rcout << " call step... " << std::endl; + } + + ll = mnre_step_sparse(fe_x, ZLam, + y, + beta_fixed, beta_random, + theta_mat, Lind); + + if (DEBUG) { + Rcpp::Rcout << " call step... DONE " << std::endl; + } + + loglk_old = (double)(ll["loglk_old"]); + lk_diff_stat = (double)(ll["loglk"]) - loglk_old; + + if (verbose >=1) { + Rcpp::Rcout << " istep " << istep << + " lk_diff_stat " << lk_diff_stat << + " lk value " << (double)ll["loglk"] << + " grad_check " << (double)ll["grad_check"] << std::endl; + } + + if (lk_diff_stat > 0) { // fit got worse + step_counter = 0; + while( (step_counter < half_steps) && ( lk_diff_stat >= 0 ) ) { + newton_step_factor = 0.5 * newton_step_factor; + lk_new = + mnre_lk(fe_x, ZLam, + beta_fixed + + newton_step_factor * as(ll["beta_fixed_diff"]), + beta_random + + newton_step_factor * as(ll["beta_random_diff"]), + y, lambda_ones, Lind); + lk_diff_stat = lk_new - loglk_old; + step_counter += 1; + + if (DEBUG) { + Rcpp::Rcout << " istep " << istep << " half-step " << step_counter << + " lk_diff_stat " << lk_diff_stat << + " lk_diff_stat / step " << lk_diff_stat / newton_step_factor << + " grad_check " << (double)ll["grad_check"] << std::endl; + } + } + } + + Rcpp::List ans; + if (lk_diff_stat >= 0) { // fit got worse + rtwr = fill_mtwm_x(ZLam, ZLam, mu); + rtwr = rtwr + sp_diag_mat; + + // TODO: either use the left factor or use rtwr + rx = Rcpp::as >(wrap(rtwr)); + solver.compute(rx); + left_factor = Rcpp::as(wrap(solver.matrixL())); + //loglk_det = arma::sum(std::log(left_factor.diag())); + //loglk_det = 0; + arma::log_det(loglk_det_val, loglk_det_sign, arma::mat(rtwr)); + loglk_det = loglk_det_val * loglk_det_sign; + + if (verbose >= 1) { + Rcpp::Rcout << "half step cycle did not find a better fit" << std::endl; + } + + ans["beta_fixed"] = beta_fixed; + ans["beta_random"] = beta_random; + ans["loglk"] = (double)(ll["loglk"]); + ans["loglk_det"] = loglk_det; + + return ans; + } else { + beta_fixed = beta_fixed + newton_step_factor * as(ll["beta_fixed_diff"]); + beta_random = beta_random + newton_step_factor * as(ll["beta_random_diff"]); + newton_step_factor = 1.0; + } + + if (std::abs(lk_diff_stat) <= tol && (double)ll["grad_check"] <= grad_tol) { + rtwr = fill_mtwm_x(ZLam, ZLam, mu); + rtwr = rtwr + sp_diag_mat; + // TODO: either use the left factor or use rtwr + rx = Rcpp::as >(wrap(rtwr)); + solver.compute(rx); + left_factor = Rcpp::as(wrap(solver.matrixL())); + //loglk_det = arma::sum(std::log(left_factor.diag())); + arma::log_det(loglk_det_val, loglk_det_sign, arma::mat(rtwr)); + loglk_det = loglk_det_val * loglk_det_sign; + + ans["beta_fixed"] = beta_fixed; + ans["beta_random"] = beta_random; + ans["loglk"] = (double)(ll["loglk"]); + ans["loglk_det"] = loglk_det; + return ans; + } + + } // steps + +} // end function + + +// [[Rcpp::export]] +double mnre_lk_penalty(const arma::mat& beta_random, + const arma::mat& theta_norm, + const arma::umat& Lind) { + + int k_class = theta_norm.n_cols; + double tmp; + double lk_penalty = 0.0; + int ind; + + for (int k1 = 0; k1 < k_class; k1++) { + tmp = 0.0; + for (int d1 = 0; d1 < beta_random.n_rows; d1++) { // dont penalize the intercept + ind = Lind(d1)-1; + tmp += beta_random(d1, k1) * beta_random(d1, k1) * theta_norm(ind, k1); + } + lk_penalty += tmp; + } + + return lk_penalty; + +} + + +// [[Rcpp::export]] +double mnre_lk_glm(const arma::sp_mat& fixed_effects, + const arma::sp_mat& random_effects, + const arma::mat& beta_fixed, + const arma::mat& beta_random, + const arma::vec& y, + const arma::umat& Lind) { + + int n_data = fixed_effects.n_rows; + int n_dim_fixed = fixed_effects.n_cols; + int n_dim_random = random_effects.n_cols; + int k_class = beta_fixed.n_cols; + + arma::mat mu = mnre_mu_x(fixed_effects, random_effects, + beta_fixed, beta_random); + + double log_lk = 0.0, lk_term; + int iy; + + for (int i=0; i < n_data; i++) { + iy = (int) y(i); + if (iy == 0) { + lk_term = std::log(1 - arma::sum(mu.row(i))); + } else { + lk_term = std::log(mu(i,iy-1)); + } + log_lk = log_lk + lk_term; + } + + return -2 * log_lk; +} + +// [[Rcpp::export]] +double mnre_lk(const arma::sp_mat& fixed_effects, + const arma::sp_mat& random_effects, + const arma::mat& beta_fixed, + const arma::mat& beta_random, + const arma::vec& y, + const arma::mat& theta_norm, + const arma::umat& Lind) { + + return mnre_lk_glm(fixed_effects, random_effects, + beta_fixed, beta_random, + y, Lind) + mnre_lk_penalty(beta_random, theta_norm, Lind); +} + + +// [[Rcpp::export]] +arma::mat mnre_mu(const arma::mat &fixed_effects, + const arma::sp_mat &random_effects, + const arma::mat &beta_fixed, + const arma::mat &beta_random) { + + int n_data = fixed_effects.n_rows; + int k_class = beta_fixed.n_cols; + arma::mat mu(n_data, k_class); + arma::mat eta = fixed_effects * beta_fixed + random_effects * beta_random; + arma::mat ee = arma::exp(eta); + + double denom; + for (int i=0; i < n_data; i++) { + denom = 1.0; + for (int k=0; k < k_class; k++) { + denom += ee(i, k); + } + + for (int k=0; k < k_class; k++) { + mu(i, k) = ee(i, k)/denom; + } + } + + return mu; +} + +/** This mu function works with the coordinates that have been expanded in the column direction + i.e. the coordinates are naturally N x D, but get expanded to N x D*k + **/ + +// [[Rcpp::export]] +arma::mat mnre_mu_x(const arma::sp_mat &fe_x, + const arma::sp_mat &re_x, + const arma::mat &beta_fixed, + const arma::mat &beta_random) { + + int n_data = fe_x.n_rows; + int k_class = beta_fixed.n_cols; + + int ndim_fe = fe_x.n_cols / k_class; + int ndim_re = re_x.n_cols / k_class; + + arma::sp_mat fixed_effects(n_data, ndim_fe); + arma::sp_mat random_effects(n_data, ndim_re); + + arma::mat mu(n_data, k_class); + arma::mat eta(n_data, k_class); + arma::mat ee(n_data, k_class); + + arma::vec denom(n_data); + + for (int k=0; k < k_class; k++) { + fixed_effects = fe_x.submat(0, ndim_fe * k, n_data-1, ndim_fe * (k+1) - 1); + random_effects = re_x.submat(0, ndim_re * k, n_data-1, ndim_re * (k+1) - 1); + eta.col(k) = fixed_effects * beta_fixed.col(k) + random_effects * beta_random.col(k); + ee.col(k) = exp(eta.col(k)); + } + + for (int i=0; i< n_data; i++) { + denom(i) = 1.0; + for (int k=0; k < k_class; k++) { + denom(i) += ee(i, k); + } + + for (int k=0; k < k_class; k++) { + mu(i, k) = ee(i, k)/denom(i); + } + } + + return mu; +} + +// [[Rcpp::export]] +Rcpp::List mnre_step_sparse(const arma::sp_mat &fixed_effects, + const arma::sp_mat &random_effects, + const arma::vec &y, + const arma::mat &beta_fixed, + const arma::mat &beta_random, + const arma::mat &lambda_norm, + const arma::uvec &Lind) { + // The left factor of the covariance should already be applied to the random effects + // TODO: if so, then lambda_norm must be 1; should that be enforced? + + arma::mat lambda_ones(lambda_norm.n_rows, lambda_norm.n_cols); + lambda_ones.fill(1); + + int log_counter = 0; + + int n_data = y.n_elem; + int k_class = arma::max(y); + + int n_dim_fixed = beta_fixed.n_rows; + int n_dim_random = beta_random.n_rows; + int n_dim = n_dim_fixed + n_dim_random; + int D_times_K = n_dim * k_class; + int Dfixed_times_K = n_dim_fixed * k_class; + int Drandom_times_K = n_dim_random * k_class; + + // check dimensions + + // if (random_effects.n_rows != n_data) { + // Rcpp::Rcout << " error: " << random_effects.n_rows << " " << n_data; + // throw std::runtime_error("random_effects.n_rows != n_data"); + // } + // + // if (random_effects.n_cols != n_dim_random * k_class) { + // Rcpp::Rcout << " error: " << random_effects.n_cols << " " << n_dim_random * k_class; + // throw std::runtime_error("random_effects.n_cols != n_dim_random * k_class"); + // } + // + // if (fixed_effects.n_rows != n_data) { + // Rcpp::Rcout << " error: " << fixed_effects.n_rows << " " << n_data; + // throw std::runtime_error("fixed_effects.n_rows != n_data"); + // } + // + // if (fixed_effects.n_cols != n_dim_fixed * k_class) { + // Rcpp::Rcout << " error: " << fixed_effects.n_cols << " " << n_dim_fixed * k_class; + // throw std::runtime_error("fixed_effects.n_cols != n_dim_fixed * k_class"); + // } + + arma::mat mu = mnre_mu_x(fixed_effects, random_effects, beta_fixed, beta_random); + + arma::mat dy(n_data, k_class); + arma::sp_mat sp_diag_mat = speye(n_dim_random * k_class, n_dim_random * k_class); + int idx1, idx2; + + arma::sp_mat ftwf = fill_mtwm_x(fixed_effects, fixed_effects, mu); + arma::sp_mat ftwr = fill_mtwm_x(fixed_effects, random_effects, mu); + arma::sp_mat rtwr = fill_mtwm_x(random_effects, random_effects, mu); + + if (DEBUG) { + Rcpp::Rcout << " mnre step " << + " ndimfixed " << n_dim_fixed << " ndimrandom " << n_dim_random << + " ftwf " << ftwf.n_rows << " " << ftwf.n_cols << + " ftwr " << ftwr.n_rows << " " << ftwr.n_cols << + " rtwr " << rtwr.n_rows << " " << rtwr.n_cols << + std::endl; + } + + // TODO: double check that the sp_diag_mat dimensions are correct, dont just fudge it like this + //rtwr = rtwr + sp_diag_mat; + rtwr = rtwr + speye(rtwr.n_rows, rtwr.n_cols); + + arma::sp_mat lhs_mat = arma::join_vert( + arma::join_horiz(ftwf, ftwr), + arma::join_horiz(ftwr.t(), rtwr) + ); + + if (DEBUG) { + Rcpp::Rcout << " mnre step " << + " mu " << mu.n_rows << " " << mu.n_cols << + " dy " << dy.n_rows << " " << dy.n_cols << std::endl; + } + + int iy; + for (int i=0; i < y.n_rows; i++) { + iy = (int)(y(i)); + for (int k=0; k < k_class; k++) { + if (iy == k+1) { + dy(i, k) = (1 - mu(i, k)); + } else { + dy(i, k) = (0 - mu(i, k)); + } + } + } + + arma::mat rhs_fe = arma::trans(fixed_effects.submat(0, 0, n_data-1, n_dim_fixed-1)) * dy; + arma::mat rhs_re = arma::trans(random_effects.submat(0, 0, n_data-1, n_dim_random-1)) * dy; + + if (DEBUG) { + Rcpp::Rcout << + " n_dim_random " << n_dim_random << + " rhs_re " << rhs_re.n_rows << " " << rhs_re.n_cols << + std::endl; + + Rcpp::Rcout << + " n_dim_fixed " << n_dim_fixed << + " rhs_fe " << rhs_fe.n_rows << " " << rhs_fe.n_cols << + std::endl; + } + + // for the penalty term. Note this assumes spherical random effects + for (int d1=0; d1 < n_dim_random; d1++) { + for (int k1=0; k1 < k_class; k1++) { + rhs_re(d1, k1) -= beta_random(d1, k1); + } + } + + arma::mat rhs_mat = arma::join_vert(arma::vectorise(rhs_fe), + arma::vectorise(rhs_re)); + + + double grad_check = arma::norm(arma::vectorise(rhs_mat % rhs_mat)); + + if (DEBUG) { + Rcpp::Rcout << " mnre step " << + " rhs_mat " << rhs_mat.n_rows << " " << rhs_mat.n_cols << + " lhs_mat " << lhs_mat.n_rows << " " << lhs_mat.n_cols << + std::endl; + } + + // SparseQR, COLAMDOrdering > solver; + + SimplicialLDLT > solver; + + SparseMatrix rx = Rcpp::as >(wrap(lhs_mat)); + MatrixXd ry = Rcpp::as(wrap(rhs_mat)); + + if (DEBUG) { + Rcpp:Rcout << "solve the linear system..." << std::endl; + } + + solver.compute(rx); + MatrixXd vsol = solver.solve(ry); + + //arma::vec vsol = spsolve(lhs_mat, rhs_mat, "lapack"); + + arma::mat db = arma::reshape( + as(wrap(vsol.block(0, 0, Dfixed_times_K, 1))), + n_dim_fixed, k_class + ); + + arma::mat du = arma::reshape( + as(wrap(vsol.block(Dfixed_times_K, 0, Drandom_times_K, 1))), + n_dim_random, k_class + ); + + double loglk_old = mnre_lk(fixed_effects, random_effects, + beta_fixed, beta_random, + y, lambda_ones, Lind); + + double loglk = mnre_lk(fixed_effects, random_effects, + beta_fixed + db, beta_random + du, + y, lambda_ones, Lind); + + if (DEBUG) { + Rcpp::Rcout << " mnre step " << + " beta_fixed " << beta_fixed.n_rows << " " << beta_fixed.n_cols << + " beta_random " << beta_random.n_rows << " " << beta_random.n_cols << + " db " << db.n_rows << " " << db.n_cols << + " du " << du.n_rows << " " << du.n_cols << + std::endl; + } + + arma::mat beta_mat_old = arma::join_vert(beta_fixed, beta_random); + arma::mat beta_mat = arma::join_vert(beta_fixed+db, beta_random+du); + Rcpp::List ans; + ans["loglk"] = loglk; + ans["loglk_old"] = loglk_old; + ans["beta_diff"] = db; + ans["beta_mat_old"] = beta_mat_old; + ans["beta_mat"] = beta_mat; + ans["beta_fixed_old"] = beta_fixed; + ans["beta_fixed"] = beta_fixed + db; + ans["beta_fixed_diff"] = db; + ans["beta_random_old"] = beta_random; + ans["beta_random"] = beta_random + du; + ans["beta_random_diff"] = du; + ans["grad_check"] = grad_check; + + return ans; +} + + + +// [[Rcpp::export]] +arma::sp_mat fill_mtwm_x(const arma::sp_mat& x1, const arma::sp_mat& x2, + const arma::mat& mu) { + + int Kclass = mu.n_cols; + + int D1 = x1.n_cols / Kclass; + int D2 = x2.n_cols / Kclass; + + int D1_times_K = D1 * Kclass; + int D2_times_K = D2 * Kclass; + arma::sp_mat mtwm(D1_times_K, D2_times_K); + arma::vec ww; + int idx1, idx2; + double aa; + arma::sp_mat tmp_mat; + int idx11, idx12, idx21, idx22; + arma::sp_mat ww_sp; + int tmp; + + arma::mat ww_con1(x1.n_rows, D1); + arma::mat ww_con2(x2.n_rows, D2); + ww_con1.fill(1); + ww_con2.fill(1); + arma::mat ww_mat(mu.n_rows, Kclass); + + for (int k1=0; k1 < Kclass; k1++) { + for (int k2=0; k2 < Kclass; k2++) { + + if (k1 == k2) { + ww = mu.col(k1) % (1-mu.col(k1)); + } else { + ww = mu.col(k1) % mu.col(k2); + } + ww = arma::sqrt(ww); + for (int i=0; i +#include + +#ifndef mnre_H +#define mnre_H + +using Eigen::ArrayXd; +using Eigen::LLT; +using Eigen::MatrixXd; +using Eigen::VectorXd; + +using Eigen::SparseMatrix; +using Eigen::SparseQR; + +using namespace Rcpp; +using namespace Eigen; + +class mnre { + int K_class; + arma::sp_mat fixed_effects; + arma::sp_mat random_effects; + arma::mat y; + arma::mat theta_mat; + +}; + +arma::sp_mat fill_mtwm_x(const arma::sp_mat& x1, const arma::sp_mat& x2, + const arma::mat& mu); + +Rcpp::List mnre_fit_sparse(const arma::sp_mat& fixed_effects, + const arma::sp_mat& random_effects, + const arma::vec& y, + const arma::mat& theta_norm, + const arma::uvec& Lind, + arma::mat beta_fixed, + arma::mat beta_random, + int verbose); + +Rcpp::List mnre_step_sparse(const arma::sp_mat &fixed_effects, + const arma::sp_mat &random_effects, + const arma::vec &y, + const arma::mat &beta_fixed, + const arma::mat &beta_random, + const arma::mat &lambda_norm, + const arma::uvec &Lind); + +double mnre_lk_glm(const arma::sp_mat& fixed_effects, + const arma::sp_mat& random_effects, + const arma::mat& beta_fixed, + const arma::mat& beta_random, + const arma::vec& y, + const arma::umat& Lind); + +double mnre_lk(const arma::sp_mat& fixed_effects, + const arma::sp_mat& random_effects, + const arma::mat& beta_fixed, + const arma::mat& beta_random, + const arma::vec &y, + const arma::mat &lambda_norm, + const arma::umat& Lind); + +arma::mat mnre_mu_x(const arma::sp_mat &fe_x, + const arma::sp_mat &re_x, + const arma::mat &beta_fixed, + const arma::mat &beta_random); + +arma::mat mnre_mu(const arma::mat &fixed_effects, + const arma::sp_mat &random_effects, + const arma::mat &beta_fixed, + const arma::mat &beta_random); + +arma::sp_mat mnre_expand_matrix(const arma::sp_mat& x1, int k_class, int directions); + +#endif