Intro2R.Rmd

---
title: "Intro2R"
author: "Jonathan Rosenblatt"
date: "March 18, 2015"
output: html_document

---
# R Basics

Tips for this introduction:
- If you are working alone, consider starting with "An Introduction to R" here:
http://cran.r-project.org/manuals.html 
- Make sure you use RStudio.
- Ctrl+return to run lines from editor.
- alt+shift+k for RStudio keyboard shortcuts.
- Ctrl+alt+j to navigate between sections
- tab for auto-completion
- Ctrl+1 to skip to editor. 
- Ctrl+2 to skip to console.
- Ctrl+8 to skip to the environment list.
- Folding:
  - alt+l collapse chunk.
  - alt+shift+l unfold chunk.
  - alt+o collapse all.
  - alt+shift+o unfold all.
  


## Simple calculator
```{r example}
10+5
70*81
2**4
2^4
log(10)       					
log(16, 2)    					
log(1000, 10)   				
```


## Probability calculator 
Wish you knew this when you did Intro To Probability class?
```{r}
dbinom(x=3, size=10, prob=0.5) 	# For X~B(n=10, p=0.5) returns P(X=3)
dbinom(3, 10, 0.5)

pbinom(q=3, size=10, prob=0.5) # For X~B(n=10, p=0.5) returns P(X<=3) 	
dbinom(x=0, size=10, prob=0.5)+dbinom(x=1, size=10, prob=0.5)+dbinom(x=2, size=10, prob=0.5)+dbinom(x=3, size=10, prob=0.5) # Same as previous

qbinom(p=0.1718, size=10, prob=0.5) # For X~B(n=10, p=0.5) returns k such that P(X<=k)=0.1718

rbinom(n=1, size=10, prob=0.5) 	
rbinom(n=10, size=10, prob=0.5)
rbinom(n=100, size=10, prob=0.5)
```


## Getting help
Get help for a particular function.
```{r, eval=FALSE}
?dbinom 
help(dbinom)
```

Search local help files for a particular string.
```{r, eval=FALSE}
??binomial
help.search('dbinom') 
```

Load a menu with several important manuals:
```{r, eval=FALSE}
help.start() 
```


## Variable asignment:
Assignments into a variable named "x":
```{r}
x = rbinom(n=1000, size=10, prob=0.5) # Works. Bad style.
x <- rbinom(n=1000, size=10, prob=0.5) # Asignments into a variable named "x"
```
More on style: http://adv-r.had.co.nz/Style.html


Print contents:
```{r}
x
print(x)  
(x <- rbinom(n=1000, size=10, prob=0.5))  # Assign and print.
```


Operate on the object
```{r}
mean(x)  
var(x)  
hist(x)  
rm(x) # remove variable
```


For more information on distributions see http://cran.r-project.org/web/views/Distributions.html


## Piping for better style and readability
```{r}
# install.packages('magrittr')
library(magrittr)
```

```{r}
x <- rbinom(n=1000, size=10, prob=0.5)

x %>% var() # Instead of var(x)
x %>% hist()  # Instead of hist(x)
x %>% mean() %>% round(2) %>% add(10) 
```

This example clearly demonstrates the benefits (from http://cran.r-project.org/web/packages/magrittr/vignettes/magrittr.html)
```{r}
# Functional (onion) style
car_data <- 
  transform(aggregate(. ~ cyl, 
                      data = subset(mtcars, hp > 100), 
                      FUN = function(x) round(mean(x, 2))), 
            kpl = mpg*0.4251)


# magrittr style
car_data <- 
  mtcars %>%
  subset(hp > 100) %>%
  aggregate(. ~ cyl, data = ., FUN = . %>% mean %>% round(2)) %>%
  transform(kpl = mpg %>% multiply_by(0.4251)) %>%
  print
```


## Vector creation and manipulation 
```{r}
c(10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21)
10:21 							
seq(from=10, to=21, by=1) 							
x <- seq(from=10, to=21, by=2) 								
x <- c(10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21) 	
x
```



You can assign AFTER the computation is finished:
```{r}
c(1,2,3)
y<- .Last.value 
y
```


Operations usually work element-wise:
```{r}
x+2
x*2    
x^2    
sqrt(x)  
log(x)   
```


## Simple plotting 
```{r}
x<- 1:100; y<- 3+sin(x) # Create arbitrary data
plot(x = x, y = y) # x,y syntax  						
plot(y ~ x) # y~x syntax (I like better)
```

Control plot appearance:   
```{r}
plot(y~x, type='l', main='Plotting a connected line')
plot(y~x, type='h', main='Sticks plot', xlab='Insert x axis label', ylab='Insert y axis label')
plot(y~x, pch=5)
plot(y~x, pch=10, type='p', col='blue', cex=4)
abline(3, 0.002)
```

Available plotting options
```{r, eval=FALSE}
example(plot)
example(points)
?plot
help(package='graphics')
```

When your plotting gets serious, move to `ggplot2` and `ggvis` as soon as possible.




___


## Data frame Manipulation 
`data.frames` extend the `matrix` class, in that they allow the binding of vectors of several classes (with same length).
```{r}
x<- 1:100; y<- 3 + sin(x) 
class(x) # R (high) level representation of an object.

# mode(x) 
# typeof(x) 
```


Create and checkout your first data frame
```{r}
frame1 <- data.frame(x=x, sin=y)	
frame1
head(frame1)
frame1 %>% head() # just print the beginning
frame1 %>% View() # Excel-like view (never edit!)

class(frame1) # the object is of type data.frame
dim(frame1)  							
dim(x)
length(frame1)
length(x)

str(frame1) # the inner structure of an object
attributes(frame1) # get the object's meta data
```

### Exctraction 
single element:
```{r}
frame1[1, 2]    						
frame1[2, 1]     						
```

Extract _column_ by index:
```{r}
frame1[1, ]      						
frame1[,1] %>% t
frame1[,1] %>% t %>% dim
```

Extract column by name:
```{r}
names(frame1)   						
frame1[, 'sin']
dim(frame1[, 'sin'])  # extract as a vector. no dim attribute.
frame1['sin'] 
dim(frame1['x',]) # extract as a data.frame. has dim attribute.
frame1[,1:2] %>% class
frame1[2] %>% class
frame1[2, ] # extract a row

frame1$sin %>% class
```

`subset()` does the same
```{r}
subset(frame1, select=sin) 
subset(frame1, select=2)
subset(frame1, select= c(2,0))
```


Sanity conservation notice!
Always think if you want to extract a vector or a frame:
- Note the difference between `[]` and `[[]]` extraction!
- Note the difference between `frame[,1]` and `frame[1]`.
```{r}
a <- frame1[1]
b <- frame1[[1]]
a==b # Seems identical. But not really:
class(a)
class(b)
# Causes different behaviour:
a[1]
b[1]
```

More about extraction: http://adv-r.had.co.nz/Subsetting.html

### dplyr package 
`dplyr` makes the manipulation of data.frames a breeze.
It is very fast, and straightforward to use.

Install the package:
```{r}
# install.packages('dplyr')
```

The following examples are taken from:
https://github.com/justmarkham/dplyr-tutorial/blob/master/dplyr-tutorial.Rmd
```{r}
# install.packages('nycflights13')
library(nycflights13)
dim(flights)
View(flights)
names(flights)
class(flights) # a tbl_df is an extension of the data.frame class
library(dplyr) # calling dplyr

filter(flights, month == 1, day == 1) #dplyr style
flights[flights$month == 1 & flights$day == 1, ] # old style
flights %>% filter(month == 1, day == 1) # dplyr with magrittr style (yes!)

filter(flights, month == 1 | month == 2)
slice(flights, 1:10) # selects rows

arrange(flights, year, month, day) # sort
arrange(flights, desc(arr_delay)) # sort descending

select(flights, year, month, day) # select columns
select(flights, year:day) # select column range
select(flights, -(year:day)) # drop columns
rename(flights, tail_num = tailnum) # rename variables
# add a new computed colume
mutate(flights,
  gain = arr_delay - dep_delay,
  speed = distance / air_time * 60) 
# you can refer to columns just created!
mutate(flights,
  gain = arr_delay - dep_delay,
  gain_per_hour = gain / (air_time / 60)
)
# keep only new variables
transmute(flights,
  gain = arr_delay - dep_delay,
  gain_per_hour = gain / (air_time / 60)
)
# simple statistics
summarise(flights,
  delay = mean(dep_delay, na.rm = TRUE)
  )

sample_n(flights, 10) # random subsample
sample_frac(flights, 0.01) # random subsample
```

Subgroup operations
```{r}
by_tailnum <- group_by(flights, tailnum)
by_tailnum %>% class # a groupping object
delay <- summarise(by_tailnum,
  count = n(),
  avg.dist = mean(distance, na.rm = TRUE),
  avg.delay = mean(arr_delay, na.rm = TRUE))
delay <- filter(delay, count > 20, avg.dist < 2000)
View(delay)

destinations <- group_by(flights, dest)
summarise(destinations,
  planes = n_distinct(tailnum),
  flights = n()
)

# Grouping works in a hirarchy. summarise() peels outer layer.
daily <- group_by(flights, year, month, day)
(per_day   <- summarise(daily, flights = n()))
(per_month <- summarise(per_day, flights = sum(flights)))
(per_year  <- summarise(per_month, flights = sum(flights)))
```





Two table operations
```{r}
airlines %>% View
flights2 <- flights %>% select(year:day, hour, origin, dest, tailnum, carrier)

flights2 %>% left_join(airlines) # join on left table with automatic matching.

flights2 %>% left_join(weather)

flights2 %>% left_join(planes, by = "tailnum") # with named matching

flights2 %>% left_join(airports, c("dest" = "faa"))

flights2 %>% left_join(airports, c("origin" = "faa"))
```

Types of join
```{r}
(df1 <- data_frame(x = c(1, 2), y = 2:1))
(df2 <- data_frame(x = c(1, 3), a = 10, b = "a"))

df1 %>% inner_join(df2) # SELECT * FROM x JOIN y ON x.a = y.a

df1 %>% left_join(df2) # SELECT * FROM x LEFT JOIN y ON x.a = y.a

df1 %>% right_join(df2) # SELECT * FROM x RIGHT JOIN y ON x.a = y.a
df2 %>% left_join(df1) 

df1 %>% full_join(df2) # SELECT * FROM x FULL JOIN y ON x.a = y.a

# return only unmatched cases
flights %>%
  anti_join(planes, by = "tailnum") %>% 
  count(tailnum, sort = TRUE) 
# SELECT * FROM x WHERE NOT EXISTS (SELECT 1 FROM y WHERE x.a = y.a)

df1 %>% semi_join(df2, by = "x")  # SELECT * FROM x WHERE EXISTS (SELECT 1 FROM y WHERE x.a = y.a)
```

Set operations
```{r}
(df1 <- data_frame(x = 1:2, y = c(1L, 1L)))
(df2 <- data_frame(x = 1:2, y = 1:2))

intersect(df1, df2) # SELECT * FROM x INTERSECT SELECT * FROM y

union(df1, df2) # SELECT * FROM x UNION SELECT * FROM y

setdiff(df1, df2) # SELECT * FROM x EXCEPT SELECT * FROM y

setdiff(df2, df1)
```

Leaving dplyr for now...


## Data Import and export 

__Note__: The [readr](https://github.com/hadley/readr) package facilitates and accelerates data importing. This section should be updated to use it.

For a complete review see:
http://cran.r-project.org/doc/manuals/R-data.html
also in  help.start() -> "Import and Export Manual" 


### Import from WEB 
`read.table()` is the main importing workhorse.
```{r}
URL <- 'http://statweb.stanford.edu/~tibs/ElemStatLearn/datasets/bone.data'
tirgul1 <- read.table(URL)
```

Always look at the imported result!
```{r}
View(tirgul1)
# hmmm... header interpreted as data. Fix with header=TRUE:
tirgul1 <- read.table(URL, header = TRUE) 
View(tirgul1)
```

### Import .csv files
Let's write a simple file so that we have something to import:
```{r}
View(airquality) #  examine the data to export
(temp.file.name <- tempfile()) # get an arbitrary file name
write.csv(x = airquality, file = temp.file.name) #export
```

Now let's import:
```{r}
# my.data<- read.csv(file='/home/jonathan/Projects/...')
my.data<- read.csv(file=temp.file.name)
View(my.data)
```

__Note__: Under MS Windows(R) you might want need '\\\' instead of '/'

### Imprt .txt files 
Tries to guess the separator
```{r, eval=FALSE}
my.data<- read.table(file='C:\\Documents and Settings\\Jonathan\\My Documents\\...') #
```
Specifies the separator explicitly
```{r, eval=FALSE}
my.data<- read.delim(file='C:\\Documents and Settings\\Jonathan\\My Documents\\...') 
```
If you care about your sanity, see ?read.table before starting imports.

### Writing Data to files

Get and set the current directory:
```{r, eval=FALSE}
getwd() #What is the working directory?
setwd() #Setting the working directory in Linux
```

```{r}
write.csv(x=tirgul1, file='/tmp/tirgul1.csv') #
```

See ?write.table for details.

### .XLS(X) files 
Strongly recommended to convert to .csv
If you still insist see:
http://cran.r-project.org/doc/manuals/R-data.html#Reading-Excel-spreadsheets

### Massive files 
Better store as matrices and not data.frames.
`scan()` is faster than `read.table()` but less convenient:

Create the example data:
```{r}
cols<- 1e3
# Note: On Windoes you might neet to change /tmp/A.txt to /temp/A.txt 
rnorm(cols^2) %>%
  matrix(ncol=cols) %>%
  write.table(file='/tmp/A.txt', col.names= F, row.names= F)
# Measure speed of import:
system.time(A<- read.table('/tmp/A.txt', header=F))
system.time(A <- scan(file='/tmp/A.txt', n = cols^2) %>%
              matrix(ncol=cols, byrow = TRUE))

file.remove('/tmp/A.txt') 
```

This matter will be revisited in the last class.

### Databases 
Start [here](https://rforanalytics.wordpress.com/useful-links-for-r/odbc-databases-for-r/)

### Hands on example (from the WEB)
```{r}
URL <- 'http://statweb.stanford.edu/~tibs/ElemStatLearn/datasets/bone.data'
tirgul1 <- read.table(URL, header = TRUE)

names(tirgul1)
tirgul1 %>% head
tirgul1 %>% tail
View(tirgul1)
dim(tirgul1)
length(tirgul1)
```

R can be object oriented (read about S3 and S4 if interested).
See how `summary()` behaves differently on different object classes:
```{r}
class(tirgul1[, 1]); class(tirgul1[, 2]); class(tirgul1[, 3]); class(tirgul1[, 4])
summary(tirgul1)
```



Matrix is more efficient than data frames. But can store only a single class of vectors.
```{r}
tirgul.matrix <- as.matrix(tirgul1) 
tirgul.matrix
class(tirgul.matrix)
# notice everything has been cast to the most general class.
class(tirgul.matrix[, 1]); class(tirgul.matrix[, 2]); class(tirgul.matrix[, 3]); class(tirgul.matrix[, 4])
summary(tirgul.matrix)
```

Note: if copy-pasting an expression bothers you (as it should!), here are some solutions:
```{r}
# The apply family of functions:
sapply(tirgul.matrix, class)

# looping
for(j in 1:ncol(tirgul.matrix)) print(class(tirgul.matrix[,j]))
```

Make sure you read `?sapply`. 
LISP fans might also like to read `?MAP`.



Operations _within_ data objects:
```{r}
plot(tirgul1$gender)
tirgul1$gender %>% plot() # 
with(tirgul1, plot(gender) ) # Same opration. Different syntax.

mean(tirgul1$age)
tirgul1$age %>% mean() # 
with(tirgul1, mean(age) ) # Same opration. Different syntax.
```


```{r}
tirgul1$age <- tirgul1$age * 365
tirgul1<- transform(tirgul1, age=age*365 )  #Age in days
with(tirgul1, mean(age) )
tirgul1<- transform(tirgul1, age=age/365 )  #Does this revert back to years?
with(tirgul1, mean(age) )
```

Then again, many of these functions are replaced by more friendly functions in the `dplyr` package (see below).


## Sorting 
```{r}
(x<- c(20, 11, 13, 23, 7, 4))
(y<- sort(x))
(ord<- order(x))
x[ord] # Exctracting along the order is the same as sorting.
ranks<- rank(x)
identical(y[ranks] , x) # Compares two objects

(z<- c('b','a','c','d','e','z'))
xz<- data.frame(x,z)
sort(xz)
xz[ord,] # Sorting a data frame using one column
```


## Looping 
For a crash course in R programming (not only data analysis) try:   
http://adv-r.had.co.nz/                            
The usual for(), while(), repeat() 
```{r}
for (i in 1:100){
    print(i)
    }
```


```{r}
for (helloeveryone in seq(10, 100, by=2) ){
    print(helloeveryone)
    }
```


## Recursion 
Typically very slow due to memory management issues.

```{r}
fib<-function(n) {
    if (n < 2) fn<-1 
    else fn<-Recall(n - 1) + Recall(n - 2) 
    return(fn)
} 
fib(30)
```


## Finding your objects 
```{r}
ls() #Lists all available objects
ls(pattern='x')

ls(pattern='[0-9]') # Search using regular expressions
ls(pattern='[A-Z]')
```

Ctrl+8 in RStudio.




# Univariate Exploratory Statistics


##  Exploring Categorical Variables  
```{r}
gender <- c(rep('Boy', 10), rep('Girl', 12))
drink <- c(rep('Coke', 5), rep('Sprite', 3), rep('Coffee', 6), rep('Tea', 7), rep('Water', 1))  
class(gender);class(drink)

cbind(gender, drink)
table1 <- table(gender, drink) 
table1										
```


              

##  Exploring Continous Variables 

Generating and exploring data
```{r}
sample1 <- rnorm(100) 							
table(sample1) 									
hist(sample1, freq=T, main='Counts')      	
hist(sample1, freq=F, main='Frequencies') 	
lines(density(sample1))                  	
rug(sample1)
```


## The Boxplot 
```{r}
boxplot(sample1)	
```



Several different visualizations:
```{r}
sample2<-rnorm(1000)     
stem(sample2)          
hist(sample2)          
plot(density(sample2))  
rug(sample2)
```



True data 
```{r}
URL <- 'http://statweb.stanford.edu/~tibs/ElemStatLearn/datasets/bone.data'
bone <- read.table(URL, header = TRUE)
names(bone)
summary(bone) 			
stripchart(bone['age'])
stem(bone[, 'age']) 									
hist(bone[, 'age'], prob=T) 							
lines(density(bone[, 'age'])) 
with(bone, rug(age))

ind<-bone[, 'gender']=='male'

boxplot(bone$age~bone$gender)
```


## Graphical parameters 
```{r}
attach(bone) 
stripchart(age)
stripchart(age~gender)
stripchart(age~gender, v=T)

boxplot(age~gender)
boxplot(age~gender, horizontal=T, col=c('pink','lightblue') )
title(main='Amazing Boxplots!')
title(sub="Well actually.. I've seen better Boxplots")

plot(density(age), main='')
plot(density(age), main='', type='h')
plot(density(age), main='', type='o')
plot(density(age), main='', type='p')
plot(density(age), main='', type='l')

?plot.default

plot(density(age),main='')
rug(age)
boxplot(age, add=T, horizontal=T, at=0.02, boxwex=0.05, col='grey')
title(expression(alpha==f[i] (beta)))
example(plotmath)

par(mfrow=c(2,1))
(males<- gender=='male')
plot(density(age[males]), main='Male') ; rug(age[males])
plot(density(age[!males]), main='Female') ; rug(age[!males])

range(age)
plot(density(age[males]), main='Male', xlim=c(9,26)) ; rug(age[males])
plot(density(age[!males]), main='Female', xlim=c(9,26)) ; rug(age[!males])
par(mfrow=c(1,2))
plot(density(age[males]), main='Male', xlim=c(9,26)) ; rug(age[males])
plot(density(age[!males]), main='Female', xlim=c(9,26)) ; rug(age[!males])

par(mfrow=c(1,1),ask=T)
plot(density(age[males]), main='Male', xlim=c(9,26)) ; rug(age[males])
plot(density(age[!males]), main='Female', xlim=c(9,26)) ; rug(age[!males])
```


## Integer data 
Integer data will most certainly produce overlaps if plotted. Either add hitter, or treat as discrete.
```{r}
r.age<-round(age)
plot(density(r.age))
rug(r.age)
plot(density(r.age, from=9))
rug(jitter(r.age))
hist(r.age)
rug(jitter(r.age))
```


## Plotting

### Preparing data for plotting
2D data can be in either _wide_ or _long_ format. 
Most R functions are designed for long formats. 
Let's start by trying to plot in the wide format.
Notice each dosage is plotted separately (yes, I could have looped).
```{r}
wide.data<-data.frame(id=1:4, age=c(40,50,60,50), dose1=c(1,2,1,2),dose2=c(2,1,2,1), dose4=c(3,3,3,3))
wide.data

plot(dose1~age, data=wide.data, ylim=range(c(dose1,dose2,dose4)), ylab='')
points(dose2~age, data=wide.data, pch=2)
points(dose4~age, data=wide.data, pch=3)
```


Plotting in long format is much easier. 
I will first convert the data manually.
```{r}
(dose.type<-c(
		rep('dose1', length(wide.data$dose1)),
		rep('dose2', length(wide.data$dose2)),
		rep('dose4', length(wide.data$dose4))))
(dose<- c(wide.data$dose1,wide.data$dose2,wide.data$dose4))
(long.id<- rep(wide.data$id,3))
(long.age<- rep(wide.data$age,3))

long.data <- data.frame(long.id, long.age, dose.type, dose)
View(long.data)

plot(dose~long.age, data=long.data, pch=as.numeric(dose.type))
```
I will now try to avoid this manual reshaping.

### Reshaping data with `tidyr` package

This is the package I recommend if you cannot reshape manually.
Example from [here](http://blog.rstudio.org/2014/07/22/introducing-tidyr/)
```{r}
library(tidyr)
library(dplyr)

# Data in wide format:
messy <- data.frame(
  name = c("Wilbur", "Petunia", "Gregory"),
  a = c(67, 80, 64),
  b = c(56, 90, 50)
)
messy

# Convert to long format:
messy %>% gather(drug, heartrate, a:b)
```

```{r}
# Another example- from wide to long:
set.seed(10)
messy <- data.frame(
  id = 1:4,
  trt = sample(rep(c('control', 'treatment'), each = 2)),
  work.T1 = runif(4),
  home.T1 = runif(4),
  work.T2 = runif(4),
  home.T2 = runif(4)
)
messy %>% head
tidier <- messy %>%  gather(key, time, -id, -trt)
tidier %>% head(8)

# From long to wide
tidy <- tidier %>%
  separate(key, into = c("location", "time"), sep = "\\.") 
tidy %>% head(8)
```

### Fancy Plotting 
```{r}
library(ggplot2)
URL <- 'http://statweb.stanford.edu/~tibs/ElemStatLearn/datasets/bone.data'
bone <- read.table(URL, header = TRUE)
qplot(spnbmd, data=bone)
qplot(x=gender, y=spnbmd, data=bone, geom='boxplot')
qplot(spnbmd, data=bone, geom='histogram')+ facet_wrap(~gender)
qplot(spnbmd, data=bone, geom='density')+ facet_wrap(~gender)
qplot(spnbmd, data=bone)+ geom_density(col='red', size=1)+ facet_wrap(~gender)
qplot(spnbmd, data=bone, fill=gender, geom='density', alpha=1)
```

Diamonds example (Taken from Wickham's web site: http://had.co.nz/stat405/)
```{r}
?diamonds
dim(diamonds)
head(diamonds)
```

```{r}
qplot(carat, data = diamonds)
qplot(carat, data = diamonds, binwidth = 1)
qplot(carat, data = diamonds, binwidth = 0.1)
qplot(carat, data = diamonds, binwidth = 0.01)
resolution(diamonds$carat)
last_plot() + xlim(0, 3)

qplot(depth, data = diamonds, binwidth = 0.2)
qplot(depth, data = diamonds, binwidth = 0.2,fill = cut) + xlim(55, 70)
qplot(depth, data = diamonds, binwidth = 0.562) +xlim(55, 70) + facet_wrap(~ cut)

qplot(table, price, data = diamonds)
qplot(table, price, data = diamonds, geom = "boxplot")
qplot(table, price, data = diamonds, geom="boxplot",group = round(table))

qplot(carat, price, data = diamonds)
qplot(carat, price, data = diamonds, alpha = I(1/10))

qplot(carat, price, data = diamonds, geom = "bin2d", main='Count Heatmap')
qplot(carat, price, data = diamonds, geom = "hex")
qplot(carat, price, data = diamonds) + geom_smooth()
```


For more information on ggplot2 see http://had.co.nz/ggplot2 


## The QQ plot 
A simple and efficient tool to compare between distributions.
```{r}
mystery.2<-function(y) {
  n<-length(y)
  y<-sort(y)
  i<-1:n
  q<-(i-0.5)/n
  x<-qnorm(q, mean(y), sqrt(var(y)))
  plot(y~x, xlab='Theoretical Quantiles', ylab='Empirical Quantiles')
}

normals.1<-rnorm(100, 0, 1); hist(normals.1)
mystery.2(normals.1); abline(0, 1)

normals.2<-rnorm(100, 0, 10); hist(normals.2)
mystery.2(normals.2); abline(0, 1)

## No need to write the function every time...
qqnorm(normals.1)   
qqnorm(normals.2)   

## How would non-normal observations look? ##
non.normals.1<-runif(100); hist(non.normals.1)
mystery.2(non.normals.1); abline(0, 1)

non.normals.2<-rexp(100, 1); hist(non.normals.2)
mystery.2(non.normals.2); abline(0, 1)

non.normals.3<-rgeom(100, 0.5); hist(non.normals.3)
mystery.2(non.normals.3); abline(0, 1)

## Adapting for a non-normal distribution: ##
qq.uniform<-function(y) {
  n<-length(y);    y<-sort(y);    i<-1:n;    q<-(i-0.5)/n
  x<-qunif(q, min=min(y), max=max(y)) #each disribution will require it's own parameters!
  plot(y~x, xlab='Theoretical Quantiles', ylab='Empirical Quantiles')
}
qq.uniform(non.normals.1);abline(0, 1)
qq.uniform(non.normals.2);abline(0, 1)
qq.uniform(normals.2);abline(0, 1)
```





# Multiple data vectors 
We now leave the single-vector world and move to the analysis of dependencies between several vectors. 

## Scatter plots
```{r}
# Sine function
x<-seq(-pi, pi, 0.01)
y<-sin(x)
plot(y~x)

#Exponent function
x<-seq(-pi, pi, 0.01)
y<-exp(x)
plot(y~x)

# Sinc function
x<-seq(-10*pi, 10*pi, 0.01)
y<-sin(x)/x
plot(y~x)

# Fancy function
x<-seq(-pi, pi, 0.01)
y<-sin(exp(x))+cos(2*x)
plot(y~x)
plot(y~x, type='l')
plot(y~x, type='o')

## Some real life data
URL <- 'http://statweb.stanford.edu/~tibs/ElemStatLearn/datasets/ozone.data'
ozone <- read.table(URL, header=T)
names(ozone)
plot(ozone)
```


## 3D plotting  
```{r}
# install.packages('rgl')
library(rgl)
plot3d(ozone[, 1:3]) 
```


## Plotting a surface 
```{r}
x <- seq(0, 1, 0.01)
y <- seq(0, 1, 0.01)
xy.grid <- expand.grid(x, y)
func1 <- function(mesh) exp(mesh[, 1]+mesh[, 2])
z <- func1(xy.grid)
xyz <- data.frame(xy.grid, z)
plot3d(xyz, xlab='x', ylab='y')  
```


##  Fitting linear lines and surfaces  
We will now try and fit linear surfaces to our data.

### Well behaved data 
```{r}
x <- 1:100
a <- 2
b <- 3.5
sigma <- 10
y <- a+b*x+rnorm(100, 0, sigma)
plot(y~x)
```

### Ordinary Least Squares 
```{r}
ols.line<-function(x, y){
    sxy<-sum( (x-mean(x) ) * (y-mean(y) ) )    
    sxx<-sum( (x-mean(x)) ^ 2 )
    b1<-sxy / sxx
    a1<-mean(y) - b1 * mean(x)
    return(list(slope=b1, intercept=a1))
}

ols<-ols.line(x, y) ; ols
abline(ols$intercept, ols$slope, lty=2, lwd=3)
predictions <-  ols$intercept + ols$slope * x
residuals<- y - predictions
plot(residuals) ; abline(h=0)
```

### Dangers of Extrapolation  
```{r}
x<-runif(1000)*5
y<-exp(x)+rnorm(1000)
plot(y~x, main='Whole relation')

rect(xleft=0, ybottom=-5, xright=2, ytop=10)

plot(y~x, main='Local relation', cex=0.5, xlim=c(0, 2), ylim=c(-5, 10));abline(v=2, lty=3)

ind<-x<=2;ind
ols.interpolating<-ols.line(x[ind], y[ind]);ols.interpolating
abline(ols.interpolating$intercept ,  ols.interpolating$slope, col='red')
text(x=0.5, y=6, labels='Interpolates Nicely', cex=2)

plot(y~x, main='Whole relation')
abline(ols.interpolating$intercept ,  ols.interpolating$slope, col='red')
abline(v=2, lty=3)
text(x=2, y=121, labels='Extrapolates Terribly!', cex=2)

# Non-linearity might be fixed with a transformation:
# Which of the following looks better (more linear)? 
plot(y~exp(x))
plot(log(y)~x)
plot(log(y)~log(x))
```

### Multivariate linear regression 
```{r}
# install.packages('rgl')
library(rgl)

xy.grid <- data.frame(x1=runif(10000), x2=runif(10000))

func1<-function(mesh, a0, a1, a2, sigma) {
    n<-nrow(mesh)
    a0 + a1 * mesh[, 1] + a2 * mesh[, 2] + rnorm(n, 0, sigma)
    }
    
# More noise hides the stucture in the data:
z<-func1(xy.grid, a0=5, a1=1, a2=3, .0); z; xyz=data.frame(xy.grid, z); plot3d(xyz, xlab='x1', ylab='x2')
z<-func1(xy.grid, a0=5, a1=1, a2=3, .4); xyz=data.frame(xy.grid, z); plot3d(xyz, xlab='x1', ylab='x2')
z<-func1(xy.grid, a0=5, a1=1, a2=3, 11); xyz=data.frame(xy.grid, z); plot3d(xyz, xlab='x1', ylab='x2')

```

`lm()` is the main workhorse for OLS solving $(X'X)^{-1} X'y$ with the QR decomposition.
```{r}
z<-func1(xy.grid, a0=5, a1=1, a2=3, .4)
xyz=data.frame(xy.grid, z)
plot3d(xyz, xlab='x1', ylab='x2')
lm(z~., xyz) # Did we exctract the correct coefficients?
```



# Date handeling 
See the `lubridate` package and manual [here](http://cran.r-project.org/web/packages/lubridate/vignettes/lubridate.html).


# String handelind 
```{r}
print("Hello\n") 	# Wrong!
show("Hello\n") 	# Wrong!
cat("Hello\n")		# Right!

# Windows directories need double escapes:
print("C:\\Program Files\\") 
cat("C:\\Program Files\\", sep="\n")

# String concatenation:
paste("Hello", "World", "!")
paste("Hello", "World", "!", sep="")
paste("Hello", " World", "!", sep="")

x <- 5
paste("x=", x)
paste("x=", x, paste="")

cat("x=", x, "\n") #Too many spaces :-(
cat("x=", x, "\n", sep="")

# Collapsing strings:
s <- c("Hello", " ", "World", "!")
paste(s)
paste(s, sep="")
paste(s, collapse="")
paste(s, collapse=" 1")


s <- c("Hello", "World!")
paste(1:3, "Hello World!")
paste(1:3, "Hello World!", sep=":")
paste(1:3, "Hello World!", sep=":", collapse="\n")
cat(paste(1:3, "Hello World!", sep=":", collapse="\n"), "\n") # cat() does not collapse :-(


# Substrings:
s <- "Hello World"
substring(s, start=4, stop=6)

# Splits:
s <- "foo, bar, baz"
strsplit(s, ", ")

s <- "foo-->bar-->baz"
strsplit(s, "-->")

# Using regular expressions (see ?regexp):
s <- "foo, bar, baz"
strsplit(s, ", *")
strsplit(s, "")

# Looking in *vectors* of strings:
(s <- apply(matrix(LETTERS[1:24], nr=4), 2, paste, collapse=""))

grep("O", s) # Returns location
grep("O", s, value=T) # Returns value


regexpr(pattern="o", text="Hello")
regexpr(pattern="o", text=c("Hello", "World!"))

s <- c("Hello", "World!")
regexpr("o", s)
s <- c("Helll ooo", "Wrld!")
regexpr("o", s)

# Fuzzy (approximate) matches:
grep ("abc", c("abbc", "jdfja", "cba")) 	# No match :-(
agrep ("abc", c("abbc", "jdfja", "cba")) 	# Match! :-)

## Note: agrep() is the function used in help.search()
s <- "foo bar baz"
gsub(pattern=" ", replacement="", s)   # Remove all the spaces
s <- "foo  bar   baz"
gsub("  ", " ", s)
gsub(" +", "", s) # Using regular expression
gsub(" +", " ", s)  # Remove multiple spaces and replace them by single spaces

s <- "foo bar baz"
sub(pattern=" ", replacement="", s) # sub() only replaces first occurance.
gsub("  ", " ", s)
```


If you use strings often, try the `stringr` package.