# PyMC modeling drug testing with truthiness

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.

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)

@pymc.deterministic
def transition_proportion_true(positive_test=positive_test,
transition_proportion=transition_proportion):
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.subplot(211)
plt.hist(true_use_fraction_samples, bins=25,
label="Posterior of drug use frequency.", normed=True)
plt.legend(loc="upper right")
plt.subplot(212)
plt.hist(transition_proportion_samples, bins=25,
label="Posterior probability to lie about use.", normed=True)
plt.legend(loc="upper right")


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.

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