@Vfeagins Hi Victor, the way you have defined your model is incorrect. When you define a Bayesian model, you are defining that using the parameters specified.
Here your model parameters are theta[1], your intercept, theta[2], your beta coefficient and X your independent variable.
If you examine and read the Winbugs model notes for Air carefully, you will see that Z is the covariate representing the NO2 concentration the bedrooms. In other words Z is the measurement of the Independent variable X with measurement error. That is the reason why you specify
mu <- alpha + beta*Z and in-turn you specify X as
X <- normal(mu, sigma)
In this manner, X has now become your stochastic, independent variable
You are modeling the probability that a child has a respiratory illness using
p <- ilogit(theta[1] + theta[2] * X)
So your model statement has to be:
mod <- model(theta[1], theta[2], X)
If you do not specify precisely the quantity and the input parameters and the covariates (the X variables), greta will not know what to model.
So your complete correct code has to be the following:
library(reticulate)
library(greta)
y <- c(21, 20, 15)
n <- c(48, 34, 21)
Z <- c(10, 30, 50)
alpha <- 4.48
beta <- 0.76
sigma2 <- 81.14
sigma <- sqrt(sigma2)
tau <- 1 / sigma2
J <- 3
theta <- normal(0, 32, dim = 2)
mu <- alpha + beta * Z # this is deterministic
X <- normal(mu, sigma) # with this you make X a random variable
p <- ilogit(theta[1] + theta[2] * X) # this is your Logistic regression
distribution(y) <- binomial(n, p) # this defines the Y, your original dependent variable
mod <- model(theta[1], theta[2], X)
draws <- mcmc(mod, n_samples = 1000)
summary(draws)
# to get a summary table of the parameters and compare to Winbugs output
# the results you get are very similar to the results in the Winbugs output
library(bayesplot)
mcmc_trace(draws)
# to visually inspect the posterior of the parameters
# to visually inspect the how well the draws mixed
# the trace of each variable / parameter should look like a catepillar
If you run correctly, the your output should look similar to what I got below (will not match exactly as this is a Bayesian model)
Iterations = 1:1000
Thinning interval = 1
Number of chains = 4
Sample size per chain = 1000
1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
theta[1] -1.04189 1.16814 0.0184699 0.236719
theta[2] 0.04977 0.03967 0.0006273 0.006324
X[1,1] 13.04112 8.90381 0.1407816 1.736653
X[2,1] 27.53257 7.38022 0.1166915 0.447116
X[3,1] 41.11623 8.56493 0.1354234 0.313373
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
theta[1] -4.728303 -1.35201 -0.71607 -0.31715 0.2348
theta[2] 0.004259 0.02476 0.03896 0.06148 0.1701
X[1,1] -6.672637 7.14481 14.51069 19.26854 28.3531
X[2,1] 13.112142 22.70398 27.43145 32.35168 42.4910
All the best with your Bayesian adventure with greta.
