Partially-within subjects designs

Author

Phillip Alday, Douglas Bates, and Reinhold Kliegl

Published

2024-06-27

Begin by loading the packages to be used.

Code

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

CairoMakie.activate!(; type="svg")

ProgressMeter.ijulia_behavior(:clear)

Code

n_subj = 40
n_item = 3
# things are expressed as "between", so "within subjects" is "between items"
item_btwn = Dict(:frequency => ["high", "medium", "low"])
design = simdat_crossed(MersenneTwister(42), n_subj, n_item;
                        item_btwn = item_btwn)
design = DataFrame(design)

120×4 DataFrame

95 rows omitted

Row	subj	item	frequency	dv
	String	String	String	Float64
1	S01	I1	high	-0.556027
2	S02	I1	high	-0.444383
3	S03	I1	high	0.0271553
4	S04	I1	high	-0.299484
5	S05	I1	high	1.77786
6	S06	I1	high	-1.1449
7	S07	I1	high	-0.468606
8	S08	I1	high	0.156143
9	S09	I1	high	-2.64199
10	S10	I1	high	1.00331
11	S11	I1	high	1.08238
12	S12	I1	high	0.187028
13	S13	I1	high	0.518149
⋮	⋮	⋮	⋮	⋮
109	S29	I3	low	-0.338763
110	S30	I3	low	-0.0953675
111	S31	I3	low	0.768972
112	S32	I3	low	1.44244
113	S33	I3	low	-0.275032
114	S34	I3	low	-0.379637
115	S35	I3	low	-0.722696
116	S36	I3	low	0.139661
117	S37	I3	low	-1.40934
118	S38	I3	low	1.05546
119	S39	I3	low	-2.23782
120	S40	I3	low	1.15915

Code

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

3×2 DataFrame

Row	item	frequency
	String	String
1	I1	high
2	I2	medium
3	I3	low

Code

m0 = let contrasts, form
    contrasts = Dict(:frequency => HelmertCoding(base="high"))
    form = @formula(dv ~ 1 + frequency +
                    (1 + frequency | subj))
    fit(MixedModel, form, design; contrasts)
end

Minimizing 73    Time: 0:00:00 ( 4.52 ms/it)

	Est.	SE	z	p	σ_subj
(Intercept)	-0.0830	0.0983	-0.84	0.3983	0.5440
frequency: low	0.0411	0.1258	0.33	0.7438	0.7051
frequency: medium	-0.0060	0.0572	-0.11	0.9160	0.2928
Residual	0.5204

Code

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

Code

θ = createθ(m0; subj=re_subj)

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

Code

σ = 1;
β = [1.0, -3, -2];

Code

fit!(simulate!(m0; θ, β, σ))

	Est.	SE	z	p	σ_subj
(Intercept)	1.2070	0.1527	7.90	<1e-14	0.6362
frequency: low	-3.2655	0.3056	-10.68	<1e-25	1.7158
frequency: medium	-2.0013	0.2044	-9.79	<1e-21	1.1859
Residual	1.2587

Code

shrinkageplot(m0)

Code

caterpillar(m0; orderby=nothing, vline_at_zero=true)

Code

design[!, :dv] .= response(m0)

120-element Vector{Float64}:
  7.34755654258727
 10.229604583695
  3.6043175855673883
  5.3566728415298925
  4.474249009070753
  7.719899856563901
  4.025122582186671
  5.209742237453524
  5.451105316896041
  5.39384999765566
  4.332314108764023
  8.376518943225815
  4.405193973744266
  ⋮
  0.07309586813780955
  1.4581187646069615
  0.7913143551731827
  1.2178125667082424
  1.5417032495112117
 -3.1374916197365295
  1.230062037484938
  0.35610481040040387
  2.3340449614726744
  1.5733946081398242
 -0.01802135858697218
 -3.7247679264549767

Code

design_partial = filter(design) do row
    subj = parse(Int, row.subj[2:end])
    item = parse(Int, row.item[2:end])
    # for even-numbered subjects, we keep all conditions
    # for odd-numbered subjects, we keep only the two "odd" items,
    # i.e. the first and last conditions
    return iseven(subj) || isodd(item)
end
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

Code

m1 = let contrasts, form
    contrasts = Dict(:frequency => HelmertCoding(base="high"))
    form = @formula(dv ~ 1 + frequency +
                    (1 + frequency | subj))
    fit(MixedModel, form, design_partial; contrasts)
end

	Est.	SE	z	p	σ_subj
(Intercept)	1.1885	0.1832	6.49	<1e-10	0.5590
frequency: low	-3.2655	0.3056	-10.68	<1e-25	1.7496
frequency: medium	-2.0198	0.2726	-7.41	<1e-12	1.4723
Residual	1.1621

Code

shrinkageplot(m1)

Code

caterpillar(m1; orderby=nothing, vline_at_zero=true)