```
library(lme4)
library(afex)
library(emmeans)
library(dplyr)
library(JuliaCall)
```

This is the **third in a series of blog posts** working through how to generate standard errors for estimated marginal means and pooled standard errors for pairwise tests of these estimates from mixed effects models. The first shows how this can be done in R by hand for models fitted with `lme4`

, largely replicating results provided by the excellent `emmeans`

package. The second shows how this can be done in Python by hand for models fitted with `statsmodels`

. This post aims to show how this can be done in Julia, while the next post will show how this can be done in Julia. This was motivated primarily by exploring methods for fitting mixed effects models in Python and Julia and finding no suitable alternative to `emmeans`

in R for estimating marginal effects and conducting pairwise tests.

While in Julia we can use RCall to interface with R to perform model fitting with `lme4`

and follow-up tests using `emmeans`

at the time of inception there were no Julia packages that recreate the functionality of `emmeans`

for mixed effects models fitted in Julia. Since then, the Effects.jl package has been released which allows you to calculate estimated marginal means and pairwise tests by defining a reference grid. However, this doesn’t fully recreate the functionality of `emmeans`

and understanding how to do so ourselves in Julia may be informative. So, we’ll do it again by hand.

Being able to calculate estimated marginal means in Julia directly is important as the MixedModels package for fitting multilevel models is a lot faster than `lme4`

and suffers from fewer issues with convergence. This might be especially important if one were to perform simulation-based power analysis where model complexity and the number of models might be high. Here, any improvement in model fitting speed is nice to have.

This post will therefore work through fitting a mixed effects model from a factorial design in Julia using the `MixedModels.jl`

package before calculating estimated marginal means and conducting pairwise tests by hand.

# Getting Started

Here’s a quick breakdown of the packages used in each program, and why they are used:

`lme4`

: used for fitting mixed effects models in R.`afex`

: used to obtain the data set with categorical predictors.`emmeans`

: used for generating estimated marginal means and pairwise tests between group estimates in R.`dplyr`

(as always!): used for general data cleaning.

Since this blog post is created using `knitr`

as the engine from R, I need to set up JuliaCall so I can pass data from R to Julia. Let’s get this set up.

`<- julia_setup() j `

From Julia we’ll need a number of packages, with their breakdown as follows:

`RCall`

: used for interfacing between Julia and R.`DataFrames`

: unsurprisingly used to work with data frames.`CSV`

: used for loading .csv files.`Query`

: used to filter, join, sort and group data.`MixedModels`

: used for fitting linear mixed effects models.`StatsModels`

: needed to use the formula notation in defining a mixed effects model.`LinearAlgebra`

: used for matrix operations, like getting the diagnonal of a variance-covariance matrix.`Combinatorics`

: for generating combinations of pairwise contrasts.`Distributions`

: used for working with distributions such as the Normal distribution for generating confidence intervals.

`using DataFrames, CSV, MixedModels, Query, StatsModels, LinearAlgebra, Combinatorics, Distributions`

We’ll first clean the data in R to show that the data set is consistent across both R and Julia methods.

```
<- get(data("fhch2010")) |> # load data from afex
fhch filter(correct == TRUE) |> # subset the data to correct answers
select(id, task, stimulus, item, rt, log_rt) # keep only used cols
```

We’ll then make the data accessible in Julia.

`julia_assign("dat", fhch) # grab data from r`

Then we can view the data in Julia to see what it looks like.

```
dat## 12960×6 DataFrame
## Row │ id task stimulus item rt log_rt
## │ Cat… Cat… Cat… Cat… Float64 Float64
## ───────┼─────────────────────────────────────────────────────
## 1 │ N1 naming word potted 1.091 0.0870947
## 2 │ N1 naming word engine 0.876 -0.132389
## 3 │ N1 naming word ideal 0.71 -0.34249
## 4 │ N1 naming nonword uares 1.21 0.19062
## 5 │ N1 naming nonword xazz 0.843 -0.170788
## 6 │ N1 naming word fill 0.785 -0.242072
## 7 │ N1 naming nonword bounge 0.662 -0.41249
## 8 │ N1 naming nonword psems 0.713 -0.338274
## ⋮ │ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮
## 12954 │ L28 lexdec word isle 1.515 0.415415
## 12955 │ L28 lexdec word reins 1.844 0.611937
## 12956 │ L28 lexdec word idiom 1.688 0.523544
## 12957 │ L28 lexdec word lawns 1.563 0.446607
## 12958 │ L28 lexdec word spiced 1.641 0.495306
## 12959 │ L28 lexdec word weighs 1.797 0.586119
## 12960 │ L28 lexdec word darted 1.531 0.425921
## 12945 rows omitted
```

# Fitting the Models

Before we fit our model, we first have to define a dictionary of contrasts. For consistency with the R and Python walkthroughs, we’ll use dummy coding throughout. As noted in the StatsModels.jl documentation on contrasts, while “the default in R and SPSS is to use the last level as the base. Here we use the first level as the base.” For further consistency with the results from the model fitting in the previous posts, we’ll manually set the base as the last level in each factor.

Typically, for factorial designs where ANOVA-like interpretations are desired, `EffectsCoding()`

[AKA sum coding], or deviation coding should be used to ensure you obtain main effects and their interactions. This relies on using the `EffectsCoding()`

function from `StatsModels`

.

```
const contrasts = Dict(
:task => EffectsCoding(base = "naming", levels = ["naming", "lexdec"]),
:stimulus => EffectsCoding(base = "word", levels = ["word", "nonword"]),
:id => EffectsCoding(),
:item => EffectsCoding()
);
```

However, for consistency with the previous posts we’ll use dummy coding like so.

```
const contrasts = Dict(
:task => DummyCoding(base = "naming", levels = ["naming", "lexdec"]),
:stimulus => DummyCoding(base = "word", levels = ["word", "nonword"]),
:id => MixedModels.Grouping(),
:item => MixedModels.Grouping()
);
```

We can then fit our model in `MixedModels`

. This has fixed effects of task, stimulus, and their interaction along with random intercepts and slopes of stimulus by participant (`id`

) and random intercepts and slopes of task by item. For R users, the fitting process is very similar to that of `lme4`

but with the minor difference of needing to formally define the formula using the `@formula()`

notation rather than as a bare formula and having to wrap it all up in a `fit()`

function to actually fit the model you’ve defined.

```
= fit(
mod
MixedModel, @formula(
~ 1 + task * stimulus +
log_rt 1 + stimulus | id) +
(1 + task | item)
(
),
dat,= contrasts
contrasts
)## Linear mixed model fit by maximum likelihood
## log_rt ~ 1 + task + stimulus + task & stimulus + (1 + stimulus | id) + (1 + task | item)
## logLik -2 logLik AIC AICc BIC
## -3438.9542 6877.9084 6899.9084 6899.9287 6982.0742
##
## Variance components:
## Column Variance Std.Dev. Corr.
## item (Intercept) 0.0085461 0.0924449
## task: lexdec 0.0267998 0.1637064 -0.78
## id (Intercept) 0.0343237 0.1852666
## stimulus: nonword 0.0075800 0.0870631 -0.02
## Residual 0.0905882 0.3009787
## Number of obs: 12960; levels of grouping factors: 600, 45
##
## Fixed-effects parameters:
## ────────────────────────────────────────────────────────────────────────
## Coef. Std. Error z Pr(>|z|)
## ────────────────────────────────────────────────────────────────────────
## (Intercept) -0.352442 0.0421441 -8.36 <1e-16
## task: lexdec 0.343417 0.0568873 6.04 <1e-08
## stimulus: nonword 0.312304 0.0223443 13.98 <1e-43
## task: lexdec & stimulus: nonword -0.255031 0.0312596 -8.16 <1e-15
## ────────────────────────────────────────────────────────────────────────
```

How does the model look in R?

```
<- lme4::lmer(
m ~ 1 + task * stimulus +
log_rt 1 + stimulus | id) +
(1 + task | item),
(data = fhch
)summary(m)
## Linear mixed model fit by REML ['lmerMod']
## Formula: log_rt ~ 1 + task * stimulus + (1 + stimulus | id) + (1 + task |
## item)
## Data: fhch
##
## REML criterion at convergence: 6899
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.087 -0.593 -0.105 0.471 5.573
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## item (Intercept) 0.00855 0.0925
## tasklexdec 0.02683 0.1638 -0.78
## id (Intercept) 0.03592 0.1895
## stimulusnonword 0.00793 0.0891 -0.02
## Residual 0.09059 0.3010
## Number of obs: 12960, groups: item, 600; id, 45
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -0.3524 0.0431 -8.18
## tasklexdec 0.3434 0.0581 5.91
## stimulusnonword 0.3123 0.0227 13.73
## tasklexdec:stimulusnonword -0.2550 0.0318 -8.03
##
## Correlation of Fixed Effects:
## (Intr) tsklxd stmlsn
## tasklexdec -0.745
## stimlsnnwrd -0.078 0.066
## tsklxdc:stm 0.063 -0.095 -0.745
```

Those results look pretty similar with some minor differences in the estimates and standard errors likely due to different algorithms being used for the model fitting process.

# Estimating Marginal Means

First, we need to set up a grid allowing us to see the unique combinations for the levels of each factor. We’ll to this by extracting the unique values in each column of our data set and then combining them by taking their product. This is then defined as a data frame we’ve named `grid`

. Note here that we’ve had to define an empty column named `y`

which we use as a placeholder for the outcome. This is just assigned with lots of zeroes. The reason for doing so is that we’re next going to define a new model formula we’ll use to define a model matrix that we can use to estimate the marginal means and their standard errors. Every formula in Julia needs an outcome, so we must insert this `y`

column to allow for these downstream calculations.

```
begin
= unique(dat[:, :task])
task = unique(dat[:, :stimulus])
stimulus = DataFrame(Iterators.product(task, stimulus))
grid rename!(grid,[:task,:stimulus])
insertcols!(grid, 3, :y => zeros(4))
end
## 4×3 DataFrame
## Row │ task stimulus y
## │ String String Float64
## ─────┼───────────────────────────
## 1 │ naming word 0.0
## 2 │ lexdec word 0.0
## 3 │ naming nonword 0.0
## 4 │ lexdec nonword 0.0
```

Next, we’ll define a new formula. In this case since we want the estimated marginal means for the fixed effects we only need to define the fixed effects as the predictors. We can use `y`

here as our placeholder for an outcome.

```
= @formula(y ~ 1 + task * stimulus)
f ## FormulaTerm
## Response:
## y(unknown)
## Predictors:
## 1
## task(unknown)
## stimulus(unknown)
## task(unknown) & stimulus(unknown)
```

Next, we’ll define a model frame using the `ModelFrame()`

function. We pass to this function our new formula (of just fixed effects), the grid of unique values, and the contrasts we used when fitting the model. This last argument is important as the deign matrix will be different if we use a different contrast coding system.

```
= ModelFrame(f, grid, contrasts = contrasts)
mf ## ModelFrame{NamedTuple{(:y, :task, :stimulus), Tuple{Vector{Float64}, Vector{String}, Vector{String}}}, StatisticalModel}(y ~ 1 + task + stimulus + task & stimulus, StatsModels.Schema with 3 entries:
## stimulus => stimulus
## task => task
## y => y
## , (y = [0.0, 0.0, 0.0, 0.0], task = ["naming", "lexdec", "naming", "lexdec"], stimulus = ["word", "word", "nonword", "nonword"]), StatisticalModel)
```

Now that we have our model frame we can finally define the model matrix by passing this input to the `ModelMatrix()`

function.

```
= ModelMatrix(mf)
mm ## ModelMatrix{Matrix{Float64}}([1.0 0.0 0.0 0.0; 1.0 1.0 0.0 0.0; 1.0 0.0 1.0 0.0; 1.0 1.0 1.0 1.0], [1, 2, 3, 4])
```

Now we can start with the process of calculating our estimated marginal means. To do this, we need the fixed effect parameter estimates, or betas.

```
= fixef(mod)
betas ## 4-element Vector{Float64}:
## -0.352441738457094
## 0.34341675691748685
## 0.312304497311043
## -0.25503093565175144
```

Then using some matrix multiplication from our grid and fixed effect parameter estimates we can calculate our estimated marginal means. Note here that we use `mm.m`

to access the underlying matrix of the model matrix object `mm`

. Matrix multiplication in Julia is simple, requiring just the `*`

operator. While we’re at it we can drop the `y`

placeholder column for the outcome as it’s no longer needed.

```
:, :emmean] = mm.m * betas
grid[## 4-element Vector{Float64}:
## -0.352441738457094
## -0.009024981539607124
## -0.040137241146050995
## 0.048248580119684414
= select!(grid, Not([:y]))
grid ## 4×3 DataFrame
## Row │ task stimulus emmean
## │ String String Float64
## ─────┼───────────────────────────────
## 1 │ naming word -0.352442
## 2 │ lexdec word -0.00902498
## 3 │ naming nonword -0.0401372
## 4 │ lexdec nonword 0.0482486
```

Then we will generate the standard errors by pulling out the variance covariance matrix.

```
= vcov(mod)
v ## 4×4 Matrix{Float64}:
## 0.00177612 -0.00178675 -7.51705e-5 8.57959e-5
## -0.00178675 0.00323617 8.57959e-5 -0.000174455
## -7.51705e-5 8.57959e-5 0.000499266 -0.0005205
## 8.57959e-5 -0.000174455 -0.0005205 0.000977165
```

We will use this to calculate the standard errors by getting the square root of the diagnonal of the model matrix multiplied by the variance covariance matrix multiplied by the transpose of the model matrix and taking its diagonal.

```
:, :SE] = sqrt(Diagonal(mm.m * v * transpose(mm.m))).diag
grid[## 4-element Vector{Float64}:
## 0.04214405136906626
## 0.03793142983653827
## 0.04609822326134779
## 0.04145067163863225
```

Let’s see how our new grid object has changed.

```
grid## 4×4 DataFrame
## Row │ task stimulus emmean SE
## │ String String Float64 Float64
## ─────┼──────────────────────────────────────────
## 1 │ naming word -0.352442 0.0421441
## 2 │ lexdec word -0.00902498 0.0379314
## 3 │ naming nonword -0.0401372 0.0460982
## 4 │ lexdec nonword 0.0482486 0.0414507
```

How does this compare to `emmeans`

in R?

```
<- emmeans(m, ~task*stimulus, lmer.df = "a")
emm
emm## task stimulus emmean SE df asymp.LCL asymp.UCL
## naming word -0.35 0.043 Inf -0.44 -0.268
## lexdec word -0.01 0.039 Inf -0.08 0.067
## naming nonword -0.04 0.047 Inf -0.13 0.052
## lexdec nonword 0.05 0.042 Inf -0.03 0.131
##
## Degrees-of-freedom method: asymptotic
## Confidence level used: 0.95
```

That’s pretty close! Minor deviations here are likely due to the different algorithms used in the underlying model.

# Pairwise tests between marginal means

## Making pairwise mean differences

We’ll skip right to our best, automated method for computing the pairwise tests between the estimated marginal means for each combination of our levels that we used in the previous posts.

We’ll first get the names for the unique combinations for each contrast from our grid. We’ll then join the values in each column together so we have the pairs of conditions.

```
= [join(join(i, "_")) for i=zip(grid[:, :task], grid[:, :stimulus])]
combo_names ## 4-element Vector{String}:
## "naming_word"
## "lexdec_word"
## "naming_nonword"
## "lexdec_nonword"
```

Then we’ll list the get the unique combinations for each pair of conditions to get each unique pairwise contrast.

```
begin
= collect(combinations(combo_names, 2))
contrast_names = [join(join(i, " - ")) for i = contrast_names]
contrast_names_final end
## 6-element Vector{String}:
## "naming_word - lexdec_word"
## "naming_word - naming_nonword"
## "naming_word - lexdec_nonword"
## "lexdec_word - naming_nonword"
## "lexdec_word - lexdec_nonword"
## "naming_nonword - lexdec_nonword"
```

Once we have these unique pairwise contrast names we’ll get their indices so we can use these values to get the estimated marginal means of the pairwise contrasts. This just requires us to get every unique combination of the numbers 1 through to 4 (the maximum row length in our estimated marginal means given Julia is 1-indexed). We create a vector that holds vectors of two values to represent the pairwise contrasts.

```
= collect(combinations(1:length(combo_names), 2))
column_combinations ## 6-element Vector{Vector{Int64}}:
## [1, 2]
## [1, 3]
## [1, 4]
## [2, 3]
## [2, 4]
## [3, 4]
```

We can then put these into a data frame with some placeholder names for the pairs of contrasts (`i1`

and `i2`

). We now know which rows in our contrast matrix we have to subtract from one another.

```
begin
= DataFrame(i1 = [], i2 = [])
combo_cols for i in 1:length(column_combinations)
push!(combo_cols.i1, column_combinations[i][1])
push!(combo_cols.i2, column_combinations[i][2])
end
end
print(combo_cols)
## 6×2 DataFrame
## Row │ i1 i2
## │ Any Any
## ─────┼──────────
## 1 │ 1 2
## 2 │ 1 3
## 3 │ 1 4
## 4 │ 2 3
## 5 │ 2 4
## 6 │ 3 4
```

To figure out which terms need to be switched on and off in our parameter estimates we can perform vectorised subtraction on the indices. We transpose this to turn the values back into a matrix which we can use for matrix multiplication with our parameter estimates.

```
= Matrix(transpose(
model_pairwise_contrasts :, :i1], :] - mm.m[combo_cols[:, :i2], :]
mm.m[combo_cols[
))## 4×6 Matrix{Float64}:
## 0.0 0.0 0.0 0.0 0.0 0.0
## -1.0 0.0 -1.0 1.0 0.0 -1.0
## 0.0 -1.0 -1.0 -1.0 -1.0 0.0
## 0.0 0.0 -1.0 0.0 -1.0 -1.0
```

To generate the estimated marginal means for each pairwise contrast we will multiply each column of the individual pairwise contrasts from the matrix above by the betas.

```
.* betas
model_pairwise_contrasts ## 4×6 Matrix{Float64}:
## -0.0 -0.0 -0.0 -0.0 -0.0 -0.0
## -0.343417 0.0 -0.343417 0.343417 0.0 -0.343417
## 0.0 -0.312304 -0.312304 -0.312304 -0.312304 0.0
## -0.0 -0.0 0.255031 -0.0 0.255031 0.255031
```

Notice that for some contrasts this just means pulling out a single fixed effect estimate, while for others there are multiple fixed effects “switched on”. So our next step is to sum these together to get a single score for each contrast. We will put this together as a single vector.

```
= vec(sum(model_pairwise_contrasts .* betas, dims = 1))
pairwise_means ## 6-element Vector{Float64}:
## -0.34341675691748685
## -0.312304497311043
## -0.4006903185767784
## 0.03111225960644387
## -0.05727356165929154
## -0.08838582126573541
```

Now we’re ready to create a data frame of pairwise contrast names and estimated marginal means for the contrasts.

```
= DataFrame(
pairwise_tests :contrast => contrast_names_final,
:est => pairwise_means
)## 6×2 DataFrame
## Row │ contrast est
## │ String Float64
## ─────┼─────────────────────────────────────────────
## 1 │ naming_word - lexdec_word -0.343417
## 2 │ naming_word - naming_nonword -0.312304
## 3 │ naming_word - lexdec_nonword -0.40069
## 4 │ lexdec_word - naming_nonword 0.0311123
## 5 │ lexdec_word - lexdec_nonword -0.0572736
## 6 │ naming_nonword - lexdec_nonword -0.0883858
```

For getting the standard errors we’ll again use our best method for doing so in the previous posts. Below we’ll append the standard errors to our table under the heading “SE”. We do this by stepping through each contrast and calculating the standard errors for that contrast using matrix multiplication.

```
begin
= Float64[]
pairwise_SE for i in 1:nrow(pairwise_tests)
= model_pairwise_contrasts[:, i]
grad_g push!(pairwise_SE, sqrt(sum(grad_g .* v .*transpose(grad_g))))
end
= pairwise_SE
pairwise_tests.SE end;
show(pairwise_tests, allcols = true)
## 6×3 DataFrame
## Row │ contrast est SE
## │ String Float64 Float64
## ─────┼────────────────────────────────────────────────────────
## 1 │ naming_word - lexdec_word -0.343417 0.0568873
## 2 │ naming_word - naming_nonword -0.312304 0.0223443
## 3 │ naming_word - lexdec_nonword -0.40069 0.0591124
## 4 │ lexdec_word - naming_nonword 0.0311123 0.0596979
## 5 │ lexdec_word - lexdec_nonword -0.0572736 0.020867
## 6 │ naming_nonword - lexdec_nonword -0.0883858 0.0621645
```

Now let’s take a look at our updated data frame.

Next we can start calculating our test statistics and 95% confidence intervals. We make the *z*-scores by dividing each mean score by the standard error. Using these scores, we can calculate the *p*-value using the normal approximation, which for a two-tailed test is 2 \(\times\) the probability of our *z*-scores on the negative scale. We then calculate the lower and upper confidence interval limits using the normal approximation. This involves taking our parameter estimate ± the standard error multiplied by the z-value that corresponds to the upper/lower tails for the normal distribution for the 2.5th and 97.5th percentiles.

```
begin
= pairwise_tests.est./pairwise_tests.SE
pairwise_tests.z = 2*cdf.(Normal(), -abs.(pairwise_tests.z))
pairwise_tests.p = pairwise_tests.est + pairwise_tests.SE * quantile(Normal(0.0, 1.0), 1-(1+0.95)/2)
pairwise_tests.lci = pairwise_tests.est - pairwise_tests.SE * quantile(Normal(0.0, 1.0), 1-(1+0.95)/2)
pairwise_tests.uci end
## 6-element Vector{Float64}:
## -0.23191969030494863
## -0.2685105454245266
## -0.28483208704136254
## 0.1481179999457679
## -0.0163750169166273
## 0.033454298695036314
```

We can then print out our final contrasts with parameter estimates, standard errors, confidence intervals, and test statistics.

```
select!(pairwise_tests, [:contrast, :est, :SE, :lci, :uci, :z, :p]);
show(pairwise_tests, allcols = true)
## 6×7 DataFrame
## Row │ contrast est SE lci uci z p
## │ String Float64 Float64 Float64 Float64 Float64 Float64
## ─────┼─────────────────────────────────────────────────────────────────────────────────────────────────────────
## 1 │ naming_word - lexdec_word -0.343417 0.0568873 -0.454914 -0.23192 -6.03679 1.57209e-9
## 2 │ naming_word - naming_nonword -0.312304 0.0223443 -0.356098 -0.268511 -13.9769 2.1555e-44
## 3 │ naming_word - lexdec_nonword -0.40069 0.0591124 -0.516549 -0.284832 -6.77844 1.21477e-11
## 4 │ lexdec_word - naming_nonword 0.0311123 0.0596979 -0.0858935 0.148118 0.521162 0.602254
## 5 │ lexdec_word - lexdec_nonword -0.0572736 0.020867 -0.0981721 -0.016375 -2.7447 0.00605668
## 6 │ naming_nonword - lexdec_nonword -0.0883858 0.0621645 -0.210226 0.0334543 -1.42181 0.155083
```

We can even round the numeric values in the columns by defining a function to do this for us.

```
function round_df(df::DataFrame, decimals::Int)
for col in eachcol(df)
if eltype(col) <: Number
.= round.(col, digits=decimals)
col end
end
return df
end
## round_df (generic function with 1 method)
= round_df(pairwise_tests, 2);
df_rounded show(df_rounded, allcols = true)
## 6×7 DataFrame
## Row │ contrast est SE lci uci z p
## │ String Float64 Float64 Float64 Float64 Float64 Float64
## ─────┼───────────────────────────────────────────────────────────────────────────────────────
## 1 │ naming_word - lexdec_word -0.34 0.06 -0.45 -0.23 -6.04 0.0
## 2 │ naming_word - naming_nonword -0.31 0.02 -0.36 -0.27 -13.98 0.0
## 3 │ naming_word - lexdec_nonword -0.4 0.06 -0.52 -0.28 -6.78 0.0
## 4 │ lexdec_word - naming_nonword 0.03 0.06 -0.09 0.15 0.52 0.6
## 5 │ lexdec_word - lexdec_nonword -0.06 0.02 -0.1 -0.02 -2.74 0.01
## 6 │ naming_nonword - lexdec_nonword -0.09 0.06 -0.21 0.03 -1.42 0.16
```

How does this match up with `emmeans`

?

```
pairs(emm, adjust = "none")
## contrast estimate SE df z.ratio p.value
## naming word - lexdec word -0.34 0.058 Inf -5.900 <.0001
## naming word - naming nonword -0.31 0.023 Inf -13.700 <.0001
## naming word - lexdec nonword -0.40 0.060 Inf -6.600 <.0001
## lexdec word - naming nonword 0.03 0.061 Inf 0.500 0.6100
## lexdec word - lexdec nonword -0.06 0.021 Inf -2.700 0.0100
## naming nonword - lexdec nonword -0.09 0.064 Inf -1.400 0.1600
##
## Degrees-of-freedom method: asymptotic
```

That looks pretty good!

But this requires the user to make a grid of unique values in the variables they want to collapse across, respecify the formula for the fixed effects, use these to generate a model frame with the original model contrasts, and then define the model matrix. This is a lot of leg work and requires a lot of new inputs if we want to make a function.

As noted earlier, at the time when I initially wrote this code there wasn’t an alternative to doing all this nicely in a package in Julia. However, in the time between getting started on a package and this blog post, the Effects.jl package was released. This requires the user to simply define a design matrix and pass the model object to calculate the estimated marginal means. It can also allow you to calculate pairwise tests on these estimated marginal means. That said, it still lacks some of simplicity of the interface and functionality found in `emmeans`

. Some time soon I’ll release a package that tries to recreate the interface `emmeans`

(if not the full functionality) of `emmeans`

by asking users to simply pass the model object and a formula, *doing the grid construction for you*, to calculate estimated marginal means and pairwise tests.