Skip to content
Permalink
Browse files

Initial commit

  • Loading branch information...
wangxj03 committed Jul 24, 2012
0 parents commit eec16a03d5cbca89674a6532a3aa7880fff401d3
676 COPYING

Large diffs are not rendered by default.

@@ -0,0 +1,18 @@
2012-03-14 Jun Yan <jun.yan@uconn.edu>

* Removed tiCox.h and tiCox.cpp. Functions there are covered
by bayesCox.

* Cleaned the style of the code under src with astyle -s2.

2012-03-10 Jun Yan <jun.yan@uconn.edu>

* Added system dependence on the boost library, which allows
removal of the boost headers files from the src directory.
The configure.ac file was adapted from that in package
RQuantLib; thank Dirk Eddelbuettel <edd@debian.org>.

* Removed std:cout lines (Xiaojing Wang).

* Replaced exit with return whenever possible. May need
more fix later.
@@ -0,0 +1,24 @@
Package: dynsurv
Version: 0.1-9
Date: 2012-07-22
Title: Dynamic models for survival data
Author: Xiaojing Wang <wangxj03@gmail.com>, Jun Yan
<jun.yan@uconn.edu>, and Ming-Hui Chen
<ming-hui.chen@uconn.edu>
Maintainer: Jun Yan <jun.yan@uconn.edu>
Depends: methods, survival, plyr, reshape, ggplot2, nleqslv
Description: Functions to fit time-varying coefficient models for
interval censored and right censored survival data. Three major
approaches are implemented: 1) Bayesian Cox model with
time-independent, time-varying or dynamic coefficients for
right censored and interval censored data; 2) Spline based
time-varying coefficient Cox model for right censored data; 3)
Transformation model with time-varying coefficients for right
censored data using estimating equations.
SystemRequirements: Boost library from http://www.boost.org
License: GPL (>= 3)
LazyLoad: Yes
LazyData: Yes
Packaged: 2012-05-23 17:03:18 UTC; jyan
Repository: CRAN
Date/Publication: 2012-05-24 06:24:14
55 MD5
@@ -0,0 +1,55 @@
f27defe1e96c2e1ecd4e0c9be8967949 *COPYING
e78327a622b68f98203a828ca2c4e031 *ChangeLog
65af30b9f009b4c47d4d55a412356cec *DESCRIPTION
b6bff42df1ad890b24ae80a2a545a51f *NAMESPACE
40833cef7d5dccd6c7ea51a79d68ddb7 *R/bayesCox.R
5a1d9158c1b0b88a55f45698a77ff71a *R/coef.R
9814ea58fa47eedacdfaeeea36a95669 *R/jump.R
47ba19846805453234e992d72f04b14c *R/plot.R
cd6ba0cd5b20db5b527d4edaab8613c9 *R/print.R
8f209f3b568ea114225c8732f8eb0eb9 *R/splineCox.R
03e386388c0fed9387b8d262ee6ce02d *R/tvTran.R
e4e3ec97e4dcbf9fca84e89ad8d5a4dc *cleanup
ec25ae93a50a125394124952475e9814 *configure
0e45691eda4e9173d00e39736fd7dec1 *configure.ac
7b2ed07bd37790cccee322eb4995406c *data/bcos.RData
f9dd9273ba88f15f2e212960c52ee91e *data/tooth.RData
dfddbbe91232148c76aa54f5c8610cd7 *man/bayesCox.Rd
88c092b2c56aa7a52155aeef423f93fc *man/bcos.Rd
f58b16042dbe007fe2aea5cdc2ea1278 *man/coef.bayesCox.Rd
7b5b7a044986e71a6059aae10ef360a2 *man/coef.splineCox.Rd
079c658157fe7418e70bf61715491294 *man/coef.tvTran.Rd
d6d0641bf5e1137e98bde58e8d35007d *man/jump.bayesCox.Rd
427e2f9da930a17747b7fe4bdbd912ae *man/nu.bayesCox.Rd
4707d4eb86718cff48f1acb555a57bd2 *man/plotCoef.Rd
6e84072606c479673000f2db0780c08e *man/plotJump.Rd
ddaf2365b830504d0e42dc14ff787fac *man/plotNu.Rd
5ec038258d3a38256c72cc0d3f52c089 *man/splineCox.Rd
2cabb5de28e590184b99ab344a7bc5d5 *man/tooth.Rd
9a37276ae6c1e153fa360fdc575d37de *man/tvTran.Rd
adace84d992dc754f1cecba279b08dfe *src/CoxModel.h
509c7d80d27562a2fcf5a329624b3820 *src/DynamicCoxModel.h
b0c53ea6b5b116db37ede16ad1b54973 *src/DynamicCoxModel_v2.h
dd4fcafd81083973c5d75dc7df8354d3 *src/DynamicGORHModel.h
8111eb602e88f2714fd94afb2de59e44 *src/DynamicModel.h
62f66b11c2627247206dde67d2d33d67 *src/GORHModel.h
7ca840acba7f5724b7801743e59af35e *src/GibbsSampler.h
0c9ee64fbf0d93bad1d0208b0b78ee62 *src/IntRegData.h
9e11d2b7d51f820fde4bf16484d67da4 *src/IntRegModel.h
3f45183e92ae7e4fdbd0516398aa1c8f *src/IntRegPar.cpp
5e6d3a37d5b48ea463d0ee9859313d36 *src/IntRegPar.h
98597848e8978eb7bc64765d4963d229 *src/IntRegPrior.h
4fabb18f7f60ea8c000684ddfb56d930 *src/Makevars.in
7f3a6be25a86974506a70cbe2161ffed *src/Makevars.win
3c43c1d75af0fa26c14a8ddcba72425b *src/TimeIndepCoxModel.h
180794ed91db7b270b05eba43f33a503 *src/TimeIndepGORHModel.h
a80a0073f930db97b6f1e4a291ba7ccd *src/TimeIndepModel.h
a30c2843f1719298c13d6603bdeda35a *src/TimeVaryingCoxModel.h
0c72193653c7cc0f2a2378dd1ea8d70a *src/TimeVaryingGORHModel.h
61cd119256529f2f279865ab40192715 *src/TimeVaryingModel.h
edd52447fb26a4b5b00ce8f085f275a3 *src/arms.c
cf3d07b92526c9e237ae072113120b91 *src/arms.h
8271952c6a683727830e6f0b3f2d2497 *src/bayesCox.cpp
ba1d1c77f25f112a66a639eb1958553f *src/rng.h
e0c72a43f3f73f7b08bfe48bd33cef81 *src/ublas.h
895ca902f5622bef0e0e99202d03b7dd *src/ublas_ext.h
@@ -0,0 +1,16 @@
useDynLib(dynsurv)
export(bayesCox, coef.bayesCox, jump, nu, jump.bayesCox, nu.bayesCox,
splineCox, coef.splineCox,
tvTran, coef.tvTran,
plotCoef, plotJumpTrace, plotJumpHist, plotNu)

S3method(print, bayesCox)
S3method(coef, bayesCox)
S3method(jump, bayesCox)
S3method(nu, bayesCox)

S3method(print, splineCox)
S3method(coef, splineCox)

S3method(print, tvTran)
S3method(coef, tvTran)
@@ -0,0 +1,249 @@
##############################################################################
##
## R package dynsurv by Xiaojing Wang, Jun Yan, and Ming-Hui Chen
## Copyright (C) 2011
##
## This file is part of the R package dynsurv.
##
## The R package dynsurv 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 3 of the License, or
## (at your option) any later version.
##
## The R package dynsurv 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 the R package dynsurv. If not, see <http://www.gnu.org/licenses/>.
##
##############################################################################

##############################################################################
# Baseline prior
##############################################################################
Gamma_fun <- function(shape=0.1, rate=0.1) {
list(shape=shape, rate=rate)
}

Const_fun <- function(value=1) {
list(value=value, value=value)
}

GammaProcess_fun <- function(mean=0.1, ctrl=1) {
list(mean=mean, ctrl=ctrl)
}

bp_fun <- function(type=c("Gamma", "Const", "GammaProcess"), ...) {
type <- match.arg(type)

if (type == "Gamma")
hyper <- Gamma_fun(...)
if (type == "Const")
hyper <- Const_fun(...)
if (type == "GammaProcess")
hyper <- GammaProcess_fun(...)

list(type=type, hyper=hyper)
}

##############################################################################
# Coefficient prior
##############################################################################
Normal_fun <- function(mean=0, sd=1) {
list(mean=mean, sd=sd)
}
AR1_fun <- function(sd=1) {
list(sd=sd, sd=sd)
}
HAR1_fun <- function(shape=2, scale=1) {
list(shape=shape, scale=scale)
}

cp_fun <- function(type=c("Normal", "AR1", "HAR1"), ...) {
type <- match.arg(type)

if (type == "Normal")
hyper <- Normal_fun(...)
if (type == "AR1")
hyper <- AR1_fun(...)
if (type == "HAR1")
hyper <- HAR1_fun(...)

list(type=type, hyper=hyper)
}

##############################################################################
# Gibbs sampler control and general control
##############################################################################
gibbs_fun <- function(iter=3000, burn=500, thin=1, verbose=TRUE, nReport=100) {
list(iter=iter, burn=burn, thin=thin, verbose=verbose, nReport=nReport)
}

control_bfun <- function(intercept=FALSE, a0=100, eps0=1) {
list(intercept=intercept, a0=a0, eps0=eps0)
}

##############################################################################
# Bayesian Cox model
##############################################################################
# grid: must be sorted with last number be finite
# base.prior:
# list(type="Gamma", shape=0.1, rate=0.1)
# coef.prior:
# list(type="Normal", mean=0, sd=1)
# list(type="AR1", sd=1)
# list(type="HAR1", shape=2, scale=1)

bayesCox <- function(formula, data, grid, out,
model=c("TimeIndep", "TimeVarying", "Dynamic"),
base.prior=list(), coef.prior=list(),
gibbs=list(), control=list()) {

Call <- match.call()

# Prepare prior information
model <- match.arg(model)
base.prior <- do.call("bp_fun", base.prior)
coef.prior <- do.call("cp_fun", coef.prior)

if (model == "TimeIndep") {
if (base.prior$type == "Gamma" && coef.prior$type == "Normal")
id <- 11
if (base.prior$type == "GammaProcess" && coef.prior$type == "Normal")
id <- 12
}
if (model == "TimeVarying") {
if (base.prior$type == "Gamma" && coef.prior$type == "AR1")
id <- 21
if (base.prior$type == "Gamma" && coef.prior$type == "HAR1")
id <- 22
if (base.prior$type == "GammaProcess" && coef.prior$type == "AR1")
id <- 23
if (base.prior$type == "GammaProcess" && coef.prior$type == "HAR1")
id <- 24

}
if (model == "Dynamic") {
if (base.prior$type == "Gamma" && coef.prior$type == "AR1")
id <- 31
if (base.prior$type == "Gamma" && coef.prior$type == "HAR1")
id <- 32
if (base.prior$type == "Const" && coef.prior$type == "AR1")
id <- 33
if (base.prior$type == "Const" && coef.prior$type == "HAR1")
id <- 34
if (base.prior$type == "GammaProcess" && coef.prior$type == "AR1")
id <- 35
if (base.prior$type == "GammaProcess" && coef.prior$type == "HAR1")
id <- 36
}

gibbs <- do.call("gibbs_fun", gibbs)
control <- do.call("control_bfun", control)

# Prepare data matrix LRX
mf <- model.frame(formula, data)
mm <- model.matrix(formula, data)

LRX <- cbind(mf[, 1][, 1:2], mm[, -1])
obsInd <- which(mf[, 1][, 3] == 1)
LRX[obsInd, 2] <- LRX[obsInd, 1]

cov.names <- colnames(mm)[-1]

LRX[LRX[, 2] == Inf, 2] <- max(tail(grid, 1), 999)
LRX[is.na(LRX[, 2]), 2] <- max(tail(grid, 1), 999)

if (control$intercept) {
LRX <- cbind(LRX[, 1:2], 1, LRX[, -c(1:2)])
cov.names <- c("intercept", cov.names)
}

colnames(LRX) <- c("L", "R", cov.names)

# Prepare results holder
K <- length(grid)
nBeta <- length(cov.names)
lambda <- rep(0, K)

if (model == "TimeIndep")
beta <- rep(0, nBeta)
else
beta <- rep(0, nBeta * K)

nu <- rep(0, nBeta)
jump <- rep(0, nBeta * K)

# Call C++ function
res <- .C("bayesCox",
as.double(LRX),
as.integer(nrow(LRX)),
as.integer(ncol(LRX) - 2),
as.double(grid),
as.integer(length(grid)),
as.character(out),
as.integer(id),
as.double(base.prior$hyper[[1]]),
as.double(base.prior$hyper[[2]]),
as.double(coef.prior$hyper[[1]]),
as.double(coef.prior$hyper[[2]]),
as.integer(gibbs$iter),
as.integer(gibbs$burn),
as.integer(gibbs$thin),
as.integer(gibbs$verbose),
as.integer(gibbs$nReport),
as.double(control$a0),
as.double(control$eps0),
lambda=as.double(lambda),
beta=as.double(beta),
nu=as.double(nu),
jump=as.integer(jump),
LPML=as.double(0),
DHat=as.double(0),
DBar=as.double(0),
pD=as.double(0),
DIC=as.double(0),
DIC3=as.double(0))

# Post fit processing
if (model != "TimeIndep")
res$beta <- matrix(res$beta, K, nBeta)

if (coef.prior$type != "HAR1")
res$nu <- NULL

if (model != "Dynamic")
res$jump <- NULL
else
res$jump <- matrix(res$jump, K, nBeta)

# Return list
rl <- list(call=Call,
grid=grid,
out=out,
model=model,
LRX=LRX,
base.prior=base.prior,
coef.prior=coef.prior,
gibbs=gibbs,
control=control,
N=nrow(LRX),
K=K,
nBeta=nBeta,
cov.names=cov.names,
est=list(lambda=res$lambda,
beta=res$beta,
nu=res$nu,
jump=res$jump),
measure=list(LPML=res$LPML,
DHat=res$DHat,
DBar=res$DBar,
pD=res$pD,
DIC=res$DIC,
DIC3=res$DIC3))

class(rl) <- "bayesCox"
rl
}

0 comments on commit eec16a0

Please sign in to comment.
You can’t perform that action at this time.