Getting Started

If you are using this code with the video, note that some slides have been added post-recording and are not shown. For this video we will be working with the USpop.csv dataset again. Make an R Notebook for this walk through tutorial to save all the code you will be learning. We will cover:

• Adding new vectors to datasets
• Numeric versus character values
• for loops in R
• transformations
• binning data
• using maps to plot points on a map
• Splitting datasets into quantiles
• Paralleling base R code to dplyr syntax to accomplish the same goals

Set up workspace

• You are working within an R project (check in the top right corner of RStudio - you should see the project name “R_Mini_Course”).

• This means that the project directory will also be set as working directory. The exception is in a R Notebook, where the working directory is where the R Notebook is saved.

• You should be saving your notebooks in the R Project and using ../.. to point at the main project directory.

• Before starting, it may be helpful to have a chunk of code that does the following:

  • clear your workspace rm(list=ls())
  • load your packages library(<package-name>)
  • check your session information sessionInfo()
  • list files in your working directory list.files(getwd())

rm(list=ls())
library(knitr)
library(maps)
library(tidyverse)

## ── Attaching packages ──────────────────────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.4.0      ✔ purrr   0.3.5 
## ✔ tibble  3.1.8      ✔ dplyr   1.0.10
## ✔ tidyr   1.2.1      ✔ stringr 1.4.1 
## ✔ readr   2.1.3      ✔ forcats 0.5.2 
## ── Conflicts ─────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ✖ purrr::map()    masks maps::map()

sessionInfo()

## R version 4.2.1 (2022-06-23)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur 11.7
## 
## Matrix products: default
## LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] forcats_0.5.2   stringr_1.4.1   dplyr_1.0.10    purrr_0.3.5     readr_2.1.3    
##  [6] tidyr_1.2.1     tibble_3.1.8    ggplot2_3.4.0   tidyverse_1.3.2 maps_3.4.1     
## [11] knitr_1.40      rmarkdown_2.17 
## 
## loaded via a namespace (and not attached):
##  [1] lubridate_1.9.0     assertthat_0.2.1    digest_0.6.30       utf8_1.2.2         
##  [5] R6_2.5.1            cellranger_1.1.0    backports_1.4.1     reprex_2.0.2       
##  [9] evaluate_0.18       httr_1.4.4          highr_0.9           pillar_1.8.1       
## [13] rlang_1.0.6         googlesheets4_1.0.1 readxl_1.4.1        rstudioapi_0.14    
## [17] jquerylib_0.1.4     googledrive_2.0.0   munsell_0.5.0       broom_1.0.1        
## [21] compiler_4.2.1      modelr_0.1.9        xfun_0.34           pkgconfig_2.0.3    
## [25] htmltools_0.5.3     tidyselect_1.2.0    fansi_1.0.3         crayon_1.5.2       
## [29] tzdb_0.3.0          dbplyr_2.2.1        withr_2.5.0         grid_4.2.1         
## [33] jsonlite_1.8.3      gtable_0.3.1        lifecycle_1.0.3     DBI_1.1.3          
## [37] magrittr_2.0.3      scales_1.2.1        cli_3.4.1           stringi_1.7.8      
## [41] cachem_1.0.6        fs_1.5.2            xml2_1.3.3          bslib_0.4.1        
## [45] ellipsis_0.3.2      generics_0.1.3      vctrs_0.5.0         tools_4.2.1        
## [49] glue_1.6.2          hms_1.1.2           fastmap_1.1.0       yaml_2.3.6         
## [53] timechange_0.1.1    colorspace_2.0-3    gargle_1.2.1        rvest_1.0.3        
## [57] haven_2.5.1         sass_0.4.2

list.files(getwd())

##  [1] "_site.yml"                                   "4.02.Basic_Summary_Statistics_in_R_files"   
##  [3] "4.02.Basic_Summary_Statistics_in_R.html"     "4.02.Basic_Summary_Statistics_in_R.Rmd"     
##  [5] "4.03.Data_Manipulation_in_R.Rmd"             "4.04.Advanced_Statistical_Concepts_in_R.Rmd"
##  [7] "4.05.R_on_CL.Rmd"                            "4.06.Programming_in_R.Rmd"                  
##  [9] "5.03.Advanced_Graphing_in_R.Rmd"             "about.Rmd"                                  
## [11] "Activity1_intro.pdf"                         "activity1.html"                             
## [13] "activity1.Rmd"                               "activity1key.html"                          
## [15] "activity1key.Rmd"                            "activity2.html"                             
## [17] "activity2.Rmd"                               "activity3.html"                             
## [19] "activity3.Rmd"                               "Body_Fat.pdf"                               
## [21] "Congrats.html"                               "Congrats.Rmd"                               
## [23] "data"                                        "docs"                                       
## [25] "images"                                      "index.html"                                 
## [27] "index.Rmd"                                   "LICENSE"                                    
## [29] "module1.html"                                "module1.Rmd"                                
## [31] "module2.html"                                "module2.Rmd"                                
## [33] "module3.html"                                "module3.Rmd"                                
## [35] "module4.html"                                "module4.Rmd"                                
## [37] "module5.html"                                "module5.Rmd"                                
## [39] "module6.html"                                "module6.Rmd"                                
## [41] "module7.html"                                "module7.Rmd"                                
## [43] "modules"                                     "R_Mini_Course.Rproj"                        
## [45] "raw_data"                                    "README.md"                                  
## [47] "README.Rmd"                                  "site_libs"

Read in dataset

You will need to add path information to the raw_data directory once you have uncompressed the data tarball. As we did in the last video, read the USpop.csv file into an object called “data”:

#data=read.csv(file="USPop.csv")
#View(data)

You may also read in the previously made object into “data”:

data <- readRDS(file = "data/USPop.rds")
head(data)

##   Year Population
## 1 1790    3929214
## 2 1800    5308483
## 3 1810    7239881
## 4 1820    9638453
## 5 1830   12866020
## 6 1840   17069453

Click the data object in the environment tab in RStudio. Notice we have two columns of data, year and population size.

Add a new column

An easy new column to add is a midpoint between two existing columns. In this way, we can add five to every item in the vector data$Year to create a new vector data$midpoint:

data$midpoint=data$Year+5
#View(data)

Click the data object in the environment tab in RStudio again. Notice there is now a third column labeled midpoint.

The “tidy” way

Occasionally, you may see alternative methods of performing some task in R using tools from the meta package tidyverse. We will introduce some basic utility from this meta package throughout the modules. Let’s add the new column using the dplyr and magrittr packages from tidyverse

tidydata <- data %>% mutate(midpoint = Year+5)

%>% from the magrittr package is a pipe! This means it is used to pass information from the left of the symbol to functions on the right.

Here we are passing the contents of data to the dplyr function mutate() creating a new column named “midpoint.” Midpoint contains the value in Year plus 5.

Initialize a vector

Another way to add to a dataset is to initialize a vector outside of a data frame.

growth=vector(mode="numeric",length=length(data$Year))

In setting up this vector, we have made it the same length as our existing vector data$Year. We also set it up to be numeric indicating we want R to treat it as numbers.

We will store population growth in this vector, and will thus need to compare the population size in year i to year i-1. What concepts have you learned in this course that will make this easy?

Fill vector via a for-loop

Similar to how we’ve done this in shell, we set a loop variable i and set the loop to start at 2 and stop at the length of data$Year. Because growth compares intervals of time, we will have a vector that is one less than the current vector length. However, R doesn’t like this, so we will set the first value as missing data.

growth[1]=NA
for (i in 2:length(data$Year)) {
    growth[i]=data$Population[i]-data$Population[i-1]+1
}

Examine the formula in the loop. Why do we add 1 to the growth for each iteration?

Another way to add columns

data=cbind(data,growth)

Here, we’ve reset the object data as a combination of the existing object and our new vector.

The “tidy” way

We combine data and growth using dplyr:

tidydata <- bind_cols(tidydata, as_tibble(growth)) %>% 
  rename(growth = value)

Here, we are using bind_cols to combining “tidydata” & “growth”. Notice that we coerced “growth” into a tibble (a tidy data frame) to match the structure of “tidydata”. We wanted the column to be called “growth”, so we renamed it with rename!

dplyr (a grammar of data manipulation) was developed with the end user in mind. The package contains a few clear, explanatory verbs that perform common data manipulations. Here are a few other manipulations and functions you should be aware about:

Getting subsets of a dataset with filter

# only rows where population size is greater than or equal to the mean
tidydata_pop_filter <- tidydata %>% filter(Population >= mean(Population))

Use the help pane to search for slice_max() or in the console type ?slice_max. Try to craft and execute a command using this function to 50% of the rows with the highest growth values. Guidance: slice_max(growth, prop = .5)

Getting specific columns or variables with select

# isolate the growth and filtered population columns
tidydata_pop_select <- tidydata_pop_filter %>% select(Population, growth)

Getting the mean population size and growth with summarise & across

# summarizing across variables
tidydata_pop_means <- tidydata_pop_select %>% summarise(across(everything(),mean))

Final words on dplyr

tidyverse was created by Hadley Wickham (Statistician and Chief Scientists at RStudio) and associates. It is meant to make the coding syntax and structure more intuitive and understandable. dplyr is one of many other packages within the tidyverse, but they all share the same underlying design and data structures (i.e., they work well together). Interested in learning more? Check out tidyverse. Check out the book R for Data Science. Check out the dplyr cheatsheet.

Visual inspection of new vector

Now, back to the data set. As we did before, let’s explore this new data visually

hist(growth)

Does this data follow a normal distribution?

Log transformation

Transformation is a nice way of converting your data into a normal distribution, an underlying assumption of many statistical tests. Below, we impose a log transformation.

data$log=log(data$growth)
hist(data$log)

How does this distribution compare to the last one?

Now, let’s look at growth over time:

#plot from last video
plot(data[,1], data[,2],pch=11, col="#66CDAA",cex=0.8)

#plot of growth (without transformation)
plot(data$Year,data$growth)

#plot of log growth
plot(data$Year,data$log)

Does the transformation change much for the plot?

Test for a statistical correlation

Let’s see how our transformation impacts a test for a statistical correlation.

cor.test(data$Year,data$growth,method="spearman")

## 
##  Spearman's rank correlation rho
## 
## data:  data$Year and data$growth
## S = 76, p-value = 3.799e-06
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.9570864

cor.test(data$Year,data$log,method="spearman")

## 
##  Spearman's rank correlation rho
## 
## data:  data$Year and data$log
## S = 76, p-value = 3.799e-06
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.9570864

Why are these values the same?

Quake Dataset

Next, let’s examine another R built-in dataset, quakes

This dataset includes 1000 records of earthquakes and their location on the globe. If we wanted to plot all of them at once, it would be overwhelming.

Basic Summary Statistics

As we did before, let’s look some basic summaries of this new dataset. The summary function breaks data into quantiles, but quantile does this a bit cleaner without mean in the middle:

summary(quakes$mag)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    4.00    4.30    4.60    4.62    4.90    6.40

quantile(quakes$mag)

##   0%  25%  50%  75% 100% 
##  4.0  4.3  4.6  4.9  6.4

As you can see, the majority of earthquakes have a magnitude below 5. If we wanted to plot earthquakes on a map, we may want a way to differentiate them based on their magnitude. We can do this via binning.

Binning data

Using the new function qunatile, let’s break up the quake data based on magnitude and add this new vector to the dataset quakes:

quakes$quant=cut(quakes$mag,breaks=quantile(quakes$mag),labels=F,include.lowest=T)
head(quakes)

##      lat   long depth mag stations quant
## 1 -20.42 181.62   562 4.8       41     3
## 2 -20.62 181.03   650 4.2       15     1
## 3 -26.00 184.10    42 5.4       43     4
## 4 -17.97 181.66   626 4.1       19     1
## 5 -20.42 181.96   649 4.0       11     1
## 6 -19.68 184.31   195 4.0       12     1

#View(quakes)

Because we split the data into quantiles, there should be roughly the same number in each group. At your table, verify this is in fact true. How many data points are in each quantile?

Binning data and adding character labels

Rather than have your data labeled 1-4, you may be interested in a more specific labeling scheme. Here, I am going to assign color values to each category with the lowest magnitude having grey value and the highest category having red value.

quakes$color=as.character(cut(quakes$mag,breaks=quantile(quakes$mag),labels=c("grey50","yellow","orange","red"),include.lowest=T))

Check out maps

Take a few minutes to read over the documentation for this new package.

??maps

Plotting quakes on a map

We will use the function map. First, let’s examine the range of latitude and longitude values in this dataset:

range(quakes$lat)

## [1] -38.59 -10.72

range(quakes$long)

## [1] 165.67 188.13

Plotting quakes on a map

Let’s use the world map to project these earthquakes. We’ll initialize a map and add points from our dataset. Note we’ve specified the limits of the plot according to our data, but added a bit of a buffer. This allows us to zoom out a bit to see where we are. Feel free to make your own adjustments here:

maps::map('world', xlim=c(min(quakes$long)-30,max(quakes$long)+30),ylim=c(min(quakes$lat)-15,max(quakes$lat)+15))
points(quakes$long,quakes$lat,pch=19)

Adjusting plot based on our bins

Finally, let’s use adjust how these data are plotted based on the quantile bins we made earlier:

#use color value to mark magnitude
maps::map('world',xlim=c(min(quakes$long)-35,max(quakes$long)+15),ylim=c(min(quakes$lat)-15,max(quakes$lat)+5))
points(quakes$long,quakes$lat,col=quakes$color,pch=19)

#use quantile value to adjust size of each point according to magnitude
maps::map('world',xlim=c(min(quakes$long)-35,max(quakes$long)+15),ylim=c(min(quakes$lat)-15,max(quakes$lat)+5))
points(quakes$long,quakes$lat,col=quakes$color,pch=19,cex=0.5*quakes$quant)