LDPDdoublyint {DPpackage} | R Documentation |
This function generates a posterior density sample for a linear dependent Poisson Dirichlet process mixture model for the analysis of doubly-invertval-censored survival data.
LDPDdoublyint(onset,failure,p,xonset,q,xfailure,xpred, grid,prior,mcmc,state,status,work.dir=NULL)
onset |
a matrix giving the interval limits for the onset times. For multiple variables the limits must be included in a sequential manner for each variable, i.e., (L1,U1,L2,U2,...). |
failure |
a matrix giving the interval limits for the time-to-event times. For multiple variables the limits must be included in a sequential manner for each variable, i.e. (L1,U1,L2,U2,...). |
p |
an integer giving the number of predictors included for each onset variable. Different design matrices are allowed for the different responses but of the same p-dimension. |
xonset |
a matrix giving the design matrices for each onset time. For multiple variables the design matrices must be included in order and includes the intercepts, i.e. (XO1,XO2,...). |
q |
an integer giving the number of predictors included for each time-to-event variable. Different design matrices are allowed for the different responses but of the same q-dimension. |
xfailure |
a matrix giving the design matrices for each time-to-event variable. For multiple variables the design matrices must be included in order and includes the intercepts, i.e. (XT1,XT2,...). |
xpred |
a matrix giving the value of the predictors for which survival and
hazard functions will be evaluated. The number of columns of xpred should be
(p+q)*nvar/2 where nvar is the number of variables. The design matrices
for the predictions must include onset predictors first and then time-to-event
predictors, i.e., (XO1,XO2,...,XT1,XT2,...). |
grid |
a matrix including the grids where survival and hazard functions are evaluated. Each row must include the grid points for different variable. |
prior |
a list giving the prior information. The list includes the following
parameter: q , a0 and b0 giving the hyperparameters for
prior distribution of the a precision parameter of the Poisson-Dirichlet process
prior, mub and sigmab giving the hyperparameters for the prior
distributions of the b precision parameter of the Poisson-Dirichlet process
prior, nu and tinv giving the hyperparameters of the
inverted Wishart prior distribution for the kernel covariance matrix,
mb and Sb giving the hyperparameters of the normal prior distribution
for the mean of the normal baseline distribution,
nub and tbinv giving the hyperparameters of the
inverted Wishart prior distribution for the for the covariance matrix
of the normal baseline distribution, and maxm giving the finite truncation
limit of the sitck-breaking representation of the Poisson-Dirichlet process. |
mcmc |
a list giving the MCMC parameters. The list must include
the following integers: nburn giving the number of burn-in
scans, nskip giving the thinning interval, nsave giving
the total number of scans to be saved, ndisplay giving
the number of saved scans to be displayed on screen (the function reports
on the screen when every ndisplay iterations have been carried
out), and tune1 and tune2 the parameters needed for the
MH candidate generating distribution for the precision parameters of
the Poisson-Dirichlet process prior. |
state |
a list giving the current value of the parameters. This list is used if the current analysis is the continuation of a previous analysis. |
status |
a logical variable indicating whether this run is new (TRUE ) or the
continuation of a previous analysis (FALSE ). In the latter case
the current value of the parameters must be specified in the
object state . |
work.dir |
working directory. |
This generic function fits a Linear Dependent Poisson-Dirichlet Process Mixture of Survival models as in Jara et al. (2009). The joint distribution of the true chronological onset times and true time-to-events, Ti, is modeled using the mixture model
log(Ti)=yi | fXi ~ fXi
fXi = int N(Xi beta, Sigma) dG(beta)
G | a, b, G0 ~ PD(a, b, G0)
where, G0 = N(beta | mb, Sb). To complete the model specification, independent hyperpriors are assumed,
Sigma | nu, T ~ IW(nu,T)
a | q, a0, b0 ~ q delta0 + (1-q) Beta(a0,b0)
b | a, mub, sigmab ~ N(mub,sigmab)I(-a,infty)
mb | m0, S0 ~ N(m0,S0)
Sb | nub, Tb ~ IW(nub,Tb)
Note that the inverted-Wishart prior is parametrized such that if A ~ IWq(nu, psi) then E(A)= psiinv/(nu-q-1).
The precision parameters are updated using a MH step. The candidate generating distribution is of the form
a' | a, tune1 ~ 0.5 delta0 + 0.5 N(a,tune1**2)
b' | b, a', tune2 ~ N(b,tune2**2)I(-a',infty)
The computational implementation of the model is based on the maxm
truncation of stick-breaking representation of
the model (see, Jara et al., 2009).
An object of class LDPDdoublyint
representing the LDPD mixture of survival models fit.
Generic functions such as print
, plot
,
summary
, and predict
have methods to show the results of the fit. The results include
mb
, Sb
, sigma
, the precision parameter
alpha
, and the number of clusters.
The list state
in the output object contains the current value of the parameters
necessary to restart the analysis. If you want to specify different starting values
to run multiple chains set status=TRUE
and create the list state based on
this starting values. In this case the list state
must include the following objects:
alpha |
a vector giving the value of the precision parameters a and b. |
b |
a matrix of dimension maxm times the number of columns in the design matrix ((p+q)*nvar/2), giving the regression coefficients for each cluster. |
sigma |
a matrix of dimension nvar times nvar giving the kernel covariance matrix. |
mb |
giving the mean of the normal baseline distributions. |
Sb |
giving the covariance matrix the normal baseline distributions. |
ncluster |
an integer giving the number of clusters. |
ss |
an interger vector defining to which of the ncluster clusters each subject belongs. |
y |
a matrix of dimension nsubject times nvar giving the value of the imputed log-survival times. |
Alejandro Jara <ajarav@udec.cl>
Jara, A., Lesaffre, E., De Iorio, M., Quintana, F. (2009). Bayesian semiparametric inference for multivariate doubly-interval-censored data. Technical Report.
## Not run: ############################################################### # HIV-AIDS Data ############################################################### data(hiv) attach(hiv) ############################################################### # Working folder ############################################################### work.dir <- NULL ############################################################### # Response and design matrices ############################################################### nsubject <- length(onsetL) onset <- cbind(onsetL,onsetU) failure <- cbind(failureL,failureU) intercept <- rep(1,nsubject) p <- 2 xonset <- cbind(IntO=intercept,trtO=trt) q <- 2 xfailure <- cbind(IntF=intercept,trtF=trt) ############################################################### # Predictions ############################################################### grid <- matrix(c(rep(seq(0,30,len=30),1), rep(seq(0,20,len=30),1)),nrow=2,byrow=T) xpred <- matrix(c(rep(c(1,0),1),rep(c(1,0),1), rep(c(1,1),1),rep(c(1,1),1)), nrow=2,byrow=T) colnames(xpred) <- colnames(cbind(xonset,xfailure)) ############################################################### # Initial state ############################################################### state <- NULL ############################################################### # Prior ############################################################### prior<-list(a0=1, b0=1, q=0.5, mub=10, sigmab=200, nu=4, tinv=diag(1,2), nub=6, tbinv=diag(1,4), m0=rep(0,4), S0=diag(100,4), maxm=40) ############################################################### # MCMC ############################################################### mcmc<-list(nburn=5000, nskip=9, ndisplay=100, nsave=5000, tune1=0.25, tune2=1) ############################################################### # Fitting the Model ############################################################### fit1 <- LDPDdoublyint(onset=onset,failure=failure, p=p,xonset=xonset, q=q,xfailure=xfailure, xpred=xpred,grid=grid, prior=prior, mcmc=mcmc, state=state, status=TRUE, work.dir=work.dir) fit1 summary(fit1) ############################################################### # Getting Information for Predictions ############################################################### # Without CI bands and intervals fit1.pred <- predict(fit1) fit1.pred plot(fit1.pred) # With CI bands and intervals fit1.pred <- predict(fit1,compute.band=TRUE) fit1.pred plot(fit1.pred) ## End(Not run)