Nothing
R/LDPDdoublyint.R
In DPpackage: Bayesian Nonparametric Modeling in R
### LDPDdoublyint.R
### Fit a semiparametric aft model for multivariate doubly interval
### censored survival data.
###
### Copyright: Alejandro Jara, 2007-2012.
###
### Last modification: 23-01-2010.
###
### This program is free software; you can redistribute it and/or modify
### it under the terms of the GNU General Public License as published by
### the Free Software Foundation; either version 2 of the License, or (at
### your option) any later version.
###
### This program is distributed in the hope that it will be useful, but
### WITHOUT ANY WARRANTY; without even the implied warranty of
### MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
### General Public License for more details.
###
### You should have received a copy of the GNU General Public License
### along with this program; if not, write to the Free Software
### Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
###
### The author's contact information:
###
### Alejandro Jara
### Department of Statistics
### Facultad de Matematicas
### Pontificia Universidad Catolica de Chile
### Casilla 306, Correo 22
### Santiago
### Chile
### Voice: +56-2-3544506 URL : http://www.mat.puc.cl/~ajara
### Fax : +56-2-3547729 Email: atjara@uc.cl
###
"LDPDdoublyint"<-
function(onset,failure,p,xonset,q,xfailure,xpred,grid,prior,mcmc,state,status,work.dir=NULL)
UseMethod("LDPDdoublyint")
"LDPDdoublyint.default"<-
function(onset,
failure,
p,
xonset,
q,
xfailure,
xpred,
grid,
prior,
mcmc,
state,
status,
work.dir=NULL)
{
#########################################################################################
# data structure
#########################################################################################
cl <- match.call()
nsubject <- nrow(onset)
nvar <- ncol(onset)/2
x <- cbind(xonset,xfailure)
y <- NULL
llower <- NULL
lupper <- NULL
j <- 0
for(i in 1:nvar)
{
j <- j+1
llower <- cbind(llower,onset[,j])
y <- cbind(y,log((onset[,j]+onset[,j+1])/2))
j <- j+1
lupper <- cbind(lupper,onset[,j])
}
j <- 0
for(i in 1:nvar)
{
j <- j+1
llower <- cbind(llower,failure[,j])
tmp <- ((failure[,j]+failure[,j+1])/2)-((onset[,j]+onset[,j+1])/2)
tmp[tmp==0] <- 0.01
y <- cbind(y,log(tmp))
j <- j+1
lupper <- cbind(lupper,failure[,j])
}
npred <- nrow(xpred)
ngrid <- ncol(grid)
#########################################################################################
# change working directory (if requested..)
#########################################################################################
if(!is.null(work.dir))
{
cat("\n Changing working directory to ",work.dir,"\n")
old.dir <- getwd() # by default work in current working directory
setwd(work.dir)
}
#########################################################################################
# Prior specification
#########################################################################################
maxm <- prior$maxm
nu <- c(prior$nu,prior$nub)
if(is.null(prior$a0))
{
a0b0 <- c(-1,-1,0)
ap <- prior$a
apr <- 0
}
else
{
a0b0 <- c(prior$a0,prior$b0,prior$q)
ap <- 0
apr <- 1
if(prior$a0<0 || prior$b0<0)
{
stop("The parameters of the Beta prior for the a parameter must be possitive.\n")
}
if(prior$q<0 || prior$q>1)
{
stop("The mixture proportion must be in (0,1).\n")
}
}
if(is.null(prior$mub))
{
a0b0 <- c(a0b0,0,-1)
bp <- prior$b
bpr <- 0
}
else
{
a0b0 <- c(a0b0,prior$mub,prior$sigmab)
bpr <- 1
bp <- abs(rnorm(1,mean=prior$mub,sd=sqrt(prior$sigmab)))
}
alpha <- c(ap,bp)
tinv1 <- prior$tinv
tinv2 <- prior$tbinv
psiinv <- solve(prior$S0)
smu <- psiinv%*%prior$m0
#########################################################################################
# mcmc specification
#########################################################################################
mcmcvec <- c(mcmc$nburn,mcmc$nskip,mcmc$ndisplay)
nsave <- mcmc$nsave
if(is.null(mcmc$tune1))
{
tune1 <- 1.1
}
else
{
tune1 <- mcmc$tune1
}
if(is.null(mcmc$tune2))
{
tune2 <- 1.1
}
else
{
tune2 <- mcmc$tune2
}
a0b0 <- c(a0b0,tune1,tune2)
#########################################################################################
# Starting values
#########################################################################################
sigma <- matrix(0,nrow=2*nvar,ncol=2*nvar)
beta <- NULL
sigma <- NULL
for(i in 1:nvar)
{
startp <- (i-1)*p+1
endp<- (i-1)*p+p
fit0 <- glm.fit(xonset[,(startp:endp)], y[,i], family=gaussian())
betatmp <- coefficients(fit0)
beta <- c(beta,betatmp)
s2 <- var(residuals(fit0))
sigma <- c(sigma,s2)
}
for(i in 1:nvar)
{
startp <- (i-1)*q+1
endp <- (i-1)*q+q
fit0 <- glm.fit(xfailure[,(startp:endp)], y[,nvar+i], family=gaussian())
betatmp <- coefficients(fit0)
beta <- c(beta,betatmp)
s2 <- var(residuals(fit0))
sigma <- c(sigma,s2)
}
sigma <- diag(sigma,2*nvar)
mub <- beta
sigmab <- diag(1,nvar*(p+q))
#########################################################################################
# output
#########################################################################################
acrate <- rep(0,2)
npar <- (2*nvar*(2*nvar+1)/2)+nvar*(p+q)+(nvar*(p+q)*(nvar*(p+q)+1)/2)+3
f <- matrix(0,nrow=npred*2*nvar,ncol=ngrid)
h <- matrix(0,nrow=npred*2*nvar,ncol=ngrid)
thetasave <- matrix(0,nrow=nsave,ncol=npar)
randsave <- matrix(0,nrow=nsave,ncol=npred*2*nvar)
#########################################################################################
# parameters depending on status
#########################################################################################
if(status==TRUE)
{
sigmainv <- solve(sigma)
sigmabinv <- solve(sigmab)
b <- matrix(0,nrow=maxm,ncol=nvar*(p+q))
for(i in 1:1)
{
b[i,] <- beta
}
for(i in 2:maxm)
{
b[i,] <- rep(0,nvar*(p+q))
}
ss <- rep(1,nsubject)
ncluster <- 1
}
if(status==FALSE)
{
alpha <- state$alpha
b <- state$b
sigma <- state$sigma
sigmainv <- solve(state$sigma)
mub <- state$mb
sigmab <- state$Sb
sigmabinv <- solve(state$Sb)
y <- state$y
ss <- state$ss
ncluster <- state$ncluster
}
#########################################################################################
# working space
#########################################################################################
cstrt <- matrix(0,nrow=maxm,ncol=nsubject)
ccluster <- rep(0,maxm)
prob <- rep(0,maxm)
prob2 <- rep(0,maxm)
weights <- rep(0,maxm)
iflagc <- rep(0,nvar*(p+q))
theta <- rep(0,nvar*(p+q))
workmc <- matrix(0,nrow=nvar*(p+q),ncol=nvar*(p+q))
workmc2 <- matrix(0,nrow=nvar*(p+q),ncol=nvar*(p+q))
workmc3 <- matrix(0,nrow=nvar*(p+q),ncol=nvar*(p+q))
workmhc <- rep(0,nvar*(p+q)*(nvar*(p+q)+1)/2)
workmhc2 <- rep(0,nvar*(p+q)*(nvar*(p+q)+1)/2)
workvc <- rep(0,nvar*(p+q))
iflagn <- rep(0,2*nvar)
workmn <- matrix(0,nrow=2*nvar,ncol=2*nvar)
workmn2 <- matrix(0,nrow=2*nvar,ncol=2*nvar)
workmn3 <- matrix(0,nrow=2*nvar,ncol=2*nvar)
workmhn <- rep(0,2*nvar*(2*nvar+1)/2)
workmhn2 <- rep(0,2*nvar*(2*nvar+1)/2)
workvn <- rep(0,2*nvar)
workvn2 <- rep(0,2*nvar)
workvn3 <- rep(0,2*nvar)
ztz <- matrix(0,nrow=nvar*(p+q),ncol=nvar*(p+q))
zty <- rep(0,nvar*(p+q))
seed <- c(sample(1:29000,1),sample(1:29000,1))
model <- matrix(0,nrow=2*nvar,ncol=(p+q))
possi <- matrix(0,nrow=2*nvar,ncol=(p+q))
fw <- rep(0,ngrid)
fw2 <- rep(0,ngrid)
varind <- rep(1,p)
for(i in 2:nvar)
{
varind <- c(varind,rep(i,p))
}
for(i in 1:nvar)
{
varind <- c(varind,rep(nvar+i,q))
}
#########################################################################################
# calling the fortran code
#########################################################################################
foo <- .Fortran("ldpddoublyintsba",
nsubject =as.integer(nsubject),
nvar =as.integer(nvar),
p =as.integer(p),
q =as.integer(q),
ngrid =as.integer(ngrid),
npred =as.integer(npred),
x =as.double(x),
llower =as.double(llower),
lupper =as.double(lupper),
grid =as.double(grid),
xpred =as.double(xpred),
maxm =as.integer(maxm),
a0b0 =as.double(a0b0),
nu =as.integer(nu),
tinv1 =as.double(tinv1),
smu =as.double(smu),
psiinv =as.double(psiinv),
tinv2 =as.double(tinv2),
mcmc =as.integer(mcmcvec),
nsave =as.integer(nsave),
ncluster =as.integer(ncluster),
ss =as.integer(ss),
alpha =as.double(alpha),
b =as.double(b),
sigma =as.double(sigma),
sigmainv =as.double(sigmainv),
mub =as.double(mub),
sigmab =as.double(sigmab),
sigmabinv =as.double(sigmabinv),
y =as.double(y),
acrate =as.double(acrate),
f =as.double(f),
h =as.double(h),
thetasave =as.double(thetasave),
randsave =as.double(randsave),
seed =as.integer(seed),
model =as.integer(model),
possi =as.integer(possi),
varind =as.integer(varind),
cstrt =as.integer(cstrt),
ccluster =as.integer(ccluster),
prob =as.double(prob),
prob2 =as.double(prob2),
weights =as.double(weights),
iflagc =as.integer(iflagc),
theta =as.double(theta),
workmc =as.double(workmc),
workmc2 =as.double(workmc2),
workmc3 =as.double(workmc3),
workmhc =as.double(workmhc),
workmhc2 =as.double(workmhc2),
workvc =as.double(workvc),
iflagn =as.integer(iflagn),
workmn =as.double(workmn),
workmn2 =as.double(workmn2),
workmn3 =as.double(workmn3),
workmhn =as.double(workmhn),
workmhn2 =as.double(workmhn2),
workvn =as.double(workvn),
workvn2 =as.double(workvn2),
workvn3 =as.double(workvn3),
ztz =as.double(ztz),
zty =as.double(zty),
fw =as.double(fw),
fw2 =as.double(fw2),
PACKAGE="DPpackage")
#########################################################################################
# save state
#########################################################################################
model.name<-"Linear Dependent Poisson-Dirichlet Survival Model for Doubly-Interval Censored Data"
nametmp<-c(paste("onset",1:nvar,sep=""),
paste("time",1:nvar,sep="")
)
pnames<-NULL
for(i in 1:(2*nvar))
{
for(j in i:(2*nvar))
{
tmp<-paste(nametmp[i],nametmp[j],sep=";")
tmp<-paste("Sigma (",tmp,")",sep="")
pnames<-c(pnames,tmp)
}
}
nametmp <- colnames(x)
pnames <- c(pnames,paste("mb",nametmp,sep="-"))
for(i in 1:(nvar*(p+q)))
{
for(j in i:(nvar*(p+q)))
{
tmp<-paste(nametmp[i],nametmp[j],sep=";")
tmp<-paste("Sb (",tmp,")",sep="")
pnames<-c(pnames,tmp)
}
}
pnames <- c(pnames,"ncluster","a","b")
thetasave <- matrix(foo$thetasave,nrow=nsave,ncol=npar)
colnames(thetasave) <- pnames
coeff <- apply(thetasave,2,mean)
nametmp<-c(paste("onset",1:nvar,sep=""),
paste("time",1:nvar,sep="")
)
nametmp <- rep(nametmp,npred)
randsave <- matrix(foo$randsave,nrow=nsave,ncol=npred*2*nvar)
colnames(randsave)<-nametmp
coeffr <- apply(randsave,2,mean)
medtable <- matrix(coeffr,nrow=npred,ncol=2*nvar,byrow=T)
nametmp<-c(paste("onset",1:nvar,sep=""),
paste("time",1:nvar,sep="")
)
colnames(medtable)<-nametmp
nametmp<-paste("prediction",1:npred)
rownames(medtable)<-nametmp
save.state <- list(thetasave=thetasave,
randsave=randsave)
state <- list( alpha=foo$alpha,
b=matrix(foo$b,nrow=maxm,ncol=nvar*(p+q)),
sigma=matrix(foo$sigma,nrow=2*nvar,ncol=2*nvar),
mb=foo$mub,
Sb=matrix(foo$sigmab,nrow=nvar*(p+q),ncol=nvar*(p+q)),
y=matrix(foo$y,nrow=nsubject,ncol=2*nvar),
ss=foo$ss,
ncluster=foo$ncluster)
acrate <- foo$acrate[1]
dptest <- foo$acrate[2]
z<-list(modelname=model.name,
coefficients=coeff,
call=cl,
acrate=acrate,
dptest=dptest,
prior=prior,
mcmc=mcmc,
state=state,
save.state=save.state,
nsubject=foo$nsubject,
nsave=nsave,
p=foo$p,
q=foo$q,
nvar=2*foo$nvar,
npred=foo$npred,
ngrid=ngrid,
grid=grid,
f=matrix(foo$f,nrow=npred*2*nvar,ncol=ngrid),
h=matrix(foo$h,nrow=npred*2*nvar,ncol=ngrid),
medtable=medtable,
acrate=foo$acrate,
work.dir=work.dir)
cat("\n\n")
class(z)<-c("LDPDdoublyint")
z
}
###
### Tools for LDPDdoublyint: print, summary, plot, predict
###
### Copyright: Alejandro Jara, 2007
### Last modification: 08-13-2007.
"print.LDPDdoublyint"<-function (x, digits = max(3, getOption("digits") - 3), ...)
{
cat("\n",x$modelname,"\n\nCall:\n", sep = "")
print(x$call)
cat("\nP(DP|data)=",x$dptest,"\n")
cat("\nPosterior Inference of Parameters:\n")
print.default(format(x$coefficients, digits = digits), print.gap = 2,
quote = FALSE)
if(!is.null(x$acrate))
{
cat("\nAcceptance Rate for Metropolis Steps = ",x$acrate)
}
cat("\n\nNumber of Subjects:",x$nsubject)
cat("\nNumber of Variables:",x$nvar)
cat("\n\n")
invisible(x)
}
"summary.LDPDdoublyint"<-function(object, hpd=TRUE, ...)
{
stde<-function(x)
{
n<-length(x)
return(sd(x)/sqrt(n))
}
hpdf<-function(x)
{
alpha<-0.05
vec<-x
n<-length(x)
alow<-rep(0,2)
aupp<-rep(0,2)
a<-.Fortran("hpd",n=as.integer(n),alpha=as.double(alpha),x=as.double(vec),
alow=as.double(alow),aupp=as.double(aupp),PACKAGE="DPpackage")
return(c(a$alow[1],a$aupp[1]))
}
pdf<-function(x)
{
alpha<-0.05
vec<-x
n<-length(x)
alow<-rep(0,2)
aupp<-rep(0,2)
a<-.Fortran("hpd",n=as.integer(n),alpha=as.double(alpha),x=as.double(vec),
alow=as.double(alow),aupp=as.double(aupp),PACKAGE="DPpackage")
return(c(a$alow[2],a$aupp[2]))
}
thetasave<-object$save.state$thetasave
ans <- c(object[c("call", "modelname")])
### Kernel
nvar <- object$nvar
dimen1 <- nvar*(nvar+1)/2
coef.p <- object$coefficients[1:dimen1]
mat <- thetasave[,1:dimen1]
if(dimen1>0)
{
coef.m <-apply(mat, 2, median)
coef.sd<-apply(mat, 2, sd)
coef.se<-apply(mat, 2, stde)
if(hpd){
limm<-apply(mat, 2, hpdf)
coef.l<-limm[1,]
coef.u<-limm[2,]
}
else
{
limm<-apply(mat, 2, pdf)
coef.l<-limm[1,]
coef.u<-limm[2,]
}
coef.table <- cbind(coef.p, coef.m, coef.sd, coef.se , coef.l , coef.u)
if(hpd)
{
dimnames(coef.table) <- list(names(coef.p), c("Mean", "Median", "Std. Dev.", "Naive Std.Error",
"95%HPD-Low","95%HPD-Upp"))
}
else
{
dimnames(coef.table) <- list(names(coef.p), c("Mean", "Median", "Std. Dev.", "Naive Std.Error",
"95%CI-Low","95%CI-Upp"))
}
ans$kernel <- coef.table
}
### Baseline
p <- object$p
q <- object$q
nt <- (p+q)*(nvar/2)
dimen2 <- nt + nt*(nt+1)/2
coef.p <- object$coefficients[(dimen1+1):(dimen1+dimen2)]
mat <- thetasave[,(dimen1+1):(dimen1+dimen2)]
if(dimen2>0)
{
coef.m <-apply(mat, 2, median)
coef.sd<-apply(mat, 2, sd)
coef.se<-apply(mat, 2, stde)
if(hpd){
limm<-apply(mat, 2, hpdf)
coef.l<-limm[1,]
coef.u<-limm[2,]
}
else
{
limm<-apply(mat, 2, pdf)
coef.l<-limm[1,]
coef.u<-limm[2,]
}
coef.table <- cbind(coef.p, coef.m, coef.sd, coef.se , coef.l , coef.u)
if(hpd)
{
dimnames(coef.table) <- list(names(coef.p), c("Mean", "Median", "Std. Dev.", "Naive Std.Error",
"95%HPD-Low","95%HPD-Upp"))
}
else
{
dimnames(coef.table) <- list(names(coef.p), c("Mean", "Median", "Std. Dev.", "Naive Std.Error",
"95%CI-Low","95%CI-Upp"))
}
ans$base <- coef.table
}
### Precision
dimen3 <- 3
coef.p <- object$coefficients[(dimen1+dimen2+1):(dimen1+dimen2+dimen3)]
mat <- thetasave[,(dimen1+dimen2+1):(dimen1+dimen2+dimen3)]
if(dimen3>0)
{
coef.m <-apply(mat, 2, median)
coef.sd<-apply(mat, 2, sd)
coef.se<-apply(mat, 2, stde)
if(hpd){
limm<-apply(mat, 2, hpdf)
coef.l<-limm[1,]
coef.u<-limm[2,]
}
else
{
limm<-apply(mat, 2, pdf)
coef.l<-limm[1,]
coef.u<-limm[2,]
}
coef.table <- cbind(coef.p, coef.m, coef.sd, coef.se , coef.l , coef.u)
if(hpd)
{
dimnames(coef.table) <- list(names(coef.p), c("Mean", "Median", "Std. Dev.", "Naive Std.Error",
"95%HPD-Low","95%HPD-Upp"))
}
else
{
dimnames(coef.table) <- list(names(coef.p), c("Mean", "Median", "Std. Dev.", "Naive Std.Error",
"95%CI-Low","95%CI-Upp"))
}
ans$prec <- coef.table
}
ans$nsubject <- object$nsubject
ans$nvar <- object$nvar
ans$acrate <- object$acrate
ans$dptest <- object$dptest
class(ans) <- "summaryLDPDdoublyint"
return(ans)
}
"print.summaryLDPDdoublyint"<-function (x, digits = max(3, getOption("digits") - 3), ...)
{
cat("\n",x$modelname,"\n\nCall:\n", sep = "")
print(x$call)
cat("\nP(DP|data)=",x$dptest,"\n")
cat("\nKernel covariance matrix:\n")
print.default(format(x$kernel, digits = digits), print.gap = 2,
quote = FALSE)
cat("\nCentering parameters:\n")
print.default(format(x$base, digits = digits), print.gap = 2,
quote = FALSE)
cat("\nPrecision parameters:\n")
print.default(format(x$prec, digits = digits), print.gap = 2,
quote = FALSE)
if(!is.null(x$acrate))
{
cat("\nAcceptance Rate for Metropolis Steps = ",x$acrate)
}
cat("\n\nNumber of Subject:",x$nsubject)
cat("\nNumber of Variable:",x$nvar)
cat("\n\n")
invisible(x)
}
"plot.LDPDdoublyint"<-function(x, hpd=TRUE, ask=TRUE, nfigr=2, nfigc=2, param=NULL, col="#bdfcc9", ...)
{
fancydensplot<-function(x, hpd=TRUE, npts=200, xlab="", ylab="", main="",col="#bdfcc9", ...)
# Author: AJV, 2006
#
{
dens <- density(x,n=npts)
densx <- dens$x
densy <- dens$y
meanvar <- mean(x)
densx1 <- max(densx[densx<=meanvar])
densx2 <- min(densx[densx>=meanvar])
densy1 <- densy[densx==densx1]
densy2 <- densy[densx==densx2]
ymean <- densy1 + ((densy2-densy1)/(densx2-densx1))*(meanvar-densx1)
if(hpd==TRUE)
{
alpha<-0.05
alow<-rep(0,2)
aupp<-rep(0,2)
n<-length(x)
a<-.Fortran("hpd",n=as.integer(n),alpha=as.double(alpha),x=as.double(x),
alow=as.double(alow),aupp=as.double(aupp),PACKAGE="DPpackage")
xlinf<-a$alow[1]
xlsup<-a$aupp[1]
}
else
{
xlinf <- quantile(x,0.025)
xlsup <- quantile(x,0.975)
}
densx1 <- max(densx[densx<=xlinf])
densx2 <- min(densx[densx>=xlinf])
densy1 <- densy[densx==densx1]
densy2 <- densy[densx==densx2]
ylinf <- densy1 + ((densy2-densy1)/(densx2-densx1))*(xlinf-densx1)
densx1 <- max(densx[densx<=xlsup])
densx2 <- min(densx[densx>=xlsup])
densy1 <- densy[densx==densx1]
densy2 <- densy[densx==densx2]
ylsup <- densy1 + ((densy2-densy1)/(densx2-densx1))*(xlsup-densx1)
plot(0.,0.,xlim = c(min(densx), max(densx)), ylim = c(min(densy), max(densy)),
axes = F,type = "n" , xlab=xlab, ylab=ylab, main=main, cex=1.2)
xpol<-c(xlinf,xlinf,densx[densx>=xlinf & densx <=xlsup],xlsup,xlsup)
ypol<-c(0,ylinf,densy[densx>=xlinf & densx <=xlsup] ,ylsup,0)
polygon(xpol, ypol, border = FALSE,col=col)
lines(c(min(densx), max(densx)),c(0,0),lwd=1.2)
segments(min(densx),0, min(densx),max(densy),lwd=1.2)
lines(densx,densy,lwd=1.2)
segments(meanvar, 0, meanvar, ymean,lwd=1.2)
segments(xlinf, 0, xlinf, ylinf,lwd=1.2)
segments(xlsup, 0, xlsup, ylsup,lwd=1.2)
axis(1., at = round(c(xlinf, meanvar,xlsup), 2.), labels = T,pos = 0.)
axis(1., at = round(seq(min(densx),max(densx),length=15), 2.), labels = F,pos = 0.)
axis(2., at = round(seq(0,max(densy),length=5), 2.), labels = T,pos =min(densx))
}
if(is(x, "LDPDdoublyint"))
{
if(is.null(param))
{
coef.p<-x$coefficients
n<-length(coef.p)
pnames<-names(coef.p)
par(ask = ask)
layout(matrix(seq(1,nfigr*nfigc,1), nrow=nfigr, ncol=nfigc, byrow = TRUE))
for(i in 1:(n-2)){
title1<-paste("Trace of",pnames[i],sep=" ")
title2<-paste("Density of",pnames[i],sep=" ")
plot(x$save.state$thetasave[,i],type='l',main=title1,xlab="MCMC scan",ylab=" ")
if(pnames[i]=="ncluster" || pnames[i]=="k")
{
hist(x$save.state$thetasave[,i],main=title2,xlab="values", ylab="probability",probability=TRUE)
}
else
{
fancydensplot(x$save.state$thetasave[,i],hpd=hpd,main=title2,xlab="values", ylab="density",col=col)
}
}
if(!is.null(x$prior$a0))
{
title1<-paste("Trace of",pnames[n-1],sep=" ")
title2<-paste("Density of",pnames[n-1],sep=" ")
plot(x$save.state$thetasave[,n-1],type='l',main=title1,xlab="MCMC scan",ylab=" ")
fancydensplot(x$save.state$thetasave[,n-1],hpd=hpd,main=title2,xlab="values", ylab="density",col=col)
}
if(!is.null(x$prior$mub))
{
title1<-paste("Trace of",pnames[n],sep=" ")
title2<-paste("Density of",pnames[n],sep=" ")
plot(x$save.state$thetasave[,n],type='l',main=title1,xlab="MCMC scan",ylab=" ")
fancydensplot(x$save.state$thetasave[,n],hpd=hpd,main=title2,xlab="values", ylab="density",col=col)
}
}
else
{
coef.p <- x$coefficients
n <- length(coef.p)
pnames<-names(coef.p)
poss<-0
for(i in 1:n)
{
if(pnames[i]==param)poss=i
}
if (poss==0 && param !="predictive")
{
stop("This parameter is not present in the original model.\n")
}
if (param !="ncluster")
{
par(ask = ask)
layout(matrix(seq(1,nfigr*nfigc,1), nrow=nfigr, ncol=nfigc, byrow = TRUE))
title1<-paste("Trace of",pnames[poss],sep=" ")
title2<-paste("Density of",pnames[poss],sep=" ")
plot(x$save.state$thetasave[,poss],type='l',main=title1,xlab="MCMC scan",ylab=" ")
fancydensplot(x$save.state$thetasave[,poss],hpd=hpd,main=title2,xlab="values", ylab="density",col=col)
}
else
{
par(ask = ask)
layout(matrix(seq(1,nfigr*nfigc,1), nrow=nfigr, ncol=nfigc, byrow = TRUE))
title1<-paste("Trace of",pnames[poss],sep=" ")
title2<-paste("Density of",pnames[poss],sep=" ")
plot(x$save.state$thetasave[,poss],type='l',main=title1,xlab="MCMC scan",ylab=" ")
hist(x$save.state$thetasave[,poss],main=title2,xlab="values", ylab="probability",probability=TRUE)
}
}
}
}
"predict.LDPDdoublyint"<-
function(object,grid,compute.band=FALSE,hpd=TRUE, ...)
{
hpdf<-function(x)
{
alpha<-0.05
vec<-x
n<-length(x)
alow<-rep(0,2)
aupp<-rep(0,2)
a<-.Fortran("hpd",n=as.integer(n),alpha=as.double(alpha),x=as.double(vec),
alow=as.double(alow),aupp=as.double(aupp),PACKAGE="DPpackage")
return(c(a$alow[1],a$aupp[1]))
}
pdf<-function(x)
{
alpha<-0.05
vec<-x
n<-length(x)
alow<-rep(0,2)
aupp<-rep(0,2)
a<-.Fortran("hpd",n=as.integer(n),alpha=as.double(alpha),x=as.double(vec),
alow=as.double(alow),aupp=as.double(aupp),PACKAGE="DPpackage")
return(c(a$alow[2],a$aupp[2]))
}
if(is(object, "LDPDdoublyint"))
{
ngrid <- object$ngrid
grid <- object$grid
nvar <- object$nvar/2
npred <- object$npred
f <- object$f
h <- object$h
med <- object$medtable
if(compute.band)
{
nsave <- object$nsave
alpha <- 0.05
tint <- 0
if(hpd)tint <- 1
llower <- NULL
lupper <- NULL
llower <- matrix(0,nrow=npred*2*nvar,ncol=ngrid)
lupper <- matrix(0,nrow=npred*2*nvar,ncol=ngrid)
llower2 <- matrix(0,nrow=npred*2*nvar,ncol=ngrid)
lupper2 <- matrix(0,nrow=npred*2*nvar,ncol=ngrid)
workv1 <- rep(0,nsave)
workv2 <- rep(0,ngrid)
foo <- .Fortran("hpdspy",
nsave =as.integer(nsave),
npred =as.integer(npred),
nvar =as.integer(nvar),
ngrid =as.integer(ngrid),
alpha =as.double(alpha),
tint =as.integer(tint),
workv1 =as.double(workv1),
workv2 =as.double(workv2),
llower =as.double(llower),
lupper =as.double(lupper),
llower2 =as.double(llower2),
lupper2 =as.double(lupper2),
PACKAGE="DPpackage")
llower <- matrix(foo$llower,nrow=npred*2*nvar,ncol=ngrid)
lupper <- matrix(foo$lupper,nrow=npred*2*nvar,ncol=ngrid)
llower2 <- matrix(foo$llower2,nrow=npred*2*nvar,ncol=ngrid)
lupper2 <- matrix(foo$lupper2,nrow=npred*2*nvar,ncol=ngrid)
if(hpd)
{
limm <- apply(object$save.state$randsave, 2, hpdf)
coef.l <- limm[1,]
coef.u <- limm[2,]
}
else
{
limm <- apply(object$save.state$randsave, 2, pdf)
coef.l <- limm[1,]
coef.u <- limm[2,]
}
med.l <- matrix(coef.l,nrow=npred,ncol=2*nvar,byrow=T)
med.u <- matrix(coef.u,nrow=npred,ncol=2*nvar,byrow=T)
colnames(med.l) <- colnames(med)
colnames(med.u) <- colnames(med)
rownames(med.l) <- rownames(med)
rownames(med.u) <- rownames(med)
}
else
{
llower <- NULL
lupper <- NULL
llower2 <- NULL
lupper2 <- NULL
med.l <- NULL
med.u <- NULL
}
out <- list(f=f,
h=h,
med=med,
med.l=med.l,
med.u=med.u,
grid=grid,
ngrid=ngrid,
nvar=nvar,
npred=npred,
llower=llower,
lupper=lupper,
llower2=llower2,
lupper2=lupper2,
modelname=object$modelname,
call=object$call,
compute.band=compute.band)
class(out) <- c("predict.LDPDdoublyint")
out
}
}
"print.predict.LDPDdoublyint"<-function (x, digits = max(3, getOption("digits") - 3), ...)
{
cat("\n",x$modelname,"\n\nCall:\n", sep = "")
print(x$call)
cat("\nPosterior Mean of Median Survival for Predictions:\n")
print.default(format(x$med, digits = digits), print.gap = 2,
quote = FALSE)
if(!is.null(x$med.l))
{
cat("\nLower CI Limit for the Median Survival for Predictions:\n")
print.default(format(x$med.l, digits = digits), print.gap = 2,
quote = FALSE)
}
if(!is.null(x$med.u))
{
cat("\nUpper CI Limit for the Median Survival for Predictions:\n")
print.default(format(x$med.u, digits = digits), print.gap = 2,
quote = FALSE)
}
cat("\n\n")
invisible(x)
}
plot.predict.LDPDdoublyint<-function(x,ask=TRUE,xlim=NULL,nfigr=1,nfigc=1, ...)
{
if(is(x, "predict.LDPDdoublyint"))
{
par(ask = ask)
layout(matrix(seq(1,nfigr*nfigc,1), nrow=nfigr, ncol=nfigc, byrow = TRUE))
compute.band <- x$compute.band
if(compute.band)
{
k <- 0
for(i in 1:x$npred)
{
for(j in 1:x$nvar)
{
tmp1 <- paste("Prediction",i,sep=":")
tmp2 <- paste("Onset",j,sep=":")
tmp3 <- paste(tmp1,tmp2,sep="-")
k <- k+1
x1 <- x$grid[j,]
y1 <- x$llower[k,]
x2 <- x$grid[j,]
y2 <- x$lupper[k,]
aa <- rbind(x2,y2)[, order(-x2, y2)]
x2 <- aa[1,]
y2 <- aa[2,]
plot(x$grid[j,],x$f[k,],type="l",lty=1,lwd=3,ylim=c(0,1),
ylab="Survival",xlab="Time",main=tmp3,bty="n")
polygon(x=c(x1,x2),y=c(y1,y2),border=NA,col="lightgray")
lines(x$grid[j,],x$f[k,],lwd=3,lty=1)
}
for(j in 1:x$nvar)
{
tmp1 <- paste("Prediction",i,sep=":")
tmp2 <- paste("Failure",j,sep=":")
tmp3 <- paste(tmp1,tmp2,sep="-")
k <- k+1
x1 <- x$grid[j,]
y1 <- x$llower[k,]
x2 <- x$grid[j,]
y2 <- x$lupper[k,]
aa <- rbind(x2,y2)[, order(-x2, y2)]
x2 <- aa[1,]
y2 <- aa[2,]
plot(x$grid[j,],x$f[k,],type="l",lty=1,lwd=3,ylim=c(0,1),
ylab="Survival",xlab="Time",main=tmp3,bty="n")
polygon(x=c(x1,x2),y=c(y1,y2),border=NA,col="lightgray")
lines(x$grid[j,],x$f[k,],lwd=3,lty=1)
}
}
k <- 0
for(i in 1:x$npred)
{
for(j in 1:x$nvar)
{
tmp1 <- paste("Prediction",i,sep=":")
tmp2 <- paste("Onset",j,sep=":")
tmp3 <- paste(tmp1,tmp2,sep="-")
k <- k+1
x1 <- grid[j,]
y1 <- x$llower2[k,]
x2 <- grid[j,]
y2 <- x$lupper2[k,]
aa <- rbind(x2,y2)[, order(-x2, y2)]
x2 <- aa[1,]
y2 <- aa[2,]
plot(x$grid[j,],x$lupper2[k,],type="l",lty=1,lwd=2,
ylab="Hazard",xlab="Time",main=tmp3,bty="n",col="lightgray")
polygon(x=c(x1,x2),y=c(y1,y2),border=NA,col="lightgray")
lines(x$grid[j,],x$h[k,],lwd=3,lty=1)
}
for(j in 1:x$nvar)
{
tmp1 <- paste("Prediction",i,sep=":")
tmp2 <- paste("Failure",j,sep=":")
tmp3 <- paste(tmp1,tmp2,sep="-")
k <- k+1
x1 <- grid[j,]
y1 <- x$llower2[k,]
x2 <- grid[j,]
y2 <- x$lupper2[k,]
aa <- rbind(x2,y2)[, order(-x2, y2)]
x2 <- aa[1,]
y2 <- aa[2,]
plot(x$grid[j,],x$lupper2[k,],type="l",lty=1,lwd=2,
ylab="Hazard",xlab="Time",main=tmp3,bty="n",col="lightgray")
polygon(x=c(x1,x2),y=c(y1,y2),border=NA,col="lightgray")
lines(x$grid[j,],x$h[k,],lwd=3,lty=1)
}
}
}
else
{
k <- 0
for(i in 1:x$npred)
{
for(j in 1:x$nvar)
{
tmp1 <- paste("Prediction",i,sep=":")
tmp2 <- paste("Onset",j,sep=":")
tmp3 <- paste(tmp1,tmp2,sep="-")
k <- k+1
plot(x$grid[j,],x$f[k,],type="l",lty=1,lwd=3,ylim=c(0,1),
ylab="Survival",xlab="Time",main=tmp3,bty="n")
}
for(j in 1:x$nvar)
{
tmp1 <- paste("Prediction",i,sep=":")
tmp2 <- paste("Failure",j,sep=":")
tmp3 <- paste(tmp1,tmp2,sep="-")
k <- k+1
plot(x$grid[j,],x$f[k,],type="l",lty=1,lwd=3,ylim=c(0,1),
ylab="Survival",xlab="Time",main=tmp3,bty="n")
}
}
k <- 0
for(i in 1:x$npred)
{
for(j in 1:x$nvar)
{
tmp1 <- paste("Prediction",i,sep=":")
tmp2 <- paste("Onset",j,sep=":")
tmp3 <- paste(tmp1,tmp2,sep="-")
k <- k+1
plot(x$grid[j,],x$h[k,],type="l",lty=1,lwd=3,
ylab="Hazard",xlab="Time",main=tmp3,bty="n")
}
for(j in 1:x$nvar)
{
tmp1 <- paste("Prediction",i,sep=":")
tmp2 <- paste("Failure",j,sep=":")
tmp3 <- paste(tmp1,tmp2,sep="-")
k <- k+1
plot(x$grid[j,],x$h[k,],type="l",lty=1,lwd=3,
ylab="Hazard",xlab="Time",main=tmp3,bty="n")
}
}
}
}
}
Try the DPpackage package in your browser
Any scripts or data that you put into this service are public.
DPpackage documentation built on May 1, 2019, 10:23 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
### LDPDdoublyint.R
### Fit a semiparametric aft model for multivariate doubly interval
### censored survival data.
###
### Copyright: Alejandro Jara, 2007-2012.
###
### Last modification: 23-01-2010.
###
### This program is free software; you can redistribute it and/or modify
### it under the terms of the GNU General Public License as published by
### the Free Software Foundation; either version 2 of the License, or (at
### your option) any later version.
###
### This program is distributed in the hope that it will be useful, but
### WITHOUT ANY WARRANTY; without even the implied warranty of
### MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
### General Public License for more details.
###
### You should have received a copy of the GNU General Public License
### along with this program; if not, write to the Free Software
### Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
###
### The author's contact information:
###
### Alejandro Jara
### Department of Statistics
### Facultad de Matematicas
### Pontificia Universidad Catolica de Chile
### Casilla 306, Correo 22
### Santiago
### Chile
### Voice: +56-2-3544506 URL : http://www.mat.puc.cl/~ajara
### Fax : +56-2-3547729 Email: atjara@uc.cl
###
"LDPDdoublyint"<-
function(onset,failure,p,xonset,q,xfailure,xpred,grid,prior,mcmc,state,status,work.dir=NULL)
UseMethod("LDPDdoublyint")
"LDPDdoublyint.default"<-
function(onset,
failure,
p,
xonset,
q,
xfailure,
xpred,
grid,
prior,
mcmc,
state,
status,
work.dir=NULL)
{
#########################################################################################
# data structure
#########################################################################################
cl <- match.call()
nsubject <- nrow(onset)
nvar <- ncol(onset)/2
x <- cbind(xonset,xfailure)
y <- NULL
llower <- NULL
lupper <- NULL
j <- 0
for(i in 1:nvar)
{
j <- j+1
llower <- cbind(llower,onset[,j])
y <- cbind(y,log((onset[,j]+onset[,j+1])/2))
j <- j+1
lupper <- cbind(lupper,onset[,j])
}
j <- 0
for(i in 1:nvar)
{
j <- j+1
llower <- cbind(llower,failure[,j])
tmp <- ((failure[,j]+failure[,j+1])/2)-((onset[,j]+onset[,j+1])/2)
tmp[tmp==0] <- 0.01
y <- cbind(y,log(tmp))
j <- j+1
lupper <- cbind(lupper,failure[,j])
}
npred <- nrow(xpred)
ngrid <- ncol(grid)
#########################################################################################
# change working directory (if requested..)
#########################################################################################
if(!is.null(work.dir))
{
cat("\n Changing working directory to ",work.dir,"\n")
old.dir <- getwd() # by default work in current working directory
setwd(work.dir)
}
#########################################################################################
# Prior specification
#########################################################################################
maxm <- prior$maxm
nu <- c(prior$nu,prior$nub)
if(is.null(prior$a0))
{
a0b0 <- c(-1,-1,0)
ap <- prior$a
apr <- 0
}
else
{
a0b0 <- c(prior$a0,prior$b0,prior$q)
ap <- 0
apr <- 1
if(prior$a0<0 || prior$b0<0)
{
stop("The parameters of the Beta prior for the a parameter must be possitive.\n")
}
if(prior$q<0 || prior$q>1)
{
stop("The mixture proportion must be in (0,1).\n")
}
}
if(is.null(prior$mub))
{
a0b0 <- c(a0b0,0,-1)
bp <- prior$b
bpr <- 0
}
else
{
a0b0 <- c(a0b0,prior$mub,prior$sigmab)
bpr <- 1
bp <- abs(rnorm(1,mean=prior$mub,sd=sqrt(prior$sigmab)))
}
alpha <- c(ap,bp)
tinv1 <- prior$tinv
tinv2 <- prior$tbinv
psiinv <- solve(prior$S0)
smu <- psiinv%*%prior$m0
#########################################################################################
# mcmc specification
#########################################################################################
mcmcvec <- c(mcmc$nburn,mcmc$nskip,mcmc$ndisplay)
nsave <- mcmc$nsave
if(is.null(mcmc$tune1))
{
tune1 <- 1.1
}
else
{
tune1 <- mcmc$tune1
}
if(is.null(mcmc$tune2))
{
tune2 <- 1.1
}
else
{
tune2 <- mcmc$tune2
}
a0b0 <- c(a0b0,tune1,tune2)
#########################################################################################
# Starting values
#########################################################################################
sigma <- matrix(0,nrow=2*nvar,ncol=2*nvar)
beta <- NULL
sigma <- NULL
for(i in 1:nvar)
{
startp <- (i-1)*p+1
endp<- (i-1)*p+p
fit0 <- glm.fit(xonset[,(startp:endp)], y[,i], family=gaussian())
betatmp <- coefficients(fit0)
beta <- c(beta,betatmp)
s2 <- var(residuals(fit0))
sigma <- c(sigma,s2)
}
for(i in 1:nvar)
{
startp <- (i-1)*q+1
endp <- (i-1)*q+q
fit0 <- glm.fit(xfailure[,(startp:endp)], y[,nvar+i], family=gaussian())
betatmp <- coefficients(fit0)
beta <- c(beta,betatmp)
s2 <- var(residuals(fit0))
sigma <- c(sigma,s2)
}
sigma <- diag(sigma,2*nvar)
mub <- beta
sigmab <- diag(1,nvar*(p+q))
#########################################################################################
# output
#########################################################################################
acrate <- rep(0,2)
npar <- (2*nvar*(2*nvar+1)/2)+nvar*(p+q)+(nvar*(p+q)*(nvar*(p+q)+1)/2)+3
f <- matrix(0,nrow=npred*2*nvar,ncol=ngrid)
h <- matrix(0,nrow=npred*2*nvar,ncol=ngrid)
thetasave <- matrix(0,nrow=nsave,ncol=npar)
randsave <- matrix(0,nrow=nsave,ncol=npred*2*nvar)
#########################################################################################
# parameters depending on status
#########################################################################################
if(status==TRUE)
{
sigmainv <- solve(sigma)
sigmabinv <- solve(sigmab)
b <- matrix(0,nrow=maxm,ncol=nvar*(p+q))
for(i in 1:1)
{
b[i,] <- beta
}
for(i in 2:maxm)
{
b[i,] <- rep(0,nvar*(p+q))
}
ss <- rep(1,nsubject)
ncluster <- 1
}
if(status==FALSE)
{
alpha <- state$alpha
b <- state$b
sigma <- state$sigma
sigmainv <- solve(state$sigma)
mub <- state$mb
sigmab <- state$Sb
sigmabinv <- solve(state$Sb)
y <- state$y
ss <- state$ss
ncluster <- state$ncluster
}
#########################################################################################
# working space
#########################################################################################
cstrt <- matrix(0,nrow=maxm,ncol=nsubject)
ccluster <- rep(0,maxm)
prob <- rep(0,maxm)
prob2 <- rep(0,maxm)
weights <- rep(0,maxm)
iflagc <- rep(0,nvar*(p+q))
theta <- rep(0,nvar*(p+q))
workmc <- matrix(0,nrow=nvar*(p+q),ncol=nvar*(p+q))
workmc2 <- matrix(0,nrow=nvar*(p+q),ncol=nvar*(p+q))
workmc3 <- matrix(0,nrow=nvar*(p+q),ncol=nvar*(p+q))
workmhc <- rep(0,nvar*(p+q)*(nvar*(p+q)+1)/2)
workmhc2 <- rep(0,nvar*(p+q)*(nvar*(p+q)+1)/2)
workvc <- rep(0,nvar*(p+q))
iflagn <- rep(0,2*nvar)
workmn <- matrix(0,nrow=2*nvar,ncol=2*nvar)
workmn2 <- matrix(0,nrow=2*nvar,ncol=2*nvar)
workmn3 <- matrix(0,nrow=2*nvar,ncol=2*nvar)
workmhn <- rep(0,2*nvar*(2*nvar+1)/2)
workmhn2 <- rep(0,2*nvar*(2*nvar+1)/2)
workvn <- rep(0,2*nvar)
workvn2 <- rep(0,2*nvar)
workvn3 <- rep(0,2*nvar)
ztz <- matrix(0,nrow=nvar*(p+q),ncol=nvar*(p+q))
zty <- rep(0,nvar*(p+q))
seed <- c(sample(1:29000,1),sample(1:29000,1))
model <- matrix(0,nrow=2*nvar,ncol=(p+q))
possi <- matrix(0,nrow=2*nvar,ncol=(p+q))
fw <- rep(0,ngrid)
fw2 <- rep(0,ngrid)
varind <- rep(1,p)
for(i in 2:nvar)
{
varind <- c(varind,rep(i,p))
}
for(i in 1:nvar)
{
varind <- c(varind,rep(nvar+i,q))
}
#########################################################################################
# calling the fortran code
#########################################################################################
foo <- .Fortran("ldpddoublyintsba",
nsubject =as.integer(nsubject),
nvar =as.integer(nvar),
p =as.integer(p),
q =as.integer(q),
ngrid =as.integer(ngrid),
npred =as.integer(npred),
x =as.double(x),
llower =as.double(llower),
lupper =as.double(lupper),
grid =as.double(grid),
xpred =as.double(xpred),
maxm =as.integer(maxm),
a0b0 =as.double(a0b0),
nu =as.integer(nu),
tinv1 =as.double(tinv1),
smu =as.double(smu),
psiinv =as.double(psiinv),
tinv2 =as.double(tinv2),
mcmc =as.integer(mcmcvec),
nsave =as.integer(nsave),
ncluster =as.integer(ncluster),
ss =as.integer(ss),
alpha =as.double(alpha),
b =as.double(b),
sigma =as.double(sigma),
sigmainv =as.double(sigmainv),
mub =as.double(mub),
sigmab =as.double(sigmab),
sigmabinv =as.double(sigmabinv),
y =as.double(y),
acrate =as.double(acrate),
f =as.double(f),
h =as.double(h),
thetasave =as.double(thetasave),
randsave =as.double(randsave),
seed =as.integer(seed),
model =as.integer(model),
possi =as.integer(possi),
varind =as.integer(varind),
cstrt =as.integer(cstrt),
ccluster =as.integer(ccluster),
prob =as.double(prob),
prob2 =as.double(prob2),
weights =as.double(weights),
iflagc =as.integer(iflagc),
theta =as.double(theta),
workmc =as.double(workmc),
workmc2 =as.double(workmc2),
workmc3 =as.double(workmc3),
workmhc =as.double(workmhc),
workmhc2 =as.double(workmhc2),
workvc =as.double(workvc),
iflagn =as.integer(iflagn),
workmn =as.double(workmn),
workmn2 =as.double(workmn2),
workmn3 =as.double(workmn3),
workmhn =as.double(workmhn),
workmhn2 =as.double(workmhn2),
workvn =as.double(workvn),
workvn2 =as.double(workvn2),
workvn3 =as.double(workvn3),
ztz =as.double(ztz),
zty =as.double(zty),
fw =as.double(fw),
fw2 =as.double(fw2),
PACKAGE="DPpackage")
#########################################################################################
# save state
#########################################################################################
model.name<-"Linear Dependent Poisson-Dirichlet Survival Model for Doubly-Interval Censored Data"
nametmp<-c(paste("onset",1:nvar,sep=""),
paste("time",1:nvar,sep="")
)
pnames<-NULL
for(i in 1:(2*nvar))
{
for(j in i:(2*nvar))
{
tmp<-paste(nametmp[i],nametmp[j],sep=";")
tmp<-paste("Sigma (",tmp,")",sep="")
pnames<-c(pnames,tmp)
}
}
nametmp <- colnames(x)
pnames <- c(pnames,paste("mb",nametmp,sep="-"))
for(i in 1:(nvar*(p+q)))
{
for(j in i:(nvar*(p+q)))
{
tmp<-paste(nametmp[i],nametmp[j],sep=";")
tmp<-paste("Sb (",tmp,")",sep="")
pnames<-c(pnames,tmp)
}
}
pnames <- c(pnames,"ncluster","a","b")
thetasave <- matrix(foo$thetasave,nrow=nsave,ncol=npar)
colnames(thetasave) <- pnames
coeff <- apply(thetasave,2,mean)
nametmp<-c(paste("onset",1:nvar,sep=""),
paste("time",1:nvar,sep="")
)
nametmp <- rep(nametmp,npred)
randsave <- matrix(foo$randsave,nrow=nsave,ncol=npred*2*nvar)
colnames(randsave)<-nametmp
coeffr <- apply(randsave,2,mean)
medtable <- matrix(coeffr,nrow=npred,ncol=2*nvar,byrow=T)
nametmp<-c(paste("onset",1:nvar,sep=""),
paste("time",1:nvar,sep="")
)
colnames(medtable)<-nametmp
nametmp<-paste("prediction",1:npred)
rownames(medtable)<-nametmp
save.state <- list(thetasave=thetasave,
randsave=randsave)
state <- list( alpha=foo$alpha,
b=matrix(foo$b,nrow=maxm,ncol=nvar*(p+q)),
sigma=matrix(foo$sigma,nrow=2*nvar,ncol=2*nvar),
mb=foo$mub,
Sb=matrix(foo$sigmab,nrow=nvar*(p+q),ncol=nvar*(p+q)),
y=matrix(foo$y,nrow=nsubject,ncol=2*nvar),
ss=foo$ss,
ncluster=foo$ncluster)
acrate <- foo$acrate[1]
dptest <- foo$acrate[2]
z<-list(modelname=model.name,
coefficients=coeff,
call=cl,
acrate=acrate,
dptest=dptest,
prior=prior,
mcmc=mcmc,
state=state,
save.state=save.state,
nsubject=foo$nsubject,
nsave=nsave,
p=foo$p,
q=foo$q,
nvar=2*foo$nvar,
npred=foo$npred,
ngrid=ngrid,
grid=grid,
f=matrix(foo$f,nrow=npred*2*nvar,ncol=ngrid),
h=matrix(foo$h,nrow=npred*2*nvar,ncol=ngrid),
medtable=medtable,
acrate=foo$acrate,
work.dir=work.dir)
cat("\n\n")
class(z)<-c("LDPDdoublyint")
z
}
###
### Tools for LDPDdoublyint: print, summary, plot, predict
###
### Copyright: Alejandro Jara, 2007
### Last modification: 08-13-2007.
"print.LDPDdoublyint"<-function (x, digits = max(3, getOption("digits") - 3), ...)
{
cat("\n",x$modelname,"\n\nCall:\n", sep = "")
print(x$call)
cat("\nP(DP|data)=",x$dptest,"\n")
cat("\nPosterior Inference of Parameters:\n")
print.default(format(x$coefficients, digits = digits), print.gap = 2,
quote = FALSE)
if(!is.null(x$acrate))
{
cat("\nAcceptance Rate for Metropolis Steps = ",x$acrate)
}
cat("\n\nNumber of Subjects:",x$nsubject)
cat("\nNumber of Variables:",x$nvar)
cat("\n\n")
invisible(x)
}
"summary.LDPDdoublyint"<-function(object, hpd=TRUE, ...)
{
stde<-function(x)
{
n<-length(x)
return(sd(x)/sqrt(n))
}
hpdf<-function(x)
{
alpha<-0.05
vec<-x
n<-length(x)
alow<-rep(0,2)
aupp<-rep(0,2)
a<-.Fortran("hpd",n=as.integer(n),alpha=as.double(alpha),x=as.double(vec),
alow=as.double(alow),aupp=as.double(aupp),PACKAGE="DPpackage")
return(c(a$alow[1],a$aupp[1]))
}
pdf<-function(x)
{
alpha<-0.05
vec<-x
n<-length(x)
alow<-rep(0,2)
aupp<-rep(0,2)
a<-.Fortran("hpd",n=as.integer(n),alpha=as.double(alpha),x=as.double(vec),
alow=as.double(alow),aupp=as.double(aupp),PACKAGE="DPpackage")
return(c(a$alow[2],a$aupp[2]))
}
thetasave<-object$save.state$thetasave
ans <- c(object[c("call", "modelname")])
### Kernel
nvar <- object$nvar
dimen1 <- nvar*(nvar+1)/2
coef.p <- object$coefficients[1:dimen1]
mat <- thetasave[,1:dimen1]
if(dimen1>0)
{
coef.m <-apply(mat, 2, median)
coef.sd<-apply(mat, 2, sd)
coef.se<-apply(mat, 2, stde)
if(hpd){
limm<-apply(mat, 2, hpdf)
coef.l<-limm[1,]
coef.u<-limm[2,]
}
else
{
limm<-apply(mat, 2, pdf)
coef.l<-limm[1,]
coef.u<-limm[2,]
}
coef.table <- cbind(coef.p, coef.m, coef.sd, coef.se , coef.l , coef.u)
if(hpd)
{
dimnames(coef.table) <- list(names(coef.p), c("Mean", "Median", "Std. Dev.", "Naive Std.Error",
"95%HPD-Low","95%HPD-Upp"))
}
else
{
dimnames(coef.table) <- list(names(coef.p), c("Mean", "Median", "Std. Dev.", "Naive Std.Error",
"95%CI-Low","95%CI-Upp"))
}
ans$kernel <- coef.table
}
### Baseline
p <- object$p
q <- object$q
nt <- (p+q)*(nvar/2)
dimen2 <- nt + nt*(nt+1)/2
coef.p <- object$coefficients[(dimen1+1):(dimen1+dimen2)]
mat <- thetasave[,(dimen1+1):(dimen1+dimen2)]
if(dimen2>0)
{
coef.m <-apply(mat, 2, median)
coef.sd<-apply(mat, 2, sd)
coef.se<-apply(mat, 2, stde)
if(hpd){
limm<-apply(mat, 2, hpdf)
coef.l<-limm[1,]
coef.u<-limm[2,]
}
else
{
limm<-apply(mat, 2, pdf)
coef.l<-limm[1,]
coef.u<-limm[2,]
}
coef.table <- cbind(coef.p, coef.m, coef.sd, coef.se , coef.l , coef.u)
if(hpd)
{
dimnames(coef.table) <- list(names(coef.p), c("Mean", "Median", "Std. Dev.", "Naive Std.Error",
"95%HPD-Low","95%HPD-Upp"))
}
else
{
dimnames(coef.table) <- list(names(coef.p), c("Mean", "Median", "Std. Dev.", "Naive Std.Error",
"95%CI-Low","95%CI-Upp"))
}
ans$base <- coef.table
}
### Precision
dimen3 <- 3
coef.p <- object$coefficients[(dimen1+dimen2+1):(dimen1+dimen2+dimen3)]
mat <- thetasave[,(dimen1+dimen2+1):(dimen1+dimen2+dimen3)]
if(dimen3>0)
{
coef.m <-apply(mat, 2, median)
coef.sd<-apply(mat, 2, sd)
coef.se<-apply(mat, 2, stde)
if(hpd){
limm<-apply(mat, 2, hpdf)
coef.l<-limm[1,]
coef.u<-limm[2,]
}
else
{
limm<-apply(mat, 2, pdf)
coef.l<-limm[1,]
coef.u<-limm[2,]
}
coef.table <- cbind(coef.p, coef.m, coef.sd, coef.se , coef.l , coef.u)
if(hpd)
{
dimnames(coef.table) <- list(names(coef.p), c("Mean", "Median", "Std. Dev.", "Naive Std.Error",
"95%HPD-Low","95%HPD-Upp"))
}
else
{
dimnames(coef.table) <- list(names(coef.p), c("Mean", "Median", "Std. Dev.", "Naive Std.Error",
"95%CI-Low","95%CI-Upp"))
}
ans$prec <- coef.table
}
ans$nsubject <- object$nsubject
ans$nvar <- object$nvar
ans$acrate <- object$acrate
ans$dptest <- object$dptest
class(ans) <- "summaryLDPDdoublyint"
return(ans)
}
"print.summaryLDPDdoublyint"<-function (x, digits = max(3, getOption("digits") - 3), ...)
{
cat("\n",x$modelname,"\n\nCall:\n", sep = "")
print(x$call)
cat("\nP(DP|data)=",x$dptest,"\n")
cat("\nKernel covariance matrix:\n")
print.default(format(x$kernel, digits = digits), print.gap = 2,
quote = FALSE)
cat("\nCentering parameters:\n")
print.default(format(x$base, digits = digits), print.gap = 2,
quote = FALSE)
cat("\nPrecision parameters:\n")
print.default(format(x$prec, digits = digits), print.gap = 2,
quote = FALSE)
if(!is.null(x$acrate))
{
cat("\nAcceptance Rate for Metropolis Steps = ",x$acrate)
}
cat("\n\nNumber of Subject:",x$nsubject)
cat("\nNumber of Variable:",x$nvar)
cat("\n\n")
invisible(x)
}
"plot.LDPDdoublyint"<-function(x, hpd=TRUE, ask=TRUE, nfigr=2, nfigc=2, param=NULL, col="#bdfcc9", ...)
{
fancydensplot<-function(x, hpd=TRUE, npts=200, xlab="", ylab="", main="",col="#bdfcc9", ...)
# Author: AJV, 2006
#
{
dens <- density(x,n=npts)
densx <- dens$x
densy <- dens$y
meanvar <- mean(x)
densx1 <- max(densx[densx<=meanvar])
densx2 <- min(densx[densx>=meanvar])
densy1 <- densy[densx==densx1]
densy2 <- densy[densx==densx2]
ymean <- densy1 + ((densy2-densy1)/(densx2-densx1))*(meanvar-densx1)
if(hpd==TRUE)
{
alpha<-0.05
alow<-rep(0,2)
aupp<-rep(0,2)
n<-length(x)
a<-.Fortran("hpd",n=as.integer(n),alpha=as.double(alpha),x=as.double(x),
alow=as.double(alow),aupp=as.double(aupp),PACKAGE="DPpackage")
xlinf<-a$alow[1]
xlsup<-a$aupp[1]
}
else
{
xlinf <- quantile(x,0.025)
xlsup <- quantile(x,0.975)
}
densx1 <- max(densx[densx<=xlinf])
densx2 <- min(densx[densx>=xlinf])
densy1 <- densy[densx==densx1]
densy2 <- densy[densx==densx2]
ylinf <- densy1 + ((densy2-densy1)/(densx2-densx1))*(xlinf-densx1)
densx1 <- max(densx[densx<=xlsup])
densx2 <- min(densx[densx>=xlsup])
densy1 <- densy[densx==densx1]
densy2 <- densy[densx==densx2]
ylsup <- densy1 + ((densy2-densy1)/(densx2-densx1))*(xlsup-densx1)
plot(0.,0.,xlim = c(min(densx), max(densx)), ylim = c(min(densy), max(densy)),
axes = F,type = "n" , xlab=xlab, ylab=ylab, main=main, cex=1.2)
xpol<-c(xlinf,xlinf,densx[densx>=xlinf & densx <=xlsup],xlsup,xlsup)
ypol<-c(0,ylinf,densy[densx>=xlinf & densx <=xlsup] ,ylsup,0)
polygon(xpol, ypol, border = FALSE,col=col)
lines(c(min(densx), max(densx)),c(0,0),lwd=1.2)
segments(min(densx),0, min(densx),max(densy),lwd=1.2)
lines(densx,densy,lwd=1.2)
segments(meanvar, 0, meanvar, ymean,lwd=1.2)
segments(xlinf, 0, xlinf, ylinf,lwd=1.2)
segments(xlsup, 0, xlsup, ylsup,lwd=1.2)
axis(1., at = round(c(xlinf, meanvar,xlsup), 2.), labels = T,pos = 0.)
axis(1., at = round(seq(min(densx),max(densx),length=15), 2.), labels = F,pos = 0.)
axis(2., at = round(seq(0,max(densy),length=5), 2.), labels = T,pos =min(densx))
}
if(is(x, "LDPDdoublyint"))
{
if(is.null(param))
{
coef.p<-x$coefficients
n<-length(coef.p)
pnames<-names(coef.p)
par(ask = ask)
layout(matrix(seq(1,nfigr*nfigc,1), nrow=nfigr, ncol=nfigc, byrow = TRUE))
for(i in 1:(n-2)){
title1<-paste("Trace of",pnames[i],sep=" ")
title2<-paste("Density of",pnames[i],sep=" ")
plot(x$save.state$thetasave[,i],type='l',main=title1,xlab="MCMC scan",ylab=" ")
if(pnames[i]=="ncluster" || pnames[i]=="k")
{
hist(x$save.state$thetasave[,i],main=title2,xlab="values", ylab="probability",probability=TRUE)
}
else
{
fancydensplot(x$save.state$thetasave[,i],hpd=hpd,main=title2,xlab="values", ylab="density",col=col)
}
}
if(!is.null(x$prior$a0))
{
title1<-paste("Trace of",pnames[n-1],sep=" ")
title2<-paste("Density of",pnames[n-1],sep=" ")
plot(x$save.state$thetasave[,n-1],type='l',main=title1,xlab="MCMC scan",ylab=" ")
fancydensplot(x$save.state$thetasave[,n-1],hpd=hpd,main=title2,xlab="values", ylab="density",col=col)
}
if(!is.null(x$prior$mub))
{
title1<-paste("Trace of",pnames[n],sep=" ")
title2<-paste("Density of",pnames[n],sep=" ")
plot(x$save.state$thetasave[,n],type='l',main=title1,xlab="MCMC scan",ylab=" ")
fancydensplot(x$save.state$thetasave[,n],hpd=hpd,main=title2,xlab="values", ylab="density",col=col)
}
}
else
{
coef.p <- x$coefficients
n <- length(coef.p)
pnames<-names(coef.p)
poss<-0
for(i in 1:n)
{
if(pnames[i]==param)poss=i
}
if (poss==0 && param !="predictive")
{
stop("This parameter is not present in the original model.\n")
}
if (param !="ncluster")
{
par(ask = ask)
layout(matrix(seq(1,nfigr*nfigc,1), nrow=nfigr, ncol=nfigc, byrow = TRUE))
title1<-paste("Trace of",pnames[poss],sep=" ")
title2<-paste("Density of",pnames[poss],sep=" ")
plot(x$save.state$thetasave[,poss],type='l',main=title1,xlab="MCMC scan",ylab=" ")
fancydensplot(x$save.state$thetasave[,poss],hpd=hpd,main=title2,xlab="values", ylab="density",col=col)
}
else
{
par(ask = ask)
layout(matrix(seq(1,nfigr*nfigc,1), nrow=nfigr, ncol=nfigc, byrow = TRUE))
title1<-paste("Trace of",pnames[poss],sep=" ")
title2<-paste("Density of",pnames[poss],sep=" ")
plot(x$save.state$thetasave[,poss],type='l',main=title1,xlab="MCMC scan",ylab=" ")
hist(x$save.state$thetasave[,poss],main=title2,xlab="values", ylab="probability",probability=TRUE)
}
}
}
}
"predict.LDPDdoublyint"<-
function(object,grid,compute.band=FALSE,hpd=TRUE, ...)
{
hpdf<-function(x)
{
alpha<-0.05
vec<-x
n<-length(x)
alow<-rep(0,2)
aupp<-rep(0,2)
a<-.Fortran("hpd",n=as.integer(n),alpha=as.double(alpha),x=as.double(vec),
alow=as.double(alow),aupp=as.double(aupp),PACKAGE="DPpackage")
return(c(a$alow[1],a$aupp[1]))
}
pdf<-function(x)
{
alpha<-0.05
vec<-x
n<-length(x)
alow<-rep(0,2)
aupp<-rep(0,2)
a<-.Fortran("hpd",n=as.integer(n),alpha=as.double(alpha),x=as.double(vec),
alow=as.double(alow),aupp=as.double(aupp),PACKAGE="DPpackage")
return(c(a$alow[2],a$aupp[2]))
}
if(is(object, "LDPDdoublyint"))
{
ngrid <- object$ngrid
grid <- object$grid
nvar <- object$nvar/2
npred <- object$npred
f <- object$f
h <- object$h
med <- object$medtable
if(compute.band)
{
nsave <- object$nsave
alpha <- 0.05
tint <- 0
if(hpd)tint <- 1
llower <- NULL
lupper <- NULL
llower <- matrix(0,nrow=npred*2*nvar,ncol=ngrid)
lupper <- matrix(0,nrow=npred*2*nvar,ncol=ngrid)
llower2 <- matrix(0,nrow=npred*2*nvar,ncol=ngrid)
lupper2 <- matrix(0,nrow=npred*2*nvar,ncol=ngrid)
workv1 <- rep(0,nsave)
workv2 <- rep(0,ngrid)
foo <- .Fortran("hpdspy",
nsave =as.integer(nsave),
npred =as.integer(npred),
nvar =as.integer(nvar),
ngrid =as.integer(ngrid),
alpha =as.double(alpha),
tint =as.integer(tint),
workv1 =as.double(workv1),
workv2 =as.double(workv2),
llower =as.double(llower),
lupper =as.double(lupper),
llower2 =as.double(llower2),
lupper2 =as.double(lupper2),
PACKAGE="DPpackage")
llower <- matrix(foo$llower,nrow=npred*2*nvar,ncol=ngrid)
lupper <- matrix(foo$lupper,nrow=npred*2*nvar,ncol=ngrid)
llower2 <- matrix(foo$llower2,nrow=npred*2*nvar,ncol=ngrid)
lupper2 <- matrix(foo$lupper2,nrow=npred*2*nvar,ncol=ngrid)
if(hpd)
{
limm <- apply(object$save.state$randsave, 2, hpdf)
coef.l <- limm[1,]
coef.u <- limm[2,]
}
else
{
limm <- apply(object$save.state$randsave, 2, pdf)
coef.l <- limm[1,]
coef.u <- limm[2,]
}
med.l <- matrix(coef.l,nrow=npred,ncol=2*nvar,byrow=T)
med.u <- matrix(coef.u,nrow=npred,ncol=2*nvar,byrow=T)
colnames(med.l) <- colnames(med)
colnames(med.u) <- colnames(med)
rownames(med.l) <- rownames(med)
rownames(med.u) <- rownames(med)
}
else
{
llower <- NULL
lupper <- NULL
llower2 <- NULL
lupper2 <- NULL
med.l <- NULL
med.u <- NULL
}
out <- list(f=f,
h=h,
med=med,
med.l=med.l,
med.u=med.u,
grid=grid,
ngrid=ngrid,
nvar=nvar,
npred=npred,
llower=llower,
lupper=lupper,
llower2=llower2,
lupper2=lupper2,
modelname=object$modelname,
call=object$call,
compute.band=compute.band)
class(out) <- c("predict.LDPDdoublyint")
out
}
}
"print.predict.LDPDdoublyint"<-function (x, digits = max(3, getOption("digits") - 3), ...)
{
cat("\n",x$modelname,"\n\nCall:\n", sep = "")
print(x$call)
cat("\nPosterior Mean of Median Survival for Predictions:\n")
print.default(format(x$med, digits = digits), print.gap = 2,
quote = FALSE)
if(!is.null(x$med.l))
{
cat("\nLower CI Limit for the Median Survival for Predictions:\n")
print.default(format(x$med.l, digits = digits), print.gap = 2,
quote = FALSE)
}
if(!is.null(x$med.u))
{
cat("\nUpper CI Limit for the Median Survival for Predictions:\n")
print.default(format(x$med.u, digits = digits), print.gap = 2,
quote = FALSE)
}
cat("\n\n")
invisible(x)
}
plot.predict.LDPDdoublyint<-function(x,ask=TRUE,xlim=NULL,nfigr=1,nfigc=1, ...)
{
if(is(x, "predict.LDPDdoublyint"))
{
par(ask = ask)
layout(matrix(seq(1,nfigr*nfigc,1), nrow=nfigr, ncol=nfigc, byrow = TRUE))
compute.band <- x$compute.band
if(compute.band)
{
k <- 0
for(i in 1:x$npred)
{
for(j in 1:x$nvar)
{
tmp1 <- paste("Prediction",i,sep=":")
tmp2 <- paste("Onset",j,sep=":")
tmp3 <- paste(tmp1,tmp2,sep="-")
k <- k+1
x1 <- x$grid[j,]
y1 <- x$llower[k,]
x2 <- x$grid[j,]
y2 <- x$lupper[k,]
aa <- rbind(x2,y2)[, order(-x2, y2)]
x2 <- aa[1,]
y2 <- aa[2,]
plot(x$grid[j,],x$f[k,],type="l",lty=1,lwd=3,ylim=c(0,1),
ylab="Survival",xlab="Time",main=tmp3,bty="n")
polygon(x=c(x1,x2),y=c(y1,y2),border=NA,col="lightgray")
lines(x$grid[j,],x$f[k,],lwd=3,lty=1)
}
for(j in 1:x$nvar)
{
tmp1 <- paste("Prediction",i,sep=":")
tmp2 <- paste("Failure",j,sep=":")
tmp3 <- paste(tmp1,tmp2,sep="-")
k <- k+1
x1 <- x$grid[j,]
y1 <- x$llower[k,]
x2 <- x$grid[j,]
y2 <- x$lupper[k,]
aa <- rbind(x2,y2)[, order(-x2, y2)]
x2 <- aa[1,]
y2 <- aa[2,]
plot(x$grid[j,],x$f[k,],type="l",lty=1,lwd=3,ylim=c(0,1),
ylab="Survival",xlab="Time",main=tmp3,bty="n")
polygon(x=c(x1,x2),y=c(y1,y2),border=NA,col="lightgray")
lines(x$grid[j,],x$f[k,],lwd=3,lty=1)
}
}
k <- 0
for(i in 1:x$npred)
{
for(j in 1:x$nvar)
{
tmp1 <- paste("Prediction",i,sep=":")
tmp2 <- paste("Onset",j,sep=":")
tmp3 <- paste(tmp1,tmp2,sep="-")
k <- k+1
x1 <- grid[j,]
y1 <- x$llower2[k,]
x2 <- grid[j,]
y2 <- x$lupper2[k,]
aa <- rbind(x2,y2)[, order(-x2, y2)]
x2 <- aa[1,]
y2 <- aa[2,]
plot(x$grid[j,],x$lupper2[k,],type="l",lty=1,lwd=2,
ylab="Hazard",xlab="Time",main=tmp3,bty="n",col="lightgray")
polygon(x=c(x1,x2),y=c(y1,y2),border=NA,col="lightgray")
lines(x$grid[j,],x$h[k,],lwd=3,lty=1)
}
for(j in 1:x$nvar)
{
tmp1 <- paste("Prediction",i,sep=":")
tmp2 <- paste("Failure",j,sep=":")
tmp3 <- paste(tmp1,tmp2,sep="-")
k <- k+1
x1 <- grid[j,]
y1 <- x$llower2[k,]
x2 <- grid[j,]
y2 <- x$lupper2[k,]
aa <- rbind(x2,y2)[, order(-x2, y2)]
x2 <- aa[1,]
y2 <- aa[2,]
plot(x$grid[j,],x$lupper2[k,],type="l",lty=1,lwd=2,
ylab="Hazard",xlab="Time",main=tmp3,bty="n",col="lightgray")
polygon(x=c(x1,x2),y=c(y1,y2),border=NA,col="lightgray")
lines(x$grid[j,],x$h[k,],lwd=3,lty=1)
}
}
}
else
{
k <- 0
for(i in 1:x$npred)
{
for(j in 1:x$nvar)
{
tmp1 <- paste("Prediction",i,sep=":")
tmp2 <- paste("Onset",j,sep=":")
tmp3 <- paste(tmp1,tmp2,sep="-")
k <- k+1
plot(x$grid[j,],x$f[k,],type="l",lty=1,lwd=3,ylim=c(0,1),
ylab="Survival",xlab="Time",main=tmp3,bty="n")
}
for(j in 1:x$nvar)
{
tmp1 <- paste("Prediction",i,sep=":")
tmp2 <- paste("Failure",j,sep=":")
tmp3 <- paste(tmp1,tmp2,sep="-")
k <- k+1
plot(x$grid[j,],x$f[k,],type="l",lty=1,lwd=3,ylim=c(0,1),
ylab="Survival",xlab="Time",main=tmp3,bty="n")
}
}
k <- 0
for(i in 1:x$npred)
{
for(j in 1:x$nvar)
{
tmp1 <- paste("Prediction",i,sep=":")
tmp2 <- paste("Onset",j,sep=":")
tmp3 <- paste(tmp1,tmp2,sep="-")
k <- k+1
plot(x$grid[j,],x$h[k,],type="l",lty=1,lwd=3,
ylab="Hazard",xlab="Time",main=tmp3,bty="n")
}
for(j in 1:x$nvar)
{
tmp1 <- paste("Prediction",i,sep=":")
tmp2 <- paste("Failure",j,sep=":")
tmp3 <- paste(tmp1,tmp2,sep="-")
k <- k+1
plot(x$grid[j,],x$h[k,],type="l",lty=1,lwd=3,
ylab="Hazard",xlab="Time",main=tmp3,bty="n")
}
}
}
}
}
Try the DPpackage package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
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.