Code
using AlgebraOfGraphics
using CairoMakie
using Chain
using DataFrames
using MixedModels
using SMLP2024: dataset
using StatsBase
using StatsModels
activate!(; type="png")
CairoMakie.
= false progress
false
Phillip Alday, Douglas Bates, and Reinhold Kliegl
2024-09-13
false
Many researchers have pointed out that contrasts should be “tested instead of, rather than as a supplement to, the ordinary ‘omnibus’ F test” (Hays, 1973, p. 601).
For a (quasi-)experimental set of data, there is (or should be) a clear a priori theoretical commitment to specific hypotheses about differences between factor levels and how these differences enter in interactions with other factors. This specification should be used in the first LMM and reported, irrespective of the outcome. If alternative theories lead to alternative a priori contrast specifications, both analyses are justified. If the observed means render the specification completely irrelevant, the comparisons originally planned could still be reported in a Supplement).
In this script, we are working through a large number of different contrasts for the same data. The purpose is to introduce both the preprogrammed (“canned”) and the general options to specify hypotheses about main effects and interactions. Obviously, we do not endorse generating a plot of the means and specifying the contrasts accordingly. This is known as the Texas sharpshooter fallacy. The link leads to an illustration and brief historical account by Wagenmakers (2018).
Irrespective of how results turn out, there is nothing wrong with specifying a set of post-hoc contrasts to gain a better understanding of what the data are trying to tell us. Of course, in an article or report about the study, the a priori and post-hoc nature of contrast specifications must be made clear. Some kind of alpha-level adjustment (e.g., Bonferroni) may be called for, too. And, of course, there are grey zones.
There is quite a bit of statistical literature on contrasts. Two “local” references are Brehm & Alday (2022) and Schad et al. (2020).
For further readings see “Further Readings” in Schad et al. (2020).
We take the KWDYZ
dataset from Kliegl et al. (2011). This is an experiment looking at three effects of visual cueing under four different cue-target relations (CTRs). Two horizontal rectangles are displayed above and below a central fixation point or they displayed in vertical orientation to the left and right of the fixation point. Subjects react to the onset of a small visual target occurring at one of the four ends of the two rectangles. The target is cued validly on 70% of trials by a brief flash of the corner of the rectangle at which it appears; it is cued invalidly at the three other locations 10% of the trials each.
We specify three contrasts for the four-level factor CTR that are derived from spatial, object-based, and attractor-like features of attention. They map onto sequential differences between appropriately ordered factor levels.
We also have a dataset from a replication and extension of this study Kliegl et al. (2015) Both data sets are available in R-package RePsychLing
dat1 = DataFrame(dataset(:kwdyz11))
cellmeans = combine(
groupby(dat1, [:CTR]),
:rt => mean,
:rt => std,
:rt => length,
:rt => (x -> std(x) / sqrt(length(x))) => :rt_semean,
)
Row | CTR | rt_mean | rt_std | rt_length | rt_semean |
---|---|---|---|---|---|
String | Float32 | Float32 | Int64 | Float64 | |
1 | val | 358.032 | 83.4579 | 20141 | 0.588067 |
2 | sod | 391.267 | 92.662 | 2863 | 1.73177 |
3 | dos | 405.146 | 92.6892 | 2843 | 1.73836 |
4 | dod | 402.3 | 95.3913 | 2863 | 1.78278 |
We use the same formula for all analyses
This is the default order of factor levels.
Controlling the ordering of levels for contrasts:
levels
to order the levelsbase
a different level can be specified.The SeqDiffCoding
contrast corresponds to MASS::contr.sdif()
in R. The assignment of random factors such as Subj
to Grouping()
is necessary when the sample size is very large. We recommend to include it always, but in this tutorial we do so only in the first example.
m1 = let levels = ["val", "sod", "dos", "dod"]
contrasts = Dict(
:CTR => SeqDiffCoding(; levels),
)
fit(MixedModel, form, dat1; contrasts, progress)
end
Est. | SE | z | p | σ_Subj | |
(Intercept) | 389.7336 | 7.0913 | 54.96 | <1e-99 | 55.1998 |
CTR: sod | 33.7817 | 3.2874 | 10.28 | <1e-24 | 23.2486 |
CTR: dos | 13.9852 | 2.3057 | 6.07 | <1e-08 | 10.7538 |
CTR: dod | -2.7470 | 2.2143 | -1.24 | 0.2148 | 9.5110 |
Residual | 69.8348 |
What does the intercept represent?
Grand Mean is mean of condition means.
HypothesisCoding
is the most general option available. We can implement all “canned” contrasts ourselves. The next example reproduces the test statistics from SeqDiffCoding
- with a minor modification illustrating the flexibility of going beyond the default version.
m1b = let levels = ["val", "sod", "dos", "dod"]
contrasts = Dict(
:CTR => HypothesisCoding(
[
-1 1 0 0
0 -1 1 0
0 0 1 -1
];
levels,
labels=["spt", "obj", "grv"],
),
)
fit(MixedModel, form, dat1; contrasts, progress)
end
Est. | SE | z | p | σ_Subj | |
(Intercept) | 389.7336 | 7.0911 | 54.96 | <1e-99 | 55.1984 |
CTR: spt | 33.7817 | 3.2876 | 10.28 | <1e-24 | 23.2501 |
CTR: obj | 13.9852 | 2.3060 | 6.06 | <1e-08 | 10.7572 |
CTR: grv | 2.7469 | 2.2146 | 1.24 | 0.2148 | 9.5151 |
Residual | 69.8348 |
The difference to the preprogrammed SeqDiffCoding
is that for the third contrast we changed the direction of the contrast such that the sign of the effect is positive when the result is in agreement with theoretical expectation, that is we subtract the fourth level from the third, not the third level from the fourth.
This contrast corresponds to contr.treatment()
in R
m2 = let
contrasts = Dict(:CTR => DummyCoding(; base="val"))
fit(MixedModel, form, dat1; contrasts, progress)
end
Est. | SE | z | p | σ_Subj | |
(Intercept) | 358.0914 | 6.1534 | 58.19 | <1e-99 | 47.9055 |
CTR: dod | 45.0200 | 4.3625 | 10.32 | <1e-24 | 32.2821 |
CTR: dos | 47.7669 | 3.5560 | 13.43 | <1e-40 | 25.5314 |
CTR: sod | 33.7817 | 3.2868 | 10.28 | <1e-24 | 23.2436 |
Residual | 69.8349 |
The DummyCoding
contrast has the disadvantage that the intercept returns the mean of the level specified as base
, default is the first level, not the GM.
The contrasts returned by DummyCoding
may be exactly what we want. Can’t we have them, but also have the intercept estimate the GM, rather than the mean of the base level? Yes, we can! We call this “You can have your cake and it eat, too”-Coding (YchycaeitCoding). And we use HypothesisCoding
to achieve this outcome.
m2b = let levels = ["val", "sod", "dos", "dod"]
contrasts = Dict(
:CTR => HypothesisCoding(
[
-1 1 0 0
-1 0 1 0
-1 0 0 1
];
levels,
labels=levels[2:end],
)
)
fit(MixedModel, form, dat1; contrasts, progress)
end
Est. | SE | z | p | σ_Subj | |
(Intercept) | 389.7336 | 7.0911 | 54.96 | <1e-99 | 55.1978 |
CTR: sod | 33.7817 | 3.2873 | 10.28 | <1e-24 | 23.2477 |
CTR: dos | 47.7669 | 3.5566 | 13.43 | <1e-40 | 25.5370 |
CTR: dod | 45.0200 | 4.3634 | 10.32 | <1e-24 | 32.2893 |
Residual | 69.8348 |
We can simply move the column with -1s for a different base.
m2c = let levels = ["val", "sod", "dos", "dod"]
contrasts = Dict(
:CTR => HypothesisCoding(
[
1 -1 0 0
0 -1 1 0
0 -1 0 1
];
levels,
labels=["val", "dos", "dod"],
)
)
fit(MixedModel, form, dat1; contrasts, progress)
end
Est. | SE | z | p | σ_Subj | |
(Intercept) | 389.7336 | 7.0907 | 54.96 | <1e-99 | 55.1949 |
CTR: val | -33.7818 | 3.2884 | -10.27 | <1e-24 | 23.2574 |
CTR: dos | 13.9852 | 2.3041 | 6.07 | <1e-08 | 10.7332 |
CTR: dod | 11.2382 | 2.3562 | 4.77 | <1e-05 | 11.4364 |
Residual | 69.8349 |
We can simply relevel the factor with a different base.
m2d = let levels = ["sod", "val", "dos", "dod"]
contrasts = Dict(
:CTR => HypothesisCoding(
[
-1 1 0 0
-1 0 1 0
-1 0 0 1
];
levels,
labels=levels[2:end],
)
)
fit(MixedModel, form, dat1; contrasts)
end
Est. | SE | z | p | σ_Subj | |
(Intercept) | 389.7336 | 7.0915 | 54.96 | <1e-99 | 55.2013 |
CTR: val | -33.7817 | 3.2874 | -10.28 | <1e-24 | 23.2483 |
CTR: dos | 13.9852 | 2.3057 | 6.07 | <1e-08 | 10.7534 |
CTR: dod | 11.2383 | 2.3585 | 4.77 | <1e-05 | 11.4644 |
Residual | 69.8348 |
EffectsCoding
estimates the difference between the Grand Mean and three of the four levels. The difference of the fourth levels can be computed from the Grand Mean and these three differences.
m3 = let levels = ["val", "sod", "dos", "dod"]
contrasts = Dict(:CTR => EffectsCoding(; levels, base="val"))
fit(MixedModel, form, dat1; contrasts, progress)
end
Est. | SE | z | p | σ_Subj | |
(Intercept) | 389.7336 | 7.0780 | 55.06 | <1e-99 | 55.0953 |
CTR: sod | 2.1396 | 1.3336 | 1.60 | 0.1086 | 6.0051 |
CTR: dos | 16.1248 | 1.4401 | 11.20 | <1e-28 | 7.3284 |
CTR: dod | 13.3778 | 1.9331 | 6.92 | <1e-11 | 12.4701 |
Residual | 69.8351 |
This contrast corresponds almost to contr.sum()
in R. The “almost” qualification refers to the fact that EffectsCoding
uses the first level as default base; contr.sum()
uses the last factor level.
m3b = let levels = ["val", "sod", "dos", "dod"]
contrasts = Dict(:CTR => EffectsCoding(; levels, base = "dod"))
fit(MixedModel, form, dat1; contrasts, progress)
end
Est. | SE | z | p | σ_Subj | |
(Intercept) | 389.7336 | 7.0913 | 54.96 | <1e-99 | 55.1994 |
CTR: val | -31.6422 | 2.6425 | -11.97 | <1e-32 | 19.9522 |
CTR: sod | 2.1396 | 1.3337 | 1.60 | 0.1087 | 6.0063 |
CTR: dos | 16.1248 | 1.4404 | 11.19 | <1e-28 | 7.3312 |
Residual | 69.8348 |
How could we achieve the default result with HypothesisCoding?
m3c = let levels = ["val", "sod", "dos", "dod"]
contrasts = Dict(
:CTR => HypothesisCoding(
[
-1/4 3/4 -1/4 -1/4 # b - GM = b - (a+b+c+d)/4 => -1/4*a + 3/4*b - 1/4*c - 1/4*d
-1/4 -1/4 3/4 -1/4 # c - GM = c - (a+b+c+d)/4 => -1/4*a - 1/4*b + 3/4*c - 1/4*d
-1/4 -1/4 -1/4 3/4 # d - GM = d - (a+b+c+d)/4 => -1/4*a - 1/4*b - 1/4*c + 3/4*d
];
levels,
labels=levels[2:end],
)
)
fit(MixedModel, form, dat1; contrasts, progress)
end
Est. | SE | z | p | σ_Subj | |
(Intercept) | 389.7336 | 7.0920 | 54.95 | <1e-99 | 55.2050 |
CTR: sod | 2.1395 | 1.3337 | 1.60 | 0.1087 | 6.0067 |
CTR: dos | 16.1248 | 1.4405 | 11.19 | <1e-28 | 7.3330 |
CTR: dod | 13.3778 | 1.9324 | 6.92 | <1e-11 | 12.4640 |
Residual | 69.8348 |
HelmertCoding
codes each level as the difference from the average of the lower levels. With the default order of CTR
levels we get the following test statistics. These contrasts are orthogonal.
m4 = let levels = ["val", "sod", "dos", "dod"]
contrasts = Dict(:CTR => HelmertCoding(; levels))
fit(MixedModel, form, dat1; contrasts, progress)
end
Est. | SE | z | p | σ_Subj | |
(Intercept) | 389.7336 | 7.0912 | 54.96 | <1e-99 | 55.1986 |
CTR: sod | 16.8909 | 1.6436 | 10.28 | <1e-24 | 11.6232 |
CTR: dos | 10.2920 | 0.8353 | 12.32 | <1e-34 | 5.2568 |
CTR: dod | 4.4593 | 0.6441 | 6.92 | <1e-11 | 4.1545 |
Residual | 69.8349 |
+ HeC1: (b - a)/2 # (391 - 358)/2 = 16.5
+ HeC2: (c - (b+a)/2)/3 # (405 - (391 + 358)/2)/3 = 10.17
+ HeC3: (d - (c+b+a)/3)/4 # (402 - (405 + 391 + 358)/3)/4 = 4.33
We can reconstruct the estimates, but they are scaled by the number of levels involved. With HypothesisCoding
we can estimate the “unscaled” differences. Also the labeling of the contrasts is not as informative as they could be.
m4b = let levels = ["val", "sod", "dos", "dod"]
contrasts = Dict(
:CTR => HypothesisCoding(
[
-1 1 0 0
-1/2 -1/2 1 0
-1/3 -1/3 -1/3 1
];
levels,
labels= ["2-1", "3-21", "4-321"]
)
)
fit(MixedModel, form, dat1; contrasts, progress)
end
Est. | SE | z | p | σ_Subj | |
(Intercept) | 389.7336 | 7.0886 | 54.98 | <1e-99 | 55.1784 |
CTR: 2-1 | 33.7817 | 3.2874 | 10.28 | <1e-24 | 23.2487 |
CTR: 3-21 | 30.8761 | 2.5065 | 12.32 | <1e-34 | 15.7751 |
CTR: 4-321 | 17.8371 | 2.5765 | 6.92 | <1e-11 | 16.6185 |
Residual | 69.8349 |
Reverse HelmertCoding
codes each level as the difference from the average of the higher levels. To estimate these effects we simply reverse the order of factor levels. Of course, the contrasts are also orthogonal.
m4c = let levels = reverse(StatsModels.levels(dat1.CTR))
contrasts = Dict(:CTR => HelmertCoding(; levels))
fit(MixedModel, form, dat1; contrasts, progress)
end
Est. | SE | z | p | σ_Subj | |
(Intercept) | 389.7336 | 7.0930 | 54.95 | <1e-99 | 55.2132 |
CTR: dos | 1.3734 | 1.1072 | 1.24 | 0.2148 | 4.7557 |
CTR: sod | -4.2039 | 0.6845 | -6.14 | <1e-09 | 3.3523 |
CTR: val | -10.5474 | 0.8808 | -11.97 | <1e-32 | 6.6508 |
Residual | 69.8347 |
+ HeC1:(c - d)/2 # (405 - 402)/2 = 1.5
+ HeC2:(b - (c+d)/2)/3 # (391 - (405 + 402)/2)/3 = -4.17
+ HeC3:(a - (b+c+d)/3/4 # (356 -(391 + 405 + 402)/3)/4 = -10.83
… and the unscaled-by-number-of-levels estimates.
m4d = let levels = ["val", "sod", "dos", "dod"]
contrasts = Dict(
:CTR => HypothesisCoding(
[
0 0 1 -1
0 1 -1/2 -1/2
1 -1/3 -1/3 -1/3
];
levels,
labels= ["3-4", "2-34", "1-234"]
)
)
fit(MixedModel, form, dat1; contrasts)
end
Est. | SE | z | p | σ_Subj | |
(Intercept) | 389.7336 | 7.0915 | 54.96 | <1e-99 | 55.2014 |
CTR: 3-4 | 2.7469 | 2.2143 | 1.24 | 0.2148 | 9.5107 |
CTR: 2-34 | -12.6117 | 2.0529 | -6.14 | <1e-09 | 10.0486 |
CTR: 1-234 | -42.1895 | 3.5231 | -11.98 | <1e-32 | 26.6007 |
Residual | 69.8348 |
For factors with more than four levels there are many options for specifying orthogonal contrasts as long as one proceeds in a top-down strictly hierarchical fashion.
Suppose you have a factor with seven levels and let’s ignore shifting columns. In this case, you have six options for the first contrast, that is 6 vs. 1, 5 vs.2 , 4 vs. 3, 3 vs. 4, 2 vs. 5, and 1 vs. 6 levels. Then, you specify orthogonal contrasts for partitions with more than 2 elements and so on. That is, you don’t specify a contrast that crosses an earlier partition line.
In the following example, after an initial 4 vs 3 partitioning of levels, we specify AnovaCoding
for the left and HelmertCoding
for the right partition.
contrasts = Dict(
:CTR => HypothesisCoding(
[
-1/4 -1/4 -1/4 -1/4 +1/3 +1/3 +1/3
-1/2 -1/2 +1/2 +1/2 0 0 0
-1/2 +1/2 -1/2 +1/2 0 0 0
+1/2 -1/2 -1/2 +1/2 0 0 0
0 0 0 0 -1 +1 0
0 0 0 0 -1/2 -1/2 1
];
levels=["A1", "A2", "A3", "A4", "A5", "A6", "A7"],
labels=["c567.1234", "B", "C", "BxC", "c6.5", "c6.56"],
),
);
There are two rules that hold for all orthogonal contrasts:
Factorial designs (i.e., lab experiments) are traditionally analyzed with analysis of variance. The test statistics of main effects and interactions are based on an orthogonal set of contrasts. We specify them with HypothesisCoding
.
An A(2) x B(2) design can be recast as an F(4) design with the levels (A1-B1, A1-B2, A2-B1, A2-B2). The following contrast specification returns estimates for the main effect of A, the main effect of B, and the interaction of A and B. In a figure With A on the x-axis and the levels of B shown as two lines, the interaction tests the null hypothesis that the two lines are parallel. A positive coefficient implies overadditivity (diverging lines toward the right) and a negative coefficient underadditivity (converging lines).
m5 = let levels = ["val", "sod", "dos", "dod"]
contrasts = Dict(
:CTR => HypothesisCoding(
[
-1 -1 +1 +1 # A
-1 +1 -1 +1 # B
+1 -1 -1 +1 # A x B
];
levels,
labels=["A", "B", "AxB"],
),
)
fit(MixedModel, form, dat1; contrasts, progress)
end
Est. | SE | z | p | σ_Subj | |
(Intercept) | 389.7336 | 7.0913 | 54.96 | <1e-99 | 55.1994 |
CTR: A | 59.0052 | 5.1822 | 11.39 | <1e-29 | 36.2045 |
CTR: B | 31.0348 | 4.6748 | 6.64 | <1e-10 | 31.7121 |
CTR: AxB | -36.5287 | 3.0926 | -11.81 | <1e-31 | 16.0038 |
Residual | 69.8348 |
It is also helpful to see the corresponding layout of the four means for the interaction of A and B (i.e., the third contrast)
B1 B2
A1 +1 -1
A2 -1 +1
Thus, interaction tests whether the difference between main diagonal and minor diagonal is different from zero.
Going beyond the four level factor; it is also helpful to see the corresponding layout of the eight means for the interaction of A and B and C.
C1 C2
B1 B2 B1 B2
A1 +1 -1 A1 -1 +1
A2 -1 +1 A2 +1 -1
Nested contrasts are often specified as follow up as post-hoc tests for ANOVA interactions. They are orthogonal. We specify them with HypothesisCoding
.
An A(2) x B(2) design can be recast as an F(4) design with the levels (A1-B1, A1-B2, A2-B1, A2-B2). The following contrast specification returns an estimate for the main effect of A and the effects of B nested in the two levels of A. In a figure With A on the x-axis and the levels of B shown as two lines, the second contrast tests whether A1-B1 is different from A1-B2 and the third contrast tests whether A2-B1 is different from A2-B2.
m6 = let levels = ["val", "sod", "dos", "dod"]
contrasts = Dict(
:CTR => HypothesisCoding(
[
-1 -1 +1 +1
-1 +1 0 0
0 0 +1 -1
];
levels,
labels=["do_so", "spt", "grv"],
),
)
fit(MixedModel, form, dat1; contrasts, progress)
end
Est. | SE | z | p | σ_Subj | |
(Intercept) | 389.7336 | 7.0937 | 54.94 | <1e-99 | 55.2184 |
CTR: do_so | 59.0052 | 5.1816 | 11.39 | <1e-29 | 36.1992 |
CTR: spt | 33.7817 | 3.2883 | 10.27 | <1e-24 | 23.2561 |
CTR: grv | 2.7469 | 2.2146 | 1.24 | 0.2148 | 9.5149 |
Residual | 69.8348 |
The three contrasts for one main effect and two nested contrasts are orthogonal. There is no test of the interaction (parallelism).
Three factors:
2 x 3 = 6 measures / subject
Arrow.Table with 396 rows, 5 columns, and schema:
:Subj String
:Group String
:A String
:B String
:dv Float64
We select an LMM supported by the data.
cntrst2 = Dict(
:Group => SeqDiffCoding(; levels=["G1", "G2"]),
:A => SeqDiffCoding(; levels=["A1", "A2"]),
:B => SeqDiffCoding(; levels=["B1", "B2", "B3"]),
)
f6_cpx = @formula dv ~ 1 + Group*A*B + (1 + A+B | Subj);
m6_cpx = fit(MixedModel, f6_cpx, dat2; contrasts=cntrst2)
issingular(m6_cpx)
f6_zcp = @formula dv ~ 1 + Group*A*B + zerocorr(1 + A+B | Subj);
m6_zcp = fit(MixedModel, f6_zcp, dat2; contrasts=cntrst2)
issingular(m6_zcp)
f6_ovi = @formula dv ~ 1 + Group*A*B + (1 | Subj);
m6_ovi = fit(MixedModel, f6_ovi, dat2; contrasts=cntrst2)
Est. | SE | z | p | σ_Subj | |
(Intercept) | 29.5404 | 1.0698 | 27.61 | <1e-99 | 8.5428 |
Group: G2 | -6.3605 | 2.1395 | -2.97 | 0.0030 | |
A: A2 | -0.0074 | 0.3879 | -0.02 | 0.9848 | |
B: B2 | 0.1856 | 0.4751 | 0.39 | 0.6960 | |
B: B3 | -0.9345 | 0.4751 | -1.97 | 0.0492 | |
Group: G2 & A: A2 | 0.2392 | 0.7758 | 0.31 | 0.7579 | |
Group: G2 & B: B2 | 1.9668 | 0.9502 | 2.07 | 0.0385 | |
Group: G2 & B: B3 | 0.0872 | 0.9502 | 0.09 | 0.9269 | |
A: A2 & B: B2 | 3.5336 | 0.9502 | 3.72 | 0.0002 | |
A: A2 & B: B3 | 0.5232 | 0.9502 | 0.55 | 0.5819 | |
Group: G2 & A: A2 & B: B2 | 2.8335 | 1.9004 | 1.49 | 0.1360 | |
Group: G2 & A: A2 & B: B3 | -2.3042 | 1.9004 | -1.21 | 0.2253 | |
Residual | 3.8580 |
Likelihood-ratio test: 3 models fitted on 396 observations
────────────────────────────────────────────────────────
DOF ΔDOF LogLik Deviance Chisq p(>Chisq)
────────────────────────────────────────────────────────
[1] 14 -1209.2528 2418.5057
[2] 17 3 -1208.8128 2417.6256 0.8800 0.8302
[3] 23 6 -1208.0644 2416.1289 1.4967 0.9597
────────────────────────────────────────────────────────
There is a significant interaction between A and the first contrast of B (i.e., B2 - B1). The interaction is not significant for A and the second contrast of B (i.e., B3 - B2). This implies that the left pair of lines in the following figure is statistically not parallel and that we do not have sufficient evidence that the right pair of lines is not parallel.
───────────────────────────────────────────────────────────────────
Coef. Std. Error z Pr(>|z|)
───────────────────────────────────────────────────────────────────
A: A2 & B: B2 3.53363 0.924789 3.82 0.0001
A: A2 & B: B3 0.523243 0.950202 0.55 0.5819
───────────────────────────────────────────────────────────────────
using Chain
tbl1 = @chain DataFrame(dat2) begin
groupby(_, [:Subj, :A, :B])
combine(_, nrow => :n, :dv => mean => :dv)
groupby(_, [:A, :B])
combine(_,
:dv => mean => :dv_M,
:dv => std => :dv_SD,
:dv => sem => :dv_SE)
end
fig1 = data(tbl1) * mapping(:B, :dv_M; color=:A) * (visual(Lines) + visual(Scatter))
draw(fig1)
StatsModels.jl provides a few commonly used contrast coding schemes, some less-commonly used schemes, and structs that allow you to manually specify your own, custom schemes.
The most commonly used contrasts are DummyCoding
and EffectsCoding
(which are similar to contr.treatment()
and contr.sum()
in R, respectively).
We also provide HelmertCoding
and SeqDiffCoding
(corresponding to base R’s contr.helmert()
and MASS::contr.sdif()
).
ContrastsCoding()
There are two ways to manually specify contrasts. First, you can specify them directly via ContrastsCoding
. If you do, it’s good practice to specify the levels corresponding to the rows of the matrix, although they can be omitted in which case they’ll be inferred from the data.
HypothesisCoding()
A better way to specify manual contrasts is via HypothesisCoding
, where each row of the matrix corresponds to the weights given to the cell means of the levels corresponding to each column (see Schad et al. (2020) for more information).