If you break up the likelihood does it have an impact on the computations

Questions v5
1

I broke my likelihood function into 4 pieces. I could probably have got it down to two, and maybe one if I used dimensions. It works with 4. However, I would like to know if this has any impact on how the sampling works, and if there is an effect on how the parameters in the different likelihoods are correlated. Would the relationship be different if I used just one likelihood.

def Su(lambda_,t):
    return pm.math.exp(-lambda_ * t)

def So(lambda_,t,p):
    return p + (1-p)*Su(lambda_,t)

def fu (lambda_,t):    
    return lambda_*pm.math.exp(-lambda_ * t)

def ho (lambda_,t,p):    
    return (1-p)*fu(lambda_,t) / (p + (1-p)*Su(lambda_,t))


# Data 
dataSOC=osrfs.query('treatment=="SOC"')

T_OS_SOC=dataSOC["os"].to_numpy()
E_OS_SOC=dataSOC["os_event"].to_numpy()

T_RFS_SOC=dataSOC["rfs"].to_numpy()
E_RFS_SOC=dataSOC["rfs_event"].to_numpy()

dataFLIN=osrfs.query('treatment=="SOC FLIN"')

T_OS_FLIN=dataFLIN["os"].to_numpy()
E_OS_FLIN=dataFLIN["os_event"].to_numpy()

T_RFS_FLIN=dataFLIN["rfs"].to_numpy()
E_RFS_FLIN=dataFLIN["rfs_event"].to_numpy()

K=2
J=2

coords={'endpoints': ['OS','RFS'], 'treatments': ['SOC', 'SOC FLIN']}

# Model

from pymc.math import exp, log
with pm.Model() as model_2:
    
    # p global 
    p_global=pm.Beta("p_global",alpha=1,beta=1)
    beta_global=pm.Deterministic("beta_global",pm.math.logit(p_global))
    
    # Priors by treatment
    vk=pm.Normal("vk",mu=beta_global, sigma=0.2,shape=(K,1))  
    sk=pm.HalfNormal("sk",2.5,shape=(K,1))

    lind_SOC=pm.Normal("lind_SOC",-1.0,sigma=5);
    lind_FLIN=pm.Normal("lind_FLIN",-1.0,sigma=5);

    d_SOC=pm.Deterministic("d_SOC",pm.math.invlogit(lind_SOC))
    d_FLIN=pm.Deterministic("d_FLIN",pm.math.invlogit(lind_FLIN))
    
    linp_SOC=pm.Normal("linp_SOC",vk[0],sk[0])
    p_SOC=pm.Deterministic("p_SOC",pm.math.invlogit(linp_SOC))

    linp_FLIN=pm.Normal("linp_FLIN",vk[1],sk[1])
    p_FLIN=pm.Deterministic("p_FLIN",pm.math.invlogit(linp_FLIN))

    beta0_OS_SOC=pm.Normal("beta0_OS_SOC",mu=0,sigma=100)
    lambda_OS_SOC=pm.math.exp(-(beta0_OS_SOC))
    OS_SOC = pm.Potential("OS_SOC",E_OS_SOC*log(ho(lambda_OS_SOC,T_OS_SOC,(p_SOC+d_SOC)))+log(So(lambda_OS_SOC,T_OS_SOC,(p_SOC+d_SOC))))

    beta0_OS_FLIN=pm.Normal("beta0_OS_FLIN",mu=0,sigma=100)
    lambda_OS_FLIN=pm.math.exp(-(beta0_OS_FLIN))
    OS_FLIN = pm.Potential("OS_FLIN",E_OS_FLIN*log(ho(lambda_OS_FLIN,T_OS_FLIN,(p_FLIN+d_FLIN)))+log(So(lambda_OS_FLIN,T_OS_FLIN,(p_FLIN+d_FLIN))))

    beta0_RFS_SOC=pm.Normal("beta0_RFS_SOC",mu=0,sigma=100)
    lambda_RFS_SOC=pm.math.exp(-(beta0_RFS_SOC))
    RFS_SOC = pm.Potential("RFS_SOC",E_RFS_SOC*log(ho(lambda_RFS_SOC,T_RFS_SOC,p_SOC))+log(So(lambda_RFS_SOC,T_RFS_SOC,p_SOC)))

    beta0_RFS_FLIN=pm.Normal("beta0_RFS_FLIN",mu=0,sigma=100)
    lambda_RFS_FLIN=pm.math.exp(-(beta0_RFS_FLIN))
    RFS_FLIN = pm.Potential("RFS_FLIN",E_RFS_FLIN*log(ho(lambda_RFS_FLIN,T_RFS_FLIN,p_FLIN))+log(So(lambda_RFS_FLIN,T_RFS_FLIN,p_FLIN))) 
    
    data = pm.sample(draws=10000,tune=5000,nuts_sampler="nutpie",cores=20,target_accept=0.99)
2

May I ask why did you do it?

Also the custom sampler settings are a bit of a red flag. If you need so extreme values there’s more often than not a problem with your model

3

I am extrapolating out of sample. I don’t want the OS and RFS curves to cross. By definition they can’t. However, they can when you estimate and extrapolate them the SOC arm does. I set up the model so SOC has the same cure fraction fraction for RFS and OS. I added a d to make OS a little higher than RFS. I did the same for FLIN. So, I’m stuck with at least two likelihoods. My data is one record for each patient, and has an OS time, OS event, RFS time, RFS event, and there is a treatment variable that is either SOC or FLIN.

I want to get the advantage of having the OS and RFS correlated by the fact that i’m reading the data one row at a time, and a row has both the OS and RFS. I don’t want to model correlation or covariance, but I want to get some benefit from the hierarchal structure.

Given that I have four sets of likelihoods I suspect that I may have lost some benefit from doing OS and RFS in the same model.

I’m not 100% clear on if what I’m talking about is a real issue. I have some small issues with divergencies, but the diagnostic plots and traces look decent and when I set up the plots using the estimated parameters they seem to track the KM curves pretty well.

I would like to replicate what the author is doing in this paper:

The author is working in R and stan and made a website, but it not clear that the website is up-to-date and there is not an easy to follow example with data.

4

Could you explain what those acronyms are?

5

SOC is standard of care for a type of cancer drug, often some type of chemotherapy. FLIN is a new drug treatment (Made up name). OS is overall survival and it is the time that you either die or drop out of the study for other reasons, or the study period ends. RFS is relapse free survival and is the time to relapse of cancer, death, of discontinuation in the study due to dropping out or the study period ending.

This is survival analysis and it looks at a population of patients. All the patients starts out both alive and relapse free. When you make a plot, you start at 1.0 or 100 percent. As you go through time, people either relapse, die, or are censored, and the curve moves down and eventually it hits zero, meaning either everyone has relapsed or died.

If you look at two curves, the higher one is better than the lower one. The lower one means that more bad things (death, relapse) are happening. A cancer drug shows that it’s better by plotting higher.

The results of the survival estimates are used in cost effectiveness models. In most cases the data spans three to six years. However, the time horizon on the cost effectiveness is often 20 or more years, even though with some cancers and other diseases patients are unlikely to live that long. The curves do extend out of sample, but the parametric estimates are the best guess of what will happen in the future.

The models generally compute discounted life years or utility from living, along with discounted costs. You compute these measures for both drugs, and find the ratio of the differences which is a statistic called the incremental cost-effectiveness ratio (ICER).

This type of analysis is used by agencies around the world to justify the purchasing of a new treatments. It is basically saying how much does it cost for each additional quality adjusted life year. If it cost less than $10,000 keep someone alive for an additional year the drug is probably a good deal, if it cost over $10,000,000 it is probably too expensive. In countries with Nationalized healthcare, they have to make a judgement on if they will pay for a new drug or not. This analysis helps them make their decision.

In my analysis when I extrapolate my survival plots my OS and RFS cross so RFS is higher than OS. However, by definition of the measure, this can’t happen, and is only happening because I’m extrapolating out of sample. I’m turning to to Bayesian analysis to keep this from happening. In my case I did it first by forcing the RFS and OS curves to have the same cure rate, and then I forced the OS to to have a slightly higher rate by adding a “d” which I forced to be a small positive number (0.02 to 0.06)