compare.saemix | R Documentation |

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).

compare.saemix(..., method = c("is", "lin", "gq"))

`...` |
Two or more objects returned by the |

`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. |

Note that the comparison between two or more models will only be valid if they are fitted to the same dataset.

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.

Maud Delattre

M Delattre, M Lavielle, MA Poursat (2014) A note on BIC in mixed effects models. Electronic Journal of Statistics 8(1) p. 456-475

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") # 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)) # 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")

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.