## Some help for R¶

In this short notebook, I will go through a few basic examples in R that you may find useful for the course.

These are just some of the things I find useful. Feel free to search around for others.

For those of you who have done some programming before, you will notice that R is a functional programming language.

### Functions¶

You might get tired of always typing http://stats191.stanford.edu/data. You could make a small function

In [1]:
useful_function = function(dataname) {
return(paste("http://stats191.stanford.edu/data/", dataname, sep=''))
}

useful_function("groundhog.table")

'http://stats191.stanford.edu/data/groundhog.table'

Let's load the heights data with less code

In [2]:
h.table = read.table(useful_function("groundhog.table"), header=TRUE, sep=',')

199024 N
199123 Y
199222 Y
199316 Y
199412 Y
199513 N

Or, for all data sets in the course directory, we might try

In [3]:
course_dataset = function(dataname, sep='', header=TRUE) {
}

199024 N
199123 Y
199222 Y
199316 Y
199412 Y
199513 N

Note that I didn't use return in the function above. By default, R will return the object in the last line of the function code.

In [4]:
test_func = function(x) {
x^2
3
}
test_func(4)

3

### Source¶

When working on a particular project or assignment, it is often easiest to type commands in a text editor and rerun them several times. The command source is an easy way to do this, and it takes either the name of a file or a URL as argument. Suppose we have a webpage http://stats191.stanford.edu/R/helper_code.R

Then, we can execute this as follows

In [5]:
source("http://stats191.stanford.edu/R/helper_code.R")

199024 N
199123 Y
199222 Y
199316 Y
199412 Y
199513 N

As you go through the course, you might copy this file to a your computer and add some other useful functions to this file.

For larger collections of functions, R allows the creation of packages that can be installed and loaded with a call to the library function. Documentation on packages can be found here.

### Concatenation, sequences¶

Many tasks involving sequences of numbers. Here are some basic examples on how to manipulate and create sequences.

The function c, concatenation, is used often in R, as are rep and seq

In [6]:
X = 3
Y = 4
c(X,Y)

1. 3
2. 4

The function rep denotes repeat.

In [7]:
print(rep(1,4))
print(rep(2,3))
c(rep(1,4), rep(2,3))

[1] 1 1 1 1
[1] 2 2 2

1. 1
2. 1
3. 1
4. 1
5. 2
6. 2
7. 2

The function seq denotes sequence. There are various ways of specifying the sequence.

In [8]:
seq(0,10,length=11)

1. 0
2. 1
3. 2
4. 3
5. 4
6. 5
7. 6
8. 7
9. 8
10. 9
11. 10
In [9]:
seq(0,10,by=2)

1. 0
2. 2
3. 4
4. 6
5. 8
6. 10

You can sort and order sequences

In [10]:
X = c(4,6,2,9)
sort(X)

1. 2
2. 4
3. 6
4. 9

Use an ordering of X to sort a list of Y in the same order

In [11]:
Y = c(1,2,3,4)
o = order(X)
X[o]
Y[o]

1. 2
2. 4
3. 6
4. 9
1. 3
2. 1
3. 2
4. 4

A word of caution. In R you can overwrite builtin functions so try not to call variables c:

In [12]:
c = 3
c

3

However, this has not overwritten the function c.

In [13]:
c(3,4,5)

1. 3
2. 4
3. 5

Other variables to be careful are the aliases T for TRUE and F for FALSE. Since we compute $t$ and $F$ statistics it is natural to also have variables named T so when you are expecting T to be TRUE you might get a surprise.

In [14]:
c(T,F)

1. TRUE
2. FALSE

For other style advice, try reading Hadley Wickham's style guide. This is part of a fairly extensive online book. Google also has a style guide.

### Indexing¶

Often times, we will want to extract a subset of rows (or columns) of a vector (or matrix). R supports using logical vectors as index objects.

In [15]:
X = c(4,5,3,6,7,9)
Y = c(4,2,65,3,5,9)
X[Y>=5]

1. 3
2. 7
3. 9

Suppose we have a data.frame and want to extract from rows or columns. Rows are the first of two indexing objects while columns correspond to the second indexing object. Suppose we want to find take the mother and daughter heights where the daughter's height is less than or equal to 62 inches. Note the "," in the square brackets below: this tells R that it is looking for a subset of rows of the data.frame.

In [16]:
library(alr3)
data(heights)
subset_heights = heights[heights$Dheight <= 62,] print(c(nrow(heights), nrow(subset_heights)))  Loading required package: car Warning message: “package ‘car’ was built under R version 3.3.2” MheightDheight 59.755.1 58.256.5 60.656.0 60.756.8 61.856.0 55.557.9 [1] 1375 354  ### Plotting¶ R has a very rich plotting library. Most of our plots will be fairly straightforward, "scatter plots". In [17]: X = c(1:40) Y = 2 + 3 * X + rnorm(40) * 10 plot(X, Y)  The plots can be made nicer by adding colors and using different symbols. See the help for function par. In [18]: plot(X, Y, pch=21, bg='red')  You can add titles, as well as change the axis labels. In [19]: plot(X, Y, pch=23, bg='red', main='A simulated data set', xlab='Predictor', ylab='Outcome')  Lines are added with abline. We'll add some lines to our previous plot: a yellow line with intercept 2, slope 3, width 3, type 2, as well as a vertical line at x=20 and horizontal line at y=60. In [20]: plot(X, Y, pch=23, bg='red', main='A simulated data set', xlab='Predictor', ylab='Outcome') abline(2, 3, lwd=3, lty=2, col='yellow') abline(h=60, col='green') abline(v=20, col='red')  You can add points and lines to existing plots. In this example, we plot the first 20 points in red in one call to plot, then add the rest in blue with an orange line connecting them. In [21]: plot(X[1:20], Y[1:20], pch=21, bg='red', xlim=c(min(X),max(X)), ylim=c(min(Y),max(Y))) points(X[21:40], Y[21:40], pch=21, bg='blue') lines(X[21:40], Y[21:40], lwd=2, lty=3, col='orange')  You can put more than one plot on each device. Here we create a 2-by-1 grid of plots In [22]: par(mfrow=c(2,1)) plot(X, Y, pch=21, bg='red') plot(Y, X, pch=23, bg='blue') par(mfrow=c(1,1))  ### Saving plots¶ Plots can be saved as pdf, png, jpg among other formats. Let's save a plot in a file called "myplot.jpg" In [23]: jpeg("myplot.jpg") plot(X, Y, pch=21, bg='red') dev.off()  pdf: 2 Several plots can be saved using pdf files. This example has two plots in it. In [24]: pdf("myplots.pdf") # make whatever plot you want # first page plot(X, Y, pch=21, bg='red') # a new call to plot will make a new page plot(Y, X, pch=23, bg='blue') # close the current "device" which is this pdf file dev.off()  pdf: 2 ### Loops¶ It is easy to use for loops in R In [25]: for (i in 1:10) { print(i^2) }  [1] 1 [1] 4 [1] 9 [1] 16 [1] 25 [1] 36 [1] 49 [1] 64 [1] 81 [1] 100  In [26]: for (w in c('red', 'blue', 'green')) { print(w) }  [1] "red" [1] "blue" [1] "green"  Note that big loops can get really slow, a drawback of many high-level languages. ### Builtin help¶ R has a builtin help system, which can be accessed and searched as follows > help(t.test) > help.search('t.test')  Many objects also have examples that show you their usage > example(lm) In [27]: help(t.test)  In [28]: example(lm)  lm> require(graphics) lm> ## Annette Dobson (1990) "An Introduction to Generalized Linear Models". lm> ## Page 9: Plant Weight Data. lm> ctl <- c(4.17,5.58,5.18,6.11,4.50,4.61,5.17,4.53,5.33,5.14) lm> trt <- c(4.81,4.17,4.41,3.59,5.87,3.83,6.03,4.89,4.32,4.69) lm> group <- gl(2, 10, 20, labels = c("Ctl","Trt")) lm> weight <- c(ctl, trt) lm> lm.D9 <- lm(weight ~ group) lm> lm.D90 <- lm(weight ~ group - 1) # omitting intercept lm> ## No test: lm> ##D anova(lm.D9) lm> ##D summary(lm.D90) lm> ## End(No test) lm> opar <- par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0)) lm> plot(lm.D9, las = 1) # Residuals, Fitted, ... lm> par(opar) lm> ## Don't show: lm> ## model frame : lm> stopifnot(identical(lm(weight ~ group, method = "model.frame"), lm+ model.frame(lm.D9))) lm> ## End(Don't show) lm> ### less simple examples in "See Also" above lm> lm> lm>  ### Distributions in R¶ In practice, we will often be using the distribution (CDF), quantile (inverse CDF) of standard random variables like the T, F, chi-squared and normal. The standard 1.96 (about 2) standard deviation rule for$\alpha=0.05$: (note that 1-0.05/2=0.975) In [29]: qnorm(0.975)  1.95996398454005 We might want the$\alpha=0.05$upper quantile for an F with 2,40 degrees of freedom: In [30]: qf(0.95, 2, 40)  3.23172699283084 So, any observed F greater than 3.23 will get rejected at the$\alpha=0.05\$ level. Alternatively, we might have observed an F of 5 with 2, 40 degrees of freedom, and want the p-value

In [31]:
1 - pf(5, 2, 40)

0.0115292150460684

Let's compare this p-value with a chi-squared with 2 degrees of freedom, which is like an F with infinite degrees of freedom in the denominator (send 40 to infinity). We also should multiply the 5 by 2 because it's divided by 2 (numerator degrees of freedom) in the F.

In [32]:
c(1 - pchisq(5*2, 2), 1 - pf(5, 2, 4000))

1. 0.00673794699908548
2. 0.00678012054849153

Other common distributions used in applied statistics are norm and t.