discreteVPC: 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)

Author

Emmanuelle Comets emmanuelle.comets@inserm.fr

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

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

Examples

Run this code

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)

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

Run the code above in your browser using DataLab

Get 50% off unlimited learning