[R] WinBUGS on survival, simple but confusing question

[tingtingzhan]

Hi All,

I'm using WinBUGS on a very simple survival model (log-normal with one
covariate "Treat"), but I cannot understand the way it handles censored
data.  I'm posting the R file which generates the data from pre-specified
parameters, as well as the .bug file.

The question is, if I use NA to denote the censored data (as suggested by
the example Mice in WinBUGS Example Vol.I), 90% of the time I get the error
message:
"value of log normal surt.BUGS[ii] must be positive"
where "ii" is the index of the first NA in my simulated survival data.  
However, there are some very rare occasions that WinBUGS performs as we
expect, and I included on such dataset in the following code.

Any help would be greatly appreciated, if someone could try out these code
and let me know if I missed anything.

Best wishes,
Tingting Zhan



###################################
## R file.   Save to "~/BUGS.survival.simple.R"
## source("~/BUGS.survival.simple.R")

rm(list=ls(all=TRUE));
library(survival);
library(R2WinBUGS);
library(BRugs);


## True parameters
alpha = c(2.31, 0.14);
scale = 0.25;

## sample size
n.Treat = n.Control = 20;
N = n.Treat + n.Control;

## covariate
Treat = c(rep(0, n.Control), rep(1, n.Treat))

## true survival time (which is censored at day 14)
surt.cens = rep(14, N);

surt = # round(c(rlnorm(n.Control, alpha[1], scale), rlnorm(n.Treat,
alpha[1]+alpha[2], scale)));
         c(9, 9, 7, 9, 8, 10, 10, 10, 7, 8, 9, 8, 9, 7, 8, 9, 9, 11, 8, 8,
13, 9, 18, 11, 8, 11, 22, 11, 14, 8, 12, 16, 13, 10, 15, 15, 8, 16, 11, 10);  
## occasionally works


## survival time and censoring time to be passed into WinBUGS
surt.BUGS = surt;
surt.BUGS[surt >= surt.cens] = NA;  
if(any(is.na(surt.BUGS))) cat("surt.BUGS taking NA (censored) value at",
which(is.na(surt.BUGS)), ".\n");

surt.cens.BUGS = surt.cens;
surt.cens.BUGS[surt < surt.cens] = 0;

## non-info prior
sd.gamma1 = sd.gamma2 = 1e-2;
pos.lim = 1e-4;
norm.sd = 1e-3;

BUGS.data = list("N", "Treat", "surt.BUGS", "surt.cens.BUGS", "sd.gamma1",
"sd.gamma2", "norm.sd", "pos.lim");

BUGS.parameters = c("alpha", "surt.sigma", "surt.BUGS");

BUGS.sim = bugs(BUGS.data, inits = NULL, BUGS.parameters, model.file =
"~\\BUGS.survival.bug", n.chains = 7, n.iter = 700, debug = T);

## end of R file
###################################


###################################
## bug file.   Save to "~/BUGS.survival.bug"

model{
  for(i in 1:N)
  {
    surt.BUGS[i] ~ dlnorm(surt.mean[i], surt.tau)I(surt.cens.BUGS[i], );
    surt.mean[i] <- alpha[1] + alpha[2]*Treat[i];
  }

  surt.tau ~ dgamma(sd.gamma1, sd.gamma2)I(pos.lim,);
  surt.sigma <- sqrt(1/surt.tau);

  alpha[1] ~ dnorm(0, norm.sd);
  alpha[2] ~ dnorm(0, norm.sd);
}

## end of bug file
###################################



--
View this message in context: http://r.789695.n4.nabble.com/WinBUGS-on-survival-simple-but-confusing-question-tp3580685p3580685.html
Sent from the R help mailing list archive at Nabble.com.