19. Paraguay Vaccination Status*

Contents 

import pymc as pm
import pandas as pd
import numpy as np
import arviz as az

rng = np.random.default_rng(6420)

19. Paraguay Vaccination Status*#

This example goes over a multilevel, or hierarchical, logistic regression model. It also shows how to use the PyMC coordinate system.

Adapted from Unit 7: paraguay.odc and paraguaynocluster.odc

Data can be found here.

This study considers factors influencing the vaccination status among 3424 children of 2552 mothers among 264 clusters in Paraguay. In this analysis, we’re specifically interested in mother-level factors related to child immunization. However, there is randomness associated with different clusters.

  • ID3: Cluster number

  • VACCODE: =1 if fully immunized, =0 otherwise

  • LB.TOT: No. of live births

  • MAGE2: mother age <20 =1, otherwise = 0

  • UN2: consensual union = 1, otherwise = 0

  • TOILET2: unsafe toilet in hh = 1, otherwise = 0

  • PR.SPOC1: spouse unskilled laborer = 1, otherwise = 0

  • SPANISH2: Spanish not hh language = 1, otherwise = 0

We need to add a random effect by cluster. This is a good use case for the PyMC coordinates system.

data = pd.read_csv("../data/paraguay.csv")
data.info()

<class 'pandas.DataFrame'>
RangeIndex: 3424 entries, 0 to 3423
Data columns (total 8 columns):
 #   Column    Non-Null Count  Dtype
---  ------    --------------  -----
 0   ID3       3424 non-null   int64
 1   VACCODE   3424 non-null   int64
 2   LB.TOT    3424 non-null   int64
 3   MAGE2     3424 non-null   int64
 4   UN2       3424 non-null   int64
 5   TOILET2   3424 non-null   int64
 6   PR.SPOC1  3424 non-null   int64
 7   SPANISH2  3424 non-null   int64
dtypes: int64(8)
memory usage: 214.1 KB

y = data["VACCODE"].to_numpy()
# separate array for clusters
clusters = data["ID3"].to_numpy()
X = data.drop(["VACCODE", "ID3"], axis=1).to_numpy()
X_aug = np.concatenate((np.ones((X.shape[0], 1)), X), axis=1)
cols = X_aug.shape[1]
y.shape, clusters.shape, X_aug.shape

((3424,), (3424,), (3424, 7))

# set up alternate coordinates, the ID3 or clusters column
cluster_idx, clusters = pd.factorize(data.ID3)
coords = {"cluster": clusters, "id": data.index.to_numpy()}

# note that the coords dict is passed to pm.Model call
with pm.Model(coords=coords) as m:
    X_data = pm.Data("X_data", X_aug)
    y_data = pm.Data("y_data", y)
    clust_idx = pm.Data("cluster_idx", cluster_idx, dims="id")

    cluster_tau = pm.Gamma("cluster_tau", 0.01, 0.01)
    cluster_variance = pm.Deterministic("cluster_variance", 1 / cluster_tau)
    beta = pm.Normal("beta", 0, tau=1e-3, shape=cols)

    cluster_effect = pm.Normal(
        "cluster_effect", 0, tau=cluster_tau, dims="cluster"
    )
    logit_p = pm.math.dot(X_data, beta) + cluster_effect[clust_idx]

    pm.Bernoulli("likelihood", logit_p=logit_p, observed=y_data)

    trace = pm.sample(3000, random_seed=rng)

Hide code cell output

 
NUTS[nutpie]: [cluster_tau, beta, cluster_effect]













az.summary(trace, var_names=["beta", "cluster_variance"], kind="stats")










  
    

      
      mean
      sd
      eti89_lb
      eti89_ub
    

  
  
    

      beta[0]
      1.4
      0.13
      1.2
      1.7
    

    

      beta[1]
      -0.07
      0.015
      -0.093
      -0.047
    

    

      beta[2]
      -0.56
      0.2
      -0.89
      -0.24
    

    

      beta[3]
      -0.2
      0.095
      -0.35
      -0.044
    

    

      beta[4]
      -0.7
      0.14
      -0.91
      -0.48
    

    

      beta[5]
      -0.29
      0.11
      -0.46
      -0.11
    

    

      beta[6]
      -0.62
      0.098
      -0.78
      -0.47
    

    

      cluster_variance
      0.53
      0.094
      0.39
      0.68
    

  





This was originally based on this old coordinates example.




No clusters#






with pm.Model() as m_nc:
    X_data = pm.Data("X_data", X_aug)
    y_data = pm.Data("y_data", y)

    beta = pm.Normal("beta", 0, tau=1e-3, shape=cols)

    logit_p = pm.math.dot(X_data, beta)

    pm.Bernoulli("likelihood", logit_p=logit_p, observed=y_data)

    trace_nc = pm.sample(3000, random_seed=rng)











Hide code cell output





NUTS[nutpie]: [beta]

















az.summary(trace_nc, var_names=["beta"], kind="stats")










  
    

      
      mean
      sd
      eti89_lb
      eti89_ub
    

  
  
    

      beta[0]
      1.4
      0.11
      1.3
      1.6
    

    

      beta[1]
      -0.064
      0.014
      -0.085
      -0.042
    

    

      beta[2]
      -0.57
      0.19
      -0.87
      -0.27
    

    

      beta[3]
      -0.19
      0.087
      -0.33
      -0.047
    

    

      beta[4]
      -0.72
      0.11
      -0.9
      -0.54
    

    

      beta[5]
      -0.45
      0.088
      -0.59
      -0.31
    

    

      beta[6]
      -0.61
      0.084
      -0.74
      -0.47
    

  









%load_ext watermark
%watermark -n -u -v -iv -p pytensor









Last updated: Sun, 18 Oct 2026

Python implementation: CPython
Python version       : 3.14.7
IPython version      : 9.16.1

pytensor: 3.3.0

arviz : 1.3.0
numpy : 2.4.6
pandas: 3.0.5
pymc  : 6.3.1













    
    

                

              

              
              
              
              
                
              
            

            
            
              
                
                



  

     Contents