Aims of this notebook

Requirements

We assume that you have:

If you struggle with basics of R, you may also find this online book useful: https://bookdown.org/ndphillips/YaRrr/


Non-parametric data analysis

Task 1: Correlation

Read the csv file data1_week5.csv from ./data. The csv file contains three variables of 1000 observations of areas surveyed and rated according to:

  • their rank in deprivation on a scale from 1 to 10 (1 highest, 10 lowest)
  • their rank in public disorder on a scale from 1 to 10 (1 highest, 10 lowest)
  • their rank in social media participation on a scale from 1 to 5 (1 highest, 5 lowest)

You are interested in testing whether these data correlate. First, test whether the assumption of normally distributed data is met and the choose the appropriate test to ascertain whether or not there are correlations between these three variables.

#your code

Task 2: Comparing two groups

Load the dataset (from the csv file data2_week5.csv) that gives you 1000 observations of male adolescents and the number of brothers they have and whether or not they had contact with the police in the past 5 years.

Assess whether those who had contact with the police differ from those who have not in the number of brothers they have. Again, make sure you test whether the assumptions are met and use an appropriate test.

#your code

Task 3: Comparing 2 (or more) groups

Read the next csv file (data3_week5.csv). This file shows the number of casualties in traffic accidents on 600 high-risk street crossings. Each crossing is either a zebra-crossing, a traffic-light crossing, or an inofficial crossing that has none of the aforementioned.

Assess whether the factor “Type of crossing” has an effect on the number of casualties.

#your code

Discrete data

Task 4: 2-by-2 tables

Let’s start with simple 2-by-2 tables.

Run the command below to create a 2-by-2 table that shows you the number of people that were refused boarding an airplane dependin whether they were drunk or not.

flight_data = array(c(100, 50, 150, 500), dim=c(2,2))
dimnames(flight_data) = list( c('Drunk', 'Not drunk')
                       , c('Denied boarding', 'Allowed boarding'))
flight_data = addmargins(flight_data, c(1,2))
flight_data

Assess whether there is an association between being drunk and being allowed to board the airplane.

Perform a Chi-Square test manually by calculating the expected cell counts and deriving the ChiSquare value.

#your code

Think for a moment: what do the expected values mean?

Now validate your results against the automatically calculated values.

#your code

Think for a moment: what does the test statistic and significance level mean?

Now follow-up on this finding if you deem this necessary. Is there an association? If so, what is it driven by?

#your code (if needed)

Think for a moment: what does each z-score mean?

Task 5: R-by-C tables

Run the code below to build another array in the simple row-by-column format but with more levels for each. The below data array expresses the count of flight delays by the season for a small airport.

flight_delays = array(c(300, 350, 200, 150, 600, 300, 200, 100, 400, 100, 50, 50, 300, 250, 200, 200), dim=c(4,4))
dimnames(flight_delays) = list( c('On time', 'Minor delay', 'Severe delay', 'Cancelled')
                       , c('Spring', 'Summer', 'Fall', 'Winter'))
flight_delays = addmargins(flight_delays, c(1,2))
flight_delays

Use the appropriate statistical test and follow-up on the findings if necessary.

#your code

Which cells contain more counts than expected, which ones contain fewer counts than expected, and for which cells does this not differ significantly?

Task 6: Multidimensional arrays

Run the code below to create a three-dimensional array that provides you with count data on whether someone was subjected to additional checks (e.g. an interview) upon entering the UK at an airport border control, which area they came from (EU vs nonEU) and whether the passenger was additionally checked by a customs office.

airport_data = array(c(1000, 500, 4000, 5000, 8000, 1000, 6000, 4000), dim=c(2,2,2))
dimnames(airport_data) = list('additional_checks' = c('yes', 'no')
                       , 'zone' = c('EU', 'nonEU')
                       , 'customs_check' = c('yes', 'no')
                       )
airport_data

To test whether there are associations between the factors “additional_checks”, “customes_check”, and “zone”, we’d want to use a loglinear model. Before we can use that model, we need to transform the representation of the data.

airport_data_long = as.data.frame(as.table(airport_data))
airport_data_long

Now the modelling process can start.

First, build and assess the independence model:

#your code

Take a look at the actual vs fitted values:

#your code

We’d like to use a better model for the data. The ideal (= perfect) model is of course the saturated (full) model. Build this model and show that the fitted values are identical to the observed ones.

#your code

The task is to find a simpler model that still is adequate. Next, build a model that in addition to the independence model contains two-way dependencies (similar to interactions) and test whether you can reject the H0 of model adequacy.

#your code

What do you conclude?

Finally, use the automated model select with the step function.

#your code

What do you conclude?

END


LS0tCnRpdGxlOiAiTm9uLXBhcmFtZXRyaWMgZGF0YSAmIGRpc2NyZXRlIGRhdGEiCmF1dGhvcjogQiBLbGVpbmJlcmcKc3VidGl0bGU6IERlcHQgb2YgU2VjdXJpdHkgYW5kIENyaW1lIFNjaWVuY2UsIFVDTApvdXRwdXQ6IGh0bWxfbm90ZWJvb2sKLS0tCgojIyBBaW1zIG9mIHRoaXMgbm90ZWJvb2sKCi0gbm9uLXBhcmFtZXRyaWMgZGF0YQogICAgLSBjb3JyZWxhdGlvbgogICAgLSBtZWFuIGNvbXBhcmlzb24gZm9yIHR3byBncm91cHMKICAgIC0gbWVhbiBjb21wYXJpc29uIGZvciB0d28gYW5kIG1vcmUgZ3JvdXBzCi0gZGlzY3JldGUgZGF0YQogICAgLSAyLWJ5LTIgdGFibGVzCiAgICAtIFItYnktQyB0YWJsZXMKICAgIC0gTXVsdGlkaW1lbnNpb25hbCBhcnJheXMKCiAgICAKIyMgUmVxdWlyZW1lbnRzCgpXZSBhc3N1bWUgdGhhdCB5b3UgaGF2ZToKCi0gcmVhZCB0aGUgcmVxdWlyZWQgbGl0ZXJhdHVyZSBmb3Igd2Vla3MgMS01Ci0gcmV2aXNlZCB0aGUgbGVjdHVyZXMKLSBjb21wbGV0ZWQgdGhlIGludHJvZHVjdG9yeSBSIHR1dG9yaWFscyAoMTIgc3RlcHMgJiBIb3cgdG8gc29sdmUgcHJvYmxlbXMpIGFzIHdlbGwgYXMgdGhlIHR1dG9yaWFsIGZyb20gd2VlayAyCi0gY29tcGxldGVkIHRoZSBob21ld29yayBpbiB0aGlzIG1vZHVsZSBzbyBmYXIKLSByZXBsaWNhdGVkIHRoZSBjb2RlIGZyb20gdGhlIGxlY3R1cmVzIChpZiBjb25jZXB0cy9SIGltcGxlbWVudGF0aW9uIGlzIHVuY2xlYXIpCgpJZiB5b3Ugc3RydWdnbGUgd2l0aCBiYXNpY3Mgb2YgUiwgeW91IG1heSBhbHNvIGZpbmQgdGhpcyBvbmxpbmUgYm9vayB1c2VmdWw6IFtodHRwczovL2Jvb2tkb3duLm9yZy9uZHBoaWxsaXBzL1lhUnJyL10oaHR0cHM6Ly9ib29rZG93bi5vcmcvbmRwaGlsbGlwcy9ZYVJyci8pCgotLS0KCiMjIE5vbi1wYXJhbWV0cmljIGRhdGEgYW5hbHlzaXMKCiMjIyBUYXNrIDE6IENvcnJlbGF0aW9uCgpSZWFkIHRoZSBjc3YgZmlsZSBgZGF0YTFfd2VlazUuY3N2YCBmcm9tIGAuL2RhdGFgLiBUaGUgY3N2IGZpbGUgY29udGFpbnMgdGhyZWUgdmFyaWFibGVzIG9mIDEwMDAgb2JzZXJ2YXRpb25zIG9mIGFyZWFzIHN1cnZleWVkIGFuZCByYXRlZCBhY2NvcmRpbmcgdG86CgotIHRoZWlyIHJhbmsgaW4gZGVwcml2YXRpb24gb24gYSBzY2FsZSBmcm9tIDEgdG8gMTAgKDEgaGlnaGVzdCwgMTAgbG93ZXN0KQotIHRoZWlyIHJhbmsgaW4gcHVibGljIGRpc29yZGVyIG9uIGEgc2NhbGUgZnJvbSAxIHRvIDEwICgxIGhpZ2hlc3QsIDEwIGxvd2VzdCkKLSB0aGVpciByYW5rIGluIHNvY2lhbCBtZWRpYSBwYXJ0aWNpcGF0aW9uIG9uIGEgc2NhbGUgZnJvbSAxIHRvIDUgKDEgaGlnaGVzdCwgNSBsb3dlc3QpCgpZb3UgYXJlIGludGVyZXN0ZWQgaW4gdGVzdGluZyB3aGV0aGVyIHRoZXNlIGRhdGEgY29ycmVsYXRlLiBGaXJzdCwgdGVzdCB3aGV0aGVyIHRoZSBhc3N1bXB0aW9uIG9mIG5vcm1hbGx5IGRpc3RyaWJ1dGVkIGRhdGEgaXMgbWV0IGFuZCB0aGUgY2hvb3NlIHRoZSBhcHByb3ByaWF0ZSB0ZXN0IHRvIGFzY2VydGFpbiB3aGV0aGVyIG9yIG5vdCB0aGVyZSBhcmUgY29ycmVsYXRpb25zIGJldHdlZW4gdGhlc2UgdGhyZWUgdmFyaWFibGVzLgoKYGBge3J9CiN5b3VyIGNvZGUKCmBgYAoKCiMjIyBUYXNrIDI6IENvbXBhcmluZyB0d28gZ3JvdXBzCgpMb2FkIHRoZSBkYXRhc2V0IChmcm9tIHRoZSBjc3YgZmlsZSBgZGF0YTJfd2VlazUuY3N2YCkgdGhhdCBnaXZlcyB5b3UgMTAwMCBvYnNlcnZhdGlvbnMgb2YgbWFsZSBhZG9sZXNjZW50cyBhbmQgdGhlIG51bWJlciBvZiBicm90aGVycyB0aGV5IGhhdmUgYW5kIHdoZXRoZXIgb3Igbm90IHRoZXkgaGFkIGNvbnRhY3Qgd2l0aCB0aGUgcG9saWNlIGluIHRoZSBwYXN0IDUgeWVhcnMuCgpBc3Nlc3Mgd2hldGhlciB0aG9zZSB3aG8gaGFkIGNvbnRhY3Qgd2l0aCB0aGUgcG9saWNlIGRpZmZlciBmcm9tIHRob3NlIHdobyBoYXZlIG5vdCBpbiB0aGUgbnVtYmVyIG9mIGJyb3RoZXJzIHRoZXkgaGF2ZS4gQWdhaW4sIG1ha2Ugc3VyZSB5b3UgdGVzdCB3aGV0aGVyIHRoZSBhc3N1bXB0aW9ucyBhcmUgbWV0IGFuZCB1c2UgYW4gYXBwcm9wcmlhdGUgdGVzdC4KCmBgYHtyfQojeW91ciBjb2RlCgpgYGAKCgojIyMgVGFzayAzOiBDb21wYXJpbmcgMiAob3IgbW9yZSkgZ3JvdXBzCgpSZWFkIHRoZSBuZXh0IGNzdiBmaWxlIChgZGF0YTNfd2VlazUuY3N2YCkuIFRoaXMgZmlsZSBzaG93cyB0aGUgbnVtYmVyIG9mIGNhc3VhbHRpZXMgaW4gdHJhZmZpYyBhY2NpZGVudHMgb24gNjAwIGhpZ2gtcmlzayBzdHJlZXQgY3Jvc3NpbmdzLiBFYWNoIGNyb3NzaW5nIGlzIGVpdGhlciBhIHplYnJhLWNyb3NzaW5nLCBhIHRyYWZmaWMtbGlnaHQgY3Jvc3NpbmcsIG9yIGFuIGlub2ZmaWNpYWwgY3Jvc3NpbmcgdGhhdCBoYXMgbm9uZSBvZiB0aGUgYWZvcmVtZW50aW9uZWQuCgpBc3Nlc3Mgd2hldGhlciB0aGUgZmFjdG9yICJUeXBlIG9mIGNyb3NzaW5nIiBoYXMgYW4gZWZmZWN0IG9uIHRoZSBudW1iZXIgb2YgY2FzdWFsdGllcy4gCgpgYGB7cn0KI3lvdXIgY29kZQoKYGBgCgojIyBEaXNjcmV0ZSBkYXRhCgojIyMgVGFzayA0OiAyLWJ5LTIgdGFibGVzCgpMZXQncyBzdGFydCB3aXRoIHNpbXBsZSAyLWJ5LTIgdGFibGVzLgoKUnVuIHRoZSBjb21tYW5kIGJlbG93IHRvIGNyZWF0ZSBhIDItYnktMiB0YWJsZSB0aGF0IHNob3dzIHlvdSB0aGUgbnVtYmVyIG9mIHBlb3BsZSB0aGF0IHdlcmUgcmVmdXNlZCBib2FyZGluZyBhbiBhaXJwbGFuZSBkZXBlbmRpbiB3aGV0aGVyIHRoZXkgd2VyZSBkcnVuayBvciBub3QuCgpgYGB7cn0KZmxpZ2h0X2RhdGEgPSBhcnJheShjKDEwMCwgNTAsIDE1MCwgNTAwKSwgZGltPWMoMiwyKSkKZGltbmFtZXMoZmxpZ2h0X2RhdGEpID0gbGlzdCggYygnRHJ1bmsnLCAnTm90IGRydW5rJykKICAgICAgICAgICAgICAgICAgICAgICAsIGMoJ0RlbmllZCBib2FyZGluZycsICdBbGxvd2VkIGJvYXJkaW5nJykpCmZsaWdodF9kYXRhID0gYWRkbWFyZ2lucyhmbGlnaHRfZGF0YSwgYygxLDIpKQpmbGlnaHRfZGF0YQpgYGAKCkFzc2VzcyB3aGV0aGVyIHRoZXJlIGlzIGFuIGFzc29jaWF0aW9uIGJldHdlZW4gYmVpbmcgZHJ1bmsgYW5kIGJlaW5nIGFsbG93ZWQgdG8gYm9hcmQgdGhlIGFpcnBsYW5lLgoKUGVyZm9ybSBhIENoaS1TcXVhcmUgdGVzdCBtYW51YWxseSBieSBjYWxjdWxhdGluZyB0aGUgZXhwZWN0ZWQgY2VsbCBjb3VudHMgYW5kIGRlcml2aW5nIHRoZSBDaGlTcXVhcmUgdmFsdWUuCgpgYGB7cn0KI3lvdXIgY29kZQoKCmBgYAoKX1RoaW5rIGZvciBhIG1vbWVudDogd2hhdCBkbyB0aGUgZXhwZWN0ZWQgdmFsdWVzIG1lYW4/XwoKTm93IHZhbGlkYXRlIHlvdXIgcmVzdWx0cyBhZ2FpbnN0IHRoZSBhdXRvbWF0aWNhbGx5IGNhbGN1bGF0ZWQgdmFsdWVzLgoKYGBge3J9CiN5b3VyIGNvZGUKCgpgYGAKCl9UaGluayBmb3IgYSBtb21lbnQ6IHdoYXQgZG9lcyB0aGUgdGVzdCBzdGF0aXN0aWMgYW5kIHNpZ25pZmljYW5jZSBsZXZlbCBtZWFuP18KCk5vdyBmb2xsb3ctdXAgb24gdGhpcyBmaW5kaW5nIGlmIHlvdSBkZWVtIHRoaXMgbmVjZXNzYXJ5LiBJcyB0aGVyZSBhbiBhc3NvY2lhdGlvbj8gSWYgc28sIHdoYXQgaXMgaXQgZHJpdmVuIGJ5PwoKYGBge3J9CiN5b3VyIGNvZGUgKGlmIG5lZWRlZCkKCmBgYAoKX1RoaW5rIGZvciBhIG1vbWVudDogd2hhdCBkb2VzIGVhY2ggei1zY29yZSBtZWFuP18KCiMjIyBUYXNrIDU6IFItYnktQyB0YWJsZXMKClJ1biB0aGUgY29kZSBiZWxvdyB0byBidWlsZCBhbm90aGVyIGFycmF5IGluIHRoZSBzaW1wbGUgcm93LWJ5LWNvbHVtbiBmb3JtYXQgYnV0IHdpdGggbW9yZSBsZXZlbHMgZm9yIGVhY2guIFRoZSBiZWxvdyBkYXRhIGFycmF5IGV4cHJlc3NlcyB0aGUgY291bnQgb2YgZmxpZ2h0IGRlbGF5cyBieSB0aGUgc2Vhc29uIGZvciBhIHNtYWxsIGFpcnBvcnQuCgpgYGB7cn0KZmxpZ2h0X2RlbGF5cyA9IGFycmF5KGMoMzAwLCAzNTAsIDIwMCwgMTUwLCA2MDAsIDMwMCwgMjAwLCAxMDAsIDQwMCwgMTAwLCA1MCwgNTAsIDMwMCwgMjUwLCAyMDAsIDIwMCksIGRpbT1jKDQsNCkpCmRpbW5hbWVzKGZsaWdodF9kZWxheXMpID0gbGlzdCggYygnT24gdGltZScsICdNaW5vciBkZWxheScsICdTZXZlcmUgZGVsYXknLCAnQ2FuY2VsbGVkJykKICAgICAgICAgICAgICAgICAgICAgICAsIGMoJ1NwcmluZycsICdTdW1tZXInLCAnRmFsbCcsICdXaW50ZXInKSkKZmxpZ2h0X2RlbGF5cyA9IGFkZG1hcmdpbnMoZmxpZ2h0X2RlbGF5cywgYygxLDIpKQpmbGlnaHRfZGVsYXlzCmBgYAoKVXNlIHRoZSBhcHByb3ByaWF0ZSBzdGF0aXN0aWNhbCB0ZXN0IGFuZCBmb2xsb3ctdXAgb24gdGhlIGZpbmRpbmdzIGlmIG5lY2Vzc2FyeS4KCmBgYHtyfQojeW91ciBjb2RlCgoKYGBgCgpfV2hpY2ggY2VsbHMgY29udGFpbiBtb3JlIGNvdW50cyB0aGFuIGV4cGVjdGVkLCB3aGljaCBvbmVzIGNvbnRhaW4gZmV3ZXIgY291bnRzIHRoYW4gZXhwZWN0ZWQsIGFuZCBmb3Igd2hpY2ggY2VsbHMgZG9lcyB0aGlzIG5vdCBkaWZmZXIgc2lnbmlmaWNhbnRseT9fCgoKIyMjIFRhc2sgNjogTXVsdGlkaW1lbnNpb25hbCBhcnJheXMKClJ1biB0aGUgY29kZSBiZWxvdyB0byBjcmVhdGUgYSB0aHJlZS1kaW1lbnNpb25hbCBhcnJheSB0aGF0IHByb3ZpZGVzIHlvdSB3aXRoIGNvdW50IGRhdGEgb24gd2hldGhlciBzb21lb25lIHdhcyBzdWJqZWN0ZWQgdG8gYWRkaXRpb25hbCBjaGVja3MgKGUuZy4gYW4gaW50ZXJ2aWV3KSB1cG9uIGVudGVyaW5nIHRoZSBVSyBhdCBhbiBhaXJwb3J0IGJvcmRlciBjb250cm9sLCB3aGljaCBhcmVhIHRoZXkgY2FtZSBmcm9tIChFVSB2cyBub25FVSkgYW5kIHdoZXRoZXIgdGhlIHBhc3NlbmdlciB3YXMgYWRkaXRpb25hbGx5IGNoZWNrZWQgYnkgYSBjdXN0b21zIG9mZmljZS4KCmBgYHtyfQphaXJwb3J0X2RhdGEgPSBhcnJheShjKDEwMDAsIDUwMCwgNDAwMCwgNTAwMCwgODAwMCwgMTAwMCwgNjAwMCwgNDAwMCksIGRpbT1jKDIsMiwyKSkKZGltbmFtZXMoYWlycG9ydF9kYXRhKSA9IGxpc3QoJ2FkZGl0aW9uYWxfY2hlY2tzJyA9IGMoJ3llcycsICdubycpCiAgICAgICAgICAgICAgICAgICAgICAgLCAnem9uZScgPSBjKCdFVScsICdub25FVScpCiAgICAgICAgICAgICAgICAgICAgICAgLCAnY3VzdG9tc19jaGVjaycgPSBjKCd5ZXMnLCAnbm8nKQogICAgICAgICAgICAgICAgICAgICAgICkKYWlycG9ydF9kYXRhCmBgYAoKVG8gdGVzdCB3aGV0aGVyIHRoZXJlIGFyZSBhc3NvY2lhdGlvbnMgYmV0d2VlbiB0aGUgZmFjdG9ycyAiYWRkaXRpb25hbF9jaGVja3MiLCAiY3VzdG9tZXNfY2hlY2siLCBhbmQgInpvbmUiLCB3ZSdkIHdhbnQgdG8gdXNlIGEgbG9nbGluZWFyIG1vZGVsLiBCZWZvcmUgd2UgY2FuIHVzZSB0aGF0IG1vZGVsLCB3ZSBuZWVkIHRvIHRyYW5zZm9ybSB0aGUgcmVwcmVzZW50YXRpb24gb2YgdGhlIGRhdGEuCgpgYGB7cn0KYWlycG9ydF9kYXRhX2xvbmcgPSBhcy5kYXRhLmZyYW1lKGFzLnRhYmxlKGFpcnBvcnRfZGF0YSkpCmFpcnBvcnRfZGF0YV9sb25nCmBgYAoKTm93IHRoZSBtb2RlbGxpbmcgcHJvY2VzcyBjYW4gc3RhcnQuCgpGaXJzdCwgYnVpbGQgYW5kIGFzc2VzcyB0aGUgaW5kZXBlbmRlbmNlIG1vZGVsOgoKYGBge3J9CiN5b3VyIGNvZGUKCgpgYGAKClRha2UgYSBsb29rIGF0IHRoZSBhY3R1YWwgdnMgZml0dGVkIHZhbHVlczoKCmBgYHtyfQojeW91ciBjb2RlCgpgYGAKCgpXZSdkIGxpa2UgdG8gdXNlIGEgYmV0dGVyIG1vZGVsIGZvciB0aGUgZGF0YS4gVGhlIGlkZWFsICg9IHBlcmZlY3QpIG1vZGVsIGlzIG9mIGNvdXJzZSB0aGUgc2F0dXJhdGVkIChmdWxsKSBtb2RlbC4gQnVpbGQgdGhpcyBtb2RlbCBhbmQgc2hvdyB0aGF0IHRoZSBmaXR0ZWQgdmFsdWVzIGFyZSBpZGVudGljYWwgdG8gdGhlIG9ic2VydmVkIG9uZXMuCgpgYGB7cn0KI3lvdXIgY29kZQoKCmBgYAoKVGhlIHRhc2sgaXMgdG8gZmluZCBhIHNpbXBsZXIgbW9kZWwgdGhhdCBzdGlsbCBpcyBhZGVxdWF0ZS4gTmV4dCwgYnVpbGQgYSBtb2RlbCB0aGF0IGluIGFkZGl0aW9uIHRvIHRoZSBpbmRlcGVuZGVuY2UgbW9kZWwgY29udGFpbnMgdHdvLXdheSBkZXBlbmRlbmNpZXMgKHNpbWlsYXIgdG8gaW50ZXJhY3Rpb25zKSBhbmQgdGVzdCB3aGV0aGVyIHlvdSBjYW4gcmVqZWN0IHRoZSBIMCBvZiBtb2RlbCBhZGVxdWFjeS4KCmBgYHtyfQojeW91ciBjb2RlCgoKYGBgCgpfV2hhdCBkbyB5b3UgY29uY2x1ZGU/XwoKRmluYWxseSwgdXNlIHRoZSBhdXRvbWF0ZWQgbW9kZWwgc2VsZWN0IHdpdGggdGhlIGBzdGVwYCBmdW5jdGlvbi4KCmBgYHtyfQojeW91ciBjb2RlCgoKYGBgCgpfV2hhdCBkbyB5b3UgY29uY2x1ZGU/XwoKIyMgRU5ECgotLS0=