In Bayesian statistics and probabilistic modeling, \(\eta\) typically represents the linear predictor or the systematic component of a regression model. In the context of the provided PyMC3 model, \(\eta\) represents the linear predictor for each observation.

Specifically, in the provided model:

\[ \begin{align*} \eta_{ij} &= \frac{{\phi_{i0j}}}{{1 + \phi_{i1j} \cdot \exp(\phi_{i2j} \cdot x_j)}} \end{align*} \]

Here, \(i\) represents the group index (from 1 to \(K\)), \(j\) represents the observation index (from 1 to \(n\)), and \(x_j\) represents the \(j\)th predictor variable.

Breaking it down:

So, \(\eta_{ij}\) is the linear predictor for the \(j\)th observation in the \(i\)th group. It’s calculated based on the predictor variable \(x_j\) and the parameters \(\phi_{i0j}\), \(\phi_{i1j}\), and \(\phi_{i2j}\).

In summary, \(\eta\) serves as the linear predictor in the hierarchical regression model, where its values are determined by the observed predictor variables and the parameters of the model.

library('rjags')
## Loading required package: coda
## Linked to JAGS 4.3.2
## Loaded modules: basemod,bugs
library('coda')

# Example: Define or load the data Y, covariate x, and number of observations n
# You can replace these with your actual data
Y <- matrix(rnorm(100), nrow = 10, ncol = 10)  # Example: 10 groups, 10 observations per group
x <- rnorm(10)  # Example: Covariate x with 10 values
n <- ncol(Y)  # Example: Number of observations per group

# Define the number of groups (K) based on your data
K <- nrow(Y)

# Mathematical Model:
# We have a hierarchical Bayesian model with parameters \(\mu\), \(\tau\), \(\theta\), and \(\phi\).
# The likelihood function is a normal distribution for each observation \(Y_{ij}\) given \(\eta_{ij}\) with precision \(\tau_C\).
# The linear predictor \(\eta_{ij}\) is calculated based on the covariate \(x_j\) and the parameters \(\phi_i\).
# The parameters \(\phi_i\) are transformed from the parameters \(\theta_i\) using exponential functions.
# Both \(\mu\) and \(\tau\) have priors as normal and gamma distributions, respectively.
# The precision \(\tau_C\) also has a gamma prior.


# Model specification
model_string <- "
model {
    for (i in 1:K) {
        for (j in 1:n) {
            Y[i, j] ~ dnorm(eta[i, j], tauC)  # Likelihood
            eta[i, j] <- phi[i, 1] / (1 + phi[i, 2] * exp(phi[i, 3] * x[j]))  # Linear predictor
        }
        phi[i, 1] <- exp(theta[i, 1])
        phi[i, 2] <- exp(theta[i, 2]) - 1
        phi[i, 3] <- -exp(theta[i, 3])
        for (k in 1:3) {
            theta[i, k] ~ dnorm(mu[k], tau[k])  # Prior for theta
        }
    }
    
    tauC ~ dgamma(1.0E-3, 1.0E-3)  # Prior for tauC
    sigmaC <- 1 / sqrt(tauC)
    varC <- 1 / tauC
    
    for (k in 1:3) {
        mu[k] ~ dnorm(0, 1.0E-4)  # Prior for mu
        tau[k] ~ dgamma(1.0E-3, 1.0E-3)  # Prior for tau
        sigma[k] <- 1 / sqrt(tau[k])
    }
}
"

# Initial values for parameters
inits <- list(
    mu = rep(0, 3),
    tau = rep(1, 3),
    theta = matrix(0, nrow = K, ncol = 3)
)

# Parameters to monitor
params <- c("mu", "tau", "theta")

# Running the model
model <- jags.model(textConnection(model_string), data = list(Y = Y, x = x, n = n, K = K), inits = inits)
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 100
##    Unobserved stochastic nodes: 37
##    Total graph size: 712
## 
## Initializing model
update(model, 1000)
samples <- coda.samples(model, variable.names = params, n.iter = 5000)

# Plotting the MCMC chains
par(mfrow=c(3,1), mar=c(4,4,2,1), oma=c(0,0,2,0))
plot(samples, ask=FALSE)

library('ggplot2')

# Define a range of values for tauC
tauC_values <- seq(0, 2, by = 0.01)

# Define the PDF for the Gamma prior
gamma_prior <- dgamma(tauC_values, shape = 0.001, rate = 0.001)

# Define the PDF for the Uniform prior
uniform_prior <- dunif(tauC_values, min = 0, max = 2)

# Create a data frame for plotting
plot_data <- data.frame(tauC = tauC_values, gamma_prior = gamma_prior, uniform_prior = uniform_prior)

# Plot the densities
ggplot(plot_data, aes(x = tauC)) +
  geom_line(aes(y = gamma_prior, color = "Gamma Prior"), size = 1) +
  geom_line(aes(y = uniform_prior, color = "Uniform Prior"), size = 1) +
  labs(x = expression(tau[C]), y = "Density", color = "Prior Distribution") +
  theme_minimal()
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

dgamma(1,1.0E-3, 1.0E-3)
## [1] 0.0009926954
dgamma(2,1.0E-3, 1.0E-3)
## [1] 0.0004961954
plot((density(dgamma(seq(0,2,0.01),1.0E-3, 1.0E-3))))