Nothing
R/saem_fit.R
In nlmixr: Nonlinear Mixed Effects Models in Population PK/PD
Defines functions
genSaemUserFunction
## saem_fit.R: population PK/PD modeling library
##
## Copyright (C) 2014 - 2016 Wenping Wang
##
## This file is part of nlmixr.
##
## nlmixr is free software: you can redistribute it and/or modify it
## under the terms of the GNU General Public License as published by
## the Free Software Foundation, either version 2 of the License, or
## (at your option) any later version.
##
## nlmixr is distributed in the hope that it will be useful, but
## WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with nlmixr. If not, see <http://www.gnu.org/licenses/>.
# genSaemUserFunction(f$rxode.pred, f$saem.pars, f$pred, f$error)
genSaemUserFunction <- function(model, PKpars = attr(model, "default.pars"), pred = NULL, err = NULL,
control = saemControl(), inPars = NULL) {
.x <- deparse(body(pred))
.len <- length(.x)
.x <- if (.x[1] == "{") .x[2:(.len - 1)] else .x
.len <- length(.x)
.nendpnt <- .len
.mod <- RxODE::RxODE(RxODE::rxGenSaem(model, function() {
return(nlmixr_pred)
}, PKpars,
sum.prod = control$sum.prod,
optExpression = control$optExpression,
loadSymengine=control$loadSymengine
))
.fnPred <- bquote(function(a, b, c) {
RxODE::rxLoad(.(.mod))
RxODE::rxLock(.(.mod))
RxODE::rxAllowUnload(FALSE)
on.exit({
RxODE::rxUnlock(.(.mod))
RxODE::rxAllowUnload(TRUE)
RxODE::rxSolveFree()
})
.Call(`_nlmixr_saem_do_pred`, a, b, c)
})
.fn <- bquote(function(a, b, c) {
RxODE::rxLoad(.(.mod))
RxODE::rxLock(.(.mod))
on.exit({
RxODE::rxUnlock(.(.mod))
RxODE::rxAllowUnload(TRUE)
RxODE::rxSolveFree()
})
if (missing(b) && missing(c)) {
.ret <- .Call(`_nlmixr_saem_fit`, a, PACKAGE = "nlmixr")
attr(.ret, "dopred") <- .(.fnPred)
return(.ret)
} else {
.curFn <- .(.fnPred)
return(.curFn(a, b, c))
}
})
.param <- RxODE::rxParam(.mod)
.inits <- names(RxODE::rxInits(.mod))
.nrhs <- length(.param) - length(.inits)
if (any(.param == "CMT")) {
inPars <- unique(c(inPars, "CMT"))
}
.parmUpdate <- sapply(.param, function(x) {
if (any(x == inPars)) {
return(0L)
} else if (any(x == .inits)) {
return(0L)
} else {
return(1L)
}
})
.fn <- eval(.fn)
attr(.fn, "form") <- "ode" ## Not sure this is necessary any more
attr(.fn, "neq") <- length(RxODE::rxState(.mod))
attr(.fn, "nlhs") <- length(RxODE::rxLhs(.mod))
attr(.fn, "nrhs") <- sum(.parmUpdate)
attr(.fn, "paramUpdate") <- .parmUpdate
attr(.fn, "rx") <- .mod
attr(.fn, "inPars") <- inPars
attr(.fn, "nendpnt") <- .nendpnt # not calculated; Is this a problem?
return(.fn)
}
#' Generate an SAEM model
#'
#' Generate an SAEM model using either closed-form solutions or ODE-based model definitions
#'
#' @param model a compiled SAEM model by gen_saem_user_fn()
#' @param PKpars PKpars function
#' @param pred pred function; This will be a focei-style pred
#' @param inPars a character vector of parameters to be read from the
#' input dataset (including time varying covariates)
#' @details Fit a generalized nonlinear mixed-effect model using the
#' Stochastic Approximation Expectation-Maximization (SAEM)
#' algorithm
#'
#' @return A user function based on the model to run the SAEM code
#'
#' @author Matthew Fidler & Wenping Wang
#' @keywords internal
#' @export
gen_saem_user_fn <- genSaemUserFunction
#' Configure an SAEM model
#'
#' Configure an SAEM model by generating an input list to the SAEM model function
#'
#' @param model a compiled saem model by gen_saem_user_fn()
#' @param data input data
#' @param inits initial values
#' @param mcmc a list of various mcmc options
#' @param ODEopt optional ODE solving options
#' @param seed seed for random number generator
#' @param distribution one of c("normal","poisson","binomial")
#' @param fixed a character vector of fixed effect only parameters (no random effects attached) to be fixed
#' @param DEBUG Integer determining if debugging is enabled
#' @param type indicates the type of optimization for the residuals; Can be one of c("nelder-mead", "newuoa")
#' @param lambdaRange This indicates the range that Box-Cox and Yeo-Johnson parameters are constrained to be; The default is 3 indicating the range (-3,3)
#' @param powRange This indicates the range that powers can take for residual errors; By default this is 10 indicating the range is c(1/10, 10) or c(0.1,10)
#' @inheritParams RxODE::rxSEinner
#' @inheritParams saemControl
#'
#' @return Returns a list neede for the saem fit procedure
#'
#' @details
#' Fit a generalized nonlinear mixed-effect model by he Stochastic
#' Approximation Expectation-Maximization (SAEM) algorithm
#'
#' @author Wenping Wang & Matthew Fidler
#' @examples
#' \donttest{
#'
#' # In this ODE system we simply specify the ODEs
#'
#' ode <- "d/dt(depot) =-KA*depot;
#' d/dt(centr) = KA*depot - KE*centr;"
#' m1 <- RxODE(ode)
#'
#'
#' # In this ode System, we also specify the concentration as C2 = centr/V
#'
#' ode <- "C2 = centr/V;
#' d/dt(depot) =-KA*depot;
#' d/dt(centr) = KA*depot - KE*centr;"
#' m2 = RxODE(ode)
#'
#' PKpars <- function() {
#' CL <- exp(lCL)
#' V <- exp(lV)
#' KA <- exp(lKA)
#' KE <- CL / V
#' }
#'
#' PRED <- function() centr / V
#' PRED2 <- function() C2
#'
#' saem_fit <- gen_saem_user_fn(model = m1, PKpars, pred = PRED)
#'
#' # Can also use PRED2
#' saem_fit <- gen_saem_user_fn(model=m2, PKpars, pred=PRED2)
#'
#'
#' # You can also use the nlmixr UI to run this model and call the lower level functions
#'
#' one.compartment <- function() {
#' ini({
#' tka <- 0.45 # Log Ka
#' tcl <- 1 # Log Cl
#' tv <- 3.45 # Log V
#' eta.ka ~ 0.6
#' eta.cl ~ 0.3
#' eta.v ~ 0.1
#' add.sd <- 0.7
#' wt.est <- 0.0
#' })
#' model({
#' ka <- exp(tka + eta.ka)
#' cl <- exp(tcl + eta.cl)
#' v <- exp(tv + eta.v + wt.est * WT)
#' d/dt(depot) = -ka * depot
#' d/dt(center) = ka * depot - cl / v * center
#' cp = center / v
#' cp ~ add(add.sd)
#' })
#' }
#' fit <- nlmixr(one.compartment, theo_sd, "saem")
#' fit
#'
#' }
#' @export
configsaem <- function(model, data, inits,
mcmc = list(niter = c(200, 300), nmc = 3, nu = c(2, 2, 2)),
ODEopt = list(atol = 1e-6, rtol = 1e-4, method = "lsoda", transitAbs = FALSE, maxeval = 100000),
distribution = c("normal", "poisson", "binomial"),
addProp = c("combined2", "combined1"),
seed = 99, fixed = NULL, DEBUG = 0,
tol = 1e-4, itmax = 100L, type = c("nelder-mead", "newuoa"),
lambdaRange = 3, powRange = 10,
odeRecalcFactor=10^(0.5),
maxOdeRecalc=5L) {
type.idx <- c("nelder-mead" = 1L, "newuoa" = 2L)
type <- match.arg(type)
type <- type.idx[type]
force(ODEopt)
names(ODEopt) <- gsub("transit_abs", "transitAbs", names(ODEopt))
ODEopt <- do.call(RxODE::rxSolve, c(list(object=NULL), ODEopt))
# mcmc=list(niter=c(200,300), nmc=3, nu=c(2,2,2));ODEopt = list(atol=1e-6, rtol=1e-4, stiff=1, transit_abs=0);distribution=c("normal","poisson","binomial");seed=99;data=dat;distribution=1;fixed=NULL
set.seed(seed)
distribution.idx <- c("normal" = 1, "poisson" = 2, "binomial" = 3)
distribution <- match.arg(distribution)
distribution <- distribution.idx[distribution]
.data <- data
## RxODE::rxTrans(data, model)
data <- list(nmdat = data)
neq <- attr(model$saem_mod, "neq")
nlhs <- attr(model$saem_mod, "nlhs")
inPars <- attr(model$saem_mod, "inPars")
ninputpars <- length(inPars)
opt <- optM <- c(list(neq = neq, nlhs = nlhs, inits = numeric(neq)),
ninputpars = ninputpars, inPars = inPars
)
model$N.eta <- attr(model$saem_mod, "nrhs")
model$nendpnt <- attr(model$saem_mod, "nendpnt")
if (is.null(model$nendpnt)) model$nendpnt <- 1
if (is.null(model$log.eta)) model$log.eta <- rep(TRUE, model$N.eta)
if (is.null(model$omega)) model$omega <- diag(model$N.eta)
if (is.null(model$res.mod)) model$res.mod <- rep(1, model$nendpnt)
if (is.null(inits$omega)) inits$omega <- rep(1, model$N.eta) * 4
if (is.null(inits$ares)) inits$ares <- 10
if (is.null(inits$bres)) inits$bres <- 1
if (is.null(inits$cres)) inits$cres <- 1
if (is.null(inits$lres)) inits$lres <- 1
if (is.null(mcmc$print)) mcmc$print <- 1
if (model$N.eta - length(inits$theta) > 0) {
inits$theta <- c(inits$theta, rep(NA_real_, model$N.eta - length(inits$theta)))
}
if (is.null(names(inits$theta))) names(inits$theta) <- rep("", length(inits$theta))
inits.save <- inits
inits$theta.fix <- matrix(names(inits$theta),
byrow = TRUE,
ncol = model$N.eta
)
inits$theta <- matrix(inits$theta, byrow = TRUE, ncol = model$N.eta)
model$cov.mod <- 1 - is.na(inits$theta)
data$N.covar <- nrow(inits$theta) - 1
inits$theta[is.na(inits$theta)] <- 0
### FIXME
mcmc$stepsize <- 0:1
mcmc$burn.in <- 300
### FIXME: chk covars as char vec
wh <- setdiff(c(model$covars, inPars), names(data$nmdat))
if (length(wh)) {
msg <- paste0("covariate(s) not found: ", paste(wh, collapse = ", "))
stop(msg)
}
s <- subset(data$nmdat, EVID == 0)
data$data <- as.matrix(s[, c("ID", "TIME", "DV", c(model$covars, inPars))])
### chk for no obs records
wh <- setdiff(unique(data$nmdat$ID), unique(data$data[, "ID"]))
if (length(wh)) {
msg <- paste0("No data with ID: ", paste(wh, collapse = ", "))
stop(msg)
}
nphi <- model$N.eta
mcov <- model$cov.mod
covstruct <- model$omega
check <- sum((covstruct - t(covstruct)) != 0)
if (check) stop("illegal covstruct")
check <- nphi - dim(covstruct)[1]
if (check) stop("nphi and covstruct dim mismatch")
check <- prod(mcov[1, ])
if (check == 0) stop("structural par(s) absent")
check <- nphi - dim(mcov)[2]
if (check) stop("nphi and ncol(mcov) mismatch")
check <- sum(dim(inits$theta) - dim(mcov) != 0)
if (check) stop("initial theta's and mcov dim mismatch")
check <- data$N.covar + 1 - dim(mcov)[1]
if (check) stop("dim mcov and N.covar mismatch")
check <- length(model$log.eta) - nphi
if (check) stop("jlog length and nphi mismatch")
check <- length(inits$omega) - nphi
if (check) stop("length of omega inits and nphi mismatch")
# check = mcmc$burn.in>sum(mcmc$niter)
# if (check) stop("#burn-in exceeds niter")
check <- prod(is.element(covstruct, c(0, 1)))
if (check == 0) warning("non-zero value(s) in covstruct set to 1")
covstruct[covstruct != 0] <- 1
check <- prod(is.element(mcov, c(0, 1)))
if (check == 0) warning("non-zero value(s) in mcov set to 1")
mcov[mcov != 0] <- 1
check <- sum(inits$theta[1, model$log.eta] <= 0)
if (check) stop("illegal initial theta's")
check <- sum(inits$omega <= 0)
if (check) stop("illegal initial omega")
# check = inits$sigma2<=0
# if (check) stop("illegal initial sigma2")
check <- sum(diag(covstruct) == 1)
if (!check) stop("0 ETA's")
y <- data$data[, "DV"]
id <- data$data[, "ID"]
check <- any(diff(unique(id)) != 1)
if (check) stop("saem classic UI needs sequential ID. check your data")
ntotal <- length(id)
N <- length(unique(id))
covariables <- if (is.null(model$covars)) NULL else unlist(stats::aggregate(.as.data.frame(data$data[, model$covars, drop = FALSE]), list(id), unique)[, -1, drop = FALSE])
if (!is.null(covariables)) dim(covariables) <- c(N, data$N.covar)
nb_measures <- table(id)
ncov <- data$N.covar + 1
nmc <- mcmc$nmc
nM <- mcmc$nmc * N
yM <- rep(y, nmc)
mlen <- max(nb_measures)
io <- t(sapply(nb_measures, function(x) rep(1:0, c(x, mlen - x))))
ix <- rep(1:dim(io)[1], nmc)
ioM <- io[ix, ]
indioM <- grep(1, t(ioM)) - 1
## mPars <- if (ninputpars == 0) NULL else unlist(stats::aggregate(.as.data.frame(data$data[, inPars]), list(id), unique)[, -1])
## if (!is.null(mPars)) {
## dim(mPars) <- c(N, ninputpars)
## opt$mPars <- mPars
## ix <- rep(1:dim(mPars)[1], nmc)
## optM$mPars <- mPars[ix, ]
## dim(optM$mPars) <- c(nmc * N, ninputpars)
## }
if (is.null(data$nmdat$CMT)) data$nmdat$CMT <- 1 ## CHECKME
if (any(is.na(data$nmdat$CMT))) {
stop("'CMT' has NA(s)")
}
## CHECKME
form <- attr(model$saem_mod, "form")
.nobs <- 0
dat <- RxODE::etTrans(data$nmdat, attr(model$saem_mod, "rx"), addCmt=TRUE, dropUnits=TRUE, allTimeVar=TRUE)
.nobs <- attr(class(dat), ".RxODE.lst")$nobs
## if(length(dat) !=7) stop("SAEM doesn't support time varying covariates yet.");
.rx <- attr(model$saem_mod, "rx")
.pars <- .rx$params
.pars <- setNames(rep(1.1, length(.pars)), .pars)
opt$.rx <- .rx
opt$.pars <- .pars
## opt$.dat <- dat;
dat <- .as.data.frame(dat[, -6])
names(dat) <- toupper(names(dat))
dat$ID <- as.integer(dat$ID)
evt <- dat
evt$ID <- evt$ID - 1
## r
evtM <- evt[rep(1:dim(evt)[1], nmc), ]
evtM$ID <- cumsum(c(FALSE, diff(evtM$ID) != 0))
i1 <- grep(1, diag(covstruct))
i0 <- grep(0, diag(covstruct))
nphi1 <- sum(diag(covstruct))
nphi0 <- nphi - nphi1
na <- length(mcmc$stepsize)
nlambda1 <- sum(mcov[, i1])
nlambda0 <- sum(mcov[, i0])
nlambda <- nlambda1 + nlambda0
nd1 <- nphi1 + nlambda1 + 1
nd2 <- nphi1 + nlambda1 + nlambda0
nb_param <- nd2 + 1
Mcovariables <- cbind(rep(1, N), covariables)[, 1:nrow(mcov)]
dim(Mcovariables) <- c(length(Mcovariables) / nrow(mcov), nrow(mcov)) # FIXME
# get fixed ix
fixed <- inits$theta.fix != ""
wh <- fixed[, i1][mcov[, i1] == 1]
len <- length(wh)
fixed.i1 <- (1:len)[wh] - 1
wh <- fixed[, i0][mcov[, i0] == 1]
len <- length(wh)
fixed.i0 <- (1:len)[wh] - 1
jlog1 <- grep(TRUE, model$log.eta)
jcov <- grep(TRUE, apply(mcov, 1, sum) > 0)
covstruct1 <- covstruct[i1, i1]
dim(covstruct1) <- c(nphi1, nphi1)
ind_cov <- grep(1, mcov[mcov > 0])
mcov1 <- matrix(mcov[, i1], ncol = length(i1))
mcov0 <- matrix(mcov[, i0], nrow = nrow(mcov), ncol = length(i0))
ind_cov1 <- grep(1, mcov1[mcov1 > 0]) - 1
ind_cov0 <- grep(1, mcov0[mcov0 > 0]) - 1
pc <- apply(mcov, 2, sum)
ipc <- cumsum(c(0, pc[1:(nphi - 1)])) + 1
ipcl1 <- ipc[jlog1]
for (x in jlog1) inits$theta[1, x] <- log(inits$theta[1, x])
idx <- as.vector(mcov1 > 0)
COV1 <- Mcovariables[, row(mcov1)[idx]]
dim(COV1) <- c(N, sum(idx))
COV21 <- crossprod(COV1)
x <- mcov1 * col(mcov1)
x <- sapply(x[idx], function(x) {
ret <- rep(0, nphi1)
ret[x] <- 1
ret
})
LCOV1 <- t(x)
dim(LCOV1) <- c(nlambda1, nphi1)
pc1 <- apply(LCOV1, 2, sum)
x1 <- diag(sum(idx))
diag(x1) <- inits$theta[, i1][idx]
MCOV1 <- x1 %*% LCOV1
jcov1 <- grep(1, LCOV1) - 1
idx <- as.vector(mcov0 > 0)
COV0 <- Mcovariables[, row(mcov0)[idx]]
dim(COV0) <- c(N, sum(idx))
COV20 <- crossprod(COV0)
x <- mcov0 * col(mcov0)
x <- sapply(x[idx], function(x) {
ret <- rep(0, nphi0)
ret[x] <- 1
ret
})
LCOV0 <- t(x)
dim(LCOV0) <- c(nlambda0, nphi0)
x1 <- diag(sum(idx))
diag(x1) <- inits$theta[, i0][idx]
if (dim(x1)[1] > 0) {
MCOV0 <- x1 %*% LCOV0
} else {
MCOV0 <- matrix(x1, nrow = 0, ncol = dim(LCOV0)[2])
}
jcov0 <- grep(1, LCOV0) - 1
mprior_phi1 <- Mcovariables %*% inits$theta[, i1]
mprior_phi0 <- Mcovariables %*% inits$theta[, i0]
Gamma2_phi1 <- diag(nphi1)
diag(Gamma2_phi1) <- inits$omega[i1]
Gamma2_phi0 <- diag(nphi0)
diag(Gamma2_phi0) <- inits$omega[i0]
phiM <- matrix(0, N, nphi)
phiM[, i1] <- mprior_phi1
phiM[, i0] <- mprior_phi0
phiM <- phiM[rep(1:N, nmc), , drop = FALSE]
.tmp <- diag(sqrt(inits$omega))
if (model$N.eta == 1) .tmp <- matrix(sqrt(inits$omega))
.mat2 <- matrix(rnorm(phiM), dim(phiM))
phiM <- phiM + .mat2 %*% .tmp
mc.idx <- rep(1:N, nmc)
statphi <- sapply(1:nphi, function(x) tapply(phiM[, x], mc.idx, mean))
statphi11 <- statphi[, i1]
dim(statphi11) <- c(N, length(i1))
statphi01 <- statphi[, i0]
dim(statphi01) <- c(N, length(i0))
statphi12 <- crossprod(phiM[, i1])
statphi02 <- crossprod(phiM[, i0])
# x = mcov *cumsum(mcov)
# x1 = cbind(x[,i1], x[,i0])
# indiphi = order(x1[x1>0]) #FINDME
niter <- sum(mcmc$niter)
niter_phi0 <- niter * .5
nb_sa <- mcmc$niter[1] * .75
nb_correl <- mcmc$niter[1] * .75
va <- mcmc$stepsize
vna <- mcmc$niter
na <- length(va)
pas <- 1 / (1:vna[1])^va[1]
for (ia in 2:na) {
end <- length(pas)
k1 <- pas[end]^(-1 / va[ia])
pas <- c(pas, 1 / ((k1 + 1):(k1 + vna[ia]))^va[ia])
}
pash <- c(rep(1, mcmc$burn.in), 1 / (1:niter))
minv <- rep(1e-20, nphi)
# preserve par order when printing iter history
mcov[mcov == 1] <- 1:nlambda
ilambda1 <- mcov[, i1]
ilambda1 <- ilambda1[ilambda1 > 0] - 1
ilambda0 <- mcov[, i0]
ilambda0 <- ilambda0[ilambda0 > 0] - 1
i1 <- i1 - 1
i0 <- i0 - 1
opt$distribution <- distribution
opt$paramUpdate <- attr(model$saem_mod, "paramUpdate")
optM$paramUpdate <- attr(model$saem_mod, "paramUpdate")
opt$ODEopt <- ODEopt
optM$ODEopt <- ODEopt
cfg <- list(
ODEopt = ODEopt,
inits = inits.save,
nu = mcmc$nu,
niter = niter,
nb_sa = nb_sa,
nb_correl = nb_correl,
niter_phi0 = niter_phi0,
nmc = nmc,
coef_phi0 = .9638, # FIXME
rmcmc = .5,
coef_sa = .95,
pas = pas,
pash = pash,
minv = minv,
N = N,
ntotal = ntotal,
y = y,
yM = yM,
phiM = phiM,
evt = as.matrix(evt),
evtM = as.matrix(evtM),
mlen = mlen,
indioM = indioM,
pc1 = pc1,
covstruct1 = covstruct1,
Mcovariables = Mcovariables,
i1 = i1,
i0 = i0,
nphi1 = nphi1,
nphi0 = nphi0,
nlambda1 = nlambda1,
nlambda0 = nlambda0,
COV0 = COV0,
COV1 = COV1,
COV20 = COV20,
COV21 = COV21,
LCOV0 = LCOV0,
LCOV1 = LCOV1,
MCOV0 = MCOV0,
MCOV1 = MCOV1,
Gamma2_phi0 = Gamma2_phi0,
Gamma2_phi1 = Gamma2_phi1,
mprior_phi0 = mprior_phi0,
mprior_phi1 = mprior_phi1,
jcov0 = jcov0,
jcov1 = jcov1,
ind_cov0 = ind_cov0,
ind_cov1 = ind_cov1,
statphi11 = statphi11,
statphi01 = statphi01,
statphi02 = statphi02,
statphi12 = statphi12,
res.mod = model$res.mod,
ares = inits$ares,
bres = inits$bres,
cres = inits$cres,
lres = inits$lres,
opt = opt,
optM = optM,
print = mcmc$print,
distribution = distribution,
par.hist = matrix(0, sum(niter), nlambda1 + nlambda0 + nphi1 + 1 + (model$res.mod == 2) + 2 * (model$res.mod == 4)),
seed = seed,
fixed.i1 = fixed.i1,
fixed.i0 = fixed.i0,
ilambda1 = as.integer(ilambda1),
ilambda0 = as.integer(ilambda0),
nobs = .nobs)
## CHECKME
s <- cfg$evt[cfg$evt[, "EVID"] == 0, "CMT"]
cfg$opt$cmt_endpnt <- cfg$optM$cmt_endpnt <- sort(unique(s))
cfg$nendpnt <- length(unique(s))
if (model$nendpnt != cfg$nendpnt) {
msg <- sprintf("mis-match in nbr endpoints in model & in data")
stop(msg)
}
t <- unlist(split(1L:length(s), s))
cfg$ysM <- rep(cfg$y[t], cfg$nmc)
cfg$ix_sorting <- t - 1 # c-index for sorting by endpnt
cfg$y_offset <- c(0, cumsum(table(s)))
s <- cfg$evtM[cfg$evtM[, "EVID"] == 0, "CMT"]
cfg$ix_endpnt <- as.integer(as.factor(s)) - 1 # to derive vecares & vecbres
s <- cfg$evtM[cfg$evtM[, "EVID"] == 0, "ID"]
t <- cumsum(c(0, table(s)))
cfg$ix_idM <- cbind(t[-length(t)], t[-1] - 1) # c-index of obs records of each subject
cfg$ares <- rep(10, cfg$nendpnt)
cfg$bres <- rep(1, cfg$nendpnt)
cfg$cres <- rep(1, cfg$nendpnt)
cfg$lres <- rep(1, cfg$nendpnt)
cfg$yj <- rep(2L, cfg$nendpnt)
cfg$propT <- rep(0L, cfg$nendpnt)
cfg$lambda <- rep(1.0, cfg$nendpnt)
cfg$low <- rep(0.0, cfg$nendpnt)
cfg$hi <- rep(1.0, cfg$nendpnt)
cfg$ares[cfg$res.mod == 2] <- 0
cfg$bres[cfg$res.mod == 1] <- 0
nres <- (1:4)[(cfg$res.mod == 10L) * 3 + (cfg$res.mod %in% c(4L, 8L, 9L)) * 2 + (cfg$res.mod %in% c(3L, 5L, 6L, 7L)) + 1]
cfg$res_offset <- cumsum(c(0, nres))
cfg$par.hist <- matrix(0, cfg$niter, nlambda1 + nlambda0 + nphi1 + sum(nres))
cfg$addProp <- c("combined1" = 1L, "combined2" = 2L)[match.arg(addProp)]
cfg$DEBUG <- cfg$opt$DEBUG <- cfg$optM$DEBUG <- DEBUG
cfg$phiMFile <- tempfile("phi-", RxODE::rxTempDir(), ".phi")
cfg$tol <- tol
cfg$itmax <- itmax
cfg$type <- type
cfg$lambdaRange <- lambdaRange
cfg$powRange <- powRange
cfg$odeRecalcFactor <- odeRecalcFactor
cfg$maxOdeRecalc <- maxOdeRecalc
cfg
}
#' Print an SAEM model fit summary
#'
#' Print an SAEM model fit summary
#'
#' @param object a saemFit object
#' @param ... others
#' @return a list
#' @export
summary.saemFit <- function(object, ...) {
fit <- object ## Rcheck hack
th <- fit$Plambda
nth <- length(th)
H <- solve(fit$Ha[1:nth, 1:nth])
se <- sqrt(diag(H))
m <- cbind(exp(th), th, se) # FIXME
## lhsVars = scan("LHS_VARS.txt", what="", quiet=TRUE)
## if (length(lhsVars)==nth) dimnames(m)[[1]] = lhsVars
dimnames(m)[[2]] <- c("th", "log(th)", "se(log_th)")
cat("THETA:\n")
print(m)
cat("\nOMEGA:\n")
print(fit$Gamma2_phi1)
if (any(fit$sig2 == 0)) {
cat("\nSIGMA:\n")
print(max(fit$sig2^2))
} else {
cat("\nARES & BRES:\n")
print(fit$sig2)
}
invisible(list(theta = th, se = se, H = H, omega = fit$Gamma2_phi1, eta = fit$mpost_phi))
}
#' Print an SAEM model fit summary
#'
#' Print an SAEM model fit summary
#'
#' @param x a saemFit object
#' @param ... others
#' @return a list
#' @export
print.saemFit <- function(x, ...) {
fit <- x ## Rcheck hack
th <- fit$Plambda
nth <- length(th)
H <- solve(fit$Ha[1:nth, 1:nth])
se <- sqrt(diag(H))
m <- cbind(exp(th), th, se) # FIXME
## lhsVars = scan("LHS_VARS.txt", what="", quiet=TRUE)
## if (length(lhsVars)==nth) dimnames(m)[[1]] = lhsVars
dimnames(m)[[2]] <- c("th", "log(th)", "se(log_th)")
cat("THETA:\n")
print(m)
cat("\nOMEGA:\n")
print(fit$Gamma2_phi1)
if (any(fit$sig2 == 0)) {
cat("\nSIGMA:\n")
print(max(fit$sig2^2))
} else {
cat("\nARES & BRES:\n")
print(fit$sig2)
}
invisible(list(theta = th, se = se, H = H, omega = fit$Gamma2_phi1, eta = fit$mpost_phi))
}
#' Fit an SAEM model
#'
#' Fit an SAEM model using either closed-form solutions or ODE-based model definitions
#'
#' @param model an RxODE model or lincmt()
#' @param data input data
#' @param inits initial values
#' @param PKpars PKpars function
#' @param pred pred function
#' @param covars Covariates in data
#' @param mcmc a list of various mcmc options
#' @param ODEopt optional ODE solving options
#' @inheritParams configsaem
#' @param seed seed for random number generator
#' @details
#' Fit a generalized nonlinear mixed-effect model using the Stochastic
#' Approximation Expectation-Maximization (SAEM) algorithm
#'
#' @return saem fit object
#'
#' @author Matthew Fidler & Wenping Wang
#' @export
saem.fit <- function(model, data, inits,
PKpars = NULL, pred = NULL,
covars = NULL,
mcmc = list(niter = c(200, 300), nmc = 3, nu = c(2, 2, 2)),
ODEopt = list(atol = 1e-06, rtol = 1e-04, method = "lsoda", transitAbs = FALSE),
distribution = c("normal", "poisson", "binomial", "lnorm"),
seed = 99) {
UseMethod("saem.fit")
}
##' @rdname saem.fit
##' @export
saem <- saem.fit
##' @rdname saem.fit
##' @export
saem.fit.nlmixr.ui.nlme <- function(model, data, inits,
PKpars = NULL, pred = NULL,
covars = NULL,
mcmc = list(niter = c(200, 300), nmc = 3, nu = c(2, 2, 2)),
ODEopt = list(atol = 1e-06, rtol = 1e-04, method = "lsoda", transitAbs = FALSE),
distribution = c("normal", "poisson", "binomial", "lnorm"),
seed = 99) {
call <- as.list(match.call(expand.dots = TRUE))[-1]
names(call)[1] <- "object"
call$est <- "saem"
return(do.call(getFromNamespace("nlmixr", "nlmixr"), call, envir = parent.frame(1)))
}
##' @rdname saem.fit
##' @export
saem.fit.function <- saem.fit.nlmixr.ui.nlme
##' @rdname saem.fit
##' @export
saem.fit.nlmixrUI <- saem.fit.nlmixr.ui.nlme
##' @rdname saem.fit
##' @export
saem.fit.RxODE <- function(model, data, inits,
PKpars = NULL, pred = NULL,
covars = NULL,
mcmc = list(niter = c(200, 300), nmc = 3, nu = c(2, 2, 2)),
ODEopt = list(atol = 1e-06, rtol = 1e-04, method = "lsoda", transitAbs = FALSE),
distribution = c("normal", "poisson", "binomial", "lnorm"),
seed = 99) {
saem_fit <- gen_saem_user_fn(model, PKpars, pred)
model <- list(saem_mod = saem_fit, covars = covars)
cfg <- configsaem(model, data, inits, mcmc, ODEopt, distribution, seed)
fit <- saem_fit(cfg)
## dyn.unload("saem_main.dll")
fit
}
##' @rdname saem.fit
##' @export
saem.fit.default <- function(model, data, inits,
PKpars = NULL, pred = NULL,
covars = NULL,
mcmc = list(niter = c(200, 300), nmc = 3, nu = c(2, 2, 2)),
ODEopt = list(atol = 1e-06, rtol = 1e-04, method = "lsoda", transitAbs = FALSE),
distribution = c("normal", "poisson", "binomial", "lnorm"),
seed = 99) {
saem_fit <- gen_saem_user_fn(model)
model <- list(saem_mod = saem_fit, covars = covars)
cfg <- configsaem(model, data, inits, mcmc, ODEopt, distribution, seed)
fit <- saem_fit(cfg)
# dyn.unload("saem_main.dll")
fit
}
##' @export
ranef.saemFit <- function(object, ...) {
return(object$eta)
}
##' @export
fixef.saemFit <- function(object, ...) {
return(object$Plambda)
}
focei.theta.saemFit <- function(object, uif, ...) {
## Get the thetas needed for FOCEi fit.
this.env <- environment()
if (inherits(uif, "function")) {
uif <- nlmixr(uif)
}
n <- uif$focei.names
thetas <- structure(rep(NA, length(n)), .Names = n)
sf <- structure(as.vector(fixed.effects(object)), .Names = uif$saem.theta.name)
for (n in names(sf)) {
thetas[n] <- sf[n]
}
.predDf <- uif$predDf
.ini <- .as.data.frame(uif$ini)
.resMat <- object$resMat
sapply(seq_along(.predDf$cond), function(i) {
x <- paste(.predDf$cond[i])
.tmp <- .ini[which(.ini$condition == x), ]
.w <- which(sapply(.tmp$err, function(x) any(x == c("prop", "propT", "pow", "powT"))))
if (length(.w) == 1) {
thetas[paste(.tmp$name[.w])] <- .resMat[i, 2]
}
.w <- which(sapply(.tmp$err, function(x) any(x == c("pow2", "powT2"))))
if (length(.w) == 1) {
thetas[paste(.tmp$name[.w])] <- .resMat[i, 3]
}
.w <- which(sapply(seq_along(.tmp$err), function(.x) {
x <- .tmp$err[.x]
if (any(x == c(
"add", "norm", "dnorm", "lnorm", "dlnorm",
"dlogn", "logn"
))) {
if (!is.na(.tmp$est[.x])) {
return(TRUE)
}
}
return(FALSE)
}))
if (length(.w) == 1) {
thetas[paste(.tmp$name[.w])] <- .resMat[i, 1]
}
.w <- which(sapply(.tmp$err, function(x) {
any(x == c("boxCox", "yeoJohnson"))
}))
if (length(.w) == 1) {
thetas[paste(.tmp$name[.w])] <- .resMat[i, 4]
}
assign("thetas", thetas, this.env)
})
return(thetas)
}
focei.eta.saemFit <- function(object, uif, ...) {
if (inherits(uif, "function")) {
uif <- nlmixr(uif)
}
## Reorder based on translation
eta.trans <- uif$saem.eta.trans
for (i in seq(1, max(eta.trans))) {
while (!(any(i == eta.trans)) && max(eta.trans) > i) {
eta.trans[eta.trans >= i] <- eta.trans[eta.trans >= i] - 1
}
}
## orig eta -> new eta
df <- .as.data.frame(uif$ini)
eta <- df[!is.na(df$neta1), ]
len <- length(eta$name)
cur.lhs <- character()
cur.rhs <- numeric()
ome <- character()
cur.ome <- object$Gamma2_phi1
for (i in seq_along(eta$name)) {
last.block <- FALSE
if (i == len) {
last.block <- TRUE
} else if (eta$neta1[i + 1] == eta$neta2[i + 1]) {
last.block <- TRUE
}
if (eta$neta1[i] == eta$neta2[i]) {
cur.lhs <- c(cur.lhs, sprintf("ETA[%d]", eta$neta1[i]))
cur.rhs <- c(cur.rhs, cur.ome[eta.trans[eta$neta1[i]], eta.trans[eta$neta2[i]]])
if (last.block) {
ome[length(ome) + 1] <- sprintf(
"%s ~ %s", paste(cur.lhs, collapse = " + "),
paste(deparse(cur.rhs), collapse = " ")
)
cur.lhs <- character()
cur.rhs <- numeric()
}
} else {
cur.rhs <- c(cur.rhs, cur.ome[eta.trans[eta$neta1[i]], eta.trans[eta$neta2[i]]])
}
}
ome <- eval(parse(text = sprintf("list(%s)", paste(ome, collapse = ","))))
return(ome)
}
as.focei.saemFit <- function(object, uif, pt = proc.time(), ..., data, calcResid = TRUE, obf = NULL,
nnodes.gq = 1, nsd.gq = 3, adjObf = TRUE,
calcCov = TRUE, covMethod = NULL,
calcCovTime = NULL, calcTables = TRUE, keep=NULL, drop=NULL, table=tableControl(),
IDlabel=NULL) {
if (is.null(calcResid)) calcResid <- FALSE
.saemCfg <- attr(object, "saem.cfg")
.saemTime <- proc.time() - pt
if (inherits(uif, "function")) {
uif <- nlmixr(uif)
}
.dist <- ""
if (any(uif$saem.distribution == c("poisson", "binomial"))) {
calcResid <- NA
.dist <- uif$saem.distribution
}
uif.new <- uif
fit <- object
mat <- random.effects(fit)
## Reorder based on translation
eta.trans <- uif$saem.eta.trans
for (i in seq(1, max(eta.trans))) {
while (!(any(i == eta.trans)) && max(eta.trans) > i) {
eta.trans[eta.trans >= i] <- eta.trans[eta.trans >= i] - 1
}
}
mat2 <- mat[, eta.trans, drop = FALSE]
th <- focei.theta(fit, uif)
for (n in names(th)) {
uif.new$est[uif.new$name == n] <- th[n]
}
ome <- focei.eta(fit, uif)
init <- list(
THTA = as.vector(th),
OMGA = ome
)
saem.time <- proc.time() - pt
if (missing(data)) {
stop("Requires Data...")
} else {
dat <- data
}
.tn <- uif$saem.theta.name
.ini <- .as.data.frame(uif$ini)
.ini <- .ini[uif$ini$name %in% .tn, ]
if (any(.ini$fix)) {
.fixed <- paste(.ini$name[.ini$fix])
.tn <- .tn[!(.tn %in% .fixed)]
}
.calcCov <- calcCov
if (uif$env$covMethod == "") {
.cov <- NULL
.addCov <- FALSE
} else if (inherits(calcCov, "matrix")) {
.cov <- calcCov
.addCov <- TRUE
} else {
.nth <- length(.tn)
.ini <- .as.data.frame(uif$ini)
.ini <- .ini[is.na(.ini$err), ]
.ini <- .ini[!is.na(.ini$ntheta), ]
.ini <- .ini[!.ini$fix, ]
.ini <- paste(.ini$name)
.calcCovTime <- proc.time()
if (calcCov) {
.covm <- object$Ha[1:.nth, 1:.nth]
.covm <- try(calc.COV(object))
.doIt <- !inherits(.covm, "try-error")
if (.doIt && dim(.covm)[1] != .nth) .doIt <- FALSE
if (.doIt) {
.tmp <- try(chol(.covm), silent = TRUE)
.addCov <- TRUE
.sqrtm <- FALSE
if (inherits(.tmp, "try-error")) {
.tmp <- .covm
.tmp <- try(sqrtm(.tmp %*% t(.tmp)), silent = FALSE)
if (inherits(.tmp, "try-error")) {
.calcCov <- FALSE
.covm <- object$Ha[1:.nth, 1:.nth]
.tmp <- try(chol(.covm), silent = TRUE)
.addCov <- TRUE
.sqrtm <- FALSE
if (inherits(.tmp, "try-error")) {
.tmp <- object$Ha[1:.nth, 1:.nth]
.tmp <- try(sqrtm(.tmp %*% t(.tmp)), silent = FALSE)
if (inherits(.tmp, "try-error")) {
.addCov <- FALSE
} else {
.sqrtm <- TRUE
}
} else {
.tmp <- object$Ha[1:.nth, 1:.nth]
}
} else {
.sqrtm <- TRUE
}
} else {
.tmp <- .covm
}
} else {
.tmp <- object$Ha[1:.nth, 1:.nth]
.tmp <- try(chol(.covm), silent = TRUE)
.calcCov <- FALSE
.addCov <- TRUE
.sqrtm <- FALSE
if (inherits(.tmp, "try-error")) {
.tmp <- object$Ha[1:.nth, 1:.nth]
.tmp <- try(sqrtm(.tmp %*% t(.tmp)), silent = FALSE)
if (inherits(.tmp, "try-error")) {
.addCov <- FALSE
} else {
.sqrtm <- TRUE
}
} else {
.tmp <- object$Ha[1:.nth, 1:.nth]
.calcCov <- FALSE
}
}
} else {
.tmp <- try(chol(.covm), silent = TRUE)
.addCov <- TRUE
.sqrtm <- FALSE
if (inherits(.tmp, "try-error")) {
.tmp <- object$Ha[1:.nth, 1:.nth]
.tmp <- try(sqrtm(.tmp %*% t(.tmp)), silent = FALSE)
if (inherits(.tmp, "try-error")) {
.addCov <- FALSE
} else {
.sqrtm <- TRUE
}
} else {
.tmp <- object$Ha[1:.nth, 1:.nth]
.calcCov <- FALSE
}
}
if (.addCov) {
if (!.calcCov) {
.cov <- RxODE::rxInv(.tmp)
} else {
.cov <- .tmp
}
attr(.cov, "dimnames") <- list(.tn, .tn)
.cov <- .cov[.ini, .ini, drop = FALSE]
}
.calcCovTime <- proc.time() - .calcCovTime
.calcCovTime <- .calcCovTime["elapsed"]
}
.ini <- .as.data.frame(uif$ini)
.ini <- .ini[!is.na(.ini$ntheta), ]
.skipCov <- !is.na(.ini$err)
.fixed <- uif$focei.fixed
.skipCov <- .skipCov | .fixed[seq_along(.skipCov)]
.covMethod <- uif$env$covMethod
if (!any(.covMethod == c("", "r", "s", "r,s"))) {
.covMethod <- ""
}
if (is.na(calcResid)) .covMethod <- ""
.allThetaNames <- c(uif$saem.theta.name, uif$saem.omega.name, uif$saem.res.name)
.m <- object$par_hist
if (ncol(.m) > length(.allThetaNames)) {
.m <- .m[, seq_along(.allThetaNames)]
}
if (.dist == "binomial") {
.dist <- "bernoulli"
}
dimnames(.m) <- list(NULL, .allThetaNames)
.fixedNames <- paste(uif$ini$name[which(uif$ini$fix)])
.rn <- ""
.likTime <- 0
if (is.na(obf)) {
.saemObf <- NA
} else if (is.null(obf)) {
.likTime <- proc.time()
.saemObf <- calc.2LL(object, nnodes.gq = nnodes.gq, nsd.gq = nsd.gq)
.likTime <- proc.time() - .likTime
.likTime <- .likTime["elapsed"]
if (nnodes.gq == 1) {
.rn <- paste0("laplace", nsd.gq)
} else {
.rn <- paste0("gauss", nnodes.gq, "_", nsd.gq)
}
} else if (is(obf, "logical")) {
if (is.na(obf)) {
.saemObf <- NA
} else if (obf) {
.likTime <- proc.time()
.saemObf <- calc.2LL(object, nnodes.gq = nnodes.gq, nsd.gq = nsd.gq)
.likTime <- proc.time() - .likTime
.likTime <- .likTime["elapsed"]
if (nnodes.gq == 1) {
.rn <- paste0("laplace", nsd.gq)
} else {
.rn <- paste0("gauss", nnodes.gq, "_", nsd.gq)
}
} else {
.saemObf <- NA
}
} else if (is(object, "numeric")) {
.saemObf <- obf
}
.notCalced <- TRUE
.cwresTime <- proc.time()
while (.notCalced) {
.env <- new.env(parent = emptyenv())
.env$table <- table
.env$IDlabel <- IDlabel
.env$nobs2 <- .saemCfg$nobs
.env$nnodes.gq <- nnodes.gq
.env$nsd.gq <- nsd.gq
.env$adjObf <- adjObf
.env$method <- "SAEM"
.env$uif <- uif
.env$saem <- object
if (.addCov) {
.env$cov <- .cov
}
.env$parHistStacked <- .data.frame(
val = as.vector(.m),
par = rep(.allThetaNames, each = nrow(.m)),
iter = rep(1:nrow(.m), ncol(.m))
)
.env$parHist <- .data.frame(iter = rep(1:nrow(.m)), .as.data.frame(.m))
if (length(.fixedNames) > 0) {
.env$parHistStacked <- .env$parHistStacked[!(.env$parHistStacked$par %in% .fixedNames), , drop = FALSE]
.env$parHist <- .env$parHist[, !(names(.env$parHist) %in% .fixedNames), drop = FALSE]
}
if (is.na(calcResid)) {
.setSaemExtra(.env, .rn)
.env$theta <- .data.frame(
lower = -Inf, theta = init$THTA, upper = Inf, fixed = .fixed[seq_along(init$THTA)],
row.names = uif$focei.names
)
.env$fullTheta <- setNames(init$THTA, uif$focei.names)
.om0 <- .genOM(.parseOM(init$OMGA))
attr(.om0, "dimnames") <- list(uif$eta.names, uif$eta.names)
.env$omega <- .om0
.env$etaObf <- .data.frame(ID = seq_along(mat2[, 1]), setNames(.as.data.frame(mat2), uif$eta.names), OBJI = NA)
.env$noLik <- FALSE
.env$objective <- .saemObf
} else if (calcResid) {
.setSaemExtra(.env, .rn)
} else {
.setSaemExtra(.env, .rn)
.env$theta <- .data.frame(
lower = -Inf, theta = init$THTA, upper = Inf, fixed = .fixed[seq_along(init$THTA)],
row.names = uif$focei.names
)
.env$fullTheta <- setNames(init$THTA, uif$focei.names)
.om0 <- .genOM(.parseOM(init$OMGA))
attr(.om0, "dimnames") <- list(uif$eta.names, uif$eta.names)
.env$omega <- .om0
.env$etaObf <- .data.frame(ID = seq_along(mat2[, 1]), setNames(.as.data.frame(mat2), uif$eta.names), OBJI = NA)
.env$noLik <- TRUE
.env$objective <- .saemObf
}
.ctl <- uif$env$ODEopt
names(.ctl) <- sub("maxsteps", "maxstepsOde", names(.ctl))
.ctl <- .ctl[names(.ctl) != "scale"]
.ctl$maxOuterIterations <- 0
.ctl$maxInnerIterations <- 0
.ctl$covMethod <- .covMethod
.ctl$sumProd <- uif$env$sum.prod
.ctl$optExpression <- uif$env$optExpression
.ctl$scaleTo <- 0
.ctl$calcTables <- calcTables
if (.saemCfg$addProp == 1L) {
.ctl$addProp <- "combined1"
} else {
.ctl$addProp <- "combined2"
}
.ctl$method <- "liblsoda"
.ctl <- do.call(foceiControl, .ctl)
if (.ctl$singleOde) {
.pars <- NULL
.mod <- uif$focei.rx1
} else {
.pars <- uif$theta.pars
.mod <- uif$rxode.pred
}
fit.f <- try(foceiFit.data.frame(
data = dat, inits = init, PKpars = .pars,
model = .mod, pred = function() {
return(nlmixr_pred)
}, err = uif$error,
lower = uif$focei.lower,
upper = uif$focei.upper,
thetaNames = uif$focei.names,
etaNames = uif$eta.names,
etaMat = mat2,
env = .env,
fixed = .fixed,
skipCov = .skipCov,
control = .ctl,
keep=keep, drop=drop
), silent = FALSE)
if (inherits(fit.f, "try-error")) {
if (is.na(calcResid)) {
warning("Error calculating nlmixr object, return classic object")
.notCalced <- FALSE
return(object)
} else if (calcResid) {
calcResid <- FALSE
} else {
calcResid <- NA
}
} else {
.notCalced <- FALSE
}
}
.cwresTime <- proc.time() - .cwresTime
if (is.na(calcResid)) {
.cwresTime <- 0
} else if (!calcResid) {
.cwresTime <- 0
}
.env <- fit.f$env
if (uif$env$covMethod == "") {
} else if (inherits(.calcCov, "matrix")) {
if (!is.null(covMethod)) {
.env$covMethod <- covMethod
}
.calcCovTime <- calcCovTime
} else if (.calcCov) {
.env$covMethod <- "linFim"
if (.addCov & .sqrtm) {
.env$covMethod <- "|linFim|"
warning("Covariance matrix non-positive definite, corrected by sqrtm(linFim %*% linFim)")
}
} else {
if (calcCov) {
warning("Linearization of FIM could not be used to calculate covariance.")
}
if (.addCov & .sqrtm) {
.env$covMethod <- "|fim|"
warning("Covariance matrix non-positive definite, corrected by sqrtm(fim %*% fim)")
} else if (!.addCov) {
warning("FIM non-positive definite and cannot be used to calculate the covariance")
}
}
if (is.null(.env$time)) {
.env$time <- .data.frame(saem = .saemTime["elapsed"], check.names = FALSE, row.names = c(""))
} else {
.time <- .env$time
.time <- .time[, !(names(.time) %in% c("optimize", "covariance"))]
.saemTime <- .saemTime["elapsed"]
if (calcResid) {
.saemTime <- .saemTime - .cwresTime["elapsed"]
.time <- .data.frame(.time, cwres = .cwresTime["elapsed"], check.names = FALSE)
}
if (.likTime > 0) {
.time <- .data.frame(.time, logLik = .likTime, check.names = FALSE)
.saemTime <- .saemTime - .likTime
}
if (uif$env$covMethod != "") {
.saemTime <- .saemTime - .calcCovTime
.time <- .data.frame(.time, covariance = .calcCovTime, check.names = FALSE)
}
.env$time <- .data.frame(saem = .saemTime, .time, check.names = FALSE, row.names = c(""))
}
.env$message <- ""
if (is.na(calcResid)) {
row.names(.env$objDf) <- .rn
} else if (calcResid) {
if (!is.na(.saemObf)) {
.llik <- -.saemObf / 2
.nobs <- .env$nobs
attr(.llik, "df") <- attr(get("logLik", .env), "df")
.objf <- ifelse(.env$adjObf, .saemObf - .nobs * log(2 * pi), .saemObf)
## for (.t in c("OBJF","objective", "objf")){
## assign(.t,.objf,.env);
## }
.tmp <- .data.frame(
OBJF = .objf,
AIC = .saemObf + 2 * attr(get("logLik", .env), "df"),
BIC = .saemObf + log(.env$nobs) * attr(get("logLik", .env), "df"),
"Log-likelihood" = as.numeric(.llik), check.names = FALSE
)
if (any(names(.env$objDf) == "Condition Number")) {
.tmp <- .data.frame(.tmp,
"Condition Number" = .env$objDf[, "Condition Number"],
check.names = FALSE
)
}
.env$objDf <- rbind(.env$objDf, .tmp)
row.names(.env$objDf) <- c("FOCEi", .rn)
}
.setSaemExtra(.env, "FOCEi")
} else {
row.names(.env$objDf) <- .rn
}
if (inherits(fit.f, "nlmixrFitData")) {
.cls <- class(fit.f)
.env <- attr(.cls, ".foceiEnv")
.cls[1] <- "nlmixrSaem"
class(fit.f) <- .cls
}
assign("uif", .syncUif(fit.f$uif, fit.f$popDf, fit.f$omega), fit.f$env)
return(fit.f)
}
.saemResidF <- function(x) {
.Call(`_saemResidF`, x)
}
.newuoa <- function(par, fn, gr, lower = -Inf, upper = Inf, control = list(), ...) {
.ctl <- control
if (is.null(.ctl$npt)) .ctl$npt <- length(par) * 2 + 1
if (is.null(.ctl$rhobeg)) .ctl$rhobeg <- 0.2
if (is.null(.ctl$rhoend)) .ctl$rhoend <- 1e-4
.ctl$iprint <- 0L
.ctl <- .ctl[names(.ctl) %in% c("npt", "rhobeg", "rhoend", "iprint", "maxfun")]
.ret <- minqa::newuoa(par, fn,
control = .ctl
)
.ret$x <- .ret$par
.ret$message <- .ret$msg
.ret$convergence <- .ret$ierr
.ret$value <- .ret$fval
return(.ret)
}
## FIXME: coef_phi0, rmcmc, coef_sa
## FIXME: Klog, rho, sa, nmc
## FIXME: N.design
## FIXME: g = gc = 1
## FIXME: ODE inits
## FIXME: Tinf for ODE
## FIXME: chk infusion poor fit
## Local Variables:
## ess-indent-level: 2
## indent-tabs-mode: nil
## End:
Try the nlmixr package in your browser
Any scripts or data that you put into this service are public.
nlmixr documentation built on March 27, 2022, 5:05 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions genSaemUserFunction
## saem_fit.R: population PK/PD modeling library
##
## Copyright (C) 2014 - 2016 Wenping Wang
##
## This file is part of nlmixr.
##
## nlmixr is free software: you can redistribute it and/or modify it
## under the terms of the GNU General Public License as published by
## the Free Software Foundation, either version 2 of the License, or
## (at your option) any later version.
##
## nlmixr is distributed in the hope that it will be useful, but
## WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with nlmixr. If not, see <http://www.gnu.org/licenses/>.
# genSaemUserFunction(f$rxode.pred, f$saem.pars, f$pred, f$error)
genSaemUserFunction <- function(model, PKpars = attr(model, "default.pars"), pred = NULL, err = NULL,
control = saemControl(), inPars = NULL) {
.x <- deparse(body(pred))
.len <- length(.x)
.x <- if (.x[1] == "{") .x[2:(.len - 1)] else .x
.len <- length(.x)
.nendpnt <- .len
.mod <- RxODE::RxODE(RxODE::rxGenSaem(model, function() {
return(nlmixr_pred)
}, PKpars,
sum.prod = control$sum.prod,
optExpression = control$optExpression,
loadSymengine=control$loadSymengine
))
.fnPred <- bquote(function(a, b, c) {
RxODE::rxLoad(.(.mod))
RxODE::rxLock(.(.mod))
RxODE::rxAllowUnload(FALSE)
on.exit({
RxODE::rxUnlock(.(.mod))
RxODE::rxAllowUnload(TRUE)
RxODE::rxSolveFree()
})
.Call(`_nlmixr_saem_do_pred`, a, b, c)
})
.fn <- bquote(function(a, b, c) {
RxODE::rxLoad(.(.mod))
RxODE::rxLock(.(.mod))
on.exit({
RxODE::rxUnlock(.(.mod))
RxODE::rxAllowUnload(TRUE)
RxODE::rxSolveFree()
})
if (missing(b) && missing(c)) {
.ret <- .Call(`_nlmixr_saem_fit`, a, PACKAGE = "nlmixr")
attr(.ret, "dopred") <- .(.fnPred)
return(.ret)
} else {
.curFn <- .(.fnPred)
return(.curFn(a, b, c))
}
})
.param <- RxODE::rxParam(.mod)
.inits <- names(RxODE::rxInits(.mod))
.nrhs <- length(.param) - length(.inits)
if (any(.param == "CMT")) {
inPars <- unique(c(inPars, "CMT"))
}
.parmUpdate <- sapply(.param, function(x) {
if (any(x == inPars)) {
return(0L)
} else if (any(x == .inits)) {
return(0L)
} else {
return(1L)
}
})
.fn <- eval(.fn)
attr(.fn, "form") <- "ode" ## Not sure this is necessary any more
attr(.fn, "neq") <- length(RxODE::rxState(.mod))
attr(.fn, "nlhs") <- length(RxODE::rxLhs(.mod))
attr(.fn, "nrhs") <- sum(.parmUpdate)
attr(.fn, "paramUpdate") <- .parmUpdate
attr(.fn, "rx") <- .mod
attr(.fn, "inPars") <- inPars
attr(.fn, "nendpnt") <- .nendpnt # not calculated; Is this a problem?
return(.fn)
}
#' Generate an SAEM model
#'
#' Generate an SAEM model using either closed-form solutions or ODE-based model definitions
#'
#' @param model a compiled SAEM model by gen_saem_user_fn()
#' @param PKpars PKpars function
#' @param pred pred function; This will be a focei-style pred
#' @param inPars a character vector of parameters to be read from the
#' input dataset (including time varying covariates)
#' @details Fit a generalized nonlinear mixed-effect model using the
#' Stochastic Approximation Expectation-Maximization (SAEM)
#' algorithm
#'
#' @return A user function based on the model to run the SAEM code
#'
#' @author Matthew Fidler & Wenping Wang
#' @keywords internal
#' @export
gen_saem_user_fn <- genSaemUserFunction
#' Configure an SAEM model
#'
#' Configure an SAEM model by generating an input list to the SAEM model function
#'
#' @param model a compiled saem model by gen_saem_user_fn()
#' @param data input data
#' @param inits initial values
#' @param mcmc a list of various mcmc options
#' @param ODEopt optional ODE solving options
#' @param seed seed for random number generator
#' @param distribution one of c("normal","poisson","binomial")
#' @param fixed a character vector of fixed effect only parameters (no random effects attached) to be fixed
#' @param DEBUG Integer determining if debugging is enabled
#' @param type indicates the type of optimization for the residuals; Can be one of c("nelder-mead", "newuoa")
#' @param lambdaRange This indicates the range that Box-Cox and Yeo-Johnson parameters are constrained to be; The default is 3 indicating the range (-3,3)
#' @param powRange This indicates the range that powers can take for residual errors; By default this is 10 indicating the range is c(1/10, 10) or c(0.1,10)
#' @inheritParams RxODE::rxSEinner
#' @inheritParams saemControl
#'
#' @return Returns a list neede for the saem fit procedure
#'
#' @details
#' Fit a generalized nonlinear mixed-effect model by he Stochastic
#' Approximation Expectation-Maximization (SAEM) algorithm
#'
#' @author Wenping Wang & Matthew Fidler
#' @examples
#' \donttest{
#'
#' # In this ODE system we simply specify the ODEs
#'
#' ode <- "d/dt(depot) =-KA*depot;
#' d/dt(centr) = KA*depot - KE*centr;"
#' m1 <- RxODE(ode)
#'
#'
#' # In this ode System, we also specify the concentration as C2 = centr/V
#'
#' ode <- "C2 = centr/V;
#' d/dt(depot) =-KA*depot;
#' d/dt(centr) = KA*depot - KE*centr;"
#' m2 = RxODE(ode)
#'
#' PKpars <- function() {
#' CL <- exp(lCL)
#' V <- exp(lV)
#' KA <- exp(lKA)
#' KE <- CL / V
#' }
#'
#' PRED <- function() centr / V
#' PRED2 <- function() C2
#'
#' saem_fit <- gen_saem_user_fn(model = m1, PKpars, pred = PRED)
#'
#' # Can also use PRED2
#' saem_fit <- gen_saem_user_fn(model=m2, PKpars, pred=PRED2)
#'
#'
#' # You can also use the nlmixr UI to run this model and call the lower level functions
#'
#' one.compartment <- function() {
#' ini({
#' tka <- 0.45 # Log Ka
#' tcl <- 1 # Log Cl
#' tv <- 3.45 # Log V
#' eta.ka ~ 0.6
#' eta.cl ~ 0.3
#' eta.v ~ 0.1
#' add.sd <- 0.7
#' wt.est <- 0.0
#' })
#' model({
#' ka <- exp(tka + eta.ka)
#' cl <- exp(tcl + eta.cl)
#' v <- exp(tv + eta.v + wt.est * WT)
#' d/dt(depot) = -ka * depot
#' d/dt(center) = ka * depot - cl / v * center
#' cp = center / v
#' cp ~ add(add.sd)
#' })
#' }
#' fit <- nlmixr(one.compartment, theo_sd, "saem")
#' fit
#'
#' }
#' @export
configsaem <- function(model, data, inits,
mcmc = list(niter = c(200, 300), nmc = 3, nu = c(2, 2, 2)),
ODEopt = list(atol = 1e-6, rtol = 1e-4, method = "lsoda", transitAbs = FALSE, maxeval = 100000),
distribution = c("normal", "poisson", "binomial"),
addProp = c("combined2", "combined1"),
seed = 99, fixed = NULL, DEBUG = 0,
tol = 1e-4, itmax = 100L, type = c("nelder-mead", "newuoa"),
lambdaRange = 3, powRange = 10,
odeRecalcFactor=10^(0.5),
maxOdeRecalc=5L) {
type.idx <- c("nelder-mead" = 1L, "newuoa" = 2L)
type <- match.arg(type)
type <- type.idx[type]
force(ODEopt)
names(ODEopt) <- gsub("transit_abs", "transitAbs", names(ODEopt))
ODEopt <- do.call(RxODE::rxSolve, c(list(object=NULL), ODEopt))
# mcmc=list(niter=c(200,300), nmc=3, nu=c(2,2,2));ODEopt = list(atol=1e-6, rtol=1e-4, stiff=1, transit_abs=0);distribution=c("normal","poisson","binomial");seed=99;data=dat;distribution=1;fixed=NULL
set.seed(seed)
distribution.idx <- c("normal" = 1, "poisson" = 2, "binomial" = 3)
distribution <- match.arg(distribution)
distribution <- distribution.idx[distribution]
.data <- data
## RxODE::rxTrans(data, model)
data <- list(nmdat = data)
neq <- attr(model$saem_mod, "neq")
nlhs <- attr(model$saem_mod, "nlhs")
inPars <- attr(model$saem_mod, "inPars")
ninputpars <- length(inPars)
opt <- optM <- c(list(neq = neq, nlhs = nlhs, inits = numeric(neq)),
ninputpars = ninputpars, inPars = inPars
)
model$N.eta <- attr(model$saem_mod, "nrhs")
model$nendpnt <- attr(model$saem_mod, "nendpnt")
if (is.null(model$nendpnt)) model$nendpnt <- 1
if (is.null(model$log.eta)) model$log.eta <- rep(TRUE, model$N.eta)
if (is.null(model$omega)) model$omega <- diag(model$N.eta)
if (is.null(model$res.mod)) model$res.mod <- rep(1, model$nendpnt)
if (is.null(inits$omega)) inits$omega <- rep(1, model$N.eta) * 4
if (is.null(inits$ares)) inits$ares <- 10
if (is.null(inits$bres)) inits$bres <- 1
if (is.null(inits$cres)) inits$cres <- 1
if (is.null(inits$lres)) inits$lres <- 1
if (is.null(mcmc$print)) mcmc$print <- 1
if (model$N.eta - length(inits$theta) > 0) {
inits$theta <- c(inits$theta, rep(NA_real_, model$N.eta - length(inits$theta)))
}
if (is.null(names(inits$theta))) names(inits$theta) <- rep("", length(inits$theta))
inits.save <- inits
inits$theta.fix <- matrix(names(inits$theta),
byrow = TRUE,
ncol = model$N.eta
)
inits$theta <- matrix(inits$theta, byrow = TRUE, ncol = model$N.eta)
model$cov.mod <- 1 - is.na(inits$theta)
data$N.covar <- nrow(inits$theta) - 1
inits$theta[is.na(inits$theta)] <- 0
### FIXME
mcmc$stepsize <- 0:1
mcmc$burn.in <- 300
### FIXME: chk covars as char vec
wh <- setdiff(c(model$covars, inPars), names(data$nmdat))
if (length(wh)) {
msg <- paste0("covariate(s) not found: ", paste(wh, collapse = ", "))
stop(msg)
}
s <- subset(data$nmdat, EVID == 0)
data$data <- as.matrix(s[, c("ID", "TIME", "DV", c(model$covars, inPars))])
### chk for no obs records
wh <- setdiff(unique(data$nmdat$ID), unique(data$data[, "ID"]))
if (length(wh)) {
msg <- paste0("No data with ID: ", paste(wh, collapse = ", "))
stop(msg)
}
nphi <- model$N.eta
mcov <- model$cov.mod
covstruct <- model$omega
check <- sum((covstruct - t(covstruct)) != 0)
if (check) stop("illegal covstruct")
check <- nphi - dim(covstruct)[1]
if (check) stop("nphi and covstruct dim mismatch")
check <- prod(mcov[1, ])
if (check == 0) stop("structural par(s) absent")
check <- nphi - dim(mcov)[2]
if (check) stop("nphi and ncol(mcov) mismatch")
check <- sum(dim(inits$theta) - dim(mcov) != 0)
if (check) stop("initial theta's and mcov dim mismatch")
check <- data$N.covar + 1 - dim(mcov)[1]
if (check) stop("dim mcov and N.covar mismatch")
check <- length(model$log.eta) - nphi
if (check) stop("jlog length and nphi mismatch")
check <- length(inits$omega) - nphi
if (check) stop("length of omega inits and nphi mismatch")
# check = mcmc$burn.in>sum(mcmc$niter)
# if (check) stop("#burn-in exceeds niter")
check <- prod(is.element(covstruct, c(0, 1)))
if (check == 0) warning("non-zero value(s) in covstruct set to 1")
covstruct[covstruct != 0] <- 1
check <- prod(is.element(mcov, c(0, 1)))
if (check == 0) warning("non-zero value(s) in mcov set to 1")
mcov[mcov != 0] <- 1
check <- sum(inits$theta[1, model$log.eta] <= 0)
if (check) stop("illegal initial theta's")
check <- sum(inits$omega <= 0)
if (check) stop("illegal initial omega")
# check = inits$sigma2<=0
# if (check) stop("illegal initial sigma2")
check <- sum(diag(covstruct) == 1)
if (!check) stop("0 ETA's")
y <- data$data[, "DV"]
id <- data$data[, "ID"]
check <- any(diff(unique(id)) != 1)
if (check) stop("saem classic UI needs sequential ID. check your data")
ntotal <- length(id)
N <- length(unique(id))
covariables <- if (is.null(model$covars)) NULL else unlist(stats::aggregate(.as.data.frame(data$data[, model$covars, drop = FALSE]), list(id), unique)[, -1, drop = FALSE])
if (!is.null(covariables)) dim(covariables) <- c(N, data$N.covar)
nb_measures <- table(id)
ncov <- data$N.covar + 1
nmc <- mcmc$nmc
nM <- mcmc$nmc * N
yM <- rep(y, nmc)
mlen <- max(nb_measures)
io <- t(sapply(nb_measures, function(x) rep(1:0, c(x, mlen - x))))
ix <- rep(1:dim(io)[1], nmc)
ioM <- io[ix, ]
indioM <- grep(1, t(ioM)) - 1
## mPars <- if (ninputpars == 0) NULL else unlist(stats::aggregate(.as.data.frame(data$data[, inPars]), list(id), unique)[, -1])
## if (!is.null(mPars)) {
## dim(mPars) <- c(N, ninputpars)
## opt$mPars <- mPars
## ix <- rep(1:dim(mPars)[1], nmc)
## optM$mPars <- mPars[ix, ]
## dim(optM$mPars) <- c(nmc * N, ninputpars)
## }
if (is.null(data$nmdat$CMT)) data$nmdat$CMT <- 1 ## CHECKME
if (any(is.na(data$nmdat$CMT))) {
stop("'CMT' has NA(s)")
}
## CHECKME
form <- attr(model$saem_mod, "form")
.nobs <- 0
dat <- RxODE::etTrans(data$nmdat, attr(model$saem_mod, "rx"), addCmt=TRUE, dropUnits=TRUE, allTimeVar=TRUE)
.nobs <- attr(class(dat), ".RxODE.lst")$nobs
## if(length(dat) !=7) stop("SAEM doesn't support time varying covariates yet.");
.rx <- attr(model$saem_mod, "rx")
.pars <- .rx$params
.pars <- setNames(rep(1.1, length(.pars)), .pars)
opt$.rx <- .rx
opt$.pars <- .pars
## opt$.dat <- dat;
dat <- .as.data.frame(dat[, -6])
names(dat) <- toupper(names(dat))
dat$ID <- as.integer(dat$ID)
evt <- dat
evt$ID <- evt$ID - 1
## r
evtM <- evt[rep(1:dim(evt)[1], nmc), ]
evtM$ID <- cumsum(c(FALSE, diff(evtM$ID) != 0))
i1 <- grep(1, diag(covstruct))
i0 <- grep(0, diag(covstruct))
nphi1 <- sum(diag(covstruct))
nphi0 <- nphi - nphi1
na <- length(mcmc$stepsize)
nlambda1 <- sum(mcov[, i1])
nlambda0 <- sum(mcov[, i0])
nlambda <- nlambda1 + nlambda0
nd1 <- nphi1 + nlambda1 + 1
nd2 <- nphi1 + nlambda1 + nlambda0
nb_param <- nd2 + 1
Mcovariables <- cbind(rep(1, N), covariables)[, 1:nrow(mcov)]
dim(Mcovariables) <- c(length(Mcovariables) / nrow(mcov), nrow(mcov)) # FIXME
# get fixed ix
fixed <- inits$theta.fix != ""
wh <- fixed[, i1][mcov[, i1] == 1]
len <- length(wh)
fixed.i1 <- (1:len)[wh] - 1
wh <- fixed[, i0][mcov[, i0] == 1]
len <- length(wh)
fixed.i0 <- (1:len)[wh] - 1
jlog1 <- grep(TRUE, model$log.eta)
jcov <- grep(TRUE, apply(mcov, 1, sum) > 0)
covstruct1 <- covstruct[i1, i1]
dim(covstruct1) <- c(nphi1, nphi1)
ind_cov <- grep(1, mcov[mcov > 0])
mcov1 <- matrix(mcov[, i1], ncol = length(i1))
mcov0 <- matrix(mcov[, i0], nrow = nrow(mcov), ncol = length(i0))
ind_cov1 <- grep(1, mcov1[mcov1 > 0]) - 1
ind_cov0 <- grep(1, mcov0[mcov0 > 0]) - 1
pc <- apply(mcov, 2, sum)
ipc <- cumsum(c(0, pc[1:(nphi - 1)])) + 1
ipcl1 <- ipc[jlog1]
for (x in jlog1) inits$theta[1, x] <- log(inits$theta[1, x])
idx <- as.vector(mcov1 > 0)
COV1 <- Mcovariables[, row(mcov1)[idx]]
dim(COV1) <- c(N, sum(idx))
COV21 <- crossprod(COV1)
x <- mcov1 * col(mcov1)
x <- sapply(x[idx], function(x) {
ret <- rep(0, nphi1)
ret[x] <- 1
ret
})
LCOV1 <- t(x)
dim(LCOV1) <- c(nlambda1, nphi1)
pc1 <- apply(LCOV1, 2, sum)
x1 <- diag(sum(idx))
diag(x1) <- inits$theta[, i1][idx]
MCOV1 <- x1 %*% LCOV1
jcov1 <- grep(1, LCOV1) - 1
idx <- as.vector(mcov0 > 0)
COV0 <- Mcovariables[, row(mcov0)[idx]]
dim(COV0) <- c(N, sum(idx))
COV20 <- crossprod(COV0)
x <- mcov0 * col(mcov0)
x <- sapply(x[idx], function(x) {
ret <- rep(0, nphi0)
ret[x] <- 1
ret
})
LCOV0 <- t(x)
dim(LCOV0) <- c(nlambda0, nphi0)
x1 <- diag(sum(idx))
diag(x1) <- inits$theta[, i0][idx]
if (dim(x1)[1] > 0) {
MCOV0 <- x1 %*% LCOV0
} else {
MCOV0 <- matrix(x1, nrow = 0, ncol = dim(LCOV0)[2])
}
jcov0 <- grep(1, LCOV0) - 1
mprior_phi1 <- Mcovariables %*% inits$theta[, i1]
mprior_phi0 <- Mcovariables %*% inits$theta[, i0]
Gamma2_phi1 <- diag(nphi1)
diag(Gamma2_phi1) <- inits$omega[i1]
Gamma2_phi0 <- diag(nphi0)
diag(Gamma2_phi0) <- inits$omega[i0]
phiM <- matrix(0, N, nphi)
phiM[, i1] <- mprior_phi1
phiM[, i0] <- mprior_phi0
phiM <- phiM[rep(1:N, nmc), , drop = FALSE]
.tmp <- diag(sqrt(inits$omega))
if (model$N.eta == 1) .tmp <- matrix(sqrt(inits$omega))
.mat2 <- matrix(rnorm(phiM), dim(phiM))
phiM <- phiM + .mat2 %*% .tmp
mc.idx <- rep(1:N, nmc)
statphi <- sapply(1:nphi, function(x) tapply(phiM[, x], mc.idx, mean))
statphi11 <- statphi[, i1]
dim(statphi11) <- c(N, length(i1))
statphi01 <- statphi[, i0]
dim(statphi01) <- c(N, length(i0))
statphi12 <- crossprod(phiM[, i1])
statphi02 <- crossprod(phiM[, i0])
# x = mcov *cumsum(mcov)
# x1 = cbind(x[,i1], x[,i0])
# indiphi = order(x1[x1>0]) #FINDME
niter <- sum(mcmc$niter)
niter_phi0 <- niter * .5
nb_sa <- mcmc$niter[1] * .75
nb_correl <- mcmc$niter[1] * .75
va <- mcmc$stepsize
vna <- mcmc$niter
na <- length(va)
pas <- 1 / (1:vna[1])^va[1]
for (ia in 2:na) {
end <- length(pas)
k1 <- pas[end]^(-1 / va[ia])
pas <- c(pas, 1 / ((k1 + 1):(k1 + vna[ia]))^va[ia])
}
pash <- c(rep(1, mcmc$burn.in), 1 / (1:niter))
minv <- rep(1e-20, nphi)
# preserve par order when printing iter history
mcov[mcov == 1] <- 1:nlambda
ilambda1 <- mcov[, i1]
ilambda1 <- ilambda1[ilambda1 > 0] - 1
ilambda0 <- mcov[, i0]
ilambda0 <- ilambda0[ilambda0 > 0] - 1
i1 <- i1 - 1
i0 <- i0 - 1
opt$distribution <- distribution
opt$paramUpdate <- attr(model$saem_mod, "paramUpdate")
optM$paramUpdate <- attr(model$saem_mod, "paramUpdate")
opt$ODEopt <- ODEopt
optM$ODEopt <- ODEopt
cfg <- list(
ODEopt = ODEopt,
inits = inits.save,
nu = mcmc$nu,
niter = niter,
nb_sa = nb_sa,
nb_correl = nb_correl,
niter_phi0 = niter_phi0,
nmc = nmc,
coef_phi0 = .9638, # FIXME
rmcmc = .5,
coef_sa = .95,
pas = pas,
pash = pash,
minv = minv,
N = N,
ntotal = ntotal,
y = y,
yM = yM,
phiM = phiM,
evt = as.matrix(evt),
evtM = as.matrix(evtM),
mlen = mlen,
indioM = indioM,
pc1 = pc1,
covstruct1 = covstruct1,
Mcovariables = Mcovariables,
i1 = i1,
i0 = i0,
nphi1 = nphi1,
nphi0 = nphi0,
nlambda1 = nlambda1,
nlambda0 = nlambda0,
COV0 = COV0,
COV1 = COV1,
COV20 = COV20,
COV21 = COV21,
LCOV0 = LCOV0,
LCOV1 = LCOV1,
MCOV0 = MCOV0,
MCOV1 = MCOV1,
Gamma2_phi0 = Gamma2_phi0,
Gamma2_phi1 = Gamma2_phi1,
mprior_phi0 = mprior_phi0,
mprior_phi1 = mprior_phi1,
jcov0 = jcov0,
jcov1 = jcov1,
ind_cov0 = ind_cov0,
ind_cov1 = ind_cov1,
statphi11 = statphi11,
statphi01 = statphi01,
statphi02 = statphi02,
statphi12 = statphi12,
res.mod = model$res.mod,
ares = inits$ares,
bres = inits$bres,
cres = inits$cres,
lres = inits$lres,
opt = opt,
optM = optM,
print = mcmc$print,
distribution = distribution,
par.hist = matrix(0, sum(niter), nlambda1 + nlambda0 + nphi1 + 1 + (model$res.mod == 2) + 2 * (model$res.mod == 4)),
seed = seed,
fixed.i1 = fixed.i1,
fixed.i0 = fixed.i0,
ilambda1 = as.integer(ilambda1),
ilambda0 = as.integer(ilambda0),
nobs = .nobs)
## CHECKME
s <- cfg$evt[cfg$evt[, "EVID"] == 0, "CMT"]
cfg$opt$cmt_endpnt <- cfg$optM$cmt_endpnt <- sort(unique(s))
cfg$nendpnt <- length(unique(s))
if (model$nendpnt != cfg$nendpnt) {
msg <- sprintf("mis-match in nbr endpoints in model & in data")
stop(msg)
}
t <- unlist(split(1L:length(s), s))
cfg$ysM <- rep(cfg$y[t], cfg$nmc)
cfg$ix_sorting <- t - 1 # c-index for sorting by endpnt
cfg$y_offset <- c(0, cumsum(table(s)))
s <- cfg$evtM[cfg$evtM[, "EVID"] == 0, "CMT"]
cfg$ix_endpnt <- as.integer(as.factor(s)) - 1 # to derive vecares & vecbres
s <- cfg$evtM[cfg$evtM[, "EVID"] == 0, "ID"]
t <- cumsum(c(0, table(s)))
cfg$ix_idM <- cbind(t[-length(t)], t[-1] - 1) # c-index of obs records of each subject
cfg$ares <- rep(10, cfg$nendpnt)
cfg$bres <- rep(1, cfg$nendpnt)
cfg$cres <- rep(1, cfg$nendpnt)
cfg$lres <- rep(1, cfg$nendpnt)
cfg$yj <- rep(2L, cfg$nendpnt)
cfg$propT <- rep(0L, cfg$nendpnt)
cfg$lambda <- rep(1.0, cfg$nendpnt)
cfg$low <- rep(0.0, cfg$nendpnt)
cfg$hi <- rep(1.0, cfg$nendpnt)
cfg$ares[cfg$res.mod == 2] <- 0
cfg$bres[cfg$res.mod == 1] <- 0
nres <- (1:4)[(cfg$res.mod == 10L) * 3 + (cfg$res.mod %in% c(4L, 8L, 9L)) * 2 + (cfg$res.mod %in% c(3L, 5L, 6L, 7L)) + 1]
cfg$res_offset <- cumsum(c(0, nres))
cfg$par.hist <- matrix(0, cfg$niter, nlambda1 + nlambda0 + nphi1 + sum(nres))
cfg$addProp <- c("combined1" = 1L, "combined2" = 2L)[match.arg(addProp)]
cfg$DEBUG <- cfg$opt$DEBUG <- cfg$optM$DEBUG <- DEBUG
cfg$phiMFile <- tempfile("phi-", RxODE::rxTempDir(), ".phi")
cfg$tol <- tol
cfg$itmax <- itmax
cfg$type <- type
cfg$lambdaRange <- lambdaRange
cfg$powRange <- powRange
cfg$odeRecalcFactor <- odeRecalcFactor
cfg$maxOdeRecalc <- maxOdeRecalc
cfg
}
#' Print an SAEM model fit summary
#'
#' Print an SAEM model fit summary
#'
#' @param object a saemFit object
#' @param ... others
#' @return a list
#' @export
summary.saemFit <- function(object, ...) {
fit <- object ## Rcheck hack
th <- fit$Plambda
nth <- length(th)
H <- solve(fit$Ha[1:nth, 1:nth])
se <- sqrt(diag(H))
m <- cbind(exp(th), th, se) # FIXME
## lhsVars = scan("LHS_VARS.txt", what="", quiet=TRUE)
## if (length(lhsVars)==nth) dimnames(m)[[1]] = lhsVars
dimnames(m)[[2]] <- c("th", "log(th)", "se(log_th)")
cat("THETA:\n")
print(m)
cat("\nOMEGA:\n")
print(fit$Gamma2_phi1)
if (any(fit$sig2 == 0)) {
cat("\nSIGMA:\n")
print(max(fit$sig2^2))
} else {
cat("\nARES & BRES:\n")
print(fit$sig2)
}
invisible(list(theta = th, se = se, H = H, omega = fit$Gamma2_phi1, eta = fit$mpost_phi))
}
#' Print an SAEM model fit summary
#'
#' Print an SAEM model fit summary
#'
#' @param x a saemFit object
#' @param ... others
#' @return a list
#' @export
print.saemFit <- function(x, ...) {
fit <- x ## Rcheck hack
th <- fit$Plambda
nth <- length(th)
H <- solve(fit$Ha[1:nth, 1:nth])
se <- sqrt(diag(H))
m <- cbind(exp(th), th, se) # FIXME
## lhsVars = scan("LHS_VARS.txt", what="", quiet=TRUE)
## if (length(lhsVars)==nth) dimnames(m)[[1]] = lhsVars
dimnames(m)[[2]] <- c("th", "log(th)", "se(log_th)")
cat("THETA:\n")
print(m)
cat("\nOMEGA:\n")
print(fit$Gamma2_phi1)
if (any(fit$sig2 == 0)) {
cat("\nSIGMA:\n")
print(max(fit$sig2^2))
} else {
cat("\nARES & BRES:\n")
print(fit$sig2)
}
invisible(list(theta = th, se = se, H = H, omega = fit$Gamma2_phi1, eta = fit$mpost_phi))
}
#' Fit an SAEM model
#'
#' Fit an SAEM model using either closed-form solutions or ODE-based model definitions
#'
#' @param model an RxODE model or lincmt()
#' @param data input data
#' @param inits initial values
#' @param PKpars PKpars function
#' @param pred pred function
#' @param covars Covariates in data
#' @param mcmc a list of various mcmc options
#' @param ODEopt optional ODE solving options
#' @inheritParams configsaem
#' @param seed seed for random number generator
#' @details
#' Fit a generalized nonlinear mixed-effect model using the Stochastic
#' Approximation Expectation-Maximization (SAEM) algorithm
#'
#' @return saem fit object
#'
#' @author Matthew Fidler & Wenping Wang
#' @export
saem.fit <- function(model, data, inits,
PKpars = NULL, pred = NULL,
covars = NULL,
mcmc = list(niter = c(200, 300), nmc = 3, nu = c(2, 2, 2)),
ODEopt = list(atol = 1e-06, rtol = 1e-04, method = "lsoda", transitAbs = FALSE),
distribution = c("normal", "poisson", "binomial", "lnorm"),
seed = 99) {
UseMethod("saem.fit")
}
##' @rdname saem.fit
##' @export
saem <- saem.fit
##' @rdname saem.fit
##' @export
saem.fit.nlmixr.ui.nlme <- function(model, data, inits,
PKpars = NULL, pred = NULL,
covars = NULL,
mcmc = list(niter = c(200, 300), nmc = 3, nu = c(2, 2, 2)),
ODEopt = list(atol = 1e-06, rtol = 1e-04, method = "lsoda", transitAbs = FALSE),
distribution = c("normal", "poisson", "binomial", "lnorm"),
seed = 99) {
call <- as.list(match.call(expand.dots = TRUE))[-1]
names(call)[1] <- "object"
call$est <- "saem"
return(do.call(getFromNamespace("nlmixr", "nlmixr"), call, envir = parent.frame(1)))
}
##' @rdname saem.fit
##' @export
saem.fit.function <- saem.fit.nlmixr.ui.nlme
##' @rdname saem.fit
##' @export
saem.fit.nlmixrUI <- saem.fit.nlmixr.ui.nlme
##' @rdname saem.fit
##' @export
saem.fit.RxODE <- function(model, data, inits,
PKpars = NULL, pred = NULL,
covars = NULL,
mcmc = list(niter = c(200, 300), nmc = 3, nu = c(2, 2, 2)),
ODEopt = list(atol = 1e-06, rtol = 1e-04, method = "lsoda", transitAbs = FALSE),
distribution = c("normal", "poisson", "binomial", "lnorm"),
seed = 99) {
saem_fit <- gen_saem_user_fn(model, PKpars, pred)
model <- list(saem_mod = saem_fit, covars = covars)
cfg <- configsaem(model, data, inits, mcmc, ODEopt, distribution, seed)
fit <- saem_fit(cfg)
## dyn.unload("saem_main.dll")
fit
}
##' @rdname saem.fit
##' @export
saem.fit.default <- function(model, data, inits,
PKpars = NULL, pred = NULL,
covars = NULL,
mcmc = list(niter = c(200, 300), nmc = 3, nu = c(2, 2, 2)),
ODEopt = list(atol = 1e-06, rtol = 1e-04, method = "lsoda", transitAbs = FALSE),
distribution = c("normal", "poisson", "binomial", "lnorm"),
seed = 99) {
saem_fit <- gen_saem_user_fn(model)
model <- list(saem_mod = saem_fit, covars = covars)
cfg <- configsaem(model, data, inits, mcmc, ODEopt, distribution, seed)
fit <- saem_fit(cfg)
# dyn.unload("saem_main.dll")
fit
}
##' @export
ranef.saemFit <- function(object, ...) {
return(object$eta)
}
##' @export
fixef.saemFit <- function(object, ...) {
return(object$Plambda)
}
focei.theta.saemFit <- function(object, uif, ...) {
## Get the thetas needed for FOCEi fit.
this.env <- environment()
if (inherits(uif, "function")) {
uif <- nlmixr(uif)
}
n <- uif$focei.names
thetas <- structure(rep(NA, length(n)), .Names = n)
sf <- structure(as.vector(fixed.effects(object)), .Names = uif$saem.theta.name)
for (n in names(sf)) {
thetas[n] <- sf[n]
}
.predDf <- uif$predDf
.ini <- .as.data.frame(uif$ini)
.resMat <- object$resMat
sapply(seq_along(.predDf$cond), function(i) {
x <- paste(.predDf$cond[i])
.tmp <- .ini[which(.ini$condition == x), ]
.w <- which(sapply(.tmp$err, function(x) any(x == c("prop", "propT", "pow", "powT"))))
if (length(.w) == 1) {
thetas[paste(.tmp$name[.w])] <- .resMat[i, 2]
}
.w <- which(sapply(.tmp$err, function(x) any(x == c("pow2", "powT2"))))
if (length(.w) == 1) {
thetas[paste(.tmp$name[.w])] <- .resMat[i, 3]
}
.w <- which(sapply(seq_along(.tmp$err), function(.x) {
x <- .tmp$err[.x]
if (any(x == c(
"add", "norm", "dnorm", "lnorm", "dlnorm",
"dlogn", "logn"
))) {
if (!is.na(.tmp$est[.x])) {
return(TRUE)
}
}
return(FALSE)
}))
if (length(.w) == 1) {
thetas[paste(.tmp$name[.w])] <- .resMat[i, 1]
}
.w <- which(sapply(.tmp$err, function(x) {
any(x == c("boxCox", "yeoJohnson"))
}))
if (length(.w) == 1) {
thetas[paste(.tmp$name[.w])] <- .resMat[i, 4]
}
assign("thetas", thetas, this.env)
})
return(thetas)
}
focei.eta.saemFit <- function(object, uif, ...) {
if (inherits(uif, "function")) {
uif <- nlmixr(uif)
}
## Reorder based on translation
eta.trans <- uif$saem.eta.trans
for (i in seq(1, max(eta.trans))) {
while (!(any(i == eta.trans)) && max(eta.trans) > i) {
eta.trans[eta.trans >= i] <- eta.trans[eta.trans >= i] - 1
}
}
## orig eta -> new eta
df <- .as.data.frame(uif$ini)
eta <- df[!is.na(df$neta1), ]
len <- length(eta$name)
cur.lhs <- character()
cur.rhs <- numeric()
ome <- character()
cur.ome <- object$Gamma2_phi1
for (i in seq_along(eta$name)) {
last.block <- FALSE
if (i == len) {
last.block <- TRUE
} else if (eta$neta1[i + 1] == eta$neta2[i + 1]) {
last.block <- TRUE
}
if (eta$neta1[i] == eta$neta2[i]) {
cur.lhs <- c(cur.lhs, sprintf("ETA[%d]", eta$neta1[i]))
cur.rhs <- c(cur.rhs, cur.ome[eta.trans[eta$neta1[i]], eta.trans[eta$neta2[i]]])
if (last.block) {
ome[length(ome) + 1] <- sprintf(
"%s ~ %s", paste(cur.lhs, collapse = " + "),
paste(deparse(cur.rhs), collapse = " ")
)
cur.lhs <- character()
cur.rhs <- numeric()
}
} else {
cur.rhs <- c(cur.rhs, cur.ome[eta.trans[eta$neta1[i]], eta.trans[eta$neta2[i]]])
}
}
ome <- eval(parse(text = sprintf("list(%s)", paste(ome, collapse = ","))))
return(ome)
}
as.focei.saemFit <- function(object, uif, pt = proc.time(), ..., data, calcResid = TRUE, obf = NULL,
nnodes.gq = 1, nsd.gq = 3, adjObf = TRUE,
calcCov = TRUE, covMethod = NULL,
calcCovTime = NULL, calcTables = TRUE, keep=NULL, drop=NULL, table=tableControl(),
IDlabel=NULL) {
if (is.null(calcResid)) calcResid <- FALSE
.saemCfg <- attr(object, "saem.cfg")
.saemTime <- proc.time() - pt
if (inherits(uif, "function")) {
uif <- nlmixr(uif)
}
.dist <- ""
if (any(uif$saem.distribution == c("poisson", "binomial"))) {
calcResid <- NA
.dist <- uif$saem.distribution
}
uif.new <- uif
fit <- object
mat <- random.effects(fit)
## Reorder based on translation
eta.trans <- uif$saem.eta.trans
for (i in seq(1, max(eta.trans))) {
while (!(any(i == eta.trans)) && max(eta.trans) > i) {
eta.trans[eta.trans >= i] <- eta.trans[eta.trans >= i] - 1
}
}
mat2 <- mat[, eta.trans, drop = FALSE]
th <- focei.theta(fit, uif)
for (n in names(th)) {
uif.new$est[uif.new$name == n] <- th[n]
}
ome <- focei.eta(fit, uif)
init <- list(
THTA = as.vector(th),
OMGA = ome
)
saem.time <- proc.time() - pt
if (missing(data)) {
stop("Requires Data...")
} else {
dat <- data
}
.tn <- uif$saem.theta.name
.ini <- .as.data.frame(uif$ini)
.ini <- .ini[uif$ini$name %in% .tn, ]
if (any(.ini$fix)) {
.fixed <- paste(.ini$name[.ini$fix])
.tn <- .tn[!(.tn %in% .fixed)]
}
.calcCov <- calcCov
if (uif$env$covMethod == "") {
.cov <- NULL
.addCov <- FALSE
} else if (inherits(calcCov, "matrix")) {
.cov <- calcCov
.addCov <- TRUE
} else {
.nth <- length(.tn)
.ini <- .as.data.frame(uif$ini)
.ini <- .ini[is.na(.ini$err), ]
.ini <- .ini[!is.na(.ini$ntheta), ]
.ini <- .ini[!.ini$fix, ]
.ini <- paste(.ini$name)
.calcCovTime <- proc.time()
if (calcCov) {
.covm <- object$Ha[1:.nth, 1:.nth]
.covm <- try(calc.COV(object))
.doIt <- !inherits(.covm, "try-error")
if (.doIt && dim(.covm)[1] != .nth) .doIt <- FALSE
if (.doIt) {
.tmp <- try(chol(.covm), silent = TRUE)
.addCov <- TRUE
.sqrtm <- FALSE
if (inherits(.tmp, "try-error")) {
.tmp <- .covm
.tmp <- try(sqrtm(.tmp %*% t(.tmp)), silent = FALSE)
if (inherits(.tmp, "try-error")) {
.calcCov <- FALSE
.covm <- object$Ha[1:.nth, 1:.nth]
.tmp <- try(chol(.covm), silent = TRUE)
.addCov <- TRUE
.sqrtm <- FALSE
if (inherits(.tmp, "try-error")) {
.tmp <- object$Ha[1:.nth, 1:.nth]
.tmp <- try(sqrtm(.tmp %*% t(.tmp)), silent = FALSE)
if (inherits(.tmp, "try-error")) {
.addCov <- FALSE
} else {
.sqrtm <- TRUE
}
} else {
.tmp <- object$Ha[1:.nth, 1:.nth]
}
} else {
.sqrtm <- TRUE
}
} else {
.tmp <- .covm
}
} else {
.tmp <- object$Ha[1:.nth, 1:.nth]
.tmp <- try(chol(.covm), silent = TRUE)
.calcCov <- FALSE
.addCov <- TRUE
.sqrtm <- FALSE
if (inherits(.tmp, "try-error")) {
.tmp <- object$Ha[1:.nth, 1:.nth]
.tmp <- try(sqrtm(.tmp %*% t(.tmp)), silent = FALSE)
if (inherits(.tmp, "try-error")) {
.addCov <- FALSE
} else {
.sqrtm <- TRUE
}
} else {
.tmp <- object$Ha[1:.nth, 1:.nth]
.calcCov <- FALSE
}
}
} else {
.tmp <- try(chol(.covm), silent = TRUE)
.addCov <- TRUE
.sqrtm <- FALSE
if (inherits(.tmp, "try-error")) {
.tmp <- object$Ha[1:.nth, 1:.nth]
.tmp <- try(sqrtm(.tmp %*% t(.tmp)), silent = FALSE)
if (inherits(.tmp, "try-error")) {
.addCov <- FALSE
} else {
.sqrtm <- TRUE
}
} else {
.tmp <- object$Ha[1:.nth, 1:.nth]
.calcCov <- FALSE
}
}
if (.addCov) {
if (!.calcCov) {
.cov <- RxODE::rxInv(.tmp)
} else {
.cov <- .tmp
}
attr(.cov, "dimnames") <- list(.tn, .tn)
.cov <- .cov[.ini, .ini, drop = FALSE]
}
.calcCovTime <- proc.time() - .calcCovTime
.calcCovTime <- .calcCovTime["elapsed"]
}
.ini <- .as.data.frame(uif$ini)
.ini <- .ini[!is.na(.ini$ntheta), ]
.skipCov <- !is.na(.ini$err)
.fixed <- uif$focei.fixed
.skipCov <- .skipCov | .fixed[seq_along(.skipCov)]
.covMethod <- uif$env$covMethod
if (!any(.covMethod == c("", "r", "s", "r,s"))) {
.covMethod <- ""
}
if (is.na(calcResid)) .covMethod <- ""
.allThetaNames <- c(uif$saem.theta.name, uif$saem.omega.name, uif$saem.res.name)
.m <- object$par_hist
if (ncol(.m) > length(.allThetaNames)) {
.m <- .m[, seq_along(.allThetaNames)]
}
if (.dist == "binomial") {
.dist <- "bernoulli"
}
dimnames(.m) <- list(NULL, .allThetaNames)
.fixedNames <- paste(uif$ini$name[which(uif$ini$fix)])
.rn <- ""
.likTime <- 0
if (is.na(obf)) {
.saemObf <- NA
} else if (is.null(obf)) {
.likTime <- proc.time()
.saemObf <- calc.2LL(object, nnodes.gq = nnodes.gq, nsd.gq = nsd.gq)
.likTime <- proc.time() - .likTime
.likTime <- .likTime["elapsed"]
if (nnodes.gq == 1) {
.rn <- paste0("laplace", nsd.gq)
} else {
.rn <- paste0("gauss", nnodes.gq, "_", nsd.gq)
}
} else if (is(obf, "logical")) {
if (is.na(obf)) {
.saemObf <- NA
} else if (obf) {
.likTime <- proc.time()
.saemObf <- calc.2LL(object, nnodes.gq = nnodes.gq, nsd.gq = nsd.gq)
.likTime <- proc.time() - .likTime
.likTime <- .likTime["elapsed"]
if (nnodes.gq == 1) {
.rn <- paste0("laplace", nsd.gq)
} else {
.rn <- paste0("gauss", nnodes.gq, "_", nsd.gq)
}
} else {
.saemObf <- NA
}
} else if (is(object, "numeric")) {
.saemObf <- obf
}
.notCalced <- TRUE
.cwresTime <- proc.time()
while (.notCalced) {
.env <- new.env(parent = emptyenv())
.env$table <- table
.env$IDlabel <- IDlabel
.env$nobs2 <- .saemCfg$nobs
.env$nnodes.gq <- nnodes.gq
.env$nsd.gq <- nsd.gq
.env$adjObf <- adjObf
.env$method <- "SAEM"
.env$uif <- uif
.env$saem <- object
if (.addCov) {
.env$cov <- .cov
}
.env$parHistStacked <- .data.frame(
val = as.vector(.m),
par = rep(.allThetaNames, each = nrow(.m)),
iter = rep(1:nrow(.m), ncol(.m))
)
.env$parHist <- .data.frame(iter = rep(1:nrow(.m)), .as.data.frame(.m))
if (length(.fixedNames) > 0) {
.env$parHistStacked <- .env$parHistStacked[!(.env$parHistStacked$par %in% .fixedNames), , drop = FALSE]
.env$parHist <- .env$parHist[, !(names(.env$parHist) %in% .fixedNames), drop = FALSE]
}
if (is.na(calcResid)) {
.setSaemExtra(.env, .rn)
.env$theta <- .data.frame(
lower = -Inf, theta = init$THTA, upper = Inf, fixed = .fixed[seq_along(init$THTA)],
row.names = uif$focei.names
)
.env$fullTheta <- setNames(init$THTA, uif$focei.names)
.om0 <- .genOM(.parseOM(init$OMGA))
attr(.om0, "dimnames") <- list(uif$eta.names, uif$eta.names)
.env$omega <- .om0
.env$etaObf <- .data.frame(ID = seq_along(mat2[, 1]), setNames(.as.data.frame(mat2), uif$eta.names), OBJI = NA)
.env$noLik <- FALSE
.env$objective <- .saemObf
} else if (calcResid) {
.setSaemExtra(.env, .rn)
} else {
.setSaemExtra(.env, .rn)
.env$theta <- .data.frame(
lower = -Inf, theta = init$THTA, upper = Inf, fixed = .fixed[seq_along(init$THTA)],
row.names = uif$focei.names
)
.env$fullTheta <- setNames(init$THTA, uif$focei.names)
.om0 <- .genOM(.parseOM(init$OMGA))
attr(.om0, "dimnames") <- list(uif$eta.names, uif$eta.names)
.env$omega <- .om0
.env$etaObf <- .data.frame(ID = seq_along(mat2[, 1]), setNames(.as.data.frame(mat2), uif$eta.names), OBJI = NA)
.env$noLik <- TRUE
.env$objective <- .saemObf
}
.ctl <- uif$env$ODEopt
names(.ctl) <- sub("maxsteps", "maxstepsOde", names(.ctl))
.ctl <- .ctl[names(.ctl) != "scale"]
.ctl$maxOuterIterations <- 0
.ctl$maxInnerIterations <- 0
.ctl$covMethod <- .covMethod
.ctl$sumProd <- uif$env$sum.prod
.ctl$optExpression <- uif$env$optExpression
.ctl$scaleTo <- 0
.ctl$calcTables <- calcTables
if (.saemCfg$addProp == 1L) {
.ctl$addProp <- "combined1"
} else {
.ctl$addProp <- "combined2"
}
.ctl$method <- "liblsoda"
.ctl <- do.call(foceiControl, .ctl)
if (.ctl$singleOde) {
.pars <- NULL
.mod <- uif$focei.rx1
} else {
.pars <- uif$theta.pars
.mod <- uif$rxode.pred
}
fit.f <- try(foceiFit.data.frame(
data = dat, inits = init, PKpars = .pars,
model = .mod, pred = function() {
return(nlmixr_pred)
}, err = uif$error,
lower = uif$focei.lower,
upper = uif$focei.upper,
thetaNames = uif$focei.names,
etaNames = uif$eta.names,
etaMat = mat2,
env = .env,
fixed = .fixed,
skipCov = .skipCov,
control = .ctl,
keep=keep, drop=drop
), silent = FALSE)
if (inherits(fit.f, "try-error")) {
if (is.na(calcResid)) {
warning("Error calculating nlmixr object, return classic object")
.notCalced <- FALSE
return(object)
} else if (calcResid) {
calcResid <- FALSE
} else {
calcResid <- NA
}
} else {
.notCalced <- FALSE
}
}
.cwresTime <- proc.time() - .cwresTime
if (is.na(calcResid)) {
.cwresTime <- 0
} else if (!calcResid) {
.cwresTime <- 0
}
.env <- fit.f$env
if (uif$env$covMethod == "") {
} else if (inherits(.calcCov, "matrix")) {
if (!is.null(covMethod)) {
.env$covMethod <- covMethod
}
.calcCovTime <- calcCovTime
} else if (.calcCov) {
.env$covMethod <- "linFim"
if (.addCov & .sqrtm) {
.env$covMethod <- "|linFim|"
warning("Covariance matrix non-positive definite, corrected by sqrtm(linFim %*% linFim)")
}
} else {
if (calcCov) {
warning("Linearization of FIM could not be used to calculate covariance.")
}
if (.addCov & .sqrtm) {
.env$covMethod <- "|fim|"
warning("Covariance matrix non-positive definite, corrected by sqrtm(fim %*% fim)")
} else if (!.addCov) {
warning("FIM non-positive definite and cannot be used to calculate the covariance")
}
}
if (is.null(.env$time)) {
.env$time <- .data.frame(saem = .saemTime["elapsed"], check.names = FALSE, row.names = c(""))
} else {
.time <- .env$time
.time <- .time[, !(names(.time) %in% c("optimize", "covariance"))]
.saemTime <- .saemTime["elapsed"]
if (calcResid) {
.saemTime <- .saemTime - .cwresTime["elapsed"]
.time <- .data.frame(.time, cwres = .cwresTime["elapsed"], check.names = FALSE)
}
if (.likTime > 0) {
.time <- .data.frame(.time, logLik = .likTime, check.names = FALSE)
.saemTime <- .saemTime - .likTime
}
if (uif$env$covMethod != "") {
.saemTime <- .saemTime - .calcCovTime
.time <- .data.frame(.time, covariance = .calcCovTime, check.names = FALSE)
}
.env$time <- .data.frame(saem = .saemTime, .time, check.names = FALSE, row.names = c(""))
}
.env$message <- ""
if (is.na(calcResid)) {
row.names(.env$objDf) <- .rn
} else if (calcResid) {
if (!is.na(.saemObf)) {
.llik <- -.saemObf / 2
.nobs <- .env$nobs
attr(.llik, "df") <- attr(get("logLik", .env), "df")
.objf <- ifelse(.env$adjObf, .saemObf - .nobs * log(2 * pi), .saemObf)
## for (.t in c("OBJF","objective", "objf")){
## assign(.t,.objf,.env);
## }
.tmp <- .data.frame(
OBJF = .objf,
AIC = .saemObf + 2 * attr(get("logLik", .env), "df"),
BIC = .saemObf + log(.env$nobs) * attr(get("logLik", .env), "df"),
"Log-likelihood" = as.numeric(.llik), check.names = FALSE
)
if (any(names(.env$objDf) == "Condition Number")) {
.tmp <- .data.frame(.tmp,
"Condition Number" = .env$objDf[, "Condition Number"],
check.names = FALSE
)
}
.env$objDf <- rbind(.env$objDf, .tmp)
row.names(.env$objDf) <- c("FOCEi", .rn)
}
.setSaemExtra(.env, "FOCEi")
} else {
row.names(.env$objDf) <- .rn
}
if (inherits(fit.f, "nlmixrFitData")) {
.cls <- class(fit.f)
.env <- attr(.cls, ".foceiEnv")
.cls[1] <- "nlmixrSaem"
class(fit.f) <- .cls
}
assign("uif", .syncUif(fit.f$uif, fit.f$popDf, fit.f$omega), fit.f$env)
return(fit.f)
}
.saemResidF <- function(x) {
.Call(`_saemResidF`, x)
}
.newuoa <- function(par, fn, gr, lower = -Inf, upper = Inf, control = list(), ...) {
.ctl <- control
if (is.null(.ctl$npt)) .ctl$npt <- length(par) * 2 + 1
if (is.null(.ctl$rhobeg)) .ctl$rhobeg <- 0.2
if (is.null(.ctl$rhoend)) .ctl$rhoend <- 1e-4
.ctl$iprint <- 0L
.ctl <- .ctl[names(.ctl) %in% c("npt", "rhobeg", "rhoend", "iprint", "maxfun")]
.ret <- minqa::newuoa(par, fn,
control = .ctl
)
.ret$x <- .ret$par
.ret$message <- .ret$msg
.ret$convergence <- .ret$ierr
.ret$value <- .ret$fval
return(.ret)
}
## FIXME: coef_phi0, rmcmc, coef_sa
## FIXME: Klog, rho, sa, nmc
## FIXME: N.design
## FIXME: g = gc = 1
## FIXME: ODE inits
## FIXME: Tinf for ODE
## FIXME: chk infusion poor fit
## Local Variables:
## ess-indent-level: 2
## indent-tabs-mode: nil
## End:
Try the nlmixr package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.