A Crazy Little Thing Called {purrr} - Part 6 : doing statistics

This is gonna be the last time I make this Queen reference in 2017.

I’ve been sharing some “stats with {purrr}” recipes on my Twitter account lately. Twitter being what it is (a source of infinite temporary content), here’s the post gathering some statistics tricks you can do with {purrr}.

Looking for normality

If your data are in a data.frame, you can simply map a shapiro.test on all the columns, keeping the one with a p.value > to 0.05 (remember, the shapiro tests for non-normality) :

library(purrr)
map(airquality, shapiro.test) %>% keep(~ .x$p.value > 0.05)

## $Wind
## 
##  Shapiro-Wilk normality test
## 
## data:  .x[[i]]
## W = 0.98575, p-value = 0.1178

Of course, you can do the same on a non-tabular list:

# rdunif is from purrr
l <- list(a = rnorm(10), b = rdunif(1000, 10))
map(l, shapiro.test) %>% keep(~ .x$p.value > 0.05)

## $a
## 
##  Shapiro-Wilk normality test
## 
## data:  .x[[i]]
## W = 0.91275, p-value = 0.3004

Also, map_if allows you to map only on numeric variables in your data.frame:

map_if(iris, is.numeric, shapiro.test)

## $Sepal.Length
## 
##  Shapiro-Wilk normality test
## 
## data:  .x[[i]]
## W = 0.97609, p-value = 0.01018
## 
## 
## $Sepal.Width
## 
##  Shapiro-Wilk normality test
## 
## data:  .x[[i]]
## W = 0.98492, p-value = 0.1012
## 
## 
## $Petal.Length
## 
##  Shapiro-Wilk normality test
## 
## data:  .x[[i]]
## W = 0.87627, p-value = 7.412e-10
## 
## 
## $Petal.Width
## 
##  Shapiro-Wilk normality test
## 
## data:  .x[[i]]
## W = 0.90183, p-value = 1.68e-08
## 
## 
## $Species
##   [1] setosa     setosa     setosa     setosa     setosa     setosa    
##   [7] setosa     setosa     setosa     setosa     setosa     setosa    
##  [13] setosa     setosa     setosa     setosa     setosa     setosa    
##  [19] setosa     setosa     setosa     setosa     setosa     setosa    
##  [25] setosa     setosa     setosa     setosa     setosa     setosa    
##  [31] setosa     setosa     setosa     setosa     setosa     setosa    
##  [37] setosa     setosa     setosa     setosa     setosa     setosa    
##  [43] setosa     setosa     setosa     setosa     setosa     setosa    
##  [49] setosa     setosa     versicolor versicolor versicolor versicolor
##  [55] versicolor versicolor versicolor versicolor versicolor versicolor
##  [61] versicolor versicolor versicolor versicolor versicolor versicolor
##  [67] versicolor versicolor versicolor versicolor versicolor versicolor
##  [73] versicolor versicolor versicolor versicolor versicolor versicolor
##  [79] versicolor versicolor versicolor versicolor versicolor versicolor
##  [85] versicolor versicolor versicolor versicolor versicolor versicolor
##  [91] versicolor versicolor versicolor versicolor versicolor versicolor
##  [97] versicolor versicolor versicolor versicolor virginica  virginica 
## [103] virginica  virginica  virginica  virginica  virginica  virginica 
## [109] virginica  virginica  virginica  virginica  virginica  virginica 
## [115] virginica  virginica  virginica  virginica  virginica  virginica 
## [121] virginica  virginica  virginica  virginica  virginica  virginica 
## [127] virginica  virginica  virginica  virginica  virginica  virginica 
## [133] virginica  virginica  virginica  virginica  virginica  virginica 
## [139] virginica  virginica  virginica  virginica  virginica  virginica 
## [145] virginica  virginica  virginica  virginica  virginica  virginica 
## Levels: setosa versicolor virginica

Bulk cor.test

As said on Twitter, this might be p.hacking, but here’s how you can do a cor.test for all columns in a data.frame, with a little help from tidystringdist and broom :

library(tidystringdist) # Works since v0.1.2 
comb <- tidy_comb_all(names(airquality))
knitr::kable(comb)

V1	V2
Ozone	Solar.R
Ozone	Wind
Ozone	Temp
Ozone	Month
Ozone	Day
Solar.R	Wind
Solar.R	Temp
Solar.R	Month
Solar.R	Day
Wind	Temp
Wind	Month
Wind	Day
Temp	Month
Temp	Day
Month	Day

bulk_cor <- pmap(comb, ~ cor.test(airquality[[.x]], airquality[[.y]])) %>% 
  map_df(broom::tidy) %>% 
  cbind(comb, .)

knitr::kable(bulk_cor, digits = 3)

V1 | V2 | estimate | statistic | p.value | parameter | conf.low | conf.high | method | alternative |
:– | :– | —: | : | –: | : | —: | : | :— | :– |
Ozone | Solar.R | 0.348 | 3.880 | 0.000 | 109 | 0.173 | 0.502 | Pearson’s product-moment correlation | two.sided |
Ozone | Wind | -0.602 | -8.040 | 0.000 | 114 | -0.706 | -0.471 | Pearson’s product-moment correlation | two.sided |
Ozone | Temp | 0.698 | 10.418 | 0.000 | 114 | 0.591 | 0.781 | Pearson’s product-moment correlation | two.sided |
Ozone | Month | 0.165 | 1.781 | 0.078 | 114 | -0.018 | 0.337 | Pearson’s product-moment correlation | two.sided |
Ozone | Day | -0.013 | -0.141 | 0.888 | 114 | -0.195 | 0.169 | Pearson’s product-moment correlation | two.sided |
Solar.R | Wind | -0.057 | -0.683 | 0.496 | 144 | -0.217 | 0.107 | Pearson’s product-moment correlation | two.sided |
Solar.R | Temp | 0.276 | 3.444 | 0.001 | 144 | 0.119 | 0.419 | Pearson’s product-moment correlation | two.sided |
Solar.R | Month | -0.075 | -0.906 | 0.366 | 144 | -0.235 | 0.088 | Pearson’s product-moment correlation | two.sided |
Solar.R | Day | -0.150 | -1.824 | 0.070 | 144 | -0.305 | 0.012 | Pearson’s product-moment correlation | two.sided |
Wind | Temp | -0.458 | -6.331 | 0.000 | 151 | -0.575 | -0.323 | Pearson’s product-moment correlation | two.sided |
Wind | Month | -0.178 | -2.227 | 0.027 | 151 | -0.328 | -0.020 | Pearson’s product-moment correlation | two.sided |
Wind | Day | 0.027 | 0.334 | 0.739 | 151 | -0.132 | 0.185 | Pearson’s product-moment correlation | two.sided |
Temp | Month | 0.421 | 5.703 | 0.000 | 151 | 0.281 | 0.543 | Pearson’s product-moment correlation | two.sided |
Temp | Day | -0.131 | -1.619 | 0.108 | 151 | -0.283 | 0.029 | Pearson’s product-moment correlation | two.sided |
Month | Day | -0.008 | -0.098 | 0.922 | 151 | -0.166 | 0.151 | Pearson’s product-moment correlation | two.sided |

Some Machine Learning

lm

Let’s build a linear model of all possible iris combinations:

res <- pmap(comb, ~ lm(airquality[[.x]] ~ airquality[[.y]]))
get_rsquared <- compose(as_mapper(~ .x$r.squared), summary)
map_dbl(res, get_rsquared)

##  [1] 1.213419e-01 3.618582e-01 4.877072e-01 2.706660e-02 1.749177e-04
##  [6] 3.225293e-03 7.608786e-02 5.670205e-03 2.258257e-02 2.097529e-01
## [11] 3.178824e-02 7.388015e-04 1.771966e-01 1.705458e-02 6.338966e-05

Any chance there’s a r.squared above 0.5 ?

map(res, get_rsquared) %>% some(~ .x > 0.5)

## [1] FALSE

rpart

We’ll build 20 rpart to find the better model. Yes, this will be better to do it with a random forest, but we’re here for the example :)

Training and validation

Let’s do the train / validation.

library(dplyr)
library(readr)
titanic <- read_csv("http://biostat.mc.vanderbilt.edu/wiki/pub/Main/DataSets/titanic3.csv")
train <- rerun(20, sample_frac(titanic, size = 0.8))
validation <- map(train, ~ anti_join(titanic, .x))

Check if the sizes of all validation datasets are the same:

map_int(validation, nrow) %>% every(~ .x == 262)

## [1] TRUE

Create a model builder

We’ll create a simple model to predict survival based on sex.

library(rpart)
rpart_pimped <- partial(rpart, formula = survived ~ sex, method = "class")
res <- map(train, ~ rpart_pimped(data = .x))

Make prediction

prediction <- map2(validation, res, ~ predict(.y, .x, type = "class"))
w_prediction <- map2(validation, prediction, ~ mutate(.x, prediction = .y))

Confusion matrix

Now let’s map a conf matrix on all these results:

library(caret)
conf_mats <- map(w_prediction, ~ confusionMatrix(.x$prediction, .x$survived))

Is the “Sensivity” for all models above 0.8?

map_dbl(conf_mats, ~ .x$byClass["Sensitivity"]) %>% every(~ .x > 0.8)

## [1] TRUE

Is the “Specificity” for all models above 0.8?

map_dbl(conf_mats, ~ .x$byClass["Specificity"]) %>% every(~ .x > 0.8)

## [1] FALSE

keep_index

Let’s modify keep a little bit so we can extract the models we need:

# Here's the keep source code
keep

## function(.x, .p, ...) {
##   sel <- probe(.x, .p, ...)
##   .x[!is.na(sel) & sel]
## }
## <environment: namespace:purrr>

# Let's tweak it a little bit
keep_index <- function(.x, .p, ...) {
  sel <- purrr:::probe(.x, .p, ...)
  which(sel)
}

So, which are the models with a sensitivity superior to 0.85?

map_dbl(conf_mats, ~ .x$byClass["Sensitivity"]) %>% keep_index(~ .x > 0.85)

##  [1]  1  3  4  6 11 13 14 15 16 18 19

And the models with a specificity superior to 0.7?

map_dbl(conf_mats, ~ .x$byClass["Specificity"]) %>% keep_index(~ .x > 0.7)

## [1]  2  7 13 14 17

Which models are in both?

sens <- map_dbl(conf_mats, ~ .x$byClass["Sensitivity"]) %>% keep_index(~ .x > 0.85)
spec <- map_dbl(conf_mats, ~ .x$byClass["Specificity"]) %>% keep_index(~ .x > 0.7)
keep(sens, map_lgl(sens, ~ .x %in% spec))

## [1] 13 14

So, I guess we’ll go for model(s) number 13, 14!