Simulating an Experiment from Scratch

Here is a worked example of simulating a partially crossed design from scratch.

First, some setup:

using DataFrames
using MixedModels, MixedModelsSim
using Random

Assemble the Design

We're going to do a 2 x 2 x 2 design. For concreteness, let's think of this as a linguistic design:

  • age old vs. young, between subjects
  • frequency high vs. low, between items
  • context matched vs. unmatched, within both.

Further, let's assume we want 40 subjects and 80 items. We can specify this design as follows:

n_subj = 40
n_item = 80
subj_btwn = Dict(:age => ["old", "young"])
item_btwn = Dict(:frequency => ["high", "low"])
both_win = Dict(:context => ["matched", "unmatched"])
Dict{Symbol, Vector{String}} with 1 entry:
  :context => ["matched", "unmatched"]

and then generate it:

rng = MersenneTwister(42)  # specify our random number generator for reproducibility
design = simdat_crossed(rng, n_subj, n_item;
                        subj_btwn = subj_btwn,
                        item_btwn = item_btwn,
                        both_win = both_win)
6400-element Vector{NamedTuple{(:subj, :age, :item, :frequency, :context, :dv), Tuple{String, String, String, String, String, Float64}}}:
 (subj = "S01", age = "old", item = "I01", frequency = "high", context = "matched", dv = -0.5560268761463861)
 (subj = "S02", age = "young", item = "I01", frequency = "high", context = "matched", dv = -0.444383357109696)
 (subj = "S03", age = "old", item = "I01", frequency = "high", context = "matched", dv = 0.027155338009193845)
 (subj = "S04", age = "young", item = "I01", frequency = "high", context = "matched", dv = -0.29948409035891055)
 (subj = "S05", age = "old", item = "I01", frequency = "high", context = "matched", dv = 1.7778610980573246)
 (subj = "S06", age = "young", item = "I01", frequency = "high", context = "matched", dv = -1.14490153172882)
 (subj = "S07", age = "old", item = "I01", frequency = "high", context = "matched", dv = -0.46860588216767457)
 (subj = "S08", age = "young", item = "I01", frequency = "high", context = "matched", dv = 0.15614346264074028)
 (subj = "S09", age = "old", item = "I01", frequency = "high", context = "matched", dv = -2.641991008076796)
 (subj = "S10", age = "young", item = "I01", frequency = "high", context = "matched", dv = 1.0033099014594844)
 ⋮
 (subj = "S32", age = "young", item = "I80", frequency = "low", context = "unmatched", dv = -0.0720635904250138)
 (subj = "S33", age = "old", item = "I80", frequency = "low", context = "unmatched", dv = -0.5636260084084584)
 (subj = "S34", age = "young", item = "I80", frequency = "low", context = "unmatched", dv = -0.48302653966207154)
 (subj = "S35", age = "old", item = "I80", frequency = "low", context = "unmatched", dv = 0.5635515836664501)
 (subj = "S36", age = "young", item = "I80", frequency = "low", context = "unmatched", dv = -1.5342115330081054)
 (subj = "S37", age = "old", item = "I80", frequency = "low", context = "unmatched", dv = -0.7925886663164663)
 (subj = "S38", age = "young", item = "I80", frequency = "low", context = "unmatched", dv = 1.1344191438635207)
 (subj = "S39", age = "old", item = "I80", frequency = "low", context = "unmatched", dv = 1.3406462534138157)
 (subj = "S40", age = "young", item = "I80", frequency = "low", context = "unmatched", dv = 0.7969024201220799)

Note that simdat_crossed returns a row table, which MixedModels.jl can process directly. For nicely displaying it, we can again use pretty_table:

pretty_table(first(design, 5))
┌────────┬────────┬────────┬───────────┬─────────┬───────────┐
│   subj │    age │   item │ frequency │ context │        dv │
│ String │ String │ String │    String │  String │   Float64 │
├────────┼────────┼────────┼───────────┼─────────┼───────────┤
│    S01 │    old │    I01 │      high │ matched │ -0.556027 │
│    S02 │  young │    I01 │      high │ matched │ -0.444383 │
│    S03 │    old │    I01 │      high │ matched │ 0.0271553 │
│    S04 │  young │    I01 │      high │ matched │ -0.299484 │
│    S05 │    old │    I01 │      high │ matched │   1.77786 │
└────────┴────────┴────────┴───────────┴─────────┴───────────┘

We can also convert it to a DataFrame and change the factors to use pooled arrays to save a bit of memory.

df = pooled!(DataFrame(design))
first(df, 5)

5 rows × 6 columns

subjageitemfrequencycontextdv
StringStringStringStringStringFloat64
1S01oldI01highmatched-0.556027
2S02youngI01highmatched-0.444383
3S03oldI01highmatched0.0271553
4S04youngI01highmatched-0.299484
5S05oldI01highmatched1.77786

Note that simdat_crossed generated a column dv for our dependent variable that has been pre-populated with noise from a standard normal distribution ($N(0,1)$). Typically, we will want to scale that, but we can do that in the simulation step. Also, this dependent variable is pure noise: we haven't yet added in effects. Adding in effects also comes in the simulation step.

But before we get to simulating, let's fit the model to the noise, just to see how things look. We're going to use effects coding for our contrasts.

contrasts = Dict(:age => EffectsCoding(base="young"),
                 :frequency => EffectsCoding(base="high"),
                 :context => EffectsCoding(base="matched"))
form = @formula(dv ~ 1 + age * frequency * context +
                    (1 + frequency + context | subj) +
                    (1 + age + context | item))
m0 = fit(MixedModel, form, design; contrasts=contrasts)
Est.SEzpσ_itemσ_subj
(Intercept)-0.01680.0128-1.310.18950.00940.0111
age: old-0.00040.0128-0.030.97340.0107
frequency: low0.01800.01311.370.1696 0.0217
context: unmatched0.00990.01320.750.45340.00620.0245
age: old & frequency: low0.00140.01310.110.9147
age: old & context: unmatched0.00680.01320.510.6089
frequency: low & context: unmatched-0.00840.0126-0.670.5043
age: old & frequency: low & context: unmatched0.00000.01260.000.9972
Residual1.0095

Assemble the Random Effects

The hard part in simulating right now is specifying the random effects. We're working on making this bit easier, but you need to specify the variance-covariance matrix of the random effects. You can see what this looks like:

vc = VarCorr(m0)
Column VarianceStd.DevCorr.
item(Intercept)0.000089100.00943915
age: old0.000113750.01066532+1.00
context: unmatched0.000038300.00618839-1.00-1.00
subj(Intercept)0.000123560.01111553
frequency: low0.000472810.02174427+1.00
context: unmatched0.000600980.02451495+1.00+1.00
Residual 1.019117631.00951356

For each grouping variable (subjects and items), there are two major components: the standard deviations ahd the correlations.

For this example, we'll just assume all the correlations and thus the covariances are 0 in order to make things simple. Then we only have to worry about the standard deviations.

Let's assume that the variability

  • between items
    • in the intercept is 1.3 times the residual variability
    • in age is 0.35 times the residual variability
    • in context is 0.75 times the residual variability
  • between subjects
    • in the intercept is 1.5 times the residual variability
    • in frequency is 0.5 times the residual variability
    • in context is 0.75 times the residual variability

Note these are always specified relative to the residual standard deviation. In other words, we think about how big the between-subject and between-item differences are relative to the between-observation differences.

We can now create the associated covariance matrices.[cholesky]

re_item = create_re(1.3, 0.35, 0.75)
3×3 LinearAlgebra.LowerTriangular{Float64, LinearAlgebra.Diagonal{Float64, Vector{Float64}}}:
 1.3   ⋅     ⋅ 
 0.0  0.35   ⋅ 
 0.0  0.0   0.75
re_subj = create_re(1.5, 0.5, 0.75)
3×3 LinearAlgebra.LowerTriangular{Float64, LinearAlgebra.Diagonal{Float64, Vector{Float64}}}:
 1.5   ⋅    ⋅ 
 0.0  0.5   ⋅ 
 0.0  0.0  0.75

We can check that we got these right by installing these parameter values into the model.

update!(m0; subj=re_subj, item=re_item)
VarCorr(m0)
Column VarianceStd.DevCorr.
item(Intercept)1.6454111.282736
age: old0.1192680.345352+0.00
context: unmatched0.5476590.740040+0.00+0.00
subj(Intercept)2.1906361.480080
frequency: low0.2434040.493360+0.00
context: unmatched0.5476590.740040+0.00+0.00
Residual 0.9736160.986720

Looks good. The values don't exactly match the values in our parameter vector because the residual standard deviation isn't exactly 1.0.

For the actual simulation, we'll need the compact form of these covariance matrices that MixedModels.jl uses internally. This compact form is the parameter vector θ and we can get it back out of the model where we just installed it:

show(m0.θ)
[1.3, 0.0, 0.0, 0.35, 0.0, 0.75, 1.5, 0.0, 0.0, 0.5, 0.0, 0.75]

Alternatively, we can also just generate θ directly from the random-effects matrices:

θ = createθ(m0; subj=re_subj, item=re_item)
show(θ)
[1.3, 0.0, 0.0, 0.35, 0.0, 0.75, 1.5, 0.0, 0.0, 0.5, 0.0, 0.75]

We could also create it directly from the covariance matrices we created, but in this case we need to make sure they're in the same order as in the VarCorr output:

θhand = vcat( flatlowertri(re_item), flatlowertri(re_subj) )
show(θhand)
[1.3, 0.0, 0.0, 0.35, 0.0, 0.75, 1.5, 0.0, 0.0, 0.5, 0.0, 0.75]
Warning

In storing the parameter vector θ, MixedModels.jl uses an ordering that yields maximum sparseness, which enables better computational efficiency. The ordering is thus dependent on the entirety of the design – both the choice of the random effects and the relative number of subjects and items. For this reason, we strongly recommend using the helper methods that allow specifying the grouping variable by name.

Assemble the Fixed Effects

The last two components we need are the residual variance and the effect sizes for the fixed effects.

σ = 5;
β = [1.0, -1.0, 2.0, -1.5, 0.3, -1.3, 1.4, 0];

The entries in the β correspond to the coefficients in the model given by

coefnames(m0)
8-element Vector{String}:
 "(Intercept)"
 "age: old"
 "frequency: low"
 "context: unmatched"
 "age: old & frequency: low"
 "age: old & context: unmatched"
 "frequency: low & context: unmatched"
 "age: old & frequency: low & context: unmatched"

Simulate

Now we're ready to actually simulate our data. We can use parametricbootstrap to do this: the parametric bootstrap actually works by simulating new data from an existing model and then looking at how the estimates fit to that new data look. In MixedModels.jl, you can specify different parameter values, such as the ones we made up for our fake data.

# typically we would use a lot more simulations
# but we want to be quick in this example
sim = parametricbootstrap(MersenneTwister(12321), 20, m0; β=β, σ=σ, θ=θ)
MixedModels.MixedModelBootstrap{Float64}(NamedTuple{(:objective, :σ, :β, :se, :θ), Tuple{Float64, Float64, NamedTuple{(Symbol("(Intercept)"), Symbol("age: old"), Symbol("frequency: low"), Symbol("context: unmatched"), Symbol("age: old & frequency: low"), Symbol("age: old & context: unmatched"), Symbol("frequency: low & context: unmatched"), Symbol("age: old & frequency: low & context: unmatched")), NTuple{8, Float64}}, StaticArrays.SVector{8, Float64}, StaticArrays.SVector{12, Float64}}}[(objective = 40018.602265232206, σ = 4.939079891966025, β = ((Intercept) = 2.4058378521152743, age: old = -1.2670883256174414, frequency: low = 2.0533714785795065, context: unmatched = -1.740107860861446, age: old & frequency: low = -0.5190328604292623, age: old & context: unmatched = 0.1451120653754165, frequency: low & context: unmatched = 2.1177660632840314, age: old & frequency: low & context: unmatched = -0.02185888186006084), se = [1.1358881306725912, 0.8670260270935046, 0.8530554501452571, 0.684703133452198, 0.43496653581431405, 0.5427927123493602, 0.4219076849206084, 0.06173849864957531], θ = [1.3802131187697046, 0.07363573065954103, -0.15131224941000407, 0.36547179743397945, -0.016468404773229334, 0.7403318686332546, 1.0755846968318812, -0.1504957254988452, -0.01632933867642655, 0.46025221334153915, -0.09160514667265254, 0.6842445398550104]), (objective = 40186.70354885261, σ = 5.001320547370504, β = ((Intercept) = -0.17367934549402927, age: old = -1.6884737399653018, frequency: low = 2.103092243374021, context: unmatched = -1.866265144821212, age: old & frequency: low = 0.44977298527607595, age: old & context: unmatched = -1.2812982019847243, frequency: low & context: unmatched = 1.0553300004934154, age: old & frequency: low & context: unmatched = -0.05650169868333586), se = [1.5482978350472023, 1.420100545024385, 0.7631451294407995, 0.7300471494057506, 0.449277042096159, 0.6266700184275188, 0.37968650487944894, 0.06251650684213131], θ = [1.1642747813325354, -0.11321712595820184, 0.05407756390151051, 0.3544440893232437, 0.027224412538752983, 0.6670148142686142, 1.7746883670892664, -0.02718461747686232, 0.011056030533929461, 0.4965640769452714, 0.15214834931870347, 0.7736222406015976]), (objective = 40293.01559980909, σ = 5.041204304375276, β = ((Intercept) = 0.6887229123586143, age: old = -1.8037166766668706, frequency: low = 0.651712588283564, context: unmatched = -1.2626074934316618, age: old & frequency: low = 0.05204279786749312, age: old & context: unmatched = -1.9319507254949027, frequency: low & context: unmatched = 2.0210983526236244, age: old & frequency: low & context: unmatched = -0.053056808219879036), se = [1.3062281948849344, 1.0992898764033403, 0.8259810827032392, 0.7586793920756463, 0.42947745486447514, 0.6168968956849005, 0.4460981249242797, 0.06301505380469095], θ = [1.2951438235640418, 0.057276897388244305, 0.04009868812191284, 0.32727649588218694, 0.017184326183367318, 0.7823302983626522, 1.3566790733136107, 0.057508425983256534, 0.06707007932360987, 0.47493608576352003, -0.002400610938860735, 0.7669626412275298]), (objective = 40086.82922469018, σ = 4.95917240448153, β = ((Intercept) = -0.013630950518911557, age: old = -0.4693672440767657, frequency: low = 2.538196869753277, context: unmatched = -2.664125115384292, age: old & frequency: low = 0.8833080201462415, age: old & context: unmatched = -0.7743307686347451, frequency: low & context: unmatched = 1.4427059716859063, age: old & frequency: low & context: unmatched = -0.05042503402891836), se = [1.3515992019984038, 1.1375452437576952, 0.8570998573960862, 0.6086845371293691, 0.44923150412889595, 0.3981812143124084, 0.46453342575439144, 0.061989655056019125], θ = [1.37055947859067, 0.04275566878224288, -0.10385162523364136, 0.37870856747467146, 0.03855219086060012, 0.822907948633788, 1.4232959260973743, 0.11124936779920248, -0.14681702625420018, 0.486808251040582, 0.13303610093568496, 0.4608335024980936]), (objective = 40383.621861971376, σ = 5.089719970565516, β = ((Intercept) = 0.27438629760169836, age: old = 0.25442213666776575, frequency: low = 1.4263281664347123, context: unmatched = -0.4025468262829919, age: old & frequency: low = -0.2741273897667077, age: old & context: unmatched = -0.7396304540571129, frequency: low & context: unmatched = 1.6351514975024959, age: old & frequency: low & context: unmatched = 0.07347458213196409), se = [1.2321365168671636, 0.9878369899542863, 0.8879052140230335, 0.70238008264155, 0.49602136205494207, 0.528166277807658, 0.4673605232514832, 0.06362149963206895], θ = [1.3230917916461893, 0.001452126811210459, -0.13098810456765728, 0.27519533049836226, 0.09222365579329991, 0.7977309485884152, 1.2093958002922522, 0.04329909005529296, 0.26132610089377445, 0.5778504303996129, -0.1053654670244077, 0.5874479970113262]), (objective = 40073.14199479938, σ = 4.958403663549425, β = ((Intercept) = 2.17318487447886, age: old = -1.6792985233710385, frequency: low = 2.1049911769832983, context: unmatched = -2.637335057515626, age: old & frequency: low = 0.34020843195517303, age: old & context: unmatched = -0.8572717880772203, frequency: low & context: unmatched = 1.2576639252661532, age: old & frequency: low & context: unmatched = -0.012957628040293462), se = [1.3399794190249927, 1.098570589475671, 0.8628381161144795, 0.6905442866024954, 0.3947177615306724, 0.580483035175793, 0.37912053981410837, 0.06198004579436781], θ = [1.4299864746348698, -0.031950139939736104, -0.062170738648142945, 0.3581917155625209, -0.036388256587482085, 0.6708231351675298, 1.3757160494786433, -0.03875772239881358, -0.03798155302747845, 0.42552345157520777, -0.16395577563020977, 0.716691311281039]), (objective = 40021.69032342658, σ = 4.931839940075771, β = ((Intercept) = 2.898462024250111, age: old = 0.340318608263236, frequency: low = 3.9239223815245685, context: unmatched = -0.8812095444111226, age: old & frequency: low = 0.7198568091002483, age: old & context: unmatched = -1.0617893533396348, frequency: low & context: unmatched = 1.88464163458346, age: old & frequency: low & context: unmatched = 0.08042828251315401), se = [1.3674322534305918, 1.1800021902631899, 0.8210589062508611, 0.71162568078969, 0.4434770894155967, 0.5934345978660238, 0.39755120855848874, 0.06164799925094714], θ = [1.3066260851695726, -0.005585488649827695, 0.05019247871190811, 0.3698999676674243, -0.12648475951524202, 0.6991484523629521, 1.4883463199104008, 0.10162625872805975, -0.0877453178603434, 0.48828944347021846, -0.0017199623849936609, 0.7517934638046927]), (objective = 40104.20130795826, σ = 4.9630825159354135, β = ((Intercept) = 0.20265194814911563, age: old = -0.20951290640163522, frequency: low = 3.497967969367615, context: unmatched = -2.4187342778288516, age: old & frequency: low = 0.009123554115819974, age: old & context: unmatched = -0.9711464983701439, frequency: low & context: unmatched = 1.8612376860127002, age: old & frequency: low & context: unmatched = 0.04828469357173061), se = [1.347339388021339, 1.137842313914241, 0.8475486365546664, 0.8232794826548814, 0.444635126983239, 0.725835539819035, 0.3934471439476498, 0.06203853144919267], θ = [1.3465801162334816, 0.03762467183631728, 0.003534877932098287, 0.34776748527330625, 0.0953854462783352, 0.693648359737749, 1.4265339873742187, -0.027754559949968775, -0.2419509671207738, 0.5028367039337461, -0.31522004467322495, 0.8314881457696149]), (objective = 40387.83092777358, σ = 5.066914461157119, β = ((Intercept) = 2.344932366372784, age: old = -1.6692821774044326, frequency: low = 1.532154249272835, context: unmatched = -1.8687050543649026, age: old & frequency: low = 0.4920182881211212, age: old & context: unmatched = -0.5723301717677096, frequency: low & context: unmatched = 0.9305139607120827, age: old & frequency: low & context: unmatched = -0.10765078579583905), se = [1.3995633393095053, 1.2170745439499269, 0.816066343916353, 0.7668691248191064, 0.4341165540005697, 0.6252812450180111, 0.4484675269087806, 0.06333643076446399], θ = [1.2883142583504812, -0.08136331430022395, -0.011972344023047707, 0.40645178684757044, 0.06632153235098799, 0.7808108282144253, 1.488518219136205, 0.14314409066762063, 0.04641651476378458, 0.4254035536763971, 0.032985437067503565, 0.7743749930172742]), (objective = 40129.371949855435, σ = 4.960165977824237, β = ((Intercept) = 0.35102886386334387, age: old = -3.1126861700728785, frequency: low = 2.2016052033209608, context: unmatched = -0.326506862346934, age: old & frequency: low = 0.0010668170621371833, age: old & context: unmatched = -0.8500950946166265, frequency: low & context: unmatched = 1.4940348325320945, age: old & frequency: low & context: unmatched = 0.02964255162290107), se = [1.6080056846293094, 1.4056408727725072, 0.9157183195479482, 0.8018799554352706, 0.47820918242995647, 0.7020009937226911, 0.39249245217106, 0.06200207472280296], θ = [1.450016853266934, -0.020855233489133138, -0.01756652599677117, 0.34511451998940135, -0.0665607033159441, 0.6954650603301508, 1.7737762567146265, 0.15943195189724113, -0.20169747610356797, 0.5294873351558672, 0.2221342877655148, 0.8396001042139347]), (objective = 40068.51334676554, σ = 4.967691811544355, β = ((Intercept) = 3.813398857199389, age: old = -0.18991731215687732, frequency: low = 3.836528578129765, context: unmatched = -0.6362405163701056, age: old & frequency: low = 0.7633862303339543, age: old & context: unmatched = -1.4975619739090837, frequency: low & context: unmatched = 1.4152386154087422, age: old & frequency: low & context: unmatched = 0.013249800568367337), se = [1.4200630547882405, 1.2404177998881354, 0.7907908417051353, 0.6480445519003972, 0.3839367057618016, 0.5312455859263734, 0.3762921740271824, 0.062096147644304434], θ = [1.28669015135694, 0.045406727804755895, -0.04327546831212057, 0.3226989735401406, -0.05820544113574143, 0.6642731293553251, 1.56031910352963, -0.02078184760139913, 0.03255902252201527, 0.42326058542975536, -0.02685338400511617, 0.6703852800921171]), (objective = 40521.07801997387, σ = 5.142231880457976, β = ((Intercept) = 1.7647791150674335, age: old = -1.1438873763298354, frequency: low = 2.4591096960544205, context: unmatched = -2.0524799615155693, age: old & frequency: low = 0.6282670328265495, age: old & context: unmatched = -1.8590186710803442, frequency: low & context: unmatched = 2.142001241301484, age: old & frequency: low & context: unmatched = -0.011610630595578453), se = [1.3212518641462356, 1.1669523658669279, 0.7318867824504897, 0.7893391607978463, 0.38952483619898315, 0.6667614882868469, 0.42733719325995745, 0.0642778985057247], θ = [1.1395816925029942, 0.02831406477525492, -0.019124476437426728, 0.3691837110438977, 0.04468684479052629, 0.7332339116939633, 1.408964345924882, 0.05358430565212392, -0.29507387829168824, 0.38968318120194906, -0.07950829364958936, 0.7568805757681676]), (objective = 40327.508774626425, σ = 5.064419899111753, β = ((Intercept) = 0.3844126094127215, age: old = -2.3522764221365335, frequency: low = 2.2576166599690586, context: unmatched = -2.494841669099275, age: old & frequency: low = -0.1543039982622944, age: old & context: unmatched = -0.724325782242761, frequency: low & context: unmatched = 0.991206333069001, age: old & frequency: low & context: unmatched = 0.07881787371977936), se = [1.1836623712097383, 0.9691093574171431, 0.7952952380763351, 0.7480977262881794, 0.41305066677512015, 0.6166186832570038, 0.4282979827806771, 0.06330524873889691], θ = [1.2524647079895463, -0.01591848195720609, -0.15315809440736142, 0.35741976048319785, 0.06200238440043158, 0.7296335195416408, 1.1808642292877651, 0.012977392062653161, -0.20108338569467768, 0.4423317569886971, -0.14789225881000628, 0.7241650442922118]), (objective = 40195.50055859437, σ = 5.017593189312361, β = ((Intercept) = -0.941988776220006, age: old = -1.9157559379803402, frequency: low = 4.019275940717263, context: unmatched = -0.03521013911328396, age: old & frequency: low = 0.5068280939465024, age: old & context: unmatched = -1.3640010050399567, frequency: low & context: unmatched = 1.9071199437809745, age: old & frequency: low & context: unmatched = 0.08520908798123417), se = [1.2280130536417706, 1.0336850487043854, 0.7911948364343102, 0.7285291614830208, 0.4318309729521333, 0.6339351065270019, 0.3644376593668428, 0.06271991486640452], θ = [1.2367475244238146, -0.01699329343154664, -0.17945141861109273, 0.3642321447712303, -0.08133092432992847, 0.6088634257702448, 1.2747207365555744, -0.10768357169640742, 0.17982351924896028, 0.4603840488589546, 0.07526901512627898, 0.7708727814837338]), (objective = 40052.15364792543, σ = 4.964592206351148, β = ((Intercept) = 0.5613955381684937, age: old = -0.9105894983644887, frequency: low = 2.436326451691786, context: unmatched = -1.4867157520005538, age: old & frequency: low = 0.2902045717910467, age: old & context: unmatched = -1.441237661744103, frequency: low & context: unmatched = 1.4912701187492605, age: old & frequency: low & context: unmatched = 0.07130975135619963), se = [1.311766823129537, 1.1371103566287935, 0.7410194159693835, 0.6360064617475678, 0.34842149722891735, 0.5100846039312992, 0.38492731447878387, 0.062057402579389354], θ = [1.2373382620178628, 0.0313572054917603, 0.052564363100947545, 0.3764982221135672, -0.03123603259465424, 0.6816814088634963, 1.4215590629216943, -0.04025803872071197, 0.06616962313137548, 0.34318908721247854, -0.13457622817078008, 0.6273102982661336]), (objective = 40363.50616806291, σ = 5.058759785402437, β = ((Intercept) = 0.002264003189998879, age: old = -2.0708710441186, frequency: low = 1.2263276427674, context: unmatched = -0.8961769799397443, age: old & frequency: low = 0.6971190204887288, age: old & context: unmatched = -0.3187108616125955, frequency: low & context: unmatched = 2.0928780924140984, age: old & frequency: low & context: unmatched = -0.02725510131242313), se = [1.4409817673431768, 1.2463565880712082, 0.8559326025320578, 0.8136295633236179, 0.457817551962977, 0.6668271412083092, 0.470460659049673, 0.06323449731753046], θ = [1.3246233350118495, -0.0609734103175935, 0.0858983630109309, 0.3404240727170446, -0.005713451899168774, 0.8197542675093232, 1.5368768895344376, 0.24247180764748488, 0.10728786846865236, 0.4502932054823649, 0.176955952172991, 0.803708704485366]), (objective = 40352.91727031978, σ = 5.069465280300661, β = ((Intercept) = 0.14618013717592385, age: old = 0.01400791614591642, frequency: low = 1.5614822145157021, context: unmatched = 0.5000970330514182, age: old & frequency: low = 0.07916639522509898, age: old & context: unmatched = -1.3076328065091218, frequency: low & context: unmatched = 1.5503568144570261, age: old & frequency: low & context: unmatched = -0.048257433116977624), se = [1.3890793768551875, 1.2166495089107168, 0.808127347029909, 0.7382208242076331, 0.45140261557441874, 0.6310061348865892, 0.3883513698992174, 0.06336831600375827], θ = [1.2340560404617549, -0.017410733379092092, 0.06528541738609503, 0.3520603321583274, 0.061935169304874055, 0.6699850570963106, 1.4951726565966297, 0.14072021172311355, 0.0446818244744674, 0.47851057245613177, 0.0592330572359787, 0.7797278039832439]), (objective = 40239.96502608555, σ = 5.022020643972092, β = ((Intercept) = -0.20636296964737968, age: old = -2.091149372787316, frequency: low = 1.8193302739687371, context: unmatched = -1.2598168770158826, age: old & frequency: low = 0.19363133423595072, age: old & context: unmatched = -0.5923372828397314, frequency: low & context: unmatched = 1.762678582880754, age: old & frequency: low & context: unmatched = -0.024342188739192975), se = [1.3946134521830145, 1.202019222594715, 0.8209372420859115, 0.7897493079023322, 0.4169430247861582, 0.6734045439245857, 0.41733801956883104, 0.06277525804965116], θ = [1.3023733044239643, -0.03407953464238272, -0.16726504455266616, 0.329711179011346, -0.010209801901895181, 0.7154637802380653, 1.4934343150133649, -0.029932998598507364, 0.10892715519770424, 0.46220240214954667, 0.29025872139442044, 0.785393983522106]), (objective = 40187.56779060092, σ = 4.9811780520504145, β = ((Intercept) = -1.1624453910717687, age: old = -1.3422147173528673, frequency: low = 2.0839571034089004, context: unmatched = -0.6215757087667217, age: old & frequency: low = 0.1717511782258361, age: old & context: unmatched = -0.8010893303336065, frequency: low & context: unmatched = 1.558257921865055, age: old & frequency: low & context: unmatched = -0.0653809260260382), se = [1.3890772425691955, 1.155834212797257, 0.8824796395740723, 0.7506891248818572, 0.43033412125508286, 0.5632738337861612, 0.5001336285912154, 0.06226472565063018], θ = [1.4321669415085623, -0.012022248133624475, 0.08251453780859445, 0.37028456149453204, 0.0955540127280059, 0.8820707980601529, 1.4418151913006478, -0.038307961307116396, -0.008980610069111728, 0.47137911560861, -0.18130643653574469, 0.6872298100950597]), (objective = 40425.43655377477, σ = 5.093668310407084, β = ((Intercept) = -1.796637112207508, age: old = -2.432607318130064, frequency: low = 1.1918804094709436, context: unmatched = -1.8377351505917991, age: old & frequency: low = 0.5170195033031822, age: old & context: unmatched = -1.0168219596159525, frequency: low & context: unmatched = 1.4621476887582376, age: old & frequency: low & context: unmatched = -0.007559495292853527), se = [1.3412599518120363, 1.1692797549142526, 0.8164667620587477, 0.7381500373841334, 0.4846182626183338, 0.6277475774420743, 0.39351294049876046, 0.06367085388008856], θ = [1.2161704982412667, 0.039356238713308636, -0.023802001527406514, 0.38239648698117146, 0.13122029849444122, 0.6687192019164322, 1.4239704194535239, -0.09263549484766839, -0.06316899763597267, 0.5228340145884798, -0.03410082642512904, 0.7720934232487263])], LinearAlgebra.LowerTriangular{Float64, Matrix{Float64}}[[1.3 0.0 0.0; 0.0 0.35 0.0; 0.0 0.0 0.75], [1.5 0.0 0.0; 0.0 0.5 0.0; 0.0 0.0 0.75]], [[1, 2, 3, 5, 6, 9], [1, 2, 3, 5, 6, 9]], [0.0, -Inf, -Inf, 0.0, -Inf, 0.0, 0.0, -Inf, -Inf, 0.0, -Inf, 0.0], (item = ("(Intercept)", "age: old", "context: unmatched"), subj = ("(Intercept)", "frequency: low", "context: unmatched")))

As mentioned above, the ordering within θ is dependent on the entire design, so if you embed the simulation code in a loop iterating over different numbers of subjects and items, it's probably better to write it as:

sim = parametricbootstrap(MersenneTwister(12321), 20, m0;
                          β=β, σ=σ, θ=createθ(m0; subj=re_subj, item=re_item))

See your power and profit!

Finally, we can turn this into a power table:

ptbl = power_table(sim)
8-element Vector{NamedTuple{(:coefname, :power), Tuple{String, Float64}}}:
 (coefname = "age: old & frequency: low", power = 0.1)
 (coefname = "frequency: low & context: unmatched", power = 1.0)
 (coefname = "age: old", power = 0.15)
 (coefname = "context: unmatched", power = 0.5)
 (coefname = "(Intercept)", power = 0.15)
 (coefname = "frequency: low", power = 0.7)
 (coefname = "age: old & context: unmatched", power = 0.35)
 (coefname = "age: old & frequency: low & context: unmatched", power = 0.0)

For nicely displaying it, we can again use pretty_table:

pretty_table(ptbl)
┌────────────────────────────────────────────────┬─────────┐
│                                       coefname │   power │
│                                         String │ Float64 │
├────────────────────────────────────────────────┼─────────┤
│                      age: old & frequency: low │     0.1 │
│            frequency: low & context: unmatched │     1.0 │
│                                       age: old │    0.15 │
│                             context: unmatched │     0.5 │
│                                    (Intercept) │    0.15 │
│                                 frequency: low │     0.7 │
│                  age: old & context: unmatched │    0.35 │
│ age: old & frequency: low & context: unmatched │     0.0 │
└────────────────────────────────────────────────┴─────────┘

We can of course convert the row table into a DataFrame:

DataFrame(ptbl)

8 rows × 2 columns

coefnamepower
StringFloat64
1age: old & frequency: low0.1
2frequency: low & context: unmatched1.0
3age: old0.15
4context: unmatched0.5
5(Intercept)0.15
6frequency: low0.7
7age: old & context: unmatched0.35
8age: old & frequency: low & context: unmatched0.0
  • choleskyTechnically, we're creating the lower Cholesky factor of these matrices, which is a bit like the matrix square root. In other words, we're creating the matrix form of standard deviations instead of the matrix form of the variances.