194 saem generalized log likelihood

DESCRIPTION

@@ -89,7 +89,7 @@ LinkingTo:
     dparser (>= 0.1.8),
     lbfgsb3c,
     Rcpp,
-    RcppArmadillo (>= 0.5.600.2.0),
+    RcppArmadillo (>= 0.11.2.3.1),
     RcppEigen (>= 0.3.3.3.0),
     rxode2 (>= 2.0.7),
     StanHeaders (>= 2.18.0)

NAMESPACE

@@ -57,7 +57,6 @@ S3method(nlmixr2Est,simulate)
 S3method(nmGetDistributionSaemLines,default)
 S3method(nmGetDistributionSaemLines,norm)
 S3method(nmGetDistributionSaemLines,rxUi)
-S3method(nmGetDistributionSaemLines,t)
 S3method(nmObjGet,atol)
 S3method(nmObjGet,coefficients)
 S3method(nmObjGet,condition)
@@ -205,6 +204,11 @@ S3method(rxUiGet,saemAres)
 S3method(rxUiGet,saemBres)
 S3method(rxUiGet,saemCovars)
 S3method(rxUiGet,saemCres)
+S3method(rxUiGet,saemDistribution)
+S3method(rxUiGet,saemErrDf)
+S3method(rxUiGet,saemErrMu)
+S3method(rxUiGet,saemErrMuFix)
+S3method(rxUiGet,saemErrMuNames)
 S3method(rxUiGet,saemEtaNames)
 S3method(rxUiGet,saemEtaTrans)
 S3method(rxUiGet,saemFixed)
@@ -225,6 +229,7 @@ S3method(rxUiGet,saemModResTotalResiduals)
 S3method(rxUiGet,saemModel)
 S3method(rxUiGet,saemModel0)
 S3method(rxUiGet,saemModelList)
+S3method(rxUiGet,saemModelNeedsLlik)
 S3method(rxUiGet,saemModelOmega)
 S3method(rxUiGet,saemModelOmegaFixed)
 S3method(rxUiGet,saemModelOmegaFixedValues)
@@ -244,6 +249,7 @@ S3method(rxUiGet,saemParamsToEstimate)
 S3method(rxUiGet,saemParamsToEstimateCov)
 S3method(rxUiGet,saemPropT)
 S3method(rxUiGet,saemResFixed)
+S3method(rxUiGet,saemResLlMod)
 S3method(rxUiGet,saemResMod)
 S3method(rxUiGet,saemResNames)
 S3method(rxUiGet,saemResValue)

R/RcppExports.R

@@ -168,6 +168,10 @@ augPredTrans <- function(pred, ipred, lambda, yjIn, low, hi) {
     .Call(`_nlmixr2est_augPredTrans`, pred, ipred, lambda, yjIn, low, hi)
 }
 
+saem_user_opt_ll_fun <- function(resParsNV) {
+    .Call(`_nlmixr2est_saem_user_opt_ll_fun`, resParsNV)
+}
+
 saem_do_pred <- function(in_phi, in_evt, in_opt) {
     .Call(`_nlmixr2est_saem_do_pred`, in_phi, in_evt, in_opt)
 }

R/focei.R

@@ -1760,7 +1760,6 @@ nlmixr2Est.fo <- function(env, ...) {
 #'@export
 nlmixr2Est.output <- function(env, ...) {
   .ui <- env$ui
-  rxode2::assertRxUiTransformNormal(.ui, " for the estimation routine 'output'", .var.name=.ui$modelName)
   rxode2::assertRxUiRandomOnIdOnly(.ui, " for the estimation routine 'output'", .var.name=.ui$modelName)
 
   .foceiFamilyControl(env, ...)

R/nmObjGet.R

@@ -131,6 +131,11 @@ nmObjGet.phiR <- function(x, ...) {
   if (is.null(.phi)) return(NULL)
   .ret <- lapply(seq_along(.phi), function(i) {
     .cov <- .phi[[i]]
+    .d <- diag(.cov)
+    if (any(.d == 0) || any(!is.finite(.d))) {
+      .d1 <- length(.d)
+      return(matrix(rep(NA, .d1*.d1), .d1, .d1))
+    }
     .sd2 <- sqrt(diag(.cov))
     .cor <- stats::cov2cor(.cov)
     dimnames(.cor) <- dimnames(.cov)
@@ -150,11 +155,13 @@ nmObjGet.phiSE <- function(x, ...) {
   if (is.null(.phi)) return(NULL)
   .ret <- as.data.frame(t(vapply(seq_along(.phi), function(i) {
     .cov <- .phi[[i]]
-    sqrt(diag(.cov))
+    suppressWarnings(sqrt(diag(.cov)))
   }, double(dim(.phi[[1]])[1]))))
   names(.ret) <- paste0("se(", names(.ret), ")")
   .id <- seq_along(.phi)
-  if (!is.null(names(.phi))) .id <- factor(.id, levels=names(.phi))
+  if (!is.null(names(.phi))) {
+    .id <- names(.phi)
+  }
   cbind(data.frame(ID=.id), .ret)
 }
 attr(nmObjGet.phiSE, "desc") <- "standard error of each individual's eta (if present)"
@@ -168,11 +175,13 @@ nmObjGet.phiRSE <- function(x, ...) {
   if (is.null(.phi)) return(NULL)
   .ret <- as.data.frame(t(vapply(seq_along(.phi), function(i) {
     .cov <- .phi[[i]]
-    sqrt(diag(.cov))/unlist(.eta[i,])*100
+    suppressWarnings(sqrt(diag(.cov))/unlist(.eta[i,])*100)
   }, double(dim(.phi[[1]])[1]))))
   names(.ret) <- paste0("rse(", names(.ret), ")%")
   .id <- seq_along(.phi)
-  if (!is.null(names(.phi))) .id <- factor(.id, levels=names(.phi))
+  if (!is.null(names(.phi))) {
+    .id <- names(.phi)
+  }
   cbind(data.frame(ID=.id), .ret)
 }
 attr(nmObjGet.phiRSE, "desc") <- "relative standard error of each individual's eta (if present)"

R/saem.R

@@ -9,6 +9,22 @@
   .opt$estimate
 }
 
+.saemLlOpt1 <- function(p1) {
+  .opt <- stats::nlm(saem_user_opt_ll_fun, p1)
+  .opt$estimate
+}
+
+.saemLlOptNewUoa <- function(par, rhoend, maxfun) {
+  .ctl <- list(npt = length(par) * 2 + 1,
+               rhobeg = 0.2,
+               rhoend = rhoend,
+               maxfun = maxfun,
+               iprint = 0L)
+  .ret <- minqa::newuoa(par, saem_user_opt_ll_fun,
+                        control = .ctl)
+  .ret$par
+}
+
 .newuoa <- function(par, fn, gr, lower = -Inf, upper = Inf, control = list(), ...) {
   .ctl <- control
   if (is.null(.ctl$npt)) .ctl$npt <- length(par) * 2 + 1
@@ -202,7 +218,6 @@
                         inits=.inits,
                         mcmc=rxode2::rxGetControl(ui, "mcmc", list(niter = c(200, 300), nmc = 3, nu = c(2, 2, 2))),
                         rxControl=rxode2::rxGetControl(ui, "rxControl", rxode2::rxControl()),
-                        distribution="normal",
                         fixedOmega=ui$saemModelOmegaFixed,
                         fixedOmegaValues=ui$saemModelOmegaFixedValues,
                         parHistThetaKeep=ui$saemParHistThetaKeep,
@@ -221,7 +236,9 @@
                         perNoCor=rxode2::rxGetControl(ui, "perNoCor", 0.75),
                         perFixOmega=rxode2::rxGetControl(ui, "perFixOmega", 0.1),
                         perFixResid=rxode2::rxGetControl(ui, "perFixResid", 0.1),
-                        resFixed=ui$saemResFixed)
+                        resFixed=ui$saemResFixed,
+                        distribution=ui$saemDistribution,
+                        resLlMod=ui$saemResLlMod)
     .print <- rxode2::rxGetControl(ui, "print", 1)
     if (inherits(.print, "numeric")) {
       .cfg$print <- as.integer(.print)
@@ -254,7 +271,7 @@
   if (is.null(.control)) {
     .control <- saemControl()
   }
-  if (!inherits(.control, "saemControl")){
+  if (!inherits(.control, "saemControl")) {
     .control <- do.call(nlmixr2est::saemControl, .control)
   }
   assign("control", .control, envir=.ui)
@@ -757,7 +774,6 @@ nmObjGetFoceiControl.saem <- function(x, ...) {
 #' @export
 nlmixr2Est.saem <- function(env, ...) {
   .ui <- env$ui
-  rxode2::assertRxUiTransformNormal(.ui, " for the estimation routine 'saem'", .var.name=.ui$modelName)
   rxode2::assertRxUiRandomOnIdOnly(.ui, " for the estimation routine 'saem'", .var.name=.ui$modelName)
   rxode2::assertRxUiEstimatedResiduals(.ui, " for the estimation routine 'saem'", .var.name=.ui$modelName)
   rxode2::assertRxUiMixedOnly(.ui, " for the estimation routine 'saem'", .var.name=.ui$modelName)

R/saemRxUiGet.R

@@ -159,7 +159,8 @@ rxUiGet.saemFixed <- function(x, ...) {
 #' @export
 rxUiGet.saemEtaTrans <- function(x, ...) {
   .ui <- x[[1]]
-  .etas <- .ui$iniDf[!is.na(.ui$iniDf$neta1), ]
+  .iniDf <- .ui$iniDf
+  .etas <- .iniDf[!is.na(.iniDf$neta1), ]
   .etas <- .etas$name[.etas$neta1 == .etas$neta2]
   .thetas <- rxUiGet.saemParamsToEstimateCov(x, ...)
   .muRefDataFrame <- .ui$muRefDataFrame
@@ -176,6 +177,7 @@ rxUiGet.saemEtaTrans <- function(x, ...) {
   }, integer(1), USE.NAMES=FALSE)
 }
 #attr(rxUiGet.saemEtaTrans, "desc") <- "Get the saem eta to theta translation"
+
 #' @export
 rxUiGet.saemOmegaTrans <- function(x, ...) {
   .etaTrans <- rxUiGet.saemEtaTrans(x, ...)
@@ -436,7 +438,8 @@ rxUiGet.saemResValue <- function(x, ...) {
 #' @export
 rxUiGet.saemEtaNames <- function(x, ...) {
   .ui <- x[[1]]
-  .etaNames <- .ui$iniDf[!is.na(.ui$iniDf$neta1), ]
+  .iniDf <- .ui$iniDf
+  .etaNames <- .iniDf[!is.na(.iniDf$neta1), ]
   .etaNames <- .etaNames[.etaNames$neta1 == .etaNames$neta2, "name"]
   ## .etaTrans <- rxUiGet.saemOmegaTrans(x, ...)
   .etaTrans <- rxUiGet.saemEtaTrans(x, ...)
@@ -453,7 +456,8 @@ rxUiGet.saemEtaNames <- function(x, ...) {
 #' @export
 rxUiGet.saemParHistOmegaKeep <- function(x, ...) {
   .ui <- x[[1]]
-  .etaNames <- .ui$iniDf[!is.na(.ui$iniDf$neta1), ]
+  .iniDf <-.ui$iniDf
+  .etaNames <- .iniDf[!is.na(.iniDf$neta1), ]
   .etaNames <- .etaNames[.etaNames$neta1 == .etaNames$neta2,]
   .names <- rxUiGet.saemEtaNames(x, ...)
   vapply(.names, function(etaName){
@@ -482,7 +486,11 @@ rxUiGet.saemParHistNames <- function(x, ...) {
   #join_cols(join_cols(Plambda, Gamma2_phi1.diag()), vcsig2).t();
   .plambda <- rxUiGet.saemParamsToEstimate(x, ...)
   .plambda <- .plambda[!rxUiGet.saemFixed(x, ...)]
-  c(.plambda, rxUiGet.saemParHistEtaNames(x, ...), rxUiGet.saemParHistResNames(x, ...))
+  if (rxUiGet.saemModelNeedsLlik(x, ...)) {
+    c(.plambda, rxUiGet.saemParHistEtaNames(x, ...))
+  } else {
+    c(.plambda, rxUiGet.saemParHistEtaNames(x, ...), rxUiGet.saemParHistResNames(x, ...))    
+  }
 }
 
 #' @export
@@ -585,7 +593,7 @@ rxUiGet.saemLogEta <- function(x, ...) {
     .w <- which(.ce$parameter == x)
     if (length(.w) == 1L) return(.ce$curEval[.w] == "exp")
     FALSE
-  }, logical(1))
+  }, logical(1),USE.NAMES=TRUE)
 }
 #attr(rxUiGet.saemLogEta, "desc") <- "saem's log.eta for saem"
 
@@ -625,7 +633,7 @@ rxUiGet.saemInitTheta <- function(x, ...) {
   .n <- vapply(.theta, function(x) ifelse(x, "FIXED", ""),
                character(1), USE.NAMES=FALSE)
   .ret <- vapply(seq_along(.logEta),
-                  function(i){
+                  function(i) {
                     .isEta <- any(.names[i] %in% .etaNames)
                     if (.logEta[i]) {
                       if (.isEta) {
@@ -695,7 +703,8 @@ rxUiGet.saemInit <- function(x, ...) {
 rxUiGet.saemThetaDataFrame <- function(x, ...) {
   .ui <- x[[1]]
   .theta <- .ui$theta
-  .fixed <- .ui$iniDf[!is.na(.ui$iniDf$ntheta), "fix"]
+  .iniDf <-.ui$iniDf
+  .fixed <- .iniDf[!is.na(.iniDf$ntheta), "fix"]
   data.frame(lower= -Inf, theta=.theta, fixed=.fixed, upper=Inf, row.names=names(.theta))
 }
 #attr(rxUiGet.saemThetaDataFrame, "desc") <- "Get theta data frame"

R/saemRxUiGetLlik.R

@@ -0,0 +1,103 @@
+#' @export
+rxUiGet.saemErrDf <- function(x, ...) {
+  if (!rxUiGet.saemModelNeedsLlik(x, ...)) return(NULL)
+  .x <- x[[1]]
+  .w <- !is.na(.x$iniDf$ntheta) & !is.na(.x$iniDf$err)
+  if (length(.w) == 0L) return(NULL)
+  .x$iniDf[.w, ]
+}
+
+#' @export
+rxUiGet.saemErrMuNames <- function(x, ...) {
+  .iniErr <- rxUiGet.saemErrDf(x, ...)
+  if (is.null(.iniErr)) return(NULL)
+  paste0("rx_err_", seq_along(.iniErr$name), "_")
+}
+
+#' @export
+rxUiGet.saemResLlMod <- function(x, ...) {
+  .iniErr <- rxUiGet.saemErrDf(x, ...)
+  if (is.null(.iniErr)) return(numeric(0))
+  setNames(
+    vapply(seq_along(.iniErr$est),
+           function(i) {
+             .low <- .iniErr$lower[i]
+             .hi <- .iniErr$upper[i]
+             .est <- .iniErr$est[i]
+             if (!is.finite(.low) && !is.finite(.hi)) {
+               return(logit(.est, .low, .hi))
+             } else if (is.finite(.low) && !is.finite(.hi)) {
+               return(log(.est - .low))
+             } else if (!is.finite(.low) && is.finite(.hi)) {
+               return(log(.low - .est))
+             } else if (!is.finite(.low) && !is.finite(.hi)) {
+               return(.est)
+             }
+           }, numeric(1), USE.NAMES=FALSE),
+    rxUiGet.saemErrMuNames(x, ...))
+}
+
+#' @export
+rxUiGet.saemErrMuFix <- function(x, ...) {
+  .iniErr <- rxUiGet.saemErrDf(x, ...)
+  if (is.null(.iniErr)) return(NULL)
+  setNames(
+    vapply(seq_along(.iniErr$est),
+           function(i) {
+             .iniErr$fix[i]
+           }, logical(1), USE.NAMES=FALSE),
+    rxUiGet.saemErrMuNames(x, ...)
+    )
+}
+
+
+#' @export
+rxUiGet.saemErrMu <- function(x, ...) {
+  .iniErr <- rxUiGet.saemErrDf(x, ...)
+  if (is.null(.iniErr)) return(NULL)
+  lapply(seq_along(.iniErr$est),
+         function(i) {
+           .name <- .iniErr$name[i]
+           .low <- .iniErr$lower[i]
+           .hi <- .iniErr$upper[i]
+           if (!is.finite(.low) && !is.finite(.hi)) {
+             return(str2lang(paste0(.name, "<- expit(rx_err_", i,
+                                    "_, ", .low, ", ", .hi, ")")))
+           } else if (is.finite(.low) && !is.finite(.hi)) {
+             return(str2lang(paste0(.name, "<- ", .low, " + exp(rx_err_", i, "_)")))
+           } else if (!is.finite(.low) && is.finite(.hi)) {
+             return(str2lang(paste0(.name, "<- ", .hi, " - exp(rx_err_", i, "_)")))
+           } else if (!is.finite(.low) && !is.finite(.hi)) {
+             return(str2lang(paste0(.name, "<- rx_err_", i, "_)")))
+           }
+         })
+}
+
+#' @export
+rxUiGet.saemModelNeedsLlik <- function(x, ...) {
+  .ui <- x[[1]]
+  .ret <- any(.ui$predDf$distribution != "norm")
+  if (.ret) return(.ret)
+  # FIXME: check for complicated error expressions
+  FALSE
+}
+
+.saemIntegrateErrMu <- function(ui, lines) {
+  if (ui$saemModelNeedsLlik) {
+    .ref <- ui$saemErrMu
+    .ret <- lines
+    .ret[[2]] <- as.call(c(list(quote(`{`)),
+                           .ref,
+                           lapply(seq_along(lines[[2]])[-1], function(i) {
+                             lines[[2]][[i]]
+                           })))
+    return(.ret)
+  }
+  lines
+}
+
+#' @export
+rxUiGet.saemDistribution <- function(x, ...) {
+  if (rxUiGet.saemModelNeedsLlik(x, ...)) return(2L)
+  1L
+}