Software for Statistics

• Computing software is essential for modern statistics

  • Large datasets

  • Visualization

  • Simulation

  • Iterative methods

• Many softwares are available

Outline

• Installing R

• Some examples

• A little bit of history

• Some thoughts on why R has been successful

Installing R

• R is most commonly used as a REPL (Read-Eval-Print-Loop)

• This is essentially the model used by a calculator:

  • Waits for user input

  • Evaluates and prints result

  • Waits for more input

Installing R

• Go to https://cran.r-project.org/ (or choose a mirror first)

• Follow instructions depending on your platform (probably Windows)

Running R

• Once installed, you can start the appropriate interface (or R directly) to get something like this:
  
  

R Under development (unstable) (2018-05-05 r74699) -- "Unsuffered Consequences"
Copyright (C) 2018 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

Loading required package: utils
> 

• The > represents a prompt indicating that R is waiting for input.

• The difficult part is to learn what to do next

Before we start, an experiment!

Question: What proportion of population sees white & gold?

• Statistics uses data to make inferences

• Model:

  • Let \(p\) be the probability of seeing white & gold

  • Assume that individuals are independent

Generally useful method: maximum likelihood

• Likelihood function: probability of observed data as function of \(p\)

\[ L(p) = P(X = 26) = {44 \choose 26} p^{26} (1-p)^{(44-26)}, p \in (0, 1) \]

• Intuition: \(p\) that gives higher \(L(p)\) is more “likely” to be correct

• Maximum likelihood estimate \(\hat{p} = \arg \max L(p)\)

How could we do this numerically?

• Pretend for the moment that we did not know how to do this.

• How could we arrive at the same solution numerically?

• Basic idea: Compute \(L(p)\) for various values of \(p\) and find minimum.

N = 44
x = 26

p = 0.5
choose(N, x) * p^x * (1-p)^(N-x)

[1] 0.05852204

“Vectorized” computations

• One distinguishing feature of R is that it operates on “vectors”

pvec = seq(0, 1, by = 0.01)
pvec

  [1] 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.10 0.11 0.12 0.13 0.14 0.15 0.16 0.17 0.18 0.19 0.20 0.21 0.22
 [24] 0.23 0.24 0.25 0.26 0.27 0.28 0.29 0.30 0.31 0.32 0.33 0.34 0.35 0.36 0.37 0.38 0.39 0.40 0.41 0.42 0.43 0.44 0.45
 [47] 0.46 0.47 0.48 0.49 0.50 0.51 0.52 0.53 0.54 0.55 0.56 0.57 0.58 0.59 0.60 0.61 0.62 0.63 0.64 0.65 0.66 0.67 0.68
 [70] 0.69 0.70 0.71 0.72 0.73 0.74 0.75 0.76 0.77 0.78 0.79 0.80 0.81 0.82 0.83 0.84 0.85 0.86 0.87 0.88 0.89 0.90 0.91
 [93] 0.92 0.93 0.94 0.95 0.96 0.97 0.98 0.99 1.00

Lvec = choose(N, x) * pvec^x * (1-pvec)^(N-x)
Lvec

  [1] 0.000000e+00 8.591575e-41 4.802734e-33 1.512457e-28 2.223726e-25 6.093745e-23 5.765981e-21 2.617468e-19
  [9] 6.936811e-18 1.218119e-16 1.545270e-15 1.506153e-14 1.180429e-13 7.700395e-13 4.294774e-12 2.091957e-11
 [17] 9.052864e-11 3.529530e-10 1.254220e-09 4.101694e-09 1.244626e-08 3.528813e-08 9.404416e-08 2.368078e-07
 [25] 5.659476e-07 1.288790e-06 2.806191e-06 5.860149e-06 1.176882e-05 2.278440e-05 4.261443e-05 7.714841e-05
 [33] 1.354251e-04 2.308597e-04 3.827207e-04 6.178014e-04 9.721737e-04 1.492843e-03 2.239047e-03 3.282888e-03
 [41] 4.708923e-03 6.612349e-03 9.095461e-03 1.226215e-02 1.621039e-02 2.102292e-02 2.675658e-02 3.343099e-02
 [49] 4.101773e-02 4.943113e-02 5.852204e-02 6.807589e-02 7.781593e-02 8.741246e-02 9.649794e-02 1.046874e-01
 [57] 1.116031e-01 1.169009e-01 1.202969e-01 1.215909e-01 1.206845e-01 1.175920e-01 1.124418e-01 1.054689e-01
 [65] 9.699819e-02 8.742011e-02 7.716176e-02 6.665536e-02 5.630807e-02 4.647572e-02 3.744302e-02 2.941171e-02
 [73] 2.249722e-02 1.673329e-02 1.208326e-02 8.455753e-03 5.722622e-03 3.736794e-03 2.348049e-03 1.415438e-03
 [81] 8.156783e-04 4.475222e-04 2.326508e-04 1.139594e-04 5.224689e-05 2.224201e-05 8.707704e-06 3.098277e-06
 [89] 9.873047e-07 2.765972e-07 6.651882e-08 1.330702e-08 2.121986e-09 2.540743e-10 2.092599e-11 1.034935e-12
 [97] 2.447773e-14 1.806704e-16 1.596089e-19 7.927831e-25 0.000000e+00

Plotting is very easy

plot(x = pvec, y = Lvec, type = "l")

Functions

• Functions can be used to encapsulate repetitive computations

• Like mathematical functions, R function also take arguments as input and “returns” an output

L = function(p) choose(N, x) * p^x * (1-p)^(N-x)
L(0.5)

[1] 0.05852204

L(x/N)

[1] 0.1216

Functions can be plotted directly

plot(L, from = 0, to = 1)

…and they can be numerically “optimized”

optimize(L, interval = c(0, 1), maximum = TRUE)

$maximum
[1] 0.5909084

$objective
[1] 0.1216

• Compare with

x / N

[1] 0.5909091

A more complicated example

• Suppose \(X_1, X_2, ..., X_n \sim Bin(N, p)\), and are independent

• Instead of observing each \(X_i\), we only get to know \(M = \max(X_1, X_2, ..., X_n)\)

• What is the maximum likelihood estimate of \(p\)? (\(N\) and \(n\) are known, \(M = m\) is observed)

A more complicated example

To compute likelihood, we need p.m.f. of \(M\) : \[ P(M \leq m) = P(X_1 \leq m, ..., X_n \leq m) = \left[ \sum_{x=0}^m {N \choose x} p^{x} (1-p)^{(N-x)} \right]^n \] and \[ P(M = m) = P(M \leq m) - P(M \leq m-1) \]

n = 10
N = 50
M = 30
F <- function(p, m)
{
    x = seq(0, m)
    (sum(choose(N, x) * p^x * (1-p)^(N-x)))^n
}
L = function(p)
{
    F(p, M) - F(p, M-1)
}

Maximum Likelihood estimate

optimize(L, interval = c(0, 1), maximum = TRUE)

$maximum
[1] 0.4996703

$objective
[1] 0.1981222

“The Dress” revisited

• What factors determine perceived color? (From 23andme.com)

Simulation: birthday problem

• R can be used to simulate random events

• Example: how likely is a common birthday in a group of 20 people?

N = 20
days = sample(365, N, rep = TRUE)
days

 [1] 112 320  19  42  66  41  73 182 314 266 154 313 351 276 218 359 257 246 195  42

length(unique(days))

[1] 19

Law of Large Numbers

• With enough replications, sample proportion should converge to probability

haveCommon = function()
{
    days = sample(365, N, rep = TRUE)
    length(unique(days)) < N
}
haveCommon()

[1] FALSE

haveCommon()

[1] FALSE

haveCommon()

[1] TRUE

haveCommon()

[1] TRUE

Law of Large Numbers

• With enough replications, sample proportion should converge to probability

• Do this sytematically:

replicate(100, haveCommon())

  [1] FALSE FALSE FALSE  TRUE FALSE FALSE  TRUE  TRUE  TRUE FALSE  TRUE FALSE FALSE FALSE  TRUE  TRUE FALSE  TRUE  TRUE
 [20]  TRUE FALSE  TRUE  TRUE  TRUE FALSE FALSE  TRUE FALSE FALSE FALSE  TRUE FALSE  TRUE FALSE  TRUE  TRUE FALSE FALSE
 [39]  TRUE FALSE FALSE  TRUE  TRUE FALSE  TRUE FALSE FALSE  TRUE  TRUE FALSE  TRUE  TRUE FALSE  TRUE FALSE FALSE FALSE
 [58]  TRUE FALSE  TRUE FALSE FALSE FALSE FALSE  TRUE FALSE  TRUE FALSE FALSE  TRUE FALSE FALSE FALSE  TRUE FALSE FALSE
 [77] FALSE  TRUE  TRUE FALSE FALSE FALSE  TRUE FALSE  TRUE FALSE FALSE FALSE FALSE  TRUE FALSE FALSE  TRUE FALSE FALSE
 [96] FALSE  TRUE FALSE FALSE FALSE

Law of Large Numbers

• With enough replications, sample proportion should converge to probability

plot(cumsum(replicate(1000, haveCommon())) / 1:1000, type = "l")
lines(cumsum(replicate(1000, haveCommon())) / 1:1000, col = "red")
lines(cumsum(replicate(1000, haveCommon())) / 1:1000, col = "blue")

A more serious example: climate change

Change in temperature (global average deviation) since 1851

library(lattice)
xyplot(Temp ~ Year, data = globalTemp, grid = TRUE)

Change in atmospheric carbon dioxide

xyplot(CO2 ~ Year, data = globalTemp, grid = TRUE)

Does change in \(CO_2\) explain temperature rise?

xyplot(Temp ~ CO2, data = globalTemp, grid = TRUE, type = c("p", "r")) # include OLS regression line

Fitting the regression model

fm = lm(Temp ~ 1 + CO2, data = globalTemp)
coef(fm) # estimated regression coefficients

 (Intercept)          CO2 
-2.836082117  0.008486628 

SSE = function(beta)
{
    with(globalTemp,
         sum((Temp - beta[1] - beta[2] * CO2)^2))
}
optim(c(0, 0), fn = SSE)

$par
[1] -2.836176636  0.008486886

$value
[1] 2.210994

$counts
function gradient 
      93       NA 

$convergence
[1] 0

$message
NULL

Fitting the regression model

• lm() gives exact solution and more statistically relevant details

summary(fm)

Call:
lm(formula = Temp ~ 1 + CO2, data = globalTemp)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.28460 -0.09004 -0.00101  0.08616  0.35926 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)
(Intercept) -2.8360821  0.1145766  -24.75   <2e-16
CO2          0.0084866  0.0003602   23.56   <2e-16

Residual standard error: 0.1218 on 149 degrees of freedom
Multiple R-squared:  0.7884,    Adjusted R-squared:  0.787 
F-statistic: 555.1 on 1 and 149 DF,  p-value: < 2.2e-16

Changing the model-fitting criteria

• Suppose we wanted to minimize sum of absolute errors instead of sum of squares

• No closed form solution any more, but general optimizer will still work:

SAE = function(beta)
{
    with(globalTemp,
         sum(abs(Temp - beta[1] - beta[2] * CO2)))
}
opt = optim(c(0, 0), fn = SAE)
opt

$par
[1] -2.832090898  0.008471257

$value
[1] 14.5602

$counts
function gradient 
     123       NA 

$convergence
[1] 0

$message
NULL

Changing the model-fitting criteria

• Compare with least squares line

coef(fm) # least squared errors

 (Intercept)          CO2 
-2.836082117  0.008486628 

opt$par  # least absolute errors

[1] -2.832090898  0.008471257

• The two lines are virtually identical in this case

• This is not always true

Another example: number of phone calls per year in Belgium

data(phones, package = "MASS")
xyplot(calls ~ year, data = phones, grid = TRUE)

Another example: number of phone calls per year in Belgium

fm2 <- lm(calls ~ year, data = phones)
SAE = function(beta)
{
    with(phones,
         sum(abs(calls - beta[1] - beta[2] * year)))
}
opt = optim(c(0, 0), fn = SAE)

coef(fm2) # least squared errors

(Intercept)        year 
-260.059246    5.041478 

opt$par  # least absolute errors

[1] -66.053297   1.353735

Another example: number of phone calls per year in Belgium

xyplot(calls ~ year, data = phones, grid = TRUE,
       panel = function(x, y, ...) {
           panel.xyplot(x, y, ...)
           panel.abline(fm2, col = "red") # least squared errors
           panel.abline(opt$par, col = "blue") # least absolute errors
       })

Summary

• Conventional statistical learning focuses on problems that can be “solved” analytically