Skip to content
Open
Changes from 1 commit
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
160 changes: 160 additions & 0 deletions docs/notebook.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,160 @@
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# ╔═╡ 40bf1d04-76d0-11eb-0d98-cf2945aa3fd4
using DataFrames, MixedModels, MixedModelsSim, Random

# ╔═╡ 09276d18-76d2-11eb-3b91-b77f2d532e82
md"""
# Simulation of linear mixed models

This notebook shows the process of simulating responses from a mixed-effects model with crossed random effects for subject (`subj`) and item (`item`).
The general approach is to create a data frame with the appropriate experimental and grouping factors and with an _inert_ dependent variable (`dv`), meaning that, at this stage `dv` is just random noise.

This model is then used in a parametric bootstrap simulation with non-inert coefficient values (β) and random-effects variances and covariances.

Begin by loading the packages to be used.
"""

# ╔═╡ 28a121e2-76d3-11eb-0896-43da695be0a4
md"""
And initialize a random number generator.
"""

# ╔═╡ b64a889c-76d0-11eb-1936-5545db219def
rng = MersenneTwister(8675309);

# ╔═╡ ce3be2ec-76e0-11eb-3fb5-67bd31816bca
md"""
Next we set the number of subjects, `sub_n`, and the number of items, `item_n`
"""

# ╔═╡ 21958b24-76ed-11eb-11a2-41223c4864f8
begin
sub_n = 200
item_n = 50
end;

# ╔═╡ b8fe4dd2-76e2-11eb-3a7b-21e2b17043f7
md"A between-subject factor, `cond`, with two levels, `easy` and `hard` is specified as"

# ╔═╡ feb89de4-76d0-11eb-115e-89e7ed2fc483
subj_bt = Dict(:cond => ["easy", "hard"]);

# ╔═╡ eb911720-76e2-11eb-002d-21a35fd7d4b2
md"These settings are used to create a Table with an _inert_ dependent variable, `dv`, as"

# ╔═╡ 286561a4-76d1-11eb-2211-61d84bc1b9a9
frm = simdat_crossed(rng, sub_n, item_n, subj_btwn = subj_bt);

# ╔═╡ aa915fd2-76ec-11eb-2cec-bff137f10b18
md"This data table is in the form of a _row table_, which is a vector of `NamedTuple`s"

# ╔═╡ cb5f962a-76ec-11eb-3b1e-6f7f1f68bca5
typeof(frm)

# ╔═╡ 6914997e-76e3-11eb-2f2b-0501a82bee66
md"It is generally best to convert this table to a `DataFrame` for display or summary"

# ╔═╡ 9f327062-76e3-11eb-16dc-7be3a2b1c09f
df = DataFrame(frm)

# ╔═╡ b0202450-76e3-11eb-3218-8bddeed1c50d
describe(df)

# ╔═╡ 0fce98d0-76e4-11eb-0ed4-430ae7972fce
md"Note that if you go back to the entry field of `sub_n` or `item_n` and change the value, the other results are automatically updated"

# ╔═╡ 34b8677c-76e4-11eb-1807-9306323a7c9f
md"""
## Fitting a model to the inert data

In the `MixedModels` package the formula language is similar to that in the `lme4` package for `R` with the exception that it must be wrapped in a call to the macro `@formula`
"""

# ╔═╡ 48e92618-76d1-11eb-053c-e77449325b9c
m0form = @formula dv ~ 1 + cond + (1|subj) + (1|item);

# ╔═╡ 7292d5fa-76e4-11eb-243f-8f7ad404d683
md"Also, contrasts are specified in the call to fit the model instead of being part of the metadata for the factor"

# ╔═╡ 761fb192-76d1-11eb-2f8e-890af7ae08d8
m0 = fit(MixedModel, m0form, frm, contrasts=Dict(:cond => EffectsCoding()))

# ╔═╡ 91169880-76d1-11eb-28da-0f4135f0e995
md"""
As this dependent variable is inert (i.e. simulated from a standard Normal distribution without any fixed- or random-effects contributions) it is not surprising the the parameter estimates for the variance components for subject and item are close to zero as are the estimates for the fixed-effects.

## Simulated model fits for non-inert data

To simulate responses from non-intert data we must specify values for the parameters in the form of a scalar, `σ`, the standard deviation of the per-observation noise, and two vectors: `β` the fixed-effects parameter vector and `θ` the vector of relative standard deviations.

The fixed-effects parameter is relative to the model matrix, `X`, which, in this case, is the overall mean DV for a typical condition and half the difference between the `hard` and `easy` conditions. It is half the difference because we use a +/- 1 coding for the `cond` factor, producing a model matrix, `X`, of
"""

# ╔═╡ d217cc16-76e8-11eb-2f64-6ff4931c87e9
Int.(m0.X') # convert to integers to save space in the display

# ╔═╡ 9d46c892-76e9-11eb-008e-67d5100eb3fa
md"Suppose we want an overall mean of 400 ms. and the difference between `hard` and `easy` to be 50 ms., then we would set"

# ╔═╡ dbd15424-76e9-11eb-2306-436929c4e756
β = [400.0, 25.0] # Note we include the decimal point to force Float64 values

# ╔═╡ fbdfaf02-76e9-11eb-2404-d9a288cf1d31
md"We also set"

# ╔═╡ 167234ae-76ea-11eb-2589-93a3e6940453
begin
σ = 25.0 # s.d. of per-observation noise
σ₁ = 100.0 # s.d. of random effects for subject
σ₂ = 50.0 # s.d. of random effects for item
θ = [σ₁, σ₂] ./ σ
end;

# ╔═╡ b75c8580-76eb-11eb-2cca-37a0702e44a5
md"And now simulate a few thousand cases."

# ╔═╡ cd924f92-76eb-11eb-3275-eb4f55282d5d
sim = parametricbootstrap(MersenneTwister(12321), 5000, m0; β=β, σ=σ, θ=θ);
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Would it make sense to look at DataFrame(first(sim.coefpvalues)) to show what the (fixed of effects of the) mixed model looks like?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sure. I was just exploring somewhat here so my choice of what to show and what not to show was somewhat arbitrary.


# ╔═╡ 3f7a0d3e-76ec-11eb-368e-17428098bc72
ptbl = power_table(sim);

# ╔═╡ 52e23b14-76ec-11eb-1787-01b2cf9888eb
DataFrame(ptbl)

# ╔═╡ Cell order:
# ╟─09276d18-76d2-11eb-3b91-b77f2d532e82
# ╠═40bf1d04-76d0-11eb-0d98-cf2945aa3fd4
# ╟─28a121e2-76d3-11eb-0896-43da695be0a4
# ╠═b64a889c-76d0-11eb-1936-5545db219def
# ╟─ce3be2ec-76e0-11eb-3fb5-67bd31816bca
# ╠═21958b24-76ed-11eb-11a2-41223c4864f8
# ╟─b8fe4dd2-76e2-11eb-3a7b-21e2b17043f7
# ╠═feb89de4-76d0-11eb-115e-89e7ed2fc483
# ╟─eb911720-76e2-11eb-002d-21a35fd7d4b2
# ╠═286561a4-76d1-11eb-2211-61d84bc1b9a9
# ╟─aa915fd2-76ec-11eb-2cec-bff137f10b18
# ╠═cb5f962a-76ec-11eb-3b1e-6f7f1f68bca5
# ╟─6914997e-76e3-11eb-2f2b-0501a82bee66
# ╠═9f327062-76e3-11eb-16dc-7be3a2b1c09f
# ╠═b0202450-76e3-11eb-3218-8bddeed1c50d
# ╟─0fce98d0-76e4-11eb-0ed4-430ae7972fce
# ╟─34b8677c-76e4-11eb-1807-9306323a7c9f
# ╠═48e92618-76d1-11eb-053c-e77449325b9c
# ╟─7292d5fa-76e4-11eb-243f-8f7ad404d683
# ╠═761fb192-76d1-11eb-2f8e-890af7ae08d8
# ╟─91169880-76d1-11eb-28da-0f4135f0e995
# ╠═d217cc16-76e8-11eb-2f64-6ff4931c87e9
# ╟─9d46c892-76e9-11eb-008e-67d5100eb3fa
# ╠═dbd15424-76e9-11eb-2306-436929c4e756
# ╟─fbdfaf02-76e9-11eb-2404-d9a288cf1d31
# ╠═167234ae-76ea-11eb-2589-93a3e6940453
# ╟─b75c8580-76eb-11eb-2cca-37a0702e44a5
# ╠═cd924f92-76eb-11eb-3275-eb4f55282d5d
# ╠═3f7a0d3e-76ec-11eb-368e-17428098bc72
# ╠═52e23b14-76ec-11eb-1787-01b2cf9888eb