discreteVPC function

VPC for non Gaussian data models

This function provides VPC plots for non Gaussian data models (work in progress)

discreteVPC(object, outcome = "categorical", verbose = FALSE, ...)

Arguments

• object: an saemixObject object returned by the saemix function. The object must include simulated data under the empirical design, using the model and estimated parameters from a fit, produced via the simulateDiscreteSaemix function.
• outcome: type of outcome (valid types are "TTE", "binary", "categorical", "count")
• verbose: whether to print messages (defaults to FALSE)
• ...: additional arguments, used to pass graphical options (to be implemented, currently not available)

Details

This function is a very rough first attempt at automatically creating VPC plots for models defined through their log-likelihood (categorical, count, or TTE data). It makes use of the new element simulate.function in the model component of the object

• for TTE data, a KM-VPC plot will be produced
• for count, categorical and binary data, a plot showing the proportion of each score/category across time will be shown along with the corresponding prediction intervals from the model These plots can be stratified over a covariate in the data set (currently only categorical covariates) by passing an argument which.cov='name' to the call

Examples

data(lung.saemix)

saemix.data<-saemixData(name.data=lung.saemix,header=TRUE,name.group=c("id"),
name.predictors=c("time","status","cens"),name.response=c("status"),
name.covariates=c("age", "sex", "ph.ecog", "ph.karno", "pat.karno", "wt.loss","meal.cal"),
units=list(x="days",y="",covariates=c("yr","","-","%","%","cal","pounds")))

weibulltte.model<-function(psi,id,xidep) {
  T<-xidep[,1]
  y<-xidep[,2] # events (1=event, 0=no event)
  cens<-which(xidep[,3]==1) # censoring times (subject specific)
  init <- which(T==0)
  lambda <- psi[id,1] # Parameters of the Weibull model
  beta <- psi[id,2]
  Nj <- length(T)
  ind <- setdiff(1:Nj, append(init,cens)) # indices of events
  hazard <- (beta/lambda)*(T/lambda)^(beta-1) # H'
  H <- (T/lambda)^beta # H
  logpdf <- rep(0,Nj) # ln(l(T=0))=0
  logpdf[cens] <- -H[cens] + H[cens-1] # ln(l(T=censoring time))
  logpdf[ind] <- -H[ind] + H[ind-1] + log(hazard[ind]) # ln(l(T=event time))
  return(logpdf)
}

simulateWeibullTTE <- function(psi,id,xidep) {
  T<-xidep[,1]
  y<-xidep[,2] # events (1=event, 0=no event)
  cens<-which(xidep[,3]==1) # censoring times (subject specific)
  init <- which(T==0)
  lambda <- psi[,1] # Parameters of the Weibull model
  beta <- psi[,2]
  tevent<-T
  Vj<-runif(dim(psi)[1])
  tsim<-lambda*(-log(Vj))^(1/beta) # nsuj events
  tevent[T>0]<-tsim
  tevent[tevent[cens]>T[cens]] <- T[tevent[cens]>T[cens]]
  return(tevent)
  }
saemix.model<-saemixModel(model=weibulltte.model,description="time model",modeltype="likelihood",
       simulate.function = simulateWeibullTTE,
                psi0=matrix(c(1,2),ncol=2,byrow=TRUE,dimnames=list(NULL,  c("lambda","beta"))),
                transform.par=c(1,1),covariance.model=matrix(c(1,0,0,0),ncol=2, byrow=TRUE))
saemix.options<-list(seed=632545,save=FALSE,save.graphs=FALSE, displayProgress=FALSE)

tte.fit<-saemix(saemix.model,saemix.data,saemix.options)
simtte.fit <- simulateDiscreteSaemix(tte.fit, nsim=100)
gpl <- discreteVPC(simtte.fit, outcome="TTE")
plot(gpl)

References

Brendel, K, Comets, E, Laffont, C, Laveille, C, Mentre, F. Metrics for external model evaluation with an application to the population pharmacokinetics of gliclazide, Pharmaceutical Research 23 (2006), 2036-2049.

Holford, N. The Visual Predictive Check: superiority to standard diagnostic (Rorschach) plots (Abstract 738), in: 14th Meeting of the Population Approach Group in Europe, Pamplona, Spain, 2005.

Ron Keizer, tutorials on VPC TODO

See Also

SaemixObject, saemix, saemix.plot.vpc, simulateDiscreteSaemix

Author(s)

Emmanuelle Comets emmanuelle.comets@inserm.fr