Power Analysis and Simulation Tutorial

This tutorial demonstrates how to conduct power analyses and data simulation using Julia and the MixedModelsSim package.

Power analysis is an important tool for planning an experimental design. Here we show how to

  1. Take existing data and calculate power by simulating new data.
  2. Adapt parameters in a given Linear Mixed Model to analyze power without changing the existing data set.
  3. Create a (simple) balanced fully crossed dataset from scratch and analyze power.
  4. Recreate a more complex dataset from scratch and analyze power for specific model parameter but various sample sizes.

Load the packages we'll be using in Julia

First, here are the packages needed in this example.

using MixedModels        # run mixed models
using MixedModelsSim     # simulation functions for mixed models
using DataFrames, Tables # work with data tables
using StableRNGs         # random number generator
using CSV                # write CSV files
using Statistics         # basic math funcions
using DataFramesMeta     # dplyr-like operations
using Gadfly             # plotting package

Define number of iterations

Here we define how many model simulations we want to do. A large number will give more reliable results, but will take longer to compute. It is useful to set it to a low number for testing, and increase it for your final analysis.

nsims = 500
500

Take existing data and calculate power by simulate new data with bootstrapping.

Build a Linear Mixed Model from existing data.

For the first example we are going to simulate bootstrapped data from an existing data set:

Experiment 2 from Kronmüller, E., & Barr, D. J. (2007). Perspective-free pragmatics: Broken precedents and the recovery-from-preemption hypothesis. Journal of Memory and Language, 56(3), 436-455.

The data we will be using through out this tutorial is a study about how in a conversation the change of a speaker or the change of precedents (which are patterns of word usage to discribe an object, e.g. one can refer to the same object "white shoes", "runners", "sneakers") affects the understanding.

Objects are presented on a screen while participants listen to instructions to move the objects around. Participants eye movements are tracked. The dependent variable is response time, defined as the latency between the onset of the test description and the moment at which the target was selected. The independent variables are speaker (old vs. new), precedents (maintain vs. break) and cognitive load (a secondary memory task).

We have to load the data and define some characteristics like the contrasts and the underlying model.

Load existing data:

kb07 = MixedModels.dataset(:kb07);
Arrow.Table: (subj = ["S030", "S030", "S030", "S030", "S030", "S030", "S030", "S030", "S030", "S030"  …  "S103", "S103", "S103", "S103", "S103", "S103", "S103", "S103", "S103", "S103"], item = ["I01", "I02", "I03", "I04", "I05", "I06", "I07", "I08", "I09", "I10"  …  "I23", "I24", "I25", "I26", "I27", "I28", "I29", "I30", "I31", "I32"], spkr = ["new", "old", "old", "new", "new", "old", "old", "new", "new", "old"  …  "old", "old", "new", "new", "old", "old", "new", "new", "old", "old"], prec = ["break", "maintain", "break", "maintain", "break", "maintain", "break", "maintain", "break", "maintain"  …  "maintain", "break", "maintain", "break", "maintain", "break", "maintain", "break", "maintain", "break"], load = ["yes", "no", "no", "no", "no", "yes", "yes", "yes", "yes", "no"  …  "yes", "yes", "yes", "yes", "no", "no", "no", "no", "yes", "yes"], rt_trunc = Int16[2267, 3856, 1567, 1732, 2660, 2763, 3528, 1741, 3692, 1949  …  2706, 4281, 2075, 3179, 1216, 2286, 1202, 1581, 1601, 1941], rt_raw = Int16[2267, 3856, 1567, 1732, 2660, 2763, 3528, 1741, 3692, 1949  …  2706, 4281, 2075, 3179, 1216, 2286, 1202, 1581, 1601, 1941])

Set contrasts:

contrasts = Dict(:spkr => HelmertCoding(),
                 :prec => HelmertCoding(),
                 :load => HelmertCoding());
Dict{Symbol, HelmertCoding} with 3 entries:
  :spkr => HelmertCoding(nothing, nothing)
  :load => HelmertCoding(nothing, nothing)
  :prec => HelmertCoding(nothing, nothing)

The chosen LMM for this dataset is defined by the following model formula:

kb07_f = @formula( rt_trunc ~ 1 + spkr+prec+load + (1|subj) + (1+prec|item) );
FormulaTerm
Response:
  rt_trunc(unknown)
Predictors:
  1
  spkr(unknown)
  prec(unknown)
  load(unknown)
  (subj)->1 | subj
  (prec,item)->(1 + prec) | item

Fit the model

kb07_m = fit(MixedModel, kb07_f, kb07; contrasts=contrasts)
Linear mixed model fit by maximum likelihood
 rt_trunc ~ 1 + spkr + prec + load + (1 | subj) + (1 + prec | item)
    logLik   -2 logLik      AIC         AICc        BIC     
 -14331.9251  28663.8501  28681.8501  28681.9513  28731.2548

Variance components:
             Column      Variance  Std.Dev.  Corr.
item     (Intercept)     133015.244 364.713
         prec: maintain   63766.936 252.521 -0.70
subj     (Intercept)      88819.436 298.026
Residual                 462443.388 680.032
 Number of obs: 1789; levels of grouping factors: 32, 56

  Fixed-effects parameters:
──────────────────────────────────────────────────────
                    Coef.  Std. Error      z  Pr(>|z|)
──────────────────────────────────────────────────────
(Intercept)     2181.85       77.4681  28.16    <1e-99
spkr: old         67.879      16.0785   4.22    <1e-04
prec: maintain  -333.791      47.4472  -7.03    <1e-11
load: yes         78.5904     16.0785   4.89    <1e-05
──────────────────────────────────────────────────────

Simulate from existing data with same model parameters

We will first look at the power of the dataset with the same parameters as in the original data set. This means that each dataset will have the exact number of observations as the original data. Here, we use the model kb07_m we fitted above to our dataset kb07.

You can use the parameparametricbootstrap() function to run nsims iterations of data sampled using the parameters from kb07_m. Set up a random seed to make the simulation reproducible. You can use your favourite number.

To use multithreading, you need to set the number of cores you want to use. E.g. in Visual Studio Code, open the settings (gear icon in the lower left corner or cmd-,) and search for "thread". Set julia.NumThreads to the number of cores you want to use (at least 1 less than your total number).

Set random seed for reproducibility:

rng = StableRNG(42);
StableRNGs.LehmerRNG(state=0x00000000000000000000000000000055)

Run nsims iterations:

kb07_sim = parametricbootstrap(rng, nsims, kb07_m; use_threads = false);
MixedModelBootstrap{Float64}(NamedTuple{(:objective, :σ, :β, :se, :θ), Tuple{Float64, Float64, NamedTuple{(Symbol("(Intercept)"), Symbol("spkr: old"), Symbol("prec: maintain"), Symbol("load: yes")), NTuple{4, Float64}}, StaticArrays.SVector{4, Float64}, StaticArrays.SVector{4, Float64}}}[(objective = 28597.686958301765, σ = 672.4742168806578, β = ((Intercept) = 2284.6127147658253, spkr: old = 29.11180730984709, prec: maintain = -390.33009138582776, load: yes = 53.14456254423481), se = [62.73035162232511, 15.899745546263143, 52.115957259123796, 15.899727273445112], θ = [0.38945863611454135, -0.3655273674264099, 0.20173044439091598, 0.4365243927487813]), (objective = 28693.335543496614, σ = 688.590203341095, β = ((Intercept) = 2123.430272091411, spkr: old = 64.88529182702143, prec: maintain = -274.9658380129044, load: yes = 83.45120453339561), se = [86.28344283642528, 16.280827881454552, 48.77442466886997, 16.280806376482772], θ = [0.6432988132115198, -0.2975099395803221, 0.2327012752355757, 0.3517873166617583]), (objective = 28679.767401170287, σ = 682.3511111198069, β = ((Intercept) = 2107.372424544438, spkr: old = 51.6008209756879, prec: maintain = -291.49513162221433, load: yes = 92.29190056137108), se = [72.75295292298387, 16.133363949049304, 43.211172669720334, 16.1333419454115], θ = [0.4427097060208063, -0.1985517356380322, 0.2664915108239762, 0.5121705841466969]), (objective = 28787.634493298312, σ = 707.037344859215, β = ((Intercept) = 2057.47333806114, spkr: old = 55.31045925394355, prec: maintain = -311.81630553494443, load: yes = 80.51174518717448), se = [70.99720815205329, 16.717017401611297, 51.06179555014658, 16.716997368167075], θ = [0.4792278411074963, -0.2223370362887258, 0.31555949817084405, 0.36255927762576184]), (objective = 28684.027373341174, σ = 689.3325316660917, β = ((Intercept) = 2197.0272881011633, spkr: old = 86.62348694185611, prec: maintain = -339.4726787051471, load: yes = 90.57445775392998), se = [69.03652170078819, 16.298384842130684, 44.17932330926803, 16.298364555381887], θ = [0.4834302731831347, -0.21511056465269132, 0.25938299391310804, 0.3484111423976133]), (objective = 28641.72686189496, σ = 680.3774659992224, β = ((Intercept) = 2216.2972774060627, spkr: old = 72.65367174561966, prec: maintain = -386.6643347672149, load: yes = 78.83130004698582), se = [63.3665779433807, 16.08665577056877, 36.70245113648232, 16.08663386167352], θ = [0.35070128691924773, -0.17290925628010323, 0.21291545842204537, 0.48908599444203127]), (objective = 28612.9086577733, σ = 670.9778647263602, β = ((Intercept) = 2188.7111918232044, spkr: old = 120.42387137964064, prec: maintain = -346.09693758195596, load: yes = 65.35910670712612), se = [77.03879394244925, 15.86445041507165, 40.818582006583455, 15.864428655497958], θ = [0.5489807867385973, -0.17648268753075882, 0.2634228604570635, 0.4236836190689666]), (objective = 28622.694996843744, σ = 668.9668293756872, β = ((Intercept) = 2208.695383143965, spkr: old = 46.16091587208839, prec: maintain = -350.6051382036844, load: yes = 83.68193326705679), se = [77.66868578097674, 15.816917741249, 44.74030164376037, 15.816895681988772], θ = [0.48907293269911223, -0.2307484149635346, 0.2683260222520525, 0.5522458171804614]), (objective = 28658.930864778413, σ = 682.4035452019535, β = ((Intercept) = 2048.2630376571633, spkr: old = 101.35055952327083, prec: maintain = -259.8256594089462, load: yes = 66.9053082634049), se = [79.5531152672955, 16.134538170346904, 48.26452874147284, 16.134515983021195], θ = [0.571992548568578, -0.3201351609994257, 0.19924790936693496, 0.396483723555056]), (objective = 28648.328109753806, σ = 675.3312691818315, β = ((Intercept) = 2186.1625524611404, spkr: old = 76.70635645471339, prec: maintain = -330.37912342177265, load: yes = 104.48896278769718), se = [82.22564573673316, 15.967386349861588, 48.71325374664613, 15.967364730790928], θ = [0.583947517030748, -0.2693202514012847, 0.275819273713866, 0.44958571365351685])  …  (objective = 28666.26900735583, σ = 687.7590494010423, β = ((Intercept) = 2081.0346693166866, spkr: old = 82.88994551351107, prec: maintain = -255.20341144884802, load: yes = 80.579325247257), se = [64.79881921514973, 16.261162552833305, 35.27794025003453, 16.26113970749452], θ = [0.41922679351611075, -0.16930998297430983, 0.19400997657230806, 0.39779010405883064]), (objective = 28653.616447037053, σ = 678.726219891093, β = ((Intercept) = 2251.881018346555, spkr: old = 99.51878671101528, prec: maintain = -387.4206859730128, load: yes = 87.21822067641467), se = [86.43486653869887, 16.04759601566194, 46.888319116472765, 16.04757184326136], θ = [0.6147519867111627, -0.32287334634669623, 0.17487699132199272, 0.46424637900488175]), (objective = 28722.61071065539, σ = 696.3334699949901, β = ((Intercept) = 2105.574486994814, spkr: old = 46.874623748600015, prec: maintain = -265.0136789956263, load: yes = 96.44671621852284), se = [59.35294033978519, 16.463936198018043, 42.60072418831378, 16.463916268627592], θ = [0.32738068103921314, -0.14133336588722986, 0.2861929200029339, 0.43357384834259505]), (objective = 28604.97231715837, σ = 669.7787176024134, β = ((Intercept) = 2060.1586282613334, spkr: old = 65.54255150440737, prec: maintain = -288.8327861598166, load: yes = 81.47871040602283), se = [69.5350143496375, 15.836091955352247, 45.54027990603889, 15.836071198125836], θ = [0.45516210649366706, -0.2395049797098539, 0.2696036282179779, 0.45795113123305237]), (objective = 28653.82247865642, σ = 678.2513135732314, β = ((Intercept) = 2305.1418570938777, spkr: old = 52.2020007138039, prec: maintain = -309.63432594337786, load: yes = 94.85995983698129), se = [87.41871083385871, 16.036400791871184, 53.452798652171595, 16.036379760837704], θ = [0.6547546103607714, -0.33715844266109435, 0.25920349798874776, 0.38567765654572983]), (objective = 28687.91150541678, σ = 686.1788806033582, β = ((Intercept) = 2077.763224970452, spkr: old = 83.5621808917557, prec: maintain = -192.4407390891706, load: yes = 83.6357986730465), se = [68.59114074209404, 16.223841384367148, 51.93831255390963, 16.22382131882905], θ = [0.4361490272812234, -0.2984333646495293, 0.2763802075257269, 0.44199844546604783]), (objective = 28678.019205632714, σ = 681.7363191444938, β = ((Intercept) = 2222.0251973850836, spkr: old = 47.85011402497766, prec: maintain = -367.94290868234367, load: yes = 72.98866859494171), se = [80.55475621971621, 16.118816928562516, 51.32175463142586, 16.11879548392605], θ = [0.5625598983346031, -0.3007519693815502, 0.2702044128793771, 0.4435549981108227]), (objective = 28606.440846947407, σ = 672.5294266924902, β = ((Intercept) = 2227.19675546103, spkr: old = 40.55671148115136, prec: maintain = -342.5952669729612, load: yes = 72.06259868418597), se = [64.52690329432359, 15.901107151137774, 45.85407447567215, 15.901087172192245], θ = [0.4124514740184291, -0.25621561333548093, 0.2553893591394055, 0.4318716198580464]), (objective = 28663.156694777088, σ = 688.0170824276888, β = ((Intercept) = 2144.3737315034095, spkr: old = 61.69102196978499, prec: maintain = -291.71255079867694, load: yes = 85.37277915376215), se = [74.07392413996861, 16.26721198145973, 38.3292342964013, 16.267187613476505], θ = [0.5277447167029292, -0.2443287316649958, 0.14740800586657923, 0.3611198742407753]), (objective = 28688.217993539572, σ = 688.1234716256314, β = ((Intercept) = 2135.159566436118, spkr: old = 87.15843669270461, prec: maintain = -341.3169839690801, load: yes = 57.958495645257756), se = [74.51542187365884, 16.26979316677296, 40.77012130162073, 16.269769808496505], θ = [0.499197484759422, -0.23484447082995188, 0.19821879083807736, 0.4350501332116482])], LowerTriangular{Float64, Matrix{Float64}}[[0.5363168233715857 0.0; -0.25981993107379 0.2653016002105174], [0.4382528181348316]], [[1, 2, 4], [1]], [0.0, -Inf, 0.0, 0.0], (item = ("(Intercept)", "prec: maintain"), subj = ("(Intercept)",)))

Try: Run the code above with or without use_threads = true.

The output DataFrame kb07_sim contains the results of the bootstrapping procedure.

df = DataFrame(kb07_sim.allpars);
first(df, 9)
nrow(df)
4500

The dataframe df has 4500 rows: 9 parameters, each from 500 iterations.

Plot some bootstrapped parameter:

σres = @where(df, :type .== "σ", :group .== "residual").value
plot(x = σres, Geom.density, Guide.xlabel("Parametric bootstrap estimates of σ"), Guide.ylabel("Density"))

βInt = @where(df, :type .== "β", :names .== "(Intercept)").value
plot(x = βInt, Geom.density, Guide.xlabel("Parametric bootstrap estimates of β (Intercept)"), Guide.ylabel("Density"))

βSpeaker = @where(df, :type .== "β", :names .== "spkr: old").value
plot(x = βSpeaker, Geom.density, Guide.xlabel("Parametric bootstrap estimates of β Speaker"), Guide.ylabel("Density"))

βPrecedents = @where(df, :type .== "β", :names .== "prec: maintain").value
plot(x = βPrecedents, Geom.density, Guide.xlabel("Parametric bootstrap estimates of β Precedents"), Guide.ylabel("Density"))

βLoad = @where(df, :type .== "β", :names .== "load: yes").value
plot(x = βLoad, Geom.density, Guide.xlabel("Parametric bootstrap estimates of β Load"), Guide.ylabel("Density"))
Parametric bootstrap estimates of β Load -200 -150 -100 -50 0 50 100 150 200 250 300 350 -150 -140 -130 -120 -110 -100 -90 -80 -70 -60 -50 -40 -30 -20 -10 0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 170 180 190 200 210 220 230 240 250 260 270 280 290 300 -200 0 200 400 -150 -145 -140 -135 -130 -125 -120 -115 -110 -105 -100 -95 -90 -85 -80 -75 -70 -65 -60 -55 -50 -45 -40 -35 -30 -25 -20 -15 -10 -5 0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95 100 105 110 115 120 125 130 135 140 145 150 155 160 165 170 175 180 185 190 195 200 205 210 215 220 225 230 235 240 245 250 255 260 265 270 275 280 285 290 295 300 h,j,k,l,arrows,drag to pan i,o,+,-,scroll,shift-drag to zoom r,dbl-click to reset c for coordinates ? for help ? -0.04 -0.03 -0.02 -0.01 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 -0.030 -0.028 -0.026 -0.024 -0.022 -0.020 -0.018 -0.016 -0.014 -0.012 -0.010 -0.008 -0.006 -0.004 -0.002 0.000 0.002 0.004 0.006 0.008 0.010 0.012 0.014 0.016 0.018 0.020 0.022 0.024 0.026 0.028 0.030 0.032 0.034 0.036 0.038 0.040 0.042 0.044 0.046 0.048 0.050 0.052 0.054 0.056 0.058 0.060 -0.03 0.00 0.03 0.06 -0.030 -0.029 -0.028 -0.027 -0.026 -0.025 -0.024 -0.023 -0.022 -0.021 -0.020 -0.019 -0.018 -0.017 -0.016 -0.015 -0.014 -0.013 -0.012 -0.011 -0.010 -0.009 -0.008 -0.007 -0.006 -0.005 -0.004 -0.003 -0.002 -0.001 0.000 0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 0.010 0.011 0.012 0.013 0.014 0.015 0.016 0.017 0.018 0.019 0.020 0.021 0.022 0.023 0.024 0.025 0.026 0.027 0.028 0.029 0.030 0.031 0.032 0.033 0.034 0.035 0.036 0.037 0.038 0.039 0.040 0.041 0.042 0.043 0.044 0.045 0.046 0.047 0.048 0.049 0.050 0.051 0.052 0.053 0.054 0.055 0.056 0.057 0.058 0.059 0.060 Density

Convert p-values of your fixed-effects parameters to dataframe

kb07_sim_df = DataFrame(kb07_sim.coefpvalues);

Have a look at your simulated data:

first(kb07_sim_df, 8)

8 rows × 6 columns

itercoefnameβsezp
Int64SymbolFloat64Float64Float64Float64
11(Intercept)2284.6162.730436.41962.08616e-290
21spkr: old29.111815.89971.830960.0671064
31prec: maintain-390.3352.116-7.489656.90594e-14
41load: yes53.144615.89973.342480.000830325
52(Intercept)2123.4386.283424.60999.88744e-134
62spkr: old64.885316.28083.985386.73721e-5
72prec: maintain-274.96648.7744-5.63751.72536e-8
82load: yes83.451216.28085.125742.96369e-7

Now that we have a bootstrapped data, we can start our power calculation.

Power calculation

The function power_table() from MixedModelsSim takes the output of parametricbootstrap() and calculates the proportion of simulations where the p-value is less than alpha for each coefficient. You can set the alpha argument to change the default value of 0.05 (justify your alpha).

ptbl = power_table(kb07_sim, 0.05)
4-element Vector{NamedTuple{(:coefname, :power), Tuple{String, Float64}}}:
 (coefname = "(Intercept)", power = 1.0)
 (coefname = "spkr: old", power = 0.99)
 (coefname = "prec: maintain", power = 1.0)
 (coefname = "load: yes", power = 1.0)

An estimated power of 1 means that in every iteration the specific parameter we are looking at was below our alpha. An estimated power of 0.5 means that in half of our iterations the specific parameter we are looking at was below our alpha. An estimated power of 0 means that for none of our iterations the specific parameter we are looking at was below our alpha.

You can also do it manually:

kb07_sim_df[kb07_sim_df.coefname .== Symbol("prec: maintain"),:]

mean(kb07_sim_df[kb07_sim_df.coefname .== Symbol("prec: maintain"),:p] .< 0.05)
1.0

For nicely displaying, you can use pretty_table:

pretty_table(ptbl)
┌────────────────┬─────────┐
│       coefname │   power │
│         String │ Float64 │
├────────────────┼─────────┤
│    (Intercept) │     1.0 │
│      spkr: old │    0.99 │
│ prec: maintain │     1.0 │
│      load: yes │     1.0 │
└────────────────┴─────────┘

Adapt parameters in a given Linear Mixed Model to analyze power without changing the existing data set.

Let's say we want to check our power to detect effects of spkr, prec, and load that are only half the size as in our pilot data. We can set a new vector of beta values with the β argument to parametricbootstrap().

Set random seed for reproducibility:

rng = StableRNG(42);
StableRNGs.LehmerRNG(state=0x00000000000000000000000000000055)

Specify β:

new_beta = kb07_m.β
new_beta[2:4] = kb07_m.β[2:4]/2
3-element Vector{Float64}:
   33.939510821966415
 -166.89531060660704
   39.295220532876364

Run nsims iterations:

kb07_sim_half = parametricbootstrap(rng, nsims, kb07_m; β = new_beta, use_threads = false);
MixedModelBootstrap{Float64}(NamedTuple{(:objective, :σ, :β, :se, :θ), Tuple{Float64, Float64, NamedTuple{(Symbol("(Intercept)"), Symbol("spkr: old"), Symbol("prec: maintain"), Symbol("load: yes")), NTuple{4, Float64}}, StaticArrays.SVector{4, Float64}, StaticArrays.SVector{4, Float64}}}[(objective = 28597.686958301798, σ = 672.4742168944547, β = ((Intercept) = 2284.6127147659076, spkr: old = -4.827703512106587, prec: maintain = -223.4347807791827, load: yes = 13.849342011224959), se = [62.73035163000962, 15.899745546588989, 52.11595731959365, 15.899727273770969], θ = [0.3894586364676791, -0.36552736817078296, 0.201730444128044, 0.4365243923112756]), (objective = 28693.335543496774, σ = 688.5902031097281, β = ((Intercept) = 2123.4302720928285, spkr: old = 30.945781003180457, prec: maintain = -108.07052740814592, load: yes = 44.15598399903543), se = [86.28344431104688, 16.280827875983537, 48.7744228213249, 16.28080637101094], θ = [0.6432988225602976, -0.29750993126729297, 0.23270125996100824, 0.35178733030045173]), (objective = 28679.767401187004, σ = 682.3514263744722, β = ((Intercept) = 2107.372425589631, spkr: old = 17.66130866381283, prec: maintain = -124.59982250550942, load: yes = 52.99667898326563), se = [72.75236107124213, 16.13337140201527, 43.21064915028541, 16.13334939837451], θ = [0.4427057268608347, -0.1985448362021727, 0.26649059459950547, 0.512165917389774]), (objective = 28787.634493298334, σ = 707.0373448738048, β = ((Intercept) = 2057.473338061341, spkr: old = 21.370948432159132, prec: maintain = -144.92099492815092, load: yes = 41.216524654074654), se = [70.99720816631935, 16.717017401956014, 51.06179555600408, 16.71699736851179], θ = [0.47922784124926554, -0.22233703710728853, 0.3155594976438599, 0.36255927757864376]), (objective = 28684.027373341156, σ = 689.3325319007566, β = ((Intercept) = 2197.027288098481, spkr: old = 52.68397611901527, prec: maintain = -172.5773680994197, load: yes = 51.27923722375926), se = [69.03652139582431, 16.298384847676765, 44.17932250966241, 16.29836456092776], θ = [0.4834302687956428, -0.2151105627171438, 0.2593829861741756, 0.348411145381006]), (objective = 28641.72686189494, σ = 680.3774660335392, β = ((Intercept) = 2216.29727740529, spkr: old = 38.714160923317934, prec: maintain = -219.76902416094478, load: yes = 39.53607951488223), se = [63.366578260154526, 16.086655771378847, 36.70245102413327, 16.086633862483442], θ = [0.35070129290688035, -0.1729092582671119, 0.2129154554475641, 0.48908599184332]), (objective = 28612.90865777327, σ = 670.9778648154439, β = ((Intercept) = 2188.711191823509, spkr: old = 86.48436055731042, prec: maintain = -179.2016269757223, load: yes = 26.063886173981256), se = [77.03879377109473, 15.864450417177677, 40.81858206236424, 15.86442865760405], θ = [0.5489807862507947, -0.1764826863542275, 0.263422861799965, 0.4236836160681758]), (objective = 28622.694996843737, σ = 668.9668293728637, β = ((Intercept) = 2208.695383144031, spkr: old = 12.22140505015063, prec: maintain = -183.7098275970499, load: yes = 44.386712734118895), se = [77.66868576815821, 15.81691774118225, 44.74030163468216, 15.816895681922027], θ = [0.489072932346809, -0.23074841485472847, 0.2683260222396354, 0.552245817506632]), (objective = 28658.930864778384, σ = 682.4035453521794, β = ((Intercept) = 2048.2630376554903, spkr: old = 67.41104870109945, prec: maintain = -92.93034880254854, load: yes = 27.61008773220655), se = [79.55311556319316, 16.13453817389694, 48.26452830436431, 16.134515986571106], θ = [0.571992552158704, -0.32013515804678266, 0.19924790665817596, 0.396483721208332]), (objective = 28648.328109753813, σ = 675.3312690093388, β = ((Intercept) = 2186.162552460998, spkr: old = 42.76684563281668, prec: maintain = -163.48381281507824, load: yes = 65.19374225491755), se = [82.22564611131742, 15.967386345784131, 48.713254128443914, 15.967364726713498], θ = [0.5839475205064165, -0.269320254492129, 0.2758192755812407, 0.44958571463695424])  …  (objective = 28666.269007355815, σ = 687.7590494965499, β = ((Intercept) = 2081.034669317012, spkr: old = 48.950434690733296, prec: maintain = -88.30810084304792, load: yes = 41.28410471404286), se = [64.7988189098356, 16.261162555091015, 35.27793999308063, 16.261139709752193], θ = [0.41922678965216487, -0.16930998032828004, 0.19400997566029493, 0.3977901051234942]), (objective = 28653.61644703708, σ = 678.7262183707182, β = ((Intercept) = 2251.8810183391815, spkr: old = 65.57927588258265, prec: maintain = -220.52537537285232, load: yes = 47.92300015085551), se = [86.43486911082705, 16.047595979731533, 46.888319842569054, 16.04757180733107], θ = [0.6147520062852723, -0.3228733476865012, 0.17487700432785466, 0.4642463962417016]), (objective = 28722.610710655375, σ = 696.3334697799232, β = ((Intercept) = 2105.5744869955697, spkr: old = 12.93511292640505, prec: maintain = -98.1183683892466, load: yes = 57.15149568488469), se = [59.35294069865901, 16.463936192933943, 42.600724493143545, 16.463916263543467], θ = [0.32738068422073585, -0.14133337048824288, 0.28619292085457043, 0.4335738501022161]), (objective = 28604.972317158357, σ = 669.7787171381538, β = ((Intercept) = 2060.1586282615144, spkr: old = 31.60304068377948, prec: maintain = -121.93747555187488, load: yes = 42.18348987300281), se = [69.53501491669756, 15.836091944379358, 45.54028061390293, 15.836071187153019], θ = [0.4551621116759513, -0.23950498169119316, 0.26960363536734533, 0.4579511338811007]), (objective = 28653.822478656402, σ = 678.2513136362046, β = ((Intercept) = 2305.141857093361, spkr: old = 18.262489891366787, prec: maintain = -142.73901533727002, load: yes = 55.564739304696), se = [87.41871067653915, 16.03640079335928, 53.45279875779163, 16.036379762325808], θ = [0.6547546088680261, -0.3371584445758443, 0.2592034969420603, 0.38567765641579754]), (objective = 28687.911505416767, σ = 686.1788805897905, β = ((Intercept) = 2077.7632249704584, spkr: old = 49.62267006981797, prec: maintain = -25.545428482539435, load: yes = 44.340578140181364), se = [68.59114077986216, 16.22384138404646, 51.938312523829, 16.223821318508353], θ = [0.4361490275262046, -0.2984333644124711, 0.27638020741062713, 0.4419984457651159]), (objective = 28678.019205632645, σ = 681.7363187743792, β = ((Intercept) = 2222.0251973845125, spkr: old = 13.910603201283372, prec: maintain = -201.04759807747627, load: yes = 33.6934480626495), se = [80.5547568674323, 16.118816919812502, 51.32175512495205, 16.118795475175943], θ = [0.5625599012189993, -0.3007519770838569, 0.2702044111247767, 0.44355500683954535]), (objective = 28606.440846947422, σ = 672.529426591652, β = ((Intercept) = 2227.196755466824, spkr: old = 6.617200665066487, prec: maintain = -175.69995636046087, load: yes = 32.767378145488706), se = [64.52690295406569, 15.901107148756894, 45.85407473847266, 15.901087169811506], θ = [0.4124514715773318, -0.2562156090758654, 0.2553893668384506, 0.431871617822354]), (objective = 28663.15669477711, σ = 688.0170827534774, β = ((Intercept) = 2144.373731502584, spkr: old = 27.75151114861961, prec: maintain = -124.8172401912601, load: yes = 46.077558621678705), se = [74.07392275078402, 16.267211989161847, 38.32923360407811, 16.26718762117864], θ = [0.5277447033586398, -0.24432872423307533, 0.1474080056965562, 0.3611198738068936]), (objective = 28688.217993539598, σ = 688.1234716400177, β = ((Intercept) = 2135.1595664362185, spkr: old = 53.218925870951715, prec: maintain = -174.4216733622509, load: yes = 18.6632751122517), se = [74.51542170624657, 16.269793167113416, 40.770121298977614, 16.269769808836998], θ = [0.4991974830697586, -0.23484447024995014, 0.19821879147666122, 0.435050133181829])], LowerTriangular{Float64, Matrix{Float64}}[[0.5363168233715857 0.0; -0.25981993107379 0.2653016002105174], [0.4382528181348316]], [[1, 2, 4], [1]], [0.0, -Inf, 0.0, 0.0], (item = ("(Intercept)", "prec: maintain"), subj = ("(Intercept)",)))

Power calculation

power_table(kb07_sim_half)
4-element Vector{NamedTuple{(:coefname, :power), Tuple{String, Float64}}}:
 (coefname = "(Intercept)", power = 1.0)
 (coefname = "spkr: old", power = 0.552)
 (coefname = "prec: maintain", power = 0.944)
 (coefname = "load: yes", power = 0.704)

Create a (simple) balanced fully crossed dataset from scratch and analyze power.

In some situations, instead of using an existing dataset it may be useful to simulate the data from scratch. This could be the case when the original data is not available but the effect sizes are known. That means that we have to:

a) specify the effect sizes manually

b) manually create an experimental design, according to which data can be simulated

If we simulate data from scratch, aside from subject and item number, we can manipulate the arguments β, σ and θ. Lets have a closer look at them, define their meaning and we will see where the corresponding values in the model output are.

Fixed Effects (βs)

β are our effect sizes. If we look again on our LMM summary from the kb07-dataset kb07_m we see our four β under fixed-effects parameters in the Coef.-column.

kb07_m
kb07_m.β
4-element Vector{Float64}:
 2181.852639291391
   67.87902164393283
 -333.7906212132141
   78.59044106575273

Residual Variance (σ)

σ is the residual-standard deviation listed under the variance components.

kb07_m
kb07_m.σ
680.0319023667461

Random Effects (θ)

The meaning of θ is a bit less intuitive. In a less complex model (one that only has intercepts for the random effects) or if we suppress the correlations in the formula with zerocorr() then θ describes the relationship between the random effects standard deviation and the standard deviation of the residual term. In our kb07_m example: The residual standard deviation is 680.032. The standard deviation of our first variance component item - (Intercept) is 364.713. Thus our first θ is the relationship: variance component devided by residual standard deviation 364.713 / 680.032 = 0.53631

kb07_m.θ
4-element Vector{Float64}:
  0.5363168233715857
 -0.25981993107379
  0.2653016002105174
  0.4382528181348316

We also can calculate the θ for variance component subj - (Intercept). The residual standard deviation is 680.032. The standard deviation of our variance component subj - (Intercept) is 298.026. Thus, the related θ is the relationship: variance component devided by residual standard deviation 298.026 / 680.032 = 0.438252

kb07_m.θ
4-element Vector{Float64}:
  0.5363168233715857
 -0.25981993107379
  0.2653016002105174
  0.4382528181348316

We can not calculate the θs for variance component item - prec: maintain this way, because it includes the correlation of item - prec: maintain and item - (Intercept). But keep in mind that the relation of item - prec: maintain-variability (252.521) and the residual-variability (680.032) is 252.521 / 680.032 = 0.3713369.

The θ vector is the flattened version of the variance-covariance matrix - a lowertrinangular matrix. The on-diagonal elements are just the standard deviations (the σ's). If all off-diagonal elements are zero, we can use our calculation above. The off-diagonal elements are covariances and correspond to the correlations (the ρ's). If they are unequal to zero, as it is in our kb07-dataset, one way to get the two missing θ-values is to take the values directly from the model we have already fitted.

See the two inner values:

kb07_m.θ
4-element Vector{Float64}:
  0.5363168233715857
 -0.25981993107379
  0.2653016002105174
  0.4382528181348316

Another way is to make use of the create_re() function. Here you have to define the relation of all random effects variabilities to the variability of the residuals, as shown above, and the correlation-matrices.

Let's start by defining the correlation matrix for the item-part.

The diagonal is always 1.0 because everything is perfectly correlated with itself. The elements below the diagonal follow the same form as the Corr. entries in the output of VarCorr(). In our example the correlation of item - prec: maintain and item - (Intercept) is -0.7. The elements above the diagonal are just a mirror image.

re_item_corr = [1.0 -0.7; -0.7 1.0]
2×2 Matrix{Float64}:
  1.0  -0.7
 -0.7   1.0

Now we put together all relations of standard deviations and the correlation-matrix for the item-group. This calculates the covariance factorization which is the theta matrix.

re_item = create_re(0.536, 0.371; corrmat = re_item_corr)
2×2 LowerTriangular{Float64, Matrix{Float64}}:
  0.536    ⋅ 
 -0.2597  0.264947

Work flow for constructing theta vector

Note: Don't be too specific with your values in create_re(). If there are rounding errors, you will get the error-message: PosDefException: matrix is not Hermitian; Cholesky factorization failed.

But you can extract the exact values like shown below:

corr_exact = VarCorr(kb07_m).σρ.item.ρ[1]
σ_residuals_exact = VarCorr(kb07_m).s
σ_1_exact = VarCorr(kb07_m).σρ.item.σ[1] / σ_residuals_exact
σ_2_exact = VarCorr(kb07_m).σρ.item.σ[2] / σ_residuals_exact

re_item_corr = [1.0 corr_exact; corr_exact 1.0]
re_item = create_re(σ_1_exact, σ_2_exact; corrmat = re_item_corr)
2×2 LowerTriangular{Float64, Matrix{Float64}}:
  0.536317   ⋅ 
 -0.25982   0.265302

Let's continue with the subj-part.

Since there the by-subject random effects have only one entry (the intercept), there are no correlations to specify and we can omit the corrmat argument.

Now we put together all relations of standard deviations and the correlation-matrix for the subj-group:

This calculates the covariance factorization which is the theta matrix.

re_subj = create_re(0.438)
1×1 LowerTriangular{Float64, Diagonal{Float64, Vector{Float64}}}:
 0.438

If you want the exact value you can use

σ_residuals_exact = VarCorr(kb07_m).s
σ_3_exact = VarCorr(kb07_m).σρ.subj.σ[1] / σ_residuals_exact
re_subj = create_re(σ_3_exact)
1×1 LowerTriangular{Float64, Diagonal{Float64, Vector{Float64}}}:
 0.4382528181348316

As mentioned above θ is the compact form of these covariance matrices:

kb07_m.θ = vcat( flatlowertri(re_item), flatlowertri(re_subj) )
4-element Vector{Float64}:
  0.5363168233715857
 -0.25981993107379
  0.2653016002105175
  0.4382528181348316

We can install these parameter in the parametricbootstrap()-function or in the model like this:

update!(kb07_m, re_item, re_subj)
Linear mixed model fit by maximum likelihood
 rt_trunc ~ 1 + spkr + prec + load + (1 | subj) + (1 + prec | item)
    logLik   -2 logLik      AIC         AICc        BIC     
 -14331.9251  28663.8501  28681.8501  28681.9513  28731.2548

Variance components:
             Column      Variance  Std.Dev.  Corr.
item     (Intercept)     133015.244 364.713
         prec: maintain   63766.936 252.521 -0.70
subj     (Intercept)      88819.436 298.026
Residual                 462443.388 680.032
 Number of obs: 1789; levels of grouping factors: 32, 56

  Fixed-effects parameters:
──────────────────────────────────────────────────────
                    Coef.  Std. Error      z  Pr(>|z|)
──────────────────────────────────────────────────────
(Intercept)     2181.85       77.4681  28.16    <1e-99
spkr: old         67.879      16.0785   4.22    <1e-04
prec: maintain  -333.791      47.4472  -7.03    <1e-11
load: yes         78.5904     16.0785   4.89    <1e-05
──────────────────────────────────────────────────────

A simple example from scratch

Having this knowledge about the parameters we can now simulate data from scratch

The simdat_crossed() function from MixedModelsSim lets you set up a data frame with a specified experimental design. For now, it only makes fully balanced crossed designs!, but you can generate an unbalanced design by simulating data for the largest cell and deleting extra rows.

Firstly we will set an easy design where subj_n subjects per age group (O or Y) respond to item_n items in each of two conditions (A or B).

Your factors need to be specified separately for between-subject, between-item, and within-subject/item factors using Dict with the name of each factor as the keys and vectors with the names of the levels as values.

We start with the between subject factors:

subj_btwn = Dict("age" => ["O", "Y"])
Dict{String, Vector{String}} with 1 entry:
  "age" => ["O", "Y"]

There are no between-item factors in this design so you can omit it or set it to nothing. Note that if you have factors which are between subject and between item, you need to put them in both dicts.

item_btwn = nothing

Next, we put within-subject/item factors in a dict:

both_win = Dict("condition" => ["A", "B"])
Dict{String, Vector{String}} with 1 entry:
  "condition" => ["A", "B"]

Define subject and item number:

subj_n = 10
item_n = 30
30

Simulate data:

dat = simdat_crossed(subj_n,
                     item_n,
                     subj_btwn = subj_btwn,
                     item_btwn = item_btwn,
                     both_win = both_win);
600-element Vector{NamedTuple{(:subj, :age, :item, :condition, :dv), Tuple{String, String, String, String, Float64}}}:
 (subj = "S01", age = "O", item = "I01", condition = "A", dv = -1.9257742340657893)
 (subj = "S02", age = "Y", item = "I01", condition = "A", dv = 0.4834989932397005)
 (subj = "S03", age = "O", item = "I01", condition = "A", dv = -0.0341260825718261)
 (subj = "S04", age = "Y", item = "I01", condition = "A", dv = 2.1983576066989046)
 (subj = "S05", age = "O", item = "I01", condition = "A", dv = 0.7033196487997226)
 (subj = "S06", age = "Y", item = "I01", condition = "A", dv = 0.9421739433138172)
 (subj = "S07", age = "O", item = "I01", condition = "A", dv = -2.195402078208497)
 (subj = "S08", age = "Y", item = "I01", condition = "A", dv = 1.5652033790201465)
 (subj = "S09", age = "O", item = "I01", condition = "A", dv = 0.13772042663634645)
 (subj = "S10", age = "Y", item = "I01", condition = "A", dv = -0.5706936158999962)
 ⋮
 (subj = "S02", age = "Y", item = "I30", condition = "B", dv = -0.007395767985365643)
 (subj = "S03", age = "O", item = "I30", condition = "B", dv = -0.0966237972492741)
 (subj = "S04", age = "Y", item = "I30", condition = "B", dv = 0.3656375243330442)
 (subj = "S05", age = "O", item = "I30", condition = "B", dv = 1.3891429220948128)
 (subj = "S06", age = "Y", item = "I30", condition = "B", dv = -0.9907756349785514)
 (subj = "S07", age = "O", item = "I30", condition = "B", dv = 0.5713132313975486)
 (subj = "S08", age = "Y", item = "I30", condition = "B", dv = -0.04845808267981579)
 (subj = "S09", age = "O", item = "I30", condition = "B", dv = -0.7222283703055851)
 (subj = "S10", age = "Y", item = "I30", condition = "B", dv = -1.2306265915264)

Have a look:

first(DataFrame(dat),8)

8 rows × 5 columns

subjageitemconditiondv
StringStringStringStringFloat64
1S01OI01A-1.92577
2S02YI01A0.483499
3S03OI01A-0.0341261
4S04YI01A2.19836
5S05OI01A0.70332
6S06YI01A0.942174
7S07OI01A-2.1954
8S08YI01A1.5652

The values we see in the column dv is just random noise.

Set contrasts:

contrasts = Dict(:age => HelmertCoding(),
                 :condition => HelmertCoding());
Dict{Symbol, HelmertCoding} with 2 entries:
  :age       => HelmertCoding(nothing, nothing)
  :condition => HelmertCoding(nothing, nothing)

Define formula:

f1 = @formula dv ~ 1 + age * condition + (1|item) + (1|subj);
FormulaTerm
Response:
  dv(unknown)
Predictors:
  1
  age(unknown)
  condition(unknown)
  (item)->1 | item
  (subj)->1 | subj
  age(unknown) & condition(unknown)

Note that we did not include condition as random slopes for item and subject. This is mainly to keep the example simple and to keep the parameter θ easier to understand (see Section 3 for the explanation of theta).

Fit the model:

m1 = fit(MixedModel, f1, dat, contrasts=contrasts)
Linear mixed model fit by maximum likelihood
 dv ~ 1 + age + condition + age & condition + (1 | item) + (1 | subj)
   logLik   -2 logLik     AIC       AICc        BIC    
  -849.7836  1699.5673  1713.5673  1713.7565  1744.3458

Variance components:
            Column    Variance   Std.Dev.  
item     (Intercept)  0.00000000 0.00000000
subj     (Intercept)  0.00005578 0.00746882
Residual              0.99469320 0.99734307
 Number of obs: 600; levels of grouping factors: 30, 10

  Fixed-effects parameters:
───────────────────────────────────────────────────────────────
                             Coef.  Std. Error      z  Pr(>|z|)
───────────────────────────────────────────────────────────────
(Intercept)             0.0210212    0.0407848   0.52    0.6063
age: Y                 -0.0398237    0.0407848  -0.98    0.3288
condition: B            0.0460315    0.0407164   1.13    0.2582
age: Y & condition: B  -0.00192799   0.0407164  -0.05    0.9622
───────────────────────────────────────────────────────────────

Because the dv is just random noise from N(0,1), there will be basically no subject or item random variance, residual variance will be near 1.0, and the estimates for all effects should be small. Don't worry, we'll specify fixed and random effects directly in parametricbootstrap().

Set random seed for reproducibility:

rng = StableRNG(42);
StableRNGs.LehmerRNG(state=0x00000000000000000000000000000055)

Specify β, σ, and θ, we just made up this parameter:

new_beta = [0., 0.25, 0.25, 0.]
new_sigma = 2.0
new_theta = [1.0, 1.0]
2-element Vector{Float64}:
 1.0
 1.0

Run nsims iterations:

sim1 = parametricbootstrap(rng, nsims, m1,
                        β = new_beta,
                        σ = new_sigma,
                        θ = new_theta,
                        use_threads = false);
MixedModelBootstrap{Float64}(NamedTuple{(:objective, :σ, :β, :se, :θ), Tuple{Float64, Float64, NamedTuple{(Symbol("(Intercept)"), Symbol("age: Y"), Symbol("condition: B"), Symbol("age: Y & condition: B")), NTuple{4, Float64}}, StaticArrays.SVector{4, Float64}, StaticArrays.SVector{2, Float64}}}[(objective = 2692.361031432655, σ = 2.047878422109851, β = ((Intercept) = -1.251469197211205, age: Y = 1.0092232885828774, condition: B = 0.22341873284284625, age: Y & condition: B = 0.0694707523830103), se = [0.7207930394871161, 0.6118632760079705, 0.083604286490418, 0.083604286490418], θ = [1.0190343381659128, 0.9359609240686053]), (objective = 2637.6876454916783, σ = 1.9664599289731421, β = ((Intercept) = 0.5933304076891701, age: Y = 0.16250401624007726, condition: B = 0.06184620640928361, age: Y & condition: B = -0.12242391752570918), se = [0.7344968193104621, 0.6637529491050558, 0.08028039042689748, 0.08028039042689748], θ = [0.8760165516143874, 1.0595496560524502]), (objective = 2584.1440148821357, σ = 1.8805353428190845, β = ((Intercept) = 0.12032043743902278, age: Y = 0.16126672420499952, condition: B = 0.2753148935121751, age: Y & condition: B = 0.06880234926264957), se = [0.5120307889801435, 0.3559679542947235, 0.07677253388627658, 0.07677253388627658], θ = [1.0719842907954271, 0.5845025504712554]), (objective = 2674.0451340881077, σ = 2.006802012678268, β = ((Intercept) = -0.08942325214774771, age: Y = 1.2445124382240857, condition: B = 0.2756056156168277, age: Y & condition: B = 0.0019319947128911019), se = [0.8064189064600313, 0.704189781105721, 0.08192734909753425, 0.08192734909753425], θ = [1.0725546028708817, 1.102112411146808]), (objective = 2687.746337155375, σ = 2.0633485113253083, β = ((Intercept) = -0.6692467458281341, age: Y = 0.510513273956098, condition: B = 0.20596770138857215, age: Y & condition: B = -0.025674164168205473), se = [0.5568391658345431, 0.44337910649084006, 0.08423585023797138, 0.08423585023797138], θ = [0.8942458505540337, 0.6671443589171886]), (objective = 2676.4559583634505, σ = 2.0373920161719585, β = ((Intercept) = 0.4394355388817507, age: Y = 0.17931804380054417, condition: B = 0.3439411667123485, age: Y & condition: B = 0.11467757327199378), se = [0.7436384903738149, 0.6778629978192658, 0.08317618076069253, 0.08317618076069253], θ = [0.8220377108840188, 1.0441744521972427]), (objective = 2662.551612757376, σ = 2.0073631997851695, β = ((Intercept) = -0.5461799640539708, age: Y = -0.020800285880995608, condition: B = 0.3215822875574087, age: Y & condition: B = -0.06355472394332838), se = [0.5695549461503429, 0.4253146251159524, 0.08195025946523654, 0.08195025946523654], θ = [1.0336204265682767, 0.6574595718319739]), (objective = 2636.7044040219603, σ = 1.9547764909530143, β = ((Intercept) = -0.3166835135941909, age: Y = 0.4757914674418475, condition: B = 0.32305158582986454, age: Y & condition: B = 0.0719370653921514), se = [0.7909222263184837, 0.7141321535806116, 0.07980341606705171, 0.07980341606705171], θ = [0.9525595938554858, 1.1480286223430511]), (objective = 2690.4276214574315, σ = 2.0452657137975776, β = ((Intercept) = 0.583469542691405, age: Y = -0.08265650422637824, condition: B = 0.22652487157529652, age: Y & condition: B = -0.08648642349416014), se = [0.7132932076066393, 0.6040754003335782, 0.08349762312022137, 0.08349762312022137], θ = [1.015807163157768, 0.9250229096081493]), (objective = 2708.5482577153466, σ = 2.083599273257543, β = ((Intercept) = 0.35898497948788627, age: Y = -0.15549678679449047, condition: B = 0.2391166597806736, age: Y & condition: B = 0.0006507184175729109), se = [0.6930264753768088, 0.5897517383666084, 0.08506258413191393, 0.08506258413191393], θ = [0.9567958140893679, 0.8857067092059125])  …  (objective = 2615.6137822050496, σ = 1.9279425653180247, β = ((Intercept) = 0.907392200377613, age: Y = -0.16011527089574332, condition: B = 0.2098416432974792, age: Y & condition: B = -0.03935004306592519), se = [0.6789558798573394, 0.5943478375176486, 0.07870792564035982, 0.07870792564035982], θ = [0.9324775018201081, 0.9662837911864878]), (objective = 2595.0147652046144, σ = 1.8770714151239487, β = ((Intercept) = -0.3097896941038077, age: Y = 0.7273109208682076, condition: B = 0.22213316933245922, age: Y & condition: B = -0.003933497669757553), se = [0.8817203524389043, 0.8096754631189734, 0.07663111963029362, 0.07663111963029362], θ = [1.0186008785175913, 1.3579266597421904]), (objective = 2612.0561177595528, σ = 1.9100663931169395, β = ((Intercept) = -0.8295393526201817, age: Y = 0.4469459973230009, condition: B = 0.21426159906915654, age: Y & condition: B = -0.0902989835244413), se = [0.7744047518346746, 0.6893301915030642, 0.07797813396624675, 0.07797813396624675], θ = [1.0119131530426086, 1.1339193876774079]), (objective = 2663.8367335004086, σ = 2.0079216479348725, β = ((Intercept) = -0.3194714891335858, age: Y = -0.7322085965324101, condition: B = 0.2928245603121291, age: Y & condition: B = -0.09158528963374925), se = [0.6057794158645861, 0.4797332425384689, 0.08197305801547944, 0.08197305801547944], θ = [1.0090115652867233, 0.7444208694973428]), (objective = 2610.387180563737, σ = 1.902111756788359, β = ((Intercept) = -0.7131620379606376, age: Y = 1.1127178367866586, condition: B = 0.10302499063161036, age: Y & condition: B = -0.02407355268197516), se = [0.8463774458815531, 0.7665134290817152, 0.07765338729800628, 0.07765338729800628], θ = [1.0334813994220318, 1.267779150103906]), (objective = 2575.2455339616836, σ = 1.852086851237604, β = ((Intercept) = -0.6544944052471983, age: Y = 0.15993095829804993, condition: B = 0.1804519240681131, age: Y & condition: B = 0.07882478295873334), se = [0.703909436663219, 0.6079202026212791, 0.07561112908083509, 0.07561112908083509], θ = [1.049419818013706, 1.0299112433771669]), (objective = 2586.2931471750444, σ = 1.8690790012136018, β = ((Intercept) = 0.7931789621247564, age: Y = -0.8998561878420233, condition: B = 0.24609980924074723, age: Y & condition: B = -0.06510498244263868), se = [0.5901276433245217, 0.42894298275896287, 0.07630483069873575, 0.07630483069873575], θ = [1.1876732237535361, 0.7141496234902017]), (objective = 2612.7028992829746, σ = 1.9315728455258265, β = ((Intercept) = 0.7342671862295446, age: Y = -0.18056171525828635, condition: B = 0.19147311460739172, age: Y & condition: B = 0.06108871204853024), se = [0.741998338980579, 0.6852775469138882, 0.0788561312092338, 0.0788561312092338], θ = [0.806816262800003, 1.1144506558328313]), (objective = 2613.5549060118287, σ = 1.9278629303074943, β = ((Intercept) = -0.21465120723189618, age: Y = 0.5783629182736614, condition: B = 0.18810412516923858, age: Y & condition: B = -0.09316412488902753), se = [0.6383987406007704, 0.5492217877657928, 0.07870467455466881, 0.07870467455466881], θ = [0.924590947263084, 0.8915915669658926]), (objective = 2704.4237783634994, σ = 2.102505843713082, β = ((Intercept) = -0.29587923132339994, age: Y = -0.42128353416205266, condition: B = 0.19992196530671996, age: Y & condition: B = 0.004263394667623395), se = [0.6590087065752558, 0.5971403439493096, 0.08583444163861477, 0.08583444163861477], θ = [0.7262369043834699, 0.8888030304982291])], LowerTriangular{Float64, Matrix{Float64}}[[0.0], [0.007488719437142797]], [[1], [1]], [0.0, 0.0], (item = ("(Intercept)",), subj = ("(Intercept)",)))

Power calculation

ptbl= power_table(sim1)
4-element Vector{NamedTuple{(:coefname, :power), Tuple{String, Float64}}}:
 (coefname = "(Intercept)", power = 0.09)
 (coefname = "age: Y", power = 0.144)
 (coefname = "condition: B", power = 0.87)
 (coefname = "age: Y & condition: B", power = 0.044)

For nicely displaying it, you can use pretty_table:

pretty_table(ptbl)
┌───────────────────────┬─────────┐
│              coefname │   power │
│                String │ Float64 │
├───────────────────────┼─────────┤
│           (Intercept) │    0.09 │
│                age: Y │   0.144 │
│          condition: B │    0.87 │
│ age: Y & condition: B │   0.044 │
└───────────────────────┴─────────┘

Recreate a more complex dataset from scratch and analyze power for specific model parameter but various sample sizes.

Recreate the kb07-dataset from scratch

For full control over all parameters in our kb07 data set we will recreate the design using the method shown above.

Define subject and item number:

subj_n = 56
item_n = 32
32

Define factors in a dict:

subj_btwn = nothing
item_btwn = nothing
both_win = Dict("spkr" => ["old", "new"],
                "prec" => ["maintain", "break"],
                "load" => ["yes", "no"]);
Dict{String, Vector{String}} with 3 entries:
  "spkr" => ["old", "new"]
  "prec" => ["maintain", "break"]
  "load" => ["yes", "no"]

Play with simdat_crossed

subj_btwn = Dict("spkr" => ["old", "new"],
                "prec" => ["maintain", "break"],
                "load" => ["yes", "no"]);
item_btwn = nothing
both_win = nothing;

subj_n = 56
item_n = 32

fake_kb07 = simdat_crossed(subj_n, item_n,
                     subj_btwn = subj_btwn,
                     item_btwn = item_btwn,
                     both_win = both_win);


fake_kb07_df = DataFrame(fake_kb07);

Simulate data:

fake_kb07 = simdat_crossed(subj_n, item_n,
                     subj_btwn = subj_btwn,
                     item_btwn = item_btwn,
                     both_win = both_win);
1792-element Vector{NamedTuple{(:subj, :spkr, :prec, :load, :item, :dv), Tuple{String, String, String, String, String, Float64}}}:
 (subj = "S01", spkr = "old", prec = "maintain", load = "yes", item = "I01", dv = 0.6917653699829488)
 (subj = "S02", spkr = "new", prec = "maintain", load = "yes", item = "I01", dv = -1.3333181995332504)
 (subj = "S03", spkr = "old", prec = "break", load = "yes", item = "I01", dv = -0.7134390925359635)
 (subj = "S04", spkr = "new", prec = "break", load = "yes", item = "I01", dv = -1.8627048468309033)
 (subj = "S05", spkr = "old", prec = "maintain", load = "no", item = "I01", dv = -0.9786246932813167)
 (subj = "S06", spkr = "new", prec = "maintain", load = "no", item = "I01", dv = -0.06604397083993412)
 (subj = "S07", spkr = "old", prec = "break", load = "no", item = "I01", dv = -0.7824314618060416)
 (subj = "S08", spkr = "new", prec = "break", load = "no", item = "I01", dv = 0.09771664348080317)
 (subj = "S09", spkr = "old", prec = "maintain", load = "yes", item = "I01", dv = 0.8900991431682521)
 (subj = "S10", spkr = "new", prec = "maintain", load = "yes", item = "I01", dv = -0.9552707865886649)
 ⋮
 (subj = "S48", spkr = "new", prec = "break", load = "no", item = "I32", dv = -0.501776189989607)
 (subj = "S49", spkr = "old", prec = "maintain", load = "yes", item = "I32", dv = -0.9629841559881522)
 (subj = "S50", spkr = "new", prec = "maintain", load = "yes", item = "I32", dv = -3.297241159819142)
 (subj = "S51", spkr = "old", prec = "break", load = "yes", item = "I32", dv = 0.6015044768817394)
 (subj = "S52", spkr = "new", prec = "break", load = "yes", item = "I32", dv = -0.002121018185968433)
 (subj = "S53", spkr = "old", prec = "maintain", load = "no", item = "I32", dv = -1.6159012315921328)
 (subj = "S54", spkr = "new", prec = "maintain", load = "no", item = "I32", dv = -0.4700845183915087)
 (subj = "S55", spkr = "old", prec = "break", load = "no", item = "I32", dv = 1.2879958397758078)
 (subj = "S56", spkr = "new", prec = "break", load = "no", item = "I32", dv = 0.05728319336087551)

Make a dataframe:

fake_kb07_df = DataFrame(fake_kb07);

Have a look:

first(fake_kb07_df,8)

8 rows × 6 columns

subjspkrprecloaditemdv
StringStringStringStringStringFloat64
1S01oldmaintainyesI010.691765
2S02newmaintainyesI01-1.33332
3S03oldbreakyesI01-0.713439
4S04newbreakyesI01-1.8627
5S05oldmaintainnoI01-0.978625
6S06newmaintainnoI01-0.066044
7S07oldbreaknoI01-0.782431
8S08newbreaknoI010.0977166

The function simdat_crossed generates a balanced fully crossed design. Unfortunately, our original design is not fully crossed. Every subject saw an image only once, thus in one of eight possible conditions. To simulate that we only keep one of every eight lines.

We sort the dataframe to enable easier selection

fake_kb07_df = sort(fake_kb07_df, [:subj, :item, :load, :prec, :spkr])

In order to select only the relevant rows of the data set we define an index which represets a random choice of one of every eight rows. First we generate a vector idx which represents which row to keep in each set of 8.

len = convert(Int64,(length(fake_kb07)/8))
idx = rand(rng, collect(1:8) , len)
224-element Vector{Int64}:
 2
 7
 3
 1
 6
 5
 5
 3
 5
 8
 ⋮
 1
 7
 1
 4
 8
 3
 7
 5
 3

Then we create an array A, of the same length that is populated multiples of the number 8. Added together A and idx give the indexes of one row from each set of 8s.

A = repeat([8], inner=len-1)
A = append!( [0], A )
A = cumsum(A)
idx = idx+A
224-element Vector{Int64}:
    2
   15
   19
   25
   38
   45
   53
   59
   69
   80
    ⋮
 1721
 1735
 1737
 1748
 1760
 1763
 1775
 1781
 1787

Reduce the fully crossed design to the original experimental design:

fake_kb07_df= fake_kb07_df[idx, :]
rename!(fake_kb07_df, :dv => :rt_trunc)

Now we can use the simulated data in the same way as above.

Set contrasts:

contrasts = Dict(:spkr => HelmertCoding(),
                 :prec => HelmertCoding(),
                 :load => HelmertCoding());
Dict{Symbol, HelmertCoding} with 3 entries:
  :spkr => HelmertCoding(nothing, nothing)
  :load => HelmertCoding(nothing, nothing)
  :prec => HelmertCoding(nothing, nothing)

Define formula, same as above:

kb07_f = @formula( rt_trunc ~ 1 + spkr+prec+load + (1|subj) + (1+prec|item) );
FormulaTerm
Response:
  rt_trunc(unknown)
Predictors:
  1
  spkr(unknown)
  prec(unknown)
  load(unknown)
  (subj)->1 | subj
  (prec,item)->(1 + prec) | item

Fit the model:

fake_kb07_m = fit(MixedModel, kb07_f, fake_kb07_df, contrasts=contrasts)
Linear mixed model fit by maximum likelihood
 rt_trunc ~ 1 + spkr + prec + load + (1 | subj) + (1 + prec | item)
   logLik   -2 logLik     AIC       AICc        BIC    
  -306.6446   613.2891   631.2891   632.1303   661.9939

Variance components:
             Column      Variance   Std.Dev.    Corr.
item     (Intercept)     0.00000000 0.00000000
         prec: maintain  0.00000000 0.00000005   NaN
subj     (Intercept)     0.00000000 0.00000000
Residual                 0.90485629 0.95123935
 Number of obs: 224; levels of grouping factors: 32, 56

  Fixed-effects parameters:
────────────────────────────────────────────────────────
                      Coef.  Std. Error      z  Pr(>|z|)
────────────────────────────────────────────────────────
(Intercept)     -0.0192356    0.0635574  -0.30    0.7622
spkr: old       -0.0443656    0.0635574  -0.70    0.4852
prec: maintain  -0.107652     0.0635574  -1.69    0.0903
load: yes       -0.00894328   0.0635574  -0.14    0.8881
────────────────────────────────────────────────────────

Set random seed for reproducibility:

rng = StableRNG(42);
StableRNGs.LehmerRNG(state=0x00000000000000000000000000000055)

Then, again, we specify β, σ, and θ. Here we use the values that we found in the model of the existing dataset:

#beta
new_beta = [2181.85, 67.879, -333.791, 78.5904] #manual
new_beta = kb07_m.β #grab from existing model

#sigma
new_sigma = 680.032 #manual
new_sigma = kb07_m.σ #grab from existing model

#theta
re_item_corr = [1.0 -0.7; -0.7 1.0]
re_item = create_re(0.536, 0.371; corrmat = re_item_corr)
re_subj = create_re(0.438)
new_theta = vcat( flatlowertri(re_item), flatlowertri(re_subj) )  #manual
new_theta = kb07_m.θ #grab from existing model
4-element Vector{Float64}:
  0.5363168233715857
 -0.25981993107379
  0.2653016002105175
  0.4382528181348316

Run nsims iterations:

fake_kb07_sim = parametricbootstrap(rng, nsims, fake_kb07_m,
                        β = new_beta,
                        σ = new_sigma,
                        θ = new_theta,
                        use_threads = false);
MixedModelBootstrap{Float64}(NamedTuple{(:objective, :σ, :β, :se, :θ), Tuple{Float64, Float64, NamedTuple{(Symbol("(Intercept)"), Symbol("spkr: old"), Symbol("prec: maintain"), Symbol("load: yes")), NTuple{4, Float64}}, StaticArrays.SVector{4, Float64}, StaticArrays.SVector{4, Float64}}}[(objective = 3636.636285794904, σ = 737.8857380598588, β = ((Intercept) = 2111.405819736409, spkr: old = 57.18418731355799, prec: maintain = -402.73230527938114, load: yes = 84.82884142923483), se = [86.4198781520241, 51.87921155401374, 69.75154228584405, 51.40482457490105], θ = [0.5174041383672098, -0.2925406860008569, 0.17923198570437404, 0.0]), (objective = 3627.1505517589667, σ = 745.9718489594835, β = ((Intercept) = 2066.3870381457955, spkr: old = 86.88158595416425, prec: maintain = -250.22254269615115, load: yes = 72.75111671504818), se = [79.33403313427756, 52.033960371049915, 54.301694889072564, 51.28062361875329], θ = [0.44849206360424093, -0.11690227391539593, 0.0, 0.0]), (objective = 3629.6460114726697, σ = 676.0321094560071, β = ((Intercept) = 2165.875554932161, spkr: old = 71.68243712507996, prec: maintain = -327.63895210119534, load: yes = 65.6877391839798), se = [101.76903256947129, 60.77137600620894, 73.84762929150146, 59.9807546356846], θ = [0.6719766827577028, -0.3395309365740736, 0.0, 0.4138892747535174]), (objective = 3614.962672019733, σ = 638.6096083714981, β = ((Intercept) = 2122.375329875919, spkr: old = 66.76109839506924, prec: maintain = -356.3902675276288, load: yes = 74.45009837055352), se = [99.03897158454768, 58.02257816878528, 77.76739862899534, 57.51874720813993], θ = [0.6970866218145167, -0.3607065382881416, 0.25397125963434714, 0.41691641982030886]), (objective = 3568.466169872345, σ = 593.8043715632481, β = ((Intercept) = 2255.2564516772013, spkr: old = 99.42944585161713, prec: maintain = -291.3938233220894, load: yes = 30.325128563838053), se = [79.96769522940897, 55.39433640750964, 66.6455581187081, 54.91397258706072], θ = [0.5368629527445188, -0.3398942528928302, 0.0, 0.45956279630211216]), (objective = 3639.5730166276935, σ = 635.0624388220865, β = ((Intercept) = 2183.5111776046515, spkr: old = 124.98699958908655, prec: maintain = -402.84445590463395, load: yes = 113.28688147963294), se = [100.46915183651359, 69.9907631548012, 92.63860588673434, 69.80802906969369], θ = [0.6257147567803769, -0.2714511643484786, 0.4458236488607296, 0.6201163497083998]), (objective = 3602.2457032191055, σ = 689.4423856680926, β = ((Intercept) = 2154.1048936678176, spkr: old = 157.04079132953518, prec: maintain = -309.845595286944, load: yes = 90.55326844499162), se = [71.26612240452253, 48.38023238131034, 64.83853476435611, 48.24976374663757], θ = [0.4169000975175746, -0.22364897189251204, 0.2568936949137387, 0.0]), (objective = 3671.9305850362794, σ = 709.1356754294544, β = ((Intercept) = 2134.180542502668, spkr: old = 33.1078583027862, prec: maintain = -218.13396655787793, load: yes = -17.898191135864092), se = [100.91510693718908, 76.44361820658908, 89.49491061658799, 76.31637867682494], θ = [0.5155653960285771, 0.024778035135990657, 0.357897114828732, 0.60075613424156]), (objective = 3611.963173433437, σ = 657.8127452056431, β = ((Intercept) = 2196.381236920855, spkr: old = 112.89298941187955, prec: maintain = -161.06181377421623, load: yes = 77.14857261946072), se = [93.03411057113848, 53.147116659937105, 73.05144271790168, 52.676819191055586], θ = [0.6432000429330725, -0.3258771068995159, 0.25564781959147387, 0.2826006560164839]), (objective = 3613.7529812886833, σ = 622.5552292324423, β = ((Intercept) = 2163.5605835363153, spkr: old = 103.1653880549227, prec: maintain = -312.5310576454006, load: yes = 91.1287774802018), se = [90.50357525207805, 67.59676754111007, 85.26010716472555, 67.36327453012777], θ = [0.5306545724380344, -0.3943820657479962, 0.22683125416919236, 0.6155402429408227])  …  (objective = 3598.4472644073526, σ = 703.7035910363653, β = ((Intercept) = 2276.906835381216, spkr: old = 202.46791609440405, prec: maintain = -412.1684722850962, load: yes = 56.07837575173253), se = [66.07157151680258, 49.856242039746604, 55.419385838291184, 49.51088781067061], θ = [0.34046660928063166, -0.15795242721395536, 0.09905896495589596, 0.11689446078407409]), (objective = 3666.818178125255, σ = 765.9819491852217, β = ((Intercept) = 2269.3215959487775, spkr: old = 56.50364838732531, prec: maintain = -432.26420967546085, load: yes = 94.99542419507878), se = [100.93038709540515, 63.111736190967655, 67.9746008145919, 62.197316049125234], θ = [0.5747922208019896, -0.18295230411142666, 0.0, 0.318248330144641]), (objective = 3641.078220540293, σ = 667.4564160755705, β = ((Intercept) = 2148.4902549511253, spkr: old = 49.69290987246406, prec: maintain = -367.3469405023088, load: yes = 116.07243788468175), se = [98.51229939254516, 65.9928713167234, 81.62281993413895, 65.65636083213344], θ = [0.6079364135221806, -0.2011680209464133, 0.33603426827640026, 0.5061696279904735]), (objective = 3644.5874918881736, σ = 748.4732756426573, β = ((Intercept) = 2271.6729862179886, spkr: old = 86.2810586609657, prec: maintain = -460.2248684740786, load: yes = 140.6954136740148), se = [80.78720445916493, 63.38602655371095, 64.17316886414869, 62.84529972061185], θ = [0.37380004408019535, -0.07830341989529281, 0.0, 0.3631903469883493]), (objective = 3639.4853122711643, σ = 728.0350055599372, β = ((Intercept) = 2242.4196711744526, spkr: old = 99.75269349229724, prec: maintain = -260.1289768856411, load: yes = 104.56434913516995), se = [98.65476894801031, 54.40551471881057, 68.181547735203, 53.539322805000424], θ = [0.6288624235658962, -0.3081038047455251, 0.0, 0.1867134395972392]), (objective = 3629.3558153012914, σ = 620.0506774392759, β = ((Intercept) = 2106.038910607525, spkr: old = 30.940425109368295, prec: maintain = -419.36807326180025, load: yes = 2.2691772244695234), se = [113.69290493638893, 61.659322805152705, 86.53513278614244, 61.23244430396673], θ = [0.8564432856064275, -0.3276687546092365, 0.4191899343568653, 0.5025205870661192]), (objective = 3651.921613995107, σ = 706.2506302237504, β = ((Intercept) = 2320.700961608762, spkr: old = -53.710769283629794, prec: maintain = -345.09126148818865, load: yes = 156.16670862694687), se = [90.70182153084089, 67.02680732177164, 90.55723313002734, 66.88375288058437], θ = [0.4718180277909702, -0.43895789352026515, 0.1681803362470318, 0.478991892983086]), (objective = 3650.3777163117425, σ = 723.6150213823325, β = ((Intercept) = 2135.205167372525, spkr: old = 14.395785830337047, prec: maintain = -334.03761632261853, load: yes = 149.8271329788226), se = [92.40663025767286, 61.71187321867084, 85.28047137183803, 61.33104245748815], θ = [0.520841031158126, -0.44255596129033714, 0.0, 0.36613783707936953]), (objective = 3645.3306146794293, σ = 688.5727614231807, β = ((Intercept) = 2155.239867156788, spkr: old = 39.581137888475425, prec: maintain = -302.28187959268064, load: yes = -43.713772168819645), se = [109.23738873517202, 66.54543903822713, 72.38700249936318, 65.64131427382078], θ = [0.7044426569258122, -0.22998110664143784, 0.0, 0.487997708024065]), (objective = 3624.994021244288, σ = 652.0841303505555, β = ((Intercept) = 2378.4837764690083, spkr: old = 125.01362481968539, prec: maintain = -431.245196927933, load: yes = 72.71395621223759), se = [90.7550542835542, 58.794066427331394, 81.80757595975109, 58.63056513180373], θ = [0.584750422175592, -0.2225281877886704, 0.4209198371906079, 0.40441880554354054])], LowerTriangular{Float64, Matrix{Float64}}[[0.0 0.0; 5.681421022013248e-8 0.0], [0.0]], [[1, 2, 4], [1]], [0.0, -Inf, 0.0, 0.0], (item = ("(Intercept)", "prec: maintain"), subj = ("(Intercept)",)))

Power calculation

power_table(fake_kb07_sim)
4-element Vector{NamedTuple{(:coefname, :power), Tuple{String, Float64}}}:
 (coefname = "(Intercept)", power = 1.0)
 (coefname = "spkr: old", power = 0.21)
 (coefname = "prec: maintain", power = 0.994)
 (coefname = "load: yes", power = 0.24)

Compare to the powertable from the existing data:

power_table(kb07_sim)
4-element Vector{NamedTuple{(:coefname, :power), Tuple{String, Float64}}}:
 (coefname = "(Intercept)", power = 1.0)
 (coefname = "spkr: old", power = 0.99)
 (coefname = "prec: maintain", power = 1.0)
 (coefname = "load: yes", power = 1.0)

We have successfully recreated the power simulation of an existing dataset from scratch. This has the advantage, that we now can iterate over different numbers of subjects and items.

Loop over subject and item sizes

When designing a study, you may be interested in trying various numbers of subjects and items to see how that affects the power of your study. To do this you can use a loop to run simulations over a range of values for any parameter.

We first define every fixed things outside the loop (same as above):

Define factors in a dict:

subj_btwn = nothing
item_btwn = nothing
both_win = Dict("spkr" => ["old", "new"],
                "prec" => ["maintain", "break"],
                "load" => ["yes", "no"]);
Dict{String, Vector{String}} with 3 entries:
  "spkr" => ["old", "new"]
  "prec" => ["maintain", "break"]
  "load" => ["yes", "no"]

Set contrasts:

contrasts = Dict(:spkr => HelmertCoding(),
                 :prec => HelmertCoding(),
                 :load => HelmertCoding());
Dict{Symbol, HelmertCoding} with 3 entries:
  :spkr => HelmertCoding(nothing, nothing)
  :load => HelmertCoding(nothing, nothing)
  :prec => HelmertCoding(nothing, nothing)

Set random seed for reproducibility:

rng = StableRNG(42);
StableRNGs.LehmerRNG(state=0x00000000000000000000000000000055)

Define formula:

kb07_f = @formula( rt_trunc ~ 1 + spkr+prec+load + (1|subj) + (1+prec|item) );
FormulaTerm
Response:
  rt_trunc(unknown)
Predictors:
  1
  spkr(unknown)
  prec(unknown)
  load(unknown)
  (subj)->1 | subj
  (prec,item)->(1 + prec) | item

Specify β, σ, and θ:

#beta
new_beta = [2181.85, 67.879, -333.791, 78.5904]
new_beta = kb07_m.β

#sigma
new_sigma = 680.032
new_sigma = kb07_m.σ

#theta
re_item_corr = [1.0 -0.7; -0.7 1.0]
re_item = create_re(0.536, 0.371; corrmat = re_item_corr)
re_subj = create_re(0.438)
new_theta = vcat( flatlowertri(re_item), flatlowertri(re_subj) )
new_theta = kb07_m.θ
4-element Vector{Float64}:
  0.5363168233715857
 -0.25981993107379
  0.2653016002105175
  0.4382528181348316

Then we define the variables that out loop will iterate over

Define subject and item numbers as arrays:

sub_ns = [20, 30, 40];
item_ns = [16, 24, 32];
3-element Vector{Int64}:
 16
 24
 32

Make an empty dataframe:

d = DataFrame();

0 rows × 0 columns

Run the loop:

for subj_n in sub_ns
    for item_n in item_ns

    # Make fully crossed data:
    fake_kb07 = simdat_crossed(subj_n, item_n,
                     subj_btwn = subj_btwn,
                     item_btwn = item_btwn,
                     both_win = both_win);
    fake_kb07_df = DataFrame(fake_kb07)

    # Reduce the fully crossed design to the original experimental design:
    fake_kb07_df = sort(fake_kb07_df, [:subj, :item, :load, :prec, :spkr])
    len = convert(Int64,(length(fake_kb07)/8))
    idx = rand(rng, collect(1:8) , len)
    A = repeat([8], inner=len-1)
    A = append!( [0], A )
    A = cumsum(A)
    idx = idx+A
    fake_kb07_df= fake_kb07_df[idx, :]
    rename!(fake_kb07_df, :dv => :rt_trunc)

    # Fit the model:
    fake_kb07_m = fit(MixedModel, kb07_f, fake_kb07_df, contrasts=contrasts);

    # Run nsims iterations:
    fake_kb07_sim = parametricbootstrap(rng, nsims, fake_kb07_m,
                        β = new_beta,
                        σ = new_sigma,
                        θ = new_theta,
                        use_threads = false);

    # Power calculation
    ptbl = power_table(fake_kb07_sim)
    ptdf = DataFrame(ptbl)
    ptdf[!, :item_n] .= item_n
    ptdf[!, :subj_n] .= subj_n
    append!(d, ptdf)

    end
end
┌ Warning: Assignment to `fake_kb07` in soft scope is ambiguous because a global variable by the same name exists: `fake_kb07` will be treated as a new local. Disambiguate by using `local fake_kb07` to suppress this warning or `global fake_kb07` to assign to the existing global variable.
└ @ none:5
┌ Warning: Assignment to `fake_kb07_df` in soft scope is ambiguous because a global variable by the same name exists: `fake_kb07_df` will be treated as a new local. Disambiguate by using `local fake_kb07_df` to suppress this warning or `global fake_kb07_df` to assign to the existing global variable.
└ @ none:9
┌ Warning: Assignment to `len` in soft scope is ambiguous because a global variable by the same name exists: `len` will be treated as a new local. Disambiguate by using `local len` to suppress this warning or `global len` to assign to the existing global variable.
└ @ none:13
┌ Warning: Assignment to `idx` in soft scope is ambiguous because a global variable by the same name exists: `idx` will be treated as a new local. Disambiguate by using `local idx` to suppress this warning or `global idx` to assign to the existing global variable.
└ @ none:14
┌ Warning: Assignment to `A` in soft scope is ambiguous because a global variable by the same name exists: `A` will be treated as a new local. Disambiguate by using `local A` to suppress this warning or `global A` to assign to the existing global variable.
└ @ none:15
┌ Warning: Assignment to `fake_kb07_m` in soft scope is ambiguous because a global variable by the same name exists: `fake_kb07_m` will be treated as a new local. Disambiguate by using `local fake_kb07_m` to suppress this warning or `global fake_kb07_m` to assign to the existing global variable.
└ @ none:23
┌ Warning: Assignment to `fake_kb07_sim` in soft scope is ambiguous because a global variable by the same name exists: `fake_kb07_sim` will be treated as a new local. Disambiguate by using `local fake_kb07_sim` to suppress this warning or `global fake_kb07_sim` to assign to the existing global variable.
└ @ none:26
┌ Warning: Assignment to `ptbl` in soft scope is ambiguous because a global variable by the same name exists: `ptbl` will be treated as a new local. Disambiguate by using `local ptbl` to suppress this warning or `global ptbl` to assign to the existing global variable.
└ @ none:33

Progress:  11%|████▍                                    |  ETA: 0:00:01
Progress:  22%|█████████                                |  ETA: 0:00:01
Progress:  33%|█████████████▌                           |  ETA: 0:00:01
Progress:  44%|█████████████████▉                       |  ETA: 0:00:01
Progress:  54%|██████████████████████▎                  |  ETA: 0:00:00
Progress:  65%|██████████████████████████▌              |  ETA: 0:00:00
Progress:  75%|██████████████████████████████▉          |  ETA: 0:00:00
Progress:  84%|██████████████████████████████████▎      |  ETA: 0:00:00
Progress:  94%|██████████████████████████████████████▊  |  ETA: 0:00:00
Progress: 100%|█████████████████████████████████████████| Time: 0:00:00

Progress:   8%|███▍                                     |  ETA: 0:00:01
Progress:  17%|██████▊                                  |  ETA: 0:00:01
Progress:  23%|█████████▍                               |  ETA: 0:00:01
Progress:  31%|████████████▉                            |  ETA: 0:00:01
Progress:  39%|████████████████▏                        |  ETA: 0:00:01
Progress:  48%|███████████████████▋                     |  ETA: 0:00:01
Progress:  56%|███████████████████████                  |  ETA: 0:00:01
Progress:  64%|██████████████████████████▎              |  ETA: 0:00:00
Progress:  72%|█████████████████████████████▋           |  ETA: 0:00:00
Progress:  80%|████████████████████████████████▉        |  ETA: 0:00:00
Progress:  88%|████████████████████████████████████▎    |  ETA: 0:00:00
Progress:  97%|███████████████████████████████████████▋ |  ETA: 0:00:00
Progress: 100%|█████████████████████████████████████████| Time: 0:00:01

Progress:   7%|██▊                                      |  ETA: 0:00:01
Progress:  13%|█████▍                                   |  ETA: 0:00:01
Progress:  20%|████████▎                                |  ETA: 0:00:01
Progress:  26%|██████████▉                              |  ETA: 0:00:01
Progress:  33%|█████████████▋                           |  ETA: 0:00:01
Progress:  40%|████████████████▍                        |  ETA: 0:00:01
Progress:  47%|███████████████████▎                     |  ETA: 0:00:01
Progress:  53%|█████████████████████▊                   |  ETA: 0:00:01
Progress:  60%|████████████████████████▌                |  ETA: 0:00:01
Progress:  66%|███████████████████████████▏             |  ETA: 0:00:01
Progress:  73%|█████████████████████████████▉           |  ETA: 0:00:00
Progress:  80%|████████████████████████████████▋        |  ETA: 0:00:00
Progress:  86%|███████████████████████████████████▍     |  ETA: 0:00:00
Progress:  93%|██████████████████████████████████████   |  ETA: 0:00:00
Progress:  99%|████████████████████████████████████████▊|  ETA: 0:00:00
Progress: 100%|█████████████████████████████████████████| Time: 0:00:01

Progress:   8%|███▌                                     |  ETA: 0:00:01
Progress:  17%|██████▊                                  |  ETA: 0:00:01
Progress:  25%|██████████▎                              |  ETA: 0:00:01
Progress:  34%|█████████████▊                           |  ETA: 0:00:01
Progress:  42%|█████████████████▎                       |  ETA: 0:00:01
Progress:  51%|████████████████████▊                    |  ETA: 0:00:01
Progress:  59%|████████████████████████▏                |  ETA: 0:00:00
Progress:  67%|███████████████████████████▋             |  ETA: 0:00:00
Progress:  76%|███████████████████████████████          |  ETA: 0:00:00
Progress:  84%|██████████████████████████████████▌      |  ETA: 0:00:00
Progress:  93%|██████████████████████████████████████   |  ETA: 0:00:00
Progress: 100%|█████████████████████████████████████████| Time: 0:00:01

Progress:   6%|██▋                                      |  ETA: 0:00:01
Progress:  13%|█████▎                                   |  ETA: 0:00:01
Progress:  19%|████████                                 |  ETA: 0:00:01
Progress:  26%|██████████▋                              |  ETA: 0:00:01
Progress:  33%|█████████████▍                           |  ETA: 0:00:01
Progress:  39%|████████████████▏                        |  ETA: 0:00:01
Progress:  46%|██████████████████▉                      |  ETA: 0:00:01
Progress:  53%|█████████████████████▋                   |  ETA: 0:00:01
Progress:  60%|████████████████████████▍                |  ETA: 0:00:01
Progress:  66%|███████████████████████████▏             |  ETA: 0:00:01
Progress:  73%|█████████████████████████████▊           |  ETA: 0:00:00
Progress:  77%|███████████████████████████████▋         |  ETA: 0:00:00
Progress:  84%|██████████████████████████████████▎      |  ETA: 0:00:00
Progress:  90%|█████████████████████████████████████    |  ETA: 0:00:00
Progress:  97%|███████████████████████████████████████▋ |  ETA: 0:00:00
Progress: 100%|█████████████████████████████████████████| Time: 0:00:01

Progress:   5%|██▏                                      |  ETA: 0:00:02
Progress:  11%|████▍                                    |  ETA: 0:00:02
Progress:  16%|██████▌                                  |  ETA: 0:00:02
Progress:  21%|████████▊                                |  ETA: 0:00:01
Progress:  27%|███████████                              |  ETA: 0:00:01
Progress:  32%|█████████████▎                           |  ETA: 0:00:01
Progress:  38%|███████████████▍                         |  ETA: 0:00:01
Progress:  43%|█████████████████▋                       |  ETA: 0:00:01
Progress:  49%|███████████████████▉                     |  ETA: 0:00:01
Progress:  54%|██████████████████████▏                  |  ETA: 0:00:01
Progress:  60%|████████████████████████▍                |  ETA: 0:00:01
Progress:  65%|██████████████████████████▋              |  ETA: 0:00:01
Progress:  70%|████████████████████████████▊            |  ETA: 0:00:01
Progress:  76%|███████████████████████████████▏         |  ETA: 0:00:00
Progress:  81%|█████████████████████████████████▎       |  ETA: 0:00:00
Progress:  87%|███████████████████████████████████▌     |  ETA: 0:00:00
Progress:  92%|█████████████████████████████████████▊   |  ETA: 0:00:00
Progress:  97%|███████████████████████████████████████▉ |  ETA: 0:00:00
Progress: 100%|█████████████████████████████████████████| Time: 0:00:01

Progress:  10%|███▉                                     |  ETA: 0:00:01
Progress:  19%|████████                                 |  ETA: 0:00:01
Progress:  28%|███████████▋                             |  ETA: 0:00:01
Progress:  38%|███████████████▋                         |  ETA: 0:00:01
Progress:  48%|███████████████████▌                     |  ETA: 0:00:01
Progress:  57%|███████████████████████▍                 |  ETA: 0:00:00
Progress:  66%|███████████████████████████▎             |  ETA: 0:00:00
Progress:  76%|███████████████████████████████▏         |  ETA: 0:00:00
Progress:  85%|███████████████████████████████████      |  ETA: 0:00:00
Progress:  95%|███████████████████████████████████████  |  ETA: 0:00:00
Progress: 100%|█████████████████████████████████████████| Time: 0:00:01

Progress:   5%|█▉                                       |  ETA: 0:00:02
Progress:   9%|███▌                                     |  ETA: 0:00:02
Progress:  13%|█████▎                                   |  ETA: 0:00:02
Progress:  17%|███████▏                                 |  ETA: 0:00:02
Progress:  22%|█████████                                |  ETA: 0:00:02
Progress:  26%|██████████▊                              |  ETA: 0:00:02
Progress:  31%|████████████▌                            |  ETA: 0:00:02
Progress:  35%|██████████████▍                          |  ETA: 0:00:02
Progress:  40%|████████████████▎                        |  ETA: 0:00:01
Progress:  44%|██████████████████▏                      |  ETA: 0:00:01
Progress:  48%|███████████████████▉                     |  ETA: 0:00:01
Progress:  53%|█████████████████████▋                   |  ETA: 0:00:01
Progress:  57%|███████████████████████▌                 |  ETA: 0:00:01
Progress:  62%|█████████████████████████▎               |  ETA: 0:00:01
Progress:  66%|███████████████████████████              |  ETA: 0:00:01
Progress:  70%|████████████████████████████▊            |  ETA: 0:00:01
Progress:  74%|██████████████████████████████▌          |  ETA: 0:00:01
Progress:  79%|████████████████████████████████▎        |  ETA: 0:00:01
Progress:  83%|██████████████████████████████████       |  ETA: 0:00:00
Progress:  87%|███████████████████████████████████▊     |  ETA: 0:00:00
Progress:  91%|█████████████████████████████████████▌   |  ETA: 0:00:00
Progress:  95%|███████████████████████████████████████▏ |  ETA: 0:00:00
Progress: 100%|████████████████████████████████████████▉|  ETA: 0:00:00
Progress: 100%|█████████████████████████████████████████| Time: 0:00:02

Progress:   4%|█▌                                       |  ETA: 0:00:03
Progress:   7%|███                                      |  ETA: 0:00:03
Progress:  11%|████▌                                    |  ETA: 0:00:02
Progress:  14%|█████▉                                   |  ETA: 0:00:02
Progress:  18%|███████▎                                 |  ETA: 0:00:02
Progress:  21%|████████▊                                |  ETA: 0:00:02
Progress:  25%|██████████▍                              |  ETA: 0:00:02
Progress:  29%|███████████▉                             |  ETA: 0:00:02
Progress:  33%|█████████████▍                           |  ETA: 0:00:02
Progress:  36%|██████████████▉                          |  ETA: 0:00:02
Progress:  40%|████████████████▍                        |  ETA: 0:00:02
Progress:  44%|█████████████████▉                       |  ETA: 0:00:02
Progress:  48%|███████████████████▌                     |  ETA: 0:00:01
Progress:  51%|█████████████████████▏                   |  ETA: 0:00:01
Progress:  55%|██████████████████████▋                  |  ETA: 0:00:01
Progress:  59%|████████████████████████▏                |  ETA: 0:00:01
Progress:  63%|█████████████████████████▋               |  ETA: 0:00:01
Progress:  66%|███████████████████████████              |  ETA: 0:00:01
Progress:  69%|████████████████████████████▌            |  ETA: 0:00:01
Progress:  73%|█████████████████████████████▉           |  ETA: 0:00:01
Progress:  77%|███████████████████████████████▌         |  ETA: 0:00:01
Progress:  81%|█████████████████████████████████        |  ETA: 0:00:01
Progress:  84%|██████████████████████████████████▋      |  ETA: 0:00:00
Progress:  88%|████████████████████████████████████▏    |  ETA: 0:00:00
Progress:  92%|█████████████████████████████████████▋   |  ETA: 0:00:00
Progress:  95%|███████████████████████████████████████▏ |  ETA: 0:00:00
Progress:  99%|████████████████████████████████████████▋|  ETA: 0:00:00
Progress: 100%|█████████████████████████████████████████| Time: 0:00:02

Our dataframe d now contains the power information for each combination of subjects and items.

print(d)
36×4 DataFrame
 Row │ coefname        power    item_n  subj_n
     │ String          Float64  Int64   Int64
─────┼─────────────────────────────────────────
   1 │ (Intercept)       1.0        16      20
   2 │ spkr: old         0.458      16      20
   3 │ prec: maintain    0.994      16      20
   4 │ load: yes         0.538      16      20
   5 │ (Intercept)       1.0        24      20
   6 │ spkr: old         0.548      24      20
   7 │ prec: maintain    1.0        24      20
   8 │ load: yes         0.696      24      20
   9 │ (Intercept)       1.0        32      20
  10 │ spkr: old         0.738      32      20
  11 │ prec: maintain    1.0        32      20
  12 │ load: yes         0.816      32      20
  13 │ (Intercept)       1.0        16      30
  14 │ spkr: old         0.564      16      30
  15 │ prec: maintain    1.0        16      30
  16 │ load: yes         0.708      16      30
  17 │ (Intercept)       1.0        24      30
  18 │ spkr: old         0.722      24      30
  19 │ prec: maintain    1.0        24      30
  20 │ load: yes         0.862      24      30
  21 │ (Intercept)       1.0        32      30
  22 │ spkr: old         0.834      32      30
  23 │ prec: maintain    1.0        32      30
  24 │ load: yes         0.938      32      30
  25 │ (Intercept)       1.0        16      40
  26 │ spkr: old         0.666      16      40
  27 │ prec: maintain    0.962      16      40
  28 │ load: yes         0.786      16      40
  29 │ (Intercept)       1.0        24      40
  30 │ spkr: old         0.86       24      40
  31 │ prec: maintain    1.0        24      40
  32 │ load: yes         0.92       24      40
  33 │ (Intercept)       1.0        32      40
  34 │ spkr: old         0.94       32      40
  35 │ prec: maintain    1.0        32      40
  36 │ load: yes         0.984      32      40

Lastly we plot our results:

categorical!(d, :item_n)

plot(d, x="subj_n", y="power",xgroup= "coefname",color="item_n", Geom.subplot_grid(Geom.point, Geom.line), Guide.xlabel("Number of subjects by parameter"), Guide.ylabel("Power"))
Number of subjects by parameter 16 24 32 item_n load: yes prec: maintain spkr: old (Intercept) -5 0 5 10 15 20 25 30 35 40 45 50 55 60 65 0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 42 44 46 48 50 52 54 56 58 60 0 20 40 60 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 -5 0 5 10 15 20 25 30 35 40 45 50 55 60 65 0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 42 44 46 48 50 52 54 56 58 60 0 20 40 60 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 -5 0 5 10 15 20 25 30 35 40 45 50 55 60 65 0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 42 44 46 48 50 52 54 56 58 60 0 20 40 60 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 -5 0 5 10 15 20 25 30 35 40 45 50 55 60 65 0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 42 44 46 48 50 52 54 56 58 60 0 20 40 60 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 -0.25 -0.20 -0.15 -0.10 -0.05 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00 1.05 1.10 1.15 1.20 1.25 1.30 1.35 1.40 1.45 1.50 1.55 1.60 -1 0 1 2 -0.20 -0.18 -0.16 -0.14 -0.12 -0.10 -0.08 -0.06 -0.04 -0.02 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.16 0.18 0.20 0.22 0.24 0.26 0.28 0.30 0.32 0.34 0.36 0.38 0.40 0.42 0.44 0.46 0.48 0.50 0.52 0.54 0.56 0.58 0.60 0.62 0.64 0.66 0.68 0.70 0.72 0.74 0.76 0.78 0.80 0.82 0.84 0.86 0.88 0.90 0.92 0.94 0.96 0.98 1.00 1.02 1.04 1.06 1.08 1.10 1.12 1.14 1.16 1.18 1.20 1.22 1.24 1.26 1.28 1.30 1.32 1.34 1.36 1.38 1.40 1.42 1.44 1.46 1.48 1.50 1.52 1.54 1.56 1.58 1.60 1.62 Power

Credit

This tutorial was conceived for ZiF research and tutorial workshop by Lisa DeBruine (Feb. 2020) presented again by Phillip Alday during the SMLP Summer School (Sep. 2020).

Updated and extended by Lisa Schwetlick & Daniel Backhaus, with the kind help of Phillip Alday, after changes to the package.