I’m new to Bayesian modeling in general and PyMC in particular. I’m trying to model in a situation in which we ask patients whether they’ve taken drugs and we further drug test a subset of these samples. I’d like to combine this data taking into account the observation that patients can (often?) lie about whether they’ve taken drugs.

I’m having trouble modeling the truthiness of the patients. I don’t quite understand how I could make their admitted use be a function of both their actual use and their tendency to lie. I’ve got some example code below. Any suggestions?

import pymc
import numpy as np

admited_use = [True, True, False, False, True, False]
test_use  = [False, True, True, True, False, False]

true_use_fraction = pymc.Uniform('true_use_fraction', 0, 1)
lying_use_fraction = pymc.Uniform('lying_use_fraction', 0, 1)

positive_test = pymc.Bernoulli('positive_test', true_use_fraction, 
                               observed=True, value=np.array(test_use))

EDIT: I want to use the admit-use data for instances in which we don't have drug testing and we'd like to stratify patients based on how often they lie about their drug use history.


1 Answer 1


I've been learning this stuff myself lately, so I took a shot at setting up a model for your situation.

I broke the setup into two pieces. First, a person either does or does not use drugs. I set this part up much the same as you have it

test_use = np.array([1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 0, 1, 1, 0, 1]) == 1
true_use_fraction = pymc.Uniform('true_use_fraction', 0, 1)
positive_test = pymc.Bernoulli('positive_test', true_use_fraction, size=len(test_use),
                               observed=True, value=test_use)

So, our prior on the proportion of the population that are drug users is uniform, and given that, whether an individual is or is not a drug user is a bernoulli trial.

Next, I assumed that there is a population frequency of lying about drug use. Like this: if you are not a user, you do not lie and tell us you do use (because, who would lie about that); if you do use drugs, there is a latent probability that you will lie about it and say you don't.

The data we have is just whether the person said they were a user, so let's derive a variable measuring if you did or did not lie

admited_use = np.array([1, 1, 0, 1, 1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 1]) == 1
user_lied = test_use*(1 - admited_use) == 1

Now,create the transition proportion variable

transition_proportion = pymc.Uniform('transition_proportion', 0, 1)

def transition_proportion_true(positive_test=positive_test,
    ttrue = np.zeros(len(positive_test))
    ttrue[positive_test] = transition_proportion
    return ttrue

The derived variable, which represents the probability you will not tell the truth, is zero if you are a non-user, and otherwise is the latent probability of lying about your drug use.

Next, let's create a bernoulli variable for if the person lied, and link it to our data

transition = pymc.Bernoulli('transition', transition_proportion_true,
                            observed=True, value=user_lied)

Finally, we create the model

model = pymc.Model([true_use_fraction, transition_proportion,
                    positive_test, transition])

and run the sampler

mcmc = pymc.MCMC(model)
mcmc.sample(500000, 25000, 1)

Now we can sample from the posteriors and draw the resulting distributions

true_use_fraction_samples = mcmc.trace('true_use_fraction')[:]
transition_proportion_samples = mcmc.trace('transition_proportion')[:]

plt.figure(figsize=(12.5, 6))
plt.hist(true_use_fraction_samples, bins=25,
         label="Posterior of drug use frequency.", normed=True)
plt.legend(loc="upper right")
plt.hist(transition_proportion_samples, bins=25,
         label="Posterior probability to lie about use.", normed=True)
plt.legend(loc="upper right")

Posterior of drug use frequency and lying frequency

One Final Note: I got hung up on one point. I had tried to set up the model with the admitted drug use (survey response) as my second training data set. This created a situation where I had a pymc.deterministic observable. Apparently pymc does not like this, and I had a frustrating hour debugging until I realized the issue.

  • $\begingroup$ Perfect, thanks for the help. I also had trouble with deterministic observables and couldn't figure out how to work around it. $\endgroup$
    – JudoWill
    Commented Aug 9, 2016 at 14:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.