How to draw a stratified sample in R

Problem

We want to draw a stratified sample from a data frame in R.

Solution

Let's look at two examples, with one or several groups.

One group

We extract 3 records from each of the species: setosa, versicolor y virginica.

Base package

set.seed(1)
iris1 <- lapply(split(iris, iris$Species), function(x) x[sample(nrow(x), 3), ])
do.call("rbind", iris1) 

              Sepal.Length Sepal.Width Petal.Length Petal.Width    Species
setosa.14              4.3         3.0          1.1         0.1     setosa
setosa.19              5.7         3.8          1.7         0.3     setosa
setosa.28              5.2         3.5          1.5         0.2     setosa
versicolor.96          5.7         3.0          4.2         1.2 versicolor
versicolor.60          5.2         2.7          3.9         1.4 versicolor
versicolor.94          5.0         2.3          3.3         1.0 versicolor
virginica.148          6.5         3.0          5.2         2.0  virginica
virginica.133          6.4         2.8          5.6         2.2  virginica
virginica.131          7.4         2.8          6.1         1.9  virginica

dplyr

library(dplyr)
set.seed(1)
iris %>%
  group_by(Species) %>%
  sample_n(., 3)

 Source: local data frame [9 x 5]
Groups: Species

  Sepal.Length Sepal.Width Petal.Length Petal.Width    Species
1          4.3         3.0          1.1         0.1     setosa
2          5.7         3.8          1.7         0.3     setosa
3          5.2         3.5          1.5         0.2     setosa
4          5.7         3.0          4.2         1.2 versicolor
5          5.2         2.7          3.9         1.4 versicolor
6          5.0         2.3          3.3         1.0 versicolor
7          6.5         3.0          5.2         2.0  virginica
8          6.4         2.8          5.6         2.2  virginica
9          7.4         2.8          6.1         1.9  virginica

Two groups

For each number of cylinders (4, 6 u 8) we will extract 2 records with automatic transmission (am = 0) and 2 with manual transmission (am = 1).

Base package

set.seed(1)
mtcars1 <- lapply(split(mtcars, list(mtcars$cyl, mtcars$am)), function(x) x[sample(nrow(x), 2), ])
do.call("rbind", mtcars1) 

                       mpg cyl  disp  hp drat    wt  qsec vs am gear carb
0.4.Merc 240D         24.4   4 146.7  62 3.69 3.190 20.00  1  0    4    2
0.4.Toyota Corona     21.5   4 120.1  97 3.70 2.465 20.01  1  0    3    1
1.4.Fiat X1-9         27.3   4  79.0  66 4.08 1.935 18.90  1  1    4    1
1.4.Lotus Europa      30.4   4  95.1 113 3.77 1.513 16.90  1  1    5    2
0.6.Hornet 4 Drive    21.4   6 258.0 110 3.08 3.215 19.44  1  0    3    1
0.6.Merc 280          19.2   6 167.6 123 3.92 3.440 18.30  1  0    4    4
1.6.Ferrari Dino      19.7   6 145.0 175 3.62 2.770 15.50  0  1    5    6
1.6.Mazda RX4 Wag     21.0   6 160.0 110 3.90 2.875 17.02  0  1    4    4
0.8.Chrysler Imperial 14.7   8 440.0 230 3.23 5.345 17.42  0  0    3    4
0.8.Hornet Sportabout 18.7   8 360.0 175 3.15 3.440 17.02  0  0    3    2
1.8.Ford Pantera L    15.8   8 351.0 264 4.22 3.170 14.50  0  1    5    4
1.8.Maserati Bora     15.0   8 301.0 335 3.54 3.570 14.60  0  1    5    8

dplyr

set.seed(1)
mtcars %>%
  group_by(cyl, am) %>%
  sample_n(., 2)

Source: local data frame [12 x 11]
Groups: cyl, am

    mpg cyl  disp  hp drat    wt  qsec vs am gear carb
1  24.4   4 146.7  62 3.69 3.190 20.00  1  0    4    2
2  21.5   4 120.1  97 3.70 2.465 20.01  1  0    3    1
3  27.3   4  79.0  66 4.08 1.935 18.90  1  1    4    1
4  30.4   4  95.1 113 3.77 1.513 16.90  1  1    5    2
5  21.4   6 258.0 110 3.08 3.215 19.44  1  0    3    1
6  19.2   6 167.6 123 3.92 3.440 18.30  1  0    4    4
7  19.7   6 145.0 175 3.62 2.770 15.50  0  1    5    6
8  21.0   6 160.0 110 3.90 2.875 17.02  0  1    4    4
9  14.7   8 440.0 230 3.23 5.345 17.42  0  0    3    4
10 18.7   8 360.0 175 3.15 3.440 17.02  0  0    3    2
11 15.8   8 351.0 264 4.22 3.170 14.50  0  1    5    4
12 15.0   8 301.0 335 3.54 3.570 14.60  0  1    5    8

Related posts

• How to select a random sample in R

References

• stackoverflow

How to draw a stratified sample in R

Problem

We want to draw a stratified sample in R. Previously, we took a random sample from a data frame. We did not control over the distribution of the subgroups. This time we will control over the distribution of each stratum keeping the same overall distribution of the original data.

Solution

Using the function createDataPartition from the caret package.

library(tidyverse)
library(caret)
set.seed(1)
planes <- as.data.frame(nycflights13::planes)
trainIndex <- createDataPartition(planes$engine,
                                  p = .5,
                                  list = FALSE,
                                  times = 1)

planesTrain <- planes[trainIndex, ]
planesTest  <- planes[-trainIndex, ]

Using the function stratified from the splitstackshape package.

library(splitstackshape)
set.seed(1)
planesTrain1 <- stratified(planes, "engine", 0.5)
planesTest1 <- planes[!(planes$tailnum %in% planesTrain$tailnum),]

Checking that we have created balanced splits of the data.

• Original data frame

planes %>% 
  group_by(engine) %>% 
    summarise(n = n()) %>%
  mutate(cum = n / sum(n))

# A tibble: 6 x 3
  engine            n      cum
               
1 4 Cycle           2 0.000602
2 Reciprocating    28 0.00843 
3 Turbo-fan      2750 0.828   
4 Turbo-jet       535 0.161   
5 Turbo-prop        2 0.000602
6 Turbo-shaft       5 0.00151  

• Training set

planesTrain %>%
  group_by(engine) %>%
  summarise(n = n()) %>%
  mutate(cum = n / sum(n))

# A tibble: 6 x 3
  engine            n      cum
               
1 4 Cycle           1 0.000602
2 Reciprocating    14 0.00842 
3 Turbo-fan      1375 0.827   
4 Turbo-jet       268 0.161   
5 Turbo-prop        1 0.000602
6 Turbo-shaft       3 0.00181 

• Test set

planesTest %>%
  group_by(engine) %>%
  summarise(n = n()) %>%
  mutate(cum = n / sum(n))

# A tibble: 6 x 3
  engine            n      cum
               
1 4 Cycle           1 0.000602
2 Reciprocating    14 0.00843 
3 Turbo-fan      1375 0.828   
4 Turbo-jet       267 0.161   
5 Turbo-prop        1 0.000602
6 Turbo-shaft       2 0.00120 

Related posts

• How to select a random sample in R
• splitstackshape

References

• The caret package

How to select a random sample in R

Problem

We want to extract a random sample from a data frame in R.

Solution

• Base package

set.seed(1)
starwars[sample(nrow(starwars), 10), ] # 10 filas
# Showing the first 5 columns
set.seed(1)
starwars[sample(nrow(starwars), 10), 1:5]

# A tibble: 10 x 5
   name            height  mass hair_color skin_color
                            
 1 Dexter Jettster    198 102   none       brown     
 2 Sebulba            112  40   none       grey, red 
 3 Luke Skywalker     172  77   blond      fair      
 4 Jar Jar Binks      196  66   none       orange    
 5 Bib Fortuna        180  NA   none       pale      
 6 Han Solo           180  80   brown      fair      
 7 Cliegg Lars        183  NA   brown      fair      
 8 Eeth Koth          171  NA   black      brown     
 9 Boba Fett          183  78.2 black      fair      
10 Yarael Poof        264  NA   none       white

• dplyr

library(tidyverse)
set.seed(1)
starwars %>%
  sample_n(10) %>%
  select(1:5)

• data.table

library(data.table)
set.seed(1)
data.table(starwars)[sample(.N, 10), 1:5]

Descriptive statistics by group in R

Problem

We'd like to report descriptive statistics in R by a grouping variable and subsetting the output statistics.

Solution

We will use the data frame iris, columns Sepal.Length and Sepal.Width and grouping by Species. In our example, we want to return the mean, the standard deviation, the skewness and kurtosis.

Subset of descriptive statistics by group

library(psych)
# Variables by index
d <- describeBy(iris[1:2], group = iris$Species)
# Two options to subset the statistics:
lapply(d, "[", , c(3, 4, 11, 12))
lapply(d, subset, , c(3, 4, 11, 12)) 

# Variables by name
i <- match(c("Sepal.Length", "Petal.Length"), names(iris))
d <- describeBy(iris[i], group = iris$Species)
lapply(d, subset, , c("mean", "sd", "skew", "kurtosis")) 

$setosa
             mean   sd skew kurtosis
Sepal.Length 5.01 0.35 0.11    -0.45
Sepal.Width  3.43 0.38 0.04     0.60

$versicolor
             mean   sd  skew kurtosis
Sepal.Length 5.94 0.52  0.10    -0.69
Sepal.Width  2.77 0.31 -0.34    -0.55

$virginica
             mean   sd skew kurtosis
Sepal.Length 6.59 0.64 0.11    -0.20
Sepal.Width  2.97 0.32 0.34     0.38

Subset of descriptive statistics without grouping

# Seleccionamos las columnas deseadas de la tabla
d <- describe(iris[1:2])
# Subsetting output statistics
d[, c(3, 4, 11, 12)]

             mean   sd skew kurtosis
Sepal.Length 5.84 0.83 0.31    -0.61
Sepal.Width  3.06 0.44 0.31     0.14

References

• Descriptive statistics in R

Descriptive statistics in R

Problem

We'd like to compute descriptive statistics in R.

Solution

The summary function returns a set of summary statistics for the input (a vector, data frame or model).

# For a variable
summary(iris$Sepal.Length)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  4.300   5.100   5.800   5.843   6.400   7.900 

# For a data frame
summary(iris)

  Sepal.Length    Sepal.Width     Petal.Length    Petal.Width          Species  
 Min.   :4.300   Min.   :2.000   Min.   :1.000   Min.   :0.100   setosa    :50  
 1st Qu.:5.100   1st Qu.:2.800   1st Qu.:1.600   1st Qu.:0.300   versicolor:50  
 Median :5.800   Median :3.000   Median :4.350   Median :1.300   virginica :50  
 Mean   :5.843   Mean   :3.057   Mean   :3.758   Mean   :1.199                  
 3rd Qu.:6.400   3rd Qu.:3.300   3rd Qu.:5.100   3rd Qu.:1.800                  
 Max.   :7.900   Max.   :4.400   Max.   :6.900   Max.   :2.500  

The fivenum function returns Tukey's five number summary (minimum, lower-hinge, median, upper-hinge, maximum) for the input data.

# For a variable
fivenum(iris$Sepal.Width)

[1] 2.0 2.8 3.0 3.3 4.4

The boxplot.stats function returns the statistics necessary for producing box plots.

boxplot.stats(iris$Sepal.Width)

$stats
[1] 2.2 2.8 3.0 3.3 4.0

$n
[1] 150

$conf
[1] 2.935497 3.064503

$out
[1] 4.4 4.1 4.2 2.0

To return a specific statistic we type boxplot.stats(iris$Sepal.Width) followed by:

$stats - vector with Tukey's five number summary.
$n - the number of non-NA observations.
$conf - the lower and upper extremes of the ‘notch’.
$out- outliers.

The psych package

For a data frame

library(psych)
describe(iris)

             vars   n mean   sd median trimmed  mad
Sepal.Length    1 150 5.84 0.83   5.80    5.81 1.04
Sepal.Width     2 150 3.06 0.44   3.00    3.04 0.44
Petal.Length    3 150 3.76 1.77   4.35    3.76 1.85
Petal.Width     4 150 1.20 0.76   1.30    1.18 1.04
Species*        5 150  NaN   NA     NA     NaN   NA
             min  max range  skew kurtosis   se
Sepal.Length 4.3  7.9   3.6  0.31    -0.61 0.07
Sepal.Width  2.0  4.4   2.4  0.31     0.14 0.04
Petal.Length 1.0  6.9   5.9 -0.27    -1.42 0.14
Petal.Width  0.1  2.5   2.4 -0.10    -1.36 0.06
Species*     Inf -Inf  -Inf    NA       NA   NA

Statistics by group

describeBy(iris, group = iris$Species)

group: setosa
             vars  n mean   sd median trimmed  mad
Sepal.Length    1 50 5.01 0.35    5.0    5.00 0.30
Sepal.Width     2 50 3.43 0.38    3.4    3.42 0.37
Petal.Length    3 50 1.46 0.17    1.5    1.46 0.15
Petal.Width     4 50 0.25 0.11    0.2    0.24 0.00
Species*        5 50  NaN   NA     NA     NaN   NA
             min  max range skew kurtosis   se
Sepal.Length 4.3  5.8   1.5 0.11    -0.45 0.05
Sepal.Width  2.3  4.4   2.1 0.04     0.60 0.05
Petal.Length 1.0  1.9   0.9 0.10     0.65 0.02
Petal.Width  0.1  0.6   0.5 1.18     1.26 0.01
Species*     Inf -Inf  -Inf   NA       NA   NA
--------------------------------------- 
group: versicolor
             vars  n mean   sd median trimmed  mad
Sepal.Length    1 50 5.94 0.52   5.90    5.94 0.52
Sepal.Width     2 50 2.77 0.31   2.80    2.78 0.30
Petal.Length    3 50 4.26 0.47   4.35    4.29 0.52
Petal.Width     4 50 1.33 0.20   1.30    1.32 0.22
Species*        5 50  NaN   NA     NA     NaN   NA
             min  max range  skew kurtosis   se
Sepal.Length 4.9  7.0   2.1  0.10    -0.69 0.07
Sepal.Width  2.0  3.4   1.4 -0.34    -0.55 0.04
Petal.Length 3.0  5.1   2.1 -0.57    -0.19 0.07
Petal.Width  1.0  1.8   0.8 -0.03    -0.59 0.03
Species*     Inf -Inf  -Inf    NA       NA   NA
--------------------------------------- 
group: virginica
             vars  n mean   sd median trimmed  mad
Sepal.Length    1 50 6.59 0.64   6.50    6.57 0.59
Sepal.Width     2 50 2.97 0.32   3.00    2.96 0.30
Petal.Length    3 50 5.55 0.55   5.55    5.51 0.67
Petal.Width     4 50 2.03 0.27   2.00    2.03 0.30
Species*        5 50  NaN   NA     NA     NaN   NA
             min  max range  skew kurtosis   se
Sepal.Length 4.9  7.9   3.0  0.11    -0.20 0.09
Sepal.Width  2.2  3.8   1.6  0.34     0.38 0.05
Petal.Length 4.5  6.9   2.4  0.52    -0.37 0.08
Petal.Width  1.4  2.5   1.1 -0.12    -0.75 0.04
Species*     Inf -Inf  -Inf    NA       NA   NA

Read compressed files in R using readr

Problem

We need to read compressed files in R.

Solution

We will use the package readr:

Files ending in .gz, .bz2, .xz, or .zip will be automatically uncompressed. Files starting with http://, https://, ftp://, or ftps:// will be automatically downloaded. Remote gz files can also be automatically downloaded and decompressed.

In our example we use the file title.ratings.tsv.gz.

library(readr)
df_ratings <- read_tsv('title.ratings.tsv.gz', na = "\\N", quote = '')
df_ratings %>% head()

We can provide the URL and it will be automatically downloaded and decompressed

df_ratings <- read_tsv('https://datasets.imdbws.com/title.ratings.tsv.gz', na = "\\N", quote = '')
df_ratings %>% head()

Results

# A tibble: 6 x 3
  tconst    averageRating numVotes
                   
1 tt0000001           5.8     1423
2 tt0000002           6.4      168
3 tt0000003           6.6     1016
4 tt0000004           6.4      100
5 tt0000005           6.2     1713
6 tt0000006           5.5       88

How to list the available data sets in R

Problem

We want to list the available data sets in R.

Solution

1. List all available data sets

data()

Data sets in package ‘datasets’:

AirPassengers           Monthly Airline Passenger Numbers 1949-1960
BJsales                 Sales Data with Leading Indicator
BJsales.lead (BJsales)
                        Sales Data with Leading Indicator
BOD                     Biochemical Oxygen Demand
CO2                     Carbon Dioxide Uptake in Grass Plants
ChickWeight             Weight versus age of chicks on different diets
DNase                   Elisa assay of DNase
EuStockMarkets          Daily Closing Prices of Major European Stock
                        Indices, 1991-1998
Formaldehyde            Determination of Formaldehyde
HairEyeColor            Hair and Eye Color of Statistics Students
...                    ...

2. List all available data sets, including those packages not currently loaded.

data(package = .packages(all.available = TRUE))

Data sets in package ‘aqp’:

amarillo                Amarillo Soils
ca630                   Soil Data from the Central Sierra Nevada Region
                        of California
munsell                 Munsell to sRGB Lookup Table for Common Soil
                        Colors
rruff.sample            Sample XRD Patterns
soil_minerals           Munsell Colors of Common Soil Minerals
sp1                     Soil Profile Data Example 1
sp2                     Honcut Creek Soil Profile Data
sp3                     Soil Profile Data Example 3
sp4                     Soil Chemical Data from Serpentinitic Soils of
                        California
sp5                     Sample Soil Database #5
sp6                     Soil Physical and Chemical Data from
                        Manganiferous Soils

Data sets in package ‘beeswarm’:

breast                  Lymph-node-negative primary breast tumors

3. List all available data sets of a specific package

data(package = "ISLR")

Data sets in package ‘ISLR’:

Auto                    Auto Data Set
Caravan                 The Insurance Company (TIC) Benchmark
Carseats                Sales of Child Car Seats
College                 U.S. News and World Report's College Data
Credit                  Credit Card Balance Data
Default                 Credit Card Default Data
Hitters                 Baseball Data
Khan                    Khan Gene Data
NCI60                   NCI 60 Data
OJ                      Orange Juice Data
Portfolio               Portfolio Data
Smarket                 S&P Stock Market Data
Wage                    Mid-Atlantic Wage Data
Weekly                  Weekly S&P Stock Market Data

References

• data{utils}