diff --git a/DESCRIPTION b/DESCRIPTION index d88d139..f47e078 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,8 +1,8 @@ -Package: GMVJM +Package: gmvjoint Type: Package Title: Joint models of survival and multivariate longitudinal data -Version: 0.5 -Date: 2022-12-02 +Version: 0.1 +Date: 2022-12-05 Author: James Murray Maintainer: James Murray Description: Fit joint models of survival and multivariate longitudinal data. The longitudinal @@ -26,4 +26,4 @@ LinkingTo: Rcpp, RcppArmadillo Encoding: UTF-8 RoxygenNote: 7.2.2 LazyData: true -URL: https://github.com/jamesmurray7/GMVJM +URL: https://github.com/jamesmurray7/gmvjoint diff --git a/NAMESPACE b/NAMESPACE index 95a6520..09291bc 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,4 +1,4 @@ -useDynLib(GMVJM, .registration=TRUE) +useDynLib(gmvjoint, .registration=TRUE) importFrom(Rcpp, evalCpp) importFrom("methods", "el") importFrom("stats", "Gamma", "as.formula", "binomial", "coef", "cov", diff --git a/NEWS.md b/NEWS.md index b19f696..c149d69 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,2 +1,2 @@ -# GMVJM version 0.1 +# gmvjoint version 0.1 * First release diff --git a/R/RcppExports.R b/R/RcppExports.R index 83c5b8b..0ef214a 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -4,90 +4,90 @@ #' Specifically for obtaining pmf in simulations (since requisite doesn't exist in R). #' @keywords internal GP1_pmf_scalar <- function(mu, phi, Y) { - .Call(`_GMVJM_GP1_pmf_scalar`, mu, phi, Y) + .Call(`_gmvjoint_GP1_pmf_scalar`, mu, phi, Y) } #' log-likelihood of Gamma #' @keywords internal ll_Gamma <- function(Y, shape, mu) { - .Call(`_GMVJM_ll_Gamma`, Y, shape, mu) + .Call(`_gmvjoint_ll_Gamma`, Y, shape, mu) } #' 'GP1' from Zamani & Ismail (2012) #' @keywords internal ll_genpois <- function(eta, phi, Y) { - .Call(`_GMVJM_ll_genpois`, eta, phi, Y) + .Call(`_gmvjoint_ll_genpois`, eta, phi, Y) } #' Survival log-density #' @keywords internal logfti <- function(b, S, SS, Fi, Fu, l0i, haz, Delta, gamma_rep, zeta) { - .Call(`_GMVJM_logfti`, b, S, SS, Fi, Fu, l0i, haz, Delta, gamma_rep, zeta) + .Call(`_gmvjoint_logfti`, b, S, SS, Fi, Fu, l0i, haz, Delta, gamma_rep, zeta) } #' The joint density #' @keywords internal joint_density <- function(b, Y, X, Z, beta, D, sigma, family, Delta, S, Fi, l0i, SS, Fu, haz, gamma_rep, zeta, beta_inds, b_inds, K) { - .Call(`_GMVJM_joint_density`, b, Y, X, Z, beta, D, sigma, family, Delta, S, Fi, l0i, SS, Fu, haz, gamma_rep, zeta, beta_inds, b_inds, K) + .Call(`_gmvjoint_joint_density`, b, Y, X, Z, beta, D, sigma, family, Delta, S, Fi, l0i, SS, Fu, haz, gamma_rep, zeta, beta_inds, b_inds, K) } #' Quadrature - standard deviation of N(mu, tau^2). #' @keywords internal maketau <- function(S, Z) { - .Call(`_GMVJM_maketau`, S, Z) + .Call(`_gmvjoint_maketau`, S, Z) } #' First derivative of the joint density with respect to b. #' @keywords internal joint_density_ddb <- function(b, Y, X, Z, beta, D, sigma, family, Delta, S, Fi, l0i, SS, Fu, haz, gamma_rep, zeta, beta_inds, b_inds, K) { - .Call(`_GMVJM_joint_density_ddb`, b, Y, X, Z, beta, D, sigma, family, Delta, S, Fi, l0i, SS, Fu, haz, gamma_rep, zeta, beta_inds, b_inds, K) + .Call(`_gmvjoint_joint_density_ddb`, b, Y, X, Z, beta, D, sigma, family, Delta, S, Fi, l0i, SS, Fu, haz, gamma_rep, zeta, beta_inds, b_inds, K) } #' Create vector of scores on fixed effects. #' @keywords internal Sbeta <- function(beta, X, Y, Z, b, sigma, family, beta_inds, K, quad, tau, w, v) { - .Call(`_GMVJM_Sbeta`, beta, X, Y, Z, b, sigma, family, beta_inds, K, quad, tau, w, v) + .Call(`_gmvjoint_Sbeta`, beta, X, Y, Z, b, sigma, family, beta_inds, K, quad, tau, w, v) } #' Hessian matrix on fixed effects. #' @keywords internal Hbeta <- function(beta, X, Y, Z, b, sigma, family, beta_inds, K, quad, tau, w, v) { - .Call(`_GMVJM_Hbeta`, beta, X, Y, Z, b, sigma, family, beta_inds, K, quad, tau, w, v) + .Call(`_gmvjoint_Hbeta`, beta, X, Y, Z, b, sigma, family, beta_inds, K, quad, tau, w, v) } #' Update for residual variance (sigma for gaussian family) #' @keywords internal vare_update <- function(X, Y, Z, b, beta, tau, w, v) { - .Call(`_GMVJM_vare_update`, X, Y, Z, b, beta, tau, w, v) + .Call(`_gmvjoint_vare_update`, X, Y, Z, b, beta, tau, w, v) } #' Update for dispersion parameter (sigma for \code{"genpois"} family) #' @keywords internal phi_update <- function(b, X, Y, Z, beta, phi, w, v, tau) { - .Call(`_GMVJM_phi_update`, b, X, Y, Z, beta, phi, w, v, tau) + .Call(`_gmvjoint_phi_update`, b, X, Y, Z, beta, phi, w, v, tau) } #' Score vector for survival parameters. #' @keywords internal Sgammazeta <- function(gammazeta, b, Sigma, S, SS, Fu, Fi, haz, Delta, w, v, b_inds, K, q, eps) { - .Call(`_GMVJM_Sgammazeta`, gammazeta, b, Sigma, S, SS, Fu, Fi, haz, Delta, w, v, b_inds, K, q, eps) + .Call(`_gmvjoint_Sgammazeta`, gammazeta, b, Sigma, S, SS, Fu, Fi, haz, Delta, w, v, b_inds, K, q, eps) } #' Hessian matrix for survival parameters. #' @keywords internal Hgammazeta <- function(gammazeta, b, Sigma, S, SS, Fu, Fi, haz, Delta, w, v, b_inds, K, q, eps) { - .Call(`_GMVJM_Hgammazeta`, gammazeta, b, Sigma, S, SS, Fu, Fi, haz, Delta, w, v, b_inds, K, q, eps) + .Call(`_gmvjoint_Hgammazeta`, gammazeta, b, Sigma, S, SS, Fu, Fi, haz, Delta, w, v, b_inds, K, q, eps) } #' Update to the baseline hazard #' @keywords internal lambdaUpdate <- function(survtimes, ft, gamma, gamma_rep, zeta, S, Sigma, b, w, v, b_inds, K, q) { - .Call(`_GMVJM_lambdaUpdate`, survtimes, ft, gamma, gamma_rep, zeta, S, Sigma, b, w, v, b_inds, K, q) + .Call(`_gmvjoint_lambdaUpdate`, survtimes, ft, gamma, gamma_rep, zeta, S, Sigma, b, w, v, b_inds, K, q) } #' Second derivative of joint density wrt b #' @keywords internal joint_density_sdb <- function(b, Y, X, Z, beta, D, sigma, family, Delta, S, Fi, l0i, SS, Fu, haz, gamma_rep, zeta, beta_inds, b_inds, K, eps) { - .Call(`_GMVJM_joint_density_sdb`, b, Y, X, Z, beta, D, sigma, family, Delta, S, Fi, l0i, SS, Fu, haz, gamma_rep, zeta, beta_inds, b_inds, K, eps) + .Call(`_gmvjoint_joint_density_sdb`, b, Y, X, Z, beta, D, sigma, family, Delta, S, Fi, l0i, SS, Fu, haz, gamma_rep, zeta, beta_inds, b_inds, K, eps) } diff --git a/README.md b/README.md index dc202bc..15daaf8 100644 --- a/README.md +++ b/README.md @@ -1,13 +1,13 @@ -# `GMVJM` +# `gmvjoint` -## What is `GMVJM`? +## What is `gmvjoint`? -* "**G**eneralised" -* **M**ulti**v**ariate -* **J**oint **M**odels. +* "**g**eneralised" +* **m**ulti**v**ariate +* **Joint** models. -GMVJM allows the user to fit joint models of survival and multivariate longitudinal data, where the +`gmvjoint` allows the user to fit joint models of survival and multivariate longitudinal data, where the longitudinal sub-models are specified by generalised linear mixed models (GLMMs). The joint models are fit via maximum likelihood using an approximate EM algorithm first proposed by Bernhardt *et al*. (2015). The GLMMs are specified using the same syntax as for package `glmmTMB` (Brooks *et @@ -15,6 +15,8 @@ al*., 2017). The joint models themselves are then the flexible extensions to th Wulfoshn and Tsiatis (1997). The user is able to simulate data under many different response types. +Currently, five families can be fit: Gaussian; Poisson; binomial; Gamma and generalised Poisson. + ## To-do list The package in current incantation is relatively skeletal, as such not a lot of post-hoc anaylses on fitted joint models is possible. As such, an immediate to-do list currently looks like @@ -26,6 +28,37 @@ which just needs to be ported over. No promises are made w.r.t timescale of these being implemented: Currently I am a PhD student and little cache is awarded for production or maintenance of R packages! +## Example +To fit a joint model, we first need to specify the longitudinal and survival sub-models. + +The longitudinal sub-model **must** be a list which contains the specification of the longitudinal process along with its random effects structure +in the same syntax as a [glmmTMB](https://cran.r-project.org/web/packages/glmmTMB/index.html) model (which itself is the same as the widely-used `lme4`). +As an example, suppose we want to fit a trivariate model on the oft-used PBC data, with a linear time-drug interaction term on albumin, a spline term on +(logged) serum bilirubin and a linear fit on spiders, we specify +```r +PBC$serBilir <- log(PBC$serBilir) +long.formulas <- list( + albumin ~ drug * time + (1 + time|id), + serBilir ~ drug * splines::ns(time, 3) + (1 + splines::ns(time, 3)|id), + spiders ~ drug * time + (1|id) +) +``` +where we note interactions and spline-time fits are possible. Currently, transformations on variables (e.g. `log(Y)`) must be done *before* this setup of formulae. + +The survival sub-model must be set-up using `Surv()` from the [survival](https://cran.r-project.org/web/packages/survival/) package e.g. +```r +surv.formula <- Surv(survtime, status) ~ drug +``` +Currently interaction terms in the survival sub-model specification are unsupported. + +Now we can do the joint model call through the main workhorse function `joint`. This notably take a *list* of family arguments which **must** match-up in the desired order as the longitudinal process +list. We call our `fit` via +```r +fit <- joint(long.formulas = long.formulas, surv.formula = surv.formula, data = PBC, + family = list("gaussian", "gaussian", "binomial")) +``` +where extra control arguments are documented in `?joint`. + ## References Bernhardt PW, Zhang D and Wang HJ. A fast EM Algorithm for Fitting Joint Models of a Binary diff --git a/GMVJM.Rproj b/gmvjoint.Rproj similarity index 100% rename from GMVJM.Rproj rename to gmvjoint.Rproj diff --git a/man/GMVJM-package.Rd b/man/gmvjoint-package.Rd similarity index 87% rename from man/GMVJM-package.Rd rename to man/gmvjoint-package.Rd index 7356948..461dce2 100644 --- a/man/GMVJM-package.Rd +++ b/man/gmvjoint-package.Rd @@ -1,12 +1,12 @@ -\name{GMVJM-package} -\alias{GMVJM-package} -\alias{GMVJM} +\name{gmvjoint-package} +\alias{gmvjoint-package} +\alias{gmvjoint} \docType{package} \title{ -\packageTitle{GMVJM} +\packageTitle{gmvjoint} } \description{ -GMVJM allows the user to fit joint models of survival and multivariate longitudinal data. The +gmvjoint allows the user to fit joint models of survival and multivariate longitudinal data. The longitudinal data is specified by generalised linear mixed models (GLMMs). The joint models are fit via maximum likelihood using an approximate EM algorithm first proposed by Bernhardt et al. (2015). The GLMMs are specified using the same syntax as for package \code{glmmTMB} Brooks et @@ -16,7 +16,7 @@ types. } \author{ - \packageAuthor{GMVJM} + \packageAuthor{gmvjoint} } \references{ diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index adfca34..e29d313 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -13,7 +13,7 @@ Rcpp::Rostream& Rcpp::Rcerr = Rcpp::Rcpp_cerr_get(); // GP1_pmf_scalar double GP1_pmf_scalar(const double mu, const double phi, const double Y); -RcppExport SEXP _GMVJM_GP1_pmf_scalar(SEXP muSEXP, SEXP phiSEXP, SEXP YSEXP) { +RcppExport SEXP _gmvjoint_GP1_pmf_scalar(SEXP muSEXP, SEXP phiSEXP, SEXP YSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -26,7 +26,7 @@ END_RCPP } // ll_Gamma double ll_Gamma(const arma::vec& Y, const double& shape, const arma::vec& mu); -RcppExport SEXP _GMVJM_ll_Gamma(SEXP YSEXP, SEXP shapeSEXP, SEXP muSEXP) { +RcppExport SEXP _gmvjoint_ll_Gamma(SEXP YSEXP, SEXP shapeSEXP, SEXP muSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -39,7 +39,7 @@ END_RCPP } // ll_genpois double ll_genpois(const arma::vec& eta, const double phi, arma::vec& Y); -RcppExport SEXP _GMVJM_ll_genpois(SEXP etaSEXP, SEXP phiSEXP, SEXP YSEXP) { +RcppExport SEXP _gmvjoint_ll_genpois(SEXP etaSEXP, SEXP phiSEXP, SEXP YSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -52,7 +52,7 @@ END_RCPP } // logfti double logfti(const arma::vec& b, const arma::rowvec& S, const arma::mat& SS, const arma::rowvec& Fi, const arma::mat& Fu, const double l0i, const arma::rowvec& haz, const int Delta, const arma::vec& gamma_rep, const arma::vec& zeta); -RcppExport SEXP _GMVJM_logfti(SEXP bSEXP, SEXP SSEXP, SEXP SSSEXP, SEXP FiSEXP, SEXP FuSEXP, SEXP l0iSEXP, SEXP hazSEXP, SEXP DeltaSEXP, SEXP gamma_repSEXP, SEXP zetaSEXP) { +RcppExport SEXP _gmvjoint_logfti(SEXP bSEXP, SEXP SSEXP, SEXP SSSEXP, SEXP FiSEXP, SEXP FuSEXP, SEXP l0iSEXP, SEXP hazSEXP, SEXP DeltaSEXP, SEXP gamma_repSEXP, SEXP zetaSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -72,7 +72,7 @@ END_RCPP } // joint_density double joint_density(const arma::vec& b, const List Y, const List X, const List Z, const arma::vec& beta, const arma::mat& D, const List sigma, const List family, const int Delta, const arma::rowvec& S, const arma::rowvec& Fi, const double l0i, const arma::mat& SS, const arma::mat& Fu, const arma::rowvec& haz, const arma::vec& gamma_rep, const arma::vec& zeta, const List beta_inds, const List b_inds, const int K); -RcppExport SEXP _GMVJM_joint_density(SEXP bSEXP, SEXP YSEXP, SEXP XSEXP, SEXP ZSEXP, SEXP betaSEXP, SEXP DSEXP, SEXP sigmaSEXP, SEXP familySEXP, SEXP DeltaSEXP, SEXP SSEXP, SEXP FiSEXP, SEXP l0iSEXP, SEXP SSSEXP, SEXP FuSEXP, SEXP hazSEXP, SEXP gamma_repSEXP, SEXP zetaSEXP, SEXP beta_indsSEXP, SEXP b_indsSEXP, SEXP KSEXP) { +RcppExport SEXP _gmvjoint_joint_density(SEXP bSEXP, SEXP YSEXP, SEXP XSEXP, SEXP ZSEXP, SEXP betaSEXP, SEXP DSEXP, SEXP sigmaSEXP, SEXP familySEXP, SEXP DeltaSEXP, SEXP SSEXP, SEXP FiSEXP, SEXP l0iSEXP, SEXP SSSEXP, SEXP FuSEXP, SEXP hazSEXP, SEXP gamma_repSEXP, SEXP zetaSEXP, SEXP beta_indsSEXP, SEXP b_indsSEXP, SEXP KSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -102,7 +102,7 @@ END_RCPP } // maketau List maketau(const List& S, const List& Z); -RcppExport SEXP _GMVJM_maketau(SEXP SSEXP, SEXP ZSEXP) { +RcppExport SEXP _gmvjoint_maketau(SEXP SSEXP, SEXP ZSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -114,7 +114,7 @@ END_RCPP } // joint_density_ddb arma::vec joint_density_ddb(const arma::vec& b, const List Y, const List X, const List Z, const arma::vec& beta, const arma::mat& D, const List sigma, const List family, const int Delta, const arma::rowvec& S, const arma::rowvec& Fi, const double l0i, const arma::mat& SS, const arma::mat& Fu, const arma::rowvec& haz, const arma::vec& gamma_rep, const arma::vec& zeta, const List beta_inds, const List b_inds, const int K); -RcppExport SEXP _GMVJM_joint_density_ddb(SEXP bSEXP, SEXP YSEXP, SEXP XSEXP, SEXP ZSEXP, SEXP betaSEXP, SEXP DSEXP, SEXP sigmaSEXP, SEXP familySEXP, SEXP DeltaSEXP, SEXP SSEXP, SEXP FiSEXP, SEXP l0iSEXP, SEXP SSSEXP, SEXP FuSEXP, SEXP hazSEXP, SEXP gamma_repSEXP, SEXP zetaSEXP, SEXP beta_indsSEXP, SEXP b_indsSEXP, SEXP KSEXP) { +RcppExport SEXP _gmvjoint_joint_density_ddb(SEXP bSEXP, SEXP YSEXP, SEXP XSEXP, SEXP ZSEXP, SEXP betaSEXP, SEXP DSEXP, SEXP sigmaSEXP, SEXP familySEXP, SEXP DeltaSEXP, SEXP SSEXP, SEXP FiSEXP, SEXP l0iSEXP, SEXP SSSEXP, SEXP FuSEXP, SEXP hazSEXP, SEXP gamma_repSEXP, SEXP zetaSEXP, SEXP beta_indsSEXP, SEXP b_indsSEXP, SEXP KSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -144,7 +144,7 @@ END_RCPP } // Sbeta arma::vec Sbeta(const arma::vec& beta, const List& X, const List& Y, const List& Z, const List& b, const List& sigma, const List& family, const List& beta_inds, const int K, const bool quad, const List& tau, const arma::vec& w, const arma::vec& v); -RcppExport SEXP _GMVJM_Sbeta(SEXP betaSEXP, SEXP XSEXP, SEXP YSEXP, SEXP ZSEXP, SEXP bSEXP, SEXP sigmaSEXP, SEXP familySEXP, SEXP beta_indsSEXP, SEXP KSEXP, SEXP quadSEXP, SEXP tauSEXP, SEXP wSEXP, SEXP vSEXP) { +RcppExport SEXP _gmvjoint_Sbeta(SEXP betaSEXP, SEXP XSEXP, SEXP YSEXP, SEXP ZSEXP, SEXP bSEXP, SEXP sigmaSEXP, SEXP familySEXP, SEXP beta_indsSEXP, SEXP KSEXP, SEXP quadSEXP, SEXP tauSEXP, SEXP wSEXP, SEXP vSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -167,7 +167,7 @@ END_RCPP } // Hbeta arma::mat Hbeta(const arma::vec& beta, const List& X, const List& Y, const List& Z, const List& b, const List& sigma, const List& family, const List& beta_inds, const int K, const bool& quad, const List& tau, const arma::vec& w, const arma::vec& v); -RcppExport SEXP _GMVJM_Hbeta(SEXP betaSEXP, SEXP XSEXP, SEXP YSEXP, SEXP ZSEXP, SEXP bSEXP, SEXP sigmaSEXP, SEXP familySEXP, SEXP beta_indsSEXP, SEXP KSEXP, SEXP quadSEXP, SEXP tauSEXP, SEXP wSEXP, SEXP vSEXP) { +RcppExport SEXP _gmvjoint_Hbeta(SEXP betaSEXP, SEXP XSEXP, SEXP YSEXP, SEXP ZSEXP, SEXP bSEXP, SEXP sigmaSEXP, SEXP familySEXP, SEXP beta_indsSEXP, SEXP KSEXP, SEXP quadSEXP, SEXP tauSEXP, SEXP wSEXP, SEXP vSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -190,7 +190,7 @@ END_RCPP } // vare_update double vare_update(const arma::mat& X, const arma::vec& Y, const arma::mat& Z, const arma::vec& b, const arma::vec& beta, const arma::vec& tau, const arma::vec& w, const arma::vec& v); -RcppExport SEXP _GMVJM_vare_update(SEXP XSEXP, SEXP YSEXP, SEXP ZSEXP, SEXP bSEXP, SEXP betaSEXP, SEXP tauSEXP, SEXP wSEXP, SEXP vSEXP) { +RcppExport SEXP _gmvjoint_vare_update(SEXP XSEXP, SEXP YSEXP, SEXP ZSEXP, SEXP bSEXP, SEXP betaSEXP, SEXP tauSEXP, SEXP wSEXP, SEXP vSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -208,7 +208,7 @@ END_RCPP } // phi_update List phi_update(const arma::vec& b, const arma::mat& X, const arma::vec& Y, const arma::mat& Z, const arma::vec& beta, const double phi, const arma::vec& w, const arma::vec& v, const arma::vec& tau); -RcppExport SEXP _GMVJM_phi_update(SEXP bSEXP, SEXP XSEXP, SEXP YSEXP, SEXP ZSEXP, SEXP betaSEXP, SEXP phiSEXP, SEXP wSEXP, SEXP vSEXP, SEXP tauSEXP) { +RcppExport SEXP _gmvjoint_phi_update(SEXP bSEXP, SEXP XSEXP, SEXP YSEXP, SEXP ZSEXP, SEXP betaSEXP, SEXP phiSEXP, SEXP wSEXP, SEXP vSEXP, SEXP tauSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -227,7 +227,7 @@ END_RCPP } // Sgammazeta arma::vec Sgammazeta(arma::vec& gammazeta, arma::vec& b, List Sigma, arma::rowvec& S, arma::mat& SS, arma::mat& Fu, arma::rowvec& Fi, arma::vec& haz, int Delta, arma::vec& w, arma::vec& v, List b_inds, int K, int q, long double eps); -RcppExport SEXP _GMVJM_Sgammazeta(SEXP gammazetaSEXP, SEXP bSEXP, SEXP SigmaSEXP, SEXP SSEXP, SEXP SSSEXP, SEXP FuSEXP, SEXP FiSEXP, SEXP hazSEXP, SEXP DeltaSEXP, SEXP wSEXP, SEXP vSEXP, SEXP b_indsSEXP, SEXP KSEXP, SEXP qSEXP, SEXP epsSEXP) { +RcppExport SEXP _gmvjoint_Sgammazeta(SEXP gammazetaSEXP, SEXP bSEXP, SEXP SigmaSEXP, SEXP SSEXP, SEXP SSSEXP, SEXP FuSEXP, SEXP FiSEXP, SEXP hazSEXP, SEXP DeltaSEXP, SEXP wSEXP, SEXP vSEXP, SEXP b_indsSEXP, SEXP KSEXP, SEXP qSEXP, SEXP epsSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -252,7 +252,7 @@ END_RCPP } // Hgammazeta arma::mat Hgammazeta(arma::vec& gammazeta, arma::vec& b, List Sigma, arma::rowvec& S, arma::mat& SS, arma::mat& Fu, arma::rowvec& Fi, arma::vec& haz, int Delta, arma::vec& w, arma::vec& v, List b_inds, int K, int q, double eps); -RcppExport SEXP _GMVJM_Hgammazeta(SEXP gammazetaSEXP, SEXP bSEXP, SEXP SigmaSEXP, SEXP SSEXP, SEXP SSSEXP, SEXP FuSEXP, SEXP FiSEXP, SEXP hazSEXP, SEXP DeltaSEXP, SEXP wSEXP, SEXP vSEXP, SEXP b_indsSEXP, SEXP KSEXP, SEXP qSEXP, SEXP epsSEXP) { +RcppExport SEXP _gmvjoint_Hgammazeta(SEXP gammazetaSEXP, SEXP bSEXP, SEXP SigmaSEXP, SEXP SSEXP, SEXP SSSEXP, SEXP FuSEXP, SEXP FiSEXP, SEXP hazSEXP, SEXP DeltaSEXP, SEXP wSEXP, SEXP vSEXP, SEXP b_indsSEXP, SEXP KSEXP, SEXP qSEXP, SEXP epsSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -277,7 +277,7 @@ END_RCPP } // lambdaUpdate arma::mat lambdaUpdate(List survtimes, arma::mat& ft, arma::vec& gamma, arma::vec& gamma_rep, arma::vec& zeta, List S, List Sigma, List b, arma::vec& w, arma::vec& v, List b_inds, int K, int q); -RcppExport SEXP _GMVJM_lambdaUpdate(SEXP survtimesSEXP, SEXP ftSEXP, SEXP gammaSEXP, SEXP gamma_repSEXP, SEXP zetaSEXP, SEXP SSEXP, SEXP SigmaSEXP, SEXP bSEXP, SEXP wSEXP, SEXP vSEXP, SEXP b_indsSEXP, SEXP KSEXP, SEXP qSEXP) { +RcppExport SEXP _gmvjoint_lambdaUpdate(SEXP survtimesSEXP, SEXP ftSEXP, SEXP gammaSEXP, SEXP gamma_repSEXP, SEXP zetaSEXP, SEXP SSEXP, SEXP SigmaSEXP, SEXP bSEXP, SEXP wSEXP, SEXP vSEXP, SEXP b_indsSEXP, SEXP KSEXP, SEXP qSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -300,7 +300,7 @@ END_RCPP } // joint_density_sdb arma::mat joint_density_sdb(const arma::vec& b, const List Y, const List X, const List Z, const arma::vec& beta, const arma::mat& D, const List sigma, const List family, const int Delta, const arma::rowvec& S, const arma::rowvec& Fi, const double l0i, const arma::mat& SS, const arma::mat& Fu, const arma::rowvec& haz, const arma::vec& gamma_rep, const arma::vec& zeta, const List beta_inds, const List b_inds, const int K, double eps); -RcppExport SEXP _GMVJM_joint_density_sdb(SEXP bSEXP, SEXP YSEXP, SEXP XSEXP, SEXP ZSEXP, SEXP betaSEXP, SEXP DSEXP, SEXP sigmaSEXP, SEXP familySEXP, SEXP DeltaSEXP, SEXP SSEXP, SEXP FiSEXP, SEXP l0iSEXP, SEXP SSSEXP, SEXP FuSEXP, SEXP hazSEXP, SEXP gamma_repSEXP, SEXP zetaSEXP, SEXP beta_indsSEXP, SEXP b_indsSEXP, SEXP KSEXP, SEXP epsSEXP) { +RcppExport SEXP _gmvjoint_joint_density_sdb(SEXP bSEXP, SEXP YSEXP, SEXP XSEXP, SEXP ZSEXP, SEXP betaSEXP, SEXP DSEXP, SEXP sigmaSEXP, SEXP familySEXP, SEXP DeltaSEXP, SEXP SSEXP, SEXP FiSEXP, SEXP l0iSEXP, SEXP SSSEXP, SEXP FuSEXP, SEXP hazSEXP, SEXP gamma_repSEXP, SEXP zetaSEXP, SEXP beta_indsSEXP, SEXP b_indsSEXP, SEXP KSEXP, SEXP epsSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; diff --git a/src/funs.cpp b/src/funs.cpp index 5f9c1be..7095923 100644 --- a/src/funs.cpp +++ b/src/funs.cpp @@ -624,7 +624,7 @@ arma::mat lambdaUpdate(List survtimes, arma::mat& ft, arma::vec& gamma, arma::ve } double mu = as_scalar(exp(S_i * zeta + Fst * rhs)); for(int l = 0; l < gh; l++){ - store(j, i) += as_scalar(w[l] * mu * exp(v[l] * 0.5 * tau)); + store(j, i) += as_scalar(w[l] * mu * exp(v[l] * sqrt(tau))); } } } diff --git a/src/init.c b/src/init.c index 312aa06..d44a16f 100644 --- a/src/init.c +++ b/src/init.c @@ -8,42 +8,42 @@ */ /* .Call calls */ -extern SEXP _GMVJM_GP1_pmf_scalar(SEXP, SEXP, SEXP); -extern SEXP _GMVJM_Hbeta(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP); -extern SEXP _GMVJM_Hgammazeta(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP); -extern SEXP _GMVJM_joint_density(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP); -extern SEXP _GMVJM_joint_density_ddb(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP); -extern SEXP _GMVJM_joint_density_sdb(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP); -extern SEXP _GMVJM_lambdaUpdate(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP); -extern SEXP _GMVJM_ll_Gamma(SEXP, SEXP, SEXP); -extern SEXP _GMVJM_ll_genpois(SEXP, SEXP, SEXP); -extern SEXP _GMVJM_logfti(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP); -extern SEXP _GMVJM_maketau(SEXP, SEXP); -extern SEXP _GMVJM_phi_update(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP); -extern SEXP _GMVJM_Sbeta(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP); -extern SEXP _GMVJM_Sgammazeta(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP); -extern SEXP _GMVJM_vare_update(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP); +extern SEXP _gmvjoint_GP1_pmf_scalar(SEXP, SEXP, SEXP); +extern SEXP _gmvjoint_Hbeta(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP); +extern SEXP _gmvjoint_Hgammazeta(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP); +extern SEXP _gmvjoint_joint_density(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP); +extern SEXP _gmvjoint_joint_density_ddb(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP); +extern SEXP _gmvjoint_joint_density_sdb(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP); +extern SEXP _gmvjoint_lambdaUpdate(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP); +extern SEXP _gmvjoint_ll_Gamma(SEXP, SEXP, SEXP); +extern SEXP _gmvjoint_ll_genpois(SEXP, SEXP, SEXP); +extern SEXP _gmvjoint_logfti(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP); +extern SEXP _gmvjoint_maketau(SEXP, SEXP); +extern SEXP _gmvjoint_phi_update(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP); +extern SEXP _gmvjoint_Sbeta(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP); +extern SEXP _gmvjoint_Sgammazeta(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP); +extern SEXP _gmvjoint_vare_update(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP); static const R_CallMethodDef CallEntries[] = { - {"_GMVJM_GP1_pmf_scalar", (DL_FUNC) &_GMVJM_GP1_pmf_scalar, 3}, - {"_GMVJM_Hbeta", (DL_FUNC) &_GMVJM_Hbeta, 13}, - {"_GMVJM_Hgammazeta", (DL_FUNC) &_GMVJM_Hgammazeta, 15}, - {"_GMVJM_joint_density", (DL_FUNC) &_GMVJM_joint_density, 20}, - {"_GMVJM_joint_density_ddb", (DL_FUNC) &_GMVJM_joint_density_ddb, 20}, - {"_GMVJM_joint_density_sdb", (DL_FUNC) &_GMVJM_joint_density_sdb, 21}, - {"_GMVJM_lambdaUpdate", (DL_FUNC) &_GMVJM_lambdaUpdate, 13}, - {"_GMVJM_ll_Gamma", (DL_FUNC) &_GMVJM_ll_Gamma, 3}, - {"_GMVJM_ll_genpois", (DL_FUNC) &_GMVJM_ll_genpois, 3}, - {"_GMVJM_logfti", (DL_FUNC) &_GMVJM_logfti, 10}, - {"_GMVJM_maketau", (DL_FUNC) &_GMVJM_maketau, 2}, - {"_GMVJM_phi_update", (DL_FUNC) &_GMVJM_phi_update, 9}, - {"_GMVJM_Sbeta", (DL_FUNC) &_GMVJM_Sbeta, 13}, - {"_GMVJM_Sgammazeta", (DL_FUNC) &_GMVJM_Sgammazeta, 15}, - {"_GMVJM_vare_update", (DL_FUNC) &_GMVJM_vare_update, 8}, + {"_gmvjoint_GP1_pmf_scalar", (DL_FUNC) &_gmvjoint_GP1_pmf_scalar, 3}, + {"_gmvjoint_Hbeta", (DL_FUNC) &_gmvjoint_Hbeta, 13}, + {"_gmvjoint_Hgammazeta", (DL_FUNC) &_gmvjoint_Hgammazeta, 15}, + {"_gmvjoint_joint_density", (DL_FUNC) &_gmvjoint_joint_density, 20}, + {"_gmvjoint_joint_density_ddb", (DL_FUNC) &_gmvjoint_joint_density_ddb, 20}, + {"_gmvjoint_joint_density_sdb", (DL_FUNC) &_gmvjoint_joint_density_sdb, 21}, + {"_gmvjoint_lambdaUpdate", (DL_FUNC) &_gmvjoint_lambdaUpdate, 13}, + {"_gmvjoint_ll_Gamma", (DL_FUNC) &_gmvjoint_ll_Gamma, 3}, + {"_gmvjoint_ll_genpois", (DL_FUNC) &_gmvjoint_ll_genpois, 3}, + {"_gmvjoint_logfti", (DL_FUNC) &_gmvjoint_logfti, 10}, + {"_gmvjoint_maketau", (DL_FUNC) &_gmvjoint_maketau, 2}, + {"_gmvjoint_phi_update", (DL_FUNC) &_gmvjoint_phi_update, 9}, + {"_gmvjoint_Sbeta", (DL_FUNC) &_gmvjoint_Sbeta, 13}, + {"_gmvjoint_Sgammazeta", (DL_FUNC) &_gmvjoint_Sgammazeta, 15}, + {"_gmvjoint_vare_update", (DL_FUNC) &_gmvjoint_vare_update, 8}, {NULL, NULL, 0} }; -void R_init_GMVJM(DllInfo *dll) +void R_init_gmvjoint(DllInfo *dll) { R_registerRoutines(dll, NULL, CallEntries, NULL, NULL); R_useDynamicSymbols(dll, FALSE);