fabricatr: what changes, what does not, and why

nest = FALSE becomes declare_level()

In 1.x, an independent level that will later be crossed or linked is written as add_level(..., nest = FALSE). The flag is doing real semantic work, and the name of the function actively fights it: you are not adding a level to the hierarchy, you are declaring a separate population.

It still works, and tells you what to write:

fabricate(
  countries = add_level(N = 3, gdp = runif(N), nest = FALSE)
)
#> Warning: `nest = FALSE` is deprecated in fabricatr 2.0.
#> Write:  declare_level(N = 3, gdp = runif(N))
#> Not:    add_level(N = 3, gdp = runif(N), nest = FALSE)
#> This warning is displayed once per session.
#> # A tibble: 3 × 2
#>   countries   gdp
#>   <chr>     <dbl>
#> 1 1         0.941
#> 2 2         0.794
#> 3 3         0.965

Nesting is never a flag in 2.0: a second add_level() nests, nest_level() nests, and declare_level() stands alone.

join_using() becomes .by

1.x’s join_using() is a quoting helper that exists to pass level names into cross_levels() and link_levels(). 2.0 takes a character vector instead, and performs the crossing with the same semantics as dplyr::cross_join(). join_using() is kept, and by = is accepted alongside .by =, so the 1.x spelling runs and warns.

fabricate(
  countries = declare_level(N = 3, gdp = runif(N, 1, 10)),
  years     = declare_level(N = 2, shock = runif(N, 0, 1)),
  obs       = cross_levels(.by = c("countries", "years"), Y = gdp + shock)
)
#> # A tibble: 6 × 6
#>   countries   gdp years shock obs       Y
#>   <chr>     <dbl> <chr> <dbl> <chr> <dbl>
#> 1 1          8.66 1     0.853 1      9.51
#> 2 2          2.44 1     0.853 2      3.29
#> 3 3          5.37 1     0.853 3      6.22
#> 4 1          8.66 2     0.642 4      9.30
#> 5 2          2.44 2     0.642 5      3.08
#> 6 3          5.37 2     0.642 6      6.01

link_levels() samples N rows from that product rather than taking all of it, and rho is now an ordinary argument rather than something smuggled inside join_using(). The 1.x spelling, link_levels(N = 2000, by = join_using(primary, secondary, rho = 0.8)), still runs, carries its rho through, and warns with the call below:

students <- fabricate(
  primary   = declare_level(N = 20, p_quality = runif(N, 1, 10)),
  secondary = declare_level(N = 15, s_quality = runif(N, 1, 10)),
  students  = link_levels(N = 2000, .by = c("primary", "secondary"), rho = 0.8,
                          score = p_quality + s_quality + rnorm(N))
)
cor(as.numeric(students$primary), as.numeric(students$secondary),
    method = "spearman")
#> [1] 0.7685407

rho and sigma sit behind the dots in the signature, link_levels(N, .by, ..., rho = 0, sigma = NULL), so a column called r, s or si is a column and not a partial match for the correlation. The cost is that they have to be named: link_levels(N, .by, 0.5) no longer sets rho, and rho = 0.5 does.

One inherited quirk surprises people. rho is the rank correlation between the two levels’ ID assignments, not between the variables defined on those levels. Students at high-numbered primary schools are matched to high-numbered secondary schools; whether that also correlates p_quality with s_quality depends entirely on whether quality happens to rise with ID, and for a runif() draw it does not. 2.0 reproduces 1.x’s behavior here exactly, down to the same draws from the same seed when mvnfast is not installed. If you want quality-to-quality sorting, define the level variables as monotone in the level ID.

For more than two levels, pass a full correlation matrix as sigma instead of the scalar rho.

by becomes .by in modify_level()

The leading dot is the tidyverse convention for an argument that controls the verb rather than naming a column, and it removes the ambiguity that arises when a dataset happens to contain a column called by. modify_level() is otherwise dplyr::mutate(), with .by giving split-apply-combine over a grouping column:

fabricate(
  N = 12,
  cluster = rep(1:3, each = 4),
  Y = rnorm(N),
  modify_level(cluster_mean = mean(Y), .by = "cluster")
)
#> # A tibble: 12 × 4
#>    ID    cluster       Y cluster_mean
#>    <chr>   <int>   <dbl>        <dbl>
#>  1 01          1 -0.896        0.0984
#>  2 02          1  0.327        0.0984
#>  3 03          1  0.929        0.0984
#>  4 04          1  0.0338       0.0984
#>  5 05          2 -1.69        -0.238 
#>  6 06          2 -0.402       -0.238 
#>  7 07          2 -0.155       -0.238 
#>  8 08          2  1.30        -0.238 
#>  9 09          3 -1.59        -0.128 
#> 10 10          3  0.505       -0.128 
#> 11 11          3 -0.140       -0.128 
#> 12 12          3  0.713       -0.128

The old spelling groups correctly rather than quietly computing an ungrouped answer, and says what to write:

fabricate(
  N = 12, cluster = rep(1:3, each = 4), Y = rnorm(N),
  modify_level(cluster_mean = mean(Y), by = "cluster")
)
#> Warning: `by =` is deprecated in fabricatr 2.0.
#> Write:  modify_level(cluster_mean = mean(Y), .by = "cluster")
#> Not:    modify_level(cluster_mean = mean(Y), by = "cluster")
#> This warning is displayed once per session.
#> # A tibble: 12 × 4
#>    ID    cluster        Y cluster_mean
#>    <chr>   <int>    <dbl>        <dbl>
#>  1 01          1 -0.551         -0.325
#>  2 02          1  0.00535       -0.325
#>  3 03          1 -0.368         -0.325
#>  4 04          1 -0.386         -0.325
#>  5 05          2 -0.178         -0.258
#>  6 06          2 -0.728         -0.258
#>  7 07          2 -0.358         -0.258
#>  8 08          2  0.231         -0.258
#>  9 09          3  0.944         -0.498
#> 10 10          3 -1.44          -0.498
#> 11 11          3  0.273         -0.498
#> 12 12          3 -1.77          -0.498

N must be named

1.x’s signature is fabricate(..., data, N, ID_label), and a bare number in the dots, fabricate(100, Y = rnorm(N)), was read as N = 100. 2.0 does not read it that way: every expression in the dots is a column and needs a name, so the same call stops with

Every column needs a name. Expression 1 in fabricate(), `100`, has none.
  If this is the number of rows, write `N = 100`: N is supplied by name.

Write fabricate(N = 100, Y = rnorm(N)). The same rule catches an unnamed column anywhere, fabricate(N = 5, rnorm(N)), which 1.x failed on with an indexing error from deep inside and which an early build of 2.0 dropped without a word.

The result is a tibble

fabricate() returns a tibble where 1.0.2 returned a plain data frame. The columns, their order, their types, and their values are the same; only the class differs, from "data.frame" to c("tbl_df", "tbl", "data.frame"). A tibble is a data frame, so anything that takes one still works.

Two things do change for a caller. A test written expect_identical(fabricate(...), some_data_frame) compares the class attribute and now fails, which is how this surfaced: six of the eight test failures in a revdepcheck of DeclareDesign 1.1.1 are that comparison and not one of them concerns a value. And single-bracket indexing does not drop to a vector, so df[, "Y"] is a one-column tibble where 1.0.2 gave a numeric vector. df$Y and df[["Y"]] are unchanged.

modify_level() at a named level works at that level

Named after a level above the current one, modify_level() evaluates its expressions once per unit of that level: N is the number of regions, the columns in view are the ones built at the regions level or above, and each result is written back to every city in its region.

fabricate(
  regions = add_level(N = 3, a = 1:3),
  cities  = add_level(N = 2, b = rnorm(N)),
  regions = modify_level(z = rnorm(N), a2 = a * 2)
)
#> # A tibble: 6 × 6
#>   regions     a cities       b       z    a2
#>   <chr>   <int> <chr>    <dbl>   <dbl> <dbl>
#> 1 1           1 1      -0.767  -2.49       2
#> 2 1           1 2       0.516  -2.49       2
#> 3 2           2 3       0.314   0.688      4
#> 4 2           2 4       0.0860  0.688      4
#> 5 3           3 5      -1.66    0.0549     6
#> 6 3           3 6      -0.898   0.0549     6

1.x evaluates the same call the same way, three draws of z and a2 from the region’s own a, but then returns the regions frame alone: three rows, with the cities and their b gone. 2.0 keeps the six rows. An early build of 2.0 ignored the level name instead and drew one z per city, which is the silent kind of wrong, and is why this section exists.

A column built at a level nested inside is out of view at the outer level, as it is in 1.x, where the same call could see nothing below the level it named. regions = modify_level(bb = mean(b)) stops and says so, because b varies within a region. To summarise it per region, leave the level name off and group: modify_level(bb = mean(b), .by = "regions").

reveal_outcomes() keeps the type of the outcome

The one place where 2.0 deliberately declines to reproduce 1.x, because what 1.x did was wrong rather than merely different.

The switching equation picks, for each unit, the potential outcome column matching that unit’s assignment. 1.x did the picking by matrix-indexing the data frame of potential outcome columns, which goes through as.matrix(), and a matrix holds one type. So a factor outcome came back as character and a Date came back as a string. For an unordered factor that is an annoyance. For an ordered one it is a wrong answer that does not announce itself:

set.seed(1)
likert <- fabricate(
  N = 200, U = rnorm(N),
  potential_outcomes(Y ~ draw_ordered(0.6 * Z + U, breaks = c(-1, 1),
                                      labels = c("lo", "mid", "hi")))
)
revealed <- fabricate(data = likert, Z = rbinom(200, 1, 0.5),
                      Yobs = reveal_outcomes(Y ~ Z))
levels(revealed$Yobs)
#> [1] "lo"  "mid" "hi"

Those are the levels the outcome was declared with. Under 1.x, Yobs was a character vector, and the levels were gone: the natural next step, making it a factor to model it, put them back in alphabetical order as hi < lo < mid. A regression on that has the wrong baseline and, if it is treated as ordered, fits polynomial contrasts against an ordering nobody chose. Nothing warns, because nothing is malformed; the numbers are simply answers to a different question.

2.0 selects column by column instead, so the revealed vector keeps whatever the potential outcomes carried: a factor stays a factor with its levels and their order, a Date stays a Date, and numeric, integer, logical, and character are as they always were.

One case has no right answer. If the conditions disagree about their levels, there is no set to keep, and inventing an order across level sets that were never meant to be compared is worse than declining. Those reveal as character, with a warning naming the outcome and saying to give every condition the same levels.

This matters beyond fabricatr, because reveal_outcomes() is re-exported by DeclareDesign and declare_measurement(Y = reveal_outcomes(Y ~ Z)) is the idiom the book runs on. Any design with a Likert or other ordered outcome was measuring it through this.

ID columns are zero-padded character, as in 1.x

IDs are zero-padded character strings, so a design with twelve clusters gets "01" through "12" and the padding keeps character sort order matching numeric order:

fabricate(
  clusters = add_level(N = 12),
  units    = add_level(N = 2)
)
#> # A tibble: 24 × 2
#>    clusters units
#>    <chr>    <chr>
#>  1 01       01   
#>  2 01       02   
#>  3 02       03   
#>  4 02       04   
#>  5 03       05   
#>  6 03       06   
#>  7 04       07   
#>  8 04       08   
#>  9 05       09   
#> 10 05       10   
#> # ℹ 14 more rows

The type is the point. An ID is a label, not a quantity, and an integer ID is a number that R will happily use as one. lm(Y ~ clusters) on integer cluster IDs silently fits a linear trend in cluster number where you meant fixed effects, and mean(clusters) returns a value rather than an error. Character IDs make both mistakes impossible.

Note that the padding width is set per level, by the number of units at that level rather than by anything global. A frame of 4 clusters holding 3 units each carries clusters running "1" to "4" alongside units running "01" to "12", as in the ICC example above. 1.x does the same.

Levels name their own ID columns after the level, exactly as in 1.x. Flat fabrication gets a column called ID, which ID_label renames and ID_label = NA suppresses:

fabricate(N = 3, Y = rnorm(N), ID_label = "respondent")
#> # A tibble: 3 × 2
#>   respondent      Y
#>   <chr>       <dbl>
#> 1 1           0.894
#> 2 2          -1.05 
#> 3 3           1.97
fabricate(N = 3, Y = rnorm(N), ID_label = NA)
#> # A tibble: 3 × 1
#>        Y
#>    <dbl>
#> 1 -0.384
#> 2  1.65 
#> 3  1.51

total_sd is a parameter, not a rescaling

This is the one change that gives different numbers from the same call, so it gets the argument on both sides rather than an assertion.

draw_normal_icc() builds each value from a cluster-level draw and an individual-level draw, so there are four ways to describe the same pair of variances: the ICC, the within-cluster standard deviation, the between-cluster standard deviation, and the total. They are related by \(\sigma_{total}^2 = \sigma_{between}^2 + \sigma_{within}^2\) and \(ICC = \sigma_{between}^2 / \sigma_{total}^2\), so any two of them pin the other two.

1.x treats total_sd differently from the other three. It draws using the ICC, then rescales the finished vector so its sample standard deviation equals total_sd. 2.0 solves for the two standard deviations instead and draws once: with \(\sigma_{total}\) set by total_sd, \(\sigma_{between} = \sigma_{total}\sqrt{ICC}\) and \(\sigma_{within} = \sigma_{total}\sqrt{1 - ICC}\).

The consequence, over 200 draws of 100 observations with ICC = 0.4 and total_sd = 2:

The case for changing it. Rescaling makes the realised standard deviation exactly total_sd in every draw. A simulation that fabricates 10,000 datasets gets 10,000 datasets whose outcome variance is identical to the fourth decimal place, which no real sampling process would produce. Any diagnosand that depends on the spread of the outcome, and any standard error with an estimated variance in it, is then computed against a DGP with one of its sources of variation switched off. For a package whose entire purpose is repeated sampling, that is the wrong default. The fix also costs nothing in expressiveness, since the closed form is exact and defined at every ICC including 0 and 1.

The case against. It is a divergence, and this vignette otherwise promises that identical code gives identical results. Anyone who has calibrated a design against 1.x’s total_sd will see their numbers move. There is also a legitimate reading of total_sd under which the old behavior is correct rather than buggy: if you want a dataset whose outcome has exactly the standard deviation you named, perhaps because you are matching a published descriptive statistic, rescaling gives you that and parameterising does not. Under that reading 2.0 has taken away a feature.

Why the change anyway. The two readings are not equally served by the two implementations. If you want the parameter, you cannot recover it from the rescaled version. If you want the exact sample standard deviation, you can still have it in one line:

y <- draw_normal_icc(clusters = rep(1:10, each = 10), ICC = 0.4)
y <- (y - mean(y)) * 2 / sd(y) + mean(y)
sd(y)
#> [1] 2

One correction to the upstream discussion, because it changes what is actually at stake. The thread worries that rescaling “will possibly distort the between group differences.” It does not. Rescaling is an affine transformation and the ICC is scale-invariant, so the realised ICC is untouched; the mean R² is the same either way. The only thing rescaling changes is the sampling variability of the variance. That is a narrower objection than the thread supposes, and still a sufficient one.

A side benefit of solving rather than rescaling is that total_sd can now be paired with any of the other three. 1.x accepts total_sd only alongside ICC, because rescaling needs an ICC to draw from first. In 2.0 any two of the ICC, sd (within), sd_between, and total_sd determine the draw, so total_sd with sd_between works:

cl <- rep(1:20, each = 10)
sd(draw_normal_icc(clusters = cl, total_sd = 3, sd_between = 2))
#> [1] 3.350732

The correlated draw does not depend on your machine

link_levels(rho = ) and link_levels(sigma = ) correlate row assignments through a Gaussian copula, which needs a draw from a multivariate normal. 1.x takes that draw one of two ways:

if (requireNamespace("mvnfast")) {
  correlated_sn <- mvnfast::rmvn(N, ncores = getOption("mc.cores", 2L), mu, sigma)
} else {
  right_chol <- chol(sigma, pivot = TRUE)
  correlated_sn <- matrix(rnorm(N * ndim), nrow = N) %*% right_chol
}

The two paths consume the random number stream differently, so the same seed gives different data depending on whether an optional package is installed, and the mvnfast branch depends on your core count as well. A simulation package should not have that property, so 2.0 always takes the second path and does not suggest mvnfast at all.

The cost is nothing you can feel. mvnfast earns its name in hundreds of dimensions; link_levels() has one dimension per level, so two or three. Measured at three levels, it saves 19 milliseconds on a million rows, and nothing measurable at the sizes anyone actually links.

The consequence is that if you have mvnfast installed, a correlated link_levels() gives different numbers here than in 1.x. Without it, the two agree exactly. Either way, 2.0 gives the same numbers as itself on every machine, which 1.x cannot promise.

On reproducing an old seed

fabricatr does not promise that a given seed produces the same fabricated data in 2.0 as it did in 1.x. Two things in this release land at a different point in the random number stream than 1.x did, and neither is treated here as a break.

The reason is what the package is for. fabricatr invents data to stand in for data you do not have yet. A draw is a draw from the distribution you specified, and which particular draw you got is not a finding: re-running with a different seed is a thing you are supposed to be able to do, and a design whose conclusion depends on the seed is a design in trouble. Compare randomizr, where the draw is the result. A real assignment of real units to real treatments is a fact about the study, recorded in its data and defensible years later, so randomizr’s stream is part of its contract and does not move across versions. fabricatr’s is not, and does.

What 2.0 does promise is the part that matters: the distribution is the one the arguments describe, and 2.0 gives the same numbers as itself on every machine, which 1.x could not say, because its correlated draw ran through mvnfast when that package happened to be installed and took a core count as an argument.

So if you have a saved analysis whose numbers came out of fabricate() under a stored seed, expect them to move, and re-run and re-store rather than hunting for the difference. The two places are:

A correlated link_levels(), if you have mvnfast installed. Covered in the section above. Without that package, 1.x and 2.0 agree exactly, verified 2,000 rows deep against a library built without it.

Any call that starts from data. fabricatr 1.x calls an internal helper, makeUU(), to invent a unique name for the workspace slot it is about to put your data frame in. That helper is sprintf("%X-%X", as.integer(Sys.time()), sample.int(.Machine$integer.max, 1)), and sample.int() consumes exactly two uniforms. So every fabricate(data = df, ...) in 1.x began two draws into the stream. 2.0 has no temporary to name, because a column expression is evaluated with one rlang::eval_tidy() against a plain named list, so it starts at the beginning:

set.seed(7)
reference <- runif(6)
set.seed(7)
fabricate(data = data.frame(u = 1:6), Z = runif(6))$Z
#> [1] 0.98890930 0.39774545 0.11569778 0.06974868 0.24374939 0.79201043
reference
#> [1] 0.98890930 0.39774545 0.11569778 0.06974868 0.24374939 0.79201043

In 1.0.2 that same call returns reference[3:8]. The offset is exactly two whatever the data frame’s size, whether or not it carries an ID column, and however many columns you then build. It is pinned by a test, not because the stream position is promised, but because a silent change in it is worth knowing about deliberately rather than discovering downstream.

N is never a column

N is injected into the evaluation mask for every column expression, so rnorm(N) works everywhere it did before, but it is not written into the output. Ask for it explicitly if you want it:

fabricate(N = 3, n_obs = N, Y = rnorm(N))
#> # A tibble: 3 × 3
#>   ID    n_obs      Y
#>   <chr> <int>  <dbl>
#> 1 1         3 -0.412
#> 2 2         3 -0.971
#> 3 3         3 -0.947

Beyond tidiness, keeping N out of the frame is what allows cross_levels() to work without a name collision between the two levels being crossed.

Existing data frames come in as levels

import_level() makes a data frame a level of the design, which is what lets two real data sets be crossed or linked. fabricate(data = df, ...) starts from one frame and builds down from it, and that is still the way to add columns to data you already have. What it cannot do is bring in a second frame, because a level is the thing cross_levels() and link_levels() work on.

The level’s ID is the frame’s own key. A column named after the level is used as it stands, the way a join matches on a shared name, and .id names the column when it is called something else:

individuals <- data.frame(individuals = c("ann", "bob", "cyd"),
                          ind_shock = c(-0.4, 0.1, 0.8))
periods <- data.frame(year = 2020:2022, period_shock = c(0.2, -0.1, 0.3))

fabricate(
  individuals = import_level(individuals),
  period      = import_level(periods, .id = "year"),
  obs         = cross_levels(.by = c("individuals", "period"),
                             Y = ind_shock + period_shock + rnorm(N))
)
#> # A tibble: 9 × 6
#>   individuals ind_shock period period_shock obs       Y
#>   <chr>           <dbl>  <int>        <dbl> <chr> <dbl>
#> 1 ann              -0.4   2020          0.2 1     0.548
#> 2 bob               0.1   2020          0.2 2     0.183
#> 3 cyd               0.8   2020          0.2 3     1.15 
#> 4 ann              -0.4   2021         -0.1 4     1.69 
#> 5 bob               0.1   2021         -0.1 5     0.357
#> 6 cyd               0.8   2021         -0.1 6     3.42 
#> 7 ann              -0.4   2022          0.3 7     2.18 
#> 8 bob               0.1   2022          0.3 8     0.724
#> 9 cyd               0.8   2022          0.3 9     3.00

.id says which incoming column is the key. The key then takes the level’s name, so year arrives as period. Its type and values are kept, where a fabricated key is a zero-padded character string, so period still holds 2020 to 2022 and joins back to periods with by = c(period = "year"). The one column name an imported frame cannot carry is N, which is the row count in every expression a level evaluates, and which is refused rather than dropped in silence. A frame with no key of its own gets numbered rows, as add_level() numbers them. A key with a repeated or a missing value is refused, because a level’s rows are its units.

Levels nest under an imported one as they nest under any other, so a real sample of villages can carry simulated citizens:

fabricate(
  villages = import_level(data.frame(villages = c("v1", "v2"),
                                     v_income = c(10, 12))),
  citizens = nest_level(N = 3, income = v_income + rnorm(N))
)
#> # A tibble: 6 × 4
#>   villages v_income citizens income
#>   <chr>       <dbl> <chr>     <dbl>
#> 1 v1             10 1         10.5 
#> 2 v1             10 2          9.11
#> 3 v1             10 3          9.69
#> 4 v2             12 4         12.0 
#> 5 v2             12 5         13.0 
#> 6 v2             12 6         12.8

n() counts the rows in hand

n() is the number of rows the level is building, and it works in every declaration: flat, add_level(), nest_level(), cross_levels(), and modify_level(), where it counts the level being worked at. N names the same quantity and still does. The reason to have both is that N is a name: a design that also carries an N in its workspace leaves a reader unsure which one an expression reads, and redesign() rebinds names, where nothing can rebind a call. A data set that carries a column called n keeps it, because R skips non-function bindings when it looks up a call, so sum(n) reads the column and n() reads the count. 1.0.2 fails with could not find function "n".

fabricate(
  villages = add_level(N = 3, wealth = rnorm(n())),
  citizens = nest_level(N = 2, income = wealth + rnorm(n())),
  villages = modify_level(n_villages = n())
)
#> # A tibble: 6 × 5
#>   villages wealth citizens income n_villages
#>   <chr>     <dbl> <chr>     <dbl>      <int>
#> 1 1         0.705 1         1.98           3
#> 2 1         0.705 2         0.890          3
#> 3 2         1.31  3         2.06           3
#> 4 2         1.31  4         1.90           3
#> 5 3        -1.39  5        -2.37           3
#> 6 3        -1.39  6        -1.66           3

Level calls and ordinary columns can be mixed in a single fabricate() call. 1.x requires that the arguments be either all level calls or none of them, so the following raises “Arguments passed to ... must either all be calls to create or modify levels, or else none of them must be”:

fabricate(
  N = 12,
  cluster = rep(1:3, each = 4),
  Y = rnorm(N),
  modify_level(cluster_sd = sd(Y), .by = "cluster")
)
#> # A tibble: 12 × 4
#>    ID    cluster       Y cluster_sd
#>    <chr>   <int>   <dbl>      <dbl>
#>  1 01          1 -0.871       0.656
#>  2 02          1  0.719       0.656
#>  3 03          1  0.111       0.656
#>  4 04          1 -0.0785      0.656
#>  5 05          2 -0.420       0.897
#>  6 06          2 -0.562       0.897
#>  7 07          2  0.998       0.897
#>  8 08          2 -1.11        0.897
#>  9 09          3 -0.142       0.811
#> 10 10          3  0.315       0.811
#> 11 11          3  1.22        0.811
#> 12 12          3 -0.699       0.811

draw_likert() accepts a manual breaks vector. The same call returns all NA under fabricatr 1.0.2:

fabricate(N = 5, x = rnorm(N), l = draw_likert(x, breaks = c(-1, 0, 1)))
#> # A tibble: 5 × 3
#>   ID         x     l
#>   <chr>  <dbl> <int>
#> 1 1     -0.285     2
#> 2 2     -1.31      1
#> 3 3     -0.391     2
#> 4 4     -0.402     2
#> 5 5      1.35      4

draw_count() takes a log link, so a Poisson rate can be built from a linear predictor the way a probability is: draw_count(latent = x, link = "log") draws at exp(x). 1.0.2 refuses this outright, answering “Count data does not accept link functions”, which left no way to write a count outcome on the scale a model would fit it on. draw_binary() and draw_binomial() also accept "logistic" as a synonym for "logit", which 1.x does not.

fabricate(N = 6, x = rnorm(N), visits = draw_count(latent = x, link = "log"))
#> # A tibble: 6 × 3
#>   ID           x visits
#>   <chr>    <dbl>  <int>
#> 1 1      0.591        4
#> 2 2      0.101        0
#> 3 3      0.931        4
#> 4 4     -0.263        0
#> 5 5     -0.00767      1
#> 6 6      0.367        1

draw_count() also takes dispersion, for counts that vary more than a Poisson allows, which is most counts outside a textbook. The default 0 is the Poisson draw as before, on the same random number stream. A positive value draws from the negative binomial with the same mean and variance mean + dispersion * mean^2, so dispersion = 0.5 at a mean of 3 gives a variance of 7.5 rather than 3. The parameter is Stata’s alpha and the reciprocal of the theta that MASS::glm.nb() reports, so a design can be simulated at the dispersion a fitted model found.

counts <- fabricate(
  N = 10000,
  poisson = draw_count(mean = 3, N = N),
  overdispersed = draw_count(mean = 3, N = N, dispersion = 0.5)
)
sapply(counts[c("poisson", "overdispersed")], function(y) round(c(mean = mean(y), var = var(y)), 2))
#>      poisson overdispersed
#> mean    3.01          3.01
#> var     3.02          7.59

And correlate() works with base R random number generators, not only with the draw_* family, by mapping each r* function to its q* counterpart:

x <- rnorm(200)
y <- correlate(rnorm, mean = 0, sd = 1, given = x, rho = -0.5)
cor(x, y, method = "spearman")
#> [1] -0.4719388

Panels with persistent shocks

draw_normal_ar() draws a first-order autoregressive process inside each unit of a panel: each period’s value is rho times the last one plus a fresh shock. It is the error structure behind the problem Bertrand, Duflo, and Mullainathan (2004) found in difference-in-differences, where shocks that persist within a unit make standard errors that ignore the persistence far too small, and clustering by unit repairs them. With the draw in hand, that failure can be diagnosed in a design before it happens in an analysis.

panel <- fabricate(
  units = declare_level(N = 100, unit_effect = rnorm(N)),
  periods = declare_level(N = 10),
  obs = cross_levels(
    .by = c("units", "periods"),
    e = draw_normal_ar(clusters = units, time = periods, rho = 0.8),
    Y = unit_effect + e,
    voted = draw_binary(latent = Y, link = "probit"),
    contacts = draw_count(mean = 2, quantile_y = pnorm(e))
  )
)
head(panel)
#> # A tibble: 6 × 8
#>   units unit_effect periods obs        e      Y voted contacts
#>   <chr>       <dbl> <chr>   <chr>  <dbl>  <dbl> <int>    <dbl>
#> 1 001         1.13  01      0001   0.722  1.86      1        3
#> 2 002        -1.77  01      0002   1.52  -0.249     1        4
#> 3 003         1.28  01      0003  -0.746  0.532     0        1
#> 4 004        -0.851 01      0004  -0.624 -1.48      0        1
#> 5 005        -1.96  01      0005   1.50  -0.464     0        4
#> 6 006         1.44  01      0006  -0.156  1.28      0        2

Four choices in the draw are worth knowing.

• Each unit’s series starts at its stationary distribution, so the first period has the same variance as the tenth and there is no burn-in to discard.

• sd is the standard deviation of the process in any one period, not of each period’s new shock. Changing rho therefore changes how persistent the shocks are and leaves the variance where it was, so designs compared across rho differ in persistence alone. stats::arima.sim() reads its sd the other way, and at rho = 0.8 the same number there gives a process 1.7 times wider.

• The correlation between two observations of a unit is rho raised to the gap in time between them. A panel with a missing wave, or waves at uneven intervals, is the same process seen at fewer points and needs no special handling. time may be a number, a Date (read in days), or the character ID that cross_levels() gives a level.

• The draw is a latent normal variable. A binary or count outcome with the same persistence comes from passing it on, as latent to draw_binary() or as quantile_y = pnorm(e) to draw_count(), as the example does. The serial correlation of the outcome is then weaker than rho, as it is with correlate().

The function draws persistent errors, not state dependence. A model in which this period’s outcome depends on last period’s outcome itself, a lagged dependent variable, is a different model, and fitting one with unit fixed effects in a short panel has a bias of its own.