Some background information

For this practical, we’ll be using open data from the following source: Wright, D. and London, K. (2011). Modern Regression Techniques Using R, London: Sage. (Chapter 8). This is available here.

This dataset contains information about a randomised controlled trial, in which primary school pupils were asked about the time they spent doing exercises (weekly). Then, their classes (within schools) were divided into 4 different groups and pupils within those classes were given different treatments:

  1. Control (no treatment at all)
  2. Leaflet (information about exercising)
  3. Leaflet+plan (the leaflet and a plan of exercise)
  4. Leaflet+plan+quiz (leaflet, plan and a quiz to do about exercising)

Pupils were asked again about the time they spent doing exercises (weekly) some time after the treatments were given.

The variable sqw2 is the outcome (time spent doing exercises after treatment - standardised). The variable sqw1 is the measure before the treatment. The variable wcond is the group or condition to which the class of the pupils was assigned

1. Define a working directory

You can use any directory in your computer. As in the example below:

setwd("C:/myfolder")

2. Read in data

This is from Wright and London’s book website companion.

mydata<-read.delim("https://us.sagepub.com/sites/default/files/upm-binaries/26934_exercise.dat")

Note: You can safely ignore the warning when reading in the data.

3. Descriptive statistics

R is very intuitive for certain basic functions. For example, you can obtain mean and standard deviation using basic functions

mean(mydata$sqw2, na.rm = TRUE)
## [1] 1.78919
mean(mydata$sqw1, na.rm = TRUE)
## [1] 1.680725
sd(mydata$sqw2, na.rm = TRUE)
## [1] 0.5633153
sd(mydata$sqw2, na.rm = TRUE)
## [1] 0.5633153

For other basic functions, you can consult the R reference card. On page 2, under the “Math” heading, you will some descriptive statistics of interest.

But then, since this is a cluster randomised controlled trial, you might want to know the means and standard deviations by class.

Task 1: Statistics by group

Create a summary dataset containing the mean and standard deviation of the outcome, grouped by class.

For this, you’ll need to load the dplyr package:

library(dplyr)

summary_data <- mydata %>%
  group_by(class) %>%
  summarise(mean_sqw2=mean(sqw2),
            stdev_sqw2=sd(sqw2))

You can browse your newly created dataset by clicking on it in the environment tab, or you can type summary_data to print it on the console.

summary_data 
## # A tibble: 22 x 3
##    class mean_sqw2 stdev_sqw2
##    <int>     <dbl>      <dbl>
##  1     1      1.54      0.685
##  2     2      2.05      0.547
##  3     3      1.79      0.542
##  4     4      1.90      0.470
##  5     5      1.78      0.583
##  6     6      2.07      0.389
##  7     7      1.68      0.602
##  8     8      1.80      0.546
##  9     9      1.53      0.572
## 10    10      1.87      0.537
## # ... with 12 more rows

If you want to export this dataset and use it in another software package, like Excel, you can use the write family of functions.

write.csv(summary_data, "summary_data.csv")

That will save a csv file (which you can open in Excel) in your working directory.

Task 2: Relationship between two variables

You can obtain the Pearson correlation by typing:

cor(mydata$sqw2, mydata$sqw1, use="comp") # "comp" stands for "complete observations"
## [1] 0.7648664

A basic bivariate plot (using ggplot) will include the name of the dataset, a y and x variable and the actual figure/symbol used to represent the values. It would have the following form:

First, we need to load the ggplot2 package

library(ggplot2)

plot1<-ggplot(mydata, 
              aes(y=sqw2, x=sqw1)) +
  geom_point()

plot1

We save plots to an object to be able to reproduce them more quickly. It is good practice, although not essential.

You can modify plot1 by including a “colour” subcommand inside the “geom_point()” command. This would be for example: geom_point(colour=“red”)

For more tips and example code, you can visit the following link: http://www.cookbook-r.com/Graphs/

ggplot2 works with “layers”. The basic plot above can be enhanced by adding a different geom (layer) that gives further instructions to ggplot. Below we’ll see some examples:

Task 3: Fitted line

Add a single fitted line to plot1

plot2<- ggplot(mydata, aes(y=sqw2, x=sqw1)) +
geom_point() +
geom_smooth(method="lm", aes(y=sqw2, x=sqw1), se=T)

plot2

Task 4: Fitted line per class

Add fitted lines per class. TO do this, you don’t need to add a different geom, but insert a group statement inside the aes subcommand.


plot3<- ggplot(mydata, aes(y=sqw2, x=sqw1, group=factor(class))) +
geom_point() +
geom_smooth(method="lm", aes(colour=factor(class)), se=F)

plot3

Bonus track

It’s easy to produce more elaborated graphics with ggplot2. Look at the example below. Try to identify the different layers that were added to plot3.


plot4<- ggplot(mydata, aes(y=sqw2, x=sqw1, by=factor(class))) +
geom_point(aes(colour=factor(class))) +
geom_smooth(method="lm", aes(y=sqw2, x=sqw1), se=F) +
facet_wrap(~class) + theme(legend.position = "none") +
ylab("After") + xlab("Before") +
theme_minimal() +
ggtitle("Relationship between exercise before and after treatment by class")

plot4

LS0tDQp0aXRsZTogIkRhdGEgbWFuaXB1bGF0aW9uIGFuZCBncmFwaGljcyBpbiBSIDxpbWcgc3JjPVwic2VtaW5hcnNfbG9nby5wbmdcIiB3aWR0aD1cIjIwMFwiIGhlaWdodD1cIjUwXCIgc3R5bGU9XCJmbG9hdDogcmlnaHQ7XCIvPiINCmF1dGhvcjogIlBhdHJpY2lvIFRyb25jb3NvIg0KZGF0ZTogIkxhdGVzdCB1cGRhdGU6IEZlYnJ1YXJ5IDIwMjAiDQpvdXRwdXQ6DQogIGh0bWxfZG9jdW1lbnQ6DQogICAgY29kZV9kb3dubG9hZDogeWVzDQogICAgaGlnaGxpZ2h0ZXI6IG51bGwNCiAgICB0aGVtZTogY29zbW8NCiAgICB0b2M6IHllcw0KICAgIHRvY19kZXB0aDogNA0KICAgIHRvY19mbG9hdDogeWVzDQogICAgZm9udHNpemU6IDEycHQNCiAgICBpbmNsdWRlczoNCiAgICAgICNpbl9oZWFkZXI6IGhlYWRlci5odG1sDQogICAgICBhZnRlcl9ib2R5OiBmb290ZXIuaHRtbA0KLS0tDQogIA0KYGBge3Igc2V0dXAsIGluY2x1ZGU9RkFMU0UsIHdhcm5pbmc9Rn0NCmtuaXRyOjpvcHRzX2NodW5rJHNldChlY2hvID0gVFJVRSkNCmBgYA0KDQojIFNvbWUgYmFja2dyb3VuZCBpbmZvcm1hdGlvbg0KDQpGb3IgdGhpcyBwcmFjdGljYWwsIHdl4oCZbGwgYmUgdXNpbmcgb3BlbiBkYXRhIGZyb20gdGhlIGZvbGxvd2luZyBzb3VyY2U6IFdyaWdodCwgRC4gYW5kIExvbmRvbiwgSy4gKDIwMTEpLiBNb2Rlcm4gUmVncmVzc2lvbiBUZWNobmlxdWVzIFVzaW5nIFIsIExvbmRvbjogIFNhZ2UuIChDaGFwdGVyIDgpLiBUaGlzIGlzIGF2YWlsYWJsZSBbKmhlcmUqXShodHRwOi8vZHguZG9pLm9yZy8xMC40MTM1Lzk3ODA4NTcwMjQ0OTcpLg0KDQpUaGlzIGRhdGFzZXQgY29udGFpbnMgaW5mb3JtYXRpb24gYWJvdXQgYSByYW5kb21pc2VkIGNvbnRyb2xsZWQgdHJpYWwsIGluIHdoaWNoIHByaW1hcnkgc2Nob29sIHB1cGlscyB3ZXJlIGFza2VkIGFib3V0IHRoZSB0aW1lIHRoZXkgc3BlbnQgZG9pbmcgZXhlcmNpc2VzICh3ZWVrbHkpLiBUaGVuLCB0aGVpciBjbGFzc2VzICh3aXRoaW4gc2Nob29scykgd2VyZSBkaXZpZGVkIGludG8gNCBkaWZmZXJlbnQgZ3JvdXBzIGFuZCBwdXBpbHMgd2l0aGluIHRob3NlIGNsYXNzZXMgd2VyZSBnaXZlbiBkaWZmZXJlbnQgdHJlYXRtZW50czoNCiAgDQoxKSBDb250cm9sIChubyB0cmVhdG1lbnQgYXQgYWxsKQ0KMikgTGVhZmxldCAoaW5mb3JtYXRpb24gYWJvdXQgZXhlcmNpc2luZykNCjMpIExlYWZsZXQrcGxhbiAodGhlIGxlYWZsZXQgYW5kIGEgcGxhbiBvZiBleGVyY2lzZSkNCjQpIExlYWZsZXQrcGxhbitxdWl6IChsZWFmbGV0LCBwbGFuIGFuZCBhIHF1aXogdG8gZG8gYWJvdXQgZXhlcmNpc2luZykNCg0KUHVwaWxzIHdlcmUgYXNrZWQgYWdhaW4gYWJvdXQgdGhlIHRpbWUgdGhleSBzcGVudCBkb2luZyBleGVyY2lzZXMgKHdlZWtseSkgc29tZSB0aW1lIGFmdGVyIHRoZSB0cmVhdG1lbnRzIHdlcmUgZ2l2ZW4uDQoNClRoZSB2YXJpYWJsZSBzcXcyIGlzIHRoZSBvdXRjb21lICh0aW1lIHNwZW50IGRvaW5nIGV4ZXJjaXNlcyBhZnRlciB0cmVhdG1lbnQgLSBzdGFuZGFyZGlzZWQpLg0KVGhlIHZhcmlhYmxlIHNxdzEgaXMgdGhlIG1lYXN1cmUgYmVmb3JlIHRoZSB0cmVhdG1lbnQuDQpUaGUgdmFyaWFibGUgd2NvbmQgaXMgdGhlIGdyb3VwIG9yIGNvbmRpdGlvbiB0byB3aGljaCB0aGUgY2xhc3Mgb2YgdGhlIHB1cGlscyB3YXMgYXNzaWduZWQNCg0KIyAxLiBEZWZpbmUgYSB3b3JraW5nIGRpcmVjdG9yeQ0KDQpZb3UgY2FuIHVzZSBhbnkgZGlyZWN0b3J5IGluIHlvdXIgY29tcHV0ZXIuIEFzIGluIHRoZSBleGFtcGxlIGJlbG93Og0KDQpgYGB7fQ0Kc2V0d2QoIkM6L215Zm9sZGVyIikNCmBgYA0KDQojIDIuIFJlYWQgaW4gZGF0YSANCg0KVGhpcyBpcyBmcm9tIFdyaWdodCBhbmQgTG9uZG9uJ3MgYm9vayB3ZWJzaXRlIGNvbXBhbmlvbi4gDQpgYGB7ciwgbWVzc2FnZT1GQUxTRSwgd2FybmluZz1GfQ0KbXlkYXRhPC1yZWFkLmRlbGltKCJodHRwczovL3VzLnNhZ2VwdWIuY29tL3NpdGVzL2RlZmF1bHQvZmlsZXMvdXBtLWJpbmFyaWVzLzI2OTM0X2V4ZXJjaXNlLmRhdCIpDQpgYGANCk5vdGU6IFlvdSBjYW4gc2FmZWx5IGlnbm9yZSB0aGUgd2FybmluZyB3aGVuIHJlYWRpbmcgaW4gdGhlIGRhdGEuDQoNCg0KIyAzLiBEZXNjcmlwdGl2ZSBzdGF0aXN0aWNzDQoNClIgaXMgdmVyeSBpbnR1aXRpdmUgZm9yIGNlcnRhaW4gYmFzaWMgZnVuY3Rpb25zLiBGb3IgZXhhbXBsZSwgeW91IGNhbiBvYnRhaW4gbWVhbiBhbmQgc3RhbmRhcmQgZGV2aWF0aW9uIHVzaW5nIGJhc2ljIGZ1bmN0aW9ucw0KDQpgYGB7cixtZXNzYWdlPUZBTFNFLCB3YXJuaW5nPUZ9DQptZWFuKG15ZGF0YSRzcXcyLCBuYS5ybSA9IFRSVUUpDQptZWFuKG15ZGF0YSRzcXcxLCBuYS5ybSA9IFRSVUUpDQpzZChteWRhdGEkc3F3MiwgbmEucm0gPSBUUlVFKQ0Kc2QobXlkYXRhJHNxdzIsIG5hLnJtID0gVFJVRSkNCmBgYA0KDQpGb3Igb3RoZXIgYmFzaWMgZnVuY3Rpb25zLCB5b3UgY2FuIGNvbnN1bHQgdGhlIFtSIHJlZmVyZW5jZSBjYXJkXShodHRwczovL2NyYW4uci1wcm9qZWN0Lm9yZy9kb2MvY29udHJpYi9TaG9ydC1yZWZjYXJkLnBkZikuIE9uIHBhZ2UgMiwgdW5kZXIgdGhlICJNYXRoIiBoZWFkaW5nLCB5b3Ugd2lsbCBzb21lIGRlc2NyaXB0aXZlIHN0YXRpc3RpY3Mgb2YgaW50ZXJlc3QuDQoNCkJ1dCB0aGVuLCBzaW5jZSB0aGlzIGlzIGEgY2x1c3RlciByYW5kb21pc2VkIGNvbnRyb2xsZWQgdHJpYWwsIHlvdSBtaWdodCB3YW50IHRvIGtub3cgdGhlIG1lYW5zIGFuZCBzdGFuZGFyZCBkZXZpYXRpb25zIGJ5IGNsYXNzLg0KDQojIyBUYXNrIDE6IFN0YXRpc3RpY3MgYnkgZ3JvdXANCg0KQ3JlYXRlIGEgc3VtbWFyeSBkYXRhc2V0IGNvbnRhaW5pbmcgdGhlIG1lYW4gYW5kIHN0YW5kYXJkIGRldmlhdGlvbiBvZiB0aGUgb3V0Y29tZSwgZ3JvdXBlZCBieSBjbGFzcy4NCg0KRm9yIHRoaXMsIHlvdSdsbCBuZWVkIHRvIGxvYWQgdGhlIGRwbHlyIHBhY2thZ2U6DQogIA0KYGBge3IsbWVzc2FnZT1GQUxTRSwgd2FybmluZz1GfQ0KbGlicmFyeShkcGx5cikNCg0Kc3VtbWFyeV9kYXRhIDwtIG15ZGF0YSAlPiUNCiAgZ3JvdXBfYnkoY2xhc3MpICU+JQ0KICBzdW1tYXJpc2UobWVhbl9zcXcyPW1lYW4oc3F3MiksDQogICAgICAgICAgICBzdGRldl9zcXcyPXNkKHNxdzIpKQ0KYGBgDQoNCllvdSBjYW4gYnJvd3NlIHlvdXIgbmV3bHkgY3JlYXRlZCBkYXRhc2V0IGJ5IGNsaWNraW5nIG9uIGl0IGluIHRoZSBlbnZpcm9ubWVudCB0YWIsIG9yIHlvdSBjYW4gdHlwZSBgc3VtbWFyeV9kYXRhYCB0byBwcmludCBpdCBvbiB0aGUgY29uc29sZS4NCg0KYGBge3IsbWVzc2FnZT1GQUxTRSwgd2FybmluZz1GfQ0Kc3VtbWFyeV9kYXRhIA0KYGBgDQoNCklmIHlvdSB3YW50IHRvIGV4cG9ydCB0aGlzIGRhdGFzZXQgYW5kIHVzZSBpdCBpbiBhbm90aGVyIHNvZnR3YXJlIHBhY2thZ2UsIGxpa2UgRXhjZWwsIHlvdSBjYW4gdXNlIHRoZSBgd3JpdGVgIGZhbWlseSBvZiBmdW5jdGlvbnMuDQoNCmBgYHtyLG1lc3NhZ2U9RkFMU0UsIHdhcm5pbmc9RiwgZXZhbD1GfQ0Kd3JpdGUuY3N2KHN1bW1hcnlfZGF0YSwgInN1bW1hcnlfZGF0YS5jc3YiKQ0KDQpgYGANCg0KVGhhdCB3aWxsIHNhdmUgYSBgY3N2YCBmaWxlICh3aGljaCB5b3UgY2FuIG9wZW4gaW4gRXhjZWwpIGluIHlvdXIgd29ya2luZyBkaXJlY3RvcnkuDQoNCiMjIFRhc2sgMjogUmVsYXRpb25zaGlwIGJldHdlZW4gdHdvIHZhcmlhYmxlcw0KICANCllvdSBjYW4gb2J0YWluIHRoZSBQZWFyc29uIGNvcnJlbGF0aW9uIGJ5IHR5cGluZzogDQogIA0KYGBge3IsbWVzc2FnZT1GQUxTRSwgd2FybmluZz1GfQ0KY29yKG15ZGF0YSRzcXcyLCBteWRhdGEkc3F3MSwgdXNlPSJjb21wIikgIyAiY29tcCIgc3RhbmRzIGZvciAiY29tcGxldGUgb2JzZXJ2YXRpb25zIg0KDQpgYGANCg0KQSBiYXNpYyBiaXZhcmlhdGUgcGxvdCAodXNpbmcgZ2dwbG90KSB3aWxsIGluY2x1ZGUgdGhlIG5hbWUgb2YgdGhlIGRhdGFzZXQsIGEgeSBhbmQgeCB2YXJpYWJsZSBhbmQgdGhlIGFjdHVhbCBmaWd1cmUvc3ltYm9sIHVzZWQgdG8gcmVwcmVzZW50IHRoZSB2YWx1ZXMuIEl0IHdvdWxkIGhhdmUgdGhlIGZvbGxvd2luZyBmb3JtOg0KICANCkZpcnN0LCB3ZSBuZWVkIHRvIGxvYWQgdGhlIGdncGxvdDIgcGFja2FnZQ0KDQpgYGB7cixtZXNzYWdlPUZBTFNFLCB3YXJuaW5nPUZ9DQpsaWJyYXJ5KGdncGxvdDIpDQoNCnBsb3QxPC1nZ3Bsb3QobXlkYXRhLCANCiAgICAgICAgICAgICAgYWVzKHk9c3F3MiwgeD1zcXcxKSkgKw0KICBnZW9tX3BvaW50KCkNCg0KcGxvdDENCmBgYA0KDQpXZSBzYXZlIHBsb3RzIHRvIGFuIG9iamVjdCB0byBiZSBhYmxlIHRvIHJlcHJvZHVjZSB0aGVtIG1vcmUgcXVpY2tseS4gSXQgaXMgZ29vZCBwcmFjdGljZSwgYWx0aG91Z2ggbm90IGVzc2VudGlhbC4NCg0KWW91IGNhbiBtb2RpZnkgcGxvdDEgYnkgaW5jbHVkaW5nIGEgImNvbG91ciIgc3ViY29tbWFuZCBpbnNpZGUgdGhlICJnZW9tX3BvaW50KCkiIGNvbW1hbmQuIFRoaXMgd291bGQgYmUgZm9yIGV4YW1wbGU6IGdlb21fcG9pbnQoY29sb3VyPSJyZWQiKQ0KDQpGb3IgbW9yZSB0aXBzIGFuZCBleGFtcGxlIGNvZGUsIHlvdSBjYW4gdmlzaXQgdGhlIGZvbGxvd2luZyBsaW5rOiBodHRwOi8vd3d3LmNvb2tib29rLXIuY29tL0dyYXBocy8NCg0KYGdncGxvdDJgIHdvcmtzIHdpdGggImxheWVycyIuIFRoZSBiYXNpYyBwbG90IGFib3ZlIGNhbiBiZSBlbmhhbmNlZCBieSBhZGRpbmcgYSBkaWZmZXJlbnQgYGdlb21gIChsYXllcikgdGhhdCBnaXZlcyBmdXJ0aGVyIGluc3RydWN0aW9ucyB0byBgZ2dwbG90YC4gQmVsb3cgd2UnbGwgc2VlIHNvbWUgZXhhbXBsZXM6DQoNCg0KIyMgVGFzayAzOiBGaXR0ZWQgbGluZQ0KDQpBZGQgYSBzaW5nbGUgZml0dGVkIGxpbmUgdG8gcGxvdDENCg0KDQpgYGB7ciwgd2FybmluZz1GfQ0KDQpwbG90MjwtIGdncGxvdChteWRhdGEsIGFlcyh5PXNxdzIsIHg9c3F3MSkpICsNCmdlb21fcG9pbnQoKSArDQpnZW9tX3Ntb290aChtZXRob2Q9ImxtIiwgYWVzKHk9c3F3MiwgeD1zcXcxKSwgc2U9VCkNCg0KcGxvdDINCmBgYA0KDQoNCg0KIyMgVGFzayA0OiBGaXR0ZWQgbGluZSBwZXIgY2xhc3MNCg0KQWRkIGZpdHRlZCBsaW5lcyBwZXIgY2xhc3MuIFRPIGRvIHRoaXMsIHlvdSBkb24ndCBuZWVkIHRvIGFkZCBhIGRpZmZlcmVudCBgZ2VvbWAsIGJ1dCBpbnNlcnQgYSBgZ3JvdXBgIHN0YXRlbWVudCBpbnNpZGUgdGhlIGBhZXNgIHN1YmNvbW1hbmQuDQoNCjxicj4NCmBgYHtyLCB3YXJuaW5nPUZ9DQpwbG90MzwtIGdncGxvdChteWRhdGEsIGFlcyh5PXNxdzIsIHg9c3F3MSwgZ3JvdXA9ZmFjdG9yKGNsYXNzKSkpICsNCmdlb21fcG9pbnQoKSArDQpnZW9tX3Ntb290aChtZXRob2Q9ImxtIiwgYWVzKGNvbG91cj1mYWN0b3IoY2xhc3MpKSwgc2U9RikNCg0KcGxvdDMNCmBgYA0KDQoNCg0KIyMgQm9udXMgdHJhY2sNCg0KSXQncyBlYXN5IHRvIHByb2R1Y2UgbW9yZSBlbGFib3JhdGVkIGdyYXBoaWNzIHdpdGggZ2dwbG90Mi4gTG9vayBhdCB0aGUgZXhhbXBsZSBiZWxvdy4gVHJ5IHRvIGlkZW50aWZ5IHRoZSBkaWZmZXJlbnQgbGF5ZXJzIHRoYXQgd2VyZSBhZGRlZCB0byBwbG90My4NCg0KPGJyPg0KYGBge3IsIHdhcm5pbmc9Rn0NCnBsb3Q0PC0gZ2dwbG90KG15ZGF0YSwgYWVzKHk9c3F3MiwgeD1zcXcxLCBieT1mYWN0b3IoY2xhc3MpKSkgKw0KZ2VvbV9wb2ludChhZXMoY29sb3VyPWZhY3RvcihjbGFzcykpKSArDQpnZW9tX3Ntb290aChtZXRob2Q9ImxtIiwgYWVzKHk9c3F3MiwgeD1zcXcxKSwgc2U9RikgKw0KZmFjZXRfd3JhcCh+Y2xhc3MpICsgdGhlbWUobGVnZW5kLnBvc2l0aW9uID0gIm5vbmUiKSArDQp5bGFiKCJBZnRlciIpICsgeGxhYigiQmVmb3JlIikgKw0KdGhlbWVfbWluaW1hbCgpICsNCmdndGl0bGUoIlJlbGF0aW9uc2hpcCBiZXR3ZWVuIGV4ZXJjaXNlIGJlZm9yZSBhbmQgYWZ0ZXIgdHJlYXRtZW50IGJ5IGNsYXNzIikNCg0KcGxvdDQNCmBgYA==
 

Patricio Troncoso

patricio.troncoso@manchester.ac.uk