A Tutorial Introduction to R

# A Tutorial Introduction to R

## Introductory Computer Programming

### Deepayan Sarkar

---

# Software for Statistics

- Computation is an essential part of modern statistics

- Handling large datasets

- Visualization

- Simulation

- Iterative methods 
--

- Many softwares are available for such computation, but we will focus on R

- Available as [Free](https://en.wikipedia.org/wiki/Free_software_movement) / 
    	  [Open Source](https://en.wikipedia.org/wiki/The_Open_Source_Definition) Software

- Popular in both academia and industry

- Easy to try out on your own
--

- Python and Julia are other good alternatives

---

<!--

# Outline

- Installing R (hopefully most of you have already done this)

- Basics of using R

- Some examples

# Installing R

* R is most commonly used as a [REPL](https://en.wikipedia.org/wiki/Read-eval-print_loop) (Read-Eval-Print-Loop)

* When it is started, R Waits for user input

* Evaluates and prints result

* Waits for more input

* There are several different _interfaces_ to do this

* R itself works on many platforms (Windows, Mac, UNIX, Linux)

* Some interfaces are platform-specific, some work on most

* R and the interface may need to be installed separately

# Installing R

* Go to <https://cran.r-project.org/> (or choose a [mirror](https://cran.r-project.org/mirrors.html) first)

* Follow instructions depending on your platform (Windows for most of you)

* This will install R, as well as a default graphical interface on Windows and Mac

* I will recommend a different interface called
  [R Studio](https://www.rstudio.com/) that needs to be installed separately
  
* I personally use yet another interface called
  [ESS](https://ess.r-project.org/) with a general-purpose
  editor called [Emacs](https://www.gnu.org/software/emacs/)

-->

<!-- R tutorial based on XLISP-STAT tutorial

This introductory tutorial is modeled after the structure,
organization, and sequence of topics presented in Luke Tierney's
XLISP-STAT Tutorial to provide a similar introduction to the R
programming environment.

-->

# What is R?

* Powerful statistical computing and graphics environment

* Widely used for statistical analysis and data science

* Implementation of the *S* statistical programming language

* Basic data types allow collections of numbers (data sets) to be manipulated as complete objects

---

# Our goal

---

* Not a complete manual

* Introduce fundamentals of R

* Basic statistical computations 
	
	* Visualizations of univariate data
	
	* Linear regression 
	
	* Defining custom functions
	
--

* The only way to learn a computer program: __use it!__

* Execute the provided examples yourself

* Make changes and see what happens

---

# Installing R

* Go to <https://cran.isid.ac.in/> (or choose another [mirror](https://cran.r-project.org/mirrors.html))

* Follow instructions depending on your platform (Windows for most of you)
--

* This will install R, as well as a default graphical interface on Windows and Mac
--

* I will recommend a different interface called
  [R Studio](https://posit.co/) that needs to be installed separately
  
* Other options:

* [Emacs](https://www.gnu.org/software/emacs/) + [ESS](https://ess.r-project.org/)
	
	* [VS Code](https://code.visualstudio.com/docs/languages/r)
	
	* [Positron](https://positron.posit.co/)

---

# Starting

* Start R (usually through an interface like RStudio)

* When ready, the command window (console) will display a prompt, typically `>`.

* Example: type the number `1` and press *return* (or *Enter*)

* R will print the result and provide a new prompt:

```r
> 1
```
```
 1
```
```r
>
```

---

# Evaluating expressions

* We can also type an 'expression' (something that requires calculation)

* R will _evaluate_ the expression and print the result, followed by a new prompt

``` r
1 + 2
```

```
[1] 3
```
```r
>
```

---

# Finishing

* To exit R, type the command `q()` or `q("no")`

```r
> q()
```

```
Save workspace image? [y/n/c]: n
```

* In subsequent examples, we will omit the `>` prompts for brevity

---

# Notation for expressions

---

* Mathematical functions use two types of standard notation:

* Prefix notation

$$
	  f: \mathbb{R}^2 \to \mathbb{R} \quad\quad\quad f: (x, y) \mapsto f(x, y)
      $$
--

* Infix notation

$$
      +: \mathbb{R}^2 \to \mathbb{R} \quad\quad\quad +: (x, y) \mapsto x + y
      $$

* Also some nonstandard notations:

$$
  \text{power} : \mathbb{R}^2 \to \mathbb{R} \quad\quad\quad ?? : (x, p) \mapsto x^p
  $$

$$
  \sin^2 : \mathbb{R} \to [0, 1] \quad\quad\quad \sin^2: x \mapsto \sin^2 x = (\sin x)^2
  $$
  
  $$
  \log : [0, \infty)^2 \to \mathbb{R} \quad\quad\quad\quad \log: (x, b) \mapsto \log_b x = \frac{\ln x}{\ln b}
  $$

---

* R mostly uses **prefix** notation or "function" notation

``` r
sin(10)
```

```
[1] -0.5440211
```

``` r
log(10)     # natural logarithm
```

```
[1] 2.302585
```

``` r
log(10, 2)  # logarithm with base 2
```

```
[1] 3.321928
```

* The hash (`#`) is recognized as a comment character in R

* Anything following it on the same line is ignored

---

* R also **infix** notation for standard arithmetic operators

``` r
10 + 2
```

```
[1] 12
```

``` r
10^2
```

```
[1] 100
```

``` r
sin(10)^2
```

```
[1] 0.295959
```

---

* But the underlying evaluation is still based on prefix notation

* We can also use the _operators_ as if they were _functions_

``` r
`+`(10, 2)
```

```
[1] 12
```

``` r
`^`(10, 2)
```

```
[1] 100
```

``` r
`^`(sin(10), 2)
```

```
[1] 0.295959
```

---

# Basic principles: Listener and Evaluator (REPL)

---

* R is an _interactive_ system

* R has two basic components:

* A **listener** (console) that waits for instructions

* An **evaluator** (engine) that executes these instructions

---

* R can handle many different kinds of data

* Basic classification: *simple data* and *compound data*

---

## Examples of **Simple Data**

* Numbers (numeric values, including integers and 'floating-point' numbers)

```r
1     # an integer  
-3.14 # a floating point ('real') number
```
--

* Logical values

```r
TRUE  # true  
FALSE # false
      ```
---

## Examples of **Simple Data**

* Strings (enclosed in single or double quotes)

```r
"This is a string 1 2 3 4"
```

* Symbols (used for naming variables or functions)

```r
x
gdp.data
this.is.a.symbol
```

---

## Examples of **Compound Data**

* **vectors** (ordered collections of elements of the same type) 
	
* **lists** (ordered collections that can contain elements of different types or even other lists/vectors)

* We commonly define compound data using the concatenate function `c()`:

``` r
c(1, 2, 3)
```

```
[1] 1 2 3
```

---

* An R session is an interaction between the user and the console (listener)

* User inputs an expression and hits *return*, the listener passes it to the evaluator

* Basic rule:

* Everything is evaluated

* The results are (usually) printed, and the listener goes back to listening

* This is known as the **Read-Eval-Print Loop** or **REPL**.
--

* Numbers and strings evaluate to themselves:

``` r
1  
```

```
[1] 1
```

``` r
"Hello"  
```

```
[1] "Hello"
```

---

* Complex expressions typically use functions (and operators) in R

* Functions are applied using parentheses, mirroring prefix notation

``` r
sqrt(1 + 2)
```

```
[1] 1.732051
```

``` r
nchar("Hello")
```

```
[1] 5
```

---

* We can sometimes wish to tell the evaluator *not* to evaluate something

* This is typically when defining expressions symbolically

* This is accomplished using the `quote()` function

``` r
quote(1 + 2)  # 1 + 2 is not evaluated
```

```
1 + 2
```

---

# Elementary Statistical Operations

---

* Statistical data usually consists of groups of numbers

* Example from _World Handbook of Political And Social Indicators_ (Bruce M. Russett et al, 1964)

* Per capita GNP (around 1957)
	
	* Adult literacy (%) (around 1960)
	
	* Enrollment in higher education per 100,000 population (around 1960)
	
* Only a subset of countries are shown (with relatively low per capita GDP)

---

country|gnppc|pctlit_adult|highered100k
:------|----:|-----------:|-----------:
Nepal|45|5|56
Afghanistan|50|2.5|12
Laos|50|17.5|4
Ethiopia|55|2.5|5
Burma|57|47.5|63
Libya|60|13|49
Sudan|60|9|34
Tanganyika|61|7.5|9
Uganda|64|27.5|14
Pakistan|70|13|165
China|73|47.5|69
India|73|19.3|220
South Vietnam|76|17.5|83
Nigeria|78|10|4
Kenya|87|22.5|5
Madagascar|88|33.5|21
Congo|92|37.5|4
Thailand|96|68|251
Bolivia|99|32.1|166
Cambodia|99|17.5|18

]

---

* For now, let us consider only the data on enrolment in higher education:

```
56 12 4 5 63 49 34 9 14 165 69 220 83 4 5 21 4 251 166 18
```

* In R, we represent this data as a vector using `c()` (combine):

``` r
c(56, 12, 4, 5, 63, 49, 34, 9, 14, 165, 69, 220, 83, 4, 5, 21, 4, 251, 166, 18)
```

```
 [1]  56  12   4   5  63  49  34   9  14 165  69 220  83   4   5  21   4 251 166  18
```

---

* The `mean()` function computes the **average** of a vector of numbers:

``` r
mean(c(56, 12, 4, 5, 63, 49, 34, 9, 14, 165, 69, 220, 83, 4, 5,
       21, 4, 251, 166, 18))
```

```
[1] 62.6
```

* The **median** of these numbers can be calculated using `median()`:

``` r
median(c(56, 12, 4, 5, 63, 49, 34, 9, 14, 165, 69, 220, 83, 4, 5,
         21, 4, 251, 166, 18))
```

```
[1] 27.5
```

---

* To avoid retyping the data, we can assign it a name

* This is done using the assignment operator `<-` or the (mostly equivalent) `=` operator:

``` r
higher.educ <- c(56, 12, 4, 5, 63, 49, 34, 9, 14, 165, 69, 220, 83, 4,
                 5, 21, 4, 251, 166, 18)
```

* Now, the symbol `higher.educ` holds the vector of 20 numbers

* If we type the symbol, R returns its value

``` r
higher.educ
```

```
 [1]  56  12   4   5  63  49  34   9  14 165  69 220  83   4   5  21   4 251 166  18
```

---

* We can easily compute numerical descriptive statistics:

``` r
mean(higher.educ)  
```

```
[1] 62.6
```

``` r
median(higher.educ)  
```

```
[1] 27.5
```

``` r
sd(higher.educ) # Standard deviation  
```

```
[1] 76.57222
```

``` r
IQR(higher.educ) # Interquartile range
```

```
[1] 64.5
```

---

* R also supports **elementwise arithmetic operations** on vectors

* We can add 1 to each value as follows:

``` r
1 + higher.educ
```

```
 [1]  57  13   5   6  64  50  35  10  15 166  70 221  84   5   6  22   5 252 167  19
```

And calculate the natural logarithms of the results:

``` r
log(1 + higher.educ)
```

```
 [1] 4.043051 2.564949 1.609438 1.791759 4.158883 3.912023 3.555348 2.302585 2.708050 5.111988
[11] 4.248495 5.398163 4.430817 1.609438 1.791759 3.091042 1.609438 5.529429 5.117994 2.944439
```

---

# Summary Statistics and Plots

---

* Another example: average monthly PM 2.5 levels

* Recorded at an air quality monitoring station in R.K.Puram (a
locality in Delhi close to ISI)

* Over a 3-year period, from January 2021 to December 2023

``` r
pm25 <- c(288, 223, 167, 156, 126, 120, 102, 106, 83, 114, 259, 282, 
          234, 183, 174, 176, 160, 139, 102, 99, 110, 173, 245, 250, 260, 
          190, 150, 164, 161, 144, 115, 138, 123, 182, 323, 280)
```

---

* The function `hist()` can be used to obtain a histogram

``` r
hist(pm25) # Produces a histogram plot
```

![plot of chunk unnamed-chunk-19](figures/roverview-0-unnamed-chunk-19-1.svg)

---

* Here are some numerical summaries.

``` r
mean(pm25)  
```

```
[1] 175.0278
```

``` r
median(pm25)  
```

```
[1] 162.5
```

``` r
sd(pm25)  
```

```
[1] 63.83796
```

``` r
IQR(pm25)
```

```
[1] 103.5
```

---

* We can look at transformations like square root and log

``` r
mean(sqrt(pm25))  
```

```
[1] 13.02359
```

``` r
median(sqrt(pm25))  
```

```
[1] 12.74741
```

``` r
mean(log(pm25))  
```

```
[1] 5.102077
```

``` r
median(log(pm25))  
```

```
[1] 5.090635
```

---

* The `boxplot()` function produces another useful summary in terms of quartiles

* We can compare the boxplots of original and transformed observations to see how they affect symmetry

``` r
par(mfrow = c(1, 3))
boxplot(pm25, main = "PM25"); boxplot(sqrt(pm25), main = "sqrt(PM25)"); boxplot(log(pm25), main = "log(PM25)")
```

![plot of chunk unnamed-chunk-22](figures/roverview-0-unnamed-chunk-22-1.svg)

---

# Scatter Plots

---

* Particularly important for us

* Shows relationship between two variables

* In R, the `plot()` function creates a scatterplot
--

* For the PM2.5 data, we may want to plot PM2.5 levels against time

* Useful to check for serial correlation or trend

* This is known as a time series plot
--

* For this, we need a sequence of integers for the time (X-axis) variable

* The function `seq()` or the shorthand `:` operator generates a list
  of consecutive integers.

---

``` r
time <- 0:35  
plot(time, pm25)
```

![plot of chunk unnamed-chunk-23](figures/roverview-0-unnamed-chunk-23-1.svg)

---

* There appears to be a clear trend

* For time series plots such as this, it is common to connect points by lines

* We can do so using the `type` argument:

* `type = "l"` (for lines only) 
	
	* `type = "o"` for points connected by lines.

---

``` r
plot(time, pm25, type = "o") # "o" stands for 'overlay'
```

![plot of chunk unnamed-chunk-24](figures/roverview-0-unnamed-chunk-24-1.svg)

---

* Scatterplots are highly useful for examining the relationship between two numerical variables

* Going back to the table of social indicators data, let us also consider adult literacy (%)

``` r
adult.lit <- c(5, 2.5, 17.5, 2.5, 47.5, 13, 9, 7.5, 27.5, 13, 47.5,
               19.3, 17.5, 10, 22.5, 33.5, 37.5, 68, 32.1, 17.5)
```

* We can then create a scatter plot of `higher.educ` against `adult.lit` using

``` r
plot(adult.lit, higher.educ)
```

---

![plot of chunk unnamed-chunk-27](figures/roverview-0-unnamed-chunk-27-1.svg)

---

# Plotting Functions

---

* The `plot()` function can also be used to graph mathematical functions

* Suppose we want to plot $\sin(x)$ from $-\pi$ to $+\pi$

* R predefines $\pi$ as `pi` (upto available precision, of course)

* So we can construct a sequence of equally spaced real numbers as follows:

``` r
x_points <- seq(-pi, pi, length.out = 50)
x_points
```

```
 [1] -3.14159265 -3.01336438 -2.88513611 -2.75690784 -2.62867957 -2.50045130 -2.37222302 -2.24399475
 [9] -2.11576648 -1.98753821 -1.85930994 -1.73108167 -1.60285339 -1.47462512 -1.34639685 -1.21816858
[17] -1.08994031 -0.96171204 -0.83348377 -0.70525549 -0.57702722 -0.44879895 -0.32057068 -0.19234241
[25] -0.06411414  0.06411414  0.19234241  0.32057068  0.44879895  0.57702722  0.70525549  0.83348377
[33]  0.96171204  1.08994031  1.21816858  1.34639685  1.47462512  1.60285339  1.73108167  1.85930994
[41]  1.98753821  2.11576648  2.24399475  2.37222302  2.50045130  2.62867957  2.75690784  2.88513611
[49]  3.01336438  3.14159265
```

---

``` r
plot(x_points, sin(x_points), type = "l")
```

![plot of chunk unnamed-chunk-29](figures/roverview-0-unnamed-chunk-29-1.svg)

---

* It is also possible to plot functions (of one argument) directly

* For an existing function like the `sin` function above, this very simple

``` r
plot(sin, from = -2 * pi, to = 2 * pi)
```

![plot of chunk unnamed-chunk-30](figures/roverview-0-unnamed-chunk-30-1.svg)

---

* We can alse define a new function to plot as follows

``` r
f <- function(x) { 2 * x + 3 * x^2 - x^3 }
plot(f, from = -10, to = 10)
```

![plot of chunk unnamed-chunk-31](figures/roverview-0-unnamed-chunk-31-1.svg)

---

# For loops and mapping

---

* A more interesting example: _loss functions_

* Let <span>$X_1, X_2, \dotsc, X_n$</span> represent observed data
  values

* Suppose we want to find $\theta$ which minimizes
  $$
  L(\theta) = \sum\limits_{i=1}^n (X_i - \theta)^2
  $$

* This is the sum of squared deviations of a every value in a dataset
  from $\theta$
  
* Suppose we want to calculate and plot $L(\theta)$

---

* Let us use the `higher.educ` data as our <span>$X_i$</span> values

* We can implement this function as follows

``` r
SSD <- function(theta) {
    dev <- higher.educ - theta
    sum(dev * dev) # value returned by function
}
```

---

* We can now try to use the earlier approach

* Evaluate this function on a generated set of input values
	
	* Plot them using `plot(..., type = "l")`
--
	
* One problem with this approach is that `SSD()` is not a vectorized function

* It can only calculate the answer for one scalar argument `theta` at a time
	
	* As a result, the following will fail with an error:

``` r
theta_vals <- seq(50, 70, length.out = 100)
plot(theta_vals, SSD(theta_vals), type = "l")
```

```
Error in xy.coords(x, y, xlabel, ylabel, log): 'x' and 'y' lengths differ
```

---

* In situations like this, we need to evaluate the answer element-wise using a for loop

* Common approach:

* Create a placeholder vector to hold the evaluated function values
	
	* Fill them in one by one in a loop

``` r
y_vals <- theta_vals # placeholder of same length as input
for (i in 1:100) y_vals[[i]] <- SSD(theta_vals[[i]])
```

* This needs one additional piece of knowledge:

* `x[[i]]` gives us the `i`-th component of a compound data vector `x`

* `x[[i]]  <- value` assigns to the `i`-th component of a compound data vector `x`

---

``` r
plot(theta_vals, y_vals, type = "l")
```

![plot of chunk unnamed-chunk-35](figures/roverview-0-unnamed-chunk-35-1.svg)

---

* Recall: the SSD function is minimized at the sample mean of a dataset

``` r
mean(higher.educ)
```

```
[1] 62.6
```

---

* For loops are generally useful, but unnecessarily complicated in this situation

* Special case: apply or map a function to all elements of a vector (of input values)

* Simpler solution: use the `sapply()` function.

---

``` r
y_vals2 <- sapply(theta_vals, SSD) # evaluates SSD(x) for each theta in theta_vals
str(y_vals2)
```

```
 num [1:100] 114578 114477 114378 114280 114184 ...
```

``` r
str(y_vals) # earlier calculation using for loop
```

```
 num [1:100] 114578 114477 114378 114280 114184 ...
```

---

``` r
plot(theta_vals, y_vals2, type = "l")
```

![plot of chunk unnamed-chunk-38](figures/roverview-0-unnamed-chunk-38-1.svg)

---

* Another option: convert "scalar" function into a vectorized function

``` r
SSD_vec <- Vectorize(SSD)
plot(SSD_vec, from = 60, to = 65)
```

![plot of chunk unnamed-chunk-39](figures/roverview-0-unnamed-chunk-39-1.svg)

---

# Generating Random and Systematic Data

---

* R provides functions for generating pseudo-random numbers

* E.g., `runif(n)` for Uniform, `rnorm(n)` for standard Normal

``` r
runif(10)
```

```
 [1] 0.20497050 0.07698767 0.68109903 0.18578231 0.47300244 0.78644715 0.15833495 0.09038630
 [9] 0.22583684 0.50226505
```

``` r
rnorm(10)
```

```
 [1] -1.0334143  0.5231131 -1.4147663  0.7842986  2.0773395 -0.7179434  0.1439253 -0.1252617
 [9] -0.1378383  0.7910844
```

* `rpois(n, lambda)` for Poisson with mean `lambda`

``` r
rpois(30, lambda = 5)
```

```
 [1]  9  3  3 11  7  4  8  7  4  7  7  5  5  3  8  8  5  1  3  7  4  5  6  4  3  6  4  2  1  1
```

---

* We have already used

* `seq(start, end)` (or `start:end`) for equally spaced sequences of integers

* `seq(a, b, length.out = n)` for equally spaced sequences of real numbers
--

* The `rep()` function is useful for generating specific patterns

* Repeat a sequence:

``` r
rep(c(1, 2, 3), 2)  
```

```
[1] 1 2 3 1 2 3
```
--

* Repeat each element a specified number of times:

``` r
rep(c(1, 2, 3), times = c(3, 2, 1))  
```

```
[1] 1 1 1 2 2 3
```

---

# Forming Subsets and Deleting Cases

---

* R uses bracket indexing `[]` to select elements from a vector or list

* **An important point to note**,  R uses 1-based indexing (not 0-based indexing)

* Define a vector `x`:

``` r
x <- c(3, 7, 5, 9, 12, 3, 14, 2)
```

* To retrieve the 1st element (index 1), we can use either of the following:

``` r
x[[5]]
```

```
[1] 12
```

``` r
x[5]
```

```
[1] 12
```

---

* To retrieve a group of elements, we must use the second form

* The "index" (inside the `[ ... ]`) is in general a vector

``` r
x[c(1, 3)]
```

```
[1] 3 5
```

---

* To exclude elements, we use negative indices

* Index can be vector or scalar

* To exclude the 5th element:

``` r
x[-5]
```

```
[1]  3  7  5  9  3 14  2
```

* To exclude multiple elements:

``` r
x[-c(1, 3, 7)]
```

```
[1]  7  9 12  3  2
```

---

* **Logical indexing**: The index can also be a logical vector

* E.g., select all elements of `x` that are greater than 3:

``` r
x[x > 3]
```

```
[1]  7  5  9 12 14
```
--

* This is actually a two-step process:

``` r
index <- x > 3
index
```

```
[1] FALSE  TRUE  TRUE  TRUE  TRUE FALSE  TRUE FALSE
```

``` r
x[index]
```

```
[1]  7  5  9 12 14
```

---

# Combining Several Lists

* Use `c()` to **combine** several short vectors into a single longer vector

``` r
z1 <- c(1, 2, 3)  
z2 <- c(4)
z3 <- c(5, 6, 7, 8)  
c(z1, z2, z3)  
```

```
[1] 1 2 3 4 5 6 7 8
```

---

# Modifying Data

---

* Goal: Replace values in a vector after it has been defined

* Standard paradigm: combine subsetting with assignment

* Suppose we want to change the the 5th element in `x` to 11:

``` r
x[[5]] <- 11  
x  
```

```
[1]  3  7  5  9 11  3 14  2
```

* The same idea works for vector indexing as well

* To change elements 1 and 3 to 15 and 16:

``` r
x[c(1, 3)] <- c(15, 16)  
x  
```

```
[1] 15  7 16  9 11  3 14  2
```

---

* Referencing versus copying: behavior of R is different from Python

* R makes "deep copies" of vectors when copying them

* This means modifying one vector does not modify other copies of the vector

``` r
x
```

```
[1] 15  7 16  9 11  3 14  2
```

``` r
y <- x # y is a copy  
x[[3]] <- 100  
x
```

```
[1]  15   7 100   9  11   3  14   2
```

``` r
y
```

```
[1] 15  7 16  9 11  3 14  2
```

---

# Getting Help

---

* Online help is available for most R functions

* You can use the `?` operator followed by the function name, or the `help()` function:

```r
?median  
help("median")
```

* The `??` operator searches the documentation for keywords.

```r
??normal
```

---

# Listing and Undefining Variables

---

* To find out which variables we have defined in the current session:

``` r
ls()
```

```
 [1] "adult.lit"   "f"           "higher.educ" "i"           "index"       "ndata"      
 [7] "plot"        "pm25"        "showCall"    "SSD"         "SSD_vec"     "theta_vals" 
[13] "time"        "x"           "x_points"    "y"           "y_vals"      "y_vals2"    
[19] "z1"          "z2"          "z3"         
```

* Sometimes useful to remove a variable to free up memory or clean up workspace

* use the `rm()` function:

``` r
rm(theta_vals, y_vals, y_vals2, SSD, SSD_vec)
ls()
```

```
 [1] "adult.lit"   "f"           "higher.educ" "i"           "index"       "ndata"      
 [7] "plot"        "pm25"        "showCall"    "time"        "x"           "x_points"   
[13] "y"           "z1"          "z2"          "z3"         
```

---

# Saving Your Work

---

* R provides mechanisms to save variables for later use:

``` r
save(higher.educ, pm25, file = "examples.rda")
```

* This saves the specified variables to a file in a special binary format

* Can be reloaded later in a different R session.

```r
load("examples.rda")
```

---

* To save your code, write code in a **file** (instead of the R prompt) and save it for later use

* RStudio and most IDEs makes this very easy

---

# Loading files

---

* Saved data files can be loaded using the `load()` function

* Similarly, R scripts can be run using the `source()` function

* This runs a previously saved R script as a series of commands:

```r
source("/path/to/script.R")
```

---

# Importing data stored in other formats

---

* Small datasets can be typed in at the R console to illustrate basic usage

* This is not practical for real world datasets which can be large

* `load()` can be used, but only when they have already been saved using R

---

* Commonly used text-based formats

* CSV (comma-separated values)
	
	* TSV (TAB-separated values)
	
	* Other simple text formats

* Easy to read directly from R using `read.csv()`, `read.table()`, etc.

---

* But other prorietary (binary) formats are also common

* `.xls` or `.xlsx` files exported by Microsoft Excel

* `.xpt` files exported by SAS

* `.sav` files exported by SPSS

* `.dta` files exported by Stata.

* Usually possible using additional packages

* See the [R Data Import/Export](https://cran.icts.res.in/doc/manuals/r-devel/R-data.html) manual

---

# Simple Linear Regression

---

* Linear regression models are usually fit using the function `lm()` (linear model)

* `lm()` fits the model, but not necessarily everything else you may want

* Usual pattern:

* Save the result of `lm()` in a variable

* Call specific functions to get coefficients, residuals, predictions, tests, etc.

---

* Example: social indicators data on `higher.educ` vs `adult.lit`

``` r
plot(adult.lit, higher.educ)
```

![plot of chunk unnamed-chunk-58](figures/roverview-0-unnamed-chunk-58-1.svg)

---

* Transform data to reduce asymmetry square roots are often useful for proportions)

``` r
plot(sqrt(adult.lit), log(higher.educ))
```

![plot of chunk unnamed-chunk-59](figures/roverview-0-unnamed-chunk-59-1.svg)

---

* Is ordinary least squares resonable for this data?

``` r
fm <- lm(log(higher.educ) ~ 1 + sqrt(adult.lit))
fm
```

```

Call:
lm(formula = log(higher.educ) ~ 1 + sqrt(adult.lit))

Coefficients:
    (Intercept)  sqrt(adult.lit)  
         2.0235           0.2901  
```

---

* The `lm()` function returns a **model object**

* The fitted coefficients can be obtained using

``` r
coef(fm)
```

```
    (Intercept) sqrt(adult.lit) 
      2.0234664       0.2900805 
```

---

* These coefficients can be used to add the regression line to the scatterplot

``` r
plot(sqrt(adult.lit), log(higher.educ))
abline(coef(fm))
```

![plot of chunk unnamed-chunk-62](figures/roverview-0-unnamed-chunk-62-1.svg)

---

* A summary of the fitted model with some basic tests

``` r
summary(fm)
```

```

Call:
lm(formula = log(higher.educ) ~ 1 + sqrt(adult.lit))

Residuals:
    Min      1Q  Median      3Q     Max 
-2.4135 -0.8809  0.0616  1.1279  2.0958

Coefficients:
                Estimate Std. Error t value Pr(>|t|)  
(Intercept)       2.0235     0.8185   2.472   0.0236 *
sqrt(adult.lit)   0.2901     0.1725   1.682   0.1099  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.367 on 18 degrees of freedom
Multiple R-squared:  0.1358,	Adjusted R-squared:  0.08779 
F-statistic: 2.829 on 1 and 18 DF,  p-value: 0.1099
```

---

* Output of `summary()` includes

* Coefficient estimates along with standard errors + $t$-tests
	
	* $F$-test for the **overall** model (compared to **mean-only** / intercept model)
	
	* These are same as there is only one linear term in addition to the intercept
--

* Can we conclude that adult literacy is independent of higher education enrolment in the population?

???

This lack of significance should not be taken too seriously, as the
regression model is based on very few data points, and small samples
rarely provide evidence for significant associations.  A larger
dataset consisting of records from 108 countries is available
[here](data/social-indicators-1964.csv), which we will consider in
another tutorial.

---

# Defining Our Own Functions

---

* A basic tool in any programming language

* In R, functions are defined using the `function()` keyword

---

* A function to compute the mean of a vector (similar to the built-in `mean()` function)

``` r
my.mean <- function(x) {
    s <- sum(x)
    n <- length(x)
	s / n
}
my.mean(runif(50))
```

```
[1] 0.4636072
```

---

* A function to delete a case from a vector, taking the vector and the 1-based index `i`

``` r
delete.case <- function(x, i) {  
    x[-i]  
 }
```
--

* Most languages have a `return` keyword, which is optional in R; this is equivalent:

``` r
delete.case <- function(x, i) {  
    return(x[-i])
}
```

* If there is no call to return, the last evaluated value is returned

---

# Anonymous Functions

---

* Functions need not be **assigned** to a variable

* The following defines a function that does not have a "name"

```r
function(n) {
    sample(1:500, n, replace = TRUE)
}
```

* This function takes a SRSWR from a population of size 500
--

* R treats functions as any other kind of data

* In particular, anonymous functions may be used as arguments wherever a function is expected

---

* Example: plot the Standard Normal density using `plot()`

``` r
plot(function(x) exp(-x^2 / 2) / sqrt(2 * pi), from = -3, to = 3)
```

![plot of chunk unnamed-chunk-67](figures/roverview-0-unnamed-chunk-67-1.svg)

---

* Mapping functions like `sapply` or `lapply` frequently use anonymous functions as arguments

* Example: take a SRSWR of various sizes from a population of size 500:

``` r
lapply(c(7, 13, 17), function(n) sample(1:500, n, replace = TRUE))
```

```
[[1]]
[1] 130 135 239 453  14 226 442

[[2]]
 [1]  57 447 354 232 286 327 110 469 432   8  29 483 181

[[3]]
 [1]  95 353 143 261  44  39 393 279 325 364 333 315  68  37 417 157 125
```

---

* Here, the expression `function(n) sample(1:500, n, replace = TRUE)` is an anonymous function

* Of course, we could name it by assigning it to a variable

``` r
srswr500 <- function(n) sample(1:500, n, replace = TRUE)
```

* and then call `lapply()` using the variable name to refer to the function

``` r
lapply(c(7, 13, 17), srswr500)
```

```
[[1]]
[1] 314   1 148 311 189  78 456

[[2]]
 [1] 447   9 386 447 160 265 268 400 429 159 423 357 262

[[3]]
 [1] 432 413 448 125 200 371 312 187 279  51  45 212 195 104  53 423  38
```

???

Note that as this is simulating a "random" experiment, the results
will be different each time the code is run.

---

# Data frames

---

* So far we have been using simple vectors as datasets

* This is impractical for any serious analysis

* Real data are typically multivariate datasets, containing several measurements on each unit
--

* This is the type of data we typically see in a spreadsheet

* Such datasets are typically stored as **data frames** in R

---

* Data frames are usually created by data import functions such as `read.csv()`

* Can also be constructed using the `data.frame()` function

* Example: we can combine the variables of the social indicators data

``` r
dsocial <-
    data.frame(country = c("Nepal", "Afghanistan", "Laos", "Ethiopia",
                           "Burma", "Libya", "Sudan", "Tanganyika",
                           "Uganda", "Pakistan", "China", "India",
                           "South Vietnam", "Nigeria", "Kenya", "Madagascar",
                           "Congo", "Thailand", "Bolivia", "Cambodia"),
               hedu = higher.educ,
               adlit = adult.lit)
```

---

* This new variable, `dsocial` is now a matrix-like object

``` r
dsocial
```

```
         country hedu adlit
1          Nepal   56   5.0
2    Afghanistan   12   2.5
3           Laos    4  17.5
4       Ethiopia    5   2.5
5          Burma   63  47.5
6          Libya   49  13.0
7          Sudan   34   9.0
8     Tanganyika    9   7.5
9         Uganda   14  27.5
10      Pakistan  165  13.0
11         China   69  47.5
12         India  220  19.3
13 South Vietnam   83  17.5
14       Nigeria    4  10.0
15         Kenya    5  22.5
16    Madagascar   21  33.5
17         Congo    4  37.5
18      Thailand  251  68.0
19       Bolivia  166  32.1
20      Cambodia   18  17.5
```

---

* Individual columns can be extracted using the `$` operator

* So we can still use the methods we have learned earlier

``` r
mean(dsocial$adlit)
```

```
[1] 22.52
```

``` r
lm(log(dsocial$hedu) ~ 1 + sqrt(dsocial$adlit))
```

```

Call:
lm(formula = log(dsocial$hedu) ~ 1 + sqrt(dsocial$adlit))

Coefficients:
        (Intercept)  sqrt(dsocial$adlit)  
             2.0235               0.2901  
```

---

* However, this is **not** the recommended way

* Most modeling and visualization functions in R are designed to work with data frames as input

* They have special features to make this easy

* A common approach: the **formula interface**

---

* We have already seen this with `lm()`

``` r
fm.social <- lm(log(hedu) ~ sqrt(adlit), data = dsocial)
summary(fm.social)
```

```

Call:
lm(formula = log(hedu) ~ sqrt(adlit), data = dsocial)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.4135 -0.8809  0.0616  1.1279  2.0958

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)   2.0235     0.8185   2.472   0.0236 *
sqrt(adlit)   0.2901     0.1725   1.682   0.1099  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.367 on 18 degrees of freedom
Multiple R-squared:  0.1358,	Adjusted R-squared:  0.08779 
F-statistic: 2.829 on 1 and 18 DF,  p-value: 0.1099
```

---

* The `plot()` function also has a similar formula interface

``` r
plot(log(hedu) ~ sqrt(adlit), data = dsocial)
abline(coef(fm.social))
```

![plot of chunk unnamed-chunk-75](figures/roverview-0-unnamed-chunk-75-1.svg)

---

* Alternative visualization packages: **lattice** and **ggplot2**

``` r
library(lattice)
xyplot(log(hedu) ~ sqrt(adlit), data = dsocial, grid = TRUE, smooth = "lm")
```

![plot of chunk unnamed-chunk-76](figures/roverview-0-unnamed-chunk-76-1.svg)

---

``` r
library(ggplot2)
ggplot(dsocial, mapping = aes(x = sqrt(adlit), y = log(hedu))) +
    geom_point() + geom_smooth(method = "lm", se = FALSE)
```

![plot of chunk unnamed-chunk-77](figures/roverview-0-unnamed-chunk-77-1.svg)

---

# Example: Full Social Indicators dataset (107 countries)

---

``` r
dsocial1964 <-
    read.csv("https://www.isid.ac.in/~deepayan/RT2026/data/social-indicators-1964.csv",
             comment.char = "#")
str(dsocial1964)
```

```
'data.frame':	107 obs. of  7 variables:
 $ Country                  : chr  "Nepal" "Afghanistan" "Laos" "Togo" ...
 $ GNP.per.Capita           : int  45 50 50 50 55 57 60 60 60 61 ...
 $ Percent.Urban            : num  4.4 7.5 4 4.5 1.7 10 4.7 18.4 5 3.3 ...
 $ Percent.Adult.Literacy   : num  5 2.5 17.5 7.5 2.5 47.5 2.5 13 9 7.5 ...
 $ Higher.Ed.per.100000     : num  56 12 4 NA 5 63 NA 49 34 9 ...
 $ Inhabitants.per.Physician: int  72000 41000 100000 58000 117000 15000 14000 5800 40000 18000 ...
 $ Radios.per.1000          : num  NA 1.7 8 4.3 4.5 5.6 11.4 63.1 0.9 3.6 ...
```

<!--

$ Country                  : chr  "Nepal" "Afghanistan" "Laos" "Togo" ...
 $ GNP.per.Capita           : int  45 50 50 50 55 57 60 60 60 61 ...
 $ Percent.Urban            : num  4.4 7.5 4 4.5 1.7 10 4.7 18.4 5 3.3 ...
 $ Percent.Adult.Literacy   : num  5 2.5 17.5 7.5 2.5 47.5 2.5 13 9 7.5 ...
 $ Higher.Ed.per.100000     : num  56 12 4 NA 5 63 NA 49 34 9 ...
 $ Inhabitants.per.Physician: int  72000 41000 100000 58000 117000 15000 14000 5800 40000 18000 ...
 $ Radios.per.1000          : num  NA 1.7 8 4.3 4.5 5.6 11.4 63.1 0.9 3.6 ...

-->

---

``` r
plot(log(Higher.Ed.per.100000) ~ sqrt(Percent.Adult.Literacy), data = dsocial1964)
abline(coef(lm(log(Higher.Ed.per.100000) ~ sqrt(Percent.Adult.Literacy), data = dsocial1964)))
```

![plot of chunk unnamed-chunk-79](figures/roverview-0-unnamed-chunk-79-1.svg)

---

``` r
xyplot(log(Higher.Ed.per.100000) ~ sqrt(Percent.Adult.Literacy), dsocial1964, grid = TRUE, smooth = "lm")
```

![plot of chunk unnamed-chunk-80](figures/roverview-0-unnamed-chunk-80-1.svg)

---

``` r
ggplot(dsocial1964, mapping = aes(x = sqrt(Percent.Adult.Literacy), y = log(Higher.Ed.per.100000))) +
    geom_point() + geom_smooth(method = "lm", se = FALSE)
```

![plot of chunk unnamed-chunk-81](figures/roverview-0-unnamed-chunk-81-1.svg)

---

``` r
fm1964 <- lm(log(Higher.Ed.per.100000) ~ sqrt(Percent.Adult.Literacy), data = dsocial1964)
coef(fm1964)
```

```
                 (Intercept) sqrt(Percent.Adult.Literacy) 
                   1.9750857                    0.4020052 
```

---

``` r
summary(fm1964)
```

```

Call:
lm(formula = log(Higher.Ed.per.100000) ~ sqrt(Percent.Adult.Literacy), 
    data = dsocial1964)

Residuals:
    Min      1Q  Median      3Q     Max 
-5.4175 -0.3814  0.2343  0.7025  2.2206

Coefficients:
                             Estimate Std. Error t value Pr(>|t|)    
(Intercept)                   1.97509    0.36691   5.383 4.93e-07 ***
sqrt(Percent.Adult.Literacy)  0.40201    0.04812   8.354 4.16e-13 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.19 on 99 degrees of freedom
  (6 observations deleted due to missingness)
Multiple R-squared:  0.4135,	Adjusted R-squared:  0.4076 
F-statistic: 69.79 on 1 and 99 DF,  p-value: 4.158e-13
```

---

# Questions?

<!--

* Interactive plot using the **plotly** package.

```r
library(plotly)
ggplotly(p)
```

-->