Clinical Trial data - request help on variable number of comparisons per trial

Questions v5
7

That is ok, from the text file I was able to get the excel and run the code. There are some things which does not make sense to start with:

You create a variable “d_” with seven elements. Two things I don’t understand is later you transform it via

d=pm.Deterministic("d",d_*reference,dims="treatments")

which only has the effect of making d[0] 0 and keeping other d the same. Due to the structure of your model later on, d_[0] therefore does not contribute to anything. On top of this you later on do:

dm=pm.Deterministic("dm",d[t])

and you go on to use. However t is an indexer that simply picks the last and first element of d and tiles it, therefore anything other than the first and last element gets discarded basically (as they are not used elsewhere). With the fact that you set d_ to 0 before, basically only thing that matters for your model is d_[-1]. And indeed this is why the posterior plots for d_ look like

posterior
posterior1191×819 60.8 KB

As you can see anything other than d_[-1] is exactly the same as your prior, very large normal distributions. If you only care about the last element, you probably don’t need to create others? For others interested, the code I am using is as follows:

import arviz as az
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
import pymc as pm
import pytensor.tensor as pt
import pytensor 

import pandas as pd
import io
blocker_string="""
study,na,r1,r2,n1,n2,t1,t2
1,2,3,3,39,38,0,1
2,2,14,7,116,114,0,1
3,2,11,5,93,69,0,1
4,2,127,102,1520,1533,0,1
5,2,27,28,365,355,0,1
6,2,6,4,52,59,0,1
7,2,152,98,939,945,0,1
8,2,48,60,471,632,0,1
9,2,37,25,282,278,0,1
10,2,188,138,1921,1916,0,1
11,2,52,64,583,873,0,1
12,2,47,45,266,263,0,1
13,2,16,9,293,291,0,1
14,2,45,57,883,858,0,1
15,2,31,25,147,154,0,1
16,2,38,33,213,207,0,1
17,2,12,28,122,251,0,1
18,2,6,8,154,151,0,1
19,2,3,6,134,174,0,1
20,2,40,32,218,209,0,1
21,2,43,27,364,391,0,1
22,2,39,22,674,680,0,1"""
data=pd.read_csv(io.StringIO(blocker_string),sep=",")

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,dims="treatments")
    d=pm.Deterministic("d",d_*reference,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(1000,tune=500)

pm.model_to_graphviz(femodel)