Explaining Recruitment to  Extremism:
A Bayesian Case-Control Approach

1]Roberto CerinaAuthor contributions: NK and CB conceived of the study. RC developed the models. RC, CB, and NK contributed to the analysis. RC and CB developed the code for the R package and the documentation. CB, NK and AZ contributed to the data collection. CB, NK and AZ developed the literature review. RC, CB, NK, and AZ contributed to the writing. We received helpful feedback and advice from Sir David Cox, Thomas Hegghammer, Bjørn Høyland, Bent Nielsen, Jacob Aasland Ravndal, and Frank Windmeijer. <cit.> reached out to us as we were finalizing our manuscript. Their analysis also uses leaked ISIS recruitment data to analyze the socio-economic correlates of joining the movement. To implement the methods described in this paper, see the associated R package, as well as documentation and vignettes: <http://extremeR.info>
2]Christopher Barrie
3]Neil Ketchley
4]Aaron Y. Zelin
[1]Maastricht University
[2]University of Edinburgh
[3]University of Oxford
[4]Brandeis University

Who joins extremist movements? Answering this question poses considerable methodological challenges. Survey techniques are infeasible and selective samples provide no counterfactual. Recruits can be assigned to contextual units, but this is vulnerable to problems of ecological inference. In this article, we take inspiration from epidemiology and elaborate a technique that combines survey and ecological approaches. The multilevel Bayesian case-control design that we propose allows us to identify individual-level and contextual factors patterning the incidence of recruitment, while accounting for rare events, contamination, and spatial autocorrelation. We validate our approach by matching a sample of Islamic State (ISIS) fighters from nine MENA countries with representative population surveys enumerated shortly before recruits joined the movement. High status individuals in their early twenties with university education were more likely to join ISIS. There is more mixed evidence for relative deprivation. We provide software for applied researchers to implement our method.

§ INTRODUCTION

Identifying who is more or less likely to join extremist movements is a pressing issue for both political science and public policy.[We understand extremists to be movements or individuals that reject or seek to subvert democracy, often through violent means <cit.>.] However, empirical research on this topic is beset by methodological challenges. Population surveys offer little insight into the phenomenon as recruits to extremism are tiny minorities in any society, and so are tiny minorities in survey samples. This is before obvious problems related to eliciting truthful responses to questions probing illicit actions. Recent innovations in survey and online digital trace methodologies have allowed researchers to obtain more accurate measures of support for extremism <cit.>. However, these approaches capture attitudes, rather than behaviour. For researchers interested in why some individuals join extremist movements and not others, the most common strategy is to collect a convenience sample of recruits after which one of two approaches are commonly preferred. A first sees researchers simply report sample proportions of a given characteristic, e.g. the percentage of recruits who have university education. While often descriptively useful, this fails to account for population baselines and other confounding factors affecting the incidence of recruitment. An alternative strategy is to assign recruits to meaningful contexts and use the characteristics of those places to explain variation in the recruitment rate. While this approach provides a counterfactual and allows for multivariate analysis, it suffers from familiar problems of ecological inference <cit.>.

The method we propose in this paper allows researchers to leverage both survey and contextual data to make robust inferences about the individual and ecological correlates of recruitment to extremism. To do so, we take inspiration from the case-control design used in epidemiology and show how it can be adapted to combine a convenience sample of cases (extremists) with controls (respondents from a representative survey). In this, we build on the recent introduction of case-control methods to political science found in the study of protest participation <cit.>. Several statistical challenges arising from the nature of extremism remain, however. Popular approaches for modelling rare events [e.g. ] do not account for spatial autocorrelation in the incidence of recruitment or the possibility of contamination between cases and controls. Further, conventional rare events approaches face potential separation problems in the absence of regularized priors for coefficients of interest. Finally, while existing approaches incorporate informative priors on the true prevalence of the outcome of interest (here: recruitment to extremism), this information is often not known.

Our approach offers a complete solution to these statistical problems and can be described as a multilevel, Bayesian case-control design that is robust to rare events, contamination, and spatial autocorrelation patterning the incidence of recruitment <cit.>. The Bayesian approach is preferable for a number of reasons. First, Bayesian probabilistic programming software provides unique flexibility in the modeling of the complex hierarchical structures characterizing recruitment into extremism. Second, it permits the use of informative priors to account for the true prevalence of the event of recruitment, as well as to regularize coefficient estimates to account for separation bias and instability when carrying out regressions <cit.>. Finally, in the absence of prior knowledge of the overall propensity of being a recruit in a given context, the model can estimate the propensity from the data <cit.>.

A great strength of our method is that applied extremism researchers can choose those parameters most relevant for their case. When sampling from national populations, the risk of contamination between cases and controls may be sufficiently low such that it does not pose a threat to inference. On the other hand, recruitment may not qualify as a rare event when comparing recruits to certain sub-populations. So too, spatial autocorrelation in recruitment may not apply if sampling from a small area or closed context. As we show, our modelling strategy if flexible to the inclusion or exclusion of these parameters, depending on the case at hand. In support of our approach, and to help guide the modelling decisions of future practitioners, we provide an extensive simulation study that compares our model to alternative frameworks, and show its robustness and superiority in predicting the true underlying probability of recruitment under various bias-inducing scenarios.

To display some of the key properties of our modeling strategy, we analyze recruitment of Sunni Muslim males in nine MENA countries to the Islamic State in Iraq and Syria (ISIS). We focus our analysis on an individual's level of education and social status — two key variables associated with recruitment to extremism found in the literature on violent Islamist movements. We show how our approach can be used to perform two types of analyses. In the first, we leverage a multilevel regression model trained on a cross-national sample of ISIS recruits and non-recruits. This provides a robust descriptive analysis about the individual-level characteristics of recruits across countries and sub-national administrative units. A second analysis focuses on two countries for which we have rich contextual information: Egypt and Tunisia. This analysis adds value by adjusting for local heterogeneity with the addition of relevant ecological covariates, allowing us to ascertain the potential sensitivity of individual-level findings to unobserved contextual confounding.

For the purposes of illustration, we implement the complete solution described above, accounting for the possibility of contamination, spatial autocorrelation, and separation in our regression coefficients. Overall, we find that high-status males with university education in their early twenties were more likely to join ISIS. We also find that relatively deprived males in Egypt were more likely to join ISIS, but not in Tunisia. This heterogeneity in the individual and contextual correlates of violent extremism demonstrates the importance of accounting for both individual and context specific factors.

While our analysis focuses on recruitment to ISIS, our hope is that this paper inspires political scientists to apply case-control methods to other instances of extremism where data on recruits and population surveys are available. Examples include participation in the 2021 attack on the Capitol Building in Washington D.C. <cit.>, recruitment to far-right movements and white supremacist groups <cit.>, as well as other examples of violent extremism <cit.>. In this spirit, the dedicated R package accompanying this paper allows applied researchers working on such cases to implement our proposed technique.[See: <http://extremeR.info>.]

§ EXPLAINING RECRUITMENT TO EXTREMISM

The principle strategy available to researchers interested in the correlates of recruitment to extremism is to sample on the dependent variable, obtaining relevant demographic information on individual extremists or members of extremist movements. In the ideal scenario, researchers are able to obtain movement membership lists, which can reveal information on tens of thousands of individuals <cit.>, although in practice such complete data is rare. Absent such lists, a well-established strategy is to leverage data from arrests or killings to generate samples of participants <cit.>. Alternatively, researchers can look to collect demographic information on extremists by either interviewing former recruits <cit.>, or by reconstructing the biographical profiles of prominent extremists from open source information <cit.>. Of course, a principle limitation of these samples is that they do not provide information on individuals outside of the subpopulation of interest, meaning that it is not possible to compare recruits to the population from which they are drawn. To remedy this, researchers either confine attention to variation amongst recruits <cit.>, or else assign individuals to meaningful contexts, e.g. universities, cities, or countries, and then use the characteristics of those units to explain cross-sectional variation in the recruitment rate <cit.>. While this latter approach is undoubtedly superior to simply analyzing sample proportions, it inevitably relies on ecological inference.[An alternative approach adopted in the extremism literature is the “life history" method, whereby researchers draw on small convenience samples of extremists in an effort to identify events and individual-level characteristics that made them more likely to participate in extremism <cit.>. As with simply reporting sample proportions, this method lacks an obvious counterfactual, meaning that we cannot rule out the possibility that factors identified as leading individuals to pursue extremist behavior are also present in individuals who did not become extremists <cit.>.]

§.§ A multilevel Bayesian case-control design

In what follows, we suggest two new methods for analyzing recruitment to extremism. The first leverages a cross-national, multilevel regression model trained on a complete sample of recruits and survey respondents. This  provides a robust descriptive analysis about the individual-level factors which characterize recruits across countries and subnational units. The model uses random effects to control for unobservable subnational heterogeneity; these are preferable to fixed effects due to potentially heavily imbalanced area-level sample sizes <cit.>. The model further uses a conditionally auto-regressive prior <cit.> to account for spatial smoothing. The second analysis focuses on single country studies where rich contextual information is available. The added value of this analysis lies in controlling for local heterogeneity in order to ascertain the robustness of any individual-level findings to contextual confounding. Taken together, our proposed setup thus plots a way forward for researchers to combine survey and ecological information for the robust analysis of recruitment to extremism.

§.§ Simple case-control set-up

We begin by describing the backbone of our model, which is a logistic regression accounting for case-control sampling protocol via an offset. Borrowing from Rota et al. , we define $r_i=\{0,1\}$ as the set of states that observation $i$ in our sample of size $n = n_0 + n_1$ can obtain, where $r_i = 1$ implies the observation is a `case', $r_i = 0$ defines a control, $n_1 = \sum_i^n \mathds{1}(r_i = 1)$ and $n_0 = \sum_i^n \mathds{1}(r_i = 0)$. In our application, a `case' would refer to a known extremist; a `control' to a survey respondent. Recall that cases are selected entirely on the dependent variable while controls come from the population that cases are drawn from. Take $N_1$ to represent the number of cases in the population of interest, and $N_0$ the number of controls. The probability of being included in the sample ($s_i = 1$) conditional on the true state of any individual can hence be understood as $P_1 = \mbox{Pr}(s_i = 1 \mid r_i = 1) = \frac{n_1}{N_1}$, while that of being sampled as a control is $P_0 = \mbox{Pr}(s_i = 1 \mid r_i = 0) = \frac{n_0}{N_0}$. The log-ratio of these sampling probabilities can then be used as an `offset' in a logistic regression, to account for the sampling protocol. The hierarchical specification of the model follows, with regression coefficients being assigned a very weakly informative prior;[See <https://github.com/stan-dev/stan/wiki/Prior-Choice-Recommendations> for prior-choice advice when using . Note further that the normal distribution in our model (and in ) is parametrized by mean and standard deviation.]

\begin{align}
r_i & \sim \mbox{Bernoulli}(\rho_i);\\
\mbox{logit}(\rho_i) & = \mbox{log}\left(\frac{P_1}{P_0}\right) + \sum_k x_{i,k} \beta_{k};\\
\beta_k & \sim N(0,10).
\end{align}

The above hierarchical model thus contains three layers: layer (1) is a model of the true state of an observation, conditional on their latent propensity $\rho$; layer (2) describes this latent propensity, by accounting for systematic variation due to heterogeneity in covariates; layer (3) models the effects of each covariate by assigning a prior probabilistic model.

§.§ Contaminated controls

Recall that that the case-control setup as described above takes known recruits and combines them with `controls' taken from survey respondents. While we know that our cases are correctly labeled, we do not know whether this is true of our controls. That is, our controls may be `contaminated' as survey respondents may have become recruits <cit.>. This is especially concerning when researchers have access to biographical information on tens of thousands of extremists  <cit.> or are comparing recruits to small sub-populations <cit.>. <cit.> outline a `latent variable' formulation of their contamination model. Below we present our version of that same model as a mixture, which we find more intuitive.

The `label' of an observation, $y_i=\{0,1\}$, is observed for all observations, while the true `state' of an observation, $r_i=\{0,1\}$, is only observed for cases. The implied probability distribution of labels conditional on being a control is:

\begin{align*}
\mbox{Pr}(y_i = 1 &\mid r_i =0 , s_i = 1 ) = 0 = \theta_{0};\\
\mbox{Pr}(y_i = 0 &\mid r_i =0 , s_i = 1 ) = 1 = (1-\theta_{0} );
\end{align*}

Due to contamination, it is possible that observations characterized by $y_i=0$ are actually in state $r_i = 1$; hence we need a probability distribution for $y \mid r_i = 1$. Let $\pi = \frac{N_1}{N_1 + N_0}$ be the prevalence of recruits in the population of interest, and let $n_u = \sum_i^n \mathds{1}(y_i = 0)$ be the number of unlabeled observations. We expect there to be $\pi n_u$ cases amongst the unlabeled observations. We can then characterize the probability distribution of labels, conditional on being a case, as:

\begin{align*}
\mbox{Pr}(y = 1 \mid r =1 , s = 1) = \frac{n_1}{n_1 + \pi n_u} = \theta_{1};\\
\mbox{Pr}(y = 0 \mid r =1 , s = 1) = \frac{\pi n_u}{n_1 + \pi n_u} = (1-\theta_{1}).
\end{align*}

Finally, our model for the latent state $r_i$ must reflect the possibility of contamination. We do this by re-defining the relative-risk of being sampled as:

\begin{align*}
\frac{P_1}{P_0} = \frac{\frac{n_1 + \pi n_u}{N_1}}{\frac{ (1-\pi)n_u}{N_0}} = \frac{n_1}{\pi n_u} + 1.
\end{align*}

The updated, hierarchical specification for the case-control model accounting for contaminated controls is then:

\begin{align}
y_i & \sim \mbox{Bernoulli}(\theta_{r_i});\\
r_i & \sim \mbox{Bernoulli}(\rho_i);\\
\mbox{logit}(\rho_i) & = \mbox{log}\left(\frac{n_1}{\pi n_u} + 1\right) + \sum_k x_{i,k} \beta_{k};\\
\beta_k & \sim N(0,10).
\end{align}

In summary, we derive our labels via two distinct data generating processes, identified by a latent state $r_i = \{1,0\}$. In the event that the latent state of a given record is that of a true control, $r_i = 0$, it is then impossible for this record to be labeled $y_i = 1$; conversely, if the latent state is that of a true case, $r_i = 1$, then it is still possible for a record to be labeled $y_i = 0$, with probability $(1-\theta_1)$. This latter model describes the issue of contamination. Note that in our application, $\theta$ is always observed, and fed to the model as data.

§.§ Area-level random effects

Survey data and information on recruits often contain information on the origin or location of residence of individuals. And we can understand individuals as nested within geographical units of increasing sizes. Generalizing, we can exploit variance at three levels: the individual, some small-area, and some large-area.

These area effects could be incorporated in the model via fixed-effects, by expanding the design matrix to include relevant dummy-variables for each area of interest. We consider this strategy unwise when trying to explain recruitment to extremism and prefer a random-effects approach. In the case of rare forms of political behaviour, our geographical units at all levels of analysis will have relatively few observations <cit.>. Additionally, for many units, we will have no cases. Finally, we know that lists of recruits are unlikely to be exhaustive; that is, we will not have data for every recruit hailing from every subnational unit or country. Here, a sample of recruitment data or similar can be treated as a non-probability sample — it is unlikely that we can have complete confidence the sample constitutes a complete or random sample of the population of interest. Given these concerns, a random effects approach is preferable as it means: 1) we are able to borrow strength across areas, which also increases efficiency, to produce more realistic estimates for the area-level coefficients <cit.>; 2) in the absence of more detailed knowledge about the data-generating process, the shrinkage effect obtained by partial pooling is more likely to shield our estimates from any systematic sampling bias among our cases <cit.>.

We can also relax some of the theoretical bias associated with the shrinkage induced by random effects via incorporating observable area-level heterogeneity in the design-matrix as fixed effects <cit.>. This is what we elect to do in single country analyses. Finally, it is worth highlighting that our goal is not to make inferences about area-level effects. Rather, we seek to strip our individual-level effects estimates of contested variance that may be associated with the provenance of the recruit. The resulting hierarchical model is as follows:

\begin{align}
y_i & \sim \mbox{Bernoulli}(\theta_{r_i});\\
r_i & \sim \mbox{Bernoulli}(\rho_i);\\
\mbox{logit}(\rho_i) & = \mbox{log}\left(\frac{n_1}{\pi n_u} + 1\right) + \sum_k x_{i,k} \beta_{k} + \phi_{l[i]} + \eta_{j[i]} ;\\
\beta_k & \sim N(0,10);\\
\phi_l & \sim N(0,\sigma_\phi);\\
\sigma_\phi = \frac{1}{\sqrt{\tau}_\phi}, \mbox{ } \tau_\phi &\sim \mbox{Gamma}(\epsilon,\epsilon);\\
\eta_j & \sim N(0,\sigma_\eta);\\
\sigma_\eta = \frac{1}{\sqrt{\tau}_\eta}, \mbox{ } \tau_\eta &\sim \mbox{Gamma}(\epsilon,\epsilon);
\end{align}

where $\epsilon$ stands for some arbitrary number, chosen as a compromise to minimize the prior information and maximise the Markov chain Monte Carlo (MCMC) convergence speed and stability.

§.§ Spatial autocorrelation

The network ties connecting actors across space play an important role in recruitment to high-risk activism and extremism <cit.>. Sometimes the ties connecting recruits will be available; more commonly this information will not be recoverable.

In the absence of detailed network information, we propose controlling for potential network effects at levels of varying scale. We work on the assumption that network ties are more likely to form between individuals who are geographically proximate. Depending on the richness of the data on recruits, we may generate distance matrices between geographical units of varying size.

To account for area-level spatial autocorrelation, we propose incorporating a version of the conditional auto-regressive (CAR) model <cit.>. This approach has been used in individual-level models of behaviour, enabling local smoothing of predictions according to behaviour observed in neighbouring areas <cit.>. The key ingredients of a CAR model are $\bm{\omega}$, a distance-weight matrix; $\alpha$, a parameter governing the degree of autocorrelation, where $\alpha =0$ implies spatial independence, and $\alpha =1$ implies an intrinsic conditional auto-regressive (ICAR) model <cit.>; and $\sigma_{\psi}$, the standard deviation of the subnational unit effects. The resulting model for spatial random effect $\psi_l \mbox{ } \forall \mbox{ }l = \{1,...,L\}$   is then:

\begin{equation*}
\psi_l \mid \psi_{l^\prime}   \sim N\left(\alpha\sum_{l^{\prime} \neq l} \omega_{ll^{\prime}} \psi_{l^\prime},\sigma_{\psi} \right).\\
\end{equation*}

In practice, we implement the ICAR specification of the model, with $\alpha = 1$, and take $\bm{\omega}$ to be the neighbourhood matrix. The neighbourhood matrix has diagonals zero (a unit cannot neighbour itself) and off-diagonal zero or one depending on whether the given units are neighbours. We choose this specification of the distance matrix because of the efficiency gains it affords in a Bayesian context  <cit.>. This leads to:

\begin{equation*}
\psi_l \mid \psi_{l^\prime}  \sim N\left(\frac{\sum_{l^{\prime} \neq l} \psi_{l^\prime}}{d_{l,l}},\frac{\sigma_{\psi}}{\sqrt{d_{l,l}}} \right),
\end{equation*}

where $d_{l,l}$ is an entry of the diagonal matrix $D$ of size $L \times L$, whose diagonal is defined as a vector of the number of neighbours of each area. The joint distribution of this model is simply a multivariate normal distribution $\bm{\phi} \sim N(0,[\tau_\psi(D - W)]^{-1})$, $\tau_\psi = \frac{1}{\sigma^2_\psi}$, which is conveniently proportional to the squared pairwise difference of neighbouring effects. Note that the sum-to-zero constraint is needed for identifiability purposes, as in its absence any constant added to the $\psi$s would cancel out in the difference.[Generally, this model has the disadvantage of being an improper-prior, as its density does not integrate to unity and is non-generative (i.e., it can only be used as a prior and not as a likelihood) though it serves our purposes within the context of a hierarchical model. Finally, we should note that this prior encodes an intrinsic dependence between subnational units. It can no longer detect the degree of spatial autocorrelation supported by the data but instead assumes that areas are explicitly dependent, and estimates coefficients accordingly.] Following <cit.>, setting the precision to $1$ and centering the model such that $\sum^L_l \psi_l = 0$, we arrive at:

\begin{equation*}
\mbox{log } p(\bm{\psi} ) \propto \mbox{exp} \left\{ -\frac{1}{2} \sum_{l^{\prime} \neq l}(\psi_l - \psi_{l^\prime})^2 \right\};
\end{equation*}

The hierarchical model we implement to incorporate the spatial component is within the Besag-York-Mollié (BYM) family <cit.>. For a given level of analysis, say the city or province in a cross-country analysis, BYM models are characterized by two random effects which explain unobserved heterogeneity: $\phi_l$ defines a non-spatial component while $\psi_l$ defines systematic variance due to spatial dependency. The typical challenge with BYM is that the two areal effects cannot be identified without imposing some structure since they are mutually dependent, meaning either component is capable of accounting for contested variance at the area-level. This leads to inefficient posterior exploration of any MCMC sample, and subsequent lack of convergence <cit.>. To overcome this, we implement a state-of-the-art solution leveraging penalized-complexity priors <cit.>, which proposes modelling the two effects as a scaled mixture such that:

\begin{align*}
\gamma_l = \sigma \left( \phi_l\sqrt{(1-\lambda)} + \psi_l\sqrt{(\lambda/s)}  \right);
\end{align*}

where $\phi$ and $\psi$ are random effects scaled to have unitary variance and $\lambda \in [0,1]$ is a mixing parameter, defining the proportion of residual variation attributable to spatial dependency. In order for the spatial and unstructured effects to share $\sigma$, they must be on the same scale <cit.>. We must therefore scale the ICAR-distributed effects, as their original scale is defined by the local neighbourhood. A proposed scaling factor is chosen such that the geometric mean of the variance parameters over the areal units is $1$, $\mbox{Var}(\psi_l)= 1$. Note that this scaling factor, $s$ in the equation above, can be calculated directly from the adjacency matrix, and hence it is not to be estimated but passed to the model as data.

The resulting hierarchical specification of our model follows:

\begin{align}
y_i & \sim \mbox{Bernoulli}(\theta_{r_i});\\
r_i & \sim \mbox{Bernoulli}(\rho_i);\\
\mbox{logit}(\rho_i) & = \mbox{log}\left(\frac{n_1}{\pi n_u} + 1\right) + \sum_k x_{i,k} \beta_{k} + \gamma_{l[i]} + \eta_{j[i]} ;\\
\beta_k & \sim N(0,10);\\
\gamma_l & = \sigma \left( \phi_l\sqrt{(1-\lambda)} + \psi_l\sqrt{(\lambda/s)}  \right);\\
\lambda & \sim \mbox{Beta}(0.5,0.5);\\
\phi_l & \sim N(0,1);\\
\psi_l \mid \psi_{l^\prime}  & \sim N\left(\frac{\sum_{l^{\prime} \neq l} \psi_{l^\prime}}{d_{l,l}},\frac{1}{\sqrt{d_{l,l}}} \right)\\
\sigma & \sim  \frac{1}{2}N(0,1);\\
\eta_j & \sim N(0,\sigma_\eta);\\
\sigma_\eta = \frac{1}{\sqrt{\tau}_\eta}, \mbox{ } \tau_\eta &\sim \mbox{Gamma}(\epsilon,\epsilon);
\end{align}

where $\frac{1}{2}N$ denotes a half-normal distribution, which is the recommended prior for the variance of BYM effects <cit.>.

§.§ Regularizing prior coefficients

Multiple contributions have highlighted problems with logistic regression coefficient estimates under rare-events <cit.>. The intuition behind these challenges is typically described as some variation on the standard separation problem where any given covariate or simple combination thereof perfectly separates cases from controls. This leads to biased and unstable point-estimates with large associated uncertainty <cit.>. A number of regularization techniques have been proposed to reduce bias and stabilize the coefficient estimates. Our preferred regularization method is that proposed by <cit.> and <cit.>. The approach assumes it should be unlikely to observe unit-changes in the (standardized) covariates that would lead to outcome changes as large as $5$ points on the logit scale. Using a slight variation on this approach to ensure sufficient regularisation, we use a Cauchy prior with scale-parameter set to $1$ for the regression coefficients, and a `looser' scale of $10$ logit points on the intercept to accomodate for the rarity of the event in the sample. The advantages of the Cauchy prior lie in its fat tails, which avoid over-shrinkage of large coefficients <cit.>. We apply this prior to our fixed effects exclusively, as the likelihood of our random effects is already structured and penalized. Our final model specification is then as follows:

\begin{align}
y_i & \sim \mbox{Bernoulli}(\theta_{r_i});\\
r_i & \sim \mbox{Bernoulli}(\rho_i);\\
\mbox{logit}(\rho_i) & = \mbox{log}\left(\frac{n_1}{\pi n_u} + 1\right) + \sum_k x_{i,k} \beta_{k} + \gamma_{l[i]} + \eta_{j[i]} ;\\
\beta_1 & \sim \mbox{Cauchy}(0,10);\\
\beta_k \mid k>1 & \sim \mbox{Cauchy}(0,1);\\
\gamma_l & = \sigma \left( \phi_l\sqrt{(1-\lambda)} + \psi_l\sqrt{(\lambda/s)}  \right);\\
\lambda & \sim \mbox{Beta}(0.5,0.5);\\
\phi_l & \sim N(0,1);\\
\psi_l \mid \psi_{l^\prime}  & \sim N\left(\frac{\sum_{l^{\prime} \neq l} \psi_{l^\prime}}{d_{l,l}},\frac{1}{\sqrt{d_{l,l}}} \right)\\    \sigma  &\sim \frac{1}{2}N(0,1);\\
\eta_j & \sim N(0,\sigma_\eta);\\
\sigma_\eta = \frac{1}{\sqrt{\tau}_\eta}, \mbox{ } \tau_\eta &\sim \mbox{Gamma}(\epsilon,\epsilon);
\end{align}

§.§ Bayesian modeling in

The model outlined above is amenable to Bayesian estimation via Monte Carlo Markov Chain (MCMC) methods. Previous contributions to the case-control literature [e.g., ] have used  <cit.> or  <cit.> as software to implement some variations on a simple Gibbs sampler. Due to the heavy computational burden imposed by the spatial prior, we propose instead to innovate by estimating this model in  <cit.>.  leverages a version of Hamiltonian Monte-Carlo (HMC) called the `No U-Turn Sampler' (NUTS) <cit.>, which dramatically improves the efficiency and speed of convergence of our Markov-Chains. A challenge we face is that  cannot handle the sampling of latent discrete parameters ($r_i$ in our hierarchical model above), posing a problem for the estimation of mixture models. The state-of-the-art solution is to marginalize the latent parameter out. In practice this means replacing our model for the observed labels $y$ with the following mixture of Bernoulli distributions:

\begin{equation}
f ( y_i \mid  \rho_i ) =  \rho_i  \mbox{ Bernoulli} ( y_i \mid \theta_{1 } )  + ( 1 - \rho_i ) \mbox{ Bernoulli} (y_i \mid \theta_{0} ).
\end{equation}

Beyond allowing for model parameters to be informed by $y_i$ according to the mixed structure above, marginalization provides significant advantages for posterior exploration and MCMC efficiency as it leverages expectations rather than sampling of discrete parameters. Listing <ref> in the Appendix presents the  code for our final model. Note that fixed-effects covariates are standardized.[We standardize both dichotomous and continuous variables as this aids convergence.] Estimates of regression coefficients on the original, unstandardised scale are computed and available in the Appendix.

§.§ Simulation study

The model that we propose extremism researchers should adopt is significantly more complex than the standard case-control design using rare-events logistic regression and requires a substantial understanding of Bayesian methods to be fully appreciated. Moreover, the model's estimation becomes roughly exponentially more computationally challenging as the sample size increases. To provide evidence that our approach is nevertheless preferable to a more straightforward case-control design, in Appendix Section <ref> we report the results of a comprehensive simulation study that compares the performance of our model against the King and Zeng model , as well as a simple fixed-effects logistic regression. We score these models according to their ability to accurately predict the underlying latent propensity of recruitment, $\mu_i = \mbox{logit}({\rho_i})$. We further investigate these models' performance in accurately estimating the intercept, regression coefficients, and residual area effects.

Simulations show that our model is robust and general. The results suggest that in a rare-event scenario, our model outperforms King and Zeng's rare-events logistic regression thanks to its ability to account for spatial auto-correlation, whilst also remaining largely unbiased to discrepancies in sample and population prevalence. As prevalence increases, our model retains a degree of robustness that neither a simple fixed-effects logistic regression, nor the `rare events logit', can offer – largely thanks to the contamination layer. This robustness extends not just to the ability to correctly estimate latent propensity $\bm{\mu}^\star$, but actively reduces bias and RMSE in the estimation of coefficients.

§ CORRELATES OF ISLAMIST EXTREMISM

Determining the demographic profile of recruits to violent extremism is a central concern in studies of terrorism and political violence. Among the most commonly studied characteristics is education <cit.>. In their review of the literature, <cit.> summarize three mechanisms linking education to reduced political violence. First, the availability of educational infrastructure might lower societal grievances and lead to a consequent decline in support for extremist groups. Second, the opportunity costs of engaging in violence are greater for more highly educated individuals. Third, education helps spread norms of tolerance and civic engagement, lowering support for forms of political extremism.

An alternative body of literature argues that greater access to higher education may increase recruitment to violent extremism, in part as a result of relative deprivation <cit.>. This argument has been popular among scholars studying political violence in the MENA where it is argued that mass educational expansion created a pool of highly educated unemployed individuals who were attracted to high risk activism <cit.>. A parallel literature focuses on the educational profiles of Islamist and jihadist groups. The key finding of this scholarship is that recruits are drawn disproportionately from the highly educated <cit.>. Again, relative deprivation arising from unrealized social status is often cited here as an explanatory factor with MENA countries having expanded higher education out of step with labor market demand, creating a “lumpen intelligentsia" vulnerable to recruitment <cit.>. An alternative interpretation is that the highly educated often have greater feelings of efficacy and a sense that they can understand and change politics <cit.>.

An adjacent body of research focuses on recruitment to transnational terrorist groups such as ISIS and Al-Qaeda. These contributions, which tend to be ecological and operate at the country level, find associations between human development indices and horizational inequalities in educational attainment as key factors patterning the flow of recruits to groups <cit.>. Using a sample of leaked border documents similar to those we exploit below, <cit.> find that recruits tend to be highly educated and have high social status. Other research stresses the local-level determinants of ISIS recruitment, comparing the demographic characteristics of fighters with local-level averages for the general population. While not multivariate, these studies conclude that fighters tend to hail from marginalized regions and areas with prehistories of contention, but that there remains significant across-context variation in the correlates of recruitment <cit.>.

In summary, there are somewhat competing predictions for why some individuals may be more likely to be recruited to violent extremism. University education, relative deprivation, and high social status are commonly advanced to explain an individual's propensity to join violent Islamist groups. Scholarship also points to important subnational, ecological factors in the incidence of recruitment. Taken together, the literature would benefit from a method that combines information at the individual and ecological levels to appraise those factors that make recruitment more or less likely.

§ DATA

Our key data source to study ISIS is a set of leaked border documents. This leak was widely covered in international news media and has been used to provide descriptive statistics on the geographical distribution and demographic characteristics of ISIS fighters from multiple MENA countries <cit.>.[Dodwell et all  suggest that around 98% of these individuals can be matched against records for ISIS recruits held by the U.S. Department of Defense.] When an individual entered Islamic State territory from Turkey they passed though an ISIS-controlled border control point. The leaked documents derive from these border control offices. They contain detailed information on the home residence of each recruit, age, education, marital status, previous employment, employment status, previous combat experience, and date of entry into ISIS-controlled territory. Figure <ref> provides an imitation of one of the border documents. We use data for nine countries in the MENA that were included in the leak. These are: Algeria, Egypt, Jordan, Kuwait, Lebanon, Libya, Morocco, Tunisia, and Yemen. In total, we have complete records for 1,051 recruits.

Example of ISIS border document (details changed)

For the case-control design, we combine individual-level ISIS recruitment data with a nationally representative sample of individuals from Wave III of the Arab Barometer  survey. The fieldwork for these surveys was completed before most recruits recorded in our border documents entered ISIS-held territory, and so may be vulnerable to contamination.[The surveys were in the field at different times for each country: December, 2012-January, 2013 for Jordan; February, 2013 for Tunisia; March-April, 2013 for Egypt and Algeria; April-June, 2013 in Morocco; July 2013 in Lebanon; November-December 2013 in Yemen; and February-March 2014 in Kuwait <cit.>.] Our choice of covariates to use from this survey is naturally constrained by the information included in the border documents. We elect to include covariates for age, marital status, university education, and student status. We combine two variables for unemployed and employment in agricultural or manual labor to create a composite variable designed to measure “low status" activity. An interaction between this variable and our university education variable is designed to capture the relative deprivation hypothesis; that is, whether highly educated individuals engaged in low status economic activity are more likely to become recruits.

As per our models above, we present a “bird's eye" and “worm's eye" view of ISIS recruits. The first approach uses just individual-level information for recruits across our nine countries. Here, we leverage a multilevel regression model trained on the complete sample of $1,051$ recruits and $5,093$ unlabeled records.[We lose $248$ cases and $225$ unlabeled records due to incompleteness. We considered incorporating a multiple imputation model in our Bayesian analysis, but decided against it: a) in the interest of clarity, as a number of layers of complexity are already present; b) in the interest of estimation stability, as incorporating further uncertainty could destabilize the estimates of coefficients; c) we have no a-priori reason to believe missingness is correlated with the outcome, and hence we do not expect that listwise deletion will bias the estimates.] The second incorporates contextual information for Egypt ($n_1 = 66, n_0 = 551$ complete records[We lose $15$ controls and $95$ cases to listwise deletion of incomplete records in the Egypt dataset.]) and Tunisia ($n_1 = 426$, $n_0 = 589$ complete records[We lose $5$ controls and $20$ cases to listwise deletion of incomplete records in the Tunisia dataset.]) at the district level. We focus on these two countries due to the availability of contextual information at the district level that is not easily accessible for the other countries in our sample. Again, the added value of this analysis lies in controlling for observable district-level heterogeneity in order to ascertain the robustness of any individual-level findings to contextual confounding. As such, for each country we include variables to capture subnational differences in demographic and labor-market composition, employment opportunities, as well as more context-specific variables designed to capture support for Islamist political organizations and prehistories of contentious politics. Our choice of contextual variables is based on existing research finding that lack of employment opportunities, prehistories of mobilization and repression, as well as support for political Islam, are predictive of ISIS recruitment <cit.>. As such, the inclusion of these variables provides insight into the potential magnitude and origin of contextual confounding when analyzing recruitment to extremism. We describe the sources and operationalization of these ecological variables in Appendix <ref>.

§ RESULTS

In this section we present the results of our two analytic approaches: the `bird's eye' view, referring to our cross-country analysis, and the `worm's eye' view, referring to our two country-specific analyses. For each, we present: i) convergence diagnostics; ii) the posterior density of fixed and random effects according to our models; 3) the posterior predictive distribution across potential recruitment profiles.

§.§ MCMC Convergence

Convergence diagnostics provide a first measure of the reliability of our parameter estimates for both the Bird's Eye and Worm's Eye models. Here, we follow Vehtari et al  and implement multiple state-of-the-art tests.

Per Appendix Figures <ref>, <ref>, and <ref>, we examine four versions of the Gelman-Rubin statistic ($\hat{R}$) to verify convergence is obtained broadly, as well as when we encounter heteroskedasticity across chains, or when these are heavy-tailed. There exist various convergence-thresholds in the literature – the most stringent requires $\hat{R}<1.01$, a medium-stringency threshold suggests $\hat{R}<1.05$ (especially if we are estimating a large number of parameters), whilst the historical recommendation was $\hat{R}<1.1$ <cit.>. Recent work demonstrates that this latter threshold prematurely diagnoses convergence in most cases <cit.>. The parameters of all of our models are broadly convergent under the harshest $1.01$ threshold for all of the measures of $\hat{R}$, with the exception of a very small number of spatial effects which are convergent under a slightly more laxed threshold, though still well below the `premature convergence' threshold $\hat{R}<1.1$. [Note that in the Bird's Eye model, this struggle is slightly exacerbated by the inclusion of governorates from Israel and Saud Arabia, for which we have no observations, and whose effects are fully interpolated via the spatial process.]

Appendix Figures <ref>, <ref>, and <ref> present five measures of Effective Sample Size (ESS), which tell us about the true number of independent draws from the joint posterior distribution after accounting for auto-correlation within chains. The measures check that the independent sample is `large enough' to ensure stability of summaries of the distribution at various moments (e.g. overall, at the median, at the tails, etc.). Mirroring the performance of the $\hat{R}$, the posterior samples for most of our estimates parameters are well above the recommended threshold ($ESS>400$) for ensuring stability of the central and tail estimates.

We further explore convergence at different quantiles of the posterior distribution of our least-convergent parameters – those with the lowest bulk and tail ESS (Figure <ref> in the Appendix presents these measures for the Bird's Eye model – the Worm's Eye analysis is available upon request).
The inference is that if these relatively low-ESS parameters showcase satisfactory ESS at every quantile, we can be reassured that the whole model has converged. These plots suggest broad reliability of estimates at every section of the distribution. Finally, we explore the mixing properties of our chains for these least-convergent parameters (Figure <ref> in the Appendix).[Here, we choose to include mixing diagnostics for the Bird's Eye model. The Worm's Eye analysis is available upon request.]
These plots broadly suggest good mixing properties of our model, even for these relatively inefficient posterior samples.

§.§ Fixed and Random Effects

Figure <ref> presents the posterior density of the individual-level fixed effects in the Bird's Eye model; Figures <ref> and <ref> present the Worm's Eye equivalent. These plots contain the main results of our models. Note that all the covariates, including dummies, are centered and scaled, hence the coefficients are to be interpreted in terms of standard deviations from the mean of each covariate (Figures <ref>, <ref> and <ref> in the Appendix are the individual-level posterior densities on the original, non-standardized scale). Since we are principally interested in the robust estimation of individual-level predictors, we display only the posterior density of individual fixed effects for all of our models.[Figures <ref> and <ref> display the standardized district-level posterior densities, whilst <ref> and <ref> present district-level coefficients on the original, non-standardized scale.] To aid with interpretation, the mean and standard deviation of each covariate are reported in the legend of each plot.

Posterior density of fixed-effect coefficients for the Bird's Eye model.

Posterior density of fixed-effect coefficients for the Worm's Eye models.

The estimated intercepts for the three models are extremely low. For the Bird's Eye model, the log-odds are in the order of negative $11$. For the Egypt Worm's Eye model, it is just over negative $13$; in Tunisia it is negative $9$. The size of the intercept is primarily driven by the size of the offset, which is in turn determined by the overall prevalence of recruitment. It is therefore not surprising that Egypt's intercept is so dramatically low, given the close-to-zero prevalence of recruitment when compared to population size ($\pi = \frac{4}{100,000}$) versus Tunisia where this prevalence is higher ($\pi = \frac{2}{1,000}$). For the Bird's Eye model, a different offset is provided for observations coming from different countries, to account for country-specific prevalence. The large and negative intercept underscores an important challenge in the explanation of why individuals join movements like ISIS: a linear combination of features capable of pushing an individual to become a recruit has to be extremely large, on the log-odds scale, to meaningfully affect the otherwise extremely low probability of recruitment.

For the purposes of illustration, we focus primarily on testing the role of education and social status in an individual's decision to join ISIS. An individual who has university education and low-status is assumed to be relatively deprived. We compare predicted log-odds, as opposed to predicted probabilities, as these are scarcely comparable due to the powerful effect of the intercept, which drags probabilities of most profiles close to zero (though see Figures <ref>,
<ref>, and <ref> for predicted probabilities of recruitment relative to the `average' profile, and Figures <ref>, <ref>, and <ref> for expected counts under different relative-deprivation scenarios). The total logit effects on probability of recruitment for different relative-deprivation profiles are shown in Figure <ref> for the Bird's Eye model, and in Figures <ref> and <ref> for the Worm's Eye. The relative deprivation hypothesis finds mixed support: at the Bird's eye level, it seems clear that being high-status plays a key role in increasing propensity of being recruited, while having university education plays a more minor role. A similar pattern is evident in Tunisia, though the effect of being high-status and having university education is starker, meaningfully increasing the propensity to join ISIS by around $3$ points on the log-odds scale compared to relatively deprived individuals. In Egypt the effects are more consistent with the traditional relative-deprivation hypothesis, though note the large prediction intervals around the total effects of relatively-deprived individuals. There is also substantial overlap between the distributions in all plots. This is largely due to the uncertainty around the intercept, which plays a role in marginalising these effects. Note further that varying prediction intervals on the effects reflect the highly unbalanced prevalence of the groups in our study. All in all, the evidence from these analyses suggest that high-status individuals were more likely to be recruited by ISIS, and that having university education on-top of a high social status further increases the likelihood of recruitment. The large prediction intervals as a result of uncertainty around the intercept underscores that much remains unknown about the underlying systematic determinants of recruitment.

Predicted propensity of recruitment for relative-deprivation profiles according the Bird's Eye model.

Predicted propensity of recruitment for relative-deprivation profiles according to the Worm's Eye models.

To fit the ICAR model we implemented the fully-connected graph shown in Figure <ref>.[In the graph, edges connect nodes identified by the centroids of governorates for each country. Minor adjustments were performed to ensure the absence of islands or sub-graphs, which would have made the analysis needlessly complicated. Note also that Israel and Saudi Arabia are included for the purpose of obtaining this fully-connected graph, but no observations were available for either of these countries in terms of recruits or Arab Barometer observations, and hence the direction of the estimates for their governorates is entirely driven by the spatial component. Figure <ref> displays the observed number of recruits per area, along with the residual for each governorate.] The spatially autocorrelated component dominates the governorate-level variance, as shown by the posterior of mixing parameter $\lambda$ in Figure <ref>, estimated via Monte-Carlo mean at close to $0.9$, suggesting around $90\%$ of the variance at the governorate level can be explained by the ICAR model.[The spatial distribution of point estimates for governorate and country-level random effects are presented in Figure <ref>.]

Fully-connected graph for the Bird’s Eye model (a) and Governorate-level variance mixing parameter - $\lambda$ (b).

We repeat these analyses for Egypt and Tunisia. Figure <ref> shows similar mixing among spatial and non-spatial components for the two countries, with around $15\%$ of the district-level variance in Egypt being explained by spatial patterns, and $19\%$ in Tunisia. It is noteworthy that very few of our contextual variables have explanatory power for predicting recruitment. Coupled with the low percentage of variance being explained by the spatial components, our analysis suggests that, in spite of our best efforts to account for observable heterogeneity, there exist a vast array of unobserved, non-spatial district-level effects, which accounts for over $80\%$ of the unexplained district-level variance in both Egypt and Tunisia. Hence this contextual variance, while properly accounted for, remains unexplained. In the Appendix we also describe Moran's I statistics for the Worm's Eye analysis as well as point estimates for the district and governorate effects in Egypt and Tunisia.

District-level variance mixing parameter - $\lambda$ - for Egypt (a) and Tunisia (b).

§.§ Predicted propensity of recruitment by profile

To conclude our analysis, we present inferences derived from the posterior predictive distribution of the out-of-sample probability of recruitment, focusing on individual-level characteristics. We produce descriptive plots on the log-odds scale, to aid interpretation and avoid the intercept-trap which would collapse all propensities to zero; plots on the relative probability and absolute probability are also presented in the Appendix, enabling the reader to appreciate the sharp non-linear `jump' which characterises the recruitment profiles of different socio-demographc groups. Note further that these predicted probabilities differ from the within-sample propensities, in that they are calculated after removing the offset, hence generalized to apply to any sampling protocol, as opposed to the specific case-control setting of our application.

What is the profile of individuals `at risk' of recruitment to ISIS according to our models? We attempt to answer this question by analyzing the predicted probabilities of all possible theoretical profiles, defined by the individual-level characteristics available in our data. Every profile is assumed to come from a hypothetical `average district'. Figure <ref> presents point estimates and prediction intervals for the log-odds of recruitment, over $160$ possible profiles in the Bird's eye model. Similar plots displaying the absolute and relative probabilities of recruitment are available in Figures <ref> and <ref> in the Appendix. Table <ref> presents the profiles of the top $10$ most likely profiles to be recruited, providing four useful metrics to interpret the results: predicted probability; predicted rate per $10,000$ people; predicted odds relative to the average profile; and log-odds.

Distribution of the predicted probabilities prediction intervals, presented on the log-odds scale to aid cross-profile comparisons. The black dotted line highlights the the zero-log-odds point, whilst the purple dotted line notes a central estimate for the median recruitment propensity across profiles.

Rank     Married     Student     University Edu     Low-Status     Age     $\hat{\mbox{P}}(r = 1 \mid X = x)$     Predicted Rate     $\frac{\hat{\mbox{P}}(r = 1 \mid X = x)}{\hat{\mbox{P}}(r = 1 \mid X = \bar{X})}$     $\mbox{logit}\left(\hat{\mbox{P}}(r = 1 \mid X = x)\right)$

1           0           0           1           0          18     0.000460     5/10000     22.561     -7.684

2           1           0           1           0          18     0.000398     4/10000     19.273     -7.828

3           0           0           0           0          18     0.000376     4/10000     18.632     -7.885

4           1           0           0           0          18     0.000324     3/10000     16.013     -8.036

5           0           0           1           0          26     0.000281     3/10000     13.745     -8.177

6           1           0           1           0          26     0.000241     2/10000     11.763     -8.331

7           0           0           0           0          26     0.000230     2/10000     11.287     -8.379

8           1           0           0           0          26     0.000197     2/10000     9.695     -8.530

9           0           1           1           0          18     0.000168     2/10000     8.295     -8.690

10           1           1           1           0          18     0.000147     1/10000     7.211     -8.828

Top 10 recruitable theoretical profiles according to the Bird's eye model. Profiles are ordered by predicted probability of recruitment net of sampling protocol. $10$ ages are evaluated, starting at $18$ (to avoid non-existent profiles) and ending at the largest observed age ($86$). The last four columns represent respectively: i. the predicted probability of recruitment; ii. the predicted rate of recruitment per $10,000$ people; iii. the predicted odds of recruitment, relative to the `average' profile; iv. the log-odds of recruitment.

From the Bird's Eye prediction intervals we notice that the predicted probability of recruitment is centered around $-15$ on the log-odds scale, again underscoring how rare the event in question is, and how few are the profiles of people who are susceptible to recruitment. A select number of profiles approach a predicted probability around $-7$, and translate to meaningful rates of recruitment; these are highlighted in the predicted-probabilities table, which show the $10$ most recruitable profiles. Looking at Table <ref>, we can say that the most likely recruitable profile (loosely characterised as a young, high-status, Sunni male with some university education who is unmarried and not currently studying), is around $23$ times as likely to be recruited as an average Sunni male from an average area in the MENA. Moreover, for every $10,000$ members of the most recruitable profile across the region, we expect $5$ to have joined ISIS.
It is worthwhile to note that, in agreement with Figure <ref>, all the most recruitable profiles are high-status individuals, and a majority of them has had some university education. Unsurprisingly, all of these profiles are under-25, and not currently studying.

The individual-level profiles highlighted at the Bird's eye level are comparable to those highlighted within Tunisia in the Worm's eye level (Figure <ref> and Table <ref>), whilst the Egypt analysis reveals more evidence for the relative-deprivation hypothesis, with a majority of the highly-recruitable profiles being relatively-deprived (Figure <ref> and Table <ref>.[For absolute and relative probabilities of recruitment from the Worm's eye models, see Appendix Figures <ref>, <ref>, <ref>, and <ref>.] The relative recruitability of the most susceptible profiles in Egypt and Tunisia is also greater. In Egypt, the most recruitable profile (loosely characterised as a young, low-status, Sunni male with some university education who is married and is currently studying) is around $157$ times as likely to be recruited as the average Egyptian Sunni male – though the expected rate of recruitment for this group is lower than for the Bird's Eye models' most recruitable profile, clocking in at $3$ recruits for every $10,000$ such individuals. The Egypt-specific propensity of recruitment is dramatically lower than that of Tunisia, again highlighting the role of contextual effects. In Tunisia, the most recruitable profile  (loosely characterised as a young, high-status, Sunni male who has university education, is unmarried and is not currently studying) has a probability of recruitment equivalent to $0.04$ on the log-odds scale. This profile is over $335$ times as likely as the average Tunisian Sunni male to be recruited, highlighting that though recruitment is still relatively rare in the population, is significantly concentrated in the top recruitment profiles. This is also supported by Figure <ref>, which shows only a handful of profiles have predicted probabilities above $\frac{1}{100}$.

Worm's Eye (Tunisia) distribution of the predicted probabilities prediction intervals, presented on the log-odds scale to aid cross-profile comparisons. The black dotted line highlights the the zero-log-odds point, whilst the purple dotted line notes a central estimate for the median recruitment propensity across profiles.

Rank     Married     Student     University Edu     Low-Status     Age     $\hat{\mbox{P}}(r = 1 \mid X = x)$     Predicted Rate     $\frac{\hat{\mbox{P}}(r = 1 \mid X = x)}{\hat{\mbox{P}}(r = 1 \mid X = \bar{X})}$     $\mbox{logit}\left(\hat{\mbox{P}}(r = 1 \mid X = x)\right)$

1           0           0           1           0          18     0.043680     437/10000     335.663     -3.086

2           1           0           1           0          18     0.031493     315/10000     242.502     -3.426

3           0           0           1           0          26     0.028925     289/10000     222.013     -3.514

4           1           0           1           0          26     0.021177     212/10000     159.752     -3.833

5           0           0           0           0          18     0.017397     174/10000     137.069     -4.034

6           1           0           0           0          18     0.012615     126/10000     98.539     -4.360

7           0           0           0           0          26     0.011647     116/10000     88.370     -4.441

8           0           0           1           0          33     0.009503     95/10000     73.103     -4.647

9           0           1           1           0          18     0.008662     87/10000     68.471     -4.740

10           1           0           0           0          26     0.008327     83/10000     64.331     -4.780

Top $10$ recruitable theoretical profiles according to the Tunisia `Worm's Eye' model. Profiles are ordered by predicted probability of recruitment net of sampling protocol. $10$ ages are evaluated, starting at $18$ (to avoid non-existent profiles) and ending at the largest observed age ($86$). The last four columns represent respectively: i. the predicted probability of recruitment; ii. the predicted rate of recruitment per $10,000$ people; iii. the predicted odds of recruitment, relative to the `average' profile; iv. the log-odds of recruitment

Worm's Eye (Egypt) distribution of the predicted probabilities prediction intervals, presented on the log-odds scale to aid cross-profile comparisons. The black dotted line highlights the the zero-log-odds point, whilst the purple dotted line notes a central estimate for the median recruitment propensity across profiles.

Rank     Married     Student     University Edu     Low-Status     Age     $\hat{\mbox{P}}(r = 1 \mid X = x)$     Predicted Rate     $\frac{\hat{\mbox{P}}(r = 1 \mid X = x)}{\hat{\mbox{P}}(r = 1 \mid X = \bar{X})}$     $\mbox{logit}\left(\hat{\mbox{P}}(r = 1 \mid X = x)\right)$

1           1           1           1           1          18     0.000287     3/10000     156.721     -8.155

2           0           1           1           1          18     0.000258     3/10000     141.021     -8.264

3           1           1           1           1          25     0.000173     2/10000     97.486     -8.660

4           0           1           1           1          25     0.000155     2/10000     87.536     -8.772

5           1           0           1           1          18     0.000113     1/10000     62.422     -9.091

6           1           1           1           0          18     0.000112     1/10000     62.529     -9.096

7           1           1           1           1          31     0.000108     1/10000     59.223     -9.130

8           0           0           1           1          18     0.000101     1/10000     56.978     -9.199

9           0           1           1           0          18     0.000098     1/10000     53.336     -9.229

10           0           1           1           1          31     0.000097     1/10000     52.977     -9.244

Top 10 recruitable theoretical profiles according to the Egypt `Worm's Eye' model. Profiles are ordered by predicted probability of recruitment net of sampling protocol. $10$ ages are evaluated, starting at $18$ (to avoid non-existent profiles) and ending at the largest observed age ($78$). The last four columns represent respectively: i. the predicted probability of recruitment; ii. the predicted rate of recruitment per $10,000$ people; iii. the predicted odds of recruitment, relative to the `average' profile; iv. the log-odds of recruitment

§ CONCLUSION

Extreme forms of political behaviour are rarely ever committed by more than a tiny subsection of any given national population. Despite their small size, these groups often have an outsized influence on state and international politics. Because of their small size, extremists are particularly hard to study using conventional statistical methods and research designs.

To address this, we propose that extremism researchers adopt a new variant of the case-control design taken from epistemology. This approach allows us to combine survey techniques with ecological forms of analysis, allowing for meaningful comparisons with the underlying populations from which recruits are drawn. To implement this, we solve a number of statistical problems when explaining rare and extreme forms of political behaviour. In particular, we build on existing contaminated case-control designs to demonstrate: 1) how best to incorporate area-level random effects when the number of recruits for a given unit is small; 2) how to account for spatial autocorrelation in this setup; 3) how to regularize coefficients to guard against separation; and 4) provide the  software package so that extremism researchers working on a range of different cases can easily apply our models. Following this approach, we are able to robustly recover individual fixed effects for relevant demographic characteristics net of context, as well as the proportion of unexplained variance attributable to unobserved spatial effects. Finally, we provide simulations to help benchmark our approach and inform the model selection of future researchers.

We demonstrate the potential of our technique with data on ISIS recruits across nine countries in the MENA. Here, the paper makes an important empirical contribution to the literature on Islamist extremism. We find confirmation that recruits tend to have a high social status and are more likely to be highly educated – and we have mixed support for the relative deprivation hypothesis. Despite the considerable advances our approach represents, we are still constrained by the sparsity of information available for certain countries in our sample. Future work might extend our analysis by introducing contextual information across all countries for which have a sufficient number of records for recruits.

Introducing finer-grain contextual data is especially appropriate given that we find a large proportion of variance in the incidence of recruitment to Islamist extremism is attributable to unexplained ecological and spatial confounding. A broader implication here is that future work should focus on obtaining precise estimates of potential local-level confounders. Often, we have only limited data on the characteristics of individual cases, and the characteristics of controls are naturally constrained by the richness of the information we have for our cases. In this regard, introducing richer ecological measures, which is often more easily accessible, could yield important insights, by illuminating the role of contextual factors net of individual-level drivers. Here, promising avenues include the measurement of sub- and cross-national recruitment networks (e.g., ), as well as the use of nontraditional data sources to derive granular, time-varying ecological data (e.g., ).