Simulating verbal fluency

Simulating verbal fluency

Verbal fluency (say as many animals as you can in a minute) produces an overdispersed count. The variance grows faster than the mean, because words come in clusters rather than as independent events. sim_fluency() generates that structure from a negative binomial, with the overdispersion pinned to a variance-to-mean ratio you choose.

sim_fluency()

sim_fluency(n, effect = 0, mu = 15, vmr = 2, zi = 0)

Argument	What it is
n	Sample size per group (the frame has 2 * n rows).
effect	Group effect as a log rate ratio. log(1.5) is a 50% higher rate in the treatment group; 0 is the null.
mu	Control-group mean count.
vmr	Variance-to-mean ratio at the control mean. 1 is Poisson; 2, 4, 8 are increasingly overdispersed.
zi	Proportion of structural zeros to inject, for a floor-heavy impaired sample. 0 for none.

It returns a data frame with group (a factor, ctrl and treat) and y (the count).

library(countkit)
set.seed(1)
d <- sim_fluency(n = 80, effect = log(1.5), mu = 15, vmr = 3)
tapply(d$y, d$group, mean)
#>  ctrl treat
#> 13.50 21.85
var(d$y) / mean(d$y)   # realized overdispersion
#> 4.23

The treatment group averages about 1.5 times the control group, as asked, and the pooled variance-to-mean ratio lands above 1, the overdispersion the negative binomial is there to handle.

A floor-heavy sample

Impaired samples often pile up at zero. Set zi to inject structural zeros on top of the count process:

set.seed(3)
d <- sim_fluency(n = 500, effect = 0, mu = 8, vmr = 2, zi = 0.25)
mean(d$y == 0)          # about a quarter are structural zeros (plus a few from sampling)
#> 0.248

A count with a spike like this is the case for a hurdle model, which separates "did they produce any words at all" from "how many, given some."

true_diff_fluency()

The estimand for bias and coverage is the true difference in expected count. This helper returns it in closed form, so you can score any model against it.

true_diff_fluency(effect = log(1.5), mu = 15)
#> 7.5

# it matches the process it describes
set.seed(9)
d <- sim_fluency(n = 1e5, effect = log(1.5), mu = 15, vmr = 2)
mean(d$y[d$group == "treat"]) - mean(d$y[d$group == "ctrl"])
#> 7.512

The analytic difference (7.5) and the empirical one from a large sample (7.512) agree, which is what lets you use it as ground truth when you measure a model's bias or confidence-interval coverage.