# R/func_compare.R In saemix: Stochastic Approximation Expectation Maximization (SAEM) Algorithm

#### Documented in compare.saemix

```###########################  Model Comparison with BIC or AIC 	#############################

#' Model comparison with information criteria (AIC, BIC).
#'
#' A specific penalty is used for BIC (BIC.cov) when the compared models have in common the
#' structural model and the covariance structure for the random effects (see Delattre et al., 2014).
#'
#' Note that the comparison between two or more models will only be valid if they are
#' fitted to the same dataset.
#'
#' @param \dots Two or more objects returned by the \code{\link{saemix}} function
#' @param method The method used for computing the likelihood : "is" (Importance Sampling),
#' "lin" (Linearisation) or "gq" (Gaussian quadrature). The default is to use importance sampling "is".
#' If the requested likelihood is not available in all model objects, the method stops with a warning.
#' @return A matrix of information criteria is returned, with at least two columns containing respectively
#' AIC and BIC values for each of the compared models. When the models have in common the structural model
#' and the covariance structure for the random effects, the matrix includes an additional column with BIC.cov
#' values that are more appropriate when the comparison only concerns the covariates.
#' @author Maud Delattre
#' @references M Delattre, M Lavielle, MA Poursat (2014) A note on BIC in mixed effects models.
#' Electronic Journal of Statistics 8(1) p. 456-475
#' @keywords model AIC BIC
#' @examples
#' data(theo.saemix)
#'
#'   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")
#'
#' # Definition of models to be compared
#' 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)
#' }
#' # Model with one covariate
#' saemix.model1<-saemixModel(model=model1cpt,modeltype="structural",
#'   description="One-compartment model, clearance dependent on weight",
#'   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,0,1,0,0,0),ncol=3,byrow=TRUE))
#' # Model with two covariates
#' saemix.model2<-saemixModel(model=model1cpt,modeltype="structural",
#'    description="One-compartment model, clearance dependent on weight, volume dependent on sex",
#'    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,0,1,0,1,0),ncol=3,byrow=TRUE))
#' # Model with three covariates
#' saemix.model3<-saemixModel(model=model1cpt,modeltype="structural",
#'   description="One-cpt model, clearance, absorption dependent on weight, volume dependent on sex",
#'   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(1,0,1,0,1,0),ncol=3,byrow=TRUE))
#'
#' \donttest{
#' # Running the main algorithm to estimate the population parameters
#' saemix.options<-list(seed=632545,save=FALSE,save.graphs=FALSE)
#' saemix.fit1<-saemix(saemix.model1,saemix.data,saemix.options)
#' saemix.fit2<-saemix(saemix.model2,saemix.data,saemix.options)
#' saemix.fit3<-saemix(saemix.model3,saemix.data,saemix.options)
#'
#' compare.saemix(saemix.fit1, saemix.fit2, saemix.fit3)
#' compare.saemix(saemix.fit1, saemix.fit2, saemix.fit3, method = "lin")
#'
#' # We need to explicitly run Gaussian Quadrature if we want to use it in
#' # the comparisons
#' saemix.fit1 <- llgq.saemix(saemix.fit1)
#' saemix.fit2 <- llgq.saemix(saemix.fit2)
#' saemix.fit3 <- llgq.saemix(saemix.fit3)
#' compare.saemix(saemix.fit1, saemix.fit2, saemix.fit3, method = "gq")
#'
#' }
#' @export

compare.saemix<-function(..., method = c("is", "lin", "gq")) {

mod.list <- list(...)
nb.mod <- length(mod.list)
if(nb.mod==1 & typeof(mod.list[])=="list") { # Models given as a list, sorting them out
mod.list1<-list()
x<-mod.list[]
i1<-1
for(i in 1:length(x)) {
if(is(x[[i]], "SaemixObject")) {
mod.list1[[i1]]<-x[[i]]
i1<-i1+1
}
}
mod.list<-mod.list1
nb.mod <- length(mod.list)
}

if (nb.mod < 2){
warning("'compare.saemix' requires at least two models.")
return()
}

list.class <- sapply(mod.list, class)
if (!all(list.class=="SaemixObject")) {
warning("All inputs should have class 'SaemixObject'.")
return()
}

method <- match.arg(method)

# This returns TRUE if the requested likelihood is available
ll_available <- function(x, method) {
length(slot(x@results, paste0("ll.", method))) > 0
}

for (k in 1:nb.mod){
x <- mod.list[[k]]
if (!ll_available( x, method)) {
warning("The likelihood calculated by 'll", method, ".saemix' is not available for model ", k)
return()
}
if(x@model@modeltype=="likelihood" & method=="lin") warning("Linearisation is not appropriate for computing likelihoods in discrete models.")
}

aic_slot <- paste0("aic.", method)
bic_slot <- paste0("bic.", method)
bic_cov_slot <- paste0("bic.covariate.", method)

info <- matrix(0,nb.mod,2)

rownames(info) <- 1:nb.mod
colnames(info) <- c("AIC","BIC")

str.model <- list()
cov.model <- list()

# Comparison of datasets : if models are estimated on different datasets, comparison does
# not make sense.

same.data <- rep(NA,(nb.mod-1))

for (k in 2:nb.mod){
same.data[k-1] <- identical(mod.list[]@data,mod.list[[k]]@data)
}

if ("FALSE" %in% same.data){
warning('Compared models should be fitted on the same data.')
return()
}

# Comparison of the statistical models : BIC.covariate only makes sense if model type,
# structural model (and residual model if appropriate), and covariance structure for the random
# effects are the same

for (k  in 1:nb.mod){
info[k,] <- c(
slot(mod.list[[k]]@results, aic_slot),
slot(mod.list[[k]]@results, bic_slot))
}

same.model.type <- rep(NA,(nb.mod-1))
same.str.model <- rep(NA,(nb.mod-1))

for (k in 2:nb.mod){
same.model.type[k-1] <- identical(mod.list[]@model@modeltype,mod.list[[k]]@model@modeltype)
same.str.model[k-1] <- identical(mod.list[]@model@model,mod.list[[k]]@model@model)
}

if (!(FALSE %in% same.model.type)){

if (!(FALSE %in% same.str.model)){

same.cov.model <- rep(NA,(nb.mod-1))

for (k in 2:nb.mod){
same.cov.model[k-1] <- identical(mod.list[]@model@covariance.model,mod.list[[k]]@model@covariance.model)
}

if (mod.list[]@model@modeltype == "structural"){

same.res.model <- rep(NA,(nb.mod-1))

for (k in 2:nb.mod){
same.res.model[k-1] <- identical(mod.list[]@model@error.model,mod.list[[k]]@model@error.model)
}

if ((!(FALSE %in% same.cov.model)) && (!(FALSE %in% same.res.model))){
bic.cov <- rep(NA,nb.mod)
for (k in 1:nb.mod){
bic.cov[k] <- slot(mod.list[[k]]@results, bic_cov_slot)
}
info <- cbind(info, bic.cov)
colnames(info) <- c("AIC","BIC","BIC.cov")

}

} else {
if (!(FALSE %in% same.cov.model)){
bic.cov <- rep(NA,nb.mod)
for (k in 1:nb.mod){
bic.cov[k] <- slot(mod.list[[k]]@results, bic_cov_slot)
}
info <- cbind(info, bic.cov)
colnames(info) <- c("AIC","BIC","BIC.cov")
}
}
}
}

message("Likelihoods calculated by ",
switch(method,
is = "importance sampling",
lin = "linearisation",
gq = "Gaussian quadrature"))
return(info)
}
```

## Try the saemix package in your browser

Any scripts or data that you put into this service are public.

saemix documentation built on Aug. 5, 2022, 5:25 p.m.