Having Trouble translating Winbugs code

I’m trying to translate the code below to pymc. It seems doable. I added a data file below. I modified the data so that I’m always looking at two arm comparisons. The tricky part for me is the SW term that is based on lagged values of W. The hard part is to make the the delta matrix, and I think it needs to be created as the output of the pytensor scan function.

It seems to me that this model code should already exist, but I can’t find it. For now, I don’t want to change the basic structure since it was created by an agency that regulates drugs submissions, so before changing it, I want to be able to replicate what they did.

If I ignore the correction for cases that have three or more treatment arms, the pymc code below produces the correct results.

Three Arm Analysis.csv (3.9 KB)

# Binomial likelihood, logit link
# Random effects model for multi-arm trials
model{                                       # *** PROGRAM STARTS
for(i in 1:ns){                              # LOOP THROUGH STUDIES
    w[i,1] <- 0                              # adjustment for multi-arm trials is zero for control arm
    delta[i,1] <- 0                          # treatment effect is zero for control arm
    mu[i] ~ dnorm(0,.0001)                   # vague priors for all trial baselines
    for (k in 1:na[i]) {                     # LOOP THROUGH ARMS
        r[i,k] ~ dbin(p[i,k],n[i,k])         # binomial likelihood
        logit(p[i,k]) <- mu[i] + delta[i,k]  # model for linear predictor
    for (k in 2:na[i]) {                                  # LOOP THROUGH ARMS
        delta[i,k] ~ dnorm(md[i,k],taud[i,k])             # trial-specific LOR distributions
        md[i,k] <-  d[t[i,k]] - d[t[i,1]] + sw[i,k]       # mean of LOR distributions (with multi-arm trial correction)
        taud[i,k] <- tau *2*(k-1)/k                       # precision of LOR distributions (with multi-arm trial correction
        w[i,k] <- (delta[i,k] - d[t[i,k]] + d[t[i,1]])    # adjustment for multi-arm RCTs
        sw[i,k] <- sum(w[i,1:k-1])/(k-1)                  # cumulative adjustment for multi-arm trials             
d[1]<-0       # treatment effect is zero for reference treatment
# vague priors for treatment effects
for (k in 2:nt){  d[k] ~ dnorm(0,.0001) }
sd ~ dunif(0,2)     # vague prior for between-trial SD
tau <- pow(sd,-2)   # between-trial precision = (1/between-trial variance)

If I ignore the third arm, the code below produced the correct results:

import arviz as az
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
import pymc as pm
from pytensor import tensor as pt
from pytensor.ifelse import ifelse
import pytensor 
%matplotlib inline
[Three Arm Analysis.csv|attachment](upload://vTXKrak4cpvJxddxOg8FRGZIzuf.csv) (3.9 KB)

N=len(trials) #
coords = {"studies": studies, "treatments":treatments}

with pm.Model(coords=coords) as femodel:

    fedata=pm.sample(10000,tune=5000, nuts_sampler="nutpie")

az.summary(fedata,var_names=["d"], hdi_prob=0.95)

	mean	sd	hdi_2.5%	hdi_97.5%	mcse_mean	mcse_sd	ess_bulk	ess_tail	r_hat
d[SK]	0.000	0.000	0.000	0.000	0.000	0.000	40000.0	40000.0	NaN
d[t-PA]	-0.003	0.031	-0.064	0.056	0.000	0.000	25876.0	30064.0	1.0
d[at-PA]	-0.158	0.044	-0.247	-0.075	0.000	0.000	10723.0	19178.0	1.0
d[SK+t-PA]	0.453	0.697	-0.926	1.825	0.005	0.004	20350.0	21758.0	1.0
d[r-PA]	-0.112	0.060	-0.228	0.006	0.001	0.000	13126.0	21698.0	1.0
d[TNK]	-0.153	0.077	-0.304	-0.001	0.001	0.000	12869.0	20943.0	1.0
d[PTCA]	-0.475	0.102	-0.677	-0.278	0.001	0.000	21433.0	28942.0	1.0

I may be heading down a dead end. The way i’m going using scan I will need to create pymc random variables inside the scan function to create the “delta” array. However, in browsing I found a comment that said you can’t do this.

It might still work if I could create the random outside of scan but then change their parameters inside of scan.

You can do that, or you can use .dist and wrap the scan in pm.CustomDist. See here.

I assume na[i] (number of arms) differs for each study. Otherwise the two loops in the winbugs can probably be vectorized easily. Maybe what you can do is you can split your dataset into groups where there are data for only two arms and data for all three arms and then write one model for each where now the loop can be vectorized and you don’t need to use scan. Sort of the like the first approach here

Thank you for the suggestions.

What led me stuck in scan was the statement: "sw[i,k] ← sum(w[i,1:k-1])/(k-1) "

In winbugs k=1,2,…number of arms (mostly 2, but could be more, but usually not more than 4 or 5). I thought that I could use scan’s ability to track an accumulating variable. However, I think that I need to do the statement: "delta[i,k] ~ dnorm(md[i,k],taud[i,k]) " in the scan loop, and it’s my understanding that I’m breaking a rule since you can’t create new pymc variables in scan.

I started down a new path and I was thinking about calculating delta + mu for the first and second treatments, and then adding a third based on the first N(number of studies with three arms) from the matrices that I already calculated. It will make the code look a little clunky but not terrible.

Thank you for the suggestions.

What led me to get stuck in scan was the statement: "sw[i,k] ← sum(w[i,1:k-1])/(k-1) "

If sw is a fully defined matrix (for every i and k) then you can probably try pytensor.tensor.extra_ops.cumsum in


if it does the same thing as numpy (I bet it does). And also create a constant matrix whose kth column is 1/(k-1) to divide by in the end.

If I did not misunderstand your problem though, it seems that if you try to do all the data together sw may not be completely defined (i.e would be like a ragged array) so you may need to split your data in to cases where only covariates and data for two arms are defined and cases where covariates and data for three arms are defined. If you are not just working with two categories that I defined but many possible categories of how these covariates and data group, then scan might indeed be better.


Thanks. I can try this. I was a little concerned that I might not be getting a new random variable for each iteration of the scan loop and maybe the same one was being re-used. I did an eval and the random variables had the same value. I will try your suggestion first. I’m concerned that the model might be doing something other than what I think and hope it should be doing. I guess that one check is if it produces a correct answer.


The tricky part for the SW variable is that it accumulates for treatments 2,3,… but always resets to zero when you when you are on a row for treatment 2. I will look at pytensor.tensor.extra_ops.cumsum.

Thanks, some very helpful suggestions.

.eval() doesn’t update the random state between runs, so you are supposed to always get the same thing. You can use pm.draw instead if you want to test the distribution.

I’m making some progress and I am close. I still have some missing pieces since I did not understand the inner outer procedure calls in the example, and what I need to feed back to the scan function. I think below that i’m calculating delta correctly, but I think I missed the part that allows the function to resample when it is used. Also, when I types out.eval(), I got numbers that made sense, but I got errors when I tried delta_.eval(). I’m not sure if this is because I have missing parts, or that eval just does not work on a custom distribution set up the way I’m doing it.

with pm.Model(coords=coords) as femodel:

    def pairs(sid,dd,W,TAU):    
        DELTA=pm.Normal.dist(mu=MD, sigma=TAUD)
        return DELTA

You’re not using CustomDist correctly. Check the docs and the example I linked. You should wrap your scan in a function and pass that function to the dist argument. Inputs to dist are args of CustomDist. Your scan also needs to track random state updates with collect_default_updates, and you probably want to return W, not DELTA.

Hi, I think I have a functioning model:

The Winbugs model produced the following output:


This output was produced by my pymc code:


I have some issues with divergences. Afterall all I am asking the model to estimate 6 fixed effect differences between drug treatments as the main estimate of interest along with 36 means, and 37 random treatment effects. The items in the model look similar to winbugs. My sd is much smaller than what winbugs computed (0.08462 winbugs and 0.008 with my code).

Here is the code. It’s looks a little rudimentary but it seems to mostly work now. I posed the data for three arms earlier.







coords = {"studies": STUDIES, "treatments":TREATMENTS}

with pm.Model(coords=coords) as remodel:

    # set priors


    # Set mu (same for all sets of comparisons)

    # Arm 1    
    # Arm 2
    # Arm 3

    redata=pm.sample(10000,tune=5000, nuts_sampler="nutpie")

You can probably improve you priors to help with divergences. Uniform priors aren’t recommend, because they introduce hard discontinuities into the posterior geometry. You also don’t need to use the tau parameterization as in winBUGS. Sometimes it’s better, but sometimes it’s not. At any rate, you can just directly define a prior over the quantity you care about. If you’re parameterizing with tau, there’s no need to first sample sigma then compute tau. I’d also make the priors more informative, your d parameter has a prior standard deviation of 10,000… why? I’m sure the prior predictive draws will look ludicrous.

General advice would be to do prior predictive sampling to make sure your priors result in scientifically plausible data. For posterior checking, I’d look at az.plot_pair with divergences=True, and look for deterministic relationships, strong correlations, corners/boundaries, and clusters of divergences.

I think the priors look silly, but I’m trying to replicate a model published by a regulatory agency in 2016. As step 1, I am trying to to replicate what they did. It would be great if the owners the PYMC would add an example page with examples taken from the publication in the link below. I don’t see anything on the site about network meta analysis. It’s another domain that they may want to cover. It’s a niche, but sometimes pharmaceuticals want to use these models when they make submissions to regulators in Canada and the UK.

I found another mistake in what I did.


Mu is an array with a shape of 36 x 1, and I want a shape of 1 x 1.

So I am replacing the code with “taud3=pm.Deterministic(“taud3”,pm.math.zeros_like(t3)+tau2(k3-1)/k3)” I will edit the last post that I made to correct it.

The new results look better, especially the sd value.


You also don’t need to wrap everything in Deterministics, only if you need to access them in the trace object

Thanks. I’m new to using pymc so I’m sure there is a better way to construct the model. I kept on getting error messages, so I created an example in excel and then tried to use it as a map. For someone who is new it’s a little challenging to know what is declarative and how to get the computational graph to do what you want it too. I will try to back out the deterministic variables so the model won’t have to track so many.

The next example that I’m going to try is similar and in the publication what I referenced. It has one difference in the the average length of follow-up is given in years. The time is used as an offset to try to normalize the time at risk in each study. The example uses a complementary log-log link, and I found code on the Bambi site, but is would be useful if it was added to pymc.math.

I tried taking out the Deterministic variables and I generated a lot of error messages. It there a general principal of when you you should create a deterministic variable and when not to.

In my case I’m not working data taking from publications and I will never be dealing with “big data”. With what I’m doing speed and efficiency are not critical.

ERROR (pytensor.graph.rewriting.basic): Rewrite failure due to: constant_folding
ERROR (pytensor.graph.rewriting.basic): node: Check{n >= 0, 0 <= p <= 1}([[nan] [na ... an] [nan]], False)
ERROR (pytensor.graph.rewriting.basic): TRACEBACK:
ERROR (pytensor.graph.rewriting.basic): Traceback (most recent call last):
  File "C:\Users\amoyn\anaconda3\Lib\site-packages\pytensor\graph\rewriting\basic.py", line 1922, in process_node
    replacements = node_rewriter.transform(fgraph, node)
  File "C:\Users\amoyn\anaconda3\Lib\site-packages\pytensor\graph\rewriting\basic.py", line 1082, in transform
    return self.fn(fgraph, node)
  File "C:\Users\amoyn\anaconda3\Lib\site-packages\pytensor\tensor\rewriting\basic.py", line 1106, in constant_folding
    required = thunk()
  File "C:\Users\amoyn\anaconda3\Lib\site-packages\pytensor\link\c\op.py", line 91, in rval
  File "C:\Users\amoyn\anaconda3\Lib\site-packages\pytensor\link\c\basic.py", line 1792, in __call__
    raise exc_value.with_traceback(exc_trace)
pymc.logprob.utils.ParameterValueError: n >= 0, 0 <= p <= 1

... omitted 

 0.00% [0/44000 00:00<?]
ValueError                                Traceback (most recent call last)
Cell In[27], line 38
     35 p3=pm.Deterministic("p3",pm.invlogit(mu[s3]+delta3))
     36 y3=pm.Binomial("y3",n=N3,p=p3,observed=R3)
---> 38 redata=pm.sample(10000,tune=1000, nuts_sampler="nutpie")

File ~\anaconda3\Lib\site-packages\pymc\sampling\mcmc.py:674, in sample(draws, tune, chains, cores, random_seed, progressbar, step, nuts_sampler, initvals, init, jitter_max_retries, n_init, trace, discard_tuned_samples, compute_convergence_checks, keep_warning_stat, return_inferencedata, idata_kwargs, nuts_sampler_kwargs, callback, mp_ctx, model, **kwargs)
    670     if not isinstance(step, NUTS):
    671         raise ValueError(
    672             "Model can not be sampled with NUTS alone. Your model is probably not continuous."
    673         )
--> 674     return _sample_external_nuts(
    675         sampler=nuts_sampler,
    676         draws=draws,
    677         tune=tune,
    678         chains=chains,
    679         target_accept=kwargs.pop("nuts", {}).get("target_accept", 0.8),
    680         random_seed=random_seed,
    681         initvals=initvals,
    682         model=model,
    683         progressbar=progressbar,
    684         idata_kwargs=idata_kwargs,
    685         nuts_sampler_kwargs=nuts_sampler_kwargs,
    686         **kwargs,
    687     )
    689 if isinstance(step, list):
    690     step = CompoundStep(step)

File ~\anaconda3\Lib\site-packages\pymc\sampling\mcmc.py:297, in _sample_external_nuts(sampler, draws, tune, chains, target_accept, random_seed, initvals, model, progressbar, idata_kwargs, nuts_sampler_kwargs, **kwargs)
    295 compiled_model = nutpie.compile_pymc_model(model)
    296 t_start = time.time()
--> 297 idata = nutpie.sample(
    298     compiled_model,
    299     draws=draws,
    300     tune=tune,
    301     chains=chains,
    302     target_accept=target_accept,
    303     seed=_get_seeds_per_chain(random_seed, 1)[0],
    304     progress_bar=progressbar,
    305     **nuts_sampler_kwargs,
    306 )
    307 t_sample = time.time() - t_start
    308 # Temporary work-around. Revert once https://github.com/pymc-devs/nutpie/issues/74 is fixed
    309 # gather observed and constant data as nutpie.sample() has no access to the PyMC model

File ~\anaconda3\Lib\site-packages\nutpie\sample.py:236, in sample(compiled_model, draws, tune, chains, cores, seed, save_warmup, progress_bar, init_mean, return_raw_trace, **kwargs)
    234         pass
    235 finally:
--> 236     results = sampler.finalize()
    238 dims = {name: list(dim) for name, dim in compiled_model.dims.items()}
    239 dims["mass_matrix_inv"] = ["unconstrained_parameter"]

ValueError: Sampling failed: All initialization points failed

Caused by:
    Logp function returned error: Logp function returned error code 1

Deterministics are never needed, I suspect you introduced an error when you removed them.

c = pm.Deterministic("x", a + b)
d = pm.Normal("d", mu=c)

Can simply be written as

c = a + b
d = pm.Normal("d", mu=c)

Which is rather less verbose.

The “rule” for knowing when to use it is: do you care about this specific intermediate computation? That means you either want to see its values in the trace or for it to show in model.to_graphviz()