No Data, No Worries : Fabricatr

Neil Shephard

n.shephard@sheffield.ac.uk

Research Software Engineer, Department of Computer Science, University of Sheffield

Scan This

ns-rse.github.io/sheffieldr-fabricatr-20240621

Motivation

• Rarely work on my own data.
• Literate programming, reproducible workflow
  • Org-Babel
  • RMarkdown
  • Quarto.
• Scripts to clean, summarise, plot and analyse data.

The Turing Way project illustration by Scriberia. Used under a CC-BY 4.0 licence. DOI: 10.5281/zenodo.3332807

I’ve been a statistician, data scientist and more recently a Research Software Engineer for a number of years now and I rarely work on my own data. I collaborate with other researchers and either analyse their data for them or help guide them through the process of undertaking the anayyses themselves. My approach to doing this is to use a literate programming approach, a term introduced by computer scientist Donald Knuth, which has seend wide adoption in statistical and data science circles whereby text and code are inter-weaved into essentially scripts that can be re-run to clean, summarise plot and analys the data.

But there is a common problem that I encounter

Problem

Rarely is the final data set available!

“To consult the statistician after an experiment is finished is often merely to ask him to conduct a post mortem examination. He can perhaps say what the experiment died of.” - R.A. Fisher

Rarely is the full data set available, which can make it challenging to get started on developing the reproducible pipeline for analysing thet data.

If I’m really lucky people will be in contact before they’ve started their research and we can collaborate on what data should be collected and avoid the post-mortem scenario RA Fisher refers to in this well known quote.

Typically though that isn’t the case they have some data that is partially collected and they want to start interrogating it.

Solution

• Database schema/dictionary?

• Simulate data!

But all is not lost, if the structure of just a subset of the data is known, for example if there is a database schema or data dictionary for the researchers data then we can get started on developing our workflow by simulating some data and developing our scripts and literate programming workflow with the sample data.

This doesn’t help with what is often the bulk of the work of any analysis which is cleaning and tidying the data but it gives us something to work with for the tables, graphs and statistical analyses we might want to undertake.

Fabricatr

• fabricatr - Imagine Your Data Before You Collect It

#| label: setup
#| warning: false
#| eval: true
#| echo: true
install.packages(c("fabricatr", "tidyverse"))
library(correlate)
library(dplyr)
library(fabricatr)
library(ggplot2)
library(tidyr)

set.seed(55138)

• Binary variables
• Categorical/factors
• Continuous
• Time-Series
• Correlated
• Intra-Class Correlation (ICC)

The fabricatr package is a relatively new package that makes it easy to simulate lots of different types of variables such as binary, categorical, continuous or count data and often with various different structures such as time-series, correlation between variables and even nested study designs with intra-class correlation.

As the tag-line on the site reads you can Imagine your data before you collect it!

I’ve setup some basic examples to show how this package works and will run through some scenarios.

Basic Example

• 100 observations.
• Binary : draw_binary().
• Binomial : draw_binomial().
• Continuous (0-1) : runif().
• Random Normal (mean = 38, sd = 16) : rnorm().

dummy_df <- fabricatr::fabricate(
    N = 100,
    binary = fabricatr::draw_binary(
        N,
        prob = runif(N)
        ),
    binomial = fabricatr::draw_binomial(
        prob = runif(N),
        trials = 2
        ),
    uniform = runif(N),
    random_normal = rnorm(N, 38, 16)
)
dummy_df |> head()

   ID binary binomial     uniform random_normal
1 001      0        2 0.241211455      14.98915
2 002      0        2 0.002083675      28.97463
3 003      0        2 0.102645113      59.94338
4 004      1        2 0.124053932      39.85355
5 005      1        0 0.863406008      45.04457
6 006      1        0 0.655067994      59.41866

The main funciton is fabricate() which needs to know the number of observations you wish to simulate in the argument N and I’ve opted to keep the dataset small here and just have 100 observations.

We create a binary variable using the draw_binary() function which takes the value N and we define the probability function we wish to use to draw samples from, in this case a uniform distribution using R’s base function runif().

We then create a binomial variable with two trials to give us what is essentially count data using the draw_binomial() function, which needs a probability distribution from which to make draws and again we use the uniform distribution but it needs to know the number of trials we wish to simulate and here we have specified 2 so the resulting values of binomial will be 0, 1 or 2.

Next we create a uniform variable sampling directly from the base R function runif().

Finall we create a random_normal variable by sampling from a normal distribution with mean of 38 and a standard deviation of 16.

Basic Example - Summary

Characteristic	N = 1001
binary	41 (41%)
binomial
0	38 (38%)
1	30 (30%)
2	32 (32%)
uniform	0.38 (0.16, 0.70)
random_normal	38 (25, 49)
1 n (%); Median (IQR)

Figure 1

Figure 2

Figure 3

Categorical Variables

• 100 observations
• Random Normal (mean = 0, sd = 1) : rnorm()
• Categorise using draw_ordered() and setting breaks.
• Label variables with break_labels

cat_df <- fabricatr::fabricate(
    N = 100,
    x = rnorm(N),
    ordered = fabricatr::draw_ordered(
        x,
        breaks = c(-Inf,
                   -1,
                   -0.5,
                   0,
                   0.5,
                   1,
                   Inf),
        break_labels = c("Group 1",
                         "Group 2",
                         "Group 3",
                         "Group 4",
                         "Group 5",
                         "Group 6")
    )
)

Fabricatr also supports generating ordered categorical variables as you might get from surveys and we can generate these using the draw_ordered() function. First we generate a random normal variable x using rnorm() for our sample size N of 100 and this is passed to the draw_ordered() function where we specify a series of break points and labels for each of the categories that have been formed.

Categorical Variables - Summary

Characteristic	N = 1001
ordered
Group 1	14 (14%)
Group 2	13 (13%)
Group 3	19 (19%)
Group 4	22 (22%)
Group 5	13 (13%)
Group 6	19 (19%)
1 n (%)

Categorical variable with six levels.

Figure 4

We can tabulate this data using the excellent gtsummary package and also plot the data to see how many there are in each category.

Correlated Continuous Variables

EXPERIMENTAL

• 100 observations
• y (mean = 38, sd = 16) : using rnorm()
• a (mean = 78, sd = 14) : rnorm() rank correlation with y of 0.8 using correlate
• b (mean = 46, sd = 18) : rnorm() rank correlation with y of -0.4 using correlate

corr_df <- fabricatr::fabricate(
    N = 100,
    y = rnorm(N, mean = 64, sd = 11),
    a = fabricatr::correlate(given = y,
                             rho = 0.8,
                             rnorm,
                             mean = 78,
                             sd = 14),
    b = fabricatr::correlate(given = y,
                             rho = -0.4,
                             rnorm,
                             mean = 46,
                             sd = 18)
)

It is also possible to generate correlated variables here we again use a sample size of 100 and simulate y from a random normal distribution with a mean of 64 and a standard deviation of 11.

We then use the correlate() function to simulate a and b which are conditional on y. A has a rank correlation coefficient of 0.8 and samples from a normal distribution with a mean of 78 and standard deviation of 14, whilst b has a negative correlation with y of -0.4 and is a continuous variable drawn from a normal distribution with mean of 46 and standard dewviation of 18. And obviously b will be negatively correlated with a, albeit with a lower correlation.

Correlated Continuous Variables - Summary

Here I’ve made some simple scatter plots of the continuous variables to show how the correlation relates. Note that it is not perfect because of the randomness involved in sample…

Correlated Continuous Variables - Summary

Characteristic	N = 1001
y	64 (57, 70)
a	78 (66, 86)
b	47 (38, 59)
1 Median (IQR)

…but as we can see from the correlation coefficients… between the variables is roughly what we specified.

Correlated Discrete Variables

EXPERIMENTAL

• 100 observations
• Binomial q1 (prob = 0.38, trials = 10) : using draw_binomial()
• Binomial q1 (prob = 0.63, trials = 5) : draw_binomial() rank correlation with q1 of 0.74 using correlate
• Binomial q3 (prob = 0.79, trials = 5) : draw_binomial() rank correlation with q2 of -0.89 using correlate

corr_discrete_df <- fabricatr::fabricate(
    N = 100,
    q1 = fabricatr::draw_binomial(prob = 0.38,
                                  trials = 10,
                                  N = N),
    q2 = fabricatr::correlate(given = q1,
                             rho = 0.74,
                             fabricatr::draw_binomial,
                             prob = 0.63,
                             trials = 10),
    q3 = fabricatr::correlate(given = q2,
                             rho = -0.89,
                             fabricatr::draw_binomial,
                             prob = 0.79,
                             trials = 10)
)

In a similar vein we can correlate the discrete variables and here we draw count data using draw_binomial and specify a positive correlation of 0.74 between q1 and q2 and a negativecorrelation of -0.89 bbetween q2 and q3.

Correlated Discrete Variables - Summary

Characteristic	N = 1001
q1	4 (3, 5)
q2
2	1 (1.0%)
3	2 (2.0%)
4	7 (7.0%)
5	21 (21%)
6	23 (23%)
7	33 (33%)
8	10 (10%)
9	1 (1.0%)
10	2 (2.0%)
q3
5	3 (3.0%)
6	9 (9.0%)
7	25 (25%)
8	27 (27%)
9	27 (27%)
10	9 (9.0%)
1 Median (IQR); n (%)

Tabulation of q1 v q2

	2	3	4	5	6	7	8	9	10
0	0	0	2	0	0	0	0	0	0
1	0	0	3	1	1	0	0	0	0
2	1	1	1	3	5	1	0	0	0
3	0	1	1	10	8	5	0	0	0
4	0	0	0	5	5	10	2	0	0
5	0	0	0	2	1	11	4	0	0
6	0	0	0	0	2	3	1	1	0
7	0	0	0	0	1	1	3	0	0
8	0	0	0	0	0	2	0	0	1
9	0	0	0	0	0	0	0	0	1

Tabulation of q2 v q3

	5	6	7	8	9	10
2	0	0	0	0	0	1
3	0	0	0	0	1	1
4	0	0	0	0	3	4
5	0	0	0	5	13	3
6	0	0	4	10	9	0
7	0	5	15	12	1	0
8	2	2	6	0	0	0
9	0	1	0	0	0	0
10	1	1	0	0	0	0

We can summarise the data, again using the excellent gtsummary package and if we tabulate q1 and q2 we can see the positive correlation as higher counts of q1 correspond to higher counts of q2.

For q2 and q3 we can see that the negative correlation we specified is borne out in the simulated data with high values of q2 (the rows) corresponding to lower values of q3 the columns.

Correlated Discrete Variables - Summary

…and as with the continuous variables we can calculate the correlation coefficients and plot them as a heatmap. The values aren’t exactly what we specified because of random sampling but the pattern is there.

Intra Class Correlation (ICC)

schools_data <- fabricatr::fabricate(
  primary_schools = fabricatr::add_level(N = 20,
                      ps_quality = runif(N, 1, 10)),
  secondary_schools = fabricatr::add_level(N = 15,
                        ss_quality = runif(N, 1, 10),
                        nest = FALSE),
  students = fabricatr::link_levels(N = 1500,
               by = fabricatr::join_using(primary_schools,
                       secondary_schools),
               SAT_score = 800 + 13 * ps_quality + 26 * ss_quality +
                           rnorm(N, 0, 50)
             )
)

• Nested study design (e.g. schools, counties)
• Lifted straight from documentation
• Students (n 1500) at primary (ps_quality; n = 20) and secondary (ss_quality; n = 15) school.
• SAT_score is baseline (800) plus an additive effect for ps_quality and ss_quality and an additional random factor.

Intra Class Correlation (ICC) - Summary

lm(SAT_score ~ ps_quality + ss_quality, data = schools_data)

Call:
lm(formula = SAT_score ~ ps_quality + ss_quality, data = schools_data)

Coefficients:
(Intercept)   ps_quality   ss_quality  
     800.03        12.89        25.82  

If we then run a simple linear regression to predict the SAT_score from the quality of the primary and secondary school we see that intercept is indeed 800 and the coefficients for the effect of primary school is 13 and secondary 26 as was specified.

Intra Class Correlation (ICC)

schools_data <- fabricatr::fabricate(
  primary_schools = fabricatr::add_level(N = 20,
                      ps_quality = runif(N, 1, 10)),
  secondary_schools = fabricatr::add_level(N = 15,
                        ss_quality = runif(N, 1, 10),
                        nest = FALSE),
  students = fabricatr::link_levels(N = 1500,
               by = fabricatr::join_using(primary_schools,
                               secondary_schools,
                               rho = 0.5),
               SAT_score = 800 + 13 * ps_quality + 26 * ss_quality +
                           rnorm(N, 0, 50)
             )
)

We can introduce correlation between the quality of the primary and secondary school by specifying a value for rho (the spearkman rank correlation) in the join_using() call.

Time Series

• Time series data with fabricatr
• Can use the forecast package to generate ARIMA data.

Time series data can also be simulated, but I’ve not had time to prepare an example of this so I refer you to the authors article. Because you can use any other data creation package with fabricatr it is possible to leverage the forecast package to generate Auto-Regressive Integrated Moving Average data.

Summary

• fabricatr powerful but easy to use.
• Simulate data and setup analytical workflow (tables/graphs/models).
• Not useful for cleaning/tidying real data (~80% of work!).
• Need good database schema/data dictionary

Links

• fabricatr
  • Getting started with fabricatr
  • Building and Importing Data
  • Common Social Science variables
  • Panel and Cross-classified data
  • Time series data with fabricatr
  • Time series data with fabricatr
• gtsummary

Quarto RevealJS Template