Aims of this notebook

The Global Terrorism Database

We will use a dataset that includes variables of 180,000 terrorist attacks between 1970-2017. Note that we excluded some variables for clarity and did some preprocessing on the data.

Suppose you are given that dataset and you’re asked to inform policy-makers about the relationship between terrorist attack details and the number of victims. This notebook asks you to perform initial descriptive statistics and then build and evaluate models of the data.

You can load the dataset as follows:

Variable codebook:

This command loads a dataframe called gtd in your notebook. You can query that dataframe as usual.

Have a look at the names (columns) in the data.frame:

names(gtd)
[1] "eventid" "iyear"   "imonth"  "iday"    "nperps"  "suicide" "ransom"  "nkill"  
[9] "nwound" 

… and show the first 10 rows:

head(gtd, 10)
         eventid iyear imonth iday nperps suicide ransom nkill nwound
1  9.733093e-313  1970      0    0      7       0      1     0      0
2  9.733144e-313  1970      1    2      3       0      0     0      0
3  9.733144e-313  1970      1    2      1       0      0     0      0
4  9.733144e-313  1970      1    3      1       0      0     0      0
5  9.733147e-313  1970      1    8      1       0      0     0      0
6  9.733148e-313  1970      1   11      1       0      0     1      0
7  9.733150e-313  1970      1   15      5       0      0     0      0
8  9.733152e-313  1970      1   19      3       0      0     0      0
9  9.733152e-313  1970      1   19      2       0      0     0      0
10 9.733153e-313  1970      1   20      1       0      0     1      0

Note that we excluded variables (for clarity) and observations (to avoid missing values), so the actual dimensions of this dataframe are:

dim(gtd)
[1] 9147    9

Understanding the data

Let’s start with understanding the data bit better. You’d want to do this to avoid modelling relationships that are not meaningful.

Task

Look at the frequencies of the number of perpetrators and subset these frequency counts by the suicide and ransom variable.

#your code comes here

Task

In which year was the number of perpetrators (on average) the highest?

#your code comes here

Task

What is the most common value of the number of persons killed and wounded?

Hint: ?hist and ?table

#your code comes here

Task

Display the relationship between the number of killed victims and the number of wounded victims in a figure.

What kind of a relationship do you expect?

Hint: ?plot

#your code comes here

Task

What is the mean number of perpetrators when the attack was a suicide attack?

Hint: ?tapply

#your code comes here

Simple regression

We’ll start with simple regression (i.e. one predictor variable):

Task

Build a simple regression model that models the number of wounded victims through the number of perpetrators.

#your code comes here

Task

How satisfied are you with your model? One way to assess the “model fit” is to calculate the root mean square error - RMSE - (residuals). Calculate that metric and think about the meaning of this fit index. What does it tell you and how satisfied are you with it?

#your code comes here

Task

Plot the fitted values (in green) and the observed values (in blue) to assess the model fit visually.

#your code comes here

Multiple regression

Now you might want to use multiple variables in your model:

Task

Build a multiple regression model that models the number of killed victims through the variables suicide and ransom. Include only the two main effects (and let the intercept in the model).

#your code comes here

Task

Have a look at a potential interaction between these two predictor variables. Use the interaction.plot function to look at the joint relationship of these two variables on the number of killed victims.

#your code comes here

What does this graph tell you? Can you identify the main effects and (potential) interaction?

Task

Now look at the interaction in a numerical manner.

Hint: ?tapply

#your code comes here

Task

Suppose you want to expand the model by adding the interaction term to it. Build that model.

#your code comes here

Task

Based on the RMSE of each of the two models above (2 main effects vs 2 main effects + 1 interaction), which one do you prefer?

#your code comes here

Task

Have a look at the distribution of the nperps and nkill column. Are there some potential outliers in there?

Hint: ?plot

#your code comes here

Task

Re-run the best fitting regression model again after exluding the potential outliers. The decision for the “best” model can be made based on the RMSE:

#your code comes here

Model selection

Now suppose you want to let the model building process be decided by a stepwise model selection procedure.

Task

Build a “null model” that only contains an intercept.

#your code comes here

Task

Build a full model for the number of wounded victims modeled through the number of perpetrators, “suicide” and “ransom”.

#your code comes here

Task

Determine the best fitting model in a backward model selection procedure.

#your code comes here

Task

Run the model selection again but this time using the forward model selection procedure.

#your code comes here

Task

Compare the RMSE of the null model, the full model and the best fitting model (if different from the full model).

Display the residuals in a graph using the colours green, black, and blue for the null model, observed values, and the full model, respectively.

Hint: ?points

#your code comes here

Did you expect the observations for the null model? See whether you can discover the reason for that relationship in the model outcome (i.e. the coefficients).

END


LS0tCnRpdGxlOiAiUmVncmVzc2lvbiBtb2RlbGxpbmcgaW4gUiwgUFNNMiIKYXV0aG9yOiAiQiBLbGVpbmJlcmciCnN1YnRpdGxlOiBEZXB0IG9mIFNlY3VyaXR5IGFuZCBDcmltZSBTY2llbmNlLCBVQ0wKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQoKIyMgQWltcyBvZiB0aGlzIG5vdGVib29rCgoKIyMgCgoKIyMgVGhlIEdsb2JhbCBUZXJyb3Jpc20gRGF0YWJhc2UKCldlIHdpbGwgdXNlIGEgZGF0YXNldCB0aGF0IGluY2x1ZGVzIHZhcmlhYmxlcyBvZiBbMTgwLDAwMCB0ZXJyb3Jpc3QgYXR0YWNrcyBiZXR3ZWVuIDE5NzAtMjAxN10oaHR0cHM6Ly93d3cua2FnZ2xlLmNvbS9TVEFSVC1VTUQvZ3RkL3ZlcnNpb24vMykuIE5vdGUgdGhhdCB3ZSBleGNsdWRlZCBzb21lIHZhcmlhYmxlcyBmb3IgY2xhcml0eSBhbmQgZGlkIHNvbWUgcHJlcHJvY2Vzc2luZyBvbiB0aGUgZGF0YS4KClN1cHBvc2UgeW91IGFyZSBnaXZlbiB0aGF0IGRhdGFzZXQgYW5kIHlvdSdyZSBhc2tlZCB0byBpbmZvcm0gcG9saWN5LW1ha2VycyBhYm91dCB0aGUgcmVsYXRpb25zaGlwIGJldHdlZW4gdGVycm9yaXN0IGF0dGFjayBkZXRhaWxzIGFuZCB0aGUgbnVtYmVyIG9mIHZpY3RpbXMuIFRoaXMgbm90ZWJvb2sgYXNrcyB5b3UgdG8gcGVyZm9ybSBpbml0aWFsIGRlc2NyaXB0aXZlIHN0YXRpc3RpY3MgYW5kIHRoZW4gYnVpbGQgYW5kIGV2YWx1YXRlIG1vZGVscyBvZiB0aGUgZGF0YS4KCllvdSBjYW4gbG9hZCB0aGUgZGF0YXNldCBhcyBmb2xsb3dzOgoKYGBge3J9CmxvYWQoJy4vZGF0YS9nbG9iYWxfdGVycm9yaXNtX2RhdGFiYXNlLlJEYXRhJykKYGBgCgpbVmFyaWFibGUgY29kZWJvb2tdKGh0dHBzOi8vd3d3LnN0YXJ0LnVtZC5lZHUvZ3RkL2Rvd25sb2Fkcy9Db2RlYm9vay5wZGYpOgoKLSBldmVudGlkOiB1bmlxdWUgZXZlbnQgaWRlbnRpZmllcgotIGl5ZWFyOiB5ZWFyCi0gaW1vbnRoOiBtb250aAotIGlkYXk6IGRheQotIG5wZXJwczogbnVtYmVyIG9mIHBlcnBldHJhdG9ycwotIHN1aWNpZGU6IHdoZXRoZXIgdGhlIGF0dGFjayB3YXMgYSBzdWljaWRlIGF0dGFjaAotIHJhbnNvbTogd2hldGhlciByYW5zb20gd2FzIGRlbWFuZGVkCi0gbmtpbGw6IG51bWJlciBvZiBraWxsZWQgdmljdGltcwotIG53b3VuZDogbnVtYmVyIG9mIHdvdW5kZWQgdmljdGltcwoKClRoaXMgY29tbWFuZCBsb2FkcyBhIGRhdGFmcmFtZSBjYWxsZWQgYGd0ZGAgaW4geW91ciBub3RlYm9vay4gWW91IGNhbiBxdWVyeSB0aGF0IGRhdGFmcmFtZSAgYXMgdXN1YWwuCgpIYXZlIGEgbG9vayBhdCB0aGUgbmFtZXMgKGNvbHVtbnMpIGluIHRoZSBkYXRhLmZyYW1lOgoKYGBge3J9Cm5hbWVzKGd0ZCkKYGBgCgouLi4gYW5kIHNob3cgdGhlIGZpcnN0IDEwIHJvd3M6CgpgYGB7cn0KaGVhZChndGQsIDEwKQpgYGAKCgpOb3RlIHRoYXQgd2UgZXhjbHVkZWQgdmFyaWFibGVzIChmb3IgY2xhcml0eSkgYW5kIG9ic2VydmF0aW9ucyAodG8gYXZvaWQgbWlzc2luZyB2YWx1ZXMpLCBzbyB0aGUgYWN0dWFsIGRpbWVuc2lvbnMgb2YgdGhpcyBkYXRhZnJhbWUgYXJlOgoKYGBge3J9CmRpbShndGQpCmBgYAoKIyMgVW5kZXJzdGFuZGluZyB0aGUgZGF0YQoKTGV0J3Mgc3RhcnQgd2l0aCB1bmRlcnN0YW5kaW5nIHRoZSBkYXRhIGJpdCBiZXR0ZXIuIFlvdSdkIHdhbnQgdG8gZG8gdGhpcyB0byBhdm9pZCBtb2RlbGxpbmcgcmVsYXRpb25zaGlwcyB0aGF0IGFyZSBub3QgbWVhbmluZ2Z1bC4KCiMjIyBUYXNrCgpMb29rIGF0IHRoZSBmcmVxdWVuY2llcyBvZiB0aGUgbnVtYmVyIG9mIHBlcnBldHJhdG9ycyBhbmQgc3Vic2V0IHRoZXNlIGZyZXF1ZW5jeSBjb3VudHMgYnkgdGhlIHN1aWNpZGUgYW5kIHJhbnNvbSB2YXJpYWJsZS4KCmBgYHtyfQojeW91ciBjb2RlIGNvbWVzIGhlcmUKYGBgCgojIyMgVGFzawoKSW4gd2hpY2ggeWVhciB3YXMgdGhlIG51bWJlciBvZiBwZXJwZXRyYXRvcnMgKG9uIGF2ZXJhZ2UpIHRoZSBoaWdoZXN0PwoKYGBge3J9CiN5b3VyIGNvZGUgY29tZXMgaGVyZQpgYGAKCiMjIyBUYXNrCgpXaGF0IGlzIHRoZSBtb3N0IGNvbW1vbiB2YWx1ZSBvZiB0aGUgbnVtYmVyIG9mIHBlcnNvbnMga2lsbGVkIGFuZCB3b3VuZGVkPwoKSGludDogYD9oaXN0YCBhbmQgYD90YWJsZWAKCmBgYHtyfQojeW91ciBjb2RlIGNvbWVzIGhlcmUKYGBgCiMjIyBUYXNrCgpEaXNwbGF5IHRoZSByZWxhdGlvbnNoaXAgYmV0d2VlbiB0aGUgbnVtYmVyIG9mIGtpbGxlZCB2aWN0aW1zIGFuZCB0aGUgbnVtYmVyIG9mIHdvdW5kZWQgdmljdGltcyBpbiBhIGZpZ3VyZS4KCldoYXQga2luZCBvZiBhIHJlbGF0aW9uc2hpcCBkbyB5b3UgZXhwZWN0PwoKSGludDogYD9wbG90YAoKYGBge3J9CiN5b3VyIGNvZGUgY29tZXMgaGVyZQpgYGAKCiMjIyBUYXNrCgpXaGF0IGlzIHRoZSBtZWFuIG51bWJlciBvZiBwZXJwZXRyYXRvcnMgd2hlbiB0aGUgYXR0YWNrIHdhcyBhIHN1aWNpZGUgYXR0YWNrPwoKSGludDogYD90YXBwbHlgCgpgYGB7cn0KI3lvdXIgY29kZSBjb21lcyBoZXJlCmBgYAoKIyMgU2ltcGxlIHJlZ3Jlc3Npb24KCldlJ2xsIHN0YXJ0IHdpdGggc2ltcGxlIHJlZ3Jlc3Npb24gKGkuZS4gb25lIHByZWRpY3RvciB2YXJpYWJsZSk6CgoKIyMjIFRhc2sKCkJ1aWxkIGEgc2ltcGxlIHJlZ3Jlc3Npb24gbW9kZWwgdGhhdCBtb2RlbHMgdGhlIG51bWJlciBvZiB3b3VuZGVkIHZpY3RpbXMgdGhyb3VnaCB0aGUgbnVtYmVyIG9mIHBlcnBldHJhdG9ycy4KCmBgYHtyfQojeW91ciBjb2RlIGNvbWVzIGhlcmUKYGBgCgojIyMgVGFzawoKSG93IHNhdGlzZmllZCBhcmUgeW91IHdpdGggeW91ciBtb2RlbD8gT25lIHdheSB0byBhc3Nlc3MgdGhlICJtb2RlbCBmaXQiIGlzIHRvIGNhbGN1bGF0ZSB0aGUgcm9vdCBtZWFuIHNxdWFyZSBlcnJvciAtIFJNU0UgLSAocmVzaWR1YWxzKS4gQ2FsY3VsYXRlIHRoYXQgbWV0cmljIGFuZCB0aGluayBhYm91dCB0aGUgbWVhbmluZyBvZiB0aGlzIGZpdCBpbmRleC4gV2hhdCBkb2VzIGl0IHRlbGwgeW91IGFuZCBob3cgc2F0aXNmaWVkIGFyZSB5b3Ugd2l0aCBpdD8KCmBgYHtyfQojeW91ciBjb2RlIGNvbWVzIGhlcmUKYGBgCgojIyMgVGFzawoKUGxvdCB0aGUgZml0dGVkIHZhbHVlcyAoaW4gZ3JlZW4pIGFuZCB0aGUgb2JzZXJ2ZWQgdmFsdWVzIChpbiBibHVlKSB0byBhc3Nlc3MgdGhlIG1vZGVsIGZpdCB2aXN1YWxseS4KCmBgYHtyfQojeW91ciBjb2RlIGNvbWVzIGhlcmUKYGBgCgoKIyMgTXVsdGlwbGUgcmVncmVzc2lvbgoKTm93IHlvdSBtaWdodCB3YW50IHRvIHVzZSBtdWx0aXBsZSB2YXJpYWJsZXMgaW4geW91ciBtb2RlbDoKCiMjIyBUYXNrCgpCdWlsZCBhIG11bHRpcGxlIHJlZ3Jlc3Npb24gbW9kZWwgdGhhdCBtb2RlbHMgdGhlIG51bWJlciBvZiBraWxsZWQgdmljdGltcyB0aHJvdWdoIHRoZSB2YXJpYWJsZXMgc3VpY2lkZSBhbmQgcmFuc29tLiBJbmNsdWRlIG9ubHkgdGhlIHR3byBtYWluIGVmZmVjdHMgKGFuZCBsZXQgdGhlIGludGVyY2VwdCBpbiB0aGUgbW9kZWwpLgoKYGBge3J9CiN5b3VyIGNvZGUgY29tZXMgaGVyZQpgYGAKCgojIyMgVGFzawoKSGF2ZSBhIGxvb2sgYXQgYSBwb3RlbnRpYWwgaW50ZXJhY3Rpb24gYmV0d2VlbiB0aGVzZSB0d28gcHJlZGljdG9yIHZhcmlhYmxlcy4gVXNlIHRoZSBgaW50ZXJhY3Rpb24ucGxvdGAgZnVuY3Rpb24gdG8gbG9vayBhdCB0aGUgam9pbnQgcmVsYXRpb25zaGlwIG9mIHRoZXNlIHR3byB2YXJpYWJsZXMgb24gdGhlIG51bWJlciBvZiBraWxsZWQgdmljdGltcy4KCmBgYHtyfQojeW91ciBjb2RlIGNvbWVzIGhlcmUKYGBgCgpXaGF0IGRvZXMgdGhpcyBncmFwaCB0ZWxsIHlvdT8gQ2FuIHlvdSBpZGVudGlmeSB0aGUgbWFpbiBlZmZlY3RzIGFuZCAocG90ZW50aWFsKSBpbnRlcmFjdGlvbj8KCiMjIyBUYXNrCgpOb3cgbG9vayBhdCB0aGUgaW50ZXJhY3Rpb24gaW4gYSBudW1lcmljYWwgbWFubmVyLgoKSGludDogYD90YXBwbHlgCgoKYGBge3J9CiN5b3VyIGNvZGUgY29tZXMgaGVyZQpgYGAKCgojIyMgVGFzawoKU3VwcG9zZSB5b3Ugd2FudCB0byBleHBhbmQgdGhlIG1vZGVsIGJ5IGFkZGluZyB0aGUgaW50ZXJhY3Rpb24gdGVybSB0byBpdC4gQnVpbGQgdGhhdCBtb2RlbC4KCgpgYGB7cn0KI3lvdXIgY29kZSBjb21lcyBoZXJlCmBgYAoKCiMjIyBUYXNrCgpCYXNlZCBvbiB0aGUgUk1TRSBvZiBlYWNoIG9mIHRoZSB0d28gbW9kZWxzIGFib3ZlICgyIG1haW4gZWZmZWN0cyB2cyAyIG1haW4gZWZmZWN0cyArIDEgaW50ZXJhY3Rpb24pLCB3aGljaCBvbmUgZG8geW91IHByZWZlcj8KCmBgYHtyfQojeW91ciBjb2RlIGNvbWVzIGhlcmUKYGBgCgojIyMgVGFzawoKSGF2ZSBhIGxvb2sgYXQgdGhlIGRpc3RyaWJ1dGlvbiBvZiB0aGUgYG5wZXJwc2AgYW5kIGBua2lsbGAgY29sdW1uLiBBcmUgdGhlcmUgc29tZSBwb3RlbnRpYWwgb3V0bGllcnMgaW4gdGhlcmU/CgpIaW50OiBgP3Bsb3RgCgpgYGB7cn0KI3lvdXIgY29kZSBjb21lcyBoZXJlCmBgYAoKIyMjIFRhc2sKClJlLXJ1biB0aGUgYmVzdCBmaXR0aW5nIHJlZ3Jlc3Npb24gbW9kZWwgYWdhaW4gYWZ0ZXIgZXhsdWRpbmcgdGhlIHBvdGVudGlhbCBvdXRsaWVycy4gVGhlIGRlY2lzaW9uIGZvciB0aGUgImJlc3QiIG1vZGVsIGNhbiBiZSBtYWRlIGJhc2VkIG9uIHRoZSBSTVNFOgoKYGBge3J9CiN5b3VyIGNvZGUgY29tZXMgaGVyZQpgYGAKCgojIyBNb2RlbCBzZWxlY3Rpb24KCk5vdyBzdXBwb3NlIHlvdSB3YW50IHRvIGxldCB0aGUgbW9kZWwgYnVpbGRpbmcgcHJvY2VzcyBiZSBkZWNpZGVkIGJ5IGEgc3RlcHdpc2UgbW9kZWwgc2VsZWN0aW9uIHByb2NlZHVyZS4KCiMjIyBUYXNrCgpCdWlsZCBhICJudWxsIG1vZGVsIiB0aGF0IG9ubHkgY29udGFpbnMgYW4gaW50ZXJjZXB0LgoKYGBge3J9CiN5b3VyIGNvZGUgY29tZXMgaGVyZQpgYGAKCgojIyMgVGFzawoKQnVpbGQgYSBmdWxsIG1vZGVsIGZvciB0aGUgbnVtYmVyIG9mIHdvdW5kZWQgdmljdGltcyBtb2RlbGVkIHRocm91Z2ggdGhlIG51bWJlciBvZiBwZXJwZXRyYXRvcnMsICJzdWljaWRlIiBhbmQgInJhbnNvbSIuCgpgYGB7cn0KI3lvdXIgY29kZSBjb21lcyBoZXJlCmBgYAoKCiMjIyBUYXNrCgpEZXRlcm1pbmUgdGhlIGJlc3QgZml0dGluZyBtb2RlbCBpbiBhIGJhY2t3YXJkIG1vZGVsIHNlbGVjdGlvbiBwcm9jZWR1cmUuCgpgYGB7cn0KI3lvdXIgY29kZSBjb21lcyBoZXJlCmBgYAoKIyMjIFRhc2sKClJ1biB0aGUgbW9kZWwgc2VsZWN0aW9uIGFnYWluIGJ1dCB0aGlzIHRpbWUgdXNpbmcgdGhlIGZvcndhcmQgbW9kZWwgc2VsZWN0aW9uIHByb2NlZHVyZS4KCmBgYHtyfQojeW91ciBjb2RlIGNvbWVzIGhlcmUKYGBgCgoKIyMjIFRhc2sKCkNvbXBhcmUgdGhlIFJNU0Ugb2YgdGhlIG51bGwgbW9kZWwsIHRoZSBmdWxsIG1vZGVsIGFuZCB0aGUgYmVzdCBmaXR0aW5nIG1vZGVsIChpZiBkaWZmZXJlbnQgZnJvbSB0aGUgZnVsbCBtb2RlbCkuCgpEaXNwbGF5IHRoZSByZXNpZHVhbHMgaW4gYSBncmFwaCB1c2luZyB0aGUgY29sb3VycyBncmVlbiwgYmxhY2ssIGFuZCBibHVlIGZvciB0aGUgbnVsbCBtb2RlbCwgb2JzZXJ2ZWQgdmFsdWVzLCBhbmQgdGhlIGZ1bGwgbW9kZWwsIHJlc3BlY3RpdmVseS4KCkhpbnQ6IGA/cG9pbnRzYAoKYGBge3J9CiN5b3VyIGNvZGUgY29tZXMgaGVyZQpgYGAKCkRpZCB5b3UgZXhwZWN0IHRoZSBvYnNlcnZhdGlvbnMgZm9yIHRoZSBudWxsIG1vZGVsPyBTZWUgd2hldGhlciB5b3UgY2FuIGRpc2NvdmVyIHRoZSByZWFzb24gZm9yIHRoYXQgcmVsYXRpb25zaGlwIGluIHRoZSBtb2RlbCBvdXRjb21lIChpLmUuIHRoZSBjb2VmZmljaWVudHMpLgoKIyMgRU5ECgotLS0=