I received many emails but below is a summary of what I was looking for.
I was able to solve the probit model by myself, but I was looking for a
logit representation.
Dr. Burn sent an example that was similar to what I was looking for.
Dr. Burn logit example.
My probit example is consistent with the results with the Bayesian
estimates from the software 'AMOS'.
model{
for(i in 1:N){
y[i]~dnorm(mu[i],tau)I(thd[rich93cat[i]],thd[rich93cat[i]+1])
mu[i]<-beta0 + beta1*rfloodmn[i]}
beta0~dnorm(0,.0001)
beta1~dnorm(0,.0001)
tau~dgamma(.001,.001)
sigma <- 1/sqrt(tau)
}
DATA list(thd=c(-5000000, -1.00313, 0.93964, 5000000),
rfloodmn=c(48.5,57,47,53.5,50,49,51,47,46.5,50,54.5,43.5,47.5,28,27,55.5,
31,30,28.5,36,39,50.5,28,21.5,21,30,26,24,26.5,39,41,54.5,45,
41.5,51.5,43,35,32.5,25,19.5,41,57.5,25.5,22.5,40,31.5,25.5,4,
4,1,18.5,59.5,25.5,42.5,38,14,19.5,20.5,18.5,19.5,21,21,22.5,
59.5,23.5,25.5,24,34.5,27.5,33,46.5,42,47,48.5,51,46.5,58.5,53,44.5,45,42.5,
36.5,45,59.5,40,39,36.5,40.5,33.5,37.5,50,48,32,
39,53,49.5,40,40.5,38,60,46.5,12,25,21,37,57.5,39,30.5,12.5,14,53.5,46.5,28,
24.5,24.5,24.5,29.5,23.5,29.5,24,50.5,47,43,57.5,
49,37.5,39,51,45.5,41,27,25.5,45,57.5,52.5,38.5,40.5,39.5,37,
34,55,59,36.5,33.5,32.5,36,55.5,35,42.5,35.5,53,58,34.5,48.5,
28.5,21.5,20,20,49,23.5,34,41.5,35.5,40,29.5,51,28,24,20.5,23,
24.5,19,20.5,17,17,23,38.5,48.5,38.5,55.5,46,54.5,41.5,35.5,37,
50,46,38.5,44.5,40.5),N=190,
rich93cat=c(1,2,1,1,2,1,1,2,2,2,1,2,1,2,2,1,2,2,2,1,2,1,2,2,2,
2,2,2,2,2,2,1,2,2,1,2,2,2,2,2,1,1,3,3,1,2,3,2,2,2,2,1,2,2,2,2,2,
2,2,3,2,3,3,2,3,3,3,2,2,2,2,2,2,2,2,1,1,1,2,2,3,3,2,1,2,2,3,3,2,
2,2,2,3,3,1,1,2,2,3,1,1,2,2,2,2,1,2,2,2,2,1,2,2,2,2,2,3,2,2,2,
2,3,3,3,2,2,3,2,2,3,3,2,2,2,1,2,2,2,2,2,1,2,2,2,1,2,2,1,2,2,2,
2,2,2,3,3,2,2,2,2,3,2,2,2,3,3,3,3,2,2,2,2,2,3,2,2,2,2,2,2,2,2,2,3,3,2,2,2,2,
2))
INITIAL VALUES list(beta0=0,beta1=0,tau=1)
The original question I was given at work was to convert the AMOS model
into WinBugs.
I hope this helps.
I have attached the WinBugs programs for whoever can use them.
If you have any questions, feel free to contact me.
Supplementary Practical 2
Pneumoconiosis in Miners
Ordered categorical outcomes
Categorical outcomes often ordered. For example:
severity of symptoms: none, mild, severe
university degree: first, upper second, lower second, third
Sometimes analysed as categorical data ignoring the ordering. A much
better analysis takes account of the ordering.
A study of lung disease in miners (Ashford, 1959) classified miners
according to severity of symptoms of pneumoconiosis and years of working in
the mines (exposure).
The exposure variable was grouped and the values are the group mid-points.
Pneumoconiosis
Exposure (yrs) Normal Mild Severe
5.8 98 0 0
15.0 51 2 1
21.5 34 6 3
27.5 35 5 8
33.5 32 10 9
39.5 23 7 8
46.0 12 6 4
51.5 4 2 5
The aim of the analysis is to assess the importance of exposure on the
severity of the disease.
Theory
Basic idea: assume an underlying continuous random variable Y whose
distribution has been grouped into the observed categories.
Y is a latent variable - i.e. unobservable.
Let be the values (called the cut-points) of Y corresponding to the
boundaries between the k observed categories.
Sometimes convenient to denote and , so that we have .
Even though Y is unobservable, the cut-points can be estimated.
Let be the probability that individual i is in category j (= 1, …, k).
Then .
The cumulative probability that individual i is in category j or below is .
So .
The model proposed by McCullagh (1980) is
where is a linear predictor representing characteristics of the ith
individual. In this case we could have
( = exposure).
This is called the proportional odds model. The reason is that if we
calculate the ratio of the odds at exposure x1 and the odds at exposure x2,
the result is just
and this is independent of the category j.
Earlier work on this example suggested that it was preferable to use log
(exposure) in this case. So the model is
In the example there are 3 categories and therefore 2 cut-points ( and )
to estimate.
The main focus of interest is the parameter which is a measure of the
effect of exposure on the severity of the disease.
Priors
Non-informative priors can be normal with low precision.
(Not much chance of informative priors for the cut-points, since they are
values of an unobservable variable!)
Sensible to impose the constraint in the priors.
This can be done by using the truncation operator I(a,b) in WinBUGS …
… for example
z ~ dnorm(mu,tau)I(a,)
constrains z to be greater than a and
z ~ dnorm(mu,tau)I(,b)
constrains z to be less than b.
So you can use the trick
theta[1] ~ dnorm(0.0, 1.0E-6)I( ,theta[2])
theta[2] ~ dnorm(0.0, 1.0E-6)I(theta[1], )
for the priors of the cut-points.
Examine the code in file supp 2 ex 1.odc and make sure that you can follow
it. Write a script to run the model for 2 chains, check convergence.
What can be concluded?
model
{
for (i in 1:N) {
mu[i] <- beta*log(period[i])
}
for (j in 1:2) {
for (i in 1:N) {
logit(gamma[i,j]) <- theta[j] - mu[i]
}
}
for (i in 1:N) {
p[i,1] <- gamma[i,1]
p[i,2] <- gamma[i,2] - gamma[i,1]
p[i,3] <- 1-gamma[i,2]
y[i] ~ dcat(p[i,])
}
theta[1] ~ dnorm(0, 1.0E-03)I(,theta[2])
theta[2] ~ dnorm(0, 1.0E-03)I(theta[1],)
beta ~ dnorm(0,0.001)
}
# Data:
list(period=c
(5.8,5.8,5.8,5.8,5.8,5.8,5.8,5.8,5.8,5.8,5.8,5.8,5.8,5.8,5.8,5.8,5.8,5.8,5.8
,5.8,
5.8,5.8,5.8,5.8,5.8,5.8,5.8,5.8,5.8,5.8,5.8,5.8,5.8,5.8,5.8,5.8,5.8,5.8,5.8,
5.8,5.8,5.8,
5.8,5.8,5.8,5.8,5.8,5.8,5.8,5.8,5.8,5.8,5.8,5.8,5.8,5.8,5.8,5.8,5.8,5.8,5.8,
5.8,5.8,5.8,
5.8,5.8,5.8,5.8,5.8,5.8,5.8,5.8,5.8,5.8,5.8,5.8,5.8,5.8,5.8,5.8,5.8,5.8,5.8,
5.8,5.8,5.8,
5.8,5.8,5.8,5.8,5.8,5.8,5.8,5.8,5.8,5.8,5.8,5.8,15,15,15,15,15,15,15,15,15,1
5,15,15,15,
15,15,15,15,15,15,15,15,15,15,15,15,15,15,15,15,15,15,15,15,15,15,15,15,15,1
5,15,15,15,
15,15,15,15,15,15,15,15,15,15,15,15,21.5,21.5,21.5,21.5,21.5,21.5,21.5,21.5,
21.5,21.5,
21.5,21.5,21.5,21.5,21.5,21.5,21.5,21.5,21.5,21.5,21.5,21.5,21.5,21.5,21.5,2
1.5,21.5,
21.5,21.5,21.5,21.5,21.5,21.5,21.5,21.5,21.5,21.5,21.5,21.5,21.5,21.5,21.5,2
1.5,27.5,
27.5,27.5,27.5,27.5,27.5,27.5,27.5,27.5,27.5,27.5,27.5,27.5,27.5,27.5,27.5,2
7.5,27.5,
27.5,27.5,27.5,27.5,27.5,27.5,27.5,27.5,27.5,27.5,27.5,27.5,27.5,27.5,27.5,2
7.5,27.5,
27.5,27.5,27.5,27.5,27.5,27.5,27.5,27.5,27.5,27.5,27.5,27.5,27.5,33.5,33.5,3
3.5,33.5,
33.5,33.5,33.5,33.5,33.5,33.5,33.5,33.5,33.5,33.5,33.5,33.5,33.5,33.5,33.5,3
3.5,33.5,
33.5,33.5,33.5,33.5,33.5,33.5,33.5,33.5,33.5,33.5,33.5,33.5,33.5,33.5,33.5,3
3.5,33.5,
33.5,33.5,33.5,33.5,33.5,33.5,33.5,33.5,33.5,33.5,33.5,33.5,33.5,39.5,39.5,3
9.5,39.5,
39.5,39.5,39.5,39.5,39.5,39.5,39.5,39.5,39.5,39.5,39.5,39.5,39.5,39.5,39.5,3
9.5,39.5,
39.5,39.5,39.5,39.5,39.5,39.5,39.5,39.5,39.5,39.5,39.5,39.5,39.5,39.5,39.5,3
9.5,39.5,
46,46,46,46,46,46,46,46,46,46,46,46,46,46,46,46,46,46,46,46,46,46,51.5,51.5,
51.5,51.5,
51.5,51.5,51.5,51.5,51.5,51.5,51.5),
y=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,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,2,2,3,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,2,2,2,2,2,2,3,3,3,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,2,2,2,2,2,3,3,3,3,3,3,3,3,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,2,2,2,2,2,2,2,2,2,2,3,3,3,
3,3,3,3,
3,3,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,3,3,3,3,3,3,
3,3,1,1,
1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,3,3,3,3,1,1,1,1,2,2,3,3,3,3,3),
N=365)
Inits
list(beta=0.0, theta=c(1,2))
Thank you to everyone who took the time to help me.
Darren Johnson
Research Scientist I
IAP World Services, Inc
USGS National Wetlands Research Center
700 Cajundome Blvd.
Lafayette, La. 70506
(337)266-8675
(337)266-8568 Fax
[log in to unmask]
-------------------------------------------------------------------
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
|