Buggsies,
I'm trying to fit a logistic regression model using WinBUGS. I'm using the
Williams rat-pup data set, which I'd like to code in a binary format as a
practice run for a problem we're working on.
I've got 303 rat pups from 32 litters. Dams received a treatment during
pregnancy, and I'm interested in quantifying the effect of treatment on pup
survival (survived = yes, survived = no). Pups are clustered within litter,
so I'd like to consider litter as a random effect.
The code I've come up with so far is shown below. When I compile this model,
BUGS tells me it is looking for a key word structure for 'live'. This makes
sense, I've coded live[i,j] so the data needs to contain an array.
This has me stumped - how do make an array when there are different numbers
of rats in each litter? Am I completely off the mark? Can someone point me
in right direction?
Regards (and thanks!),
Mark Stevenson
Code ...
=========================
model
{for (i in 1:litter) {
for(j in offset[i]:(offset[i+1]-1))
r[i,j] ~ dbin(live[i,j], denom[i,j])
logit(live[i,j]) <- beta0 + (beta1*treat[i,j]) + h[i]
}
h[i] ~ dnorm(0.0, tau)
}
beta0 ~ dnorm(0.0, 1.0E-6)
beta1 ~ dnorm(0.0, 1.0E-6)
beta2 ~ dnorm(0.0, 1.0E-6)
tau ~ dgamma(0.001, 0.001)
sigma <- 1 / sqrt(tau)
}
=========================
DATA
list(litter = 32,
offset =
c(1,13,24,34,43,54,64,74,83,92,97,106,113,123,129,139,146,159,171,180,189,
197,205,218,230,240,250,259,272,277,284,294,303),
# Tells you which pup belongs to which litter. Note 33 elements in this
vector.
treat =
c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1 ... to 303),
live =
c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,1,0,
0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,
1,0,0,0,0,0,0,0,0,1,0,0,0,0,0 ... to 303),
denom =
c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1 ... to 303))
=========================
INITS
list(beta0 = 0, beta1 = 0, beta2 = 0, tau = 1)
=========================
*******************************************
Mark Stevenson
EpiCentre, IVABS
Massey University
Private Bag 11-222 Palmerston North
NEW ZEALAND
Ph: +64 6 350 6149
Fx: +64 6 350 5716
[log in to unmask]
http://www.dairywin.co.nz/
*******************************************
|