Evaluating A Program On Hr With Spatially Structured Performance Data

The Data

The time period under study comprises 8 years from 2008 to 2015. In the summer of 2011, a company’s HR management has introduced a special training programme (with voluntary participation) to improve sales skills and ability to approach potential customers of sales managers. The study area consists of 32 London boroughs, each assigned to a sales manager. For each manager, the following key variables have been measured: Gender (Female or Male), Training (No=if the manager didn’t participate in the programme; Yes=if the manager participate in the programme), Age (in years). For each district (associated with each manager), the annual Sales for the 8 years 2008-2015 have been recorded.

The Statistical Model

The observed output Sales is a counting variable and we assume that it has a Poisson distribution:

$$Sales~Poisson\left(E\bullet \rho \right)$$

where the mean is defined in terms of the rate $\rho$ and the expected number $E$ for each specific areal unit. $E$ plays the role of an offset in the model and is a known quantity for each areal unity that should not be estimated. For more information about ways for computing $E$ the interested reader are referred to Congdon (2017), Elliott et al., (2000), and Lesaffre and Lawson (2012).

In the first considered model, from now on called model1 , an additive structure including the categorical covariates Training , Gender and Age and a nonlinear effect of Year is specified:

$$\log \left(\rho \right)={\beta }_0+{\beta }_1\bullet Training+{\beta }_2\bullet Gender+{\beta }_3\bullet Age+f_1\left(Year\right).$$

The log link is the natural choice for the Poisson model. The $\beta$’s are the regression coefficients of the linear effect and $f_1$ is the unknown function of the nonlinear effect in Year. ${\beta }_0$ is thought as the average outcome rate in the entire study region. We note that model1 does not take into account the spatial dependence amongst the areal units.

To take into account the spatial piece of information, in a second considered model, model2 , an unknown function $f_2$ of the spatial effect has been added:

$$\log \left(\rho \right)={\beta }_0+{\beta }_1\bullet Training+{\beta }_2\bullet Gender+{\beta }_3\bullet Age+f_1\left(Year\right)+f_2\left(ID\right)$$

where $ID$ is a vector identifying each of the 32 London boroughs. We take a Besag prior on $f_2$ as specified in the following subsection;

The two considered models, model1 and model2 , have been compared on the basis of their DIC and WAIC scores (see subsection 3.4 below).

A third model, model3 , has been considered to evaluate the effectiveness of the training activity. The total $Salesaftertraining$ has been modelled as a Poisson distribution:

$Salesaftertraining\mathrm{}~Poisson\left(4E\bullet \rho \right)$

where the mean is defined in terms of the rate $\rho$ and the expected total number $4E$ for each specific areal unit for the four considered years after training. A log link with additive structure is specified:

$$\log \left(\rho \right)={\beta }_0+{\beta }_1\bullet Training+{\beta }_2\bullet MPS+{\beta }_3\bullet Training\bullet MPS+f\left(ID\right)$$

where the log rate $\log \left(\rho \right)$ is a linear function of $Training$ , $MPS$ (Mean Proportional Sales before training) and their interaction. An unknown function $f$ of the spatial effect with Besag spatial structure has been also added to the model.

The Besag Spatial Model

Our data set presents an areal spatial structure, where the observations are related to geographical regions (London boroughs) with adjacency information. We assume a Besag spatial structure as defined below (Wang, Yue, & Faraway, 2018).

A Gaussian increment is defined between neighbouring regions (regions that share a common border) i and j as:

$$f\left({\mathit{x}}_i\right)-\mathit{}f\left({\mathit{x}}_j\right)\mathit{}~\mathit{}N\left(0,\frac{{\sigma }_f^2}{w_{ij}}\right)$$

where ${\mathit{x}}_i$ and ${\mathit{x}}_j$ are the centroids of the regions and $w_{ij}$ are the positive and symmetric weights.

Assuming the increments are independent, it can be shown that the full conditional distribution of $f\left({\mathit{x}}_i\right)$ is normal:

$f\left({\mathit{x}}_i\right)|\mathit{}f\left({\mathit{x}}_{-i}\right),\mathit{}\mathit{\tau }\mathit{}~N\left(\frac{{\sum }_{j~i}{w}_{ij}f\left({\mathit{x}}_j\right)}{{\sum }_{j~i}{w}_{ij}},\frac{{\sigma }_f^2}{{\sum }_{j~i}{w}_{ij}}\right)$.

This prior is called Besag model because it is a special case of the intrinsic autoregressive models introduced by Besag and Kooperg (1995).

More specific or complex spatial structures can be adopted, such as the Besag-York-Molliè model. More details can be found in Besag, York & Mollie (1991) and Blangiardo and Cameletti (2015).

The DIC and WAIC

Model comparison have been based on DIC and WAIC. DIC stands for Deviance Information Criteria (Spiegelhalter, Best, Carlin, & Van der Linde, 2002). It is a score assigned to a model that take into account how good is the model fit and the effective number of parameters, that is the complexity of the model. The less the DIC the best the model. WAIC stand for Watanabe-Akaike Information Criteria (or Widely Applicable Information Criteria) and it has been introduced in Watanabe (2010). In some sense, WAIC is an improvement of DIC (Gelman et al., 2013), but for simple models they give very similar scores.

Results of preliminary Explorative Data Analysis.

In the preliminary Explorative Data Analysis no outlier where found in the data set. From the comparative violin plot in Figure 1 , it clearly appears that yearly Proportional Sales does not depend on Gender but does depend on Training.

See Full Size >

Results of the Bayesian Statistical data analysis.

Instead of MCMC simulations, the INLA approach used in this study provides approximations to the posterior marginal distributions of the parameters. Summaries of these posterior distributions include posterior means and 95% Credible Intervals (CI), which play in Bayesian statistics the same role of the maximum likelihood estimates and 95% Confidence Intervals in classical Frequentist statistics.

The estimated distribution of the parameters of model1 are summarized as follows for each categorical covariate in the model: $Training$ (mean=0.364, 95% CI= [0.346,0.382]), $Gender$ (mean=0.079, 95% CI= [0.061,0.097]), $Age$ (mean=-0.021, 95% CI= [-0.025,-0.018]). Note that all CIs don’t contain 0. It means that all the related coefficients should be considered significantly different from 0. That is all covariates in model1 appears to be relevant drivers of Sales.

From the results of model2 fitting, only the CI related to $Training$ does not contain 0 (mean=0.426, 95% CI= [0.260,0.591]), that is $Training$ is the only significant driver of Sales.

In model comparison, $model2$ (DIC=4619.337, WAIC=5038.290) appears to be better than $model1$ (DIC=6468.970, WAIC=6673.552).

Maps may be a very useful tools for representing areal data. In Figure 2 maps of model2 fitted values are shown for each year. In the middle map the white units are those whose manager has not adopted the new approach. The non white areal units are related to managers that followed the training programme and adopted the new approach. It is easy to find units in the maps with low performance (green area) in the first four years (2008-2011), whose sales manager didn’t adopted the programme and with low performance also in the following 4 years (2012-2015). There are also units that improved their performance from green to orange because the related manager followed the training programme and adopted the new approach. This kind of maps may be very effective in showing the impact of the training programme on each areal unit.

The estimated distribution of the parameters of model3 are summarized as follows for each covariate in the model: $MSB$ (mean=10.288, CI=[9.307,11.289]), $Training$ (mean=0.734, CI=[0.602,0.861]), $Interaction$ (mean=2.506, CI=[1.007,4.076]). Note that all CIs don’t contain 0.

See Full Size >