I attach herewith, the code and a sample of the dataset.
Any suggestion will be greatly appreciated.
Chris Guure
Biostatistics Department
University of Ghana
library(R2WinBUGS)
library(coda)model1<-function(){
for (i in 1:N) {
# likelihood function
ms[i] ~ dbin( p [ i ], N )
logit(p [ i ] ) <- alpha + bage*age[ i ] + bpam*pam[i ] + bpah*pah[i]
}
### prior for intercept
alpha ~ dnorm(0,0.0001)
# prior for slopes
bage ~ dnorm(0,0.0001)
bpam ~ dnorm(0,0.0001)
bpah ~ dnorm(0,0.0001)
# OR for alpha
or.age<-exp(bage)
# OR for hbp
or.pam <- exp(bpam)
# OR for fdm
or.pah <- exp(bpah)
}
data=cbind(ms=c(0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0),
age=c(77, 83, 75, 78, 75, 83, 85, 80, 80, 85, 76, 77, 80, 76, 88, 77, 81, 78,
85, 81),
pam=c(0, 0, 1, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1),
pah=c(1, 1, 0, 0, 1, 0, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0))
N=20
#### initial values
###
lineinits <- function(){
list(alpha=1,bage= 0.05, bpam=0.05, bpah=0.05)
}
## CODA output
lineout1 <- bugs(data, lineinits, c("bage", "alpha","bpam", "bpah", "or.age", "or.pam",
"or.pah"), model1,n.iter = 11000, n.burnin = 1000, n.chains = 2, codaPkg = T, DIC = TRUE)
### Posterior summaries
line.coda <- read.bugs(lineout1)
summary(line.coda)
______________________________________________
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.