In [1]:
import arviz as az
import bambi as bmb
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
In [2]:
az.style.use('arviz-darkgrid')
Here we will analyze a dataset from experimental psychology in which a sample of 36 human participants engaged in what is called the shooter task, yielding 3600 responses and reaction times (100 from each subject). The link above gives some more information about the shooter task, but basically it is a sort of crude first-person-shooter video game in which the subject plays the role of a police officer. The subject views a variety of urban scenes, and in each round or "trial" a person or "target" appears on the screen after some random interval. This person is either Black or White (with 50% probability), and they are holding some object that is either a gun or some other object like a phone or wallet (with 50% probability). When a target appears, the subject has a very brief response window -- 0.85 seconds in this particular experiment -- within which to press one of two keyboard buttons indicating a "shoot" or "don't shoot" response. Subjects receive points for correct and timely responses in each trial; subjects' scores are penalized for incorrect reponses (i.e., shooting an unarmed person or failing to shoot an armed person) or if they don't respond within the 0.85 response window. The goal of the task, from the subject's perspective, is to maximize their score.
The typical findings in the shooter task are that (a) subjects are quicker to respond to armed targets than to unarmed targets, but are especially quick toward armed black targets and especially slow toward unarmed black targets, and (b) subjects are more likely to shoot black targets than white targets, whether they are armed or not.
In [3]:
shooter = pd.read_csv('data/shooter.csv', na_values='.')
shooter.head(10)
Out[3]:
The design of the experiment is such that the subject, target, and object (i.e., gun vs. no gun) factors are fully crossed: each subject views each target twice, once with a gun and once without a gun.
In [4]:
pd.crosstab(shooter['subject'], [shooter['target'], shooter['object']])
Out[4]:
The response speeds on each trial are recorded given as reaction times (milliseconds per response), but here we invert them to and multiply by 1000 so that we are analyzing response rates response rates (responses per second). There is no theoretical reason to prefer one of these metrics over the other, but it turns out that response rates tend to have nicer distributional properties than reaction times (i.e., deviate less strongly from the standard Gaussian assumptions), so response rates will be a little more convenient for us by allowing us to use some fairly standard distributional models.
In [5]:
shooter['rate'] = 1000.0/shooter['time']
In [6]:
plt.hist(shooter['rate'].dropna());
Our first model is analogous to how the data from the shooter task are usually analyzed: incorporating all subject-level sources of variability, but ignoring the sampling variability due to the sample of 50 targets. This is a Bayesian generalized linear mixed model (GLMM) with a Normal response and with intercepts and slopes that vary randomly across subjects. Of note here is the C(x, Sum) syntax, which is from the Patsy library that we use to parse formulae. This instructs Bambi to use contrast codes of -1 and +1 for the two levels of each of the fixed factors of race (black vs. white) and object (gun vs. no gun), so that the race and object coefficients can be interpreted as simple effects on average across the levels of the other factor (directly analogous, but not quite equivalent, to the main effects). This is the standard coding used in ANOVA.
In [7]:
subj_model = bmb.Model(shooter, dropna=True)
subj_fitted = subj_model.fit('rate ~ C(race, Sum)*C(object, Sum)',
random=['C(race, Sum)*C(object, Sum)|subject'],
samples=1000, init=None)
First let's visualize the default priors that Bambi automatically decided on for each of the parameters. We do this by calling the .plot() method of the Model object.
In [8]:
subj_model.plot();
The priors on the fixed effects seem quite reasonable. Recall that because of the -1 vs +1 contrast coding, the coefficients correspond to half the difference betwen the two levels of each factor. So the priors on the fixed effects essentially say that the black vs. white and gun vs. no gun (and their interaction) response rate differences are very unlikely to be as large as a full response per second. The priors on the standard deviations of the random effects are yoked to the corresponding fixed effect priors (see our technical paper for more info).
Now let's visualize the model estimates. We do this by calling the .plot() method of the MCMCResults object that resulted from our call to Model.fit().
In [9]:
subj_fitted.posterior
m_vars = [str(v) for v in subj_model.backend.model.vars]
plot_vars = []
# is not str(v).endswith(offset)]
In [10]:
az.plot_trace(subj_fitted, compact=True);
Each distribution in the plots above has 4 densities because we used 4 MCMC chains, so we are viewing the results of all 4 chains prior to their aggregation. The main message from the plot above is that the chains all seem to have converged well and the resulting posterior distributions all look quite reasonable. It's a bit easier to digest all this information in a concise, tabular form, which we can get by calling the .summary() method of the MCMCResults object.
In [11]:
# by default the variables plotted by subj_model.plot() are different than those from:
az.summary(subj_fitted)
Out[11]:
The take-home message from the analysis seems to be that we do find evidence for the usual finding that subjects are especially quick to respond (presumably with a shoot response) to armed black targets and especially slow to respond to unarmed black targets (while unarmed white targets receive "don't shoot" responses with less hesitation). We see this in the fact that the marginal posterior for the C(race, Sum)[S.black]:C(object, Sum)[S.gun] interaction coefficient is concentrated strongly away from 0.
A major flaw in the analysis above is that random stimulus effects are ignored. The model does include random effects for subjects, reflecting the fact that the subjects we observed are but a sample from the broader population of subjects we are interested in and that potentially could have appeared in our study. But the targets we observed -- the 50 photographs of white and black men that subjets responded to -- are also but a sample from the broader theoretical population of targets we are interested in talking about, and that we could have just as easily and justifiably used as the experimental stimuli in the study. Since the stimuli comprise a random sample, they are subject to sampling variability, and this sampling variability should be accounted in the analysis by including random stimulus effects. For some more information on this, see here, particularly pages 62-63.
To account for this, we let the intercept and object slope vary randomly across targets. Random object slopes across targets are possible because, if you recall, the design of the study was such that each target gets viewed twice by each subject, once with a gun and once without a gun. However, because each target is always either white or black, it's not possible to add random slopes for the race factor or the interaction.
In [12]:
stim_model = bmb.Model(shooter, dropna=True)
stim_fitted = stim_model.fit('rate ~ C(race, Sum)*C(object, Sum)',
random=['C(race, Sum)*C(object, Sum)|subject',
'C(object, Sum)|target'],
samples=1000)
Now let's look at the results...
In [13]:
az.plot_trace(stim_fitted, var_names=stim_model.term_names, compact=True);
In [14]:
az.summary(stim_fitted, var_names=stim_model.term_names)
Out[14]:
There are two interesting things to note here. The first is that the key interaction effect, C(race, Sum)[S.black]:C(object, Sum)[S.gun] is much less clear now. The marginal posterior is still mostly concentrated away from 0, but there's certainly a nontrivial part that overlaps with 0; 3.2% of the distribution, to be exact.
In [15]:
(stim_fitted.posterior['C(race, Sum):C(object, Sum)'] < 0).mean()
Out[15]:
The second interesting thing is that the two new variance components in the model, those associated with the random stimulus effects, are actually rather large. This actually largely explains the first fact above, since if these where estimated to be close to 0 anyway, the model estimates wouldn't be much different than they were in the subj_model. It makes sense that there is a strong tendency for different targets to elicit difference reaction times on average, which leads to a large estimate of 1|target_sd. Less obviously, the large estimate of C(object, Sum)[S.gun]|target_sd (targets tend to vary a lot in their response rate differences when they have a gun vs. some other object) also makes sense, because in this experiment, different targets were pictured with different non-gun objects. Some of these objects, such as a bright red can of Coca-Cola, are not easily confused with a gun, so subjects are able to quickly decide on the correct response. Other objects, such as a black cell phone, are possibly easier to confuse with a gun, so subjects take longer to decide on the correct response when confronted with this object. Since each target is yoked to a particular non-gun object, there is good reason to expect large target-to-target variability in the object effect, which is indeed what we see in the model estimates.
Here we seek evidence of the second traditional finding, that subjects are more likely to response 'shoot' toward black targets than toward white targets, regardless of whether they are armed or not. Currently the dataset just records whether the given response was correct or not, so first we transformed this into whether the response was 'shoot' or 'dontshoot'.
In [16]:
shooter['shoot_or_not'] = shooter['response'].astype(str)
# armed targets
new_vals = {'correct': 'shoot', 'incorrect':'dontshoot', 'timeout': np.nan}
shooter['shoot_or_not'][shooter['object']=='gun'] = \
shooter['response'][shooter['object']=='gun'].astype(str).replace(new_vals)
# unarmed targets
new_vals = {'correct': 'dontshoot', 'incorrect':'shoot', 'timeout': np.nan}
shooter['shoot_or_not'][shooter['object']=='nogun'] = \
shooter['response'][shooter['object']=='nogun'].astype(str).replace(new_vals)
# view result
shooter.head(20)
Out[16]:
Let's skip straight to the correct model that includes random stimulus effects. This looks quite similiar to the stim_model from above except that we change the response to the new shoot_or_not variable -- notice the [shoot] syntax indicating that we wish to model the prbability that shoot_or_not=='shoot', not shoot_or_not=='dontshoot' -- and then change to family='bernoulli' to indicate a mixed effects logistic regression.
In [17]:
stim_response_model = bmb.Model(shooter.dropna())
stim_response_fitted = stim_response_model.fit(
'shoot_or_not[shoot] ~ C(race, Sum)*C(object, Sum)',
random=['C(race, Sum)*C(object, Sum)|subject',
'C(object, Sum)|target'],
samples=1500, family='bernoulli')
Show the trace plot
In [18]:
az.plot_trace(stim_response_fitted, compact=True);
Looks pretty good! Now for the more concise summary.
In [19]:
az.summary(stim_response_fitted)
Out[19]:
There is some slight evidence here for the hypothesis that subjects are more likely to shoot the black targets, regardless of whether they are armed or not, but the evidence is a little weak. The marginal posterior for the C(race, Sum) coefficient is mostly concentrated away from 0, but it overlaps even more in this case with 0 than did the key interaction effect in the previous model.
In [20]:
(stim_response_fitted.posterior['C(race, Sum)'] < 0).mean()
Out[20]: