library(R2OpenBUGS)
library(BRugs)
library(coda)
logist.mod <- function(){
for (i in 1:N) {
r[i] ~ dbin(p[i], n[i]) #asumsi1
b[i] ~ dnorm(0, tau) #asumsi2
logit(p[i]) <- alpha0 + alpha1 * x1[i] + alpha2 * x2[i]
+ alpha12 * x1[i] * x2[i] + b[i]
}
#distribusi prior
alpha0 ~ dnorm(0, 1.0E-6)
alpha1 ~ dnorm(0, 1.0E-6)
alpha2 ~ dnorm(0, 1.0E-6)
alpha12 ~ dnorm(0, 1.0E-6)
tau ~ dgamma(0.001, 0.001)
sigma <- 1 / sqrt(tau)
}
data.1 <- list(r = c(10, 23, 23, 26, 17, 5, 53, 55, 32, 46, 10, 8, 10, 8, 23, 0, 3, 22, 15, 32, 3), #r = germinated
n = c(39, 62, 81, 51, 39, 6, 74, 72, 51, 79, 13, 16, 30, 28, 45, 4, 12, 41, 30, 51, 7), #jumlah total percobaan r
x1 = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1),
x2 = c(0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1),
N = 21, tau=.5)
parameters <- c("alpha0", "alpha1", "alpha2", "alpha12")
inits <- function(){
list(alpha0 = .5, alpha1 = 1,
alpha2 = .2, alpha12 = .3)
}
fit <- bugs(data.1, inits, logist.mod,
parameters = c("alpha0", "alpha1", "alpha2", "alpha12"),
n.chains = 3, n.iter = 1000, codaPkg = T)
codaobject <- read.bugs(fit)
plot(codaobject)
alpha0 <- fit2$sims.list$alpha0
alpha1 <- fit2$sims.list$alpha1
fit2$summary
plot(alpha0,alpha1, pch=c(16,17), col=c('red', 'blue'))
codaobject <- read.bugs(fit)
plot(codaobject)
summary(codaobject)
ld50 <- -alpha0/alpha1
hist(ld50, breaks=50)
quantile(ld50, c(.025, .5, .975))
library(BRugs)
library(coda)
logist.mod <- function(){
for (i in 1:N) {
r[i] ~ dbin(p[i], n[i]) #asumsi1
b[i] ~ dnorm(0, tau) #asumsi2
logit(p[i]) <- alpha0 + alpha1 * x1[i] + alpha2 * x2[i]
+ alpha12 * x1[i] * x2[i] + b[i]
}
#distribusi prior
alpha0 ~ dnorm(0, 1.0E-6)
alpha1 ~ dnorm(0, 1.0E-6)
alpha2 ~ dnorm(0, 1.0E-6)
alpha12 ~ dnorm(0, 1.0E-6)
tau ~ dgamma(0.001, 0.001)
sigma <- 1 / sqrt(tau)
}
data.1 <- list(r = c(10, 23, 23, 26, 17, 5, 53, 55, 32, 46, 10, 8, 10, 8, 23, 0, 3, 22, 15, 32, 3), #r = germinated
n = c(39, 62, 81, 51, 39, 6, 74, 72, 51, 79, 13, 16, 30, 28, 45, 4, 12, 41, 30, 51, 7), #jumlah total percobaan r
x1 = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1),
x2 = c(0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1),
N = 21, tau=.5)
parameters <- c("alpha0", "alpha1", "alpha2", "alpha12")
inits <- function(){
list(alpha0 = .5, alpha1 = 1,
alpha2 = .2, alpha12 = .3)
}
fit <- bugs(data.1, inits, logist.mod,
parameters = c("alpha0", "alpha1", "alpha2", "alpha12"),
n.chains = 3, n.iter = 1000, codaPkg = T)
codaobject <- read.bugs(fit)
plot(codaobject)
alpha0 <- fit2$sims.list$alpha0
alpha1 <- fit2$sims.list$alpha1
fit2$summary
plot(alpha0,alpha1, pch=c(16,17), col=c('red', 'blue'))
codaobject <- read.bugs(fit)
plot(codaobject)
summary(codaobject)
ld50 <- -alpha0/alpha1
hist(ld50, breaks=50)
quantile(ld50, c(.025, .5, .975))
Komentar
Posting Komentar