Regression Discontinuity Design (RDD) with stochtree

Introduction

We study conditional average treatment effect (CATE) estimation for regression discontinuity designs (RDD), in which treatment assignment is based on whether a particular covariate --- referred to as the running variable --- lies above or below a known value, referred to as the cutoff value. Because treatment is deterministically assigned as a known function of the running variable, RDDs are trivially deconfounded: treatment assignment is independent of the outcome variable, given the running variable (because treatment is conditionally constant). However, estimation of treatment effects in RDDs is more complicated than simply controlling for the running variable, because doing so introduces a complete lack of overlap, which is the other key condition needed to justify regression adjustment for causal inference. Nonetheless, the CATE at the cutoff, $X=c$, may still be identified provided the conditional expectation $E[Y \mid X,W]$ is continuous at that point for all $W=w$. We exploit this assumption with the leaf regression BART model implemented in Stochtree, which allows us to define an explicit prior on the CATE. We now describe the RDD setup and our model in more detail, and provide code to implement our approach.

Regression Discontinuity Design

We conceptualize the treatment effect estimation problem via a quartet of random variables $(Y, X, Z, U)$. The variable $Y$ is the outcome variable; the variable $X$ is the running variable; the variable $Z$ is the treatment assignment indicator variable; and the variable $U$ represents additional, possibly unobserved, causal factors. What specifically makes this correspond to an RDD is that we stipulate that $Z = I(X > c)$, for cutoff $c$. We assume $c = 0$ without loss of generality.
The following figure depicts a causal diagram representing the assumed causal relationships between these variables. Two key features of this diagram are one, that $X$ blocks the impact of $U$ on $Z$: in other words, $X$ satisfies the back-door criterion for learning causal effects of $Z$ on $Y$. And two, $X$ and $U$ are not descendants of $Z$.

Using this causal diagram, we may express $Y$ as some function of its graph parents, the random variables $(X,Z,U)$: $$Y = F(X,Z,U).$$ In principle, we may obtain draws of $Y$ by first drawing $(X,Z,U)$ according to their joint distribution and then applying the function $F$. Similarly, we may relate this formulation to the potential outcomes framework straightforwardly: \begin{equation} \begin{split} Y^1 &= F(X,1,U),\\ Y^0 &= F(X,0,U). \end{split} \end{equation} Here, draws of $(Y^1, Y^0)$ may be obtained (in principle) by drawing $(X,Z,U)$ from their joint distribution and using only the $(X,U)$ elements as arguments in the above two equations, "discarding" the drawn value of $Z$. Note that this construction implies the consistency condition: $Y = Y^1 Z + Y^0 ( 1 - Z)$. Likewise, this construction implies the no interference condition because each $Y_i$ is considered to be produced with arguments ($X_i, Z_i, U_i)$ and not those from other units $j$; in particular, in constructing $Y_i$, $F$ does not take $Z_j$ for $j \neq i$ as an argument.

Next, we define the following conditional expectations \begin{equation} \begin{split} \mu_1(x) &= E[ F(x, 1, U) \mid X = x] ,\\ \mu_0(x) &= E[ F(x, 0, U) \mid X = x], \end{split} \end{equation} with which we can define the treatment effect function $$\tau(x) = \mu_1(x) - \mu_0(x).$$ Because $X$ satisfies the back-door criterion, $\mu_1$ and $\mu_0$ are estimable from the data, meaning that \begin{equation} \begin{split} \mu_1(x) &= E[ F(x, 1, U) \mid X = x] = E[Y \mid X=x, Z=1],\\ \mu_0(x) &= E[ F(x, 0, U) \mid X = x] = E[Y \mid X=x, Z=0], \end{split} \end{equation} the right-hand-sides of which can be estimated from sample data, which we supposed to be independent and identically distributed realizations of $(Y_i, X_i, Z_i)$ for $i = 1, \dots, n$. However, because $Z = I(X >0)$ we can in fact only learn $\mu_1(x)$ for $X > 0$ and $\mu_0(x)$ for $X < 0$. In potential outcomes terminology, conditioning on $X$ satisfies ignorability, $$(Y^1, Y^0) \perp \!\!\! \perp Z \mid X,$$ but not strong ignorability, because overlap is violated. Overlap would require that $$0 < P(Z = 1 \mid X=x) < 1 \;\;\;\; \forall x,$$ which is clearly violated by the RDD assumption that $Z = I(X > 0)$. Consequently, the overall ATE, $\bar{\tau} = E(\tau(X)),$ is unidentified, and we must content ourselves with estimating $\tau(0)$, the conditional average effect at the point $x = 0$, which we estimate as the difference between $\mu_1(0) - \mu_0(0)$. This is possible for continuous $X$ so long as one is willing to assume that $\mu_1(x)$ and $\mu_0(x)$ are both suitably smooth functions of $x$: any inferred discontinuity at $x = 0$ must therefore be attributable to treatment effect.

Conditional average treatment effects in RDD

We are concerned with learning not only $\tau(0)$, the "RDD ATE" (e.g. the CATE at $x = 0$), but also RDD CATEs, $\tau(0, \mathrm{w})$ for some covariate vector $\mathrm{w}$. Incorporating additional covariates in the above framework turns out to be straightforward, simply by defining $W = \varphi(U)$ to be an observable function of the (possibly unobservable) causal factors $U$. We may then define our potential outcome means as \begin{equation} \begin{split} \mu_1(x,\mathrm{w}) &= E[ F(x, 1, U) \mid X = x, W = \mathrm{w}] = E[Y \mid X=x, W=\mathrm{w}, Z=1],\\ \mu_0(x,\mathrm{w}) &= E[ F(x, 0, U) \mid X = x, W = \mathrm{w}] = E[Y \mid X=x, W =\mathrm{w}, Z=0], \end{split} \end{equation} and our treatment effect function as \begin{equation} \tau(x,\mathrm{w}) = \mu_1(x,\mathrm{w}) - \mu_0(x,\mathrm{w}) \end{equation} We consider our data to be independent and identically distributed realizations $(Y_i, X_i, Z_i, W_i)$ for $i = 1, \dots, n$. Furthermore, we must assume that $\mu_1(x,\mathrm{w})$ and $\mu_0(x,\mathrm{w})$ are suitably smooth functions of $x$, {\em for every} $\mathrm{w}$; in other words, for each value of $\mathrm{w}$ the usual continuity-based identification assumptions must hold.

With this framework and notation established, CATE estimation in RDDs boils down to estimation of condition expectation functions $E[Y \mid X=x, W=\mathrm{w}, Z=z]$, for which we turn to BART models.

The BARDDT Model

We propose a BART model where the trees are allowed to split on $(x,\mathrm{w})$ but where each leaf node parameter is a vector of regression coefficients tailored to the RDD context (rather than a scalar constant as in default BART). In one sense, such a model can be seen as implying distinct RDD ATE regressions for each subgroup determined by a given tree; however, this intuition is only heuristic, as the entire model is fit jointly as an ensemble of such trees. Instead, we motivate this model as a way to estimate the necessary conditional expectations via a parametrization where the conditional treatment effect function can be explicitly regularized, as follows.

Let $\psi$ denote the following basis vector: \begin{equation} \psi(x,z) = \begin{bmatrix} 1 & z x & (1-z) x & z \end{bmatrix}. \end{equation} To generalize the original BART model, we define $g_j(x, \mathrm{w}, z)$ as a piecewise linear function as follows. Let $b_j(x, \mathrm{w})$ denote the node in the $j$th tree which contains the point $(x, \mathrm{w})$; then the prediction function for tree $j$ is defined to be: \begin{equation} g_j(x, \mathrm{w}, z) = \psi(x, z) \Gamma_{b_j(x, \mathrm{w})} \end{equation} for a leaf-specific regression vector $\Gamma_{b_j} = (\eta_{b_j}, \lambda_{b_j}, \theta_{b_j}, \Delta_{b_j})^t$. Therefore, letting $n_{b_j}$ denote the number of data points allocated to node $b$ in the $j$th tree and $\Psi_{b_j}$ denote the $n_{b_j} \times 4$ matrix, with rows equal to $\psi(x,z)$ for all $(x_i,z_i) \in b_j$, the model for observations assigned to leaf $b_j$, can be expressed in matrix notation as: \begin{equation} \begin{split} \mathbf{Y}_{b_j} \mid \Gamma_{b_j}, \sigma^2 &\sim \mathrm{N}(\Psi_{b_j} \Gamma_{b_j},\sigma^2)\\ \Gamma_{b_j} &\sim \mathrm{N} (0, \Sigma_0), \end{split} \end{equation} where we set $\Sigma_0 = \frac{0.033}{J} \mathbf{I}$ as a default (for $x$ vectors standardized to have unit variance in-sample).

This choice of basis entails that the RDD CATE at $\mathrm{w}$, $\tau(0, \mathrm{w})$, is a sum of the $\Delta_{b_j(0, \mathrm{w})}$ elements across all trees $j = 1, \dots, J$: \begin{equation} \begin{split} \tau(0, \mathrm{w}) &= E[Y^1 \mid X=0, W = \mathrm{w}] - E[Y^0 \mid X = 0, W = \mathrm{w}]\\ & = E[Y \mid X=0, W = \mathrm{w}, Z = 1] - E[Y \mid X = 0, W = \mathrm{w}, Z = 0]\\ &= \sum_{j = 1}^J g_j(0, \mathrm{w}, 1) - \sum_{j = 1}^J g_j(0, \mathrm{w}, 0)\\ &= \sum_{j = 1}^J \psi(0, 1) \Gamma_{b_j(0, \mathrm{w})} - \sum_{j = 1}^J \psi(0, 0) \Gamma_{b_j(0, \mathrm{w})} \\ & = \sum_{j = 1}^J \Bigl( \psi(0, 1) - \psi(0, 0) \Bigr) \Gamma_{b_j(0, \mathrm{w})} \\ & = \sum_{j = 1}^J \Bigl( (1,0,0,1) - (1,0,0,0) \Bigr) \Gamma_{b_j(0, \mathrm{w})} \\ &= \sum_{j=1}^J \Delta_{b_j(0, \mathrm{w})}. \end{split} \end{equation} As a result, the priors on the $\Delta$ coefficients directly regularize the treatment effect. We set the tree and error variance priors as in the original BART model.

The following figures provide a graphical depiction of how the BARDDT model fits a response surface and thereby estimates CATEs for distinct values of $\mathrm{w}$. For simplicity only two trees are used in the illustration, while in practice dozens or hundreds of trees may be used (in our simulations and empirical example, we use 150 trees).

Two regression trees with splits in x and a single scalar w. Node images depict the g(x,w,z) function (in x) defined by that node's coefficients. The vertical gap between the two line segments in a node that contain x=0 is that node's contribution to the CATE at X = 0. Note that only such nodes contribute for CATE prediction at x=0

The two top figures show the same two regression trees as in the preceding figure, now represented as a partition of the x-w plane. Labels in each partition correspond to the leaf nodes depicted in the previous picture. The bottom figure shows the partition of the x-w plane implied by the sum of the two trees; the red dashed line marks point W=w* and the combination of nodes that include this point

Left: The function fit at W = w* for the two trees shown in the previous two figures, shown superimposed. Right: The aggregated fit achieved by summing the contributes of two regression tree fits shown at left. The magnitude of the discontinuity at x = 0 (located at the dashed gray vertical line) represents the treatment effect at that point. Different values of w will produce distinct fits; for the two trees shown, there can be three distinct fits based on the value of w.

An interesting property of BARDDT can be seen in this small illustration --- by letting the regression trees split on the running variable, there is no need to separately define a 'bandwidth' as is used in the polynomial approach to RDD. Instead, the regression trees automatically determine (in the course of posterior sampling) when to 'prune' away regions away from the cutoff value. There are two notable features of this approach. One, different trees in the ensemble are effectively using different local bandwidths and these fits are then blended together. For example, in the bottom panel of the second figure, we obtain one bandwidth for the region $d+i$, and a different one for regions $a+g$ and $d+g$. Two, for cells in the tree partition that do not span the cutoff, the regression within that partition contains no causal contrasts --- all observations either have $Z = 1$ or $Z = 0$. For those cells, the treatment effect coefficient is ill-posed and in those cases the posterior sampling is effectively a draw from the prior; however, such draws correspond to points where the treatment effect is unidentified and none of these draws contribute to the estimation of $\tau(0, \mathrm{w})$ --- for example, only nodes $a+g$, $d+g$, and $d+i$ provide any contribution. This implies that draws of $\Delta$ corresponding to nodes not predicting at $X=0$ will always be draws from the prior, which has some intuitive appeal.

Demo

In this section, we provide code for implementing our model in stochtree on a popular RDD dataset. First, let us load stochtree and all the necessary libraries for our posterior analysis.

import matplotlib.pyplot as plt
import seaborn as sns
import numpy as np
import pandas as pd
from sklearn.tree import DecisionTreeRegressor, plot_tree
from stochtree import BARTModel

Dataset

The data comes from Lindo et al (2010), who analyze data on college students enrolled in a large Canadian university in order to evaluate the effectiveness of an academic probation policy. Students who present a grade point average (GPA) lower than a certain threshold at the end of each term are placed on academic probation and must improve their GPA in the subsequent term or else face suspension. We are interested in how being put on probation or not, $Z$, affects students' GPA, $Y$, at the end of the current term. The running variable, $X$, is the negative distance between a student's previous-term GPA and the probation threshold, so that students placed on probation ($Z = 1$) have a positive score and the cutoff is 0. Potential moderators, $W$, are:

• gender (male),
• age upon entering university (age_at_entry)
• a dummy for being born in North America (bpl_north_america),
• the number of credits taken in the first year (totcredits_year1)
• an indicator designating each of three campuses (loc_campus 1, 2 and 3), and
• high school GPA as a quantile w.r.t the university's incoming class (hsgrade_pct).

# Load and organize data
data = pd.read_csv("https://raw.githubusercontent.com/rdpackages-replication/CIT_2024_CUP/refs/heads/main/CIT_2024_CUP_discrete.csv")

y = data.loc[:,"nextGPA"].to_numpy()
x = data.loc[:,"X"].to_numpy()

x = x/np.std(x)
w = data.iloc[:,3:11]

ordered_cat = pd.api.types.CategoricalDtype(ordered=True)
unordered_cat = pd.api.types.CategoricalDtype(ordered=False)
w.loc[:,"totcredits_year1"] = w.loc[:,"totcredits_year1"].astype(ordered_cat)
w.loc[:,"male"] = w.loc[:,"male"].astype(unordered_cat)
w.loc[:,"bpl_north_america"] = w.loc[:,"bpl_north_america"].astype(unordered_cat)
w.loc[:,"loc_campus1"] = w.loc[:,"loc_campus1"].astype(unordered_cat)
w.loc[:,"loc_campus2"] = w.loc[:,"loc_campus2"].astype(unordered_cat)
w.loc[:,"loc_campus3"] = w.loc[:,"loc_campus3"].astype(unordered_cat)
c = 0
n = data.shape[0]
z = np.where(x > c, 1.0, 0.0)
# Window for prediction sample
h = 0.1
test = (x > -h) & (x < h)
ntest = np.sum(np.where(test, 1, 0))

Target estimand

Generically, our estimand is the CATE function at $x = 0$; i.e. $\tau(0, \mathrm{w})$. The key practical question is which values of $\mathrm{w}$ to consider. Some values of $\mathrm{w}$ will not be well-represented near $x=0$ and so no estimation technique will be able to estimate those points effectively. As such, to focus on feasible points --- which will lead to interesting comparisons between methods --- we recommend restricting the evaluation points to the observed $\mathrm{w}_i$ such that $|x_i| \leq \delta$, for some $\delta > 0$. In our example, we use $\delta = 0.1$ for a standardized $x$ variable. Therefore, our estimand of interest is a vector of treatment effects: $$\tau(0, \mathrm{w}_i) \;\;\; \forall i \;\text{ such that }\; |x_i| \leq \delta$$

Implementing BARDDT

In order to implement our model, we write the Psi vector, as defined before: Psi = np.column_stack([np.ones(n), z * x, (1 - z) * x, z]). The training matrix for the model is np.column_stack([x, w]), which we feed into the BARTModel sampler via the X_train parameter. The basis vector Psi is fed into the function via the leaf_basis_train parameter. The parameter list barddt_mean_params defines options for the mean forest (a different list can be defined for a variance forest in the case of heteroscedastic BART, which we do not consider here). Importantly, in this list we define parameter sigma2_leaf_init = np.diag(np.repeat(0.1/150, 4)), which sets $\Sigma_0$ as described above. Now, we can fit the model, which is saved in object barddt_model.

Once the model is fit, we need 3 elements to obtain the CATE predictions: the basis vectors at the cutoff for $z=1$ and $z=0$, the test matrix $[X \quad W]$ at the cutoff, and the testing sample. We define the prediction basis vectors $\psi_1 = [1 \quad 0 \quad 0 \quad 1]$ and $\psi_0 = [1 \quad 0 \quad 0 \quad 0]$, which correspond to $\psi$ at $(x=0,z=1)$, and $(x=0,z=0)$, respectively. These vectors are written into Python as Psi1 = np.column_stack([np.ones(n), np.repeat(c, n), np.zeros(n), np.ones(n)]) and Psi0 = np.column_stack([np.ones(n), np.zeros(n), np.repeat(c, n), np.zeros(n)]). Then, we write the test matrix at $(x=0,\mathrm{w})$ as xmat_test = np.column_stack([np.zeros(n), w])[test,:]. Finally, we must define the testing window. As discussed previously, our window is set such that $|x| \leq 0.1$, which can be set in Python as test = (x > -h) & (x < h).

Once all of these elements are set, we can obtain the outcome predictions at the cutoff by running barddt_model.predict(xmat_test, Psi1) (resp. barddt_model.predict(xmat_test, Psi0)). Each of these calls returns a list, from which we can extract element y_hat to obtain the posterior distribution for the outcome. In the code below, the treated and control outcome predictions are saved in the matrix objects pred1 and pred0, respectively. Now, we can obtain draws from the CATE posterior by simply subtracting these matrices. The function below outlines how to perform each of these steps in Python.

def estimate_barddt(y,x,w,z,test,c,num_gfr=10,num_mcmc=100,seed=None):
    ## Lists of parameters for the Stochtree BART function
    barddt_global_params = {
        "standardize": True,
        "sample_sigma_global": True,
        "sigma2_global_init": 0.1
    }
    if seed is not None:
        barddt_global_params["random_seed"] = seed
    barddt_mean_params = {
        "num_trees": 50,
        "min_samples_leaf": 20,
        "alpha": 0.95,
        "beta": 2,
        "max_depth": 20,
        "sample_sigma2_leaf": False,
        "sigma2_leaf_init": np.diag(np.repeat(0.1/150, 4))
    }
    ## Set basis vector for leaf regressions
    n = y.shape[0]
    Psi = np.column_stack([np.ones(n), z * x, (1 - z) * x, z])
    covariates = np.column_stack([x, w])
    ## Model fit
    barddt_model = BARTModel()
    barddt_model.sample(
        X_train=covariates,
        y_train=y,
        leaf_basis_train=Psi,
        num_gfr=num_gfr,
        num_mcmc=num_mcmc,
        general_params=barddt_global_params,
        mean_forest_params=barddt_mean_params
    )
    ## Define basis vectors and test matrix for outcome predictions at X=c
    Psi1 = np.column_stack([np.ones(n), np.repeat(c, n), np.zeros(n), np.ones(n)])
    Psi0 = np.column_stack([np.ones(n), np.zeros(n), np.repeat(c, n), np.zeros(n)])
    Psi1 = Psi1[test,:]
    Psi0 = Psi0[test,:]
    xmat_test = np.column_stack([np.zeros(n), w])[test,:]
    ## Obtain outcome predictions
    pred1 = barddt_model.predict(xmat_test, Psi1)
    pred0 = barddt_model.predict(xmat_test, Psi0)
    ## Obtain CATE posterior
    return pred1 - pred0

Now, we proceed to fit the BARDDT model.

num_chains = 4
num_gfr = 2
num_mcmc = 100
cate_result = np.empty((ntest, num_chains*num_mcmc))
for i in range(num_chains):
    cate_rdd = estimate_barddt(y,x,w,z,test,c,num_gfr=2,num_mcmc=100,seed=i)
    cate_result[:,(i*num_mcmc):((i+1)*num_mcmc)] = cate_rdd

We now proceed to analyze the CATE posterior. The figure produced below presents a summary of the CATE posterior produced by BARDDT for this application. This picture is produced fitting a regression tree, using $W$ as the predictors, to the individual posterior mean CATEs: \begin{equation} \bar{\tau}_i = \frac{1}{M} \sum_{h = 1}^M \tau^{(h)}(0, \mathrm{w}_i), \end{equation} where $h$ indexes each of $M$ total posterior samples. As in our simulation studies, we restrict our posterior analysis to use $\mathrm{w}_i$ values of observations with $|x_i| \leq \delta = 0.1$ (after normalizing $X$ to have standard deviation 1 in-sample). For the Lindo et al (2010) data, this means that BARDDT was trained on $n = 40,582$ observations, of which 1,602 satisfy $x_i \leq 0.1$, which were used to generate the effect moderation tree.

## Fit regression tree
y_surrogate = np.mean(cate_rdd, axis=1)
X_surrogate = w.iloc[test,:]
cate_surrogate = DecisionTreeRegressor(min_impurity_decrease=0.0001)
cate_surrogate.fit(X=X_surrogate, y=y_surrogate)
plot_tree(cate_surrogate, impurity=False, filled=True, feature_names=w.columns, proportion=False, label='root', node_ids=True)
plt.show()

The resulting effect moderation tree indicates that course load (credits attempted) in the academic term leading to their probation is a strong moderator. Contextually, this result is plausible, both because course load could relate to latent character attributes that influence a student's responsiveness to sanctions and also because it could predict course load in the current term, which would in turn have implications for the GPA (i.e. it is harder to get a high GPA while taking more credit hours). The tree also suggests that effects differ by age and gender of the student. These findings are all prima facie plausible as well.

To gauge how strong these findings are statistically, we can zoom in on isolated subgroups and compare the posteriors of their subgroup average treatment effects. This approach is valid because in fitting the effect moderation tree to the posterior mean CATEs we in no way altered the posterior itself; the effect moderation tree is a posterior summary tool and not any additional inferential approach; the posterior is obtained once and can be explored freely using a variety of techniques without vitiating its statistical validity. Investigating the most extreme differences is a good place to start: consider the two groups of students at opposite ends of the treatment effect range discovered by the effect moderation tree:

• Group A a male student that attempted more than 4.8 credits in their first year (rightmost leaf node, colored red, comprising 211 individuals)
• Group B a female student of any gender who entered college younger than 19 (leftmost leaf node, colored deep orange, comprising 369 individuals).

Subgroup CATEs are obtained by aggregating CATEs across the observed $\mathrm{w}_i$ values for individuals in each group; this can be done for individual posterior samples, yielding a posterior distribution over the subgroup CATE: \begin{equation} \bar{\tau}_A^{(h)} = \frac{1}{n_A} \sum_{i : \mathrm{w}_i} \tau^{(h)}(0, \mathrm{w}_i), \end{equation} where $h$ indexes a posterior draw and $n_A$ denotes the number of individuals in the group A.

The code below produces a contour plot for a bivariate kernel density estimate of the joint CATE posterior distribution for subgroups A and B. The contour lines are nearly all above the $45^{\circ}$ line, indicating that the preponderance of posterior probability falls in the region where the treatment effect for Group A is greater than that of Group B, meaning that the difference in the subgroup treatment effects flagged by the effect moderation tree persist even after accounting for estimation uncertainty in the underlying CATE function.

predicted_nodes = cate_surrogate.apply(X=X_surrogate)
posterior_group_a = np.mean(cate_result[predicted_nodes==2,:],axis=0)
posterior_group_b = np.mean(cate_result[predicted_nodes==6,:],axis=0)
posterior_df = pd.DataFrame({'group_a': posterior_group_a, 'group_b': posterior_group_b})
sns.kdeplot(data=posterior_df, x="group_b", y="group_a")
plt.axline((0, 0), slope=1, color="black", linestyle=(0, (3, 3)))
plt.show()

As always, CATEs that vary with observable factors do not necessarily represent a causal moderating relationship. Here, if the treatment effect of academic probation is seen to vary with the number of credits, that does not imply that this association is causal: prescribing students to take a certain number of credits will not necessarily lead to a more effective probation policy, it may simply be that the type of student to naturally enroll for fewer credit hours is more likely to be responsive to academic probation. An entirely distinct set of causal assumptions are required to interpret the CATE variations themselves as causal. All the same, uncovering these patterns of treatment effect variability are crucial to suggesting causal mechanism to be investigated in future studies.

References

Lindo, Jason M., Nicholas J. Sanders, and Philip Oreopoulos. "Ability, gender, and performance standards: Evidence from academic probation." American economic journal: Applied economics 2, no. 2 (2010): 95-117.