r.rmd

---
title: "R cheatsheet"
author: Emil Rehnberg
bibliography: refs.bib
date: "`r Sys.Date()`"
csl: shiki.csl
output:
  pdf_document:
    highlight: zenburn
  html_document:
    toc: TRUE
    toc_float: TRUE
    css: styles.css
---

```{r set-options, echo=FALSE, cache=FALSE}
options(width = 200)
```

```{r echo=FALSE, eval=FALSE}

library(rmarkdown); require(shiny)
rmdFilePath <- "r.rmd"
allexamples <- FALSE
rmarkdown::render(rmdFilePath, output_format="html_document") # "all_document"

# TODO
# - try out Bundler for R -- https://github.com/robertzk/lockbox
```

```{r, eval=FALSE, echo = FALSE}

devtools::install_github('cran/colorout')

```

```{r echo=FALSE, message=FALSE}
set.seed(308)
require(colorout)
require(dplyr)
require(tidyr)
require(magrittr)
require(ggplot2)
```

```{r "initiate polymorphic exclude", echo = FALSE}
include <- function(obj, elms) UseMethod("include")
exclude <- function(obj, elms) UseMethod("exclude")
```

# 目的

cheatsheet for R. it's a place to dump typical annotated R-code.

# Datasets {.tabset .tabset-fade .tabset-pills}

## Overview

`R` has several open datasets

## numerical

```{r}
library(magrittr)
airquality %>% str
mtcars %>% str
USArrests %>% str
```

## Time-series

```{r}
AirPassengers %>% str
Nile %>% str
freeny %>% str # with numerical variables also!
```

## factors

```{r}
CO2 %>% str
ChickWeight %>% str # Ord.factor ordinal factor!
esoph %>% str # multiple Ord.factor
iris %>% str
npk %>% str
warpbreaks %>% str
# library(randomForest)
# data(imports85)
# imports85 %>% str # factors & ord.factor

data(quine, package = "MASS")
quine %>% str
```

## lists

```{r}
Harman23.cor %>% str
```

## tables

multidimensional tables!

```{r}
Titanic %>% str
```

## Survival data suitable data

```{r}
data("GBSG2", package = "TH.data")
GBSG2 %>% str
with(GBSG2, { survival::Surv(time, cens) }) %>% head(10)
```

# R-language design {.tabset .tabset-fade .tabset-pills}

## Overview

Section for R-language design and concepts. Packages, code, data loading

## loading / unloading packages

in order to attach / load / require packages use `library` (or `require` if you like problems).

to detach / unload / remove a package from a session you can use `detach`

```{r}
exists("is_in") # from magrittr package
detach(package:magrittr)
exists("is_in")
library(magrittr)
exists("is_in")
```

## error handling

use `tryCatch`

with error handler

```{r "error / exception handling", echo = TRUE}

library(magrittr)
myError <- simpleError("エッラー！")

tryCatch({
  print(1);
  stop(myError)
  print(2)
}, error = function(err){
  message("error function messages:")
  (err %>% print)
}, finally = {
  print("always finish run with finally()")
  message("finished tryCatch")
})

```

with warning handler

```{r}

tryCatch({
  print(1)
  warning("ACHTUNG! ACHTUNG!")
  print(2)
}, warning = function(war){
  message("warning function messages:")
  (war %>% print)
}, finally = {
  print("always finish run with finally()")
  message("finished tryCatch")
})

```

`stop` returns exit code 1 (check with `$?` in ruby and `$?.exitstatus` especially)

you can also use the `tools` package. but `pskill` does not returns a 0 exit status. `pskill` seems to be rather rough around the edges, use with care.

```{r}
# tools::pskill(Sys.getpid(), tools::SIGINT)
```

## Control Flow

E.g. if / else statements

```{r "control flow, all the below can be parsed"}

if(1 == 2)
  print("1 == 2") else
if(2 == 2)
  print("2 == 2") else
  print("ALL FALSE")

if(1 == 2) {
  print("1 == 2")
} else if(2 == 2) {
    print("2 == 2")
  } else {
    print("ALL FALSE")
  }

if(1 == 2) {
  print("1 == 2")
} else {
  if(2 == 2) {
    print("2 == 2")
  } else {
    print("ALL FALSE")
  }
}

```

## ellipsis

`...` / ellipsis or "dot-dot-dot"

use an arbitrary number of arguments to a function

```{r}
a <- c(1,  2,  NA, 4, NA)
b <- c(NA, NA, NA, 5, 6)
c <- c(7,  8,  NA, 9, 10)
```

catch the arguments in a list

```{r}
elTrail1 <- function(...) list(...)
```
catch the arguments in a data frame

```{r}
elTrail2 <- function(...) data.frame(...)
```

extract the first argument

```{r}
elTrail3 <- function(...) ..1
```

extract the second argument

```{r}
elTrail4 <- function(...) ..2
```

## Environments

```{r}
exists("abc")
abc <- 2
exists("abc")
exists("pi", envir = emptyenv())
exists("glm", envir = as.environment("package:base"))
exists("glm", envir = as.environment("package:stats"))
```

## dump object to executable R code

`dput` is your friend. NOTE: not suitable for S4 objects. try `writeRDS` or `save` for that.

there's some gotchas though. When you're dumping a factor - the output might end up crazy long. Use `as.character` to get what you're looking for.

```{r "dput gotcha"}
library(dplyr)
(USArrests %<>% mutate(state = USArrests %>% row.names %>% as.factor)) %>% str
USArrests$state %>% head(2) %>% dput
USArrests$state %>% head(2) %>% as.character %>% dput
```

## parse CL args

parsing command line arguments / parameters

- use `#!/usr/bin/env Rscript` on top of the executable file (`tca.r` let's say)
- call `./tca.r --args a=1 name=\'Billy\'`
- use the `commandArgs`, `parse` and `eval` functions to read the parameters into your script
- sometimes if `trailingOnly` is set, `--args` is still in there

```{r "commandArgs"}
library(magrittr)
arguments <- c("a=1", "name=\'Billy\'") # commandArgs(trailingOnly = TRUE)
for(argument in arguments){
  eval(parse(text = argument))
}
```

into a `data.frame`. This works but then all columns are strings

```{r "command line arguments into data.frame"}
arguments %>% strsplit("=") %>% sapply(function(strings) purrr::set_names(strings[2], strings[1]) )
```

use `jsonlite` instead by sending in JSON data. E.g. executing:
`./predict-co2uptake '[{"Type": "Quebec", "Treatment": "nonchilled", "conc": 350}]'`
on the command line

```{r "command line arguments into data.frame using jsonlite"}
jsonlite::fromJSON('[{"Type": "Quebec", "Treatment": "nonchilled", "conc": 350}]') %>% str
```


## Object system {.tabset .tabset-fade .tabset-pills}

### Overview

R has a S3 and S4 Object system. S3 is preferred using polymorphism.
Ref: [adv-R](http://adv-r.had.co.nz/OO-essentials.html)

### Base types

R "base types" are like a C struct, has 1) object contents and 2) a type. This is not an object system as you can't add objects to it.

```{r}

f <- function() NULL
typeof(f) # [1] "closure"
is.function(f)

typeof(sum) # [1] "builtin" # for primitive functions
is.primitive(sum)

```

use `pryr::otype` to understand what kind of system an object follows

```{r}

pryr::otype(character()) # [1] "base"
pryr::otype(integer()) # [1] "base"
pryr::otype(numeric()) # [1] "base"

pryr::otype(matrix()) # [1] "base"
pryr::otype(list()) # [1] "base"
pryr::otype(c()) # [1] "base"

pryr::otype(factor()) # [1] "S3"
pryr::otype(ts()) # [1] "S3"
pryr::otype(data.frame()) # [1] "S3"
as.Date("1999-01-01") %>% pryr::otype() # [1] "S3"

df <- data.frame(x = 1:5, y = LETTERS[1:5])
pryr::otype(df)
pryr::otype(df$x)
pryr::otype(df$y) # factors are S3

pryr::otype(mean)

```

use `pryr::ftype`

```{r}

pryr::ftype(`%in%`)
pryr::ftype(sum)
pryr::ftype(t.data.frame)
pryr::ftype(t.test) # Tricky!
pryr::ftype(writeLines)
pryr::ftype(unlist)

pryr::ftype(mean)
methods(mean)
pryr::ftype(mean.Date)

```


### S3

#### methods

- query all methods for a S3 generic function
- query all methods S3 function available for a class

```{r, ""}

pryr::ftype(mean)
methods(mean)
methods(class = "Date")

```


#### generic function "UseMethod"

TODO: add example where you pass more arguments

base for `R`'s S3 object system using polymorphism, switching function (function / method dispatch) on the class of the called object.

```{r, "UseMethod"}

asEmpty <- function(object) UseMethod("asEmpty")
asEmpty.integer <- function(object) integer(0)
asEmpty.character <- function(object) character(0)

methods(asEmpty)

1:3 %>% asEmpty
letters %>% asEmpty

```

#### constructor function

example. think the factor and `data.frame` functions

```{r}

foo <- function(x) {
  if (!is.numeric(x)) stop("X must be numeric")
  structure(list(x), class = "foo")
  }

```

#### Use-case

adding custom classes to objects and using polymorphism

```{r}

mkNameFeatures <- function(obj) UseMethod("mkNameFeatures")
mkNameFeatures.se <- function(names) data.frame(name = names, namesCount = names %>% strsplit(" ") %>% purrr::map_dbl(length))
mkNameFeatures.ja <- function(names) data.frame(name = names, hasSpace = names %>% grepl(pattern = " ") )

initClass <- function(obj, cls) obj %>% structure(class = cls %>% append(class(obj)) %>% unique)

jaNames <- c("田中角栄", "杉本 レイ", "上杉謙信") %>% initClass("ja")
seNames <- c("Emil Erik Alexis Rehnberg", "stella", "balle AF raggare") %>% initClass("se")

jaNames %>% mkNameFeatures
seNames %>% mkNameFeatures

jaNames %>% mkNameFeatures %>% initClass("ja")

```



# Standard Libraries {.tabset .tabset-fade .tabset-pills}

## Overview

Looking at functions from e.g. `base` and `stats` packages that are loaded by default R.

## Allocation

as in allocation of empty objects / place holders

```{r "place holder allocation"}
vector("integer", 6)
vector("double", 7)
vector("character", 8)
vector("logical", 3)
```

## File Manipulation

Working with files. Checking existance, copy `cp`, removing `rm`, appending etc.

```{r "file-manipulation"}
library(magrittr)
tmpDir <- file.path("", "tmp", "r-file-manipulation")
dir.exists(tmpDir)
dir.create(tmpDir)
dir.exists(tmpDir)
list.files(tmpDir) # return empty character atomic vector

filePath1 <- tmpDir %>% file.path("A.md")
filePath2 <- tmpDir %>% file.path("B.md")
filePath3 <- tmpDir %>% file.path("C.md")
cat("file A\n", file = filePath1)
file.exists(filePath2)
cat("file B\n", file = filePath2)
file.exists(filePath2)
file.append(filePath1, filePath2)
file.info(filePath1)
# file.show(filePath1) # opens file, only run if in an interactive session
file.copy(filePath1, filePath3)
file.remove(filePath1, filePath2, filePath3)
file.remove(tmpDir) # works on EMPTY directories
```

## Extract

aka `[`. Use `drop` to keep the matrix.

```{r "base:extract"}
(m <- matrix(1:9, nrow = 3, dimnames = list(c("a", "b", "c"), LETTERS[1:3])))
m[1, ] # first row of m (atomic vector)
m[1, , drop = FALSE] # first row of m still as matrix
```

## replications

replicate elements of vectors and lists

```{r "base::rep"}
rep(1:4, 2)
rep(1:4, each = 2)       # not the same.
rep(1:4, c(2,2,2,2))     # same as second.
rep(1:4, c(2,1,2,1))
rep(1:4, each = 2, len = 4)    # first 4 only.
rep(1:4, each = 2, len = 10)   # 8 integers plus two recycled 1's.
rep(1:4, each = 2, times = 3)  # length 24, 3 complete replications

(x <- LETTERS[1:4] %>% setNames(letters[1:4]) %>% as.factor)
rep(x, times = 2)
rep(x, each = 2)
rep.int(x, times = 2)  # no names
rep(x, len = 10)
rep_len(x, length.out = 10) # no names
```

## iterate over data frame variables

use `base::seq_along` in a for loop to iterate over a data frame's variables / columns

```{r}
(scaledArrests <- USArrests) %>% head
for(i in seq_along(USArrests)){
  scaledArrests[[i]] <- scale(USArrests[[i]] %>% as.numeric) %>% round(3)
}
scaledArrests %>% head
```

be aware that the same thing (minus the row names) can be acheived much more nicely with `sapply`

```{r}
USArrests %>% dplyr::select(-state) %>% sapply(scale) %>% as.data.frame %>% head
```

## reduce function

`base::Reduce`

Common Higher-Order Functions in Functional Programming Languages (similar to Common Lisp's `(reduce)`)

Simple example

```{r}
add <- function(x) Reduce("+", x)
add(1:4)
```

Custom cumulative function. Approximation for pi. `right` argument determines if process starts from the right (or from the left).

```{r}
c(3, 7, 15, 1, 292) %>%
  Reduce(f = function(incr, acc) incr + 1 / acc, right = TRUE)
```

## pattern matching

`base::switch`

exact pattern matching in R

```{r}
centre <- function(x, type) {
  switch(type,
         mean = mean(x),
         median = median(x),
         trimmed = mean(x, trim = .1))
}
x <- rcauchy(10)
centre(x, "mean")
centre(x, "median")
centre(x, "trimmed")
```

## split df to list

e.g. split a `data.frame`, rows / observations to list. NOTE: second argument is a "group" e.g. `1:10`. If you in this case have a data frame with ten rows, each row will be in an entry in the returned list.

```{r "listify observations to iterate over rows"}
data.frame(a = 1:10
           ,b = 11:20
           ,c = 51:60) %>%
  split(seq(nrow(.)))
```

or by some other grouping

```{r}
split(airquality, airquality$Month)
```

## work inside closed data environment

`base::with` and `base::within`

```{r}
with(airquality, { cor(Temp, Ozone, use = "complete") })

aq <- within(airquality, {     # Notice that multiple vars can be changed
    lOzone <- log(Ozone)
    Month <- factor(month.abb[Month])
    cTemp <- round((Temp - 32) * 5/9, 1) # From Fahrenheit to Celsius
    S.cT <- Solar.R / cTemp  # using the newly created variable
    rm(Day, Temp)
})
head(aq)
```

## system information

invoke system calls / command line. check system info and library / package info.

```{r "OS calls"}
system("date")
system("echo $USER")

Sys.getenv()

Sys.info() %>% as.matrix

.libPaths()
```

useful for installing packages to a specific library / folder: `install.packages("magrittr", lib = .libPaths()[2])`

default lib parameter is `.libPaths()[1]` which tends to be user specific

`R_LIBS_USER` sets the personal library for a user.

`Sys.unsetenv("R_LIBS_USER")` can unset this so you use a global library path instead

R home

```{r}
R.home()
Sys.getenv("R_HOME")
```

## Citation

```{r message=FALSE}
require(magrittr)
```

```{r}
citation() %>% toBibtex # defaults to R-documentation
citation("magrittr") %>% toBibtex
```

## package version

check the version of an available package

```{r "packageVersion"}
packageVersion("magrittr")
```

## set reference level

`stats:relevel`

to set reference/base/stardard/default level for e.g. a logistic regression use `relevel` with the `ref` flag

```{r}
warpbreaks %>% str
warpbreaks$tension <- relevel(warpbreaks$tension, ref = "M")
warpbreaks %>% str
```

## scientific notation

print in scientific notation

```{r}
library(magrittr)
1234567890123
123456
123456 %>% format(scientific = TRUE)
```

you can also tinker with the `options(scipen=)` option

## print declared variable

use parentheses to return / print result during variable declaration

```{r}
(a <- 1)
a
```

## dicotomize

discretize continuous variable. (turn continuous variable discrete)

```{r}
Z <- stats::rnorm(10000)
table(cut(Z, breaks = -6:6))

breaks <- c(-3, 1, 5) # one more break than labels
labels <- c("low", "high")
U <- runif(10, min=-3, max=5)
cut(U,
    breaks = breaks,
    labels = labels,
    ordered_result = TRUE, # for ordered factors (not always well supported)
    include.lowest = TRUE)
```

## matching

`match` et. al match on the first character forward.

```{r}
1:10 %in% c(1,3,5,9)
charmatch("m",   c("mean", "median", "mode"))
charmatch("med", c("mean", "median", "mode"))
```

#### using grep

matching character strings/vectors using regexes or character strings. use `value` flag to return the actual matches and not the index.

```{r}
txt <- c("arm", "foot", "lefroo", "bafoobar")
grep("foo", txt)
grep("foo", txt, value = TRUE)
```

matching two exact expressions, either or

```{r}
library(magrittr)
txt <-
  c( "The", "licenses", "for", "most", "software", "are", "designed", "to", "take", "away", "your", "freedom", "to", "share"
    ,"and", "change", "it.", "", "By", "contrast,", "the", "GNU", "General", "Public", "License", "is", "intended", "to"
    ,"guarantee", "your", "freedom", "to", "share", "and", "change", "free", "software", "--", "to", "make", "sure"
    ,"the", "software", "is", "free", "for", "all", "its", "users")
```

returns index

```{r}
txt %>% grep(pattern = "ak")
txt %>% grep(pattern = "ak|an")
```

returns value

```{r}
txt %>% grep(pattern = "ak", value = TRUE)
txt %>% grep(pattern = "ak|an", value = TRUE)
```

returns logical (grep logical)

```{r}
txt %>% grepl(pattern = "ak")
txt %>% grepl(pattern = "ak|an")
```

matching on multiple character strings inside the same character string

```{r}
txt %>% grep(pattern = "s.*ar", value = TRUE)

txt %>% grep(pattern = "ra", value = TRUE)
txt %>% grep(pattern = "n.*ra", value = TRUE)
```

matching on end of string by using dollar sign `$` at the end of the pattern

```{r}
strs <-
  c("Product.Age", "Salesforce.Center", "Category", "Standard.Months",
    "FY15.10", "FY15.11", "FY15.12", "FY16.01", "FY16.02", "FY16.03",
    "FY16.04", "ofc-FY16.05", "ofc-FY16.06", "FY16.07", "ofc-FY16.08", "FY16.09",
    "FY16.10", "FY16.11", "FY16.12", "FY17.01", "FY17.02", "FY17.03",
    "FY17.04", "FY17.05", "FY17.06", "FY15.09.1",
    "FY15.10.1", "FY15.11.1", "FY15.12.1", "FY16.01.1", "FY16.02.1",
    "FY16.03.1", "FY16.04.1", "ofc-FY16.05.1", "ofc-FY16.06.1", "FY16.07.1",
    "FY16.08.1", "FY16.09.1", "FY16.10.1", "ofc-FY16.11.1", "ofc-FY16.12.1",
    "FY17.01.1", "FY17.02.1", "FY17.03.1", "FY17.04.1", "FY17.05.1",
    "FY17.06.1")
strs %>% grep(pattern = "(FY\\d{2}).(\\d{2})$", value = TRUE)
```

matching on start of string by using caret `^` at the start of the pattern

```{r}
strs %>% grep(pattern = "^(FY\\d{2}).(\\d{2})$", value = TRUE)
```

## padding

padding with zeros

```{r}
sprintf("%04d", 20)
```

or less stable

```{r}
sprintf("%02s", c("10", "8", "12"))
```

case with hours and minutes, similar to a clock

```{r "sprintf clock"}
sprintf("%02s:%02s", 0:5, 15)
sprintf("%02s:%02s", 0:5, c(5,25,55)) # not expected
sprintf("%02s:%02s", rep(0:5, each = 3), c(5,25,55)) # expected

minutesAsClock <-
  function(minutes)
    sprintf( "%02s:%02s"
            ,minutes %>% magrittr::divide_by_int(60)
            ,minutes %>% magrittr::mod(60)) # modulo `%%`
(0:12*60+40) %>% minutesAsClock
```

Windows

```{r "sprintf on windows"}

sprintf("%02s", c("10", "8", "12")) # doesn't return what you'd expect on windows >.<
# sprintf("%02d", c("10", "8", "12")) # however this does! using `d`
# not sure it works as intended on Unix systems though

```

## substr

can extract characters based on position

```{r}
c('201611', '201304') %>% substr(5,6)
c('201611', '201304') %>% substr(1,4)
```

Can be working similar to T-SQL's `LEFT`-function. No need to use regex

Can also be thought of as deleting parts of strings / characters.

## numerical characteristics of the machine

`.Machine` has a lot of numerical characteristics

e.g. the rounding error

```{r}
.Machine$double.eps
.Machine
```

## setting names for arrays

use `stats::setNames`

```{r}
setNames( 1:3, c("foo", "bar", "baz") )
setNames( 1:3, nm = c("foo", "bar", "baz") )
setNames( nm = c("foo", "bar", "baz") ) # Special case
```

or `base::structure`; NOTE: you need to name any arguments, unlike `setNames`.

```{r}
structure( 1:3, nm = c("foo", "bar", "baz") ) # don't use `nm`
structure( 1:3, names = c("foo", "bar", "baz") )
# structure( 1:3, c("foo", "bar", "baz") ) # Error in structure(1:3, c("foo", "bar", "baz")) : attributes must be named
```

## scale

for normalizing / standardizing variables in `data.frame`

`scale` takes three arguments, the data, center and scale.

- `center` what you subtract from the values (unless set to `FALSE`). Default is the mean of the column in question.
- `scale` is the demoninator (unless set to `FALSE`). Default is the standard deviation.

```{r "scale function"}

airquality %>% summary

airquality %>% scale %>% summary # NOTE: the means and symmetrics between min/max

```

function to scale columns to `[0,1]` unit range

TODO: add scale function to go between `[-1,1]`

```{r}

unitScale <- function(dta)
  scale(dta
        ,center = dta %>% purrr::map_dbl(min, na.rm = TRUE)
        ,scale = dta %>% purrr::map(range, na.rm = TRUE) %>% purrr::map_dbl(diff)
        )
airquality %>% unitScale %>% summary

```

# Helpers {.tabset .tabset-fade .tabset-pills}

### Overview

Section for good helper functions.

### exclude

delete / remove / exclude / reject elements from a vector

```{r "vector remove element(s)"}
exclude.default <- function(vctr, elements) vctr[! vctr %in% elements]
a <- c("master", "of", "puppets")
a %>% exclude(c("of", "master"))
```

### Curry

perform currying in R with the `functional` package. Or just use

```{r}
Curry <- function(FUN, ...){
  .orig <- list(...)
  function(...) do.call(FUN, c(.orig, list(...)))
}
```

### isEmpty

checks if list / atomic vector (array) is empty

```{r}
library(magrittr)
isEmpty <- function(obj) UseMethod("isEmpty")
isEmpty.default <- function(obj) obj %>% length %>% equals(0)
isEmpty.data.frame <- function(dtf) dtf %>% nrow %>% equals(0)
isEmpty.tbl <- function(obj) obj %>% as.data.frame %>% isEmpty
isEmpty.character <- function(chrs) chrs %>% `==`("")

integer(0) %>% isEmpty
1:3 %>% isEmpty
list() %>% isEmpty
1:3 %>% as.list %>% isEmpty
data.frame(price = integer(0)) %>% isEmpty
data.frame(price = 1:3) %>% isEmpty
"" %>% isEmpty
'' %>% isEmpty
'a' %>% isEmpty
```

### isNonEmpty

checks if list / atomic vector (array) is non-empty

```{r}
library(magrittr)
isNonEmpty <- function(obj) UseMethod("isNonEmpty")
isNonEmpty.default <- function(obj) obj %>% length %>% is_greater_than(0)
isNonEmpty.data.frame <- function(obj) obj %>% nrow %>% is_greater_than(0)
isNonEmpty.tbl <- function(obj) obj %>% as.data.frame %>% isNonEmpty

integer(0) %>% isNonEmpty
1:3 %>% isNonEmpty
list() %>% isNonEmpty
1:3 %>% as.list %>% isNonEmpty
data.frame(price = integer(0)) %>% isNonEmpty
data.frame(price = 1:3) %>% isNonEmpty
```

### asQuoted

put single quotes around the argument. This is useful mainly for the underscore variants of `dplyr` functions.

```{r "asQuoted"}
asQuoted <- function(arg) paste0("'", arg, "'")
set.seed(12)
threshold <- .7
data.frame(a = runif(10)) %>% filter_(paste("a > ", threshold %>% asQuoted))
```

### compact

similar to ruby's `compact`.

```{r}

library(magrittr)
compact <- function(object) { UseMethod("compact") }
compact.default <- function(obj) obj[which(!is.na(obj))]
compact.character <- function(obj) obj[which(obj != "")] %>% compact.default
compact.list <- function(lst){
  lst[lapply(lst,
             function(elm){
               if(elm %>% is.atomic){
                 elm %>% length
               }else{
                 elm %>% nrow
               }
             }) > 0]
}

( list(a = 1,
       b = 2,
       c = integer(0),
       d = character(0),
       e = 3,
       f = data.frame(ii = integer(0), cc = character(0))
       ,g = NA
       ) %>% compact
  )

```

### prepend

add `prepend` (for nice naming) since `append` exists

```{r}
prepend <- function(x, values) { append(values, x) }
append(1:2, 6:8)
prepend(1:2, 6:8)
```

### append / prepend class

```{r "prepend class"}

appendClass <- function(obj, cls) {
  class(obj) %<>% append(cls)
  obj
}

prependClass <- function(obj, cls) {
  class(obj) %<>% prepend(cls)
  obj
}

```

### wordArray

inspired by ruby's `%w`

```{r}
library(magrittr)
wordArray <- function(string){ string %>% strsplit("[[:space:]]+") %>% unlist }
"carat color value" %>% wordArray
```

### dataDivision

Divide by zero with set default

```{r}
library(magrittr)
dataDivision <- function(numerator, denominator, divZero = NA){
  ( numerator / denominator ) %>%
   { ifelse(is.finite(.), ., divZero) }
}
dataDivision(10,2)
dataDivision(10,0)
```

### Diagonal proportion of a matrix

returns if object is square shaped i.e. number of rows and columns are the equal

```{r}
isSquareShaped <- function(object) { equals(ncol(object), nrow(object)) }
```

calculates the diagonal proportion in a square-shaped object

```{r}
library(magrittr)
diagonalProportion <-
  function(object){
    if(object %>% isSquareShaped %>% not){
      warning("calculating the diagonal proportion only applies to squared objects")
      return(0)
    }
    divide_by(object %>% diag %>% sum,
              object %>% sum)}
(diag(3) + 0.5) %T>% print %>% diagonalProportion
```

### Interval

function for creating an interval sequence (integers)

```{r}
library(magrittr)
interval <- function(mid, radius) { mid %>% { seq(. - radius, . + radius) } }
```

### mkdir -p

function for creating deep folder seemlessly

```{r}
library(magrittr)
mkdirp <- function(filepath) { filepath %>% dir.create(recursive = TRUE, showWarnings = FALSE) }
```

### rest

function similar to `rest`/`next`/`pop` in closure. I.e. returns object without the first element.
works with e.g. atomic vectors and lists. (will not work with a `data.frame`)

```{r}
library(magrittr)
rest <- function(stk){
  if(stk %>% length %>% is_greater_than(1) ){
    stk[2:length(stk)]
  }else{
    NULL
  }
}

1:4 %>% rest
1:4 %>% as.list %>% rest
```

### truncateLast

similar to the above `rest` function but for the last object in a stack
works with e.g. atomic vectors and lists. (will not work with a `data.frame`)

```{r}
library(magrittr)
truncateLast <- function(stk){
  if(stk %>% length %>% is_greater_than(1) ){
    stk[1:(length(stk) - 1)]
  }else{
    NULL
  }
}

4:1 %>% truncateLast
4:1 %>% as.list %>% truncateLast
```

### coalesce function

Using `Reduce`, I don't know *why* this works >.<

```{r}
a <- c(1,  2,  NA, 4, NA)
b <- c(NA, NA, NA, 5, 6)
c <- c(7,  8,  NA, 9, 10)

coalesce <- function(...) {
  list(...) %>%
    Reduce(f = function(nxt, acc) {
      i <- which(is.na(nxt))
      nxt[i] <- acc[i]
      nxt})
}
coalesce(a, b, c)
```

based on [stack overflow](http://stackoverflow.com/questions/19253820/how-to-implement-coalesce-efficiently-in-r) the manual index version below is faster.

```{r}
coalesce <- function(...) {
  output <- ..1
  for (element in list(...)[-1]) {
    naIndexes <- which(is.na(output))
    output[naIndexes] <- element[naIndexes]
  }
  output
}
coalesce(a, b, c)
```

### multiple plots

using `ggplot` get multiple plots in the same window

TODO: add to package, would be nice, actually maybe make a helper package w documentation etc.

Multiple plot function

ggplot objects can be passed in `...`, or to plotlist (as a list of ggplot objects)
- `cols`:   Number of columns in layout
- `layout`: A matrix specifying the layout. If present, 'cols' is ignored.

If the layout is something like `matrix(c(1,2,3,3), nrow=2, byrow=TRUE)`,
then plot 1 will go in the upper left, 2 will go in the upper right, and
3 will go all the way across the bottom.


```{r "multiplot function"}

multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
  library(grid)

  # Make a list from the ... arguments and plotlist
  plots <- c(list(...), plotlist)

  numPlots = length(plots)

  # If layout is NULL, then use 'cols' to determine layout
  if (is.null(layout)) {
    # Make the panel
    # ncol: Number of columns of plots
    # nrow: Number of rows needed, calculated from # of cols
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                    ncol = cols, nrow = ceiling(numPlots/cols))
  }

 if (numPlots==1) {
    print(plots[[1]])

  } else {
    # Set up the page
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))

    # Make each plot, in the correct location
    for (i in 1:numPlots) {
      # Get the i,j matrix positions of the regions that contain this subplot
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))

      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                      layout.pos.col = matchidx$col))
    }
  }
}
```

`multiplot` demo

```{r}
library(ggplot2)

p1 <- ggplot(ChickWeight, aes(x=Time, y=weight, colour=Diet, group=Chick)) +
	geom_line() +
	ggtitle("Growth curve for individual chicks")

p2 <- ggplot(ChickWeight, aes(x=Time, y=weight, colour=Diet)) +
	geom_point(alpha=.3) +
	geom_smooth(alpha=.2, size=1) +
	ggtitle("Fitted growth curve per diet")

p3 <- ggplot(subset(ChickWeight, Time==21), aes(x=weight, colour=Diet)) +
	geom_density() +
	ggtitle("Final weight, by diet")

p4 <- ggplot(subset(ChickWeight, Time==21), aes(x=weight, fill=Diet)) +
	geom_histogram(colour="black", binwidth=50) +
	facet_grid(Diet ~ .) +
	ggtitle("Final weight, by diet") +
	theme(legend.position="none") # No legend (redundant in this graph)

multiplot(p1, p2, p3, p4, cols = 2)
```

taken from [R cookbook](http://www.cookbook-r.com/Graphs/Multiple_graphs_on_one_page_(ggplot2)/)


# Object Types {.tabset .tabset-fade .tabset-pills}

## Overview

Section for types like lists, data.frame, factors etc.

## Dates {.tabset .tabset-fade .tabset-pills}

### Overview

section for working with dates

### parsing

parsing dates. `strftime` is a wrapper for `format.POSIXct` and `strptime` for `POSIXct`

```{r date_parsing}
strptime("2011-03-27 01:30:00", "%Y-%m-%d %H:%M:%S") # creates a `POSIXlt POSIXt`
strftime("2011-03-27 01:30:00", "%Y-%m-%d %H:%M:%S") # creates a `character`
strptime("01:30:00", "%H:%M:%S")
strptime("20160927", "%Y%m%d")
strptime(20160927, "%Y%m%d")
strptime(20160927, "%Y%m%d") %>% weekdays
as.POSIXct("2011-03-27 01:30:00")
as.POSIXct("2011-03-27 01:30:00") %>% class
as.POSIXct("2011-03-27 01:30:00") %>% class
as.POSIXlt("2011-03-27 01:30:00") %>% unclass %>% unlist
```

### extracting data

```{r "date data extracting"}
strptime("20160927", "%Y%m%d") %>% weekdays
as.POSIXct("2011-03-27 01:30:00") %>% months
strptime("2011-03-27 01:30:00", "%Y-%m-%d %H:%M:%S") %>% format("%H")
```

### combining date & time

```{r}
library(magrittr)
d1 <- as.Date("1970-10-12")
t1 <- "10:40"
paste(d1, t1) %>% as.POSIXct
```

### weeks

```{r}
today <- Sys.Date()
format(today, format = "%Y-%W")
```

## Factors {.tabset .tabset-fade .tabset-pills}

### Overview

section for working with factors

### generate factor levels

```{r}
gl(2, 3, labels = c("Control", "Treatment"))
```

### drop factor levels

Use the `droplevels` function. Returns the full argument object with pruned factor levels.

```{r}
df <-
  data.frame(letters=letters[1:5],
             numbers=seq(1:5))
levels(df$letters)
subdf <- df %>% subset(numbers <= 3)
levels(subdf$letters)
subdf$letters <- droplevels(subdf$letters)
levels(subdf$letters)
```

### Replace NA with value in factor

Replace NA with value in a factor variables.

```{r}
replaceFactorNAs <- function(fctr, naLevel = "Unknown"){
  ifelse(fctr %>% is.na, naLevel, fctr %>% as.character) %>%
    factor(fctr %>% levels %>% append(naLevel))
}
```

### getFactors

extract factor columns from data frame.

```{r}
library(magrittr)
getFactors <- function(dataFrame){
  dataFrame %>%
    names %>%
    sapply(function(colName){
             dataFrame %>%
               getElement(colName) %>%
               is.factor
    }) %>%
    dataFrame[,.]
}
exams <-
  expand.grid(res = c("Y", "N"),
              score = 1:4,
              id = LETTERS[1:5])
exams %>% head %T>% print %>% getFactors
```

### ordered factors

with ordered / ordinal factors you can filter using less / greater / larger / smaller than.

```{r "ordered factors"}
mp <- structure(list(pa = c("M 1", "M 2", "M 3")
                     ,fy = c("FY17.05", "FY17.04", "FY17.03"))
                , .Names = c("pa", "fy")
                , row.names = c(NA, -3L)
                , class = "data.frame")
mp %<>%
  dplyr::mutate( pa = pa %>% ordered
                ,fy = fy %>% as.ordered)
mp %>% str
mp %>% dplyr::filter(pa > "M 1")
```

### transforming factors

```{r}
library(magrittr)
levels <- 1:4
factorVariable <- factor(levels,
                         labels=levels+10)
factorVariable %>% as.integer # returns underlying levels
factorVariable %>% as.character # returns labels
```

### compact factor

to weed out empty factors

```{r "compact factor"}
asCompactFactor <- function(fctr) fctr %>% as.character %>% as.factor
(nonCompactFactor <- factor(c("A", "C", "A"), levels = LETTERS[1:5]))
nonCompactFactor %>% asCompactFactor

```

### edit factors

alter / remap factor-levels

- use a look-up table
- `remapFactor` function is MUCH faster for large factors
- NOTE: the mapping is reversed in the two approaches
- NOTE: use `as.character` in remapFactor2 to avoid serious bugs

remapFactor code taken from [vignette](https://cran.r-project.org/web/packages/gdata/vignettes/mapLevels.pdf)

```{r}
(fvar <-
  gl(n = 3, # number of levels
     k = 2, # number of replications
     labels = c("Control", "Treatment1", "Treatment2")))

remapFactor <- function(fctr, fctrMap){ levels(fctr) <- fctrMap ; fctr }
remapFactor2 <- function(fctr, fctrMap){ fctr %>% as.character %>% fctrMap[.] %>% unname %>% as.factor }

factorLookUpMap <-
  list(
        "trmt" = "Treatment1"
       ,"trmt" = "Treatment2"
       ,"ctrl" = "Control"
       )
fvar %>% remapFactor(factorLookUpMap)

factorLookUpMap2 <-
  c(
     "Treatment1" = "trmt"
    ,"Treatment2" = "trmt"
    ,"Control"    = "ctrl"
    )
fvar %>% remapFactor2(factorLookUpMap2)

```

or use `factor` to reorder factor levels

```{r}
set.seed(12)
ist <- InsectSprays$spray %>% sample %>% head(10) %>% as.character
ist %>% factor(levels = c("F", "D", "E", "A", "B"))
.Last.value %>% table
```

or when it's more complicated

```{r}
c(1, 0, 0, NA, 1) %>%
  factor(labels = c("N", "Y")) %>%
  factor(c("N", "Y", "U"))
```

### rm one level factors

from a data frame you might want to remove columns that are factors with only one level (since there's no variance). (compact factor / sparse factor)

```{r "rmOneLevelFactors"}
isOneLevelFactor <- function(vctr) {(vctr %>% is.factor) & (vctr %>% levels %>% length %>% equals(1))}
rmOneLevelFactors <- function(dta)
  dta %>%
    purrr::map(function(col) if(col %>% isOneLevelFactor) NULL else col ) %>%
    compact %>% as.data.frame

data.frame(a = gl(2,2), b = gl(1,4), d = 1:4) %>% rmOneLevelFactors

```

## Data Frames {.tabset .tabset-fade .tabset-pills}

### Overview

section for working with `data.frame`

### combining

to concatenating / combining / binding data frames

```{r "combining data frames"}

d1 <- data.frame(a = 1, b = 2)
d2 <- data.frame(b = 1, a = 2)
rbind(d1, d2) # returns expected combined data.frames, names are intact

```


## Lists {.tabset .tabset-fade .tabset-pills}

### Overview

section for working with lists

### Subsetting

Subset / Select specific elements in a list.

Use single square brackets and use numerical index >.<

TODO: Add Helper for this??

```{r "list select"}
Harman23.cor[c(2,3)] %>% str
```

using names to exclude / include list elements

```{r "named list subsetting"}
include.list <- function(vctr, elements) vctr[vctr %>% names %in% elements]
exclude.list <- function(vctr, elements) vctr[! vctr %>% names %in% elements]
l1 <-
  list(
     qsec = 2
    ,mpg = 3
    ,disp = 4
    ,vs = 5
    ,wt = 6
    ,am = 7
   )

l1 %>% include(c("qsec", "vs", "am")) %>% str
l1 %>% exclude(c("qsec", "vs", "am")) %>% str
```

### combining

adding / appending / pushing / concatenating lists can be done with `append`

```{r "combing lists"}

la <- list(a = 1:2, rr = runif(4))
lb <- list(b = 3)
append(la, lb)

```

same with tibble

```{r}

library(tibble)
la <- lst(a = 1:2, rr = runif(4))
lb <- lst(b = 3)
append(la, lb)

```

lists with overlapping object names will still merge (not overwrite) and extraction gets iffy

```{r "combining list gotcha"}

la <- list(a = 1:2, rr = runif(4))
lr <- list(rr = 3)
(ll <- append(la, lr))
ll[["rr"]] # grabs the first object with corresponding name. use index to access the later "rr"

```

tibble is consistent with list

```{r}

la <- tibble::lst(a = 1:2, rr = runif(4))
lr <- tibble::lst(rr = 3)
(ll <- append(la, lr))
ll[["rr"]]

purrr::update_list(la, !!! lr) # use purrr::update_list with splicing !!! / UQS()

```



### multi-layer lists

build larger lists and iterate with them

```{r "multi-layer lists"}

library(magrittr)
dates <-
  list( day1 = "2015-12-12"
       ,day2 = "2016-07-12"
       ,day3 = "2014-01-03") %>%
    sapply(as.Date, simplify = FALSE)
wkdays <- dates %>% sapply(weekdays, simplify = FALSE)
mnths <- dates %>% sapply(months, simplify = FALSE)

```

How do we combine these in one named list?

```{r}

mapply(FUN = function(dta, name){
         list( dta = dta
              ,weekday = wkdays[[name]]
              ,month = mnths[[name]]
              )
       }
  ,dates
  ,names(dates)
  ,SIMPLIFY = FALSE
  )

```

or even better

```{r}

dd <- dates %>% sapply(function(dta){list(date = dta)}, simplify = FALSE)
dd %>%
  sapply(function(lst){
           append(lst, list(weekday = lst$date %>% weekdays))
          }, simplify = FALSE)

```

or more succinct if it makes sense

```{r}

dates %>%
  sapply(function(dta){
           list( date = dta
                ,weekday = dta %>% weekdays)
          }, simplify = FALSE)

```

use tibble and purrr

```{r "multi-layers with purrr and tibble"}
library(tibble)
library(purrr)
dates %>%
  as.tibble %>%
  map(function(date){
        lst( date = date
            ,weekday = date %>% weekdays)
       }) # %>% as.tibble

```

multiple arguments to map does NOT work as expected

```{r}

dates %>%
  as.tibble %>%
  map(dates %>% names
      ,.f = function(date, dayLabel){
        lst( date = date
            ,month = date %>% months
            ,label = dayLabel)
      })

```

mix nested lists and vectors

```{r}
datesList <-
  dates %>%
    map(function(date){
      lst( date = date
          ,meta = lst( weekday = date %>% weekdays)
         )
     })

datesList %>% map(function(ddt){
    ddt$meta %<>% append(list(month = ddt$date %>% months))
    ddt
   })

```

### iteration and names

preserving names in the returning list

```{r "lists and preserving names"}
library(magrittr)
dayList <- list(Day1 = 12, Day2 = 34, `Day 3` = 23)
sapply(dayList, sqrt, simplify = FALSE)
```

accessing list item name inside function

```{r "list item access from inside function"}

dayPrinter <- function(dta, listName) paste0(listName, ": ", dta)
Map(dayPrinter
    ,dayList
    ,names(dayList))

```

or with `mapply`, since `Map` is directly using `mapply`

```{r}

mapply(FUN = dayPrinter, dayList, names(dayList), SIMPLIFY = FALSE)

```

use a function to remove the boiler plate

```{r}

mapplyWithNames <- function(lst, FUN) mapply(FUN = FUN, lst, names(lst), SIMPLIFY = FALSE)
dayList %>% mapplyWithNames(dayPrinter)

```


# Package Usages {.tabset .tabset-fade .tabset-pills}

## Overview

Section for different packages. Giving simple or typical use-cases.

## yaml {.tabset .tabset-fade .tabset-pills}

Overview of the yaml pacakge

### Usage

Very useful for config files e.g. with credentials. The default argument to `read_yaml` is a file path to yaml.

```{r }

yaml::read_yaml(text = "{ name: John Smith, Age: 33 }")

```


## glue {.tabset .tabset-fade .tabset-pills}

### Overview

Section for string interpolation package

### glue

main function for string interpolation similar to how it works in Ruby.

```{r}
abc = "DADDY"
glue::glue("Who's your {abc}?")

name <- "Fred"
age <- 50
anniversary <- as.Date("1991-10-12")
glue::glue(
  'My name is {name},',
  'my age next year is {age + 1},',
  'my anniversary is {format(anniversary, "%A, %B %d, %Y")}.')

# single braces can be inserted by doubling them
glue::glue("My name is {name}, not {{name}}.")

# Named arguments can also be supplied
glue::glue(
  'My name is {name},',
  ' my age next year is {age + 1},',
  ' my anniversary is {format(anniversary, "%A, %B %d, %Y")}.',
  name = "Joe",
  age = 40,
  anniversary = as.Date("2001-10-12"))

# `glue_data()` is useful in magrittr pipes
mtcars %>% glue::glue_data("{rownames(.)} has {hp} hp")

# Alternative delimiters can also be used if needed
one <- "1"
glue::glue("The value of $e^{2\\pi i}$ is $<<one>>$.", .open = "<<", .close = ">>")
```

### quoting

there's `single_quote`, `double_quote` and `backtick` functions for surrounding the argument with the select quote. Also there's `collapse` that works similar to how it works in `paste`.

```{r}
x <- 1:5
glue::glue('Values of x: {glue::collapse(glue::backtick(x), sep = ", ", last = " and ")}')
```

## devtools {.tabset .tabset-fade .tabset-pills}

### Overview

Section for devtools basics and qwirks. Use it for creating packages, downloading them etc.

### install.github

install packages from github

```{r "devtools::install_github", eval = FALSE}
devtools::install_github("cran/purrr")
devtools::install_github("hadley/stringr")
devtools::install_github("rstudio/rmarkdown")
devtools::install_github("rstudio/rmarkdown", ref = '00b58f6') # v1.5
devtools::install_github(c("cran/cli", "cran/utf8", "cran/pillar", "cran/tibble"))
devtools::install_github("rstudio/rmarkdown", ref = '79a741b') # v1.1
```

### create

generates a R-package template

```{r "devtools::create", eval = FALSE}
devtools::create("path/to/package/pkgname")
devtools::create("~/code/R/packages/niftytools")
```

### add imports

adds `magrittr` dependency to import field to DESCRIPTION file in a package

```{r "devtools::use_package", eval = FALSE}
devtools::use_package("magrittr")
```

### document

writes roxygen2 documentation, writing NAMESPACE and .Rd (documentation) files.

```{r "devtools::document", eval = FALSE}
devtools::load_all()
devtools::document()
```

### load package contents

Use `devtools::load_all()` function to load the package during development.

### using testthat

Run `devtools::use_testthat()` to set-up testing during package development.

### run tests

Run `devtools::test()` to run the package tests.

## testthat {.tabset .tabset-fade .tabset-pills}

### Overview

Section for all things hadley's testing package `testthat`.

### tests

No news is good news. successful test silently return the first argument.

- `expect_identical` tests with `identical`
- `expect_equal` tests with `all.equal`
- `expect_equivalent` tests with `all.equal` and `check.attributes = FALSE`
- with multiple assertions in the same `it`-block, all assertions are tested and any false assertion gets reported

```{r "testthat tests"}
# library(testthat)

# testing equality
a <- 10
testthat::expect_equal(a, 10)
testthat::expect_equal(10L, 10)
# Use expect_equal() when testing for numeric equality
testthat::expect_equal(sqrt(2) ^ 2, 2)

# testing logical return
testthat::expect_true(2 == 2)
testthat::expect_false(2 != 2)

# greater/less than tests
a <- 9
testthat::expect_lt(a, 10)
testthat::expect_lte(a, 9)
testthat::expect_gt(a, 8)
testthat::expect_gte(a, 9)

# regex match test
testthat::expect_match("Testing is fun", "f.n")

# type and class tests
testthat::expect_null(NULL)
testthat::expect_type(1, "double") # type test
testthat::expect_is(1, "numeric")  # class test

testthat::expect_length(1:10, 10)

# output expectations
testthat::expect_output(str(mtcars), "32 obs")
testthat::expect_output(str(mtcars), "11 variables")
testthat::expect_output(str(mtcars), "11 VARIABLES", ignore.case = TRUE) # extra arguments are sent to grepl
# -- Messages
f <- function(x) {
  if (x < 0) message("*x* is already negative")
  -x
}
testthat::expect_message(f(-1))
testthat::expect_message(f(-1), "already negative")
testthat::expect_message(f(1), NA)
# -- Warnings
f <- function(x) {
  if (x < 0) warning("*x* is already negative")
  -x
}
testthat::expect_warning(f(-1))
testthat::expect_warning(f(-1), "already negative")
testthat::expect_warning(f(1), NA)
## -- Errors
f <- function() stop("My error!")
testthat::expect_error(f())
testthat::expect_error(f(), "My error!")
## -- Silence
# doesn't check silent/invisible return, but rather if there's messages/warnings/errors being output
testthat::expect_silent("123")

# failure expectation
testthat::expect_failure(testthat::expect_identical(1, 1L))
```

### stubs

implemented with `testthat::with_mock`

```{r "with_mock; NOTE: Can't mock functions in base packages anymore", eval = FALSE}
testthat::with_mock(
   all.equal = function(x, y, ...) TRUE
  ,testthat::expect_equal(2 * 3, 4) # silently returns the return of the first argument
  ,.env = "base"
) # return the output of the last expression

testthat::with_mock(
   `base::identical` = function(x, y, ...) TRUE
  ,testthat::expect_identical(3, 5)
  ,testthat::expect_identical(TRUE, 6)
)
```

## dplyr {.tabset .tabset-fade .tabset-pills}

### Overview

Section for dplyr basics and qwirks

### programming

evaluating variable names and values. This is for dplyr 6.0+

```{r "dplyr programming"}

dta <- data.frame(
   v1 = sample(5,10, replace = TRUE)
  ,v2 = sample(5,10, replace = TRUE)
)

```

dplyr filter is a little bit qwirkier.

Note the parentheses around the quo eval, it's required in this situation.

```{r}

myVar <- rlang::quo(v1)
dta %>% dplyr::filter((!!myVar) == 1)

```

while not required in this situation

```{r}

myVal <- 1
dta %>% dplyr::filter(v1 == !!myVal)

```

not that col is a string / character

```{r}

which_col <- "v1"
which_val <- 1

dta %>% dplyr::filter((!!rlang::sym(which_col))==which_val)

```

OR

```{r}

dta %>% dplyr::filter(UQ(rlang::sym(which_col))==which_val)

```

old syntax still works

```{r}

dta %>% dplyr::filter_(.dots= paste0(which_col, "== ", which_val))

```

### subsetting rows

observations

```{r "dplyr filter"}
dplyr::filter(iris, Sepal.Length > 7)
```

remove duplicate rows

```{r "dplyr distinct"}
iris %>% str
dplyr::distinct(iris) %>% str
```

Randomly select fraction of rows.

```{r "dplyr sample_frac"}
dplyr::sample_frac(iris, 0.05, replace = TRUE)
```

Randomly select n rows.

```{r "dplyr sample_n"}
dplyr::sample_n(iris, 10, replace = TRUE)
```

Select rows by position.

```{r "dplyr slice"}
dplyr::slice(iris, 10:15)
```

Select top n rows based on some column

```{r "dplyr top_n"}
CO2 %>% top_n(2, uptake)  # select top 2 observations with highest uptake
CO2 %>% top_n(2, conc)    # conc has a lot of ties
CO2 %>% top_n(2, -uptake) # select top 2 observations with LOWest uptake
```

### selects

- different ways of using `dplyr::select`
- use array with `one_of` to select columns from data frame

```{r}
library(magrittr); library(dplyr)
airquality %>% head %>% select(Temp, Wind)
airquality %>% head %>% select_("Temp", "Wind")
var1 <- "Temp"; var2 <- "Wind"
airquality %>% head %>% select_(var1, var2)

wordArray <- function(string){ string %>% strsplit("[[:space:]]+") %>% unlist }
columns <- "Temp Wind" %>% wordArray
airquality %>% head %>% select(one_of(columns))
# airquality %>% head %>% select_(one_of(columns)) # ERRORS: can not find function `one_of`
```

### filter

If you get weird time-series results then standard library `filter` is most likely first in the namespace.

Use `%in%` or `is_in` for filtering on multiple items.

```{r}
library(dplyr)
set.seed(12)
threshold <- 0.7
dts <- data.frame(a1 = runif(20) %>% round(1), b1 = runif(20) %>% round(1))
dts %<>% mutate(b1Lvl = b1 %>% cut(breaks =c(0, threshold, 1), include.lowest = TRUE, labels = c("low", "high")) )

dts %>% filter(a1 %>% is_in(c(0, 0.2)))
dts %>% filter(a1 %in% c(0, 0.2))

dts %>% filter_("b1 > .7")
dts %>% filter_(paste("b1 >", threshold))
dts %>% filter_("b1Lvl == 'high'")
asQuoted <- function(arg) paste0("'", arg, "'")
b1Level <- "high"
dts %>% filter_(paste("b1Lvl == ", b1Level %>% asQuoted))
```

## pwr {.tabset .tabset-fade .tabset-pills}

### Overview

Overview for the `pwr` package, package for calculating statistical power.

### pwr.chisq.test

Chi-square test power calculations.

```{r "chi-square test power", eval = FALSE}
pwr::pwr.chisq.test(
  w = NULL,  # Effect size - 0.1, 0.3, and 0.5 represent small, medium, and large effect sizes respectively
  N = NULL,  # Sample size
  df = NULL, # degrees of freedom
  sig.level = 0.05,
 power = NULL
  )
```

```{r "chi-square test power cont."}
pwr::pwr.chisq.test(w = .5, N = 145, df = 1)
```

## rlang {.tabset .tabset-fade .tabset-pills}

### Overview

Overview for the `rlang` package. `quo` function is also imported by `dplyr`.

### quo

```{r "dplyr and quo"}
rlang::quo(paste(letters[1:3], collapse = "+")) # doesn't work
rlang::quo(!!paste(letters[1:3], collapse = "+"))
```

no magrittr piping

```{r, eval = FALSE}
paste(letters[1:3], collapse = "+") %>% quo(!!.)
```

sending

```{r}
sortArg <- rlang::quo(cyl)
mtcars %>% dplyr::select(mpg, cyl, disp) %>% dplyr::arrange(!!sortArg)

sortArg <- rlang::quo(cyl %>% (dplyr::desc))
mtcars %>% dplyr::select(mpg, cyl, disp) %>% dplyr::arrange(!!sortArg)
```

doesn't work

```{r, eval = FALSE}
sortArg <- "cyl %>% (dplyr::desc)"
mtcars %>% dplyr::select(mpg, cyl, disp) %>% dplyr::arrange(!!rlang::quo(sortArg))
```

## RSQLServer {.tabset .tabset-fade .tabset-pills}

### Overview

Overview for the `RSQLServer` package, package used for connecting R and MS SQL Server

### workflow

```{r "RSQLServer workflow", eval = FALSE}
oboedbConnect <- function() {
  DBI::dbConnect(RSQLServer::SQLServer()
                 ,server = "OBOE"
                 ,database = "Oboe"
                 ,file = "~/sql.yml" # file for sql settings and credentials
                 )
}
```

SELECTS and grab data

```{r, eval = FALSE}
n <- 5L
statusQuery <- dbplyr::build_sql("SELECT TOP ", n, " * FROM [Oboe].[dbo].[BookingStatus_lkp]")
DBI::dbGetQuery(oboedbConnect(), statusQuery)
```

complex SELECTS

```{r, eval = FALSE}
complicatedSelectQuery <- "
DECLARE @nstat INTEGER = (SELECT 7);

SELECT TOP (@nstat) *
FROM [Oboe].[dbo].[BookingStatus_lkp]
"
DBI::dbGetQuery(oboedbConnect(), complicatedSelectQuery)
```

INSERTs with `dbExecute`

```{r, eval = FALSE}
DBO::dbExecute(oboedbConnect(), "INSERT INTO cars (speed, dist) VALUES (1, 1), (2, 2), (3, 3);")
```

example sql settings file

```
OBOE:
    server: 10.XX.Y.ZZZ
    port: &port 1433
    type: &type sqlserver
    user: etr
    password: &pass blabla
DW:
    server: 10.RR.S.TTT
    port: *port
    type: *type
    user: etr2
    password: *pass
```

## tibble {.tabset .tabset-fade .tabset-pills}

### Overview

Overview of tibble package.

### use-case

use instead of `data.frame` since it's more performant and less qwirky

- A `tibble` never changes the input type.
    - No more worry of characters being automatically turned into strings.
- A `tibble` can have columns that are lists.
- A `tibble` can have non-standard variable names.
    - can start with a number or contain spaces.
    - To use this refer to these in a backtick.
- It only recycles vectors of length 1.
- It never creates row names.

```{r "tibble"}
library(tibble)

a <- 1:5
tibble(a, b = a * 2) # similar to data.frame

lst(n = 5, x = runif(n))
```

`tibble` never coerces its inputs

```{r}
str(tibble(letters))
str(tibble(x = list(diag(1), diag(2))))
```

or munges column names

```{r}
tibble(`a + b` = 1:5)
```

With the SE version, you give it a list of formulas/expressions

```{r}
tibble_(list( x = ~1:10
             ,y = quote(x * 2)))
```

## purrr {.tabset .tabset-fade .tabset-pills}

### Overview

Overview of purrr package.

### use-case

- `map()` makes a list.
- `map_lgl()` makes a logical vector.
- `map_int()` makes an integer vector.
- `map_dbl()` makes a double vector.
- `map_chr()` makes a character vector.

```{r "purrr use case"}
library(purrr)

mtcars %>%
  split(.$cyl) %>% # from base R
  map(~ lm(mpg ~ wt, data = .)) %>%
  map(summary) %>%
  map_dbl("r.squared")
```

```{r "two level split, 2nd split only on half"}
dff <-
  data.frame(F1 = gl(2,3, labels = c("a", "b"))
             ,F2 = gl(3,2, labels = LETTERS[1:3])
             ,r1 = rnorm(6)
             )
stratDff <- dff %>% split(.$F1)
stratDff$b %<>% split(.$F2)

```

### map

```{r "purrr::map"}
library(purrr)

sums <-
  mtcars %>%
    split(.$cyl) %>% # from base R
    map(~ lm(mpg ~ wt, data = .)) %>%
    map(summary)

sums %>% map_dbl("r.squared")
sums %>% map_chr("r.squared")
```

to add names use pipe `set_names` after

```{r "purrr::map with names"}
sums %>% map("call") %>% purrr::set_names(paste0("CylinderCnt_", c(4,6,8)))
```

`.id` argument is available for `purrr::map_df`

```{r ".id argument only works with map_df"}

purrr::rerun(5, rnorm(100)) %>%
  set_names(paste0("sim", 1:5)) %>%
  purrr::map_df(~ dplyr::data_frame(value = .x, step = 1:100), .id = "simulation")

```

### imap - map with index / names

iterate with names for `map`, `walk` etc.
`imap(lst, .f)` is short for `map(lst, names(list), .f)`

```{r}

c(mohammad = "^moh$|^muhamed$|^mohamed$",ahmad = "^ahmed$|^ahmat$") %>% purrr::imap(~ paste0("regex from ", .y, ": ", .x))

c(mohammad = "^moh$|^muhamed$|^mohamed$",ahmad = "^ahmed$|^ahmat$") %>% purrr::imap_chr(~ paste0("regex from ", .y, ": ", .x))

```



### partial

Curry a function. i.e. prefill arguments for a function. by default the prefilled arguments are evaluated lazily.

```{r "purrr::partial"}
compact1 <- function(x) discard(x, is.null)
compact2 <- partial(discard, .p = is.null)
```

and the generated source code is very similar to what we made by hand

```{r}
compact1
compact2

f <- partial(runif, n = rpois(1, 5))
f
f()
f()
```

You can override this by saying .lazy = FALSE

```{r}
f <- partial(runif, n = rpois(1, 5), .lazy = FALSE)
f
f()
f()
```

### prepend

companion to `append`

```{r purrr::prepend}
x <- as.list(1:3)
x %>% append("a")
x %>% prepend("a")
```

### generate random discrete numbers

```{r "purrr::rdunif"}
rdunif(1e3, 10) %>% table
```

### reduce

similar to `base::Reduce`. except the argument order makes sense :D

```{r "purrr::reduce"}
1:10 %>% reduce(`*`)
```

### rerun

rerun some code and `n` times and output the result in a list of length `n`

```{r "purrr::rerun"}

3 %>% rerun(rnorm(5))
2 %>% rerun(nrm = rnorm(3), unif = runif(3))

```

### dealing with failure

use `safely` function.
usage is very restrictive!

```{r "purrr::safely"}
safe_log <- safely(log)
safe_log(10)
safe_log("a")

loggedList <- list(1, 10, "a") %>% map(safely(log))
loggedList %>% str

loggedList %>% transpose %>% str
```

### scalar predicates

predicates for an object being of a type and being scalar (i.e. length 1)

- `is_scalar_list(x)`
- `is_scalar_atomic(x)`
- `is_scalar_vector(x)`
- `is_scalar_numeric(x)`
- `is_scalar_integer(x)`
- `is_scalar_double(x)`
- `is_scalar_character(x)`
- `is_scalar_logical(x)`

### setting names

snake case wrapper for `setNames`, tweaked defaults and stricter argument checks

```{r "purrr::set_names"}
set_names(1:4, letters[1:4])
set_names(letters[1:4]) # names is defaulting to the vector itself
```

### splice

function for appending to a list, or `splicing` all arguments into a list.

probably `update_list` is preferred though due to the edge-case shown below

```{r}
library(purrr)
```

want to an arg3 to inputs

```{r}
inputs <- list(arg1 = "a", arg2 = "b")
```

below is not sufficient

```{r}
list(inputs, arg3 = c("c1", "c2")) %>% str()
c(inputs, arg3 = c("c1", "c2")) %>% str()
```

use `splice`. `splice()` concatenates the elements of inputs with arg3

```{r}
splice(inputs, arg3 = c("c1", "c2")) %>% str()
```

edge-case when two elements get named the same

```{r}
splice(inputs, arg2 = c("c1", "c2")) %>% str()
```

use update_list instead

```{r}
update_list(inputs, arg2 = c("c1", "c2")) %>% str()
```

you can append but it's clunky

```{r}
append(inputs, list(arg3 = c("c1", "c2")))
```

or worse with a complete statement, mutating inputs

```{r}
inputs[["arg3"]] <- c("c1", "c2")
inputs
```

### split/order/sort by list components

DEPRECATED see [github release](https://github.com/tidyverse/purrr/tree/v0.2.3)

```{r "purrr::split_by", eval = FALSE}

(l1 <- transpose(list(x = sample(4), y = letters[1:4])))
l1 %>% order_by("x") # strings can't be used anymore
l1 %>% sort_by("x")

(l2 <- rerun(5, g = sample(2, 1), y = rdunif(5, 10)))
l2 %>% split_by("g") %>% str()
l2 %>% split_by("g") %>% map(. %>% map("y"))

```

### transpose

transpose a list to change order of indexing of a symmetric list

```{r "purrr::transpose"}
safe_log <- safely(log)
loggedList <- list(1, 10, "a") %>% map(safe_log)
loggedList %>% str
loggedList %>% transpose %>% str
```

### type predicates

- `is_list(x)`
- `is_atomic(x)`
- `is_vector(x)`
- `is_numeric(x)`
- `is_integer(x)`
- `is_double(x)`
- `is_character(x)`
- `is_logical(x)`
- `is_null(x)`
- `is_function(x)`

### update a list

modification of a list. similar to `splice` but can overwrite elements with same name.

```{r "purrr::update_list"}

inputs <- list(arg1 = "a", arg2 = "b")
purrr::update_list(inputs, arg2 = c("c1", "c2")) %>% str()

```

instead of using splice

```{r}
splice(inputs, arg2 = c("c1", "c2")) %>% str()

purrr::update_list(inputs, arg3 = c("c1", "c2")) %>% str()

SpecialArguments <- list(arg1 = c("c1", "c2"))
```

make it a deeper list

```{r}
purrr::update_list(inputs, SpecialArguments = SpecialArguments) %>% str
```

update list contents with overwrites using splicing !!! or `UQS()`

```{r}
purrr::update_list(inputs, !!! SpecialArguments) %>% str
```

if you add a list to list though use good old append

```{r}
spArgList <- list(specialArgs = c("c2", "c3"))
```

below broke

```{r, eval = FALSE}
inputs %>% purrr::update_list(spArgList) # will spit out inputs
```

```{r}
inputs %>% append(spArgList)
```

### when

pattern match with conditions using `when`. Matching the 2nd condition.

```{r "purrr::when"}
1:10 %>%
  when( sum(.) <=  50 ~ sum(.)
       ,sum(.) <= 100 ~ sum(.)/2
       ,~ 0
       )
```

matching both the first and 2nd conditions > first condition wins!

```{r}
1:10 %>%
  when( sum(.) <=   x ~ sum(.)
       ,sum(.) <= 2*x ~ sum(.)/2
       ,~ 0
       ,x = 60
       )
```

doesn't match the condition so run default

```{r}
iris %>%
  subset(Sepal.Length > 10) %>%
  when( nrow(.) > 0 ~ .
       ,~ iris %>% head(10)
       )
```

### iterating over multiple vectors simultaneously

use `pmap` to iterate over an arbitrary number of arguments / vectors at the same time.

```{r "purrr::pmap"}
library(purrr)
mu  <- list(5, 10, -3)
sgm <- list(1, 5, 10)
n   <- list(1, 3, 5)
```

with two arguments, you can use `map2`

```{r}
map2(mu, sgm, rnorm, n = 5) %>% str()
```

for more arguments, use a list and `pmap`.

no names means order is significant!!

```{r}
argsList <- list(n, mu, sgm)
argsList %>% pmap(rnorm)
```

use a named list for more comprehensable code

```{r}
argsList <- list(mean = mu, sd = sgm, n = n)
argsList %>% pmap(rnorm)

mu  <- list(5, 10, -3)
sigma <- list(1, 5, 10)
n <- list(1, 3, 5)
args1 <- list(n, mu, sigma)
args1 %>%
  pmap(rnorm) %>%
  str()
```

### walk

use `walk`/`pwalk` to iterate but you're interested in the side-effect instead of the returning value.

e.g. you have a list of plot objects and a list of filenames you want the plots saved to

```{r "purrr:walk"}
library(ggplot2)
plots <- mtcars %>%
  split(.$cyl) %>%
  map(~ggplot(., aes(mpg, wt)) + geom_point())
paths <- stringr::str_c(names(plots), ".pdf")

pwalk(list(paths, plots), ggsave, path = tempdir())
```

### get-attr

read attributes infix

```{r "purrr::get-attr"}
factor(1:3) %@% "levels"
mtcars %@% "class"
```

### empty?

`is_empty` function for vectors or lists

```{r purrr::is_empty}
is_empty(NULL)
is_empty(list())
is_empty(list(NULL))
```

### keep / discard / compact

usually called `filter`/`select`/`find_all` and `reject`/`drop` in other languages. (these names are taken however)

`compact` works well with both lists and vectors

```{r "purrr::keep"}
rep(10, 10) %>%
  map(sample, 5) %>%
  keep(function(x) mean(x) > 6)
```

Or use a formula

```{r}
rep(10, 10) %>%
  map(sample, 5) %>%
  keep(~ mean(.x) > 6)
```

Using a string instead of a function will select all list elements.
Where that subelement is `TRUE.`

```{r}
(x <- rerun(5, a = rbernoulli(1), b = sample(10)))
x %>% keep("a")
x %>% discard("a")
```

filter list on is.integer

```{r}
list(a = 1:3
     ,b = gl(4,3)
     ,d = integer()
     ,e = letters
     ) %>% purrr::keep(is.integer)
```

compact keeps all non-NULL and empty items

```{r}
list( a = 1
     ,b = NULL
     ,c = FALSE
     ,d = "FALSE"
     ,e = integer()
     ,f = character()
     ,g = list(NULL)
     ,h = NA
     ) %>% (purrr::compact)

c(1, NULL, 3, NA, 5) %>% compact
c(1, NULL, 3, NA, 5) %>% purrr::compact()
```

### conditional map

use `map_if` with predicate function and `map_at` with names or index

```{r "purrr::map_if"}
iris %>%
  purrr::map_if(is.factor, as.character) %>%
  str()
```

Specify which columns to map with a numeric vector of positions:

```{r}
mtcars %>% map_at(c(1, 4, 5), as.character) %>% str()
```

Or with a vector of names:

```{r}
mtcars %>% map_at(c("cyl", "am"), as.character) %>% str()
```

### lmap

apply a function to list-elements of a list, and must return a list. the length of the list going in can be of a different size to that one going out.

Let's write a function that returns a larger list or an empty list
depending on some condition. This function also uses the names
metadata available in the attributes of the list-element.

```{r "purrr::lmap"}
library(purrr)
maybe_rep <- function(e) {
  n <- rpois(1, 2)
  out <- rep_len(e, n)
  if (length(out) > 0) {
    names(out) <- paste0(names(e), seq_len(n))
  }
  out
}

x <- list(a = 1:4, b = letters[5:7], c = 8:9, d = letters[10])
```

The output size varies each time we map f()

```{r}
x %>% lmap(maybe_rep)
```

We can apply f() on a selected subset of x

```{r}
x %>% lmap_at(c("a", "d"), maybe_rep)
```

Or only where a condition is satisfied

```{r}
x %>% lmap_if(is.character, maybe_rep)
```

each element in the list is sent in as a list

```{r}
list(a = 1:2, b = letters[1:3]) %>% purrr::lmap(function(eal) { class(eal) %>% print; eal}) %>% invisible
```

making it possible to access the names of the list entries

```{r}
x %>% purrr::lmap(function(e) { names(e) %<>% paste0("_3"); e})
```

and you'll have to dig deeper to e.g. alter the names of the original list entries

```
list(a = c(B = 2, V = 3)) %>% lmap(function(eal) {names(eal[[1]]) %<>% paste0(0); eal})
```

### at depth

apply a function at a set depth of a nested list

`x %>% at_depth(0, fun)` is equivalent to `fun(x)`.

`x %>% at_depth(1, fun)` is equivalent to `map(x, fun)`.

`x %>% at_depth(2, fun)` is equivalent to `map(x, . %>% map(fun))`.

```{r "purrr::modify_depth prep"}

l1 <- list(
  obj1 = list(
    prop1 = list(param1 = 1:2, param2 = 3:4),
    prop2 = list(param1 = 5:6, param2 = 7:8)
    ),
  obj2 = list(
    prop1 = list(param1 = 9:10, param2 = 11:12),
    prop2 = list(param1 = 13:14, param2 = 15:16)
    )
  )

```

In the above list, "obj" is level 1, "prop" is level 2 and "param"
is level 3. To apply `sum()` on all params, we map it at depth 3:

```{r "purrr::at_depth DEPRECATED", eval = FALSE}

l1 %>% at_depth(3, sum)
```

map() lets us pluck the elements prop1/param2 in obj1 and obj2:

```{r}
l1 %>% map(c("prop1", "param2")) %>% str()
```

But what if we want to pluck all param2 elements? Then we need to act at a lower level:

```{r}
l1 %>% at_depth(2, "param2") %>% str()

```

```{r "purrr::modify_depth example"}

l1 %>% purrr::modify_depth(3, sum )
l1 %>% purrr::map(c("prop1", "param2"))
l1 %>% purrr::modify_depth(2, "param2")
l1 %>% purrr::modify_depth(2, ~ log(.$param2) )

```


### has_element

predicate function for checking inclussion / containing of objects in another.

works with lists and vectors. but there's some gotchas.

```{r "purrr::has_element"}

x <- list(1:10, 5, 9.9)
x %>% has_element(1:10)
x %>% has_element(3)
list(1, 2, 3) %>% has_element(2)
c(1, 2, 3) %>% has_element(2)
```

gotchas!

```{r}
1:3 %>% as.list %>% has_element(2)
1:3 %>% has_element(2)
```

instead use

```{r}
1:3 %>% as.list %>% has_element(2L)
1:3 %>% has_element(2L)

```

### invoke

Set of functions facilitating invokation / execution of combinations of functions and parameters

```{r "purrr::invoke"}
# Invoke a function with a list of arguments
invoke(runif, list(n = 10))
# Invoke a function with named arguments
invoke(runif, n = 10)

list("01a", "01b") %>%
  invoke(paste, ., sep = ".")
```

invoke map

```{r "purrr::invoke_map"}
# Invoke a list of functions, each with different arguments
invoke_map(list(runif, rnorm), list(list(n = 10), list(n = 5)))
# Or with the same inputs:
invoke_map(list(runif, rnorm), list(list(n = 5)))
invoke_map(list(runif, rnorm), n = 5)
# Or the same function with different inputs:
invoke_map("runif", list(list(n = 5), list(n = 10)))

# Or as a pipeline
list(m1 = mean, m2 = median) %>% invoke_map(x = rcauchy(100))
list(m1 = mean, m2 = median) %>% invoke_map_dbl(x = rcauchy(100))
```

### lift

composition helper. helps with sending parameters to functions.

- `d`: dots `...`
- `l`: list
- `v`: vector

`lift_dl` means lift dots to list.

```{r "purrr::lift functions"}
### Lifting from ... to list(...) or c(...)
x <- list(x = c(1:100, NA, 1000), na.rm = TRUE, trim = 0.9)
lift_dl(mean)(x)
lift(mean)(x) # lift is alias for lift_dl

# Or in a pipe:
mean %>% lift_dl %>% invoke(x)

# Default arguments can also be specified directly in lift_dl()
list(c(1:100, NA, 1000)) %>% lift_dl(mean, na.rm = TRUE)()

# lift_dl() and lift_ld() are inverse of each other.
# Here we transform sum() so that it takes a list
fun <- sum %>% lift_dl()
fun(list(3, NA, 4, na.rm = TRUE))
# Now we transform it back to a variadic function
fun2 <- fun %>% lift_ld()
fun2(3, NA, 4, na.rm = TRUE)
```

### combinations

read the `cross` help for more complex situations.

```{r "purrr::cross"}
dta <-
  list(id = c("John", "Jane"),
       greeting = c("Hello.", "Bonjour."),
       sep = c("! ", "... ")
       )

combinations <- dta %>% cross()
combinations %>% str
combinations %>% map(lift(paste))
```

### detect

`find`/`detect` in ruby. returns FIRST result based on some data and a predicate.

works for vectors or lists.

```{r "purrr::detect"}
is_even <- function(x) x %% 2 == 0

3:10 %>% detect(is_even)
3:10 %>% detect_index(is_even)

3:10 %>% detect(is_even, .right = TRUE)
3:10 %>% detect_index(is_even, .right = TRUE)

l1 <-
  list( l11 = list(a = 1, b = 2, mes = "1st item")
       ,l12 = list(a = 3, b = 4, mes = "2nd item"))
l1 %>% detect(~ .$b == 4)
```

### every / some

`all`/`any` in Ruby. predicate functions.

```{r "purrr::every"}
x <- list(0, 1, TRUE)
x %>% every(identity)
x %>% some(identity) # returns TRUE on the last predicate check

y <- list(0:10, 5.5)
y %>% every(is.numeric)
y %>% every(is.integer)
```

### flatten

similar to `unlist` but only flattens one layer at a time.

```{r "purrr::flatten"}
(x <- list( rerun(2, sample(4))
           ,rerun(2, sample(4))))
x %>% flatten()
rerun(2, sample(4)) %>% flatten_int()
```

### negate

negate a predicate function

```{r "purrr::negate"}
set.seed(12)
x <- transpose(list(x = 1:10, y = rbernoulli(10)))
x %>% keep("y") %>% length() # keep the list entries with y being TRUE
x %>% keep(negate("y")) %>% length() # same as above but y being FALSE
```

### null-default

inspired from ruby's `||`

```{r "purrr::null-default"}
1 %||% 2
NULL %||% 2
NA %||% 2
integer() %||% 2
```

## tidyr {.tabset .tabset-fade .tabset-pills}

### Overview

Usage basics for this data-munging library

### gather

unpivot a data frame. i.e. turn columns into rows. Will use the WHOLE data-frame, so remove columns you don't want in the final result.

```{r tidyr::gather}
library(magrittr)
stocks <- data.frame(
  time = as.Date('2009-01-01') + 0:7,
  dayNum = 1:8,
  X = rnorm(8, 0, 1),
  Y = rnorm(8, 0, 2),
  Z = rnorm(8, 0, 4)
)
stocks %>% str

stocks %>%
  tidyr::gather( stock # name of the unpivoted data label
                ,price # name of the unpivoted data value
                ,-c(time, dayNum)) # REMOVE the columns you still want aggregates for
```

### spread

pivot a data frame. i.e. turn a factor column with numeric column and return numeric columns for each factor

```{r tidyr::spread}
library(magrittr)
stocks <- data.frame(
   time = as.Date('2009-01-01') + 0:7
  ,dayNum = 1:8
  ,stock = c("X", "Y", "Z")
  ,price = rnorm(3*8)
)
stocks %>% str

stocks %>%
  tidyr::spread( stock  # name of the factor column to pivot
                ,price) # name of the numeric column to pivot
```

## XLConnect {.tabset .tabset-fade .tabset-pills}

### Overview

Section for XLConnect basic workflows.

### installation

```{r "install XLconnect", eval = FALSE}
devtools::install_github('cran/xlconnectjars')
devtools::install_github('cran/xlconnect')
```

### writing excel sheets

```{r "writing xl sheets"}
library(magrittr)
library(XLConnect)
irisTbl <-
  iris %>% with({table(Petal.Width, Species)}) %>% as.data.frame
wbPath <- "/tmp/ss.xlsx"
wb <- loadWorkbook(wbPath, create = TRUE)
wb %>% createSheet(name = "irisTable")
wb %>% appendWorksheet(irisTbl, sheet = "irisTable")
wb %>% saveWorkbook
file.remove(wbPath)
```


## tidytext {.tabset .tabset-fade .tabset-pills}

```{r "installing tidytext", eval = FALSE}

paste0("cran/",
       c("hunspell", "mnormt", "psych", "RWekajars", "RWeka", "Snowball", "SnowballC", "tokenizers"
         ,"janeaustenr", "ISOcodes", "stopwords")) %>%
  c(.,"tidyverse/broom", "juliasilge/tidytext") %>%
  (devtools::install_github)

```

Small example working with text data. In this case Jane Austin books.

add row and chapter numbers.

```{r "small example of working with text data"}

janeaustenr::austen_books() %>% str

(crudeBooks <-
  janeaustenr::austen_books() %>%
    dplyr::group_by(book) %>%
    dplyr::mutate(
      line = row_number()
      ,chapter =
        cumsum(stringr::str_detect(text, stringr::regex("^chapter [\\divxlc]", ignore_case = TRUE)))
     ) %>%
    (dplyr::ungroup)
)

```

### Tokenize your data

restructure data as *one-token-per-row* instead.

```{r "tidytext::unnest_tokens"}

(books <- crudeBooks %>% tidytext::unnest_tokens(word, text))

```


## magrittr {.tabset .tabset-fade .tabset-pills}

### Overview

magrittr stuff. TODO: Add Anki stuff

```{r echo = FALSE}
library(magrittr)
```

### pipe anonymous function

wrap the function in paratheses.

```{r}
f1 <- gl(2,4) %>% sample
f2 <- gl(2,4) %>% sample
table(f1, f2) %>%
  (function(table){ divide_by(table %>% diag %>% sum,
                              table %>% sum) })
# or use special magrittr syntax-sugar
table(f1, f2) %>%
  { divide_by(diag(.) %>% sum,
              sum(.)) }
```

### full package function lookup

to use full package path to function, use paratheses to avoid an exeception

```{r "magrittr package function lookup qwirk"}
c(1, NULL, 3, NA, 5) %>% (purrr::compact)

# there's no problems if you're using paratheses after the function though, so if there's extra params, piping works as usual
c(1, NULL, 3, NA, 5) %>% purrr::compact(.)

# or with extra parameters
zeroPaster <- function(e1, e2) paste0(e2, e1)
list(1:3, letters[1:3]) %>% purrr::pmap(.f = zeroPaster)
# list(1:3, letters[1:3]) %>% purrr::pmap(.l = ., .f = function(e1, e2) paste0(e2, e1)) # also works
1:3 %>% purrr::map2(letters[1:3], zeroPaster)
```

## regexes {.tabset .tabset-fade .tabset-pills}

### Overview

- see `?regex` for more details on Regular Expressions
- can be used for many things. e.g. `strsplit`
- also see `grep` section

### classes

- `[:alnum:]` Alphanumeric characters: `[:alpha:]` and `[:digit:]`.
- `[:alpha:]` Alphabetic characters: `[:lower:]` and `[:upper:]`.
- `[:blank:]` Blank characters: space and tab, and possibly other locale-dependent characters such as non-breaking space.
- `[:cntrl:]` Control characters.  In ASCII, these characters have octal codes 000 through 037, and 177 (‘DEL’).  In another character set, these are the equivalent characters, if any.
- `[:digit:]` Digits: `0 1 2 3 4 5 6 7 8 9`.
- `[:graph:]` Graphical characters: `[:alnum:]` and `[:punct:]`.
- `[:lower:]` Lower-case letters in the current locale.
- `[:print:]` Printable characters: `[:alnum:]`, `[:punct:]` and space.
- `[:punct:]` Punctuation characters: `! " # $ % & ' ( ) * + , - . / : ; < = > ? @ [ \ ] ^ _ \` { | } ~`.
- `[:space:]` Space characters: tab, newline, vertical tab, form feed, carriage return, space and possibly other locale-dependent characters.
- `[:upper:]` Upper-case letters in the current locale.
- `[:xdigit:]` Hexadecimal digits: `0 1 2 3 4 5 6 7 8 9 A B C D E F a b c d e f`.

### Reverse logic

use e.g. `[^0-9]`

```{r}

phoneNums <-
  c("081398664455", "083114605084", "089606983350", "0878 83282686",
"085770714422", "085237141280", "+6285348551672", "0215522687",
"+60109445138", "085791400805", "085795044447", "0811296113",
"082136714653", "082137487870", "081380380959", "081225120000",
"081291551721", "082225633200", "0895348595086", "085604306828",
"082376342537", "081336618616", "081316564350", "08161368240",
"082346633851", "082284374650", "082157122629", "08126063500",
"087821102719", "089618876742", "085888292599", "082292175370",
"082259505304", "085230018600", "081247883268", "081230068486",
"+6282285553134", "087881492265", "087882320948", "081219630831",
"089530803830", "085640996911", "081363930762", "0217891234",
"0895361710777", "085819650926", "0815 27967189", "08567704404",
"081216802488", "085259810046", "081315259249", "082161568728",
"081221410683", "087753142839", "082151785551", "081318762621",
"082226329094", "08121808035", "0878830-15351", "081221626240",
"081214858079", "081802222650", "089514562860", "082346907926",
"082344655309", "085891433990", "081217068606", "081219048706",
"085892290313", "085264074000", "081283644779", "083167604590",
"081372558794", "085340731850", "+6283831611122", "081381034778",
"08129209872", "082130763333", "081310063397", "08138507777",
"0816605861", "085292434342", "081213457989", "089505709376",
"081282341665", "082216363055", "+6285719223582", "081288402780",
"02812345678910", "081916590999", "082232166668", "082111947219",
"085749183589", "0851-0503-6982", "085692342663", "082187471400",
"081296840701", "089517733237", "085228881033", "089505153702",
"0271633565")

phoneNums %>% gsub("[0-9]+", "\\1", .)
phoneNums %>% gsub("[^0-9]+", "\\1", .)

phoneNums %>% gsub("\\d+", "\\1", .)
phoneNums %>% gsub("\\D+", "\\1", .)

```

or with classes

```{r}

phoneNums %>% gsub("[[:digit:]]+", "\\1", .)
phoneNums %>% gsub("[^[:digit:]]+", "\\1", .)

```

### splitting a string

split a string into a atomic vector

```{r}
library(magrittr)
"name bad  good \t heaven" %>% strsplit("[[:space:]]+") %>% unlist
```


### extracting based on pattern

example for extracting numbers from strings

- double escape regex-backslashes
- use `.*` to match before and after the `()`-group

```{r}
files <-
  c("data/attendance/workshop-2016-13-Attendance.tsv"
   ,"data/attendance/workshop-2016-19-Attendance.tsv"
   ,"data/attendance/workshop-2016-18-Attendance.tsv"
   ,"data/attendance/workshop-2015-34-Attendance.tsv"
   ,"data/attendance/workshop-2016-15-Attendance.tsv"
   ,"data/attendance/workshop-2016-04-Attendance.tsv"
   ,"data/attendance/workshop-2016-01-Attendance.tsv"
)
rs <- ".*(\\d{4}[-]\\d+).*"
files %>% sub(rs, "\\1", .)
```

regex in `dir` to list specific files from file system in a directory. similar to `ls` in Unix.

Below:

- Say that `dir` returns the files below.
- Use the `pattern` flag like below to only extract `\d{4}-\d{2}` file names.

```{r, eval = FALSE}

files <-
  c( "total-2016-25-attendance.tsv", "total-2016-26-attendance.tsv"
    ,"total-2016-27-attendance.tsv", "total-2016-28-attendance.tsv"
    ,"total-2016-29-attendance.tsv", "total-2016-30-attendance.tsv"
    ,"total-2016-31-attendance.tsv", "total-2016-32-attendance.tsv"
    ,"total-2016-33-attendance.tsv", "total-2016-34-attendance.tsv"
    ,"total-2016-35-attendance.tsv", "total-2016-36-attendance.tsv"
    ,"total-2016-37-attendance.tsv", "total-holiday-attendance.tsv"
  )

dir(attendanceDir(center)
    ,pattern = paste(classType, "[-](\\d{4}[-]\\d{2}).*", sep = ""))

```

```{r}
timeslots <-
  c("13:05", "14:40", "18:40", "19:40", "12:40", "13:40",
    "18:40", "19:40", "12:40", "14:40", "18:05", "18:40",
    "19:40", "11:40", "12:40", "14:40", "12:40", "13:40",
    "14:40", "16:40", "12:40", "13:40", "18:40", "20:40")
```

Extracting hours

```{r}
timeslots %>%
  gsub("(\\d{2}).*", "\\1", .)
```

Extracting minutes

```{r}
timeslots %>%
  gsub(".*(\\d{2})", "\\1", .)
```

### mulitple patterns

single character

```{r}
FirstName <- c("Ben","Brck","Adam","Molly","Eve")
FirstName %>% grepl(pattern = "[AaEeIi]")
```

multiple character

```{r}
MobileTel <- c("170766666", "18132452345", "138789", "153213", "111342000"
               ,"1772341", "18452342", "19123123", "17089991", "12188123")

MobileTel %>%
  grepl(pattern = "1340|1341|1342|1346|136|150|152|183|184|187|188")
```

### pattern begins with

match on the beginning of the string

```{r}
MobileTel <- c("170766666", "18132452345", "138789", "153213", "111342000"
               ,"1772341", "18452342", "19123123", "17089991", "12188123")

MobileTel %>%
  grepl(pattern = "^1340|^1341|^1342|^1346|^136|^150|^152|^183|^184|^187|^188")
```

# Cheats {.tabset .tabset-fade .tabset-pills}

## Overview

Section for "tricks".

## saving files

`save` saves the R-object(s) to disc and keeps the naming

```{r, eval = FALSE}
save(foo,file="data.Rda")
load("data.Rda")
```

`saveRDS` saves one object and does not keep the name for the object.

```{r, eval = FALSE}
saveRDS(foo, file="data.Rda") # even for S4 objects
bar <- readRDS(file="data.Rda")
```

```{r, eval = FALSE}
write
write.table
write.csv
```

```{r, eval = FALSE}
library(MASS)
mat <- matrix(1:100,nrow=20)
write.matrix(mat,'file.prn',sep = "\t")

mat2 <- as.matrix(read.table("file.prn", as.is = TRUE))
mat2 <- unname(mat2) # make mat2 a true matrix by removing the names
all.equal(mat, mat2) # [1] TRUE
```

```{r, eval = FALSE}
dput; dget # not for S4 objects
```

```{r, eval = FALSE}
dump; source # not for S4 objects
```

## remove names

```{r}
require(magrittr)
dt <- data.frame(V1 = runif(100) %>% round,
                 V2 = runif(100) %>% round)
(crossMatrix <-
  ( with(dt, table(V1,V2))
    %>% as.data.frame(row.names = LETTERS[1:4])
    %>% as.matrix
    )
  )
crossMatrix %>% unname
```

or NULLify either row or column names

```{r}

colnames(crossMatrix) <- NULL
crossMatrix
rownames(crossMatrix) <- NULL
crossMatrix

```

## remove named columns

```{r message=FALSE}
require(magrittr)
```

```{r}
df <- data.frame(x=1:5, y=2:6, z=3:7, u=4:8)
df
df %>% subset(select=-c(z,u))
df[ , -which(names(df) %in% c("z","u"))] # also works but not as nice IMO
```

## insertion in characters

or strings. how to find something similar to the use case of Ruby's `insert`

```{r "character insertion"}
"shjat" %>% gsub('^([a-z]{2})([a-z]+)$',"\\1_\\2",.)

insert <- function(obj, ...) UseMethod("insert")
insert.character <- function(chr, position, insertion){ # only works with \w characters
  chr %>%
    gsub( paste0("^(\\w{", position - 1, "})(\\w+)$")
         ,paste0("\\1", insertion, "\\2")
         ,.)
}
"shjat" %>% insert(3, "_")
```

## time-series convertions

```{r}

devtools::install_github("cran/xts")
library(xts)
price <- structure(list(
  date = c(20070103L, 20070104L, 20070105L, 20070108L, 20070109L,
           20070110L, 20070111L, 20070112L, 20070115L),
  close = c(54.7, 54.77, 55.12, 54.87, 54.86, 54.27, 54.77, 55.36, 55.76)),
  .Names = c("date", "close"), class = "data.frame",
  row.names = c("1", "2", "3", "4", "5", "6", "7", "8", "9"))
price$date %<>% as.character %>% as.Date(format="%Y%m%d")
xts(price$close, price$date)

```

```{r}

library(zoo)
pricez <- read.zoo(text="   DATE  CLOSE
 1    20070103 54.700
 2    20070104 54.770
 3    20070105 55.120
 4    20070108 54.870
 5    20070109 54.860
 6    20070110 54.270
 7    20070111 54.770
 8    20070112 55.360
 9    20070115 55.760
 ")
index(pricez) <- as.Date(as.character(index(pricez)), format="%Y%m%d")
pricez

```

`ts` (time-series) object to `data.frame`

```{r}
library(zoo)
UKgas %>%
  { data.frame(date = as.Date(as.yearqtr(time(.))),
               gas = as.matrix(.)) }
```

## Tables

```{r, eval=FALSE, echo=FALSE}

dependencies <- c("gtools", "gdata", "gmodels")
devtools::install_github(paste("cran", dependencies, sep = "/"))

```

```{r}

library(gmodels)
data(infert, package = "datasets")

CrossTable(warpbreaks$wool,
           warpbreaks$tension,
           chisq = TRUE,
           prop.t = TRUE,
           digits = 2,
           dnn = c("Wool", "Tension"))

CrossTable(infert$education, infert$induced, expected = TRUE)
CrossTable(infert$education, infert$induced, expected = TRUE, format="SAS")
CrossTable(infert$education, infert$induced, expected = TRUE, format="SPSS")
```

## cumulatively sum within a group

Cumulativily sum the `carb` column within each `am`-group

```{r}
library(magrittr)
library(dplyr)
( mtcars
  %>% filter(cyl == "6")
  %>% select(carb, gear, am)
  %>% group_by(am)
  %>% arrange(am)
  %T>% print
  %>% do(cumsum(.["carb"]))
  )
```

## list to data.frame

given a list of data frames. use `Reduce` to bind them togetherq

```{r}
library(magrittr)
library(dplyr)
c(3, 6, 12) %>%
  lapply(function(day){ airquality %>% filter(Day == day) }) %>%
  Reduce(rbind, .)
```

# Comparing objects {.tabset .tabset-fade .tabset-pills}

### Overview

section for comparing objects

### comparing scalars

```{r "comparing scalars"}
1 == 1L           # not strict
identical(1, 1L)  # strict

x1 <- 0.5 - 0.3
x2 <- 0.3 - 0.1
x1 == x2          # FALSE on most machines
all.equal(x1, x2) # TRUE everywhere
```

### comparing vectors

```{r "comparing vectors"}
v1 <- 1:3               # integer
v2 <- c(2-1, 3-1, 4-1)  # double
v1 == v2
all.equal(v1, v2)
identical(v1, v2)

v1s <- c(2,3,1) # order significant, naturally
v1 == v1s
all.equal(v1, v1s)

```

comparing named vectors

```{r "comparing named vectors"}

v3n <- c(`a` = 1, b = 2)
v3ns <- c(b = 2, `a` = 1)
v3n == v3ns
all.equal(v3n, v3ns)

```

### comparing lists

```{r "comparing lists"}

l1 <- list(aa = 1, bb = 2)
l2 <- list(bb = 2, aa = 1)
all.equal(l1, l2, check.attributes = FALSE) # l1 == l2 # not implemented

```

order insignificant list comparison

```{r }

isEquivalentList <- function(list1, list2){
  if(length(list1) != length(list2)) return(FALSE)
  setdiff(list1, list2) %>% isEmpty
}
l3 <- list(bb = 2, aa = 1, cc = 3)
l4 <- list(bb = 2, aa = 31)
l5 <- list(2, 1)
isEquivalentList(l1, l2)
isEquivalentList(l2, l1)
isEquivalentList(l3, l1)
isEquivalentList(l1, l3)

isEquivalentList(l1, l4) # checks values

isEquivalentList(l1, l5) # does not check names

isEquivalentList(v1, v2) # works on vectors too

```

# Plots {.tabset .tabset-fade .tabset-pills}

### Overview

- `ggplot2` will have precedence to base standard library `graphics::plot`

### helpers

`cut_width` is a helper function for easier cutting of numerical values.

```{r ggplot_cut_width}
library(magrittr)
set.seed(12)
(rr <- runif(10, max = 9) %>% round)
rr %>% ggplot2::cut_width(3)
```

### ggplot in functions

with `ggplot` in a function, if you're passing arguments for `aes` use `aes_string`.

```{r}
library(ggplot2)
cars <- within(mtcars, { cyl <- as.factor(cyl) })
renderBoxplot <- function(yvar){
  ggplot(cars, aes_string(x = "cyl", y = yvar)) + geom_boxplot()
}
renderBoxplot("mpg")
renderBoxplot("hp")
```

### Typical

```{r}
library(ggplot2)
(ggplot(airquality,
        aes(x = paste(Month, Day, sep="-"), y = Wind))
  + geom_point(size = .7)
  + geom_line(size = .5, colour = "blue")
  + ylab("Wind (mph)")
  + xlab("Days")
  + ylim(0, 25)   # sets limit on y-axis
  + ggtitle("Wind over time")
  + geom_hline(aes(yintercept = median(Wind)), colour = "red")
  + scale_x_discrete(breaks=NULL)
)
```

### line plot

line plots for e.g. a grouped (group by) data doesn't work since there's only one group, but ggplot expects there to be more than one group and it draws lines within the group. Hence, create a dummy group for ggplot.

typical error: `geom_path: Each group consist of only one observation. Do you need to adjust the group aesthetic?`

```{r}
library(ggplot2)
library(magrittr)
library(dplyr)
( airquality
  %>% group_by(Month)
  %>% summarise(wind = mean(Wind))
  %>% ggplot(aes(y = wind, x = Month, group = 1))
    + geom_point()
    + geom_line()
 )
```

with line groups

```{r}
library(ggplot2)
( ggplot( within(airquality, {Month %<>% as.factor})
         ,aes( x = Day
              ,y = Temp
              ,group = Month
              ,colour = Month)
       )
  + geom_line(size = .5)
  + geom_point(size = .7, fill = "white")
)
```

### non-numeric argument

typical error: `non-numeric argument to binary operator`

this can be due to many reasons! most common is some syntax error, check trailing `+` and paratheses.

### Histogram

```{r}
ggplot(airquality, aes(Wind)) + geom_histogram()
```

### Boxplot

```{r}
library(ggplot2)
(ggplot(mpg, aes(class, hwy))
 + geom_boxplot(outlier.size = 0.5)
)
```

to deal with horizontal labels being to many, you can angle (diagonal, twist, turn, transpose, flip) the labels

```{r}
aq <-
  within(airquality, {
           mthStg = Day %>% cut(breaks = c(0, 10, 20,32), labels = c("A-ge","B-chuu","C-jou"))
           yrStg = paste(Month, mthStg, sep = "-")
  })
( ggplot(aq, aes(yrStg, Temp))
  + geom_boxplot(outlier.size = .5)
  + ylab("Temperature")
  + xlab("Year Stage")
  + theme(axis.text.x = element_text(angle = 90, hjust = 1))
  )
```

### Dendrogram

See hierarchical clustering for base `plot` and `hclust`

```{r "dendrogram_ggplot"}

library(ggplot2)
library(ggdendro)
arrestsHc <- USArrests %>% dist %>% hclust("ave")
arrestsHc %>% ggdendrogram
arrestsHc %>% ggdendrogram(rotate = TRUE)

```

```{r "hc dendrogram with graphics::plot"}

arrestsHc %>% plot
arrestsHc %>% rect.hclust(5)
arrestsHc %>% plot
arrestsHc %>% rect.hclust(h = 40, which = c(2, 6), border = 3:4) # border is for colors

arrestsHc %>% plot
arrestsHc %>% rect.hclust(3)
arrestsHc %>% rect.hclust(h = 30, which = c(3, 7), border = 3:4)

```

extracting dendrogram plot data

```{r "dendrogram data extraction"}

library(ggplot2)
arrestsDendroDt <- arrestsHc %>% as.dendrogram %>% dendro_data(type = "rectangle")
ggplot(segment(arrestsDendroDt)) +
  geom_segment(aes(x = x, y = y, xend = xend, yend = yend)) +
  coord_flip() +
  scale_y_reverse(expand = c(.2, 0))
.Last.value + theme_dendro() # remove axises and background

ggplot(segment(arrestsHc %>% as.dendrogram %>% dendro_data(type = "triangle"))) +
  geom_segment(aes(x = x, y = y, xend = xend, yend = yend)) +
  coord_flip() +
  scale_y_reverse(expand = c(.2, 0))

```

regression tree diagrams

```{r ""}

library(tree)
data(cpus, package = "MASS")
model <- tree(log10(perf) ~ syct + mmin + mmax + cach + chmin + chmax
              ,data = cpus)
tree_data <- dendro_data(model)
ggplot(segment(tree_data)) +
  geom_segment(aes(x = x, y = y, xend = xend, yend = yend, size = n)
               ,colour = "blue", alpha = 0.5) +
  scale_size("n") +
  geom_text(data = label(tree_data)
            ,aes(x = x, y = y, label = label), vjust = -0.5, size = 3) +
  geom_text(data = leaf_label(tree_data)
            ,aes(x = x, y = y, label = label), vjust = 0.5, size = 2) +
  theme_dendro()

ggplot(segment(tree_data)) +
  geom_segment(aes(x = x, y = y, xend = xend, yend = yend) ,colour = "blue") +
  scale_size("n") +
  geom_text(data = label(tree_data)
            ,aes(x = x, y = y, label = label), vjust = -0.5, size = 3) +
  geom_text(data = leaf_label(tree_data)
            ,aes(x = x, y = y, label = label), vjust = 0.5, size = 2) +
  theme_dendro()

ggplot() +
  geom_segment(data = tree_data$segments
               ,aes(x = x, y = y, xend = xend, yend = yend)) +
  geom_text(data = tree_data$labels
            ,aes(x = x, y = y, label = label), size = 3, vjust = 0) +
  geom_text(data = tree_data$leaf_labels
            ,aes(x = x, y = y, label = label), size = 3, vjust = 1) +
  theme_dendro()

```

```{r "trees with rpart"}

library(rpart)
rpart_data <-
  rpart(Kyphosis ~ Age + Number + Start
        ,method = "class"
        ,data = kyphosis) %>% dendro_data

ggplot() +
  geom_segment(data = rpart_data$segments
               ,aes(x = x, y = y, xend = xend, yend = yend)
               ,colour = "blue") +
  geom_text(data = rpart_data$labels
            ,aes(x = x, y = y, label = label), size = 3, vjust = 0) +
  geom_text(data = rpart_data$leaf_labels
            ,aes(x = x, y = y, label = label), size = 3, vjust = 1) +
  theme_dendro()

co2_rpart_data <-
  rpart(Treatment ~ Type + conc + uptake
        ,method = "class"
        ,data = CO2) %>% dendro_data

ggplot() +
  geom_segment(data = co2_rpart_data$segments
               ,aes(x = x, y = y, xend = xend, yend = yend)) +
  geom_text(data = co2_rpart_data$labels
            ,aes(x = x, y = y, label = label), size = 3, vjust = 0) +
  geom_text(data = co2_rpart_data$leaf_labels
            ,aes(x = x, y = y, label = label), size = 3, vjust = 1) +
  theme_dendro()

```

### Facet grid

- use the formula to tell how and what you want to stratify the data on.
- Plot is split differently if `Species` is on the x-side
- some added smoothing
- use y and x side if you want to split on multiple categorical variables

```{r}
( ggplot(iris, aes(y = Sepal.Length,
                   x = Sepal.Width))
  + facet_grid(Species ~ .)
  + geom_smooth()
  )
```

### Themes

```{r ggplot_themes}
library(dplyr)
library(ggplot2)
diamondsTrunc <- diamonds[1:5000,]

( ggplot(data=diamondsTrunc, aes(carat,price ))
  + geom_point(aes(colour= color))
  + theme_bw()
  )

( ggplot(data=diamondsTrunc, aes(carat,price ))
  + geom_point(aes(colour= color))
  + theme_classic()
  )

( ggplot(data=diamondsTrunc, aes(carat,price ))
  + geom_point(aes(colour= color))
  + theme_bw()
  )

( ggplot(data=diamondsTrunc, aes(carat,price ))
  + geom_point(aes(colour= color))
  + theme_dark()
  )

( ggplot(data=diamondsTrunc, aes(carat,price ))
  + geom_point(aes(colour= color))
  + theme_gray()
  )

( ggplot(data=diamondsTrunc, aes(carat,price ))
  + geom_point(aes(colour= color))
  + theme_light()
  )

( ggplot(data=diamondsTrunc, aes(carat,price ))
  + geom_point(aes(colour= color))
  + theme_linedraw()
  )

( ggplot(data=diamondsTrunc, aes(carat,price ))
  + geom_point(aes(colour= color))
  + theme_minimal()
  )
```

### Matrices {.tabset .tabset-fade .tabset-pills}

#### Overview

useful to see pairwise correlation of variables.

#### graphics::pairs

pairs with added loess smoother in lower and correlation in upper

```{r}
panel.cor <- function(x, y, digits=2, prefix="", cex.cor, ...)
{
    usr <- par("usr"); on.exit(par(usr))
    par(usr = c(0, 1, 0, 1))
    r <- abs(cor(x, y))
    txt <- format(c(r, 0.123456789), digits=digits)[1]
    txt <- paste(prefix, txt, sep="")
    if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
    text(0.5, 0.5, txt, cex = cex.cor * r)
}

pairs(~Sepal.Length+Sepal.Width+Petal.Length+Petal.Width, data=iris,
      lower.panel=panel.smooth, upper.panel=panel.cor,
      pch=20, main="Iris Scatterplot Matrix")
```

#### lattice

Has options to condition the scatterplot matrix on a factor.

```{r}
library(lattice)
super.sym <- trellis.par.get("superpose.symbol")
splom(~iris[1:4],
      groups = Species,
      data = iris,
      panel = panel.superpose,
      key = list(title = "Three Varieties of Iris",
                 columns = 3,
                 points = list(pch = super.sym$pch[1:3],
                 col = super.sym$col[1:3]),
                 text = list(c("Setosa", "Versicolor", "Virginica"))))
```

#### car

The `car` package can condition the scatterplot matrix on a factor, and optionally include lowess and linear best fit lines, and boxplot, densities, or histograms in the principal diagonal, as well as rug plots in the margins of the cells.

```{r, eval=FALSE, echo=FALSE}

dependencies <- c("minqa", "nloptr", "Rcpp", "RcppEigen", "lme4", "pbkrtest", "SparseM", "MatrixModels", "quantreg", "car")
devtools::install_github(paste("cran", dependencies, sep = "/"))

```

```{r}
library(car)
scatterplot.matrix(~mpg+disp+drat+wt|cyl, data=mtcars,
  	main="Three Cylinder Options")
```

#### gclus

The `gclus` package provides options to rearrange the variables so that those with higher correlations are closer to the principal diagonal. It can also color code the cells to reflect the size of the correlations.

```{r, eval=FALSE, echo=FALSE}
devtools::install_github("cran/gclus")
```

```{r}

library(gclus)
dta <- mtcars[c(1,3,5,6)] # get data
dta.r <- abs(cor(dta)) # get correlations
dta.col <- dmat.color(dta.r) # get colors

```

reorder variables so those with highest correlation are closest to the diagonal

```{r}

dta.o <- order.single(dta.r)
cpairs(dta, dta.o, panel.colors=dta.col, gap=.5,
main="Variables Ordered and Colored by Correlation" )

```

#### gpairs

```{r, eval = FALSE}

dependencies <- c("barcode", "lmtest", "vcd", "gpairs")
devtools::install_github(paste("cran", dependencies, sep = "/"))

```

```{r}
library(gpairs)

if (allexamples) {
  y <- data.frame(A=c(rep("red", 100), rep("blue", 100)),
                  B=c(rnorm(100),round(rnorm(100,5,1),1)),
                  C=runif(200),
                  D=c(rep("big", 150), rep("small", 50)),
                  E=rnorm(200))
  gpairs(y)
}

if (allexamples) {
  data(iris)
  gpairs(iris)
  gpairs(iris, upper.pars = list(scatter = 'stats'),
         scatter.pars = list(pch = substr(as.character(iris$Species), 1, 1),
                             col = as.numeric(iris$Species)),
         stat.pars = list(verbose = FALSE))
  gpairs(iris, lower.pars = list(scatter = 'corrgram'),
         upper.pars = list(conditional = 'boxplot', scatter = 'loess'),
         scatter.pars = list(pch = 20))
}

if (allexamples) {
  data(Leaves)
  gpairs(Leaves[1:10], lower.pars = list(scatter = 'loess'))
  gpairs(Leaves[1:10], upper.pars = list(scatter = 'stats'),
         lower.pars = list(scatter = 'corrgram'),
         stat.pars = list(verbose = FALSE), gap = 0)
  corrgram(Leaves[,-33])
}
```

#### GGally

```{r, eval=FALSE, echo=FALSE}

dependencies <- c("prettyunits", "progress", "reshape", "GGally")
devtools::install_github(paste("cran", dependencies, sep = "/"))

```

```{r ggally_matrix_plot }
library(GGally)
library(dplyr)
if (allexamples){
  ds <- read.csv("for-gally.csv")
  ds$sex <- ifelse(ds$female==1, "female", "male") %>% as.factor
  ds$housing <- ifelse(ds$homeless==1, "homeless", "housed") %>% as.factor
  ggpairs(ds,
  				columns=c("i1", "cesd", "housing", "sex"), # list the factor variables first to get better boxplots
  				diag=list(continuous="density",
  									discrete="bar"),
  				axisLabels="show")
}

if (allexamples){
  data(tips, package = "reshape")
  ggpairs(tips[, 1:3])
  ggpairs(tips, 1:3, columnLabels = c("Total Bill", "Tip", "Sex"))
  ggpairs(tips, 1:3, upper = "blank")
  ggpairs(tips[, 1:3], axisLabels="internal") # Only Variable Labels on the diagonal (no axis labels)
  ggpairs(tips[, 1:3], axisLabels="none") # Only Variable Labels on the outside (no axis labels)

  ggpairs(
    tips[, c(1, 3, 4, 2)],
    upper = list(continuous = "density", combo = "box"),
    lower = list(continuous = "points", combo = "dot")
  )
}

if (allexamples){
  data(diamonds, package="ggplot2")
  diamonds.samp <- diamonds %>% sample_n(30)

  ggpairs(
    diamonds.samp,
    columns = "carat cut clarity depth" %>% wordArray,
    mapping = ggplot2::aes(color = cut),
    upper = list(continuous = wrap("density", alpha = 0.5),
                 combo = "box"),
    lower = list(continuous = wrap("points", alpha = 0.3),
                 combo = wrap("dot", alpha = 0.4)),
    title = "Diamonds"
  )
}

custom_car <- ggpairs(mtcars[, c("mpg", "wt", "cyl")], upper = "blank", title = "Custom Example")

```

ggplot example taken from example(geom_text)

```{r}
plot <- ggplot2::ggplot(mtcars, ggplot2::aes(x=wt, y=mpg, label=rownames(mtcars)))
plot <- plot +
  ggplot2::geom_text(ggplot2::aes(colour=factor(cyl)), size = 3) +
  ggplot2::scale_colour_discrete(l=40)
custom_car[1, 2] <- plot
personal_plot <- ggally_text(
  "ggpairs allows you\nto put in your\nown plot.\nLike that one.\n <---"
)
custom_car[1, 3] <- personal_plot
custom_car
```

#### ggplot2::plotmatrix

Depricated.

```{r, eval = FALSE}
library(ggplot2)
plotmatrix(with(iris, data.frame(Sepal.Length, Sepal.Width, Petal.Length, Petal.Width)))
```

### ggplot2 with multiple datasets and colors

```{r}
library(dplyr)
library(magrittr)
library(ggplot2)
dpassed <- rep(c("Yes", "No"), each=1000)
dspeed <- c(rnorm(1000, 2),
           rnorm(1000, 1))
lvlData <- data.frame(passed = dpassed, speed = dspeed)
samples <- factor(c("Passed", "Total"))

```

TODO: FIX ME

```{r , eval = FALSE}
(
  ggplot(lvlData %>% filter(passed == "Yes"),
         aes(x=speed, color="Sample", shape="Sample", linetype="Sample"))
  + stat_ecdf(data=lvlData,
              aes(x=speed, color="Total", shape="Total", linetype="Total"))
  + stat_ecdf(aes(color="Passed", shape="Passed", linetype="Passed"))
  + xlab("")
  + ylab("inclusion")
  + labs(title="Average days to completion")
  + scale_color_manual(breaks=samples, values=c("green", "red"))
  + scale_shape_manual(breaks=samples, values=c(16, 16))
  + scale_linetype_manual(breaks=samples, values=c(1, 1))
  + labs(color = "dataset")
  + theme(legend.title=element_blank())
)
```

### add linear regression line

```{r}
dlmDtExample <-
  structure(list(date =
                 structure(c(3407, 3438, 3468, 3499, 3530, 3560, 3591, 3621, 3652, 3683, 3712, 3743, 3773, 3804, 3834,
                             3865, 3896, 3926, 3957, 3987, 4018, 4049, 4077, 4108, 4138, 4169, 4199, 4230, 4261, 4291,
                             4322, 4352, 4383, 4414, 4442, 4473, 4503, 4534, 4564, 4595, 4626, 4656, 4687, 4717, 4748,
                             4779, 4807, 4838, 4868, 4899, 4929, 4960, 4991, 5021, 5052, 5082, 5113, 5144, 5173, 5204,
                             5234, 5265, 5295, 5326, 5357), class = "Date"),
                 gas = c(4.9, 4.5, 3.1, 3.1, 4.6, 4.8, 6.7, 8.6, 8.3, 7.2, 9.2, 5.5, 4.7, 4.4,
                         3.4, 2.8, 4, 5.1, 6.5, 9.2, 7.7, 7.7, 8.9, 5.7, 5, 4.5, 3.3,
                         2.8, 4, 5.6, 6.6, 10.3, 8.5, 7.9, 8.9, 5.4, 4.4, 4, 3, 3.1, 4.4,
                         5.5, 6.5, 10.1, 7.7, 9, 9, 6.5, 5.1, 4.3, 2.7, 2.8, 4.6, 5.5,
                         6.9, 9.5, 8.8, 8.7, 10.1, 6.1, 5, 4.5, 3.1, 2.9, 4.8)),
            .Names = c("date", "gas"),
            row.names = c(NA, -65L),
            class = "data.frame")

( ggplot(dlmDtExample, aes(x = date, y = gas))
  + geom_line()
  + ggtitle("UK gas consumption: amount in 106 tonnes coal equivalent")
  + geom_smooth(method = "lm", se = FALSE)
)

```

or a polynomial

```{r}
( ggplot(dlmDtExample, aes(x = date, y = gas))
  + geom_line()
  + ggtitle("UK gas consumption: amount in 106 tonnes coal equivalent")
  + geom_smooth(method = "lm", se = FALSE, formula = y ~ poly(x, 4))
)

```

# Statistics {.tabset .tabset-fade .tabset-pills}

## Overview

Section for statistics. Distributions, testing etc.

## quantiles

quantiles for different probabilities, use `stats:quantile`

```{r "quantile"}

randomNormals <- rnorm(1001)
randomNormals %>% quantile
randomNormals %>% quantile(.4)
randomNormals %>% quantile(probs = c(.1, .5, 1, 2, 5, 10, 50, NA)/ 100)

```

## Distributions {.tabset .tabset-fade .tabset-pills}

### Overview

Section for probability distributions, e.g. Normal/Gaussian, Binomial, Beta

### Beta distribution

- Typically used to model probabilities.
- Conjugate Prior to a Binomial distribution.
- Say that we are modelling a baseball batter
  - alpha parameter can be interpreted as number of hits
  - beta parameter can be interpreted as number of misses
- Ref: [wiki](https://en.wikipedia.org/wiki/Beta_distribution)

Let's look at a Beta distribution with alpha = 81 and beta = 219.

```{r}

alpha0 <- 81
beta0 <- 219

qbeta(.5, alpha0, beta0)
alpha0 / (alpha0 + beta0) # mean for this probability

```

standard beta distribution function. shape1/2 being alpha/beta.

```{r, eval = FALSE}
dbeta(x, shape1, shape2, ncp = 0, log = FALSE)
pbeta(q, shape1, shape2, ncp = 0, lower.tail = TRUE, log.p = FALSE)
qbeta(p, shape1, shape2, ncp = 0, lower.tail = TRUE, log.p = FALSE)
rbeta(n, shape1, shape2, ncp = 0)
```

- `dbeta` - returns density. high values for `x` close to the mean. E.g. usable for plotting the distribution.
- `qbeta` - returns quantile. E.g. 0.5 quantile is the mean.
- `pbeta` - returns probability. E.g. the probability of hitting 0.25 or less is 22% but ~88% to be under 0.33 (obviously this is pure theory for a probability in a one event trial)
- `rbeta` - returns random variables outcomes from a Beta(alpha0, beta0)

```{r}

n <- 21
x <- seq(0, 1, length = n)
data.frame(x
           ,d = dbeta(x, alpha0, beta0)
           ,quantile = qbeta(x, alpha0, beta0) # cumulative quantile of attaning probability x
           ,prob = pbeta(x, alpha0, beta0) # cumulative probability of attaning quantile x
           ,random = rbeta(n, alpha0, beta0) # random outcomes of n beta distributions with alpha0/beta0
           ) %>% round(3)

```

Bayesian setting. New data comes in.

- The batter has 100 more hits and out of 300 tries.
- update in bayesian way using bayes rule
  - Likelihood follows a Binomial distribution.
  - The posterior is another Beta with new alpha and beta.

```{r}

alpha1 <- alpha0 + 100 # add the hits
beta1 <- beta0 + 300 - 100 # add the misses

data.frame(x
           ,q0 = qbeta(x, alpha0, beta0)
           ,q1 = qbeta(x, alpha1, beta1)
           )

```

Plotting this distributions

```{r}

library(ggplot2)

sim <- data.frame(a = c(81, 81 + 100),
                  b = c(219, 219 + 200)) %>%
  dplyr::group_by(a, b) %>%
  dplyr::do(dplyr::data_frame(x = seq(0, 1, .001), y = dbeta(x, .$a, .$b))) %>%
  dplyr::mutate(Parameters = paste0("\u03B1 = ", a, ", \u03B2 = ", b)) %>%
  dplyr::ungroup() %>%
  dplyr::mutate(Parameters = factor(Parameters, levels = unique(Parameters)))

sim %>% dplyr::filter(a == 81) %>%
  ggplot(aes(x, y, color = Parameters)) + geom_line() +
  xlim(0, .5) + ylab("Density of beta")

sim %>%
  ggplot(aes(x, y, color = Parameters)) + geom_line() +
  xlim(0, .5) + ylab("Density of beta")

```


## Statistical Tests {.tabset .tabset-fade .tabset-pills}

### Overview

Section for test. Parametric and non-paratmetric.

### Mann-Whitney-Wilcoxon Test

In statistics, the Mann-Whitney U test (also called the Mann-Whitney-Wilcoxon (MWW), Wilcoxon rank-sum test, or Wilcoxon-Mann–Whitney test) is a nonparametric test of the null hypothesis that two samples come from the same population against an alternative hypothesis, especially that a particular population tends to have larger values than the other.[@wikiMannWhitneyUTest]

Unlike the t-test it does not require the assumption of normal distributions. It is nearly as efficient as the t-test on normal distributions.[@wikiMannWhitneyUTest]

Use the Wilcoxon signed-rank test when samples are related / paired.

Example code [@wilcoxTest]

Two-sample test.

```{r}
library(graphics)
x <- c(0.80, 0.83, 1.89, 1.04, 1.45, 1.38, 1.91, 1.64, 0.73, 1.46)
y <- c(1.15, 0.88, 0.90, 0.74, 1.21)
wilcox.test(x, y, alternative = "greater")
wilcox.test(x, y, alternative = "greater",
            exact = FALSE, correct = FALSE) # H&W large sample approximation

wilcox.test(rnorm(10), rnorm(10, 2), conf.int = TRUE)

```

Formula interface.

```{r}

boxplot(Ozone ~ Month, data = airquality)
wilcox.test(Ozone ~ Month, data = airquality,
            subset = Month %in% c(5, 8))

```

### Wilcoxon signed-rank test

The Wilcoxon signed-rank test is a non-parametric statistical hypothesis test used when comparing two related samples, matched samples, or repeated measurements on a single sample to assess whether their population mean ranks differ (i.e. it is a paired difference test). It can be used as an alternative to the paired Student's t-test, t-test for matched pairs, or the t-test for dependent samples when the population cannot be assumed to be normally distributed. [@lowry11]

Example code [@wilcoxTest], [@rTutorWilcoxonSignedRankTest]

```{r}
x <- c(1.83,  0.50,  1.62,  2.48, 1.68, 1.88, 1.55, 3.06, 1.30)
y <- c(0.878, 0.647, 0.598, 2.05, 1.06, 1.29, 1.06, 3.14, 1.29)

wilcox.test(x, y, paired = TRUE, alternative = "greater")
```

One-sample test version.

```{r }
wilcox.test(y - x, alternative = "less",
            exact = FALSE, correct = FALSE) # H&W large sample approximation

library(MASS)
wilcox.test(immer$Y1, immer$Y2, paired=TRUE)
```

### Contingency table test

look into the `stats` cheatsheet for usage of these tests.

example data

```{r "data for contingency tables test"}
library(magrittr)
partyGenderTable <-
  structure( c(762, 484, 327, 239, 468, 477)
            ,.Dim = 2:3
            ,.Dimnames = structure(list( gender = c("F", "M")
                                        ,party = c("Democrat", "Independent", "Republican"))
            ,.Names = c("gender", "party"))
            ,class = "table")

set.seed(12)
experimentDt <-
  data.frame( exposure = gl(2, 20, labels = c("Control", "Treatment"))
             ,outcome = gl(2, 20, labels = c("Sick", "Healhty")) %>% sample(replace = TRUE)
             )
```

#### Pearson chi-square test

can be applied to tables.

```{r}
(XsqTest <- partyGenderTable %>% chisq.test)
XsqTest$observed   # observed counts (same as M)
XsqTest$expected   # expected counts under the null
XsqTest$residuals  # Pearson residuals
XsqTest$stdres     # standardized residuals
```

or to variables directly.

```{r}
with(experimentDt, { chisq.test(exposure, outcome) })
```

#### Fisher's exact test

use function `fisher.test` and apply similar to `chisq.test`.

#### McNemar's test

for paired observations. use function `mcnemar.test`.

```{r "McNemars test"}
(Performance <-
  matrix(c(794, 86, 150, 570),
         nrow = 2,
         dimnames = list("1st Survey" = c("Approve", "Disapprove"),
                         "2nd Survey" = c("Approve", "Disapprove")))
 )

Performance %>% mcnemar.test
```

Here there's a significant difference in the approval from the 1st to the 2nd survey.

### ANOVA

```{r "ANOVA and Tukey HSD pairwise comparisons"}
library(magrittr)
warpbreaks %T>%
  { message("Data Summary"); print(summary(.)); message() } %>%
  aov(breaks ~ wool + tension, data = .) %T>%
  { message("ANOVA Summary"); print(summary(.)); message(); message("Pairwise comparisons") } %>%
  TukeyHSD("tension", ordered = TRUE)
```

# Supervised {.tabset .tabset-fade .tabset-pills}

### Notations {#supervised-notations}

- $p$ number of variables / parameters to fit
- $m$ number of variable to try each split (Random Forest)

### Model selection

Linear regression. BSS - Best subset selection

```{r leaps regsubsets}
library(magrittr)
data(Hitters, package="ISLR")

Hitters %>% summary
Hitters %<>% na.omit
regfit.full <- leaps::regsubsets(Salary ~ ., data = Hitters)
reg.summary <- regfit.full %>% summary
reg.summary %>% names

summariseRegsubset <- function(regsubset, icLabel){
  ic <- regsubset %>% summary %>% .[[icLabel]]
  plot(ic, xlab = "#features", ylab = icLabel)
  points(ic %>% which.min, ic %>% min, pch = 20, col = "red")
  coef(regsubset, ic %>% which.min)
}

"Mallow's CP" %>% message
reg.summary$cp
summariseRegsubset(regfit.full, "cp")
plot(regfit.full, scale = "Cp") # see ?leaps::plot.regsubsets for more scale options "r2", "adjr2" and "bic" are also available

"BIC" %>% message
reg.summary$bic

summariseRegsubset(regfit.full, "bic")
```

Forward / Backward selection

```{r "leabs forward / backward selection with leaps", eval = FALSE}

regfitFwd <- leaps::regsubsets(Salary ~ ., data = Hitters, method = "forward")
summariseRegsubset(regfitFwd, "bic")
plot(regfitFwd, scale = "bic")

regfitBwd <- leaps::regsubsets(Salary ~ ., data = Hitters, method = "backward")
summariseRegsubset(regfitBwd, "bic")
plot(regfitBwd, scale = "Cp")

```

GLM using `bestglm`

### Logistic Regression

```{r "install ROCR", eval = FALSE}

dependencies <- c("gplots", "rocr")
devtools::install_github(paste("cran", dependencies, sep = "/"))

```

In R, `glm` with `family` set is how logistic regression is implemented.

```{r "logistic regression"}
logRegMod <-
  glm(wool ~ .,
      data = warpbreaks,
      family = binomial(link = "logit"))

logRegMod %>% summary

```

calculate model BIC

```{r "bic calc"}

BIC <- function(model, nSamples) model %>% AIC(k = log(nSamples))
BIC(logRegMod, warpbreaks %>% nrow)

```

ROC

```{r , "log reg ROC"}

logRegMod %>% fitted %>%
  ROCR::prediction(warpbreaks$wool) %>%
  ROCR::performance(measure = "tpr",
                    x.measure = "fpr") %>%
  (ROCR::plot)

```

### GLM gotcha

with variables that has too many labels compared to data observations (i.e. we have no DF)

```
>     glm(poisson2 ~.
+         ,data = strata %>% toStrataTrainingDta(outcomeVar = "poisson2")
+         ,family = "poisson")
Error: inner loop 1; cannot correct step size
In addition: Warning message:
step size truncated due to divergence
```

### Negative Binomial

with GLM

```{r, "negative binomial example"}
data(quine, package = "MASS")
quine %>% str
quine.nb1 <- MASS::glm.nb(Days ~ Sex/(Age + Eth*Lrn), data = quine)
quine.nb2 <- update(quine.nb1, . ~ . + Sex:Age:Lrn)
quine.nb3 <- update(quine.nb2, Days ~ .^4)
anova(quine.nb1, quine.nb2, quine.nb3)
quine.nb1 %>% summary
```

### SVM

```{r "svm basic example"}

model <- e1071::svm(Species ~ ., data = iris, probability = TRUE)
summary(model)

```

model fitted vs real value (no CV, cheating of course)

```{r}

fitted(model) %>% table(iris$Species)

```

obtain label probabilities for each observation

```{r}
newObservations <- subset(iris %>% head, select = -Species)
predict(model, newObservations, decision.values = TRUE, probability = TRUE)
```

### recursive partitioning

#### CART

clustering with CART. see `ggdendro` for plotting with `ggplot2`

```{r "rpart"}
library(rpart)
rpart(Kyphosis ~ Age + Number + Start
      ,method = "class"
      ,data = kyphosis) %T>% plot %>% text(use.n = TRUE)
```

#### Conditional inference trees

Classification in Kyphosis dataset like we did with `rpart`

```{r "party::ctree on Kyphosis data"}
library(rpart)
library(partykit) # library(party)
ctree(Kyphosis ~ Age + Number + Start
      ,data = kyphosis) %T>% print %>% plot(main = "Conditional Inference of Kyphosis data")
```

regression

```{r "party::ctree regression"}
set.seed(290875)

ozoneDt <- airquality %>% dplyr::filter(!is.na(Ozone))
airct <-
  ctree( Ozone ~ .
        ,data = ozoneDt)
airct %T>% print %>% plot
```

you can grab the data in the different segments clusters using `party::where`. using this you can extract cluster data or aggregates.

```{r "party::ctree segments data/aggregates"}
library(dplyr)
ozoneDt %>%
  mutate(segment = predict(airct)) %>%
  group_by(segment) %>%
  summarise(avgOzone = mean(Ozone, na.rm = TRUE))

```

Mean squared error

```{r, eval = FALSE}
mean((airq$Ozone - predict(airct))^2)
```

extract observation terminal node ID using either `predict` or `where`

```{r "party::ctree regression data extraction"}
airct %>% predict(type = "node") # airct %>% party::where # only in party
```

```{r "party::ctree BinaryTree handling", eval = FALSE}
airct %>% nodes(4) # extract nodes # party only
airct %>% response %>% summary # extract the response variable from the formula. # party only
```

nominal variable classification

```{r "party::ctree classification"}
irisct <- ctree(Species ~ .,data = iris)
irisct %T>% print %>% plot
table(predict(irisct), iris$Species)
```

response of tree for the given observations (probabilities of levels in Species)

```{r, eval = FALSE}
treeresponse(irisct, newdata = iris[c(1,61, 120),]) # party only
```

ordinal regression. `ME` from the `mammoexp` data frame is an ordinal factor.

```{r "party::ctree ordinal regression"}
data("mammoexp", package = "TH.data")
mammoct <- ctree(ME ~ ., data = mammoexp)
plot(mammoct)
```

response of tree for the given observations (probabilities of levels in ME)

```{r, eval = FALSE}
treeresponse(mammoct, newdata = mammoexp[1:4,]) # party only
```

survival analysis.

```{r "party::ctree survival"}
library("survival")
data("GBSG2", package = "TH.data")
GBSG2 %>% str
GBSG2ct <- ctree(Surv(time, cens) ~ .,data = GBSG2)
plot(GBSG2ct)

```

the tree response here is a whole survival fit (survfit object, e.g. can be plotted)

```{r, eval = FALSE}
treeresponse(GBSG2ct, newdata = GBSG2[1:2,]) # party only
```

options (see `?ctree_control` for options)

- `mincriterion` for significance level
- `testtype` for multiple testing method
- `minbucket` for smallest possiple bucket
- `maxsurrogate` unsure!

```{r}
data(plantTraits, package = "cluster")
ctree( height ~ .
      ,data = plantTraits
      ,mincriterion = .9
      # ,testtype = "Teststatistic"
      ) %>% plot

ctree( height ~ .
      ,data = plantTraits
      ,testtype = "Bonferroni"
      ) %>% plot

ctree( height ~ .
      ,data = plantTraits
      ,mincriterion = .8
      ,testtype = "Bonferroni"
      ) %>% plot

ctree( height ~ .
      ,data = plantTraits
      ,minbucket = 10
      # ,testtype = "Teststatistic"
      ) %>% plot

ctree( height ~ .
      ,data = plantTraits
      ,maxdepth = 3 #,testtype = "Teststatistic"
      ) %>% plot

ctree( height ~ .
      ,data = plantTraits
      ,mincriterion = .9 #,testtype = "Teststatistic"
      ) %>% plot

ctree( height ~ .
      ,data = plantTraits
      ,mincriterion = .9 ,maxsurrogate = 2 # ,testtype = "teststatistic"
      ) %>% plot
```

extracting node splits conditions for each terminal node.

[SO link](https://stackoverflow.com/questions/21443203/ctree-how-to-get-the-list-of-splitting-conditions-for-each-terminal-node)

```{r "install emilrehnberg/party.readpaths", eval = FALSE}
devtools::install_github("emilrehnberg/party.readpaths")
```

```{r "generated data and ctree"}
shiftFirstPart <- function(vctr, divideBy, proportion = .5){
  vctr[vctr %>% length %>% multiply_by(proportion) %>% round %>% seq] %<>% divide_by(divideBy)
  vctr
}
set.seed(11)
n <- 13000
gdt <-
  data.frame( is_buyer = runif(n) %>% shiftFirstPart(1.5) %>% round %>% factor(labels = c("no", "yes"))
             ,age = runif(n) %>% shiftFirstPart(1.5) %>%
               cut(breaks = c(0, .3, .6, 1), include_lowest = TRUE, ordered_result = TRUE, labels = c("low", "mid", "high"))
             ,city = runif(n) %>% shiftFirstPart(1.5) %>%
               cut(breaks = seq(0,1,.2), include_lowest = TRUE, labels = c("Chigaco", "Boston", "Memphis", "LA", "ATL"))
             ,point = runif(n) %>% shiftFirstPart(1.2)
             )
gct <- party::ctree( is_buyer ~ ., data = gdt)

party.readpaths::readCtreePaths(gct, gdt) # dts <- gdt; ct <- gct
gct %>% plot

```

In `partykit`

```{r, eval = FALSE}
partykit:::.list.rules.party(ct)
```

```{r}
detach(package:partykit)
airq <- subset(airquality, !is.na(Ozone))
act <-
  party::ctree( Ozone ~ .
        ,data = airq # ,maxsurrogate = 3 # partykit
        ,controls = party::ctree_control(maxsurrogate = 3) # party
        )
party.readpaths::readCtreePaths(act, airq) # party only

data(plantTraits, package = "cluster")
pct <- party::ctree( height ~ .  ,data = plantTraits)
party.readpaths::readCtreePaths(pct, plantTraits)
```

prediction for the `party` package.

terminal nodes & response data prediction

```{r "prediction on party::ctree objects"}
pct %>% party::where(newdata = plantTraits[1:10,])
pct %>% party::treeresponse(newdata = plantTraits[1:10,]) # not sure what this is.. >.<
```

response prediction

```{r}
pct %>% predict(newdata = plantTraits[1:10,])
```


### PLS - Partial Least Squares
```{r}
```

### Random Forests {.tabset .tabset-fade .tabset-pills}

#### Overall notes

- After each tree is built, all of the data are run down the tree, and proximities are computed for each pair of cases. [@rfWorkings]
    - If two cases occupy the same terminal node, their proximity is increased by one.
    - At the end of the run, the proximities are normalized by dividing by the number of trees.
    - Proximities are used in replacing missing data, locating outliers, and producing illuminating low-dimensional views of the data.

- to compare random forests, perform wilcoxon test on mean differences of absolute errors. [@esl2]

In case of level errors like:

```
Error in randomForest.default(m, y, ...) :
  New factor levels in xtest not present in x
```

1. It can be due to more labels being available in test data.
2. The Function also do not like ordered factors. Original factors only!
3. To delve deeper into the actual levels in the `randomForest`-object use `rfObj$forest$xlevels`

If you are using `predict`, you might stuble upon the following error if you forgot the `keep.forest = TRUE` argument.:

```
Error in predict.randomForest(model, newdata = CV) :  No forest component in the object
```

#### Classification

- Categorical (e.g. Binary) Response
- authors suggest parameters should be (See [Notations](#supervised-notations))[@esl2rfDetails]
    - $\min(m) = 1$
    - $m = \sqrt(p)$

Example code [@rfPackage] [@sll10a]

```{r, eval=FALSE, echo=FALSE}
devtools::install_github('cran/randomForest')
```

```{r}
library(randomForest)
set.seed(71)
(iris.rf <-
	randomForest(
		Species ~ .,
		data=iris,
		importance=TRUE,
    do.trace = 100, # adds the extra output per 100 tree iteration
		proximity=TRUE))

```

Look at variable importance:

```{r}
round(importance(iris.rf), 2)
```

- proximity: refers to a proximity matrix which is the distance $\in[0,1]$ between the observations.
    - A matrix of proximity measures among the input (based on the frequency that pairs of data points are in the same terminal nodes).
- By performing MDS on the proximity we can get a feeling for how the observations group between each-other.
- $1 - prox(n, k)$ are squared distances, hence we can do MDS on those distances to see groupings [@rfScaling]

```{r}
## Do MDS on 1 - proximity:
iris.mds <- cmdscale(1 - iris.rf$proximity, eig=TRUE) # Classical (Metric) Multidimensional Scaling
op <- par(pty="s")
pairs(cbind(iris[,1:4], # Scatterplot Matrices
            iris.mds$points),
      cex=0.6,
      gap=0,
      col=c("red", "green", "blue")[as.numeric(iris$Species)],
      main="Iris Data: Predictors and MDS of Proximity Based on RandomForest")
par(op)
print(iris.mds$GOF)
```

include test data in the `randomForest` call to get the full proximity matrix between test and train data [@rfProximityForAll]

```{r}
set.seed(71)
ind <- sample(1:150,140,replace = FALSE)
train <- iris[ind,]
test <- iris[-ind,]

(iris.rf1 <- randomForest(x = train[,1:4],
                         y = train[,5],
                         xtest = test[,1:4],
                         ytest = test[,5],
                         importance=TRUE,
                         proximity=TRUE))
dim(iris.rf1$test$prox)
```

"x" can be a matrix instead of a data frame:

```{r}
set.seed(17)
x <- matrix(runif(5e2), 100)
y <- gl(2, 50) # Generate Factor Levels: 2 levels, 50 of each level
(myrf <- randomForest(x, y)) # parentheses leads to printing of the declaration
predict(myrf, x)

```

stratified sampling: draw 20, 30, and 20 of the species to grow each tree.

```{r}

randomForest(Species ~ .,
             data = iris,
             sampsize=c(20, 30, 20))
```

#### Regression

- continuous response
- authors suggest parameters should be (See [Notations](#supervised-notations))[@esl2rfDetails]
    - $\min(m) = 5$
    - $m = \sqrt(p)$

Example code [@rfPackage] [@sll10b]

```{r message=FALSE}

library(MASS)
library(randomForest)
library(magrittr)
library(ggplot2)

```

```{r}

set.seed(101)
Boston %>% str # ?Boston # to see more data info. it's housing values in Boston suburbs.

train <- 1:nrow(Boston) %>% sample(300)

randomForest(medv ~ .,
             data = Boston,
             subset = train)

nTreesToFit <- 400
maxVarsToTry <- 4 # 13
forests <-
  lapply(1:maxVarsToTry, function(varsToTry){
    randomForest(medv ~ .,
                 data = Boston,
                 subset = train,
                 mtry = varsToTry, # number of variables tried at each split
                 importance = TRUE,
                 ntree = nTreesToFit)
  })

```

Variables of highest "Importance"

```{r}

forests[[maxVarsToTry]] %>% importance %>% round(2) %>% as.data.frame %>% arrange(desc(`%IncMSE`)) # Removes row names >.<

testData <- Boston[-train,]
errors <-
  data.frame(varsPerSplit = 1:maxVarsToTry,
             oob = forests %>% sapply(function(forest){ forest$mse[nTreesToFit] }),
             test = forests %>%
               sapply(function(forest){
                 pred <- predict(forest, testData)
                 with(testData,
                      mean((medv - pred)^2))
               })) %>%
  tidyr::gather( src
                ,errors
                ,-c(varsPerSplit))

ggplot(errors, aes(x = varsPerSplit %>% as.factor, y = errors, group = src, color = src)) +
  geom_line() +
  geom_point() +
  xlab("Variables per split") +
  ylab("Mean Squared Error")

set.seed(131)
(ozone.rf <- randomForest(Ozone ~ ., data=airquality, mtry=3,
                         importance=TRUE, na.action=na.omit))

```

Show "importance" of variables: higher value mean more important:

```{r}

ozone.rf %>% importance %>% round(2)

```

randomForest call with test data

```{r}

set.seed(131)
trainObservations <- seq(1:nrow(airquality)) %>% sample(133)
testData <- airquality[-trainObservations, ] %>% na.omit # no NAs in the test data!
randomForest(Ozone ~ .,
             data=airquality,
             subset=trainObservations,
             xtest=testData %>% dplyr::select(-Ozone),
             ytest=testData %>% use_series("Ozone"), # atomic vector needed
             mtry=3,
             importance=TRUE,
             na.action=na.omit)

```

"complicated" formula:

```{r}

(swiss.rf <- randomForest(sqrt(Fertility) ~ . - Catholic + I(Catholic < 50),
                          data=swiss))
predict(swiss.rf, swiss)

```

Test use of 53-level factor as a predictor:

```{r}

set.seed(1)
x <- data.frame(x1=gl(53, 10),
                x2=runif(530),
                y=rnorm(530))
(rf1 <- randomForest(y ~ ., data=x, ntree=10))

```

Grow no more than 4 nodes per tree:

```{r}

randomForest(Species ~ ., data=iris, maxnodes=4, ntree=30) %>% treesize

```

test proximity in regression

```{r}

(iris.rrf <-
  randomForest(Sepal.Width ~ .,
               data=iris,
               ntree=101,
               proximity=TRUE,
               oob.prox=FALSE))
iris.rrf$proximity %>% str

```

#### Imputation

Random forests has two ways of replacing missing values.

1. fast; less performant. If the $m$-th variable is
    - continuous: compute the median of all values of this variable in class $j$. Use value to replace all missing values of the $m$-th variable in class $j$.
    - categorical: replace the most frequent non-missing value in class $j$. These replacement values are called fills.
2. computationally more expensive but has more performant (even with large amounts of missing data). Replaces missing values only in the training set.
    1. begins with a rough and inaccurate filling in of the missing values.
    2. Perform forest run and computes proximities.
    3. If $x(m,n)$ is a
        - continuous: estimate its fill as an average over the non-missing values of the $m$-th variables weighted by the proximities between the $n$-th case and the non-missing value case.
        - categorical: replace it by the most frequent non-missing value where frequency is weighted by proximity.
    4. Now iterate-construct a forest again using these newly filled in values, find new fills and iterate again. Our experience is that 4-6 iterations are enough.


`randomForest::na.roughfix` implements imputation #1.

```{r}
iris.na <- iris

set.seed(111)
for (i in 1:4) iris.na[sample(150, sample(20)), i] <- NA ## artificially drop some data values.

iris.na %>% summary
na.roughfix(iris.na) %>% summary

randomForest(Species ~ ., iris.na, na.action = na.omit)
randomForest(Species ~ ., iris.na, na.action = na.roughfix)
```

`randomForest::rfImpute` implements imputation #2.

```{r}
iris.na <- iris

set.seed(111)
for (i in 1:4) iris.na[sample(150, sample(20)), i] <- NA ## artificially drop some data values.

set.seed(222)
iris.imputed <- rfImpute(Species ~ ., iris.na)

set.seed(333)
randomForest(Species ~ ., iris.na, na.action = na.omit)
randomForest(Species ~ ., iris.imputed)
```

### Boosting { .tabset .tabset-fade .tabset-pills }

#### Regression

random forest boosting with continous response.

- authors suggest $m = \sqrt(p)$ (See [Notations](#supervised-notations))
- see latter half of `mpv Movies/statistical-learning/08-tree-based-methods/StatsLearning Lect10 R trees B 111213-IY7oWGXb77o.mp4`

```{r, eval=FALSE, echo=FALSE}
devtools::install_github("cran/gbm")
```

```{r message=FALSE}

library(MASS)
library(magrittr)
library(gbm) # Gradient Boosted Machines

```

```{r}

set.seed(101)
Boston %>% str # ?Boston # to see more data info. it's housing values in Boston suburbs.

train <- 1:nrow(Boston) %>% sample(300)

boost.boston <-
  gbm(medv ~ .,
      data = Boston[train,],
      distribution = "gaussian",
      n.trees = 10000,
      shrinkage = 0.01,
      interaction.depth = 4)
boost.boston %>% summary # gives variance importance plot.

```

variables lstat (lower status) and rm (rooms) are impacting the most.

```{r}

boost.boston %>% plot(i = "lstat") # partical dependency plot for the two most important plots.
boost.boston %>% plot(i = "rm")

n.trees <-
  seq(from = 100,
      to = 10000,
      by = 100)
testData <- Boston[-train, ]
predmat <-
  predict(boost.boston,
          newdata = testData,
          n.trees = n.trees)
predmat %>% dim
berr <-
  with(testData,
       apply((predmat - medv)^2, 2, mean))
plot(n.trees,
     berr,
     pch = 19,
     ylab = "Mean Squared Error",
     xlab = "# Trees",
     main = "Boosting Test Error")
```

From the RF w/o boosting.
Run that example first if you want to compare the boost vs non-boost results.

```{r "boosting plot: doesn't work", eval = FALSE}

abline(h = min(errors$test),
       col="red")
```

#### Classification

```{r message=FALSE}
library(magrittr)
library(gbm) # Gradient Boosted Machines
```

```{r}
set.seed(101)
trainObservations <- sample(1:nrow(iris), 75)
testData <- iris[-trainObservations, ]

( fit <-
	gbm(Species~.,
			data=iris[trainObservations, ],
			distribution="multinomial"))
nTrees <- 4^(1:3)

predictions <-
  predict(fit,
          newdata = testData,
          n.trees = nTrees)

( classifiationTables <-
  sapply(nTrees %>% as.character,
         function(treeSize){
           ( predictions[,,treeSize]
             %>% apply(1, which.max)
             %>% factor(labels = predictions %>% dimnames %>% .[[2]])
             %>% table(testData$Species)
            )
         }, simplify = FALSE))

```

prediction accuracy

```{r}

( classifiationTables
  %>% sapply(diagonalProportion)
 )

```

only thing required for one n.trees

```{r, eval =FALSE}

(
  predictions
  %>% apply(1, which.max)
  %>% factor(labels = predictions %>% dimnames %>% .[[2]])
  %>% table(testData$Species)
 )

```

### HMM

Hidden markov models. With `depmixS4` package

```{r eval = FALSE}

devtools::install_github(paste("cran", c("truncnorm", "Rsolnp", "depmixS4"), sep="/"))

```

data `speed` consists of three time series with three variables:

1. response time `rt`
2. accuracy `corr` (and previous `corr` in `prev`)
3. covariate, `Pacc`, the relative pay-off for speeded versus accurate responding.

```{r "simple depmix model"}

data("speed", package = "depmixS4")
speed %>% str
speed$Pacc %>% table %>% as.matrix

speedSimpleModel <- list(mixModel =
  depmixS4::depmix(
    response = rt ~ 1
    ,data = speed
    ,nstates = 2 # number of model states
    ,trstart = runif(4) # transition parameters
   )
 )

speedSimpleModel$fitted <- depmixS4::fit(speedSimpleModel$mixModel, emc = depmixS4::em.control(rand = FALSE))

speedSimpleModel %>% print

```

Covariates on the transition probabilities

```{r, "depmix with covariates in transitions"}

speedTransCovarModel <- list(unfitted =
  depmixS4::depmix(
    response = rt ~ 1
    ,data = speed
    ,nstates = 2
    ,family = gaussian()
    ,transition = ~ scale(Pacc)
    ,instart = runif(2)
   )
 )
speedTransCovarModel$fitted <- depmixS4::fit(speedTransCovarModel$unfitted, emc = depmixS4::em.control(rand = FALSE) )

```


### Classification

Other classification methods are e.g.

- C4.5 (weka package?) can be used for classification with more labels
    - categorical inputs splits in multiple nodes, one for each level value


# Unsupervised { .tabset .tabset-fade .tabset-pills }

### Overview

methods for unsupervised learning.

### PCA/PCR

principal component analysis

```{r "PCA USArrests"}

prcomp(USArrests %>% dplyr::select(-state))  # inappropriate, due to lack of scaling

```

prints standard deviations & loadings

```{r}

(usarrestsPcaFit <- prcomp(USArrests %>% dplyr::select(-state), scale = TRUE))
usarrestsPcaFit %>% summary   # prints Importance of components
usarrestsPcaFit %>% loadings  # prints loadings?
usarrestsPcaFit$scores        # prints principal components
usarrestsPcaFit %>% biplot    # plots a PC1/PC2 with loadings

```

using formula

```{r}

prcomp(~ Murder + Assault + Rape, data = USArrests, scale = TRUE)
prcomp(~ . - state, data = USArrests, scale = TRUE)

```

only works with numerical variables

```{r, eval = FALSE}
prcomp(warpbreaks %>% dplyr::select(-breaks), scale = TRUE) # doesn't work
```

PCA components with color from factor variable.

```{r "PCA iris"}
prcompIris <- prcomp(iris %>% dplyr::select(-Species), scale = TRUE)
prcompIris$x %>% as.data.frame %>%
  ggplot2::ggplot(ggplot2::aes(x = PC1, y = PC2, color = iris$Species)) + ggplot2::geom_point()
```

to fortify these simple PCA plots look up [fortify PCA ploting](https://cran.r-project.org/web/packages/ggfortify/vignettes/plot_pca.html)

### Factor Analysis

Exploratory Factor Analysis

[ref](https://www.statmethods.net/advstats/factor.html)

```{r, "efa"}

harmanFA <-
  factanal(covmat = Harman74.cor
           ,factors = 3
           ,rotation = "varimax"
           )

harmanFA %>% print(digits = 2, cutoff = .3, sort = TRUE)

firstLoadings <- harmanFA$loadings[,1:2]
firstLoadings %>% plot(type = "n")
firstLoadings %>% text(labels = Harman74.cor$cov %>% row.names, cex = .7)

```


### hierarchical clustering

calculate a distance matrix and run it through `hclust`. The `hclust` object can be plotted through `graphics::plot`

```{r "hierarchical clustering"}
USArrests %>% dist %>% hclust("average") %>% plot
USArrests %>% dist %>% hclust("complete") %>% plot(hang = -1)
USArrests %>% dist %>% hclust(method = "average")
USArrests %>% dist %>% hclust(method = "complete")
```

```{r "hclust plot w labels"}
npk %>% dplyr::select(-block) %>% dist %>% hclust("complete") %>% plot(hang = -1)
npk %>% dplyr::select(-block) %>% dist %>% hclust("complete") %>% plot(hang = -1, labels = npk$block)
```

prune tree, compare 2 and 4 group pruning

```{r "simple prune tree"}
USArrests %>% dist %>% hclust %>%
  cutree(k = c(2,4)) %>% as.data.frame %>%
  with({ table(get("2"), get("4")) })
```

1. Do centroid clustering and *squared* Euclidean distance,
2. cut the tree into ten clusters
3. reconstruct the upper part of the tree from the cluster centers.

```{r "hclust cut tree example"}
hc <- hclust(dist(USArrests)^2, "cen")
memb <- cutree(hc, k = 10)
cent <-
  sapply(1:10, function(groupIndex){
    colMeans(USArrests[memb == groupIndex, 1:4, drop = FALSE])
   }, simplify = FALSE) %>%
    Reduce(f = rbind) %>% as.matrix

hc1 <- dist(cent)^2 %>% hclust(method = "cen", members = table(memb))
opar <- par(mfrow = c(1, 2))
plot(hc,  labels = FALSE, hang = -1, main = "Original Tree")
plot(hc1, labels = FALSE, hang = -1, main = "Re-start from 10 clusters")
par(opar)
```

### Random Forest

Example code [@rfPackage].

The `unsupervised' case:

```{r "MDS with RF"}
library(randomForest)
set.seed(17)
iris.urf <- randomForest(iris %>% dplyr::select(-Species))
MDSplot(iris.urf, iris$Species) # Species sets the colors.
```

Works with a combination of numeric and factor variables.

```{r "RF MDS w factors & numericals"}
co2.rf <- randomForest(CO2)
MDSplot(co2.rf, CO2$Type)

co2.rf <- randomForest(CO2 %>% dplyr::select(-Treatment, -Type))
MDSplot(co2.rf, CO2$Type)
```

only factors also works

```{r "RF MDS w factors only"}
warp.rf <- randomForest(warpbreaks %>% dplyr::select(-breaks))
MDSplot(warp.rf, 1:2)
```

### Dissimilarity Matrix Calculation

`cluster::daisy`

Calculates pairwise dissimilarities (distances) between observations. Variables may be of mixed types.

```{r "daisy setup"}
library(cluster)
library(magrittr)
data(agriculture)
agriculture %>% str
```

daisy with euclidean distance (non-standardized variables i.e. default)

```{r "daisy euclidean distance"}
( d.agr <-
  agriculture
  %>% daisy(metric = "euclidean", stand = FALSE)
  %T>% print
  %>% as.matrix %>% .[, "DK"]
 )
```

compare with gower metric

```{r "daisy gower metric"}
agriculture %>% daisy(metric = "gower") # %>% hclust %>% plot
```

```{r "daisy flower example setup"}
data(flower)
flower %>% str
flower %>%
  daisy(type = list(asymm = "V3")) %>%
  summary
flower %>%
  daisy(type = list( asymm = c("V1", "V3")
                    ,ordratio = "V7")) %>%
  summary

# hierarchical clustering
flower %>% daisy(type = list(asymm = 3)) %>% hclust %>% plot
```

### Agglomerative nesting (hierarchical clustering)

`cluster::agnes`

Computes agglomerative hierarchical clustering of the dataset.

```{r "agnes examples setup"}
data(votes.repub)
```

```{r "agnes example 1"}
par(mfrow = c(2, 1))
( votes.repub
  %>% agnes(metric = "manhattan", stand = TRUE)
  %T>% print
  %>% plot
 )
```

### Forecasting { .tabset .tabset-fade .tabset-pills }

Section for forecasting

```{r, eval=FALSE, echo=FALSE}

dependencies <- c("timeDate", "quadprog", "TTR", "quantmod", "tseries", "fracdiff", "RcppArmadillo", "forecast")
devtools::install_github(paste("cran", dependencies, sep = "/"))

```

```{r}
library(magrittr)
library(mgcv)
library(forecast)
x <-
  rpois(100, 1 + sin(seq(0, 3*pi, l=100))) %>%
  ts(f = 12)
tt <- 1:100
(season <- x %>% seasonaldummy) %>% head
# x <- ts(rpois(100,1+sin(seq(0,3*pi,l=100))),f=12)
fit <- gam(x ~ s(tt, k=5) + season,
           family="poisson")
fit %>% plot
fit %>% summary

fcast <-
  predict(fit,
          se.fit=TRUE,
          newdata=list(tt=101:112,
                       season=seasonaldummy(x,h=12)))

x %>% plot(xlim = c(0, 10.5))
fcast$fit %>%
  exp %>%
  ts(f = 12, s = 112/12) %>%
  lines(col = 2)
fcast %>%
  { .$fit - 2*.$se } %>%
  exp %>%
  ts(f = 12, s = 112/12) %>%
  lines(col = 2, lty = 2)
fcast %>%
  { .$fit + 2*.$se } %>%
  exp %>%
  ts(f = 12, s = 112/12) %>%
  lines(col = 2, lty = 2)
```

# Neural networks { .tabset .tabset-fade .tabset-pills }

### Overview

Section for neural nets. Doesn't fit in neither Supervised nor Unsupervised sections since NNs can be both. AKA Artificial Neural Networks.

Packages

- `neuralnet`
- `nnet`
- `DARCH`
- `deepnet` [cran](https://cran.r-project.org/package=deepnet)
- `h2o`

### Theory

#### Basics

NN is a network of perceptrons. I.e. functions that take inputs & bias into an "activation function" and that yields an output.

The perceptron includes:

1. linear function
2. Activation function, transforming the output from the linear function to the desired outcome range.

The outputs are compared to known labels and coefficients/weights are adjusted (via gradient decent) repeatedly until either 1. max number of allowed iterations or 2. an acceptable error rate is hits.

The process of feeding inputs to multiple perceptrons yielding an output is called *feed-forward*. I.e. the units do not form a cycle!

#### Deep networks

aka DNN - deep neural networks

Vladmir Perervenko articles

- [THIRD GENERATION NEURAL NETWORKS: DEEP NETWORKS](https://www.mql5.com/en/articles/1103)
- [EVALUATION AND SELECTION OF VARIABLES FOR MACHINE LEARNING MODELS](https://www.mql5.com/en/articles/2029)
- [DEEP NEURAL NETWORK WITH STACKED RBM. SELF-TRAINING, SELF-CONTROL](https://www.mql5.com/en/articles/1628)

### h2o

- is Java setup requiring JVM
- the R installation is for an h2o API
- interface is on the web

- [ui](http://localhost:54321) - after running `h2o::h2o.init`
- [vignette code examples](https://github.com/h2oai/h2o-3/tree/master/h2o-docs/src/booklets/v2_2015/source/R_Vignette_code_examples)
- [cran](https://CRAN.R-project.org/package=h2o)
- [github](https://github.com/h2oai/h2o-3)
- [DNN use case, handwritten digits](http://didericksen.github.io/deeplearning-r-h2o/)

supports

- Deep Learning
- Naive Bayes
- Principal Components Analysis (PCA)
- K-means
- Stacked Ensembles
- Generalized Linear Models (GLM)
- Gradient Boosting Machine (GBM)
- Generalized Low Rank Model (GLRM)
- Distributed Random Forest (DRF)
- Word2ve

#### Installation

Check version @ [download HP](https://www.h2o.ai/download/)

```{r "install latest h2o", eval = FALSE}

"cran" %>% paste(c("RCurl", "jsonlite"), sep = "/") %>% (devtools::install_github)

# look at the \\d{4} in the repos, exchange when there's a new version. See new versions at the link above.
install.packages("h2o", type="source", repos="http://h2o-release.s3.amazonaws.com/h2o/master/4208/R")

```

#### getting started

```{r "h2o get started", eval = FALSE}

h2o::h2o.init(nthreads = 2, max_mem_size = '4g') # use -1 to use all

```

#### importing data

to work with data on h2o you need to import/upload data. Import is to a single node e.g. locally, but if you have a cluster you'll have to upload (`h2o.uploadFile()`).

```{r "h2o importing data", eval = FALSE}

#To import small iris data file from H2O’s package:
iris.hex <- h2o.importFile(
  path = system.file("extdata", "iris.csv", package="h2o")
  ,destination_frame = "iris.hex"
 )
iris.hex %>% (h2o::h2o.str)
iris.hex %>% as.data.frame %>% str
iris %>% str

#To import an entire folder of files as one data object:
# pathToFolder = "/Users/data/airlines/"
# airlines.hex = h2o.importFile(path = pathToFolder, destination_frame = "airlines.hex")

#To import from HDFS and connect to H2O in R using the IP and port of an H2O instance running on your Hadoop cluster:
 # h2o.init(ip= <IPAddress>, port =54321, nthreads = -1)
 # pathToData = "hdfs://mr-0xd6.h2oai.loc/datasets/airlines_all.csv"
 # airlines.hex = h2o.importFile(path = pathToData, destination_frame = "airlines.hex")

# Converts R object "iris" into H2O object "iris.hex"
iris.hex <- as.h2o(iris, destination_frame= "iris.hex")
head(iris.hex)

```

#### working w data

```{r "h2o working with data", eval = FALSE}

h2o::h2o.init()
airlines.hex <-
  h2o::h2o.importFile(path = "/tmp/allyears2k.csv" # From "https://s3.amazonaws.com/h2o-airlinesunpacked/allyears2k.csv"
                      ,destination_frame = "airlines.hex"
                      )
airlines.hex %>% class
airlines.hex %>% str
airlines.hex %>% summary

# View quantiles and histograms
quantile(x = airlines.hex$ArrDelay, na.rm = TRUE)
h2o::h2o.hist(airlines.hex$ArrDelay)

# Find number of flights by airport
originFlights <- airlines.hex %>%
  h2o::h2o.group_by(by = "Origin"
                    ,nrow("Origin")
                    ,gb.control=list(na.methods="rm")
                    )
# make the H2OFrame into a standard data.frame
# H2OFrames are worked on in the h2o framework, i.e. not your regular R-session environment.
originFlights %>% as.data.frame %>% str

# Find number of flights per month
(flightsByMonth <- airlines.hex %>%
  h2o::h2o.group_by(by = "Month"
                    ,nrow("Month")
                    ,gb.control=list(na.methods="rm"))
 )

# Find months with the highest cancellation ratio
cancellationsByMonth <- airlines.hex %>%
  h2o::h2o.group_by(by = "Month"
                    ,sum("Cancelled")
                    ,gb.control=list(na.methods="rm"))
cancellation_rate <- cancellationsByMonth$sum_Cancelled / flightsByMonth$nrow
rates_table <- h2o::h2o.cbind(flightsByMonth$Month, cancellation_rate)
rates_table %>% as.data.frame

# CONTRIBUTE
# >>>>.<<<<
# This is disgusting, how can the output list not be named >.<
# Construct test and train sets using sampling
airlines.split <- airlines.hex %>% h2o.splitFrame(ratios = 0.85)
airlines.train <- airlines.split[[1]]
airlines.test <- airlines.split[[2]]

# Display a summary using table-like functions
h2o::h2o.table(airlines.train$Cancelled)
h2o::h2o.table(airlines.test$Cancelled)

# CONTRIBUTE: the 'converting into factors' is atrocious, help out with this.

```

GLM in `h2o`

```{r "h2o.glm", eval = FALSE}

# Define the data for the model and display the results
airlines.glm <- h2o::h2o.glm(
  training_frame=airlines.train
  , x=c("Origin", "Dest", "DayofMonth", "Year", "UniqueCarrier", "DayOfWeek", "Month", "DepTime", "ArrTime", "Distance")
  , y="IsDepDelayed"
  , family = "binomial"
  , alpha = 0.5
  )
# includes confusion matrix and a lot of stuff
print(airlines.glm)
# View model information: training statistics, performance, important variables
summary(airlines.glm)

# Predict using GLM model
airlines.test[,X]
pred <- h2o::h2o.predict(object = airlines.glm, newdata = airlines.test)
# Look at summary of predictions: probability of TRUE class (p1)
# CONTRIBUTE: the Booklet is incorrect
summary(pred)

# TODO: try to read results, they look funny.
# Questions: does the h2o.glm require a numeric outcome?
# the above uses regularization due to setting alpha, is standardization of the data required beforehand?

```

```
> print(airlines.glm)
Model Details:
==============

H2OBinomialModel: glm
Model ID:  GLM_model_R_1519091204698_3
GLM Model: summary
    family  link                                regularization number_of_predictors_total number_of_active_predictors
1 binomial logit Elastic Net (alpha = 0.5, lambda = 1.482E-4 )                        283                         174
  number_of_iterations   training_frame
1                    6 RTMP_sid_ac2f_11

Coefficients: glm coefficients
      names coefficients standardized_coefficients
1 Intercept   124.608326                 -0.083354
2  Dest.ABE    -0.038939                 -0.038939
3  Dest.ABQ     0.619992                  0.619992
4  Dest.ACY     0.000000                  0.000000
5  Dest.ALB     0.000000                  0.000000

---
         names coefficients standardized_coefficients
279      Month     0.043643                  0.082146
280 DayofMonth    -0.030185                 -0.276507
281  DayOfWeek     0.026462                  0.050342
282    DepTime     0.000799                  0.372481
283    ArrTime    -0.000070                 -0.033979
284   Distance     0.000269                  0.155981

H2OBinomialMetrics: glm
** Reported on training data. **

MSE:  0.2137609
RMSE:  0.4623429
LogLoss:  0.6163326
Mean Per-Class Error:  0.3841098
AUC:  0.7192129
Gini:  0.4384258
R^2:  0.1426972
Residual Deviance:  45984.58
AIC:  46334.58

Confusion Matrix (vertical: actual; across: predicted) for F1-optimal threshold:
         NO   YES    Error          Rate
NO     6119 11576 0.654196  =11576/17695
YES    2236 17374 0.114023   =2236/19610
Totals 8355 28950 0.370245  =13812/37305

Maximum Metrics: Maximum metrics at their respective thresholds
                        metric threshold    value idx
1                       max f1  0.383832 0.715568 279
2                       max f2  0.111481 0.847229 388
3                 max f0point5  0.549731 0.681752 185
4                 max accuracy  0.529522 0.661520 197
5                max precision  0.978396 1.000000   0
6                   max recall  0.046384 1.000000 399
7              max specificity  0.978396 1.000000   0
8             max absolute_mcc  0.549731 0.324555 185
9   max min_per_class_accuracy  0.527486 0.660639 198
10 max mean_per_class_accuracy  0.549731 0.662039 185

Gains/Lift Table: Extract with `h2o.gainsLift(<model>, <data>)` or `h2o.gainsLift(<model>, valid=<T/F>, xval=<T/F>)`

```



### deepnet

- package is SUPER bare-bones
    - only works with matrices
    - does not use S3/S4 objects
    - no formulas
- autoencoder part seems to be supervised

```{r, "install deepnet", eval = FALSE}

devtools::install_github("cran/deepnet")

```

```{r "deepnet common setup"}

# only works with matrices!!
# obsEach <- 50
# x <- data.frame(var1 = c(rnorm(obsEach, 1, 0.5), rnorm(obsEach, -0.6, 0.2))
#                 ,var2 = c(rnorm(obsEach, -0.8, 0.2), rnorm(obsEach, 2, 1)))
# y <- c(rep(1, obsEach), rep(0, obsEach))

Var1 <- c(rnorm(50, 1, 0.5), rnorm(50, -0.6, 0.2))
Var2 <- c(rnorm(50, -0.8, 0.2), rnorm(50, 2, 1))
x <- matrix(c(Var1, Var2), nrow = 100, ncol = 2)
y <- c(rep(1, 50), rep(0, 50))

test_Var1 <- c(rnorm(50, 1, 0.5), rnorm(50, -0.6, 0.2))
test_Var2 <- c(rnorm(50, -0.8, 0.2), rnorm(50, 2, 1))
test_x <- matrix(c(test_Var1, test_Var2), nrow = 100, ncol = 2)
test_y <- y

```

feed-forward neural network

```{r "deepnet Training Neural Network"}

# ?deepnet::nn.train
nn <- deepnet::nn.train(x, y, hidden = c(5))
deepnet::nn.test(nn, test_x, test_y) # returns error rate for classification

```

deep belief neural network. This function internally uses `rbm.train()` to train a restricted Boltzmann machine (which can also be used individually)

The difference is mainly in the contrastive divergence algorithm that trains the restricted Boltzmann machines.
It is set via cd, giving the number of iterations for Gibbs sampling inside the learning algorithm.
[ref](http://www.rblog.uni-freiburg.de/2017/02/07/deep-learning-in-r/)

```{r "deepnet Training Deep Belief Neural Network"}

# ?deepnet::dbn.dnn.train
dbDnn <- deepnet::dbn.dnn.train(x, y, hidden = c(5, 5))

## predict by dbDnn
# Test new samples by Trainded NN,return error rate for classification
deepnet::nn.test(dbDnn, test_x, test_y)

```

Training a Deep neural network with weights initialized by Stacked AutoEncoder

```{r "deepnet Training Stacked Autoencoder Neural Network"}

# ?deepnet::sae.dnn.train
saeDnn <- deepnet::sae.dnn.train(x, y, hidden = c(5, 5))

## predict by saeDnn
# Test new samples by Trainded NN,return error rate for classification
deepnet::nn.test(saeDnn, test_x, test_y)

```

use `nn.predict` for predicting using the created from the train functions

```{r "deepnet predict"}

nnYPred <- deepnet::nn.predict(nn, test_x) %>% round(3)
dbDnnYPred <- deepnet::nn.predict(dbDnn, test_x) %>% round(4)
saeDnnYPred <- deepnet::nn.predict(saeDnn, test_x) %>% round(5)

# weight initialzer prediction comparisons
table(nnYPred, dbDnnYPred)
table(nnYPred, saeDnnYPred)
table(dbDnnYPred, saeDnnYPred)

```

### neuralnet

- does basic NNs
- overall janky user experience
    - requires numeric output
    - no dot formulas

```{r, "install neuralnet & caTools & ISLR", eval = FALSE}

devtools::install_github("cran/ISLR")
devtools::install_github("cran/caTools")
devtools::install_github("cran/neuralnet")

```

we will be using the `College` data set and try to predict the `Private` column based on the other columns / features.

NN needs scaled features. Here scale the features to [0,1]

```{r, "create data set with scaled features"}

data(College, package = "ISLR")
str(College)

unitScale <- function(dta)
  scale(dta
        ,center = dta %>% purrr::map_dbl(min, na.rm = TRUE)
        ,scale = dta %>% purrr::map(range, na.rm = TRUE) %>% purrr::map_dbl(diff)
        )

scaledCollegeFeatures <-
  College %>% dplyr::select(-Private) %>% unitScale %>% as.data.frame %>%
    dplyr::mutate(Private = College$Private %>% as.numeric %>% subtract(1))

trainObsNumbers <- scaledCollegeFeatures %>% nrow %>% sample(., .*.2, replace = FALSE)

```

Applying `neuralnet`. Here with `hidden` set to 10 neurons for each of the 3 hidden layers.

```{r "neural net function"}

frml <-
  paste("Private ~"
        ,scaledCollegeFeatures %>% names %>% setdiff("Private") %>% paste(collapse = "+")
        ) %>% as.formula

# ?neuralnet::neuralnet
nn <-
  neuralnet::neuralnet(frml
                       ,data = scaledCollegeFeatures[trainObsNumbers,]
                       ,hidden = c(10,10,10) # hidden neurons at each layer.
                       ,linear.output = FALSE
                       )


```

if want to continue and read [NN blog post](https://www.kdnuggets.com/2016/08/begineers-guide-neural-networks-r.html/2)

Plotting the output. We see the 2 hidden layers, each with 4,3 neurons / nodes. The blue text is the "bias" (intercept for the perceptron linear function).

```{r "ploting with neuralnet"}

neuralnet::neuralnet("Private ~ Enroll + Top10perc + Outstate + Books + PhD + Terminal + Expend"
                     ,data = scaledCollegeFeatures[trainObsNumbers,]
                     ,hidden = c(4,3) # hidden neurons at each layer.
                     ,linear.output = FALSE
                     ) %>% plot

```

### autoencoder NN

An autoencoder neural network is an unsupervised learning algorithm that applies backpropagation to adjust its weights, attempting to learn to make its target values (outputs) to be equal to its inputs.
In other words, it is trying to learn an approximation to the identity function, so as its output is similar to its input, for all training examples.
With the sparsity constraint enforced (requiring that the average, over training set, activation of hidden units be small), such autoencoder automatically learns useful features of the unlabeled training data, which can be used for, e.g., data compression (with losses), or as features in deep belief networks.

The training is performed by optimizing the autoencoder's cost function $J(W,b)$ that depends on the autoencoder's weights $W$ and biases $b$.
The optimization (searching for a local minimum) is performed with the `optim` function using one of the three methods: `BFGS`, `L-BFGS-B`, or `CG` (see details in `help(optim)`).

After the optimization converges, the mean squared error between the output and input matrix (either the training matrix, or a test matrix) is evaluated as a measure of goodness of fit of the autoencoder.

see package `autoencoder`

[cran](https://CRAN.R-project.org/package=autoencoder)

[wiki](https://en.wikipedia.org/wiki/Autoencoder)

```{r 'install autoencoder', eval = FALSE}

devtools::install_github("cran/autoencoder")

?autoencoder::autoencode

```

it exports 3 functions

- `autoencoder::visualize.hidden.units`
- `autoencoder::autoencode`
- `autoencoder::predict.autoencoder`

Train the autoencoder on unlabeled set of 5000 image patches of
size `Nx.patch` by `Ny.patch`, randomly cropped from 10 nature photos:
Load a training matrix with rows corresponding to training examples,
and columns corresponding to input channels (e.g., pixels in images):

```{r "autoencoder::autoencode example", eval = FALSE}

data('training_matrix_N=5e3_Ninput=100' ## the matrix contains 5e3 image
     ,package = "autoencoder")          ## patches of 10 by 10 pixels

```

Set up the autoencoder architecture:

```{r }

nl=3                          ## number of layers (default is 3: input, hidden, output)
unit.type = "logistic"        ## specify the network unit type, i.e., the unit's
                              ## activation function ("logistic" or "tanh")
Nx.patch=10                   ## width of training image patches, in pixels
Ny.patch=10                   ## height of training image patches, in pixels
N.input = Nx.patch*Ny.patch   ## number of units (neurons) in the input layer (one unit per pixel)
N.hidden = 10*10              ## number of units in the hidden layer
lambda = 0.0002               ## weight decay parameter
beta = 6                      ## weight of sparsity penalty term
rho = 0.01                    ## desired sparsity parameter
epsilon <- 0.001              ## a small parameter for initialization of weights
                              ## as small gaussian random numbers sampled from N(0,epsilon^2)
max.iterations = 2000         ## number of iterations in optimizer

```

Train the autoencoder on training.matrix using BFGS optimization method
(see help('optim') for details):
`WARNING: the training can take a long time (~1 hour) for this dataset!`


```{r "don't run!", eval = FALSE }
autoencoder.object <-
  autoencoder::autoencode(
    X.train=training.matrix
    ,nl=nl
    ,N.hidden=N.hidden
    ,unit.type=unit.type
    ,lambda=lambda
    ,beta=beta
    ,rho=rho
    ,epsilon=epsilon
    ,optim.method="BFGS"
    ,max.iterations=max.iterations
    ,rescale.flag=TRUE
    ,rescaling.offset=0.001
    )
```

N.B.: Training this autoencoder takes a long time, so in this example we do not run the above
`autoencode` function, but instead load the corresponding pre-trained `autoencoder.object`.


Report mean squared error for training and test sets:

```{r, eval = FALSE}

cat("autoencode(): mean squared error for training set: ",
round(autoencoder.object$mean.error.training.set,3),"\n")

```

Extract weights W and biases b from autoencoder.object:

```{r, eval = FALSE}

W <- autoencoder.object$W
b <- autoencoder.object$b

```

Visualize hidden units' learned features:

```{r, eval = FALSE}

visualize.hidden.units(autoencoder.object,Nx.patch,Ny.patch)

```

# Best practices {.tabset .tabset-fade .tabset-pills}

### Overview

Section for listing best programming practices in R.

### function checks

Given

```{r}
wt_mean <- function(x, w) sum(x * w) / sum(x)
wt_mean(1:6, 1:3)
```

there's gotchas with R's vector recycling rules. Use checks

```{r}
wt_mean <- function(x, w) {
  if (length(x) != length(w)) {
    stop("`x` and `w` must be the same length", call. = FALSE)
  }
  sum(w * x) / sum(x)
}
```

`x` and `w` must be the same length

```{r, eval = FALSE}
wt_mean(1:6, 1:3) # > Error: `x` and `w` must be the same length
```

you can use `stopifnot` function to do multiple checks and get informative messages back

```{r}
wt_mean <- function(x, w, na.rm = FALSE) {
  stopifnot(is.logical(na.rm), length(na.rm) == 1)
  stopifnot(length(x) == length(w))

  if (na.rm) {
    miss <- is.na(x) | is.na(w)
    x <- x[!miss]
    w <- w[!miss]
  }
  sum(w * x) / sum(x)
}
```

`is.logical(na.rm)` is not TRUE

```{r, eval = FALSE}
wt_mean(1:6, 6:1, na.rm = "foo") #> Error: is.logical(na.rm) is not TRUE
```

Note that when using `stopifnot` you assert what should be true rather than checking for what might be wrong.

# installation of packages  {.tabset .tabset-fade .tabset-pills}

## Overview

notes for installing packages

## after R upgrade

try to run `sudo R CMD javareconf` in the terminal.

```{r, eval = FALSE}
install.packages(c("devtools", "magrittr", "rJava", "dplyr"))
library(magrittr)
installFromCran <- function(pkgs) devtools::install_github(pkgs %>% paste("cran", ., sep = "/"))
partyPkgs <- c('mvtnorm', 'modeltools', 'zoo', 'sandwich', 'strucchange', 'TH.data', 'survival', 'multcomp', 'coin', 'party')
tibblePkgs <- c('rlang', 'tibble')
roxygen2Pkgs <- c('stringi', 'stringr', 'brew', 'backports', 'rprojroot', 'desc', 'commonmark', 'xml2', 'roxygen2')
testthatPkgs <- c('crayon', 'praise', 'testthat')
rattlePkgs <- c("RGtk2", "rattle")
lintrPkgs <- c( "registry", "iterators", "pkgmaker", "registry", "rngtools", "gridBase", "foreach", "doParallel"
               ,"NMF", "irlba", "rex", "stringdist", "igraph", "lintr")
rmdPkgs <- c("htmltools", "bitops", "caTools", "base64enc", "rmarkdown")
xlconnectPkgs <- c('rJava','XLConnectJars', 'XLConnect')
dplyrPkgs <- c( "bindr", "plogr", "bindrcpp", "glue", "pkgconfig", "rlang", "dplyr" )
rsqlPkgs <- c("dplyr", "DBI", "dbplyr", "RSQLServer")
knitrPkgs <- c("evaluate", "highr", "markdown", "yaml", "knitr")
c( partyPkgs, tibblePkgs, roxygen2Pkgs, testthatPkgs, rattlePkgs, lintrPkgs, rmdPkgs
  ,xlconnectPkgs, dplyrPkgs, rsqlPkgs, knitrPkgs, "tidyr", "matrixcalc", "dlm", "dbplyr") %>% installFromCran

devtools::install_github(c("emilrehnberg/party.readpaths"))
devtools::install_github(c("ramnathv/slidify"))
devtools::install_github(c("ramnathv/slidifyLibraries"))
```

No longer needed

```{r, eval = FALSE}
devtools::install_github(c( 'emilrehnberg/party.readpaths', 'jmp75/rClr'
                           ,"bescoto/RMSSQL", "bescoto/dplyr.mssql"))
```

copy below for running on server to install packages

```{r, eval = FALSE}
library(magrittr)
installFromCran <- function(pkgs) devtools::install_github(pkgs %>% paste("cran", ., sep = "/"))
dplyrPkgs <- c("bindr", "bindrcpp", "glue", "pkgconfig", "rlang", "dplyr" )
installFromCran(dplyrPkgs)
```

## tidyverse

```{r, "tidyverse installation", eval = FALSE}

c("r-dbi/DBI"
  ,"r-lib/rlang"
  ,"tidyverse/dplyr"
  ,"tidyverse/tibble"
  ,"tidyverse/purrr"
  ,"tidyverse/dbplyr"
  ,"cran/rematch" ,"cran/cellranger","tidyverse/readxl"
 ) %>% devtools::install_github()

c("r-lib/cli","r-lib/crayon"
  ,"rstudio/rstudioapi"
  ,"tidyverse/magrittr", "tidyverse/hms"
  ,"cran/BH", "mdlincoln/clipr" "tidyverse/readr"
  ,"tidyverse/forcats", "tidyverse/haven", "tidyverse/lubridate"
  ,"tidyverse/modelr"
  ,"r-lib/evaluate", "cran/yaml", "cran/markdown", "yihui/xfun", "yihui/highr", "yihui/knitr"
  ,"cran/base64enc", "eddelbuettel/digest", "cran/htmltools", "cran/jsonlite", "yihui/mime", "cran/tinytex","rstudio/rmarkdown"
  ,"r-lib/processx" , "r-lib/callr", "r-lib/fs", "edwindj/whisker", "r-lib/withr", "tidyverse/reprex"
  ) %>% devtools::install_github()

devtools::install_github( "tidyverse/tidyverse" )

```

### tibble

installing `tibble`

```{r, eval = FALSE}

depsLoc <- paste("r-lib", c("crayon", "pillar", "rlang"), sep = "/")
repoLoc <- "tidyverse/tibble"
c(depsLoc, repoLoc) %>% devtools::install_github()

```


### ggplot2

installation of `ggplot2` using `devtools`

```{r, "ggplot2 installation", eval = FALSE}

devtools::install_github(file.path("cran", c("RColorBrewer", "dichromat", "plyr", "colorspace", "munsell", "labeling", "viridisLite", "gtable", "lazyeval", "reshape2", "scales")))
devtools::install_github("tidyverse/ggplot2")

```

add-on packages

#### gifski

```{r "gifski installation", eval = FALSE}

depLoc <- "cran/png"
pkgLoc <- "r-rust/gifski"

c(depLoc, pkgLoc) %>% devtools::install_github(force = TRUE)

```


#### gganimate

```{r "gganimate installation", eval = FALSE}

tweenrLocs <- c("cran/farver", "cran/tweenr")
depsLoc <- c("r-lib/rlang", "r-lib/progress", "tidyverse/glue", tweenrLocs, "tidyverse/ggplot2")
pkgLoc <- "thomasp85/gganimate"
c(depsLoc, pkgLoc) %>% devtools::install_github()

```

#### plotly

```{r "plotly installation", eval = FALSE}

depLoc <-
  c("RcppCore/Rcpp", "rstudio/promises", "rstudio/httpuv", "rstudio/crosstalk", "Rdatatable/data.table", "cran/hexbin"
    ,"tidyverse/tibble", "r-lib/rlang", "tidyverse/purrr"
    ,"ramnathv/htmlwidgets", "r-lib/gtable", "tidyverse/ggplot2", "tidyverse/dplyr", "eddelbuettel/digest")
pkgLoc <- "ropensci/plotly"
c(depLoc, pkgLoc) %>% devtools::install_github(dependecies=TRUE)

```

## rvest

website scraper package. Needs `xml2` which looks for `xml2-config` and `pkg-config` both which has some shenanigans. Look into them and environment variables if you have issues.

```{r, "rvest installation", eval = FALSE}

c( "r-lib/xml2", "cran/httr", "hadley/rvest") %>% devtools::install_github()

```{r, eval = FALSE }
c("httpuv", "xtable", "sourcetools", "shiny") %>% paste0("cran/",.) %>% devtools::install_github()
```

## RODBC

and `RODBCext`. perform before trying to install on mac. look for other solutions on linux/mac.

```
brew update && brew install unixODBC
```

## party

install `party`.

```{r, eval = FALSE }
devtools::install_github('cran/party')
```

installing `party` dependencies.

```{r, eval = FALSE }
devtools::install_github('cran/mvtnorm')
devtools::install_github('cran/modeltools')
devtools::install_github('cran/zoo')      # sandwich dependency
devtools::install_github('cran/sandwich')
devtools::install_github('cran/strucchange')
devtools::install_github('cran/TH.data')  # multcomp dependency
devtools::install_github('cran/multcomp') # coin dependency
devtools::install_github('cran/coin')
devtools::install_github('cran/party')
```

## partykit

install `partykit`.

```{r, eval = FALSE }
devtools::install_github('cran/partykit')
```

installing `partykit` dependencies.

```{r, eval = FALSE }
devtools::install_github('cran/Formula')
```

## roxygen2

```{r, eval = FALSE }
devtools::install_github('cran/brew')
devtools::install_github('cran/backports') # rprojroot dependency
devtools::install_github('cran/rprojroot') # desc dependency
devtools::install_github('cran/desc')
devtools::install_github('cran/commonmark')
devtools::install_github('cran/xml2')
devtools::install_github('cran/roxygen2')
```

## testthat

```{r, eval = FALSE }

depsLoc <- c('r-lib/crayon', 'cran/praise')
pkgLoc <- 'r-lib/testthat'
c(depsLoc, pkgLoc) %>% devtools::install_github()

```

## rattle

```{r, eval = FALSE }
devtools::install_github('cran/rattle')
```

installing `rattle` dependencies. install `gtk+` via brew.

```{r, eval = FALSE }
devtools::install_github('cran/RGtk2')
```

## rJava

try to run `sudo R CMD javareconf` in the terminal.

```{r, eval=FALSE , echo=FALSE}
install.packages("rJava")
```

## lintr

installing lintr

```{r, eval = FALSE, echo = FALSE}
devtools::install_github("cran/registry") # pkgmaker dependency
devtools::install_github("cran/iterators") # foreach dependency
devtools::install_github( file.path("cran", c( "pkgmaker", "registry", "rngtools", "gridBase", "foreach", "doParallel" ))) # NMF dependencies
devtools::install_github( file.path("cran", c("NMF", "irlba"))) # igraph dependencies
devtools::install_github( file.path("cran", c("rex", "stringdist", "igraph")))
devtools::install_github("cran/lintr")

devtools::install_github( file.path("cran", c("rex", "stringdist", "igraph")))
devtools::install_github("cran/lintr")

```

## bestglm

installing `bestglm` package, similar to `leaps` but not limited to linear regression.

```{R "bestglm installation", eval = FALSE}

devtools::install_github(file.path("cran", c("glmnet", "grpreg")))
devtools::install_github("cran/bestglm")
library(bestglm)

```

## RSQLServer

install `rJava`.

CAREFUL! Do NOT install: `devtools::install_github('agstudy/rsqlserver')`

```{r, eval=FALSE, echo=FALSE }
install.packages("RSQLServer")
devtools::install_github('jmp75/rClr') # brew install mono # for this
devtools::install_github('imanuelcostigan/RSQLServer')
```

## dplyr.mssql

```{r, eval=FALSE , echo=FALSE}
devtools::install_github("bescoto/RMSSQL")
devtools::install_github("bescoto/dplyr.mssql") # this package
```

Also you'll need these packages: `DBI`, `dplyr`, `RJDBC`, `assertthat`, `dplyr`, `stringr`

If there's a timeout-problem, try changing download method.

```{r, eval=FALSE , echo=FALSE}
options(download.file.method = "wininet")

devtools::install_github("bescoto/dplyr.mssql", method = "wininet") # works? instead of setting the session options?
```

if there's problems with stringi, try below (the latter argument is for dealing with the case where the package has no matching version for R)

```{r, eval=FALSE , echo=FALSE}
R CMD INSTALL stringi_1.1.1.tar.gz --configure-args='--disable-pkg-config'
```

### Dependency hell

`dplyr`, `RSQLServer` and `dplyr.mssql` has this dependency hell going on ATM (161121). `dplyr.mssql` depends on `dplyr` version 4 but
latest `RSQLServer` requires >5.0. Downgrade `RSQLServer` and `dplyr` to make `dplyr.mssql` happy.

```{sh, eval=FALSE, echo=FALSE}
wget https://github.com/imanuelcostigan/RSQLServer/archive/v0.2.0.tar.gz
R CMD INSTALL v0.2.0.tar.gz

wget https://cran.rstudio.com/src/contrib/Archive/dplyr/dplyr_0.4.3.tar.gz
R CMD INSTALL dplyr_0.4.3.tar.gz

```

# workflows

### Logging {.tabset .tabset-fade .tabset-pills}

#### Overview

Section for understanding logging and directing output streams.

defaults:

```{r "default direction"}
1:2 # goes to stdout()
message("messaging") # goes to stderr()
warning("warning") # goes to stderr()
```

to redirect message to standard output stream

```{r "redirect message to stdout", eval = FALSE}
message("more messaging") # going to stderr()
sink(stdout(), type = "message")
message("messaging to std out!") # goes to stdout()
```

redirecting regular output to standard error

```{r "redirecting output stderr", eval = FALSE}
1:2 # to stdout()
sink(stderr(), type = "output")
3:4 # to stderr()
```

#### redirecting IO

all IO, stdout, stderr and warnings are sent to standard output stream.

start with `#!/bin/usr/env Rscript`

```{r, eval = FALSE}
sink(stdout(), type="message") # MAIN command for redirection
message("message: testing 12")
warning("achtung: deutchen furher")
error("E$ failure") # doesn't work, using `stop` if you want to send error messages and stop the execution.
```

[source](http://mazamascience.com/WorkingWithData/?p=888)

### variable relations {.tabset .tabset-fade .tabset-pills}

#### Overview

Section for investigating relations between variables

See plots section for pairwise correlation plots.

#### pairwise correlations

`stats::cor `

only works on numerical variables

```{r stats_cor_pairwise}
cor(airquality[,c("Ozone", "Wind", "Temp")],
    use = "complete")
```

#### descriptives over strata

stratum / labels / levels

```{r strata_descriptives}
with(CO2, {
  aggregate(uptake
            ,by = list(Type = Type, Treatment = Treatment)
            ,FUN = summary)
})

with(CO2, {
  aggregate(CO2[,c("uptake", "conc")]
            ,by = list(Type = Type, Treatment = Treatment)
            ,FUN = median)
})
```

`dplyr` version might be more performant for single values, but you can't use `summary` like above.

```{r strata_descriptives_dplyr}
library(dplyr)
( CO2
  %>% group_by(Type, Treatment)
  %>% summarise(uptakeMean = mean(uptake), updateMedian = median(uptake), concMean = mean(conc))
  )
```

#### fill out missing combinations / gaps

if there's gaps in the data where there should data, maybe 0 data, you can join in a dataset with full combinations.

```{r}
library(dplyr)
incompleteDt <-
  data.frame( cdate = c(201609, 201610, 201610, 201611, 201612, 201612)
             ,status = c("Bad", "OK", "Bad", "OK", "OK", "Bad")
             ,metric = c(runif(6)))
completeLabels <- expand.grid(cdate = c(201609, 201610, 201611, 201612), status = c("Bad", "OK"))
incompleteDt %>% right_join(completeLabels) %>% mutate(metric = ifelse(is.na(metric), 0, metric))
```

#### clustering with mixed variable types

taken from [article](https://www.r-bloggers.com/clustering-mixed-data-types-in-r/)

```{r}
library(ISLR)
College %>% str
```

### A/B testing {.tabset .tabset-fade .tabset-pills}

#### Overview

Section for A/B testing workflow

#### Bayesian A/B testing

main function is th `bayesTest` function. Which takes two samples (A and B), conjugate prior parameters (see `?bayesTest` for expected named vectors of parameters) and the distribution for likelihood.

```{r "plots setup"}
opar <- par(mfrow = c(2, 2))
```

##### with Bernoulli posterior

```{r, "bayesAB installation", eval = FALSE}

devtools::install_github("cran/bayesAB")

```

```{r "bayesAB package Bernoulli posterior"}

nSamples <- 250
A_binom <- rbinom(nSamples, size = 1, prob = .25) # size = 1 makes sure it's a bernoulli trial, so the outcome is either 0 or 1
B_binom <- rbinom(nSamples, size = 1, prob = .2)

```

- let's ignore that we know the samples in this case
- we know the Bernoulli distribution probability is ~0.2-0.3.
- with this in mind, we want prior parameters for a Beta that roughly matces these (Beta dist is conjugate prior to Binomial)

```{r}

bayesAB::plotBeta(100, 200) # looks a little off
bayesAB::plotBeta(65, 200) # looks way better

(AB1 <-
  bayesAB::bayesTest( A_binom
                     ,B_binom
                     ,priors = c('alpha' = 65, 'beta' = 200) # params for a Beta (conjugate to Binomial)
                     ,distribution = 'bernoulli'
                     ))
AB1 %>% summary
AB1 %>% plot

```

##### with Poisson posterior

```{r "bayesAB package Poisson posterior"}

A_pois <- rpois(nSamples, 6.5)
B_pois <- rpois(nSamples, 5.5)

bayesAB::plotGamma(30, 5) # 5-7 seem likely enough

(AB2 <-
  bayesAB::bayesTest( A_pois
                     ,B_pois
                     ,priors = c(shape = 30, rate = 5)
                     ,distribution = 'poisson'
                     ))
AB2 %>% summary
AB2 %>% plot

```

##### combining distribution

The combination here refers to combining the events in the Bernoulli and Poisson examples. The example being the first example is the click-through event and the second one being the number of clicks on the page. But the second example is based on the situation that the lead is already on the page. What about the expected number of clicks from the get-go? Then you'd have to combine the first and second example somehow and luckily `bayesAB::combine` does this.

```{r "bayesAB::combine"}
(AB3 <-
  bayesAB::combine( AB1
                   ,AB2
                   ,f = multiply_by
                   ,params = c('Probability', 'Lambda') # the kind of parameters for the posteriors in AB1/2
                   ,newName = 'Expectation'
                   )
 )
AB3 %>% summary # note the lower expecations on clicks here
AB3 %>% plot
```

```{r "finish plot"}
par(opar)
```

# References