Running examples in ‘saemix-Ex.R’ failed
The error most likely occurred in:

> base::assign(".ptime", proc.time(), pos = "CheckExEnv")
> ### Name: default.saemix.plots
> ### Title: Wrapper functions to produce certain sets of default plots
> ### Aliases: default.saemix.plots basic.gof advanced.gof covariate.fits
> ###   individual.fits
> ### Keywords: plots
> 
> ### ** Examples
> 
> 
> data(theo.saemix)
> 
> saemix.data<-saemixData(name.data=theo.saemix,header=TRUE,sep=" ",na=NA, 
+   name.group=c("Id"),name.predictors=c("Dose","Time"),
+   name.response=c("Concentration"),name.covariates=c("Weight","Sex"),
+   units=list(x="hr",y="mg/L",covariates=c("kg","-")), name.X="Time")
Using the object called theo.saemix in this R session as the data.
[1] "Weight" "Sex"   


The following SaemixData object was successfully created:

Object of class SaemixData
    longitudinal data for use with the SAEM algorithm
Dataset theo.saemix 
    Structured data: Concentration ~ Dose + Time | Id 
    X variable for graphs: Time (hr) 
    covariates: Weight (kg), Sex (-) 
      reference class for covariate Sex :  0 
> 
> model1cpt<-function(psi,id,xidep) { 
+ 	  dose<-xidep[,1]
+ 	  tim<-xidep[,2]  
+ 	  ka<-psi[id,1]
+ 	  V<-psi[id,2]
+ 	  CL<-psi[id,3]
+ 	  k<-CL/V
+ 	  ypred<-dose*ka/(V*(ka-k))*(exp(-k*tim)-exp(-ka*tim))
+ 	  return(ypred)
+ }
> 
> saemix.model<-saemixModel(model=model1cpt,
+   description="One-compartment model with first-order absorption", 
+   psi0=matrix(c(1.,20,0.5,0.1,0,-0.01),ncol=3, byrow=TRUE,
+   dimnames=list(NULL, c("ka","V","CL"))),transform.par=c(1,1,1),
+   covariate.model=matrix(c(0,1,0,0,0,0),ncol=3,byrow=TRUE),fixed.estim=c(1,1,1),
+   covariance.model=matrix(c(1,0,0,0,1,0,0,0,1),ncol=3,byrow=TRUE),
+   omega.init=matrix(c(1,0,0,0,1,0,0,0,1),ncol=3,byrow=TRUE),error.model="constant")


The following SaemixModel object was successfully created:

Nonlinear mixed-effects model
  Model function:  One-compartment model with first-order absorption
function (psi, id, xidep) 
{
    dose <- xidep[, 1]
    tim <- xidep[, 2]
    ka <- psi[id, 1]
    V <- psi[id, 2]
    CL <- psi[id, 3]
    k <- CL/V
    ypred <- dose * ka/(V * (ka - k)) * (exp(-k * tim) - exp(-ka * 
        tim))
    return(ypred)
}
  Nb of parameters: 3 
      parameter names:  ka V CL 
      distribution:
     Parameter Distribution Estimated
[1,] ka        log-normal   Estimated
[2,] V         log-normal   Estimated
[3,] CL        log-normal   Estimated
  Variance-covariance matrix:
   ka V CL
ka  1 0  0
V   0 1  0
CL  0 0  1
  Error model: constant , initial values: a.=1 
  Covariate model:
     ka V CL
[1,]  0 1  0
[2,]  0 0  0
    Initial values
              ka  V    CL
Pop.CondInit 1.0 20  0.50
Cov.CondInit 0.1  0 -0.01
> 
> # Reducing the number of iterations due to time constraints for CRAN
> saemix.options<-list(seed=632545,save=FALSE,save.graphs=FALSE,nbiter.saemix=c(100,100))
> 
> saemix.fit<-saemix(saemix.model,saemix.data,saemix.options)
 ----------- FAILURE REPORT -------------- 
 --- failure: the condition has length > 1 ---
 --- srcref --- 
: 
 --- package (from environment) --- 
saemix
 --- call from context --- 
initialiseMainAlgo(saemix.data, saemix.model, saemix.options)
 --- call from argument --- 
if (class(chol.omega) == "try-error") {
    chol.omega <- saemix.model["omega.init"][ind.eta, ind.eta] <- omega[ind.eta, 
        ind.eta] <- mydiag(nrow = length(ind.eta), ncol = length(ind.eta))
    cat("Problem inverting covariance matrix, setting initial Omega to diagonal.\n")
}
 --- R stacktrace ---
where 1: initialiseMainAlgo(saemix.data, saemix.model, saemix.options)
where 2: saemix(saemix.model, saemix.data, saemix.options)

 --- value of length: 2 type: logical ---
[1] FALSE FALSE
 --- function from context --- 
function (saemix.data, saemix.model, saemix.options) 
{
    structural.model <- saemix.model["model"]
    nb.parameters <- saemix.model["nb.parameters"]
    N <- saemix.data["N"]
    pres <- saemix.model["error.init"]
    i0.omega2 <- which((1 - mydiag(saemix.model["covariance.model"])) > 
        0)
    indest.omega <- which(saemix.model["covariance.model"] > 
        0)
    i1.omega2 <- saemix.model@indx.omega
    ind.res <- saemix.model["indx.res"]
    id <- saemix.data["data"][, saemix.data["name.group"]]
    if (length(saemix.data["name.covariates"]) == 0) 
        tab <- data.frame(id = id)
    else tab <- data.frame(id = id, saemix.data["data"][, saemix.data["name.covariates", 
        drop = FALSE]])
    temp2 <- unique(tab)
    temp <- tab[!duplicated(id), , drop = FALSE]
    if (dim(temp)[1] != dim(temp2)[1]) 
        cat("Some covariates have time-varying values; only the first is taken into account in the current version of the algorithm.\n")
    if (length(saemix.data["name.covariates"]) > 0) {
        Mcovariates <- data.frame(id = rep(1, N), temp[, 2:dim(temp)[2]])
    }
    else {
        Mcovariates <- data.frame(id = rep(1, N))
    }
    j.cov <- which(rowSums(saemix.model["betaest.model"]) > 0)
    betaest.model <- saemix.model["betaest.model"][j.cov, , drop = FALSE]
    Mcovariates <- Mcovariates[, j.cov, drop = FALSE]
    for (icol in dim(Mcovariates)[2]) if (is.factor(Mcovariates[, 
        icol])) 
        Mcovariates[, icol] <- as.numeric(Mcovariates[, icol]) - 
            1
    saemix.model["betaest.model"] <- betaest.model
    temp1 <- betaest.model[-c(1), , drop = FALSE]
    saemix.model["covariate.model"] <- temp1
    fixedpsi.ini <- saemix.model["psi0"][1, ]
    betaI.ini <- transpsi(matrix(fixedpsi.ini, nrow = 1), saemix.model["transform.par"])
    fixed.ini <- saemix.model["betaest.model"] * 0
    fixed.ini[1, ] <- betaI.ini
    nr.psi0 <- dim(saemix.model["psi0"])[1]
    nr.cov <- dim(saemix.model["betaest.model"])[1]
    if (nr.psi0 > nr.cov) {
        saemix.model["psi0"] <- saemix.model["psi0"][1:nr.cov, 
            , drop = FALSE]
        nr.psi0 <- dim(saemix.model["psi0"])[1]
    }
    if (nr.psi0 < nr.cov) {
        psi1 <- saemix.model["psi0"][nr.psi0, ]
        for (j in (nr.psi0 + 1):nr.cov) saemix.model["psi0"] <- rbind(saemix.model["psi0"], 
            psi1)
        nr.psi0 <- dim(saemix.model["psi0"])[1]
    }
    if (nr.psi0 > 1) 
        fixed.ini[2:nr.psi0, ] <- saemix.model["psi0"][2:nr.psi0, 
            ]
    covariate.estim <- matrix(rep(saemix.model["fixed.estim"], 
        nr.psi0), byrow = TRUE, ncol = length(saemix.model["fixed.estim"]))
    covariate.estim <- covariate.estim * saemix.model["betaest.model"]
    betas.ini <- fixed.ini[which(saemix.model["betaest.model"] > 
        0)]
    betas.ini <- matrix(betas.ini, ncol = 1)
    nb.betas <- sum(saemix.model["betaest.model"])
    ind.covariate <- which(saemix.model["betaest.model"] == 1)
    LCOV <- MCOV <- matrix(data = 0, nrow = nb.betas, ncol = nb.parameters)
    j1 <- 1
    COV <- matrix(nrow = dim(Mcovariates)[1], ncol = 0)
    pfix <- matrix(data = 0, nrow = 1, ncol = nb.parameters)
    mean.phi <- matrix(data = 0, nrow = N, ncol = nb.parameters)
    for (j in 1:nb.parameters) {
        jcov <- which(saemix.model["betaest.model"][, j] == 1)
        lambdaj <- fixed.ini[jcov, j]
        aj <- as.matrix(Mcovariates[, jcov])
        COV <- cbind(COV, aj)
        nlj <- length(lambdaj)
        j2 <- j1 + nlj - 1
        LCOV[j1:j2, j] <- matrix(data = 1, nrow = nlj, ncol = 1)
        j1 <- j2 + 1
        if (length(jcov) <= 1) 
            mean.phi[, j] <- aj * lambdaj
        else mean.phi[, j] <- aj %*% lambdaj
        pfix[j] <- length(lambdaj)
    }
    indx.betaI <- cumsum(c(0, pfix[1:(nb.parameters - 1)])) + 
        1
    idx <- 1:nb.betas
    indx.betaC <- idx[is.na(match(idx, indx.betaI))]
    saemix.model["indx.fix"] <- indx.betaI
    saemix.model["indx.cov"] <- indx.betaC
    COV2 <- t(COV) %*% COV
    j.covariate <- which(LCOV == 1)
    MCOV[j.covariate] <- betas.ini
    betas <- betas.ini
    ind.fix1 <- which(covariate.estim[ind.covariate] == 1)
    ind.fix0 <- which(covariate.estim[ind.covariate] == 0)
    COV1 <- COV[, ind.fix1]
    dstatCOV <- COV[, ind.fix0, drop = FALSE] %*% MCOV[ind.fix0, 
        ]
    covariate.estim1 <- covariate.estim
    covariate.estim1[, i0.omega2] <- 0
    ind.fix11 <- which(covariate.estim1[ind.covariate] == 1)
    covariate.estim0 <- covariate.estim
    covariate.estim0[, i1.omega2] <- 0
    ind.fix10 <- which(covariate.estim0[ind.covariate] == 1)
    MCOV0 <- MCOV[ind.fix10, i0.omega2, drop = FALSE]
    COV0 <- COV[, ind.fix10]
    j0.covariate <- which(LCOV[ind.fix10, i0.omega2] == 1)
    flag.fmin <- as.integer(sum(covariate.estim0[1, ]) > 0)
    chdat <- new(Class = "SaemixRepData", data = saemix.data, 
        nb.chains = saemix.options$nb.chains)
    NM <- chdat["NM"]
    IdM <- chdat["dataM"]$IdM
    yM <- chdat["dataM"]$yM
    XM <- chdat["dataM"][, c(saemix.data["name.predictors"], 
        saemix.data["name.cens"], saemix.data["name.mdv"], saemix.data["name.ytype"]), 
        drop = FALSE]
    io <- matrix(data = 0, nrow = N, ncol = max(saemix.data["nind.obs"]))
    for (i in 1:N) io[i, 1:saemix.data["nind.obs"][i]] <- 1
    ioM <- do.call(rbind, rep(list(io), saemix.options$nb.chains))
    ind.ioM <- which(t(ioM) != 0)
    DYF <- matrix(data = 0, nrow = dim(ioM)[2], ncol = dim(ioM)[1])
    if (length(i0.omega2) > 0) {
        xmat <- covariate.estim[, i0.omega2]
        if (is.null(dim(xmat))) 
            xmat <- matrix(xmat, ncol = length(i0.omega2))
        i0.temp <- which(colSums(xmat) == 0)
        ind0.eta <- i0.omega2[i0.temp]
    }
    else ind0.eta <- c()
    if (length(ind0.eta) > 0) {
        idx <- 1:nb.parameters
        ind.eta <- idx[-c(ind0.eta)]
    }
    else ind.eta <- 1:nb.parameters
    nb.etas <- length(ind.eta)
    itest.phi <- 1:NM
    ltest.phi <- length(itest.phi)
    phiM <- matrix(data = 0, nrow = NM, ncol = nb.parameters)
    etaM <- matrix(data = 0, nrow = NM, ncol = nb.etas)
    mean.phiM <- do.call(rbind, rep(list(mean.phi), saemix.options$nb.chains))
    kt <- 0
    omega <- saemix.model["omega.init"]
    chol.omega <- try(chol(omega[ind.eta, ind.eta]), silent = TRUE)
    if (class(chol.omega) == "try-error") {
        chol.omega <- saemix.model["omega.init"][ind.eta, ind.eta] <- omega[ind.eta, 
            ind.eta] <- mydiag(nrow = length(ind.eta), ncol = length(ind.eta))
        cat("Problem inverting covariance matrix, setting initial Omega to diagonal.\n")
    }
    phiMc <- mean.phiM
    while (ltest.phi > 0) {
        kt <- kt + 1
        if (kt == 100) 
            stop("stats:fit.saemix:FailedInitialParameterGuess\nFailed to find a valid initial parameter guess\n")
        end
        etaMc <- 0.5 * matrix(rnorm(NM * nb.etas), ncol = nb.etas) %*% 
            chol.omega
        phiMc[, ind.eta] <- mean.phiM[, ind.eta] + etaMc
        etaM[itest.phi, ] <- etaMc[itest.phi, ]
        phiM[itest.phi, ] <- phiMc[itest.phi, ]
        psiM <- transphi(phiM, saemix.model["transform.par"])
        fpred <- structural.model(psiM, IdM, XM)
        inan <- (is.na(fpred) + is.infinite(fpred) + (Im(fpred) != 
            0))
        itest.phi <- unique(IdM[inan])
        ltest.phi <- length(itest.phi)
    }
    var.eta <- mydiag(saemix.model["omega.init"])
    theta0 <- c(fixedpsi.ini, var.eta[i1.omega2], pres[saemix.model["indx.res"]])
    l1 <- betas.ini
    l1[indx.betaI] <- transphi(matrix(l1[indx.betaI], nrow = 1), 
        saemix.model["transform.par"])
    allpar0 <- c(l1, var.eta[i1.omega2], pres[ind.res])
    Dargs <- list(IdM = IdM, XM = XM, yM = yM, NM = NM, N = N, 
        nobs = saemix.data["ntot.obs"], yobs = saemix.data["data"][, 
            saemix.data["name.response"]], transform.par = saemix.model["transform.par"], 
        error.model = saemix.model["error.model"], structural.model = structural.model, 
        etype.exp = which(saemix.model["error.model"] == "exponential"))
    nb.parest <- sum(covariate.estim) + sum(saemix.model["covariance.model"][upper.tri(saemix.model["covariance.model"], 
        diag = TRUE)]) + 1 + as.integer(saemix.model["error.model"] == 
        "combined")
    Uargs <- list(nchains = saemix.options$nb.chains, nb.parameters = nb.parameters, 
        nb.betas = nb.betas, nb.etas = nb.etas, nb.parest = nb.parest, 
        indx.betaC = indx.betaC, indx.betaI = indx.betaI, ind.res = ind.res, 
        indest.omega = indest.omega, i0.omega2 = i0.omega2, i1.omega2 = i1.omega2, 
        j.covariate = j.covariate, j0.covariate = j0.covariate, 
        ind.fix10 = ind.fix10, ind.fix11 = ind.fix11, ind.fix1 = ind.fix1, 
        ind.fix0 = ind.fix0, MCOV0 = MCOV0, COV = COV, COV0 = COV0, 
        COV1 = COV1, LCOV = LCOV, COV2 = COV2, dstatCOV = dstatCOV, 
        Mcovariates = Mcovariates, ind.ioM = ind.ioM)
    omega.eta <- omega[ind.eta, ind.eta]
    varList <- list(pres = pres, ind0.eta = ind0.eta, ind.eta = ind.eta, 
        omega = omega, MCOV = MCOV, domega2 = do.call(cbind, 
            rep(list((sqrt(mydiag(omega.eta))) * saemix.options$rw.ini), 
                nb.etas)), diag.omega = mydiag(omega))
    stepsize <- rep(1, saemix.options$nbiter.tot)
    stepsize[(saemix.options$nbiter.saemix[1] + 1):saemix.options$nbiter.tot] <- 1/(1:saemix.options$nbiter.saemix[2])
    stepsize[1:saemix.options$nbiter.burn] <- 0
    opt <- list(stepsize.rw = saemix.options$stepsize.rw, stepsize = stepsize, 
        proba.mcmc = saemix.options$proba.mcmc, nbiter.mcmc = saemix.options$nbiter.mcmc, 
        nbiter.sa = saemix.options$nbiter.sa, alpha1.sa = saemix.options$alpha.sa, 
        alpha0.sa = 10^(-3/saemix.options$nbiter.sa), nbiter.saemix = saemix.options$nbiter.saemix, 
        maxim.maxiter = saemix.options$maxim.maxiter, flag.fmin = flag.fmin)
    return(list(saemix.model = saemix.model, Dargs = Dargs, Uargs = Uargs, 
        varList = varList, opt = opt, DYF = DYF, phiM = phiM, 
        mean.phi = mean.phi, betas = betas, fixedpsi.ini = fixedpsi.ini, 
        allpar0 = allpar0))
}
<bytecode: 0x559574db4ea8>
<environment: namespace:saemix>
 --- function search by body ---
Function initialiseMainAlgo in namespace saemix has this body.
 ----------- END OF FAILURE REPORT -------------- 
Fatal error: the condition has length > 1