This document describes your second assignment for Statistical Computing and Empirical Methods (Unit EMATM0061) on the MSc in Data Science. Before starting the assignment it is recommend that you first watch video lectures 3, 4 and 5.
You are encouraged to discuss these questions with your colleagues.
Begin by creating an Rmarkdown document with html output. You are not expected to hand in this piece of work, but it is a good idea to get used to using Rmarkdown.
You will need to install the Tidyverse set of packages if you have not already done so:
install.packages("tidyverse")
Next load the Tidyverse library which is a superset of the ggplot2 library which we shall use for this assignment.
library(tidyverse)
For the purpose of this assignment we shall use the Hawks data set from the Stat2Data package. This is a brilliant data set containing information about 908 Hawks collected by students and faculty at Cornell College in Mount Vernon, Iowa. First install the “Stat2Data” library:
install.packages("Stat2Data")
We can then load the data set as follows.
This will load a data frame called “Hawks” into your environment. We will use a slightly smaller dataframe which we construct as follows:
hawksSmall<-drop_na(select(Hawks,Age,Day,Month,Year,CaptureTime,Species,Wing,Weight,Tail))
1 Visualisation
We begin by studying some of the concepts from Lecture 3. 1
library(Stat2Data)
data("Hawks")
1.1 Types of variables
Check how many rows and columns hawksSmall using the dim() function. Use the head() function to display the top 5 rows of the hawksSmall data frame. Your result should look something like this:
## Age Day Month Year CaptureTime Species Wing Weight Tail ##1I19 91992 13:30 RT385 920219 ##2I22 91992 10:30 RT376 930221 ##3I23 91992 12:45 RT381 990235 ##4I23 91992 10:50 CH265 470220 ##5I27 91992 11:15 SS205 170157
For each of the following variables say whether they continuous, discrete or categorical. Discuss this with your colleagues.
1. Month 2. Species 3. Age
4. Wing 5. Weight
The information available at https://rdrr.io/cran/Stat2Data/man/Hawks.html may help. 1.2 What’s wrong with this plot?
Write down some problems with the plot displayed below.
500
400
300
200
100
150 200 250
Culmen 10
20 30
Weight 2000
1500 1000 500
Hopefully your report doesn’t contain plots like this!
2
1.3 Generate a histogram
Next use a combination of the functions ggplot() and geom_histogram to create a histogram plot of the weights of the Hawks within the hawksSmall data frame with bin widths of 100 grams. Your result should look something like this:
150
100
50
0
0 500
1000 1500 2000
Weight (gm)
Describe the aesthetic used within this plot.
Which term best describes the shape of the data distribution of Hawk weights: “Unimodal”, “Bimodal” or “Trimodal”?
1.4 Generate a density plot
Use a combination of the functions ggplot() and geom_density() to create a density plot of the tail lengths of the Hawks within the hawksSmall data frame. Your result should something like this:
3
Count
0.015
0.010
0.005
0.000
Create the following plots for yourself:
150
200 250
Tail (mm)
Recreate your plot with the argument adjust = 0.5 and adjust = 1. Describe the role played by the adjust argument within the geom_density() function. How many modes does the data distribution of Hawk tail lengths have?
0.03
0.02
0.01
0.00
Species CH
RT SS
Density Density
150
200 250
Tail (mm)
4
SS
RT
CH
Species CH
RT SS
1.5 Scatter plots
150 200 250
Tail (mm)
How many aesthetics are present within the following plot? What are the glyphs within this plot?
2000
1500
1000
500
0
Species CH
RT SS
Weight (gm) Species
150
200 250
Tail (mm)
Generate a plot similar to the above plot using the ggplot() and geom_point() functions. 1.6 Trend lines and facet wraps
Generate the following plot using the ggplot(), geom_point(), geom_smooth() and facet_wrap() functions.
5
CH
RT
SS
900
600
300
2000
1500
1000
500
900
600
300
0
Species CH
RT SS
160 180 200 220
150 200 250
Tail (mm)
125 150 175 200 225
What are the visual cues being used within this plot? Based on the above plot, what can we say about the relationship between the weight of the hawks and their tail lengths?
As an additional challenge, if you have time, you could try adding an annotation to your plot which highlights the weight of the heaviest hawk in the data set.
If you want to learn more about ggplot2 take a look at the ggplot2 gallery https://exts.ggplot2.tidyverse. org/gallery.
2 Data wrangling
We next turn to the data wrangling concepts from Lecture 4.
2.1 Select and filter functions
Use a combination of the select() and filter() functions to generate a data frame called “hSF” which is a sub-table of the original Hawks data frame with the following characteristics:
1. Your data frame should include the columns:
a) “Wing”,
b) “Weight” c) “Tail”.
2. Your data frame should contain a row for every hawk such that:
a) They belong to the species of Red-Tailed hawks b) They have weight at least 1kg.
6
Weight (gm)
Make use of the pipe operator to simplify your code.
How many variables does the dataframe hSF have? What would you say to communicate this information to a Machine Learning practitioner?
How many examples does the dataframe hSF have? How many observations? How many cases?
2.2 The arrange function
Use the arrange() function to sort the hSF data frame created in the previous section so that the rows appear in order of increasing wing span.
Use the head command to print out the top five rows of your sorted data frame. Your results should look something like this:
## Wing Weight Tail ## 1 37.2 1180 210 ## 2 111.0 1340 226 ## 3 199.0 1290 222 ## 4 241.0 1320 235 ## 5 262.0 1020 200
2.3 Join and rename functions
The species of Hawks within the data frame have been indicated via a two letter code. The correspondence between these codes and the full names is given by the following data frame.
## ## 1 ## 2 ## 3
species_code species_name_full
CH RT SS
Cooper’s Red-tailed Sharp-shinned
Recreate the above data frame containing the correspondence between codes and the full species names and give your data frame an appropriate name.
Use a combination of the functions left_join(), the rename() and the **select()* functions to create a new data frame called “hawksFullName” which is the same as the “Hawks” data frame except that the Species column contains the full names rather than the two letter codes.
## Joining, by = "Species"
Use a combination of the head() and select() functions to print out the top seven rows of the columns “Species”, “Wing” and “Weight” of the data frame called “hawksFullName”. Do this without modifying the data frame you just created. Your result should something like this:
## ## 1 ## 2 ## 3 ## 4 ## 5 Sharp-shinned 205 170 ## 6 Red-tailed 412 1090 ## 7 Red-tailed 370 960
Species Wing Weight Red-tailed 385 920 Red-tailed 376 930 Red-tailed 381 990
Cooper’s 265 470
Does it matter what type of join function you use here? In what situations would it make a difference?
7
2.4 The mutate function
Suppose that the fictitious “Healthy Hawks Society”1 has proposed a new measure called the “bird BMI” which attempts to measure mass of a hawk standardized by their wing span. The bird BMI is equal to the weight of the hawk (in grams) divided by their wing span (in millimeters) squared. That is,
Bird-BMI := 1000 × Weight/Wing-span2.
Use the mutate(), select() and arange() functions to create a new data frame called “hawksWithBMI” which has the same number of rows as the original Hawks data frame but only two columns – one with their Species and one with their “bird BMI”. The rows should appear in descending order of “bird BMI”. The top 8 rows of your data frame should look something like this:
## Species bird_BMI
- ##1 RT 852.69973
- ##2 RT 108.75741
- ##3 RT 32.57493
- ##4 RT 22.72688
- ##5 CH 22.40818
- ##6 RT 19.54932
- ##7 CH 15.21998
- ##8 RT 14.85927
Use the filter() function to remove those cases where the bird BMI exceeds 100 from your data frame. Then generate a violin plot of your data which shows the distribution of “bird BMIs” broken down by species. Your result should look something like this:
SS
RT
CH
Species CH
RT SS
Species
1 4 16
Bird BMI
2.5 Summarize and group-by functions
Using the dataframe “hawksFullName”, from problem 3 above, in combination with the summarize() and the groupby functions, create a summary table, broken down by Hawk species, which contains the following
1Both the “Healthy Hawks Society” and the concept of “bird BMI” were made up purely for this assignment. 8
summary quantities:
1. The number of rows;
2. The mean average wing span in centimeters;
3. The median wing span in centimeters;
4. The trimmed mean average wing span in centimeters (trim=0.1); 5. The mean average of the ratio between wing span and tail length.
Your final result should look something like this:
## # A tibble: 3 x 6 ## Species num_rows mn_wing md_wing t_mn_wing tail_wing_ratio
## <chr> <int> ## 1 Cooper’s 70 ## 2 Red-tailed 577 ## 3 Sharp-shinned 261
<dbl> <dbl> <dbl> <dbl> 244. 240 243. 1.22 383. 384 385. 1.73 185. 191 184. 1.26
Next create a summary table of the following form: Your summary table will show the number of missing values, broken down by species, for the columns Wing, Weight, Culmen, Hallux, Tail, StandardTail, Tarsus and Crop. You can complete this task by combining the select(), group_by(), summarize(), across(), everything(), sum() and is.na() functions. You should end with a summary table of the following form:
## # A tibble: 3 x 9
## Species ## <chr> ## 1 Cooper’s 1 ## 2 Red-tailed 0 ## 3 Sharp-shinned 0
Wing Weight Culmen Hallux Tail StandardTail Tarsus Crop
<int> <int> <int> <int> <int>
<int> <int> <int>
0 00 5 43 5 33
0 0 0
19 250 68
62 21 538 254 233 68
You can learn more about dplyr from the Tidyverse website https://dplyr.tidyverse.org/index.html.
3 Exploratory data analysis
We shall now illustrate some concepts from Lecture 5 on Exploratory Data Analysis in the context of our Hawks data set.
3.1 Combining location estimators with the summarise function
Use a combination of the summarise(), mean() and median() to compute the sample mean, sample median and trimmed sample mean (with q = 0.1) of the Hawk’s wing length and Hawk’s weight. Your result should look something like this:
## Wing_mean Wing_t_mean Wing_med Weight_mean Weight_t_mean Weight_med ## 1 315.6375 322.2297 370 772.0802 779.3681 970
Combine with the group_by() function to obtain a break down by species. Your result should look something like this:
9
## # A tibble: 3 x 7 ## Species Wing_mean Wing_t_mean Wing_med Weight_mean Weight_t_mean Weight_med
## <fct> <dbl> ## 1 CH 244. ## 2 RT 383. ## 3 SS 185.
<dbl> <dbl> <dbl> 243. 240 420. 385. 384 1094. 184. 191 148.
<dbl> <dbl> 410. 378.
1089. 1070 140. 155
3.2 Location and dispersion estimatiors under linear transformations
Suppose that a variable of interest Xi has values X1, . . . , Xn. Suppose that X1, . . . , Xn has sample mean μˆ. Let a, b ∈ R be real numbers and define a new variable X ̃i with values X ̃1,…,X ̃n defined by X ̃i = a·Xi +b fori=1,···,n. ShowthatX ̃1,…,X ̃n hassamplemeana·μˆ+b.
Suppose further that X1, . . . , Xn has sample variance SX2 . What is the sample variance of X ̃1, . . . , X ̃n? What is the sample standard deviation of X ̃1, . . . , X ̃n?
3.3 Robustness of location estimators
In this exercise we shall investigate the robustness of several location estimators: The sample mean, sample median and trimmed mean.
We begin by extracting a vector called “hal” consisting of the talon lengths of all the hawks with any missing values removed.
To investigate the effect of outliers on estimates of location we generate a new vector called “corrupted_hall” with 10 outliers each of value 100 created as follows:
We can then compute the mean of the original sample and the corrupted sample as follows.
mean(hal) ## [1] 26.41086
mean(corrupted_hal)
## [1] 27.21776
Now let’s investigate what happens as the number of outliers changes from 0 to 1000. The code below generates a vector called “means_vect” which gives the sample means of corrupted samples with different numbers of outliers. More precisely, means_vect is a vector of length 1001 with the i-th entry equal to the mean of a sample with i − 1 outliers.
hal<-Hawks$Hallux # Extract the vector of hallux lengths hal<-hal[!is.na(hal)] # Remove any nans
outlier_val<-100 num_outliers<-10 corrupted_hal<-c(hal,rep(outlier_val,times=num_outliers))
10
num_outliers_vect<-seq(0,1000) means_vect<-c()
for(num_outliers in num_outliers_vect){ corrupted_hal<-c(hal,rep(outlier_val,times=num_outliers)) means_vect<-c(means_vect,mean(corrupted_hal))
}
Copy and modify the above code to create an additional vector called “medians_vect” of length 1001 with the i-th entry equal to the median of a sample “corrupted_hal” with i − 1 outliers.
Ammend the code further to add an additional vector called “t_means_vect” of length 1001 with the i-th entry equal to the trimmed mean of a sample with i−1 outliers, where the trimmed mean has a trim fraction q = 0.1.
You should now have four vectors: “num_outliers_vect”, “means_vect”, “medians_vect” and “t_means_vect”. Combine these vectors into a data frame with the following code.
df_means_medians<-data.frame(num_outliers=num_outliers_vect,
mean=means_vect,t_mean=t_means_vect,
median=medians_vect)
Now use the code below to reshape and plot the data. The function pivot_longer() below is used to
reshape the data. Don’t worry if this operation is unclear at this stage. Its use will be explained soon.
df_means_medians%>% pivot_longer(!num_outliers, names_to = "Estimator", values_to = "Value")%>% ggplot(aes(x=num_outliers,color=Estimator,
linetype=Estimator,y=Value))+
geom_line()+xlab("Number of outliers")
The output of your code should look as follows:
100
80
60
40
Estimator mean
median t_mean
Value
0 250
500 750
Number of outliers
1000
11
3.4 Box plots and outliers
Use the functions ggplot() and geom_boxplot() to create a box plot which summarises the distribution of hawk weights broken down by species. Your plot should look as follows:
2000
1500
1000
500
0
CH RT SS
Species
Note the outliers displayed as individual dots.
Suppose we have a sample X1, · · · , Xn. Let q25 denote the 0.25-quantile of the sample and let q75 denote the 0.75-quantile of the sample. We can then define the interquartile range, denoted IQR by IQR := q75 − q25. In the context of boxplots and outlier Xi is any numerical value such that the following holds if either of the following holds:
Xi <q25−1.5×IQR Xi >q75+1.5×IQR.
Create a function called “num_outliers” which computes the number of outliers within a sample (with missing values excluded).
Now combine your function num_outliers() with the functions group_by() and summarise() to com- pute the number of outlier for the three samples of hawk weights broken down by specied. Your result should look as follows:
## # A tibble: 3 x 2
-
## Species num_outliers_weight
-
## <fct> <int> ## 1 CH 3 ## 2 RT 13 ## 3 SS 4
3.5 Covariance and correlation under linear transformations
Suppose that we have a pair of variables: Xi with values X1, . . . , Xn and Yi with values Y1, . . . , Yn. Suppose that X1,…,Xn and Y1,…,Yn have sample covariance ΣX,Y . Let a, b ∈ R be real numbers and define a new variable X ̃i with values X ̃1,…,X ̃n defined by X ̃i = a·Xi +b for i = 1,··· ,n. In addition, let c, d ∈ R
12
Weight
be real numbers and define a new variable Y ̃i with values Y ̃1,…,Y ̃n defined by Y ̃i = c·Yi +d for i = 1,··· ,n. What is the covariance between the pair of variables X ̃1, . . . , X ̃n and Y ̃1, . . . , Y ̃n?
Suppose that X1, . . . , Xn and Y1, . . . , Yn have correlation ρX,Y . What is the correlation between the pair of variables X ̃1,…,X ̃n and Y ̃1,…,Y ̃n?
13
