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)
trials=data[["n1","n2"]]
events=data[["r1","r2"]]
trt=data[["t1","t2"]].to_numpy()-1
N=len(trials) #
studies=data["study"]
treatments=["SK","t-PA","at-PA","SK+t-PA","r-PA","TNK","PTCA"]
coords = {"studies": studies, "treatments":treatments}
reference=np.array([0,1,1,1,1,1,1])
with pm.Model(coords=coords) as femodel:
r=pm.MutableData('r',events)
n=pm.MutableData('n',trials)
t=pm.MutableData('t',trt)
d_=pm.Normal("d_",mu=0,sigma=100,shape=(6,))
d=pm.Deterministic("d",pm.math.concatenate([[pt.constant(0)],d_]),dims="treatments")
mu=pm.Normal("mu",mu=0,sigma=100,dims="studies",shape=(N,1))
dm=pm.Deterministic("dm",d[t])
linpred=pm.Deterministic("linpred",mu+dm)
p=pm.Deterministic("p",pm.invlogit(linpred))
y=pm.Binomial("y",n=n,p=p,observed=r)
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
