---
title: "SISMID 2021 R Notes: Introduction to R"
author: |
  | Jon Wakefield
  | University of Washington
date: "`r Sys.Date()`"
output:
  beamer_presentation:
    theme: "AnnArbor"
    colortheme: "dolphin"
    keep_tex: yes
  slidy_presentation: default
---

```{r setup, echo=FALSE, message=FALSE, warning=FALSE}
library(knitr)
opts_chunk$set(collapse=TRUE, fig.align='center', tidy=TRUE, tidy.opts=list(blank=TRUE, width.cutoff=40,strip.white=TRUE), warning=FALSE,message=FALSE,cache=TRUE)
```

# Preliminaries

## Objectives


In this set of ``R`` notes you will:

- See the basic types of object that ``R`` uses, including spatial objects

- Learn how to install packages

- See how to read in data from a website 

- Learn about basic plotting and regression models

Spatial data:

- Lancashire case-control point data

- Washington counties

- Scottish districts

## Accessing the ``R`` code

To get ``R`` code alone, load the ``knitr`` library as below and then use the ``purl`` command -- I comment out these lines to save time
\vspace{.2in}
\small
```{r}
#install.packages("knitr")
#library(knitr)
#purl("2020-SISMID-Lec1-R.Rmd")
```

You need to be in the directory that contains the ``.Rmd`` file -- the output file is ``2020-SISMID-Lec1-R.R`` and is also in the current directory, which can be found by typing ``getwd()`` 

## Introduction to ``R``

\small
This document is an ``R`` Markdown file. It is a relatively easy way to write a report (or notes!) that contains data analysis -- great for reproducibility. To learn more about R Markdown, check out the online book at

[https://bookdown.org/yihui/rmarkdown/](https://bookdown.org/yihui/rmarkdown/)

We describe various types of objects (vectors, matrices, arrays, data frames, lists). This can be confusing but the ``class`` command is useful.

Some functions require certain types of objects as inputs, so if a function doesn't do as expected, it can be that the wrong type of input is being used.

``R``'s basic data types are character, numeric, integer, complex, and logical. 

``R``'s basic data structures include vector, list, matrix, data frame, and factor.

## Vectors

In ``R``, we use ``<-`` or ``=`` to assign the value on the right hand side to an object on the left hand side. 

For example, 
\vspace{.2in}
\small
```{r}
x1 <- "Hello"
x1
class(x1)
y1 <- c(1, 2, 3)
y1
class(y1)
```

## Vector

The ``c(...)`` notation defines a vector, i.e., a collection of elements of the same type. For some special vectors, there are faster ways to define them

For example, 
\vspace{.2in}
\small
```{r, message=FALSE, output=20, collapse=TRUE, tidy=TRUE,tidy.opts=list(width.cutoff=20)}
x2 <- 1:10
x2 
class(x2)
x2a <- seq(1,5,.5)
x2a
class(x2a)
```

---

Similarly, there are two-dimensional extensions of a vector (a matrix), and higher-dimensional extensions of a vector (an array). We give some examples of subsetting, with more to follow.
\vspace{.2in}
\small
```{r}
x3 <- matrix(1:24, nrow = 3)
x3
x3[2,3]
class(x3)
```

## Objects of type matrix and array

Similarly, there are two-dimensional extensions of a vector (a matrix), and higher-dimensional extensions of a vector (an array). When data is input the leftmost subscript moves fastest
\vspace{.2in}
\scriptsize
```{r, message=FALSE, collapse=TRUE, tidy=TRUE,tidy.opts=list(width.cutoff=50)}
y2 <- array(1:24, dim = c(2, 3, 4))
y2[,,1]
y2[,,4]
y2[1,2,2]
class(y2)
```

## Objects of type matrix and array

Many operators in ``R`` are vectorized, meaning that operations occur in parallel in certain ``R`` objects.
\vspace{.2in}
\scriptsize
```{r}
x4 <- 3
x4 <- x4^2 + 1
x4
y3 <- c(2,3,0)
y3 <- log(y3)
y3
```


## Accessing elements

There are three operators that can be used to extract subsets of ``R`` objects.

The ``[`` operator always returns an object of the same class as the original. It can be used to select multiple elements of an object

The ``[[`` operator is used to extract elements of a list or a data frame. It can only be used to extract a single element and the class of the returned object will not necessarily be a list or data frame.

The ``$`` operator is used to extract elements of a list or data frame by literal name. Its semantics are similar to that of ``[[``

## Data frame
A tabular form of data where each column may contain different types of values (i.e., similar to a spreadsheet), is called a data frame in ``R``.
\vspace{.2in}
\scriptsize
```{r, tidy.opts=list(width.cutoff=50)}
x5 <- data.frame(ID = 1:4,name = c("A", "B", "C", "D"), 
                color = factor(c("red", "blue", "blue", "green")))
class(x5)
sapply(x5,class)
x5$ID
x5[,1]
x5[1,]
```

##


\scriptsize
```{r, message=FALSE, collapse=TRUE, tidy=TRUE,tidy.opts=list(width.cutoff=50)}
x5[,2]
x5$color
x5[1:2, ]
x5[1:2, 2:3]
x5[, c("ID", "name")]
```

## Factors

Notice that although factors look similar to characters, they can behave very differently in action. 

We encourage new users of ``R`` to read the section on 'R Nuts and Bolts' of this online book
[\textcolor{red}{R Programming for Data Science}](https://bookdown.org/rdpeng/rprogdatascience/r-nuts-and-bolts.html#factors) by Roger Peng,
(if not the whole book) for more details.



## Lists

Another commonly used data structure is 'list'. It creates a list of elements that are not necessarily of the same type, dimension, or have the same attributes. For example, 
\vspace{.2in}
\scriptsize
```{r, message=FALSE, collapse=TRUE, tidy=TRUE,tidy.opts=list(width.cutoff=50)}
y4 <- list(apple = 1:20, orange = x5)
class(y4)
y4
```

## Lists

To subset a list: 
\vspace{.2in}
\scriptsize
```{r, message=FALSE, collapse=TRUE, tidy=TRUE,tidy.opts=list(width.cutoff=50)}
y4[[1]] 
y4[["apple"]]
y4$apple
```

## Lists

To subset a list: 
\vspace{.2in}
\scriptsize
```{r, message=FALSE, collapse=TRUE, tidy=TRUE,tidy.opts=list(width.cutoff=50)}
y4[["orange"]]
y4$orange
```

## Coercion

Objects can be explicitly coerced from one class to another using the ``as.*`` functions, if available.
\vspace{.2in}
\scriptsize
```{r, message=FALSE, collapse=TRUE, tidy=TRUE,tidy.opts=list(width.cutoff=50)}
x <- 0:6
class(x)
as.numeric(x)
as.logical(x)
as.character(x)
```

## ``R`` packages: Installation

To install packages that are hosted on [CRAN](https://cran.r-project.org), we only need to use the ``install.packages`` function.
\vspace{.2in}
\scriptsize
```{r eval = FALSE}
#install.packages("ggplot2", dep = TRUE)
#install.packages("rgdal", dep = TRUE)
```
\normalsize
Some packages not hosted on CRAN. For example, the INLA package that we will be using for space-time smoothing is hosted on its own [website](http://www.r-inla.org).
\vspace{.2in}
\scriptsize
```{r eval = FALSE}
# Commented out to save time
# install.packages("INLA", 
#                  repos=c(getOption("repos"),
#                          INLA="https://inla.r-inla-download.org/R/stable"))
```

## Load ``R`` packages

Once a package is installed, you can load it into your workspace. Then you can use all the functions and data provided in that package. 
\vspace{.2in}
\scriptsize
```{r}
library(ggplot2)
library(rgdal)
library(SUMMER)
```
\normalsize

You can check your loaded packages in the workspace by (I haven't evaluated here, as I have lots of packages installed!).
\vspace{.2in}
\scriptsize
```{r eval = TRUE}
#sessionInfo()
```

# Spatial ``R``

## ``R`` for Spatial Analysis

``R`` has extensive spatial capabilities, the 
Spatial task view is [\textcolor{red}{here}](http://cran.r-project.org/web/views/Spatial.html) 

Bivand et al (2013), is **the** reference book for GIS in ``R``, though unfortunately it's quite hard to follow if you're an ``R`` novice.



## Representing Spatial Data

Spatial classes were defined to represent and handle spatial data, so that data can be exchanged between different classes.

Class definitions are objects that contain the formal definition of a class of ``R`` objects, usually referred to as an S4 class.

Just as components of a list are accessed using ``$``, slots of an object are accessed using ``@``

## Representing Spatial Data

The ``sp`` library is the workhorse for representing spatial data.

The most basic spatial object is a 2d or 3d point: a set of coordinates may be used to define a ``SpatialPoints`` object.

From the help function:

``SpatialPoints(coords, proj4string=CRS(as.character(NA)),bbox = NULL)``

\normalsize

- PROJ.4 is a library (not an ``R`` package) for performing conversions between cartographic projections.

- The points in a ``SpatialPoints`` object may be associated with a set of attributes to give a ``SpatialPointsDataFrame`` object.

## Creating a Spatial Object

As an example, the `splancs` library was pre-``sp`` and so does not use spatial objects.

`splancs` contains a number of useful functions for analyzing spatial referenced point data.
\vspace{.2in}
\scriptsize
```{r, message=FALSE, collapse=TRUE, tidy=TRUE,tidy.opts=list(width.cutoff=50)}
library(sp)
library(splancs)
data(southlancs) # case control data
summary(southlancs)
```

## Creating a Spatial Object

Let's look at the case control points as regular ``R`` data
\vspace{.2in}
\scriptsize
```{r, message=FALSE, collapse=TRUE, tidy=TRUE,tidy.opts=list(width.cutoff=50)}
class(southlancs)
names(southlancs)
points <- southlancs[,c("x","y")]
summary(points)
```

## Creating a Spatial Object

We convert into a ``SpatialPoints`` object and then create a ``SpatialPointsDataFrame`` data frame.
\vspace{.2in}
\scriptsize
```{r, message=FALSE, collapse=TRUE, tidy=TRUE,tidy.opts=list(width.cutoff=50)}
casepoints <- SpatialPoints(southlancs[,c("x","y")])
summary(casepoints)
```


## Creating a Spatial Object


\scriptsize
```{r, message=FALSE, collapse=TRUE, tidy=TRUE,tidy.opts=list(width.cutoff=50)}
class(casepoints)
proj4string(casepoints)
```

\normalsize
No info on where these points are on planet Earth...


## Examining a Spatial Object


Try out some commands
\vspace{.2in}
\scriptsize
```{r, message=FALSE, collapse=TRUE, tidy=TRUE,tidy.opts=list(width.cutoff=50)}
head(coordinates(casepoints),3)
bbox(casepoints)
casepoints@coords[1:3,]
```

## Examining a Spatial Object

\scriptsize
```{r, message=FALSE, collapse=TRUE, tidy=TRUE,tidy.opts=list(width.cutoff=50)}
str(casepoints)
```

## Creating a Spatial Object

Now create a ``SpatialPointsDataFrame`` data frame.
\vspace{.2in}
\scriptsize

```{r, message=FALSE, collapse=TRUE, tidy=TRUE,tidy.opts=list(width.cutoff=50)}
caseDF <- SpatialPointsDataFrame(coords=casepoints,
                   data=as.data.frame(southlancs$cc))
class(caseDF)
names(caseDF)
```

## Creating a Spatial Object

Now create a ``SpatialPointsDataFrame`` data frame.
\vspace{.2in}
\scriptsize

```{r, message=FALSE, collapse=TRUE, tidy=TRUE,tidy.opts=list(width.cutoff=50)}
summary(caseDF)
```

## Visualizaing Spatial Data

Plotting spatial data can be provided in a variety of ways, see Chapter 3 of Bivand et al. (2013).

The most obvious is to use the regular plotting functions, by converting ``Spatial`` dataframes to regular dataframes, for example using ``as.data.frame``.




## Reading Shapefiles

ESRI (a company one of whose products is ArcGIS) shapefiles consist of three files, and this is a common form. 

+ The first file (*.shp) contains the geography of each shape. 

+ The second file (*.shx) is an index file which contains record offsets. 

+ The third file (*.dbf) contains feature attributes with one record per feature.


## Reading Shapefiles

Many repositories host data:

Globally, the Database of Global Administrative Areas [\textcolor{red}{here}](https://gadm.org) provides maps and spatial data for all countries and their subdivisions (admin-0, admin-1, admin-2,...)


Locally: The UW Geospatial Data Resources: 
[\textcolor{red}{here}](https://guides.lib.uw.edu/friendly.php?s=research/gis). 

As an example, consider  Washington county data that was downloaded from this site.

The data 
consists of the three files: wacounty.shp, wacounty.shx, wacounty.dbf.

The following code reads in these data and then draws a county level
map of 1990 populations.


## Reading Shapefiles

First load the libraries.

```{r, message=FALSE, collapse=TRUE, tidy=TRUE,tidy.opts=list(width.cutoff=50)}
library(maps)
library(maptools)
```

The three required data files (wacounty.*) can be downloaded from here:

\small
http://faculty.washington.edu/jonno/SISMIDmaterial/
 
\normalsize 
and can be put in a directory ``R-examples`` (which you'll need to create). 

You can put in any old directory, including the current directory, so long as you use the correct name when you call.

The following code does this for you



## Reading Shapefiles

\scriptsize

```{r, message=FALSE, collapse=TRUE,echo=F, tidy=TRUE,tidy.opts=list(width.cutoff=50)}
if (!("R-examples" %in% list.files())) {
  dir.create("R-examples")
}
```

```{r, message=FALSE, collapse=TRUE, tidy=TRUE,tidy.opts=list(width.cutoff=40),eval=T}

download.file("http://faculty.washington.edu/jonno/SISMIDmaterial/wacounty.dbf",
              destfile = "R-examples/wacounty.dbf")
download.file("http://faculty.washington.edu/jonno/SISMIDmaterial/wacounty.shp",destfile = "R-examples/wacounty.shp")
download.file("http://faculty.washington.edu/jonno/SISMIDmaterial/wacounty.shx",destfile = "R-examples/wacounty.shx")

```


## Reading Shapefiles

\scriptsize

```{r, message=FALSE, collapse=TRUE, tidy=TRUE,tidy.opts=list(width.cutoff=40),eval=T}

# Now read all the files into a single data frame
wacounty <- rgdal::readOGR(dsn="R-examples",layer="wacounty")
class(wacounty)

# If the chunk throws an error, comment out the PREVIOUS automated downloading chunk (or changes knit option 'eval=T' to  'eval=F' for the PREVIOUS chunk). And manually download the .shp, .shx and .dbf files from the website. 


```


```{r, message=FALSE, collapse=TRUE, tidy=TRUE,tidy.opts=list(width.cutoff=30),echo=FALSE}
wacounty <- rgdal::readOGR(dsn="R-examples",layer="wacounty")
class(wacounty)
```

## Looking at the data

\footnotesize
```{r, message=FALSE, collapse=TRUE, tidy=TRUE,tidy.opts=list(width.cutoff=30)}
names(wacounty)
# Let's see what these variables look like
wacounty$AreaName[1:3] # county names
wacounty$AreaKey[1:3] # FIPS codes
```

## Examining ``SpatialPolygonDataFrames``

\scriptsize
```{r }
proj4string(wacounty)
head(wacounty)
```

## Examining the data 

Let's look in the ``data`` slot
\vspace{.2in}
\scriptsize

```{r }
head(wacounty@data)
```


## Drawing a map

We look at some variables.
\vspace{.2in}
\scriptsize

```{r, message=FALSE, collapse=TRUE, tidy=TRUE,tidy.opts=list(width.cutoff=40)}
wacounty$INTPTLAT[1:3] # latitude
wacounty$INTPTLNG[1:3] # longitude
wacounty$CNTY[1:3]    
wacounty$TotPop90[1:3]
```

## Drawing a map

We set up the colors to map.
\vspace{.2in}
\scriptsize

```{r, message=FALSE, collapse=TRUE, tidy=TRUE,tidy.opts=list(width.cutoff=40),fig.show="hide"}
library(RColorBrewer)
plotvar <- wacounty$TotPop90 # variable we want to map
summary(plotvar)
nclr <- 8 # next few lines set up the color scheme for plotting
plotclr <- brewer.pal(nclr,"Oranges")
brks <- round(quantile(plotvar,probs=seq(0,1,1/(nclr))),digits=1)
colornum <- findInterval(plotvar,brks,all.inside=T)
colcode <- plotclr[colornum]
plot(wacounty,col=colcode)
legend(-119,46,legend=leglabs(round(brks,digits=1)),fill=plotclr,cex=0.35,bty="n")
```

## Drawing a map

```{r,collapse=TRUE, fig.height=4.5,fig.width=4.5,echo=FALSE,fig.cap="1990 Washington population counts by census tracts"}

plot(wacounty,col=colcode)
legend(-119,46,legend=leglabs(round(brks,digits=1)),fill=plotclr,cex=0.35,bty="n")

```

## Drawing a map

As an alternative we can use the ``spplot`` function, which uses lattice (trellis) plot methods for spatial data with attributes.


## Drawing a map: an alternative way

\scriptsize

```{r,collapse=TRUE, fig.height=4.8,fig.width=4.8,fig.asp=1.0,echo=TRUE,fig.cap="1990 Washington population counts by county"}

spplot(wacounty,zcol="TotPop90",col.regions=brewer.pal(9,"Purples"),cuts=8)

```

## Highlight a county

\scriptsize

```{r highlight, collapse=TRUE,  echo=T, fig.show="hide"}
# identify counties of interest
xx=which(wacounty$CNTY==33)
xx2=which(wacounty$CNTY==63)
# plot the whole state
plot(wacounty)
# highlight counties of interest
plot(wacounty[xx, ], col="pink", border="red", add=T,lwd=2.5)
plot(wacounty[xx2, ], col="lightblue", border="darkblue", add=T, lwd=2.5)
# Add some labels
text(coordinates(wacounty[xx, ]), "King", cex=0.75, pos=3, offset=0.15)
text(coordinates(wacounty[xx2, ]), "Spokane", cex=0.7, pos=1, offset=0.25)
points(coordinates(wacounty[c(xx, xx2), ]), pch=16, cex=0.75)

```


## Highlight a county
\scriptsize
```{r , ref.label="highlight", collapse=TRUE, fig.height=5,fig.width=5, echo=F}

plot(wacounty, border="#00000075")

# highlight counties of interest
plot(wacounty[xx, ], col="#ff000075", add=T)
plot(wacounty[xx2, ], col=NA, border="#0000ff75", add=T, lwd=2.5)

# Add some labels
text(coordinates(wacounty[xx, ]), "King", cex=0.75, pos=3, offset=0.25)
text(coordinates(wacounty[xx2, ]), "Spokane", cex=0.7, pos=1, offset=0.25)
points(coordinates(wacounty[c(xx, xx2), ]), pch=16, cex=0.75)

```




## Scottish lip cancer data

We will first fit a number of models to the famous Scottish lip cancer data.

We have counts of disease, expected numbers and an area-based covariate (proportion in agriculture, fishing and farming) in each of 56 areas.

\vspace{.2in}
\scriptsize

```{r, message=FALSE, collapse=TRUE, tidy=TRUE,tidy.opts=list(width.cutoff=50),eval=T}
library(SpatialEpi)
data(scotland)
Y <- scotland$data$cases
X <- scotland$data$AFF
E <- scotland$data$expected
# Relative risk estimates
smr <- Y/E
summary(E)
summary(smr)
```

## Distribution of SMRs

\scriptsize
```{r, echo=TRUE, collapse=TRUE,fig.height=2.5,fig.width=3,tidy.opts=list(width.cutoff=50),eval=T}

ggplot(data.frame(smr),aes(x=smr))+geom_histogram(color="blue",fill="lightblue",bins=20)+labs(x="SMR",y="Disease Counts")

```

## Mapping

\scriptsize

```{r, echo=TRUE, collapse=TRUE,fig.height=2.8,fig.width=4, tidy.opts=list(width.cutoff=50),eval=T}

smap <- scotland$spatial.polygon
scotd <- scotland$data[,c("county.names", "cases", "expected", "AFF")]
scotd$SIR <- scotd$cases/scotd$expected

```

## Mapping

\scriptsize
```{r, echo=TRUE, collapse=TRUE,fig.height=2.8,fig.width=4, tidy.opts=list(width.cutoff=50),eval=T}

# set the polygons ID slot values to correspond to county names
sapply(slot(smap, "polygons"), function(x){slot(x, "ID")})
rownames(scotd) <- scotd$county

```

## Mapping

\scriptsize
```{r, echo=TRUE, collapse=TRUE,fig.height=2.8,fig.width=4, tidy.opts=list(width.cutoff=50),eval=T}

smap <- SpatialPolygonsDataFrame(smap,scotd, match.ID = TRUE)
spplot(smap,zcol="SIR",col.regions=brewer.pal(9,"Purples"),cuts=8)

```

## Exercises

+ Experiment with the ``spplot`` command, changing colors, for example

+ Map the expected numbers from the Scottish lip cancer example

