TABLE OF CONTENTS
00main/splinesurv.agdata [ Functions ]
NAME
splinesurv.agdata --- main estimation function
FUNCTION
This is the main fitting function for the splinesurv routine. It accepts a data frame in a specified format, conducts the initialization procedure, then either starts the MCMC loop within R or calls compiled C code that does the same thing.
See also the call graph linked from 00main.
This function should not be called directly, instead the interface provided by splinesurv.formula should be used.
SYNOPSIS
4246 splinesurv.agdata <- function(x, hazard = NULL, frailty = NULL, regression = NULL, 4247 control = NULL, coda = FALSE, initial = NULL, verbose = 3, usec = TRUE, ...)
INPUTS
x a data frame with the following columns:
i cluster index
j subject index
time event time
delta event indicator
... covariates
See the package documentation for detail on the remaining options
OUTPUTS
an object of class splinesurv, see package documentation.
SOURCE
4250 { 4251 4252 if(verbose >= 1) cat("Initializing...\n") 4253 4254 agdata <- x 4255 rm(x) 4256 call <- match.call() 4257 m <- length(unique(agdata$i)) 4258 if(m == 1) warning("Single cluster: frailty density estimate is meaningless") 4259 Ji <- table(agdata$i) 4260 4261 if(verbose >= 2) cat("\tSetting initial parameters...\n") 4262 4263 # Parse input (control) 4264 control.in <- control 4265 # default control options 4266 control.default <- list( 4267 burnin = 500, # Length of the burn - in period 4268 maxiter = 1000, # Max number of iterations 4269 thin = 1, # Degree of thinning 4270 tun.auto = TRUE, # Auto - calibrate tuning parameters 4271 tun.int = 100 # Interval for calibration of the acceptance rate 4272 ) 4273 control <- control.default 4274 controlnames <- names(control) 4275 innames <- names(control.in) 4276 if(!is.null(control.in)){ 4277 # replace default control settings by any specified in input 4278 for(n in innames) eval(parse(text = paste("control$", 4279 match.arg(n, controlnames), " <- control.in$", n, sep = ""))) 4280 } 4281 if(control$burnin > control$maxiter) stop("Burnin cannot be greater than maxiter") 4282 4283 # Parse input (frailty) 4284 frailty.in <- frailty 4285 # default settings for frailty RCurve options 4286 frailty.default <- list( 4287 type = "spline", 4288 spline.ord = 4, 4289 spline.adaptive = TRUE, 4290 spline.knots = NULL, 4291 spline.knotspacing = "equal", 4292 spline.nknots = NULL, 4293 spline.nknots.prior = NULL, 4294 spline.nknots.hyper = NULL, 4295 spline.ncandknots = 100, 4296 spline.bdmconst=.4, 4297 spline.maxoccknots = 35, 4298 spline.maxknot = 5, 4299 spline.par = NULL, 4300 spline.min=-100, 4301 spline.penalty = "none", 4302 spline.penaltyfactor = 1, 4303 spline.meanpenalty = 1e10, 4304 spline.priorvar = 0.1, 4305 spline.hyper = c(0.01, 0.01), 4306 spline.tun = 1, 4307 spline.accept = 0, 4308 param.dist = "none", 4309 param.par = NULL, 4310 param.priorvar = 0.1, 4311 param.hyper = c(0.01, 0.01), 4312 param.tun = 1, 4313 param.accept = 0, 4314 weight = 0.5, 4315 weight.priorvar = 0.1, 4316 weight.hyper = c(1, 2), 4317 weight.tun = 0.01, 4318 weight.accept = 0, 4319 accept = 0 4320 ) 4321 frailty <- frailty.default 4322 frailtynames <- names(frailty) 4323 if(!is.null(frailty.in)){ 4324 # replace default frailty options by input 4325 for(n in names(frailty.in)) eval(parse(text = paste("frailty$", 4326 match.arg(n, frailtynames), " <- frailty.in$", n, sep = ""))) 4327 } 4328 # set additional RCurve settings 4329 frailty$type <- match.arg(frailty$type, c("spline", "parametric", "both")) 4330 frailty$hasspline <- frailty$type == "spline" | frailty$type == "both" 4331 frailty$haspar <- frailty$type == "parametric" | frailty$type == "both" 4332 frailty$spline.knotspacing <- match.arg(frailty$spline.knotspacing, c("equal")) 4333 frailty$spline.penalty <- match.arg(frailty$spline.penalty, 4334 c("2diff", "2deriv", "log2deriv", "none")) 4335 frailty$param.dist <- match.arg(frailty$param.dist, c("none", "gamma", "lognormal")) 4336 if(m == 1 & frailty$haspar) stop("parametric component not allowed for single cluster") 4337 if(frailty$haspar & frailty$param.dist == "none") { 4338 warning(paste("no distribution specified for frailty parametric component", 4339 "-- setting to gamma")) 4340 frailty$param.dist <- "gamma" 4341 } 4342 # match prior settings and set default hyperparameters 4343 frailty$spline.nknots.prior <- match.arg(frailty$spline.nknots.prior, 4344 c("poisson", "geometric", "poissonmix", "negbin", "power")) 4345 if(frailty$spline.nknots.prior == "poisson" & is.null(frailty$spline.nknots.hyper)) 4346 frailty$spline.nknots.hyper <- 10 4347 if(frailty$spline.nknots.prior == "geometric" & is.null(frailty$spline.nknots.hyper)) 4348 frailty$spline.nknots.hyper <- 0.1 4349 if(frailty$spline.nknots.prior == "poissonmix" & is.null(frailty$spline.nknots.hyper)) 4350 frailty$spline.nknots.hyper <- c(10, 30) 4351 if(frailty$spline.nknots.prior == "negbin" & is.null(frailty$spline.nknots.hyper)) 4352 frailty$spline.nknots.hyper <- c(2, .1) 4353 if(frailty$spline.nknots.prior == "power" & is.null(frailty$spline.nknots.hyper)) 4354 frailty$spline.nknots.hyper <- -1 / 2 4355 4356 frailty$spline.norm <- TRUE 4357 frailty$name <- "frailty" 4358 4359 # Parse input (hazard) 4360 hazard.in <- hazard 4361 # default hazard RCurve 4362 hazard.default <- list( 4363 type = "spline", 4364 spline.ord = 4, 4365 spline.adaptive = TRUE, 4366 spline.knotspacing = "mixed", 4367 spline.nknots = NULL, 4368 spline.nknots.prior = NULL, 4369 spline.nknots.hyper = NULL, 4370 spline.ncandknots = 100, 4371 spline.maxoccknots = 35, 4372 spline.bdmconst=.4, 4373 spline.knots = NULL, 4374 spline.par = NULL, 4375 spline.min=-100, 4376 spline.penalty = "none", 4377 spline.penaltyfactor = 1, 4378 spline.priorvar = 0.1, 4379 spline.hyper = c(0.01, 0.01), 4380 spline.tun = 1, 4381 spline.accept = 0, 4382 param.dist = "none", 4383 param.par = NULL, 4384 param.priorvar = 0.1, 4385 param.hyper = c(0.01, 0.01), 4386 param.tun = 1, 4387 param.accept = 0, 4388 weight = 0.5, 4389 weight.priorvar = 0.1, 4390 weight.hyper = c(1, 2), 4391 weight.tun = 0.01, 4392 weight.accept = 0 4393 ) 4394 hazard <- hazard.default 4395 haznames <- names(hazard) 4396 if(!is.null(hazard.in)){ 4397 for(n in names(hazard.in)) eval(parse(text = paste("hazard$", 4398 match.arg(n, haznames), " <- hazard.in$", n, sep = ""))) 4399 } 4400 # other hazard settings 4401 hazard$type <- match.arg(hazard$type, c("spline", "parametric", "both")) 4402 hazard$hasspline <- hazard$type == "spline" | hazard$type == "both" 4403 hazard$haspar <- hazard$type == "parametric" | hazard$type == "both" 4404 hazard$spline.knotspacing <- match.arg(hazard$spline.knotspacing, 4405 c("quantile", "equal", "mixed")) 4406 hazard$spline.penalty <- match.arg(hazard$spline.penalty, 4407 c("2diff", "2deriv", "log2deriv", "none")) 4408 hazard$param.dist <- match.arg(hazard$param.dist, 4409 c("none", "exponential", "weibull", "lognormal")) 4410 if(hazard$haspar & hazard$param.dist == "none") { 4411 warning(paste("no distribution specified for hazard parametric component", 4412 "-- setting to weibull")) 4413 hazard$param.dist <- "weibull" 4414 } 4415 # priors and hyperparameters for the number of knots 4416 hazard$spline.nknots.prior <- match.arg(hazard$spline.nknots.prior, 4417 c("poisson", "geometric", "poissonmix", "negbin", "power")) 4418 if(hazard$spline.nknots.prior == "poisson" & is.null(hazard$spline.nknots.hyper)) 4419 hazard$spline.nknots.hyper <- 10 4420 if(hazard$spline.nknots.prior == "geometric" & is.null(hazard$spline.nknots.hyper)) 4421 hazard$spline.nknots.hyper <- 0.1 4422 if(hazard$spline.nknots.prior == "poissonmix" & is.null(hazard$spline.nknots.hyper)) 4423 hazard$spline.nknots.hyper <- c(10, 30) 4424 if(hazard$spline.nknots.prior == "negbin" & is.null(hazard$spline.nknots.hyper)) 4425 hazard$spline.nknots.hyper <- c(2, .1) 4426 if(hazard$spline.nknots.prior == "power" & is.null(hazard$spline.nknots.hyper)) 4427 hazard$spline.nknots.hyper <- -1 / 2 4428 4429 # default weight settings 4430 if(!hazard$haspar) hazard$weight <- 1 4431 if(!hazard$hasspline) hazard$weight <- 0 4432 if(!frailty$haspar) frailty$weight <- 1 4433 if(!frailty$hasspline) frailty$weight <- 0 4434 4435 hazard$spline.norm <- FALSE 4436 hazard$name <- "hazard" 4437 hazard$x <- agdata$time 4438 4439 # Parse input (regression) 4440 reg.in <- regression 4441 reg.default <- list( 4442 priorvar = 0.1, 4443 hyper = c(0.01, 0.01), 4444 tun = 1, 4445 accept = 0 4446 ) 4447 regression <- reg.default 4448 regnames <- names(regression) 4449 if(!is.null(reg.in)) for(n in names(reg.in)) 4450 eval(parse(text = paste("regression$", match.arg(n, regnames), 4451 " <- reg.in$", n, sep = ""))) 4452 4453 # Default settings for the initial number of knots 4454 if(is.null(hazard$spline.nknots)) hazard$spline.nknots <- if(hazard$spline.adaptive) 4455 min(nknotsPriorMean(hazard), hazard$spline.maxoccknots) else 4456 max(min(round(sum(Ji) / 4), 35), 1) 4457 if(is.null(frailty$spline.nknots)) frailty$spline.nknots <- if(frailty$spline.adaptive) 4458 min(nknotsPriorMean(frailty), frailty$spline.maxoccknots) else 4459 max(1, min(round(m / 4), 35), 1) 4460 4461 4462 if(verbose >= 2) cat("\tFitting Cox survival models...\n") 4463 4464 # Cox fit with gamma frailties for initial values of Ui and beta 4465 varnames <- colnames(agdata)[ - (1:4)] 4466 qvarnames <- paste("`", varnames, "`", sep = "") 4467 if(m > 1){ 4468 coxfit <- coxph(as.formula(paste("Surv(time, delta)~", 4469 paste(qvarnames, collapse = " + "), " + frailty(i)")), data = agdata) 4470 Ui <- exp(coxfit$frail) 4471 if(var(Ui) < 1e-5) { 4472 Ui <- 2 * runif(length(Ui)) 4473 Ui <- 1 + (Ui - mean(Ui)) 4474 Ui[Ui < 0] <- 1 4475 } 4476 }else{ 4477 coxfit <- coxph(as.formula(paste("Surv(time, delta)~", 4478 paste(qvarnames, collapse = " + "))), data = agdata) 4479 Ui <- 1 4480 } 4481 # set initial frailies and regression structure 4482 frailty$x <- Ui 4483 beta <- coxfit$coef 4484 regression$m <- m 4485 regression$Ji <- Ji 4486 regression$covariates <- as.matrix(agdata[, -(1:4)], sum(Ji), length(beta)) 4487 regression$time <- agdata$time 4488 regression$status <- agdata$delta 4489 regression$cluster <- as.integer(agdata$i) 4490 regression <- updateregression(regression, beta) 4491 4492 rm(coxfit) 4493 4494 # Parametric fits 4495 if(verbose >= 2 & (frailty$haspar | hazard$haspar)) 4496 cat("\tFitting parametric components...\n") 4497 4498 hazard <- fitparametric(hazard, agdata) 4499 frailty <- fitparametric(frailty, Ui) 4500 4501 # Spline knots 4502 if(verbose >= 2 & (frailty$hasspline | hazard$hasspline)) 4503 cat("\tComputing spline knots...\n") 4504 4505 hazbounds <- range(agdata$time) + c(-.1, .1) * diff(range(agdata$time)) 4506 hazbounds[hazbounds < 0] <- 0 4507 hazard <- makeknots(hazard, agdata$time[agdata$delta == 1], bounds = hazbounds) 4508 frailty <- makeknots(frailty, Ui, bounds = c(0, max(max(Ui), min(2 * max(Ui), 4509 frailty$spline.maxknot)))) 4510 4511 # Evaluate the splines and integrals 4512 if(verbose >= 2 & (frailty$hasspline | hazard$hasspline)) 4513 cat("\tConstructing spline basis functions...\n") 4514 4515 hazard <- makesplinebasis(hazard, usec = usec) 4516 frailty <- makesplinebasis(frailty, usec = usec) 4517 4518 # Penalty matrices 4519 if(verbose >= 2 & (frailty$hasspline | hazard$hasspline)) 4520 cat("\tInitializing penalty matrices...\n") 4521 4522 hazard <- makepenalty(hazard, usec = usec) 4523 frailty <- makepenalty(frailty, usec = usec) 4524 4525 4526 if(verbose >= 2 & (frailty$hasspline | hazard$hasspline)) 4527 cat("\tObtaining initial values for spline parameters...\n") 4528 4529 {{{# Initial values for the theta vectors 4530 4531 if(hazard$haspar & hazard$hasspline){ 4532 oldhazweight <- hazard$weight 4533 hazard$weight <- 1; 4534 hazard <- weightcurve(hazard) 4535 } 4536 if(frailty$haspar & frailty$hasspline){ 4537 oldfrailweight <- frailty$weight 4538 frailty$weight <- 1; 4539 frailty <- weightcurve(frailty) 4540 } 4541 4542 # Initial values for hazard parameters 4543 if(hazard$hasspline){ 4544 theta.haz <- rep(0, hazard$spline.nknots + hazard$spline.ord) 4545 if(usec){ # use fast C code to compute likelihoods 4546 par <- as.double(theta.haz); status <- as.double(regression$status); 4547 lp <- as.double(regression$lp); frailrep <- as.double(rep(frailty$x, Ji)); 4548 hazParY <- as.double(if(hazard$haspar) hazard$param.y else rep(0, length(lp))); 4549 hazParYcum <- as.double(if(hazard$haspar) hazard$param.ycum else 4550 rep(0, length(lp))); 4551 weight <- hazard$weight; B <- as.double(hazard$spline.basis); 4552 C <- as.double(hazard$spline.basiscum); 4553 P <- as.double(hazard$spline.penaltyfactor * hazard$spline.penaltymatrix); 4554 penaltyType <- as.integer(pmatch(hazard$spline.penalty, 4555 c("none", "2diff", "2deriv", "log2deriv")) - 1); 4556 splinemin <- as.double(hazard$spline.min) 4557 sigma2 <- hazard$spline.priorvar 4558 sigma2 <- 100; sigma2target <- hazard$spline.priorvar 4559 #compute initial values by slowly decreasing the prior variance for stability 4560 while(sigma2 > sigma2target){ 4561 sigma2 <- as.double(sigma2 / 10) 4562 opt.theta.haz <- optim(par, fn = cmklik.spline.haz, gr = cmkgr.spline.haz, 4563 status = status, lp = lp, frailrep = frailrep, hazParY = hazParY, 4564 hazParYcum = hazParYcum, weight = weight, B = B, C = C, P = P, 4565 penaltyType = penaltyType, sigma2 = sigma2, min = splinemin, 4566 method = "BFGS", control = list(fnscale=-1)) 4567 par <- as.double(opt.theta.haz$par) 4568 } 4569 opt.theta.haz <- optim(par, fn = cmklik.spline.haz, gr = cmkgr.spline.haz, 4570 status = status, lp = lp, frailrep = frailrep, hazParY = hazParY, 4571 hazParYcum = hazParYcum, weight = weight, B = B, C = C, P = P, 4572 penaltyType = penaltyType, sigma2 = sigma2, min = splinemin, 4573 method = "BFGS", control = list(fnscale=-1, maxit = 1), hessian = TRUE) 4574 rm(par, status, lp, hazParY, hazParYcum, weight, B, C, P, penaltyType, 4575 splinemin, sigma2) 4576 gcout <- gc() 4577 }else{ # usec=FALSE 4578 hazard <- updatespline(hazard, theta.haz) 4579 gcout <- gc() 4580 sigma2 <- 100; sigma2target <- hazard$spline.priorvar 4581 #compute initial values by slowly decreasing the prior variance 4582 while(sigma2 > sigma2target){ 4583 sigma2 <- sigma2 / 10; 4584 hazard$spline.priorvar <- sigma2 4585 opt.theta.haz <- optim(hazard$spline.par, 4586 fn = mklik.spline.haz, 4587 gr = mkgr.spline.haz, 4588 method = "BFGS", 4589 control = list(fnscale=-1), 4590 hazard = hazard, 4591 frailty = frailty, 4592 regression = regression, 4593 hessian = FALSE) 4594 hazard <- updatespline(hazard, opt.theta.haz$par) 4595 } 4596 opt.theta.haz <- optim(hazard$spline.par, 4597 fn = mklik.spline.haz, 4598 gr = mkgr.spline.haz, 4599 method = "BFGS", 4600 control = list(fnscale=-1), 4601 hazard = hazard, 4602 frailty = frailty, 4603 regression = regression, 4604 hessian = TRUE) 4605 } 4606 gcout <- gc() 4607 hazard <- updatespline(hazard, opt.theta.haz$par) 4608 4609 } 4610 4611 # Initial values for frailty parameters 4612 if(frailty$hasspline){ 4613 frailty$spline.fixedind <- 1 4614 # parameter vector not including the fixed index 4615 theta.frail <- rep(0, frailty$spline.nknots + frailty$spline.ord - 1) 4616 { 4617 opt.theta.frail <- optim(theta.frail, 4618 fn = mklik.spline.frail.init, 4619 method = "BFGS", 4620 control = list(fnscale=-1), 4621 hazard = hazard, 4622 frailty = frailty, 4623 regression = regression, 4624 hessian = TRUE) 4625 } 4626 gcout <- gc() 4627 # add the fixed index back in 4628 theta.frail <- repairfrailtypar(opt.theta.frail$par, frailty$spline.fixedind) 4629 badtheta <- TRUE; j <- 0; 4630 # set the mean of the estimated frailty density to be exactly 1 4631 while(badtheta){ 4632 j <- j + 1 4633 thetaj <- suppressWarnings(log(-sum(frailty$spline.basisexp[ - j] * exp(theta.frail[ - j])) / frailty$spline.basisexp[j])) 4634 if(!is.nan(thetaj)) badtheta <- FALSE 4635 } 4636 theta.frail[j] <- thetaj; 4637 frailty <- updatespline(frailty, theta.frail) 4638 } 4639 4640 gcout <- gc() 4641 # reweight the inital curves 4642 if(hazard$hasspline & hazard$haspar){ 4643 hazard$weight <- oldhazweight 4644 hazard <- weightcurve(hazard); 4645 } 4646 if(frailty$hasspline & frailty$haspar){ 4647 frailty$weight <- oldfrailweight 4648 frailty <- weightcurve(frailty) 4649 } 4650 4651 gcout <- gc() 4652 # Evaluate variances and hessians for candidate generation 4653 frailty$tun <- diff(range(Ui))^2 / 6 4654 hess.coef <- mkhess.coef(regression$coefficients, hazard, frailty, regression) 4655 Sigma.coef <- solve(-hess.coef) 4656 regression$candcov <- Sigma.coef 4657 regression$cholcandcov <- chol(Sigma.coef, pivot = TRUE) 4658 if(hazard$hasspline){ 4659 hess.haz <- opt.theta.haz$hess 4660 Sigma.haz <- inverthessian(hess.haz) 4661 hazard$spline.candcov <- Sigma.haz 4662 hazard$spline.candsd <- rep(1, hazard$spline.maxoccknots + hazard$spline.ord) 4663 hazard$spline.cholcandcov <- chol(Sigma.haz, pivot = TRUE) 4664 rm(hess.haz, Sigma.haz) 4665 } 4666 if(frailty$hasspline){ 4667 hess.frail <- opt.theta.frail$hess 4668 Sigma.frail <- inverthessian(hess.frail) 4669 frailty$spline.candcov <- Sigma.frail 4670 frailty$spline.candsd <- rep(1, frailty$spline.maxoccknots + frailty$spline.ord) 4671 frailty$spline.cholcandcov <- chol(Sigma.frail, pivot = TRUE) 4672 rm(hess.frail, Sigma.frail) 4673 } 4674 # construct a numeric hessian for the parametric parameters 4675 if(hazard$haspar){ 4676 # (I don't remember why it doesn't just use numHess.par like for the frailty, but 4677 # there must be a reason for this inelegant setup) 4678 temphaz <- list(haspar = hazard$haspar, hasspline = hazard$hasspline, 4679 weight = hazard$weight, spline.y = hazard$spline.y, spline.ycum = hazard$spline.ycum, 4680 name = hazard$name, param.dist = hazard$param.dist, x = hazard$x, y = hazard$y, 4681 ycum = hazard$ycum, param.y = hazard$param.y, param.ycum = hazard$param.ycum, 4682 param.par = hazard$param.par, param.priorvar = hazard$param.priorvar) 4683 eps <- 1e-5 4684 par <- temphaz$param.par 4685 hess <- matrix(0, length(par), length(par)) 4686 for(i in 1:length(par)){ 4687 # this constructs numerical gradients by finite differences 4688 for(j in 1:length(par)){ 4689 par1 <- par;par2 <- par;par3 <- par; 4690 if(i == j) {par1[i] <- par[i] + eps;par2 <- par1;par3[i] <- par[i] + 2 * eps} 4691 else {par1[i] <- par[i] + eps; par2[j] <- par[j] + eps; par3 <- par1; 4692 par3[j] <- par[j] + eps} 4693 g1 <- (mklik.param.haz(par1, temphaz, frailty, regression) - 4694 mklik.param.haz(par, temphaz, frailty, regression)) / eps 4695 g2 <- (mklik.param.haz(par3, temphaz, frailty, regression) - 4696 mklik.param.haz(par2, temphaz, frailty, regression)) / eps 4697 hess[i, j] <- (g2 - g1) / eps 4698 } 4699 } 4700 Sigma.par.haz <- inverthessian(hess) 4701 hazard$param.candcov <- Sigma.par.haz 4702 hazard$param.cholcandcov <- chol(Sigma.par.haz, pivot = TRUE) 4703 rm(Sigma.par.haz, hess, temphaz) 4704 } 4705 if(frailty$haspar){ 4706 Sigma.par.frail <- inverthessian(numHess.par(frailty$param.par, mklik.param.frail, 4707 hazard = hazard, frailty = frailty, regression = regression)) 4708 frailty$param.candcov <- Sigma.par.frail 4709 frailty$param.cholcandcov <- chol(Sigma.par.frail, pivot = TRUE) 4710 rm(Sigma.par.frail) 4711 } 4712 4713 4714 }}} 4715 4716 if(frailty$hasspline) frailty$spline.fvar <- frailtysplinefvar(frailty) 4717 #browser() 4718 gcout <- gc() 4719 # Store initial values in parameter history 4720 history <- inithistory(hazard, frailty, regression, control) 4721 avg.tunhist <- NULL 4722 avg.accepthist <- NULL 4723 4724 4725 # an internal function to prepare the curve structures for calling the C code 4726 prepforCall <- function(curve){ 4727 if(curve$type == "parametric") return(curve) 4728 if(curve$spline.nknots < curve$spline.maxoccknots){ 4729 addcols <- curve$spline.maxoccknots - curve$spline.nknots 4730 curve$spline.knots <- c(curve$spline.knots, rep(-Inf, addcols)) 4731 curve$spline.par <- c(curve$spline.par, rep(-Inf, addcols)) 4732 curve$spline.candsd <- c(curve$spline.candsd, rep(-Inf, addcols)) 4733 curve$spline.basis <- cbind(curve$spline.basis, matrix(-Inf, 4734 dim(curve$spline.basis)[1], addcols)) 4735 curve$spline.basisint <- c(curve$spline.basisint, rep(-Inf, addcols)) 4736 if(curve$name == "hazard") curve$spline.basiscum <- cbind(curve$spline.basiscum, 4737 matrix(-Inf, dim(curve$spline.basiscum)[1], addcols)) 4738 if(curve$name == "frailty") curve$spline.basisexp <- c(curve$spline.basisexp, 4739 rep(-Inf, addcols)) 4740 } 4741 curve$spline.candocc <- attr(curve$spline.candknots, "occupied") 4742 return(curve) 4743 } 4744 if(usec){ 4745 hazard <- prepforCall(hazard) 4746 frailty <- prepforCall(frailty) 4747 } 4748 4749 # this does nothing, it just creates a useful marker 4750 main <- function() {} 4751 4752 if(verbose >= 1) cat("Starting MCMC...\n") 4753 4754 iter <- 1 # Counts the recorded iterations 4755 iter.el <- 0 # Counts elapsed iterations in each thinning cycle 4756 4757 if(verbose >= 3) cat(iter, " ") 4758 4759 4760 while(iter < control$maxiter) 4761 { 4762 if(usec){ 4763 # C version of the main loop 4764 gcout <- gc() 4765 nexttunint <- iter - iter%%control$tun.int + control$tun.int 4766 enditer <- min(nexttunint, control$maxiter) 4767 out <- .Call("SplineSurvMainLoop", hazard, frailty, regression, history, iter, 4768 enditer, control$thin, verbose) 4769 iter <- enditer 4770 }else{ 4771 4772 # R version of the main loop 4773 4774 iter.el <- iter.el + 1 4775 4776 if(verbose >= 4) cat(iter.el) 4777 4778 # MH update of frailties 4779 frailty <- mh.frail(hazard, frailty, regression) 4780 4781 # MH update of regression parameters 4782 regression <- mh.coef(hazard, frailty, regression) 4783 4784 # MH update of baseline parameters 4785 hazard <- mh.hazard.spline(hazard, frailty, regression) 4786 4787 # MH update of frailty density parameters 4788 frailty <- mh.frailty.spline(hazard, frailty, regression) 4789 4790 # MH update of parametric baseline parameters 4791 hazard <- mh.hazard.param(hazard, frailty, regression) 4792 4793 # MH update of parametric frailty parameters 4794 frailty <- mh.frailty.param(hazard, frailty, regression) 4795 4796 # MH update of weights 4797 hazard <- mh.weight("hazard", hazard, frailty, regression) 4798 frailty <- mh.weight("frailty", hazard, frailty, regression) 4799 4800 # Update of the sigmas / taus 4801 hazard <- updatepostvar.curve(hazard) 4802 frailty <- updatepostvar.curve(frailty) 4803 regression <- updatepostvar.coef(regression) 4804 4805 # Birth - death - move 4806 if(hazard$spline.adaptive) 4807 hazard <- mh.bdm("hazard", hazard, frailty, regression) 4808 if(frailty$spline.adaptive) 4809 frailty <- mh.bdm("frailty", hazard, frailty, regression) 4810 4811 # record history 4812 if(iter.el == control$thin){ 4813 iter <- iter + 1 4814 history <- updatehistory(history, iter, hazard, frailty, regression) 4815 iter.el <- 0 4816 if(verbose >= 3) cat(" ", iter, " ", sep = "") 4817 } 4818 } 4819 4820 {{{ # Periodic calibration check 4821 if(iter%%control$tun.int == 0 & iter < control$maxiter){ 4822 4823 if(verbose == 1 | verbose == 2) cat(iter, " ") 4824 4825 calinds <- (iter - control$tun.int + 1):iter 4826 4827 # Calibration of the tuning parameters for acceptance rate 4828 if(control$tun.auto & iter <= control$burnin){ 4829 if(verbose >= 3) cat("\n Calibration ...\n") 4830 4831 4832 tunnames <- c("regression$tun", 4833 "hazard$spline.tun", 4834 "frailty$spline.tun", 4835 "hazard$param.tun", 4836 "frailty$param.tun", 4837 "hazard$weight.tun", 4838 "frailty$weight.tun", 4839 "frailty$tun") 4840 alltun <- rep(0, length(tunnames)) 4841 for(i in 1:length(alltun)) eval(parse(text = paste("alltun[", i, "] <- ", 4842 tunnames[i]))) 4843 avg.tunhist <- rbind(avg.tunhist, alltun) 4844 avg.accepthist <- rbind(avg.accepthist, 4845 apply(history$accept[calinds, ], 2, mean)) 4846 for(g in 1:length(alltun)){ 4847 if(all(avg.accepthist[, g]>.25)) alltun[g] <- alltun[g] * 2 4848 if(all(avg.accepthist[, g]<.25)) alltun[g] <- alltun[g] / 2 4849 if(any(avg.accepthist[, g]>.25) & any(avg.accepthist[, g]<.25)){ 4850 # build linear regression model for tuning parameters 4851 fit <- lm(y~x, data.frame(x = avg.tunhist[, g], 4852 y = avg.accepthist[, g])) 4853 out <- try(max(1e-3, nlm(accrate.predict.lm, alltun[g], m = fit)$est)) 4854 if(inherits(out, "try-error")) 4855 alltun[g] <- avg.tunhist[which.min((avg.accepthist-.25)^2)] 4856 else 4857 alltun[g] <- out 4858 } 4859 } 4860 for(i in 1:length(alltun)) eval(parse(text = paste(tunnames[i], 4861 " <- alltun[", i, "]"))) 4862 4863 if(verbose >= 4){ 4864 outmat <- cbind(avg.accepthist, alltun); 4865 rownames(outmat) <- tunnames; 4866 colnames(outmat) <- c("acceptance", "tuning") 4867 print(outmat) 4868 } 4869 } 4870 # Print full iteration info 4871 if(verbose >= 5){ 4872 cat("Frailties: ", submean(history$frailty, calinds), "\n") 4873 cat("Coefficients: ", submean(history$coefficients, calinds), "\n") 4874 cat("Hazard.spline: ", submean(history$hazard.spline.par, calinds), "\n") 4875 cat("Frailty.spline: ", submean(history$frailty.spline.par, calinds), "\n") 4876 cat("Hazard.param: ", submean(history$hazard.param.par, calinds), "\n") 4877 cat("Frailty.param: ", submean(history$frailty.param.par, calinds), "\n") 4878 cat("Prior variances: ", submean(history$priorvar, calinds), "\n") 4879 } 4880 4881 #browser() 4882 4883 } 4884 }}} 4885 4886 } # end main loop 4887 4888 gcout <- gc() 4889 if(verbose > 0) cat("Done!\n") 4890 hazard <- makeoutputcurve(hazard) 4891 frailty <- makeoutputcurve(frailty) 4892 4893 gcout <- gc() 4894 {{{ # Construct output 4895 if(control$burnin < iter){ 4896 if(hasspline(frailty)){ 4897 frail.par.rs <- rowSums(exp(history$frailty.spline.par)) 4898 history$frailty.spline.par <- exp(history$frailty.spline.par) / frail.par.rs 4899 } 4900 sub <- (1:(iter)) > (control$burnin) 4901 posterior.mean <- list(coefficients = submean(history$coefficients, sub), 4902 frailty = submean(history$frailty, sub), 4903 hazard.spline.par = submean(history$hazard.spline.par, sub), 4904 hazard.param.par = submean(history$hazard.param.par, sub), 4905 hazard.weight = submean(history$hazard.weight, sub), 4906 frailty.spline.par = submean(history$frailty.spline.par, sub), 4907 frailty.param.par = submean(history$frailty.param.par, sub), 4908 frailty.weight = submean(history$frailty.weight, sub) 4909 ) 4910 # save mean number of knots in spline.nknots for output 4911 if(hasspline(hazard)) 4912 hazard$spline.nknots <- mean(apply(history$hazard.spline.knots[sub, ], 1, 4913 function(x) sum(x>-Inf))) - 2 * hazard$spline.ord 4914 if(hasspline(frailty)) { 4915 frailty$spline.nknots <- mean(apply(history$frailty.spline.knots[sub, ], 1, 4916 function(x) sum(x>-Inf))) - 2 * frailty$spline.ord 4917 history$frailty.spline.par <- log(history$frailty.spline.par) 4918 posterior.mean$frailty.spline.par <- log(posterior.mean$frailty.spline.par) 4919 } 4920 colnames(history$coefficients) <- varnames 4921 names(posterior.mean$coefficients) <- varnames 4922 }else{ 4923 postmean = NULL 4924 } 4925 if(coda){ 4926 # convert history to mcmc objects 4927 library(coda) 4928 for(i in 1:length(history)) history[[i]] <- as.mcmc(history[[i]]) 4929 } 4930 # clear history rownames 4931 rownames(history$frailty) <- rownames(history$coefficients) <- 4932 rownames(history$splinepar.haz) <- rownames(history$splinepar.frail) <- NULL 4933 control$iter <- iter 4934 }}} 4935 4936 gcout <- gc() 4937 # output object 4938 out <- list(call = call, history = history, posterior.mean = posterior.mean, 4939 hazard = hazard, frailty = frailty, control = control, data = agdata) 4940 class(out) <- "splinesurv" 4941 return(out) 4942 }