Partially-within subjects designs

After working through this page you will be able to:

• distinguish within-unit from between-unit factors and identify a partially-within design;
• simulate data for such a design with MixedModelsSim;
• construct designs of arbitrary crossing and nesting from tidy tables with join operations;
• specify the random-effects structure appropriate to crossed subjects and items.

NoteBefore you start

Prerequisites: Analysis of the sleepstudy data.

Datasets used: simulated data (no external dataset).

1 Within-unit and between-unit factors

An experimental factor is within a grouping factor (subjects, items, schools, …) if each level of the grouping factor is observed at more than one level of the experimental factor, and between if each level of the grouping factor is observed at exactly one level. The same factor can be within one grouping factor and between another: in a typical psycholinguistic experiment, word frequency is a property of the words, so it is between items, but every subject reads both high- and low-frequency words, so it is within subjects.

This distinction matters because it determines the possible random-effects structure. A random slope for a factor by a grouping factor is only meaningful if the factor varies within that grouping factor: a subject who only ever sees one level of a factor provides no information about how their response changes across levels, so there is no subject-specific slope to estimate. In the word-frequency example, the model can include a by-subject random slope for frequency, but only a random intercept for items.

We call a design partially within when at least one experimental factor is within one grouping factor but between another, or more generally when different units are observed under different subsets of the conditions. Almost every design with crossed subjects and items is partially within in this sense; see the discussion of crossed and nested random factors in the MRK17 model-selection page for the related distinction between crossing and nesting of the grouping factors themselves.

Because no real dataset ever has exactly the design you are contemplating, this page works entirely with simulated data. Simulation is also the honest way to answer design questions (“can I even estimate this random effect?”) before collecting data — see the power simulation page for the full workflow.

Begin by loading the packages to be used.

using AlgebraOfGraphics
using CairoMakie
using DataFrames
using MixedModels
using MixedModelsSim
using MixedModelsMakie
using ProgressMeter
using Random
using StableRNGs

const progress = isinteractive()

We will use the same contrasts throughout: Helmert coding for the three-level frequency factor and effects coding for the two-level priming factor introduced later.

contrasts = Dict(:frequency => HelmertCoding(base="high"),
                 :prime => EffectsCoding(base="unrelated"))

Dict{Symbol, StatsModels.AbstractContrasts} with 2 entries:
  :prime     => EffectsCoding("unrelated", nothing)
  :frequency => HelmertCoding("high", nothing)

2 Simulating a crossed design

MixedModelsSim provides simdat_crossed, which generates a fully crossed subject-by-item design in one call. We simulate 40 subjects crossed with 30 items, with a three-level frequency factor manipulated between items.

Note the naming convention: simdat_crossed expresses everything as “between”, so a factor that is within subjects but between items is passed as item_btwn. (A factor manipulated within both subjects and items would be passed as both_win.)

n_subj = 40
n_item = 30
# a between-item factor is (here) a within-subject factor
item_btwn = Dict(:frequency => ["high", "medium", "low"])
design = DataFrame(simdat_crossed(StableRNG(42), n_subj, n_item;
                                  item_btwn=item_btwn))
first(design, 6)

6×4 DataFrame

Row	subj	item	frequency	dv
	String	String	String	Float64
1	S01	I01	high	-0.670252
2	S02	I01	high	0.447122
3	S03	I01	high	1.37363
4	S04	I01	high	1.30954
5	S05	I01	high	0.12607
6	S06	I01	high	0.683948

Every subject sees every item, and the frequency levels cycle across items, so each subject sees 10 items at each frequency level:

unique!(select(design, :item, :frequency))

30×2 DataFrame

5 rows omitted

Row	item	frequency
	String	String
1	I01	high
2	I02	medium
3	I03	low
4	I04	high
5	I05	medium
6	I06	low
7	I07	high
8	I08	medium
9	I09	low
10	I10	high
11	I11	medium
12	I12	low
13	I13	high
⋮	⋮	⋮
19	I19	high
20	I20	medium
21	I21	low
22	I22	high
23	I23	medium
24	I24	low
25	I25	high
26	I26	medium
27	I27	low
28	I28	high
29	I29	medium
30	I30	low

The dv column is filled with standard normal noise — a placeholder, not something we would ever analyze. To simulate a response with a known random-effects structure, we first construct (without fitting) a model whose formula encodes the design: frequency is within subjects, so it gets a by-subject random slope, but it is between items, so items get only a random intercept.

m0 = let form
    form = @formula(dv ~ 1 + frequency +
                    (1 + frequency | subj) + (1 | item))
    MixedModel(form, design; contrasts)
end;

Next we choose the population parameters to simulate from. create_re builds the (relative) covariance factor for one grouping factor from the standard deviations of its random effects (in the order they appear in the formula: intercept, then the two frequency contrasts) and optionally a correlation matrix:

corrmat = [ 1    0.1 -0.2
            0.1  1    0.1
           -0.2  0.1  1 ]
re_subj = create_re(1.2, 1.5, 1.5; corrmat)

3×3 LinearAlgebra.LowerTriangular{Float64, Matrix{Float64}}:
  1.2    ⋅         ⋅ 
  0.15  1.49248    ⋅ 
 -0.3   0.180907  1.45852

re_item = create_re(0.8)

1×1 LinearAlgebra.LowerTriangular{Float64, LinearAlgebra.Diagonal{Float64, Vector{Float64}}}:
 0.8

createθ assembles these into the parameter vector θ in the form the model expects:

θ = createθ(m0; subj=re_subj, item=re_item)

7-element Vector{Float64}:
  1.2
  0.15000000000000002
 -0.30000000000000004
  1.49248115565993
  0.18090680674665818
  1.4585173044131932
  0.8

Together with fixed effects β and a residual standard deviation σ, we can now overwrite the placeholder response with a simulated one and fit the model to it. Note the use of a seeded StableRNG so that the simulated data — and hence the fit — are reproducible.

σ = 1;
β = [1.0, -3, -2];
simulate!(StableRNG(1), m0; θ, β, σ);

fit!(m0; progress)

println(m0)

Linear mixed model fit by maximum likelihood
 dv ~ 1 + frequency + (1 + frequency | subj) + (1 | item)
   logLik   -2 logLik     AIC       AICc        BIC    
 -1969.7518  3939.5035  3961.5035  3961.7258  4017.4944

Variance components:
               Column      Variance Std.Dev.   Corr.
subj     (Intercept)        1.330905 1.153649
         frequency: low     2.376447 1.541573 +0.14
         frequency: medium  2.499404 1.580950 -0.26 +0.21
item     (Intercept)        0.555417 0.745263
Residual                    0.953603 0.976526
 Number of obs: 1200; levels of grouping factors: 40, 30

  Fixed-effects parameters:
─────────────────────────────────────────────────────────
                       Coef.  Std. Error      z  Pr(>|z|)
─────────────────────────────────────────────────────────
(Intercept)         0.964472    0.229306   4.21    <1e-04
frequency: low     -2.64553     0.297278  -8.90    <1e-18
frequency: medium  -1.70629     0.268588  -6.35    <1e-09
─────────────────────────────────────────────────────────

The estimates recover the simulation parameters reasonably well. The shrinkage plot shows how the by-subject predictions are pulled from the within-subject estimates (the spread-out cloud) toward the population estimate:

shrinkageplot!(Figure(; size=(500, 500)), m0, :subj)

caterpillar!(Figure(; size=(600, 650)), m0, :subj; orderby=nothing, vline_at_zero=true)

3 Building designs with tidy tables and joins

simdat_crossed is convenient, but it only generates fully crossed designs, and its keyword conventions have to be memorized. There is a more general recipe that scales to designs of arbitrary complexity, based on the tidy data principle: give each kind of experimental unit its own table, with one row per unit and one column per property of that unit. The full design then falls out of combining these small tables with join operations. This is exactly what simdat_crossed does internally — we are just doing it in the open.

3.1 One table per kind of unit

Our experiment has two kinds of units. Subjects have an identifier and (to be used shortly) an assignment to a counterbalancing list:

using SMLP2026: tagpad
subj_table = DataFrame(subj=tagpad(1:n_subj, "S"),
                       list=repeat(["A", "B"], n_subj ÷ 2))
first(subj_table, 4)

4×2 DataFrame

Row	subj	list
	String	String
1	S01	A
2	S02	B
3	S03	A
4	S04	B

Items have an identifier and a frequency class — frequency lives in the item table precisely because it is a between-item factor:

item_table = DataFrame(item=tagpad(1:n_item, "I"),
                       frequency=repeat(["high", "medium", "low"], n_item ÷ 3))
first(item_table, 4)

4×2 DataFrame

Row	item	frequency
	String	String
1	I01	high
2	I02	medium
3	I03	low
4	I04	high

(The helper tagpad produces zero-padded labels such as "S01" so that string sorting matches numeric order.)

3.2 Full crossing with crossjoin

A fully crossed design is the Cartesian product of the unit tables, which is what crossjoin computes: every row of the first table paired with every row of the second.

crossed = crossjoin(select(subj_table, :subj), item_table)
nrow(crossed)

1200

This is the same design that simdat_crossed produced, which we can verify with a join: antijoin returns the rows of its first argument that have no match in the second, so an empty anti-join means every design cell we built is present in design.

isempty(antijoin(crossed, design; on=[:subj, :item, :frequency]))

true

3.3 Attaching covariates with leftjoin

Because unit-level properties live in unit-level tables, adding a covariate to the design is a one-line join. leftjoin keeps all rows of its first argument and copies in the matching columns of the second; each item’s word length is automatically repeated for every subject who sees that item.

item_covariates = DataFrame(item=item_table.item,
                            wordlength=rand(StableRNG(2026), 3:9, n_item))
leftjoin!(crossed, item_covariates; on=:item)
first(crossed, 4)

4×4 DataFrame

Row	subj	item	frequency	wordlength
	String	String	String	Int64?
1	S01	I01	high	3
2	S01	I02	medium	3
3	S01	I03	low	8
4	S01	I04	high	6

The same pattern adds subject covariates (join on :subj), and it works just as well on real data as on simulated designs — see the power simulation page for it in action.

3.4 Counterbalancing with a Latin square

Now for a factor that simdat_crossed cannot express: a priming manipulation (related vs. unrelated prime) that is within both subjects and items, but where each subject sees each item only once — so each subject sees each item in only one of the two conditions. The standard solution is a Latin square: items are rotated through the conditions across counterbalancing lists, and each subject is assigned to one list.

The rotation is itself just a small tidy table. In wide form, one row per item and one column per list:

square = DataFrame(item=item_table.item,
                   A=repeat(["related", "unrelated"], n_item ÷ 2),
                   B=repeat(["unrelated", "related"], n_item ÷ 2))
first(square, 4)

4×3 DataFrame

Row	item	A	B
	String	String	String
1	I01	related	unrelated
2	I02	unrelated	related
3	I03	related	unrelated
4	I04	unrelated	related

stack reshapes it to long (“tidy”) form — one row per item-list combination:

rotation = stack(square, [:A, :B]; variable_name=:list, value_name=:prime)
first(rotation, 4)

4×3 DataFrame

Row	item	list	prime
	String	String	String
1	I01	A	related
2	I02	A	unrelated
3	I03	A	related
4	I04	A	unrelated

The full design is then: cross subjects with items, and use the (list, item) pair to look up each trial’s condition in the rotation table.

design_ls = leftjoin(crossjoin(subj_table, item_table), rotation;
                     on=[:list, :item])
disallowmissing!(design_ls)
nrow(design_ls)

1200

One technical note: leftjoin cannot know in advance that every row will find a match, so the joined columns allow missing values. disallowmissing! asserts that no missing values are actually present and narrows the column types accordingly, which keeps the modeling code downstream simple.

The counterbalancing worked: each subject sees every combination of frequency and priming, in equal numbers, even though no subject sees any item twice.

combine(groupby(design_ls, [:list, :frequency, :prime]), nrow => :n)

12×4 DataFrame

Row	list	frequency	prime	n
	String	String	String	Int64
1	A	high	related	100
2	A	medium	unrelated	100
3	A	low	related	100
4	A	high	unrelated	100
5	A	medium	related	100
6	A	low	unrelated	100
7	B	high	unrelated	100
8	B	medium	related	100
9	B	low	unrelated	100
10	B	high	related	100
11	B	medium	unrelated	100
12	B	low	related	100

The general recipe, then, is:

• crossjoin crosses: units are combined factorially.
• leftjoin on a shared key nests and annotates: rows inherit the properties of the unit (or unit combination) they belong to.
• antijoin (or subset) removes cells: deliberately, for incomplete designs such as Latin squares beyond two conditions, or to mimic anticipated data loss.

Any crossing or nesting structure — students within classes within schools, items within sublists within lists, trials within blocks — can be built by composing these operations on small per-unit tables.

4 Fitting the partially-within design

The design now licenses a richer random-effects structure. Priming is within subjects and within items, so it can (and should) get random slopes for both. Frequency remains between items:

design_ls.dv = zeros(nrow(design_ls));
m1 = let form
    form = @formula(dv ~ 1 + frequency * prime +
                    (1 + frequency + prime | subj) +
                    (1 + prime | item))
    MixedModel(form, design_ls; contrasts)
end;

As before, we choose population parameters and simulate a response. The by-subject structure now has four dimensions (intercept, two frequency contrasts, one priming contrast) and the by-item structure two (intercept, priming contrast); for simplicity we simulate with uncorrelated random effects.

re_subj1 = create_re(1.2, 1.5, 1.5, 1.0)
re_item1 = create_re(0.8, 0.5)
θ1 = createθ(m1; subj=re_subj1, item=re_item1)
β1 = [1.0, -3, -2, 1.5, 0.5, -0.5];
simulate!(StableRNG(1), m1; θ=θ1, β=β1, σ);

fit!(m1; progress)

println(m1)

Linear mixed model fit by maximum likelihood
 dv ~ 1 + frequency + prime + frequency & prime + (1 + frequency + prime | subj) + (1 + prime | item)
   logLik   -2 logLik     AIC       AICc        BIC    
 -2083.1022  4166.2043  4206.2043  4206.9168  4308.0058

Variance components:
               Column      Variance Std.Dev.   Corr.
subj     (Intercept)        1.091728 1.044858
         frequency: low     1.849664 1.360024 +0.37
         frequency: medium  2.950267 1.717634 -0.09 +0.34
         prime: related     0.967740 0.983738 +0.23 +0.02 -0.30
item     (Intercept)        0.965141 0.982416
         prime: related     0.250507 0.500507 -0.32
Residual                    0.990055 0.995015
 Number of obs: 1200; levels of grouping factors: 40, 30

  Fixed-effects parameters:
──────────────────────────────────────────────────────────────────────────
                                        Coef.  Std. Error      z  Pr(>|z|)
──────────────────────────────────────────────────────────────────────────
(Intercept)                          1.09416    0.245539    4.46    <1e-05
frequency: low                      -3.08266    0.309413   -9.96    <1e-22
frequency: medium                   -2.10415    0.300424   -7.00    <1e-11
prime: related                       1.61106    0.182671    8.82    <1e-17
frequency: low & prime: related      0.436779   0.117316    3.72    0.0002
frequency: medium & prime: related  -0.434119   0.0677321  -6.41    <1e-09
──────────────────────────────────────────────────────────────────────────

We keep the simulated response for reuse below:

design_ls.dv .= response(m1);

The shrinkage plots now come in two flavors, one per grouping factor. Note that the by-item panels involve only the intercept and the priming slope — there is no by-item frequency slope to plot, because the design does not support one.

shrinkageplot!(Figure(; size=(700, 700)), m1, :subj)

shrinkageplot!(Figure(; size=(400, 400)), m1, :item)

caterpillar!(Figure(; size=(450, 550)), m1, :item; orderby=nothing, vline_at_zero=true)

5 When the design degrades: missing cells

Real datasets rarely arrive as balanced as they were designed. Suppose the medium-frequency items turn out to be unusable for the odd-numbered subjects (an equipment problem, a bad session — the reason does not matter). Deliberately removing design cells is again a join: build a table of the lost trials and antijoin it away.

lost = crossjoin(subj_table[1:2:end, [:subj]],
                 subset(item_table, :frequency => ByRow(==("medium")))[:, [:item]])
design_partial = antijoin(design_ls, lost; on=[:subj, :item])
nrow(design_partial)

1000

For half of the subjects, frequency is now effectively a two-level factor:

sort!(unique!(select(design_partial, :subj, :frequency)), :subj)

100×2 DataFrame

75 rows omitted

Row	subj	frequency
	String	String
1	S01	high
2	S01	low
3	S02	high
4	S02	medium
5	S02	low
6	S03	high
7	S03	low
8	S04	high
9	S04	medium
10	S04	low
11	S05	high
12	S05	low
13	S06	high
⋮	⋮	⋮
89	S36	medium
90	S36	low
91	S37	high
92	S37	low
93	S38	high
94	S38	medium
95	S38	low
96	S39	high
97	S39	low
98	S40	high
99	S40	medium
100	S40	low

Crucially, the same model can still be fit. The mixed model does not require every subject to appear in every cell; subjects simply contribute information only about the contrasts their data speak to, and their predicted effects for the missing cells are informed by the population estimates (partial pooling).

m2 = let form
    form = @formula(dv ~ 1 + frequency * prime +
                    (1 + frequency + prime | subj) +
                    (1 + prime | item))
    fit(MixedModel, form, design_partial; contrasts, progress)
end

println(m2)

Linear mixed model fit by maximum likelihood
 dv ~ 1 + frequency + prime + frequency & prime + (1 + frequency + prime | subj) + (1 + prime | item)
   logLik   -2 logLik     AIC       AICc        BIC    
 -1746.1921  3492.3842  3532.3842  3533.2422  3630.5393

Variance components:
               Column      Variance Std.Dev.   Corr.
subj     (Intercept)        1.086595 1.042399
         frequency: low     1.863285 1.365022 +0.19
         frequency: medium  2.224399 1.491442 -0.33 +0.26
         prime: related     0.938871 0.968954 +0.16 +0.06 -0.38
item     (Intercept)        0.968029 0.983884
         prime: related     0.225554 0.474926 -0.55
Residual                    1.008795 1.004388
 Number of obs: 1000; levels of grouping factors: 40, 30

  Fixed-effects parameters:
──────────────────────────────────────────────────────────────────────────
                                        Coef.  Std. Error      z  Pr(>|z|)
──────────────────────────────────────────────────────────────────────────
(Intercept)                          1.119       0.272349   4.11    <1e-04
frequency: low                      -3.08266     0.310233  -9.94    <1e-22
frequency: medium                   -2.0793      0.293368  -7.09    <1e-11
prime: related                       1.59614     0.21069    7.58    <1e-13
frequency: low & prime: related      0.436779    0.111976   3.90    <1e-04
frequency: medium & prime: related  -0.449038    0.12938   -3.47    0.0005
──────────────────────────────────────────────────────────────────────────

Compare the by-subject shrinkage with the complete design above: the predictions for the subjects with missing cells are pulled more strongly toward the population values, because less data means the population estimate carries relatively more weight.

shrinkageplot!(Figure(; size=(700, 700)), m2, :subj)

caterpillar(m2, :subj; orderby=nothing, vline_at_zero=true)

This robustness has limits: if a factor becomes between subjects for all subjects, the by-subject slope is no longer identified at all, and heavily unbalanced designs yield noisier variance-component estimates. Simulating the anticipated imbalance before running the study — exactly as done here — is the way to find out whether your design still supports the model you intend to fit.

6 See Also

• The multiple random-effects terms chapter of Embrace Uncertainty — crossed, nested, and partially crossed random-effects structures in depth.
• The MRK17 model-specification page — between/within factors and crossed vs. nested grouping factors in a real dataset.
• The power simulation page — using simulated designs (including join-based covariates) to plan sample sizes.

7 Exercises

1. Classify the factors. For the Latin-square design simulated here, which experimental factors are within subjects, which are within items, and which are within both? Why does that distinction determine where each factor can appear in the random-effects structure?

NoteSolution

Frequency is between items (each item has one frequency class) but within subjects (each subject sees all three classes); priming is within both subjects and items (each subject sees both conditions, and each item appears in both conditions — across different subjects). A factor can have a by-subject random slope only if it varies within subjects, and likewise for items: so frequency gets a by-subject slope but no by-item slope, while priming gets both. A factor that is between a grouping factor cannot have a random slope for it — there is no within-unit variation for the slope to capture.

2. Change the design size. Re-run the simulation with fewer items (say n_item = 6). How does decreasing n_item affect your ability to estimate the by-item random effects?

NoteSolution

With very few items the by-item variance components and any item-level correlations are poorly estimated and the fit is prone to singularity; increasing n_item gives the model more information per item-level effect and generally stabilizes those estimates. Try n_item = 6, 30, 90 and compare issingular and the estimated item variances.

3. Build a nested design. Using only small tidy tables and joins, construct the design table for a study with 4 schools, 3 classes within each school, and 5 students within each class, where each student takes the same 2 tests. (Hint: nesting means that a class “belongs to” its school, so the class table should carry the school identifier.)

NoteSolution

Nesting is expressed by building each lower-level table from the table above it, so that every row carries its parent’s key; crossing is expressed with crossjoin.

schools = DataFrame(school=tagpad(1:4, "H"))
classes = transform(crossjoin(schools, DataFrame(classnum=1:3)),
                    [:school, :classnum] => ByRow((s, c) -> string(s, "c", c)) => :class)
students = transform(crossjoin(classes, DataFrame(studnum=1:5)),
                     [:class, :studnum] => ByRow((c, s) -> string(c, "s", s)) => :student)
design = crossjoin(select(students, :school, :class, :student),
                   DataFrame(test=["math", "reading"]))
nrow(design)  # 4 * 3 * 5 * 2 = 120

Because each class label contains its school (and each student label its class), the identifiers are globally unique and the grouping factors can be used directly in a model formula such as (1 | school) + (1 | class) + (1 | student). See the discussion of nested factors on the MRK17 page for why unique labels matter.

---

This page was rendered from git revision 29e8d33 using Quarto 1.10.18 and Julia 1.12.7.