Here is one more example. On the first example there was only one study with three arms. This one has 4 studies with three arms and the estimates are a lot cleaner. In the example there is a variable for time to adjust for differential follow-up, but I need a inverse complimentary log log link to implement. Ignoring time, the output produces the plots below that look pretty clean.
Code below:
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
import pandas as pd
import io
data_string="""studyid,time,t1,t2,t3,r1,r2,r3,n1,n2,n3,study
0,5.8,1,2,3,43,34,37,1081,2213,1102,1.MRC-E
1,4.7,1,2,,29,20,,416,424,,2.EWPH
3,3,1,2,,140,118,,1631,1578,,3.SHEP
4,3.8,1,3,,75,86,,3272,3297,,4.HAPPHY
5,4,1,4,5,302,154,119,6766,3954,4096,5.ALLHAT
6,3,1,4,,176,136,,2511,2508,,6.INSIGHT
7,4.1,1,5,,200,138,,2826,2800,,7.ANBP-2
8,1,1,6,,8,1,,196,196,,8.ALPINE
9,3.3,2,4,,154,177,,4870,4841,,9.FEVER
10,3,2,5,,489,449,,2646,2623,,10.DREAM
11,4.5,2,5,,155,102,,2883,2837,,11.HOPE
12,4.8,2,5,,399,335,,3472,3432,,12.PEACE
13,3.1,2,6,,202,163,,2721,2715,,13.CHARM
14,3.7,2,6,,115,93,,2175,2167,,14.SCOPE
15,3.8,3,4,5,70,32,45,405,202,410,15.AASK
16,4,3,4,5,97,95,93,1960,1965,1970,16.STOP-2
17,5.5,3,4,,799,567,,7040,7072,,17.ASCOT
18,4.5,3,4,,251,216,,5059,5095,,18.NORDIL
19,4,3,4,,665,569,,8078,8098,,19.INVEST
20,6.1,3,5,,380,337,,5230,5183,,20.CAPPP
21,4.8,3,6,,320,242,,3979,4020,,21.LIFE
22,4.2,4,6,,845,690,,5074,5087,,22.VALUE"""
data=pd.read_csv(io.StringIO(data_string),sep=",")
N=len(data)
N1=data[["n1"]].astype(dtype=np.int32).to_numpy()
N2=data[["n2"]].astype(dtype=np.int32).to_numpy()
N3=data[["n3"]].dropna().astype(dtype=np.int32).to_numpy()
R1=data[["r1"]].astype(dtype=np.int32).to_numpy()
R2=data[["r2"]].astype(dtype=np.int32).to_numpy()
R3=data[["r3"]].dropna().astype(dtype=np.int32).to_numpy()
T1=data[["t1"]].astype(dtype=np.int32).to_numpy()-1
T2=data[["t2"]].astype(dtype=np.int32).to_numpy()-1
T3=data[["t3"]].dropna().astype(dtype=np.int32).to_numpy()-1
TIME=data["time"].to_numpy()
S3=data.dropna(subset=["t3"])["studyid"]
K=np.array([1,2,3])
STUDIES=data["study"]
TREATMENTS=["Diuretic","Placebo","Beta Blocker","CCB","ACE inhibitor","ARB"]
coords = {"studies": STUDIES, "treatments":TREATMENTS}
with pm.Model(coords=coords) as remodel:
# set priors
d_=pm.Normal("d_",mu=0,tau=0.0001,shape=(5,))
d=pm.Deterministic("d",pm.math.concatenate([[pt.constant(0)],d_]),dims="treatments")
sd=pm.Uniform("sd",0,5)
tau=1/(sd*sd)
# Set mu (same for all sets of comparisons)
mu=pm.Normal("mu",mu=0,tau=0.0001,dims="studies",shape=(N,1))
# Arm 1
r1=pm.ConstantData('r1',R1)
n1=pm.ConstantData('n1',N1)
t1=pm.ConstantData('t1',T1)
p1=pm.invlogit(mu)
y1=pm.Binomial("y1",n=n1,p=p1,observed=r1)
# Arm 2
r2=pm.ConstantData('r2',R2)
n2=pm.ConstantData('n2',N2)
t2=pm.ConstantData('t2',T2)
md2=d[t2]-d[t1]
taud2=pm.math.zeros_like(mu)+tau
delta2=pm.Normal("delta2",mu=md2,tau=taud2)
p2=pm.invlogit(mu+delta2)
y2=pm.Binomial("y2",n=n2,p=p2,observed=r2)
w2=delta2-d[t2]+d[t1]
# Arm 3
r3=pm.ConstantData('r3',R3)
n3=pm.ConstantData('n3',N3)
t3=pm.ConstantData('t3',T3)
s3=pm.ConstantData('s3',S3)
k3=pt.constant(3)
sw3=(w2[s3])/(3-1)
taud3=pm.math.zeros_like(t3)+tau*2*(3-1)/3
md3=d[t3]-d[t1[s3]]+sw3
delta3=pm.Normal("delta3",mu=md3,tau=taud3);
p3=pm.invlogit(mu[s3]+delta3)
y3=pm.Binomial("y3",n=n3,p=p3,observed=r3)
redata=pm.sample(100000,tune=10000, nuts_sampler="nutpie", random_seed=8675309)
az.plot_trace(redata,var_names=["d","sd"])
az.summary(redata,var_names=["d","sd"], hdi_prob=0.95)
