Dear all,
My dataset is a series of readings of blood pressure taken from
a cohort of ~2000 patients. I also have data on how long they survived
before they had a stroke.
I am trying to model the longitudinal blood pressure and stroke
survival processes jointly - giving each patient their own
intercept for blood pressure (a01) and within-individual error
for readings of blood pressure (precision=tau or standard
deviation=sigma). My model specifies that both the intercept,
a01, and standard deviation, sigma, predict stroke risk. That is
individuals j's survival time ~ Weibull(r, hr[j]) where
-log(hr[j]) = b0 + b1*a01[j] + b2*sigma[j].
The trouble is that sigma isn't centred on zero and consequently
the correlation between b0 and b2 is very high (>90%) and the
autocorrelation for both parameters is substantial (>90%) even
at large lags (eg. 50).
If I try to calculate the mean of sigma[]
(ie. sigma.mean <- mean(sigma[])), there is no problem,
but if I then try to centre sigma (ie. centred.sigma[j] <- sigma[j] -
sigma.mean), and use that in the survival model,
then the model suddenly takes at least half an hour to compile
and I haven't yet waited long enough for it to complete a
single update. Without the attempted centering, it takes maybe 15
seconds to compile, and updates fine.
Does anybody have any suggestions about why this could be happening?
Thanks for your time,
Robert
model {
for (j in 1:J) { # For each individual
# Specify the distribution of the random intercept
# and within-individual error for blood pressure
# model
a01[j] ~ dnorm(0, tau2)
tau[j] ~ dgamma(g1, g2)
# The survival model uses the standard deviation of
# each individual's blood pressure, not the precision.
sigma[j] <- 1/sqrt(tau[j])
# Try to centre the standard deviation to reduce the
# correlation between it's survival coefficient (b2)
# and the intercept.
# Adding this line causes the problem described
# above
centred.sigma[j] <- sigma[j] - sigma.mean
# Survival model
lp[j] <- b0 + b1*a01[j] + b2*centred.sigma[j]
hr[j] <- exp(-lp[j])
t[j] ~ dweib(r,hr[j])I(t.cen[j],)
}
for( i in 1:N){ # For each reading of blood pressure
# Set up an random effects model for the progression
# of blood pressure with time (called "start" here).
# The dataset is unbalanced, so I use the vector idno
# to match each reading to the patient from whom it comes.
lsbp[i] ~ dnorm( mu[i], tau[idno[i]] )
mu[i] <- a0 + a01[idno[i]] + a1*start[i]
}
# Calculate the mean of sigma for use in centring
sigma.mean <- mean(sigma[])
a0 ~ dnorm(4.5, 4)
a1 ~ dnorm(-.1, 25)
b0 ~ dnorm(3.5, .25)
b1 ~ dnorm(0, .0001)
b2 ~ dnorm(0, .0001)
r ~ dgamma(1,.001)
g1 ~ dexp(.3)
g2 ~ dgamma(1,.1)
tau2 ~ dgamma(1, .1)
sigma2 <- 1/sqrt(tau2)
}
-------------------------------------------------------------------
This list is for discussion of modelling issues and the BUGS software.
For help with crashes and error messages, first mail [log in to unmask]
To mail the BUGS list, mail to [log in to unmask]
Before mailing, please check the archive at www.jiscmail.ac.uk/lists/bugs.html
Please do not mail attachments to the list.
To leave the BUGS list, send LEAVE BUGS to [log in to unmask]
If this fails, mail [log in to unmask], NOT the whole list
|