Nothing
#' Automatic statistical model building
#'
#' Automatic statistical model building is available directly in the lixoftConnectors package using the function
#' \code{runModelBuilding}.
#'
#' buildmlx uses SAMBA (Stochastic Approximation for Model Building Algorithm), an iterative procedure to accelerate and optimize the process of model building by identifying
#' at each step how best to improve some of the model components. This method allows to find
#' the optimal statistical model which minimizes some information criterion in very few steps.
#'
#' Penalization criterion can be either a custom penalization of the form gamma*(number of parameters),
#' AIC (gamma=2) or BIC (gamma=log(N)).
#'
#' Several strategies can be used for building the covariate model at each iteration of the algorithm:
#' \code{direction="full"} means that all the possible models are compared (default when the number of covariates
#' is less than 10). Othrwise, \code{direction} is the mode of stepwise search of \code{stepAIC {MASS}},
#' can be one of "both", "backward", or "forward", with a default of "both" when there are at least 10 covariates.
#lambda='cv', glmnet.settings=NULL)
#' See https://monolix.lixoft.com/rsmlx/ for more details.
#' @param project a string: the initial Monolix project
#' @param final.project a string: the final Monolix project (default adds "_built" to the original project)
#' @param model components of the model to optimize c("residualError", "covariate", "correlation"), (default="all")
#' @param prior list of prior probabilities for each component of the model (default=NULL)
#' @param weight list of penalty weights for each component of the model (default=NULL)
#' @param coef.w1 multiplicative weight coefficient used for the first iteration only (default=0.5)
#' @param paramToUse list of parameters possibly function of covariates (default="all")
#' @param covToTest components of the covariate model that can be modified (default="all")
#' @param covToTransform list of (continuous) covariates to be log-transformed (default="none")
#' @param center.covariate TRUE/FALSE center the covariates of the final model (default=FALSE)
#' @param criterion penalization criterion to optimize c("AIC", "BIC", "BICc", gamma) (default=BICc)
#' @param test TRUE/FALSE perform additional statistical tests for building the model (default=TRUE)
#' @param ll TRUE/FALSE compute the observe likelihood and the criterion to optimize at each iteration
#' @param linearization TRUE/FALSE whether the computation of the likelihood is based on a linearization of the model (default=FALSE)
#' @param fError.min minimum fraction of residual variance for combined error model (default = 1e-3)
#' @param seq.cov TRUE/FALSE whether the covariate model is built before the correlation model
#' @param seq.cov.iter number of iterations before building the correlation model (only when seq.cov=F, default=0)
#' @param seq.corr TRUE/FALSE whether the correlation model is built iteratively (default=TRUE)
#' @param p.max maximum p-value used for removing non significant relationships between covariates and individual parameters (default=0.1)
#' @param p.min vector of 3 minimum p-values used for testing the components of a new model (default=c(0.075, 0.05, 0.1))
#' @param direction method for covariate search c("full", "both", "backward", "forward"), (default="full" or "both")
#' @param steps maximum number of iteration for stepAIC (default=1000)
#' @param n.full maximum number of covariates for an exhaustive comparison of all possible covariate models (default=10)
#' @param max.iter maximum number of iterations (default=20)
#' @param explor.iter number of iterations during the exploratory phase (default=2)
#' @param print TRUE/FALSE display the results (default=TRUE)
#' @param nb.model number of models to display at each iteration (default=1)
#' @return a new Monolix project with a new statistical model.
#' @seealso \code{getModelBuildingSettings} settings for model building with lixoftConnectors \cr
#' \code{runModelBuilding} run model building with lixoftConnectors \cr
#' \code{getModelBuildingResults} results for model building with lixoftConnectors
#' @examples
#' \dontrun{
#' # RsmlxDemo1.mlxtran is a Monolix project for modelling the pharmacokinetics (PK) of warfarin
#' # using a PK model with parameters ka, V, Cl.
#'
#' # By default, buildmlx will compute the best statistical model in term of BIC, i.e ,
#' # the best covariate model, the best correlation model for the three random effects and the best
#' # residual error model in terms of BIC.
#' # In this example, three covariates (wt, age, sex) are available with the data and will be used
#' # for building the covariate model for the three PK parameters:
#' r1 <- buildmlx(project="RsmlxDemo1.mlxtran")
#'
#' # Here, the covariate model will be built for V and Cl only and log-transformation of all
#' # continuous covariates will also be considered:
#' r2 <- buildmlx(project="RsmlxDemo1.mlxtran", paramToUse=c("V", "Cl"), covToTransform="all")
#'
#' # Only the covariate model will be built, using AIC instead of BIC:
#' r3 <- buildmlx(project="RsmlxDemo1.mlxtran", model="covariate", criterion="AIC")
#' }
#'
#' # See http://monolix.lixoft.com/rsmlx/buildmlx/ for detailed examples of use of buildmlx
#' # Download the demo examples here: http://monolix.lixoft.com/rsmlx/installation
#'
#' @importFrom MASS addterm dropterm
#' @importFrom stats coef as.formula model.matrix deviance formula extractAIC factor.scope nobs terms update update.formula
#' @importFrom utils data write.csv packageVersion
#' @importFrom dplyr filter select rename arrange bind_rows mutate
#' @importFrom dplyr %>%
#' @export
buildmlx <- function(project=NULL, final.project=NULL, model="all", prior=NULL, weight=NULL, coef.w1=0.5,
paramToUse="all", covToTest="all", covToTransform="none", center.covariate=FALSE,
criterion="BICc", linearization=FALSE, ll=T, test=T,
direction=NULL, steps=1000, n.full=10,
max.iter=20, explor.iter=2, fError.min=1e-3,
seq.cov=FALSE, seq.cov.iter=0, seq.corr=TRUE,
p.max=0.1, p.min=c(0.075, 0.05, 0.1),
print=TRUE, nb.model=1)
{
message("This functionality is directly available in the lixoftConnectors package with the function runModelBuilding.", call. = FALSE)
ptm <- proc.time()
initRsmlx()
dashed.line <- "--------------------------------------------------\n"
plain.line <- "__________________________________________________\n"
dashed.short <- "-----------------------\n"
plain.short <- "_______________________\n"
op.original <- options()
op.new <- options()
op.new$lixoft_notificationOptions$warnings <- 1 #hide the warning messages
options(op.new)
RsmlxDemo1.project <- RsmlxDemo2.project <- warfarin.data <- resMonolix <- NULL
pi <- 4*atan(1)
if (!is.null(project)) {
r <- prcheck(project, f="build", paramToUse=paramToUse, model=model)
if (r$demo)
return(r$res)
project <- r$project
} else {
project <- mlx.getProjectSettings()$project
}
method.ll <- iop.ll <- pen.coef <- NULL
r <- buildmlx.check(project, final.project, model, paramToUse, covToTest, covToTransform, center.covariate,
criterion, linearization, ll, test, direction, steps, max.iter, explor.iter,
seq.cov, seq.cov.iter, seq.corr, p.max, p.min, print, nb.model, prior, weight, n.full)
if (!is.null(r$change)) return(list(change=F))
for (j in 1:length(r))
eval(parse(text=paste0(names(r)[j],"= r[[j]]")))
r <- def.variable(weight=weight, prior=prior, criterion=criterion)
for (j in 1:length(r))
eval(parse(text=paste0(names(r)[j],"= r[[j]]")))
is.weight <- weight$is.weight
is.prior <- NULL
final.dir <- sub(pattern = "(.*)\\..*$", replacement = "\\1", final.project)
if (dir.exists(final.dir))
unlink(final.dir, recursive=TRUE)
project.dir <- mlx.getProjectSettings()$directory
if (!dir.exists(project.dir))
dir.create(project.dir)
buildmlx.dir <- file.path(mlx.getProjectSettings()$directory,"buildmlx")
Sys.sleep(0.1)
if (dir.exists(buildmlx.dir))
unlink(buildmlx.dir, recursive=TRUE)
Sys.sleep(0.1)
dir.create(buildmlx.dir)
summary.file = file.path(buildmlx.dir,"summary.txt")
Sys.sleep(0.1)
if (!dir.exists(final.dir))
dir.create(final.dir, recursive=T)
#-------------------------------------------------
to.cat <- paste0("\n", dashed.line, "\nBuilding:\n")
if (model$covariate) to.cat <- c(to.cat, " - The covariate model\n")
if (model$correlation) to.cat <- c(to.cat, " - The correlation model\n")
if (model$residualError) to.cat <- c(to.cat, " - The residual error model\n")
to.cat <- c(to.cat, "\n")
print_result(print, summary.file, to.cat=to.cat, to.print=NULL)
print.line <- F
# if (!launched.tasks[["populationParameterEstimation"]]) {
# to.cat <- paste0(plain.line,"\nEstimation of the population parameters using the initial model ... \n")
# print_result(print, summary.file, to.cat=to.cat, to.print=NULL)
# print.line <- T
# mlx.runPopulationParameterEstimation()
# }
#-------------------------------------------------
p.ini <- mlx.getPopulationParameterInformation()
rownames(p.ini) <- p.ini$name
ind.omega <- grep("omega_",p.ini[['name']])
omega <- p.ini$name[ind.omega]
omega.ini <- p.ini[ind.omega,]
error.model <- mlx.getContinuousObservationModel()$errorModel
obs.dist <- mlx.getContinuousObservationModel()$distribution
covariate.model <- mlx.getIndividualParameterModel()$covariateModel
cov.ini <- names(covariate.model[[1]])
correlation.model <- lapply(mlx.getIndividualParameterModel()$correlationBlocks$id, sort)
if (length(correlation.model)==0)
correlation.model <- NULL
error.model.ini <- error.model
covariate.model.ini <- covariate.model
correlation.model.ini <- correlation.model
to.cat <- ("- - - Initialization - - -\n")
if (!print.line)
to.cat <- paste0(plain.line, to.cat)
print_result(print, summary.file, to.cat=to.cat, to.print=NULL)
if (model$covariate) {
to.cat <- "\nCovariate model:\n"
to.print <- formatCovariateModel(covariate.model, cov.ini)
print_result(print, summary.file, to.cat=to.cat, to.print=to.print)
}
if (model$correlation) {
to.cat <- "\nCorrelation model:\n"
to.print <- ifelse(!is.null(correlation.model), correlation.model, "NULL")
print_result(print, summary.file, to.cat=to.cat, to.print=to.print)
}
if (model$residualError) {
to.cat <- "\nResidual error model:\n"
to.print <- formatErrorModel(error.model)
print_result(print, summary.file, to.cat=to.cat, to.print=to.print)
}
if (!mlx.getLaunchedTasks()[["populationParameterEstimation"]]) {
to.cat <- "\nEstimation of the population parameters using the initial model ... \n"
print_result(print, summary.file, to.cat=to.cat, to.print=NULL)
mlx.runPopulationParameterEstimation()
}
if (!mlx.getLaunchedTasks()[["conditionalDistributionSampling"]]) {
to.cat <- "Sampling of the conditional distribution using the initial model ... \n"
print_result(print, summary.file, to.cat=to.cat, to.print=NULL)
mlx.runConditionalDistributionSampling()
}
gi <- mlx.getSimulatedIndividualParameters()
gi <- gi %>% filter(rep==gi$rep[nrow(gi)]) %>% select(-rep)
lin.ll <- method.ll=="linearization"
launched.tasks <- mlx.getLaunchedTasks()
if (iop.ll) {
if (!(method.ll %in% launched.tasks[["logLikelihoodEstimation"]])) {
if (lin.ll & !launched.tasks[["conditionalModeEstimation"]])
mlx.runConditionalModeEstimation()
to.cat <- "Estimation of the log-likelihood of the initial model ... \n"
print_result(print, summary.file, to.cat=to.cat, to.print=NULL)
mlx.runLogLikelihoodEstimation(linearization = lin.ll)
}
ll.ini <- compute.criterion(criterion, method.ll, weight, pen.coef)
ll <- formatLL(mlx.getEstimatedLogLikelihood()[[method.ll]], criterion, ll.ini, is.weight, is.prior)
list.criterion <- ll.ini
to.cat <- paste0("\nEstimated criteria (",method.ll,"):\n")
to.print <- round(ll,2)
print_result(print, summary.file, to.cat=to.cat, to.print=to.print)
# print_result(print, summary.file, to.cat=plain.line, to.print=NULL)
}
mlx.saveProject(final.project)
#-----------------------------------------
if (!is.null(covToTransform)) {
covariate <- mlx.getCovariateInformation()$covariate
for (cov.name in covToTransform) {
covkj <- covariate[[cov.name]]
lckj <- paste0("logt",toupper(substr(cov.name, 1, 1)), substr(cov.name, 2, nchar(cov.name)))
if (!(lckj %in% mlx.getCovariateInformation()$name)) {
tr.str <- paste0(lckj,' = "log(',cov.name,"/",signif(mean(covkj),digits=2),')"')
trs <- paste0("lixoftConnectors::addContinuousTransformedCovariate(",tr.str,")")
eval(parse(text=trs))
# colnames(weight$covariate) <- gsub(cov.name, lckj, colnames(weight$covariate))
foo <- colnames(weight$covariate)
weight$covariate <- cbind(weight$covariate, weight$covariate[,cov.name])
colnames(weight$covariate) <- c(foo, lckj)
}
covToTest <- c(covToTest, lckj)
}
covToTest <- unique(setdiff(covToTest, covToTransform))
covToTransform <- NULL
mlx.saveProject(final.project)
if (!mlx.getLaunchedTasks()[["populationParameterEstimation"]]) {
to.cat <- "\nEstimation of the population parameters using the transformed covariates ... \n"
print_result(print, summary.file, to.cat=to.cat, to.print=NULL)
mlx.runPopulationParameterEstimation()
}
if (!mlx.getLaunchedTasks()[["conditionalDistributionSampling"]]) {
to.cat <- "Sampling of the conditional distribution using the the transformed covariates ... \n"
print_result(print, summary.file, to.cat=to.cat, to.print=NULL)
mlx.runConditionalDistributionSampling()
}
}
#--------------
stop.test <- FALSE
corr.test <- FALSE
iter <- 0
cov.names0 <- cov.names <- NULL
if (identical(covToTest,"all"))
covFix = NULL
else
covFix <- setdiff(mlx.getCovariateInformation()$name, covToTest)
if (iop.ll) {
ll.prev <- Inf
ll.new <- ll.ini
}
sp0 <- NULL
cov.test <- NULL
e <- NULL
while (!stop.test) {
iter <- iter + 1
if (iter==1){
weight.covariate <- weight$covariate*coef.w1
weight.correlation <- weight$correlation*coef.w1
} else {
weight.covariate <- weight$covariate
weight.correlation <- weight$correlation
}
to.cat <- paste0(plain.line,"- - - Iteration ",iter," - - -\n")
print_result(print, summary.file, to.cat=to.cat, to.print=NULL)
obs.dist0 <- obs.dist
error.model0 <- error.model
covariate.model0 <- covariate.model
correlation.model0 <- correlation.model
if (iop.ll)
ll0 <- ll
if (model$residualError) {
res.error <- errorModelSelection(pen.coef=pen.coef[-c(1, 2)], nb.model=nb.model, f.min=fError.min)
if (nb.model==1)
error.model <- res.error
else {
error.model <- mlx.getContinuousObservationModel()$errorModel
for (k in (1:length(error.model)))
error.model[[k]] <- as.character(res.error[[k]]$error.model[1])
}
}
pmax.cov <- p.max
if (model$covariate) {
# pmax.cov <- ifelse(iter <= 1, 1, p.max)
res.covariate <- covariateModelSelection(pen.coef=pen.coef[1], nb.model=nb.model, weight=weight.covariate,
covToTransform=covToTransform, covFix=covFix, direction=direction,
steps=steps, p.max=pmax.cov, paramToUse=paramToUse, sp0=sp0, iter=iter,
correlation.model = correlation.model, n.full=n.full, eta=e)
res.covariate$res <- sortCov(res.covariate$res, cov.ini)
if (iter>explor.iter) sp0 <- res.covariate$sp
covToTransform <- setdiff(covToTransform, res.covariate$tr0)
covariate.model <- res.covariate$model
e <- res.covariate$residuals
if (nb.model==1)
cov.select <- rownames(res.covariate$res)
else
cov.select <- rownames(res.covariate$res[[1]])
cov.names <- lapply(covariate.model[cov.select], function(x) {sort(names(which(x)))})
cov.names0 <- lapply(covariate.model0[cov.select], function(x) {sort(names(which(x)))})
cov.test <- covariate.test(cov.test, covToTest, covToTransform, paramToUse)
} else {
e <- mlx.getSimulatedRandomEffects()
}
if (model$correlation & !corr.test) {
if (isTRUE(all.equal(cov.names0,cov.names))) # & isTRUE(all.equal(error.model0,error.model)))
corr.test <- TRUE
if (!seq.cov & iter>seq.cov.iter)
corr.test <- TRUE
if (corr.test & (seq.cov==T | seq.cov.iter>0)) {
to.cat <- "Start building correlation model too ... \n"
print_result(print, summary.file, to.cat=to.cat, to.print=NULL)
}
}
if (model$correlation & corr.test) {
pen.corr <- ifelse(iter <= 1, pen.coef[1], pen.coef[1])
res.correlation <- correlationModelSelection(e0=e, pen.coef=pen.corr, nb.model=nb.model,
corr0=correlation.model0, seqmod=seq.corr, weight=weight.correlation)
if (nb.model==1)
correlation.model <- res.correlation
else
correlation.model <- res.correlation[[1]]$block
if (length(correlation.model)==0)
correlation.model <- NULL
} else {
res.correlation <- lapply(mlx.getIndividualParameterModel()$correlationBlocks$id, sort)
}
if (length(res.correlation)==0)
res.correlation <- NULL
#-------------------------------
if (max.iter>0 | nb.model>1) {
eq.cov <- isTRUE(all.equal(cov.names0,cov.names)) & (pmax.cov == p.max)
eq.err <- isTRUE(all.equal(error.model0,error.model))
eq.corr <- isTRUE(all.equal(correlation.model0,correlation.model))
eq.dist <- isTRUE(all.equal(obs.dist0,obs.dist))
if (!model$correlation | corr.test) {
if ( eq.cov & eq.err & eq.dist & eq.corr )
stop.test <- TRUE
if ( model$covariate & eq.cov & eq.dist & eq.corr )
stop.test <- TRUE
}
if (stop.test) {
to.cat <- "\nNo difference between two successive iterations\n"
print_result(print, summary.file, to.cat=to.cat, to.print=NULL)
}
if (!stop.test | nb.model>1) {
if (model$covariate) {
to.cat <- "\nCovariate model:\n"
to.print <- res.covariate$res
print_result(print, summary.file, to.cat=to.cat, to.print=to.print)
}
if (model$correlation) {
to.cat <- "\nCorrelation model:\n"
if (!is.null(res.correlation))
to.print <- res.correlation
else
to.print <- "NULL"
print_result(print, summary.file, to.cat=to.cat, to.print=to.print)
}
if (model$residualError) {
to.cat <- "\nResidual error model:\n"
to.print <- res.error
print_result(print, summary.file, to.cat=to.cat, to.print=to.print)
}
}
}
if (!stop.test) {
p.est <- mlx.getEstimatedPopulationParameters()
mlx.setInitialEstimatesToLastEstimates(fixedEffectsOnly = T)
p.ini <- mlx.getPopulationParameterInformation()
rownames(p.ini) <- p.ini$name
i.omega <- which(grepl("omega_", p.ini$name) & (!identical(p.ini$method,"FIXED")))
p.ini$initialValue[i.omega] <- p.est[p.ini$name[i.omega]]*3
# p.ini[omega,] <- omega.ini
jcor <- grep("corr_",p.ini$name)
if (length(jcor)>0) p.ini <- p.ini[-jcor,]
mlx.setPopulationParameterInformation(p.ini)
if (model$residualError) {
emodel <- error.model
odist <- mlx.getContinuousObservationModel()$distribution
for (k in (1:length(emodel))) {
if (identical(emodel[[k]],"exponential")) {
emodel[[k]] <- "constant"
odist[[k]] <- "lognormal"
} else {
odist[[k]] <- "normal"
}
}
mlx.setErrorModel(emodel)
mlx.setObservationDistribution(odist)
}
if (model$covariate) {
if (length(res.covariate$add.covariate) >0) {
for (k in 1:length(res.covariate$add.covariate))
eval(parse(text=res.covariate$add.covariate[[k]]))
}
mlx.setCovariateModel (covariate.model)
}
if (model$correlation & corr.test)
mlx.setCorrelationBlocks(correlation.model)
#-------------------------------
if (max.iter>0) {
if (iop.ll) {
if (ll.new > ll.prev) {
g <- mlx.getGeneralSettings()
g$autochains <- FALSE
g$nbchains <- g$nbchains+1
mlx.setGeneralSettings(g)
}
ll.prev <- ll.new
}
g=mlx.getConditionalDistributionSamplingSettings()
g$nbminiterations <- max(100, g$nbminiterations)
# if (!is.null(g$enablemaxiterations)) {
# g$enablemaxiterations <- T
# g$nbmaxiterations <- 500
# g$nbminiterations <- 150
# g$ratio=0.03
# }
mlx.setConditionalDistributionSamplingSettings(g)
# g <- getPopulationParameterEstimationSettings()
# g$nbexploratoryiterations <- 400
# g$nbsmoothingiterations <- 200
# g$smoothingautostop <- g$exploratoryautostop <- F
# g$simulatedannealing <- T
# setPopulationParameterEstimationSettings(g)
}
mlx.saveProject(final.project)
if (dir.exists(final.dir))
unlink(final.dir, recursive=TRUE)
mlx.loadProject(final.project)
if (max.iter>0) {
to.cat <- paste0("\nRun scenario for model ",iter," ... \nEstimation of the population parameters... \n")
print_result(print, summary.file, to.cat=to.cat, to.print=NULL)
# mlx.runPopulationParameterEstimation(parameter=gi)
mlx.runPopulationParameterEstimation()
to.cat <- "Sampling from the conditional distribution... \n"
print_result(print, summary.file, to.cat=to.cat, to.print=NULL)
mlx.runConditionalDistributionSampling()
gi <- mlx.getSimulatedIndividualParameters()
gi <- gi %>% filter(rep==gi$rep[nrow(gi)]) %>% select(-rep)
if (iop.ll) {
to.cat <- "Estimation of the log-likelihood... \n"
print_result(print, summary.file, to.cat=to.cat, to.print=NULL)
if (lin.ll & !launched.tasks[["conditionalModeEstimation"]])
mlx.runConditionalModeEstimation()
mlx.runLogLikelihoodEstimation(linearization = lin.ll)
ll.new <- compute.criterion(criterion, method.ll, weight, pen.coef)
list.criterion <- c(list.criterion, ll.new)
}
buildmlx.project.iter <- file.path(buildmlx.dir,paste0("iteration",iter,".mlxtran"))
mlx.saveProject(buildmlx.project.iter)
if (iop.ll) {
if (stop.test)
ll <- ll0
else {
ll <- formatLL(mlx.getEstimatedLogLikelihood()[[method.ll]], criterion, ll.new, is.weight)
}
to.cat <- paste0("\nEstimated criteria (",method.ll,"):\n")
to.print <- round(ll,2)
print_result(print, summary.file, to.cat=to.cat, to.print=to.print)
}
if (iter >= max.iter) {
stop.test <- TRUE
to.cat <- "Maximum number of iterations reached\n"
print_result(print, summary.file, to.cat=to.cat, to.print=NULL)
}
}
}
if (max.iter==0)
stop.test <- TRUE
mlx.saveProject(final.project)
}
change.error.model <- NULL
if (iop.ll) {
ll.min <- min(list.criterion)
iter.opt <- which.min(list.criterion)
if (iter.opt==1)
buildmlx.project.iter <- project
else {
buildmlx.project.iter <- file.path(buildmlx.dir,paste0("iteration",iter.opt-1,".mlxtran"))
}
mlx.loadProject(buildmlx.project.iter)
mlx.saveProject(final.project)
} else {
iter.opt <- iter
}
if (test) {
if (!iop.ll) {
g <- as.list(mlx.getLaunchedTasks())$logLikelihoodEstimation
if (!linearization & !("importanceSampling" %in% g))
mlx.runLogLikelihoodEstimation(linearization=F)
if (linearization & !("linearization" %in% g))
mlx.runLogLikelihoodEstimation(linearization=T)
ll.min <- compute.criterion(criterion, method.ll, weight, pen.coef)
}
to.cat <- paste0(plain.line,"- - - Further tests - - -\n")
print_result(print, summary.file, to.cat=to.cat, to.print=NULL)
if (model$covariate & sum(mlx.getIndividualParameterModel()$variability$id)>1) {
g0 <- mlx.getIndividualParameterModel()
covariate <- random.effect <- p.ttest <- p.lrt <- in.model <- p.value <- NULL
r.test <- covariate.test(cov.test, covToTest, covToTransform, paramToUse)
r.test <- r.test %>% filter(!in.model) %>% select(-in.model)
if (is.weight) {
w.cov <- weight$covariate[cbind(r.test[['parameter']], r.test[['covariate']])]
r.test <- r.test %>% mutate(p.value = p.weight(p.value, w.cov, pen.coef[1]))
}
r.cov0 <- res.covariate$r.cov0
for (j in 1:nrow(r.test)) {
pj <- r.test$parameter[j]
if (!is.null(r.cov0[[pj]])) {
if (r.test$covariate[j] %in% r.cov0[[pj]])
r.test$p.value[j] <- 1
}
}
i.min <- which(as.numeric(r.test$p.value) < p.min[1])
if (length(i.min)>0) {
to.cat <- paste0(plain.short,"Add parameters/covariates relationships:\n")
to.print <- r.test[i.min,]
print_result(print, summary.file, to.cat=to.cat, to.print=to.print)
g <- mlx.getIndividualParameterModel()
stop.test <- F
for (i in i.min) {
param.i <- r.test$parameter[i]
cov.i <- r.test$covariate[i]
g$covariateModel[[param.i]][cov.i] <- T
}
mlx.setIndividualParameterModel(g)
iter <- iter+1
to.cat <- paste0("\nRun scenario for model ",iter," ... \nEstimation of the population parameters... \n")
print_result(print, summary.file, to.cat=to.cat, to.print=NULL)
buildmlx.project.iter <- file.path(buildmlx.dir,paste0("iteration",iter,".mlxtran"))
mlx.saveProject(buildmlx.project.iter)
mlx.runPopulationParameterEstimation()
} else {
stop.test <- F
}
if (any(mlx.getObservationInformation()$type != "continuous"))
mlx.runStandardErrorEstimation(linearization=F)
else
mlx.runStandardErrorEstimation(linearization=T)
#mlx.runStandardErrorEstimation(linearization = lin.ll)
# if (is.weight | is.prior) {
# w.cov <- weight$covariate[cbind(r.test[['parameter']], r.test[['covariate']])]
# r.test <- r.test %>% mutate(p.value = p.weight(p.value, w.cov, pen.coef[1]))
# }
g <- mlx.getIndividualParameterModel()
n.param <- g$name
n.cov <- names(g$covariateModel[[1]])
r.test <- mlx.getTests()$wald
pv <- as.numeric(gsub("<", "", r.test$p.value))
pv[which(is.nan(pv))] <- 0.99999
list.ipc <- NULL
for (np in n.param) {
gp <- g$covariateModel[[np]]
ngp <- names(which(gp))
if (length(ngp) > 0) {
for (nc in ngp) {
g$covariateModel[[np]][nc] <- F
ipc <- grep(paste0("beta_",np,"_",nc), r.test$parameter)
pv[ipc] <- p.weight(pv[ipc], weight$covariate[np, nc], pen.coef[1])
if (min(pv[ipc]) < p.min[2])
g$covariateModel[[np]][nc] <- T
else
list.ipc <- c(list.ipc, ipc)
}
}
}
if (identical(g$covariateModel, g0$covariateModel))
stop.test <- T
if (length(list.ipc) >0) {
to.cat <- paste0(plain.short,"Remove parameters/covariates relationships:\n")
method <- statistics <- parameter <- NULL
to.print <- (r.test %>% select(-c(method, statistics)) %>%
rename(coefficient=parameter))[list.ipc,]
print_result(print, summary.file, to.cat=to.cat, to.print=to.print)
}
g1 <- g
if (!stop.test) {
if (length(list.ipc) >0) {
mlx.setIndividualParameterModel(g)
iter <- iter+1
to.cat <- paste0("\nRun scenario for model ",iter," ... \nEstimation of the population parameters... \n")
print_result(print, summary.file, to.cat=to.cat, to.print=NULL)
buildmlx.project.iter <- file.path(buildmlx.dir,paste0("iteration",iter,".mlxtran"))
mlx.saveProject(buildmlx.project.iter)
mlx.runPopulationParameterEstimation()
}
if (lin.ll) {
if(!launched.tasks[["conditionalModeEstimation"]])
mlx.runConditionalModeEstimation()
} else {
to.cat <- "Sampling from the conditional distribution... \n"
print_result(print, summary.file, to.cat=to.cat, to.print=NULL)
mlx.runConditionalDistributionSampling()
}
to.cat <- "Estimation of the log-likelihood... \n"
print_result(print, summary.file, to.cat=to.cat, to.print=NULL)
mlx.runLogLikelihoodEstimation(linearization = lin.ll)
ll.new <- compute.criterion(criterion, method.ll, weight, pen.coef)
ll <- formatLL(mlx.getEstimatedLogLikelihood()[[method.ll]], criterion, ll.new, is.weight, is.prior)
to.cat <- paste0("\nEstimated criteria (",method.ll,"):\n")
to.print <- round(ll,2)
print_result(print, summary.file, to.cat=to.cat, to.print=to.print)
if (ll.new < ll.min) {
ll.min <- ll.new
mlx.saveProject(final.project)
} else {
mlx.loadProject(final.project)
}
}
# --- Wald tests
g <- as.list(mlx.getLaunchedTasks())$standardErrorEstimation
if (!linearization & !("stochasticApproximation" %in% g))
mlx.runStandardErrorEstimation(linearization=F)
if (!("linearization" %in% g))
mlx.runStandardErrorEstimation(linearization=T)
r.test <- mlx.getTests()$wald
g <- mlx.getIndividualParameterModel()
n.param <- g$name
n.cov <- names(g$covariateModel[[1]])
pv <- as.numeric(gsub("<", "", r.test$p.value))
pv[which(is.nan(pv))] <- 0
list.ipc <- NULL
for (np in n.param) {
gp <- g$covariateModel[[np]]
ngp <- names(which(gp))
if (length(ngp) > 0) {
for (nc in ngp) {
ipc <- grep(paste0("beta_",np,"_",nc), r.test$parameter)
pv[ipc] <- p.weight(pv[ipc], weight$covariate[np, nc], pen.coef[1])
if (max(pv[ipc]) > p.min[2]) {
g$covariateModel[[np]][nc] <- F
list.ipc <- c(list.ipc, ipc)
}
}
}
}
if (length(list.ipc) > 0 & !identical(g$covariateModel, g0$covariateModel) & !identical(g$covariateModel, g1$covariateModel)) {
to.cat <- paste0(plain.short,"Remove parameters/covariates relationships:\n")
method <- statistics <- parameter <- NULL
to.print <- (r.test %>% select(-c(method, statistics)) %>%
rename(coefficient=parameter))[list.ipc,]
print_result(print, summary.file, to.cat=to.cat, to.print=to.print)
mlx.setIndividualParameterModel(g)
iter <- iter+1
to.cat <- paste0("\nRun scenario for model ",iter," ... \nEstimation of the population parameters... \n")
print_result(print, summary.file, to.cat=to.cat, to.print=NULL)
buildmlx.project.iter <- file.path(buildmlx.dir,paste0("iteration",iter,".mlxtran"))
mlx.saveProject(buildmlx.project.iter)
mlx.runPopulationParameterEstimation()
if (lin.ll) {
if(!launched.tasks[["conditionalModeEstimation"]])
mlx.runConditionalModeEstimation()
} else {
to.cat <- "Sampling from the conditional distribution... \n"
print_result(print, summary.file, to.cat=to.cat, to.print=NULL)
mlx.runConditionalDistributionSampling()
}
to.cat <- "Estimation of the log-likelihood... \n"
print_result(print, summary.file, to.cat=to.cat, to.print=NULL)
mlx.runLogLikelihoodEstimation(linearization = lin.ll)
ll.new <- compute.criterion(criterion, method.ll, weight, pen.coef)
ll <- formatLL(mlx.getEstimatedLogLikelihood()[[method.ll]], criterion, ll.new, is.weight, is.prior)
to.cat <- paste0("\nEstimated criteria (",method.ll,"):\n")
to.print <- round(ll,2)
print_result(print, summary.file, to.cat=to.cat, to.print=to.print)
if (ll.new < ll.min) {
ll.min <- ll.new
mlx.saveProject(final.project)
}
}
# --- end Wald tests
}
if (model$correlation) {
test.cor <- T
cor.block0 <- mlx.getIndividualParameterModel()$correlationBlocks$id
while (test.cor) {
mlx.loadProject(final.project)
p.cortest <- NULL
if (!mlx.getLaunchedTasks()$conditionalDistributionSampling)
mlx.runConditionalDistributionSampling()
r.test <- correlationTest()$p.value %>% filter(!in.model) %>% rename(p.value=p.cortest)
param1 <- gsub("eta_","",r.test$randomEffect.1)
param2 <- gsub("eta_","",r.test$randomEffect.2)
w.cor <- weight$correlation[cbind(param1, param2)]+weight$correlation[cbind(param2, param1)]
r.test <- r.test %>% mutate(p.value = p.weight(p.value, w.cor, pen.coef[1]))
i.min <- which(as.numeric(r.test$p.value) < p.min[3])
g <- mlx.getIndividualParameterModel()
param.list <- unlist(g$correlationBlocks$id)
if (length(i.min)>0) {
i.min <- i.min[which.min(r.test$p.value[i.min])]
param1 <- gsub("eta_","",r.test$randomEffect.1[i.min])
param2 <- gsub("eta_","",r.test$randomEffect.2[i.min])
test.cor <- F
if ( !(param1 %in% param.list) & !(param2 %in% param.list) ) {
l.block <- length(g$correlationBlocks$id)+1
g$correlationBlocks$id[[l.block]] <- c(param1, param2)
test.cor <- T
}
if ( !(param1 %in% param.list) & (param2 %in% param.list) ) {
l.block <- grep(param2, g$correlationBlocks$id)
g$correlationBlocks$id[[l.block]] <- c(g$correlationBlocks$id[[l.block]], param1)
test.cor <- T
}
if ( (param1 %in% param.list) & !(param2 %in% param.list) ) {
l.block <- grep(param1, g$correlationBlocks$id)
g$correlationBlocks$id[[l.block]] <- c(g$correlationBlocks$id[[l.block]], param2)
test.cor <- T
}
if (test.cor) {
to.cat <- paste0(plain.short, "Add correlation:\n")
to.print <- r.test[i.min,]
print_result(print, summary.file, to.cat=to.cat, to.print=to.print)
to.print <- g$correlationBlocks$id
print_result(print, summary.file, to.cat=NULL, to.print=to.print)
gi <- mlx.getPopulationParameterInformation()
gi$initialValue[which(gi$name==paste0("omega_",param1))] <- 3*gi$initialValue[which(gi$name==paste0("omega_",param1))]
gi$initialValue[which(gi$name==paste0("omega_",param2))] <- 3*gi$initialValue[which(gi$name==paste0("omega_",param2))]
mlx.setPopulationParameterInformation(gi)
mlx.setIndividualParameterModel(g)
###
iter <- iter+1
buildmlx.project.iter <- file.path(buildmlx.dir,paste0("iteration",iter,".mlxtran"))
mlx.saveProject(buildmlx.project.iter)
to.cat <- paste0("\nRun scenario for model ",iter," ... \nEstimation of the population parameters... \n")
print_result(print, summary.file, to.cat=to.cat, to.print=NULL)
mlx.runPopulationParameterEstimation()
if (lin.ll) {
if(!launched.tasks[["conditionalModeEstimation"]])
mlx.runConditionalModeEstimation()
} else {
to.cat <- "Sampling from the conditional distribution... \n"
print_result(print, summary.file, to.cat=to.cat, to.print=NULL)
mlx.runConditionalDistributionSampling()
}
to.cat <- "Estimation of the log-likelihood... \n"
print_result(print, summary.file, to.cat=to.cat, to.print=NULL)
mlx.runLogLikelihoodEstimation(linearization = lin.ll)
ll.new <- compute.criterion(criterion, method.ll, weight, pen.coef)
ll.disp <- formatLL(mlx.getEstimatedLogLikelihood()[[method.ll]], criterion, ll.new, is.weight, is.prior)
to.cat <- paste0("\nEstimated criteria (",method.ll,"):\n")
to.print <- round(ll.disp,2)
print_result(print, summary.file, to.cat=to.cat, to.print=to.print)
if (ll.new < ll.min) {
ll.min <- ll.new
mlx.saveProject(final.project)
} else {
test.cor <- F
}
} else {
test.cor <- F
}
} else {
test.cor <- F
}
}
mlx.loadProject(final.project)
p <- mlx.getEstimatedPopulationParameters()
if (any(mlx.getObservationInformation()$type != "continuous")) {
mlx.runStandardErrorEstimation(linearization=F)
se <- mlx.getEstimatedStandardErrors()$stochasticApproximation$se
names(se) <- mlx.getEstimatedStandardErrors()$stochasticApproximation$parameter
} else {
mlx.runStandardErrorEstimation(linearization=T)
se <- mlx.getEstimatedStandardErrors()$linearization$se
names(se) <- mlx.getEstimatedStandardErrors()$linearization$parameter
}
z <- as.numeric(p)/as.numeric(se[names(p)])
names(z) <- names(p)
pv <- pnorm(-abs(z))*2
pv.corr <- pv[grep("corr_", names(pv))]
if (length(which(pv.corr>p.min[2])) > 0) {
mlx.saveProject(buildmlx.project.iter)
ind <- mlx.getEstimatedIndividualParameters()$saem
pv.block <- strsplit(gsub("corr_", "", names(pv.corr)),"_")
ind.mod <- mlx.getIndividualParameterModel()
cb <- ind.mod$correlationBlocks$id
cl <- unlist(cb)
pv.tot <- rep(0,length(cl))
names(pv.tot) <- cl
for (k in 1: length(pv.corr))
pv.tot[pv.block[[k]]] <- pv.tot[pv.block[[k]]] + log(pv.corr[k])
param0 <- names(which.max(pv.tot))
cb <- lapply(cb, function(x) setdiff(x, param0))
i.cb <- which(lapply(cb, length)==1)
if (length(i.cb)>0)
cb[[i.cb]] <- NULL
if (!identical(cor.block0, cb)) {
ind.mod$correlationBlocks$id <- cb
to.cat <- paste0(plain.short,"Test correlation model:\n")
to.print <- cb
print_result(print, summary.file, to.cat=to.cat, to.print=to.print)
mlx.setInitialEstimatesToLastEstimates(fixedEffectsOnly = TRUE)
mlx.setIndividualParameterModel(ind.mod)
iter <- iter+1
buildmlx.project.iter <- file.path(buildmlx.dir,paste0("iteration",iter,".mlxtran"))
mlx.saveProject(buildmlx.project.iter)
to.cat <- paste0("Run scenario for model ",iter," ... \nEstimation of the population parameters... \n")
print_result(print, summary.file, to.cat=to.cat, to.print=NULL)
mlx.runPopulationParameterEstimation(parameters=ind)
if (lin.ll) {
if(!launched.tasks[["conditionalModeEstimation"]])
mlx.runConditionalModeEstimation()
} else {
to.cat <- "Sampling from the conditional distribution... \n"
print_result(print, summary.file, to.cat=to.cat, to.print=NULL)
mlx.runConditionalDistributionSampling()
}
to.cat <- "Estimation of the log-likelihood... \n"
print_result(print, summary.file, to.cat=to.cat, to.print=NULL)
mlx.runLogLikelihoodEstimation(linearization = lin.ll)
ll.new <- compute.criterion(criterion, method.ll, weight, pen.coef)
ll <- formatLL(mlx.getEstimatedLogLikelihood()[[method.ll]], criterion, ll.new, is.weight, is.prior)
to.cat <- paste0("\nEstimated criteria (",method.ll,"):\n")
to.print <- round(ll,2)
print_result(print, summary.file, to.cat=to.cat, to.print=to.print)
if (ll.new < ll.min) {
ll.min <- ll.new
mlx.saveProject(final.project)
}
}
}
}
}
if (model$covariate & nb.model>1)
res.covariate$res <- sortCov(res.covariate$res[[1]], cov.ini)
mlx.loadProject(final.project)
if (model$covariate)
covariate.model.print <- formatCovariateModel(mlx.getIndividualParameterModel()$covariateModel)
if (model$correlation) {
correlation.model.print <- lapply(mlx.getIndividualParameterModel()$correlationBlocks$id, sort)
if (length(correlation.model.print)==0)
correlation.model.print <- NULL
}
if (model$residualError)
error.model.print <- formatErrorModel(mlx.getContinuousObservationModel()$errorModel)
if (iop.ll) {
ll.final <- compute.criterion(criterion, method.ll, weight, pen.coef)
ll <- formatLL(mlx.getEstimatedLogLikelihood()[[method.ll]], criterion, ll.final, is.weight, is.prior)
}
to.cat <- paste0(plain.line,"\nFinal statistical model:\n")
print_result(print, summary.file, to.cat=to.cat, to.print=NULL)
if (model$covariate) {
to.cat <- "\nCovariate model:\n"
to.print <- covariate.model.print
print_result(print, summary.file, to.cat=to.cat, to.print=to.print)
}
if (model$correlation) {
to.cat <- "\nCorrelation model:\n"
if (!is.null(correlation.model.print))
to.print <- correlation.model.print
else
to.print <- "NULL"
print_result(print, summary.file, to.cat=to.cat, to.print=to.print)
}
if (model$residualError) {
to.cat <- "\nResidual error model:\n"
to.print <- error.model.print
print_result(print, summary.file, to.cat=to.cat, to.print=to.print)
}
if (iop.ll & max.iter>0) {
to.cat <- paste0("\nEstimated criteria (",method.ll,"):\n")
to.print <- round(ll,2)
print_result(print, summary.file, to.cat=to.cat, to.print=to.print)
}
test.del <- FALSE
if (model$covariate & center.covariate) {
foo <- lapply(res.covariate$model,function(x) {which(x)})
cov.model <- unique(unlist(lapply(foo,function(x) {names(x)})))
cov.type <- mlx.getCovariateInformation()$type[cov.model]
cov.cont <- names(cov.type[cov.type=="continuous"])
for (ck in cov.cont) {
cck <- paste0("c",ck)
covk <- mlx.getCovariateInformation()$covariate[[ck]]
tr.str <- paste0(cck,' = "',ck,"-",signif(mean(covk),digits=2),'"')
tr.str <- paste0("lixoftConnectors::addContinuousTransformedCovariate(",tr.str,")")
eval(parse(text=tr.str))
g=mlx.getIndividualParameterModel()$covariateModel
cg <- lapply(g, function(x) {foo <- x[ck]; x[ck]<-x[cck]; x[cck]=foo; return(x)})
mlx.setCovariateModel (cg)
test.del <- TRUE
covariate.model <- cg
}
}
if (test.del & dir.exists(final.dir)) {
unlink(final.dir, recursive=TRUE)
}
g=mlx.getScenario()
g$tasks[[2]]=TRUE
mlx.setScenario(g)
mlx.saveProject(final.project)
mlx.loadProject(final.project)
error.model <- mlx.getContinuousObservationModel()$errorModel
covariate.model <- mlx.getIndividualParameterModel()$covariateModel
correlation.model <- lapply(mlx.getIndividualParameterModel()$correlationBlocks$id, sort)
if (length(correlation.model)==0)
correlation.model <- NULL
dt <- proc.time() - ptm
res <- list(project=final.project, niter=iter.opt, time=dt["elapsed"])
if (model$covariate)
res <- c(res, list(covariate.model=covariate.model))
if (model$correlation)
res <- c(res, list(correlation.model=correlation.model))
if (model$residualError)
res <- c(res, list(error.model=error.model))
to.cat <- paste0("\ntotal time: ", round(dt["elapsed"], digits=1),"s\n", plain.line)
print_result(print, summary.file, to.cat=to.cat, to.print=NULL)
res$change <- !(identical(error.model,error.model.ini) &
identical(covariate.model,covariate.model.ini) &
identical(correlation.model, correlation.model.ini))
res$change.error.model <- change.error.model
res$weight <- weight
res$covToTest <- covToTest
options(op.original)
return(res)
}
buildmlx.check <- function(project, final.project, model, paramToUse, covToTest, covToTransform, center.covariate,
criterion, linearization, ll, test, direction, steps, max.iter, explor.iter,
seq.cov, seq.cov.iter, seq.corr, p.max, p.min, print, nb.model, prior, weight, n.full) {
if (length(mlx.getIndividualParameterModel()$variability)>1)
stop("Multiple levels of variability are not supported in this version of buildmlx", call.=FALSE)
if (is.null(final.project))
final.project <- paste0(sub(pattern = "(.*)\\..*$", replacement = "\\1", project),"_built.mlxtran")
if (!grepl("\\.",final.project))
final.project <- paste0(final.project,".mlxtran")
if (!grepl("\\.mlxtran",final.project))
stop(paste0(final.project, " is not a valid name for a Monolix project (use the .mlxtran extension)"), call.=FALSE)
if (!is.logical(test))
stop(" 'test' should be boolean", call.=FALSE)
if (!is.logical(center.covariate))
stop(" 'center.covariate' should be boolean", call.=FALSE)
if (!is.logical(linearization))
stop(" 'linearization' should be boolean", call.=FALSE)
if (!is.logical(ll))
stop(" 'll' should be boolean", call.=FALSE)
if (!is.logical(print))
stop(" 'print' should be boolean", call.=FALSE)
if (!is.numeric(criterion) && !(criterion %in% c("AIC", "BIC", "BICc")))
stop(" 'criterion' should be in {'AIC', 'BIC', 'BICc'} or be numerical > 0", call.=FALSE)
if (is.numeric(criterion) && criterion<=0)
stop(" 'criterion' should be in {'AIC', 'BIC', 'BICc'} or be numerical > 0", call.=FALSE)
if (!is.numeric(steps) | steps<=0)
stop(" 'steps' should be numerical > 0", call.=FALSE)
if ((round(max.iter)!=max.iter) | max.iter<=0)
stop(" 'max.iter' should be an integer > 0", call.=FALSE)
if ((round(explor.iter)!=explor.iter) | explor.iter<=0)
stop(" 'explor.iter' should be an integer > 0", call.=FALSE)
if ((round(n.full)!=n.full) | n.full<=0)
stop(" 'n.full' should be an integer > 0", call.=FALSE)
if (!is.logical(seq.cov))
stop(" 'seq.cov' should be boolean", call.=FALSE)
if ((round(seq.cov.iter)!=seq.cov.iter) | seq.cov.iter<0)
stop(" 'seq.cov.iter' should be an integer >= 0", call.=FALSE)
if (!is.logical(seq.corr))
stop(" 'seq.corr' should be boolean", call.=FALSE)
if (!is.numeric(p.max) | p.max<=0 | p.max>1)
stop(" 'p.max' should be probability > 0", call.=FALSE)
if (!is.numeric(p.min) | (length(p.min) != 3))
stop(" 'p.min' should be a 3-vector of probabilities", call.=FALSE)
if (min(p.min)<=0 | max(p.min) >1)
stop(" 'p.min' should be a 3-vector of probabilities", call.=FALSE)
if ((round(nb.model)!=nb.model) | nb.model<=0)
stop(" 'nb.model' should be an integer > 0", call.=FALSE)
if (!is.logical(print))
stop(" 'print' should be boolean", call.=FALSE)
cov.info <- mlx.getCovariateInformation()
cov.names <- cov.info$name
cov.types <- cov.info$type
j.trans <- grep("transformed",cov.types)
j.cont <- grep("continuous",cov.types)
cont.cov <- cov.names[j.cont]
# j.strat <- grep("stratification",cov.types)
# strat.cov <- cov.names[j.strat]
if (identical(covToTransform,"all"))
covToTransform = cov.names[j.cont]
else if (identical(covToTransform,"none"))
covToTransform=NULL
if (!is.null(covToTransform)) {
ncat0 <- names(which(cov.types[covToTransform]=="categorical"))
if (length(ncat0)>0) {
warning(paste0(ncat0, " is a categorical covariate and will not be transformed..."), call.=FALSE)
covToTransform <- setdiff(covToTransform, ncat0)
}
ncov0 <- covToTransform[(!(covToTransform %in% cov.names))]
if (length(ncov0)>0) {
warning(paste0(ncov0, " is not a valid covariate and will be ignored"), call.=FALSE)
covToTransform <- setdiff(covToTransform, ncov0)
}
}
if (identical(covToTest,"all"))
covToTest = cov.names
else if (identical(covToTest,"none"))
covToTest=NULL
if (!is.null(covToTest)) {
ncov0 <- covToTest[(!(covToTest %in% cov.names))]
if (length(ncov0)>0) stop(paste0(ncov0, " is not a valid covariate"), call.=FALSE)
}
#covToTest <- setdiff(covToTest, strat.cov)
covToTransform <- intersect(covToTest, covToTransform)
if (length(covToTransform)==0)
covToTransform <- NULL
ind.dist <- mlx.getIndividualParameterModel()$distribution
param.names <- names(ind.dist)
if (identical(paramToUse,"all"))
paramToUse <- param.names
nparam0 <- paramToUse[(!(paramToUse %in% param.names))]
if (length(nparam0)>0)
stop(paste0(nparam0, " is not a valid parameter"), call.=FALSE)
model.names <- c("residualError", "covariate", "correlation")
if (identical(model,"all")) model <- model.names
# if (is.null(cov.names)) model <- setdiff(model, "covariate")
# if (sum(mlx.getIndividualParameterModel()$variability$id) < 2) model <- setdiff(model, "correlation")
mod0 <- model[(!(model %in% c(model.names, "variance")))]
if (length(mod0)>0)
stop(paste0(mod0, " is not a valid model component"), call.=FALSE)
foo <- model
model <- list(F, F, F)
names(model) <- model.names
model[foo] <- T
#------------------------------
if (!any(mlx.getIndividualParameterModel()$variability$id))
stop("\nA least one parameter with random effects is required\n", call.=FALSE)
if (is.null(mlx.getContinuousObservationModel()))
model$residualError <- FALSE
if (is.null(mlx.getCovariateInformation()))
model$covariate <- FALSE
if (sum(mlx.getIndividualParameterModel()$variability$id) < 2)
model$correlation <- FALSE
if (!any(unlist(model))) {
return(list(change=F))
warning("\nThere is no statistical model to build...\n", call.=FALSE)
}
idir <- NULL
if (model$covariate) {
dir.names <- c("full", "both", "backward", "forward")
if (is.null(direction)) {
nbcov <- length(mlx.getCovariateInformation()$name)
direction <- ifelse(nbcov<=n.full,"full","both")
idir <- direction
}
dir0 <- direction[(!(direction %in% dir.names))]
if (length(dir0)>0)
stop(paste0(dir0, " is not a valid direction name"), call.=FALSE)
}
if (ll=="none") {
warning(" Use 'll=FALSE' instead of 'll='none'' ", call.=FALSE)
ll = FALSE
}
if (ll=="importanceSampling") {
warning(" Use 'linearization=FALSE' instead of 'll='importanceSampling'' ", call.=FALSE)
ll = TRUE
linearization=FALSE
}
if (ll=="linearization") {
warning(" Use 'linearization=TRUE' instead of 'll='linearization'' ", call.=FALSE)
ll = TRUE
linearization=TRUE
}
if (linearization)
method.ll<-"linearization"
else
method.ll <- "importanceSampling"
iop.ll <- ll
if (model$covariate) {
p.cov <- prior$covariate
w.cov <- weight$covariate
cov.model <- do.call(rbind, mlx.getIndividualParameterModel()$covariateModel)
if (!is.null(p.cov) & !is.null(w.cov)) {
warning("Covariate model: only 'weight' or 'prior' can be defined, not both. 'weight' will be used and prior will be ignored", call.=FALSE)
p.cov <- NULL
}
#if (!is.null(p.cov) & is.null(w.cov)) {
# if (length(p.cov)==1) {
# foo <- p.cov
# p.cov <- cov.model
# p.cov[is.logical(cov.model)] <- foo
if (length(p.cov) > 1) {
if (!identical(sort(colnames(cov.model)), sort(colnames(p.cov))) | !identical(sort(rownames(cov.model)), sort(rownames(p.cov))))
stop("prior$covariate should be a matrix whose column names are the names of the covariates and whose row names are the names of the parameters", call.=FALSE)
else
prior$covariate <- p.cov[,colnames(cov.model)]
}
# if (is.null(p.cov) & is.null(w.cov))
# w.cov <- 1
# if (!is.null(w.cov)) {
# if (length(w.cov)==1) {
# foo <- w.cov
# w.cov <- cov.model
# w.cov[is.logical(cov.model)] <- foo
# } else {
if (length(w.cov) > 1) {
if (!identical(sort(colnames(cov.model)), sort(colnames(w.cov))) | !identical(sort(rownames(cov.model)), sort(rownames(w.cov))))
stop("weight$covariate should be a matrix whose column names are the names of the covariates and whose row names are the names of the parameters", call.=FALSE)
else
weight$covariate <- w.cov[,colnames(cov.model)]
}
} else {
seq.cov <- F
}
if (model$correlation) {
p.cor <- prior$correlation
w.cor <- weight$correlation
if (!is.null(p.cor) & !is.null(w.cor)) {
warning("Correlation model: only 'weight' or 'prior' can be defined, not both. 'weight' will be used and prior will be ignored", call.=FALSE)
p.cor <- NULL
}
var.model <- mlx.getIndividualParameterModel()$variability$id
n.param <- names(which(var.model))
d.param <- length(n.param)
if (!is.null(p.cor)) {
if (length(p.cor) > 1)
if (!all(n.param %in% colnames(p.cor)) | !all(n.param %in% rownames(p.cor)) |
(!identical(p.cor, t(p.cor)) & !identical(p.cor, lower.tri(p.cor)*p.cor))) {
print(p.cor)
stop("prior$correlation should be a symetrical or triangular inferior square matrix whose column names and row names are the names of the parameters with variability", call.=FALSE)
}
}
if (!is.null(w.cor)) {
if (length(w.cor)==1)
if (length(p.cor) > 1)
if (!all(n.param %in% colnames(w.cor)) | !all(n.param %in% rownames(w.cor)) |
(!identical(w.cor, t(w.cor)) & !identical(w.cor, lower.tri(w.cor)*w.cor))) {
print(w.cor)
stop("weight$correlation should be a symetrical or triangular inferior square matrix whose column names and row names are the names of the parameters with variability", call.=FALSE)
}
}
weight$correlation <- w.cor
prior$correlation <- p.cor
}
return(list(covToTransform=covToTransform, paramToUse=paramToUse, covToTest=covToTest,
final.project=final.project, model=model, direction=direction, idir=idir,
seq.cov=seq.cov, seq.corr=seq.corr, iop.ll=iop.ll, method.ll=method.ll, weight=weight, prior=prior))
}
#------------------------
formatCovariateModel <- function(m, cov.ini=NULL) {
i0 <- which(unlist(lapply(m,function(x) {identical(x,"none")})))
if (length(i0)>0)
m <- m[-i0]
param.names <- names(m)
cov.names <- names(m[[1]])
if (length(cov.names)>0) {
if (is.vector(m[[1]]) | (!is.null(nrow(m[[1]])) && nrow(m[[1]])==1)) {
mr <- matrix(unlist(m),ncol=length(cov.names),byrow=TRUE)*1
row.names(mr) <- param.names
colnames(mr) <- cov.names
if (length(cov.names)==1)
mr <- as.data.frame(mr)
else {
if (length(param.names)>1)
mr <- as.data.frame(mr[,order(colnames(mr))])
else {
mr <- data.frame(mr)
}
}
} else {
mr <- list()
nr <- nrow(m[[1]])
for (k in 1:nr) {
mk <- lapply(m, function(x) x[k,])
mrk <- matrix(unlist(mk),ncol=length(cov.names),byrow=TRUE)*1
row.names(mrk) <- param.names
colnames(mrk) <- cov.names
mrk <- mrk[,order(colnames(mrk))]
mr[[k]] <- as.data.frame(mrk)
}
}
if (!is.null(cov.ini))
mr <- sortCov(mr, cov.ini)
return(mr)
}
}
sortCov <- function(r, cov.ini) {
n1 <- cov.ini
n3 <- c("ll", "df", "criterion")
if (is.data.frame(r)) {
n2 <- setdiff(names(r), c(n1,n3))
rs <- r[,c(n1,n2)]
} else {
n2 <- setdiff(names(r[[1]]), c(n1,n3))
rs <- lapply(r, function(x) x[c(n1,n2,n3)])
}
return(rs)
}
formatErrorModel <- function(m) {
i.ln <- which(mlx.getContinuousObservationModel()$distribution=="logNormal")
if (length(i.ln) > 0)
m[i.ln] <- "exponential"
return(m)
}
formatLL <- function(ll, criterion, cr, is.weight, is.prior=F) {
is.prior <- F
llr <- ll[c('AIC', 'BIC', 'BICc')]
if (is.prior) {
llr['criterion'] <- cr
} else {
if (is.numeric(criterion)) {
if (is.weight)
llr['w.criterion'] <- cr
else
llr['criterion'] <- cr
} else if (is.weight){
llr[paste0("w",criterion)] <- cr
}
}
if (!is.null(ll) && !is.na(ll["standardError"]))
llr["s.e."] <- ll["standardError"]
return(llr)
}
compute.criterion <- function(criterion, method.ll, weight=NULL, pen.coef=NULL) {
ofv <- mlx.getEstimatedLogLikelihood()[[method.ll]][["OFV"]]
ind.model <- mlx.getIndividualParameterModel()
cov.model <- do.call(rbind, ind.model$covariateModel)*1
if (ncol(cov.model) > 0) {
cov.info <- mlx.getCovariateInformation()
cat.name <- cov.info$name[grep("categorical",cov.info$type)]
if (length(cat.name)>0) {
cat.cov <- cov.info$covariate[cat.name]
cat.cov[] <- lapply( cat.cov, factor)
cat.nb <- unlist(lapply(cat.cov, nlevels)) -1
for (j in cat.name)
cov.model[, j] <- cov.model[, j]*cat.nb[j]
}
i1 <- names(which(ind.model$variability$id))
i0 <- names(which(!ind.model$variability$id))
pen.covariate <- sum((cov.model*weight$covariate)[i1,])*pen.coef[1] + sum((cov.model*weight$covariate)[i0,])*pen.coef[2]
} else {
pen.covariate <- 0
}
cB <- ind.model$correlationBlocks$id
cor.model <- weight$correlation*0
for (k in seq_along(cB))
cor.model[cB[[k]], cB[[k]]] <- 1
pen.correlation <- sum((lower.tri(cor.model)*cor.model)*weight$correlation)*pen.coef[1]
v <- ind.model$variability$id
pen.variance <- sum(v*weight$variance[names(v)])*pen.coef[1]
pen.pop <- length(v)*pen.coef[2]
error.model <- mlx.getContinuousObservationModel()$errorModel
pen.error <- 0
for (k in seq_along(error.model))
pen.error <- pen.error + pen.coef[2+k]*(1 + grepl("combined",error.model[[k]]))
cr <- ofv + pen.pop + pen.covariate + pen.correlation + pen.variance + pen.error
return(cr)
}
compute.criterion.old <- function(criterion, method.ll, weight=NULL, pen.coef=NULL) {
ll <- mlx.getEstimatedLogLikelihood()[[method.ll]]
if (identical(criterion,"AIC")) cr <- ll[["AIC"]]
else if (identical(criterion,"BIC")) cr <- ll[["BIC"]]
else if (identical(criterion,"BICc")) cr <- ll[["BICc"]]
else {
d <- (ll[["AIC"]] - ll[["OFV"]])/2
cr <- ll[["OFV"]]+d*criterion
}
cov.model <- do.call(rbind, mlx.getIndividualParameterModel()$covariateModel)
w.covariate <- sum(cov.model*(weight$covariate - 1))*pen.coef[1]
cB <- mlx.getIndividualParameterModel()$correlationBlocks$id
cor.model <- weight$correlation*0
for (k in seq_along(cB))
cor.model[cB[[k]], cB[[k]]] <- 1
w.correlation <- sum((lower.tri(cor.model)*cor.model)*(weight$correlation-1))*pen.coef[1]
v <- mlx.getIndividualParameterModel()$variability$id
w.variance <- sum(v*(weight$variance[names(v)]-1))*pen.coef[1]
return(cr + w.covariate + w.correlation + w.variance)
}
print_result <- function(print, summary.file, to.cat=NULL,to.print=NULL) {
if (file.exists(summary.file))
sink(summary.file, append=TRUE)
else
sink(summary.file)
if (!is.null(to.cat))
cat(to.cat)
if (!is.null(to.print))
print(to.print)
sink()
if (print) {
if (!is.null(to.cat))
cat(to.cat)
if (!is.null(to.print))
print(to.print)
}
}
covariate.test <- function(cov.test, covToTest, covToTransform, paramToUse) {
p.ttest <- random.effect <- covariate <- parameter <- id <- NULL
r.test <- covariateTest()$p.value.randomEffects %>%
rename(p.value=p.ttest) %>%
mutate(parameter=gsub("eta_","", random.effect)) %>%
select(-random.effect)
r.test <- r.test %>%
filter(covariate %in% covToTest & parameter %in% paramToUse)
if (!is.null(cov.test))
r.test$p.value <- pmin(cov.test$p.value, r.test$p.value)
r.test <- r.test[,c("parameter", "covariate", "p.value", "in.model")]
return(r.test)
}
p.weight <- function(p, pw, coef) {
p <- pmax(p, 2.2e-16)
p <- pmin(p, 0.99999)
A <- pmax(p/(1-p)*(exp(pw*coef/2)-1)/(exp(coef/2)-1) , 0)
p <- A/(1+A)
# r <- exp(coef*(pw-1)/2)
# p <- p*r/(1 - p + p*r)
return(p)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.