Dear BUGSusers,
sorry for a late reporting of the solution. I was out of the office for some
time.
I have received several suggestions. Thanks for all of them. Probably nicest
solutions can be found below.
========================================================================
ORIGINAL PROBLEM
========================================================================
Dear BUGSusers,
does anybody have an experience with doubly interval censored data in Bugs?
My data consists of pairs of intervals giving the onset (time when the
individual starts to be at risk) and death times. I.e. it is known that a
onset time A lies between u[1] and u[2] and a death time C lies between u[3]
and u[4], u[1] < u[2] and u[3] < u[4]. The quantity of a main interest is then
a difference D = C - A. There is no problem to assign a parametric
distribution to both time quantities A and C, and to use Bugs function I() to
indicate interval censored data.
e.g.
model{
for (i in 1:N) { ## loop over persons
a[i] ~ dweib(gammaa, lambdaa)I(leftem[i], aupper[i])
blower[i] <- max(a[i], leftcar[i])
c[i] ~ dweib(gammab, lambdab)I(blower[i], rightcar[i])
}
where leftem = u[1]
rightem = u[2]
leftcar = u[3]
rightcar = u[4]
Program of this type works if for all individuals u[1] < u[2] < u[3] < u[4].
The problem occurs if some individual experiences both onset and death during
one interval, i.e. if u[1]=u[3] and u[2]=u[4]. In that case Bugs always
crashes.
==========================================================================
Solution by Andrew Millard
--------------------------
I should have seen this first time - your censoring was not symmetrical.
This updates successfully:
model{
for (i in 1:N) {
a[i] ~ dweib(gammaa, lambdaa)I(leftem[i], aupper[i])
aupper[i] <- min(rightem[i], c[i]) ## ADD THIS LINE
blower[i] <- max(a[i], leftcar[i])
c[i] ~ dweib(gammab, lambdab)I(blower[i], rightcar[i])
}
You had forced new cvalues of c to be more than a, but new values of a did
not have to be less than c.
The Jama Valley example on my website is rather more complicated, but has
the same multiple censoring of unknowns by one another.
Andrew
============================================================================
Same solution by Andrew Thomas
------------------------------
Hello Arnost,
I can not exactly reproduce your problem, but then I am using a
later version of WinBUGS. There were also problems with my version of
WinBUGS. Like you say there are no problems if rightem < leftcar.
In this case the max function does not do anything. If I let the two
intervals overlap WinBUGS was not able to generate samples, one of the
intervals shrunk to zero. When I look at your code I see the treatment of
a[i] and c[i] are not symmetric. So I introduced
aupper[i] <- min(c[i], rightem[i])
a[i] ~ dweib(gammaa, lambdaa)I(leftem[i], aupper[i])
things then seemed to work.
Interval censored data somewhat like you have crops up in archoeology with
radio carbon data and stratification (see Andrew Millards web page)
Regards
Andrew
===============================================================================
=
Additional comments
-------------------
Andrews' solutions work almost perfectly. The problem only occurs if both
sampled A and C are too close each other. However then "epsilon trick" helps.
e.g.
model{
for (i in 1:N) {
a[i] ~ dweib(gammaa, lambdaa)I(leftem[i], aupper[i])
aupper[i] <- min(rightem[i], c[i]) - 0.05
blower[i] <- max(a[i], leftcar[i]) + 0.05
c[i] ~ dweib(gammab, lambdab)I(blower[i], rightcar[i])
}
David Spiegelhalter also suggested to use
a[i] ~ dweib(gammaA, lambdaA)I(leftem[i],b[i])
b[i] ~ dweib(gammaB, lambdaB)I(a[i],rightcar[i])
for the cases where u[1]=u[3] and u[2]=u[4]. This would probably also work
but it would be more complicated to indicate these cases where u[1]=u[3]
and u[2]=u[4]. So that I think that Andrews' solution is better.
Arnost
-------------------------------------------
Arnost Komarek
Biostatistisch Centrum K.U. Leuven
U.Z. St. Rafael
Kapucijnenvoer 35
B-3000 Leuven
Belgie
Tel : +32 - 16 - 33 68 86
Fax: +32 - 16 - 33 69 00
Web: http://www.med.kuleuven.ac.be/biostat/
-------------------------------------------------------------------
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
|