discreteVPC {saemix}

VPC for non Gaussian data models

Description

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

Usage

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

Author(s)

Emmanuelle Comets [e-mail]

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

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)