Categorical Predictors (ANOVA family)

Categorical Predictors (ANOVA family)#

In the previous tutorial we met the lm() class for estimating univariate and multiple regression models. Here we’ll build upon this by demonstrating how you can work with models that have categorical or factor predictors. All models include the following methods to make it easy to work with categorical predictors which we’ll explore in this tutorial:

Method

Description

.set_factors()

Set a variable as a factor with alphabetaically sorted levels by default, or specify the exact level ordering; this changes the data-type to categorical/enum

.unset_factors()

Unset all factor variables to their original datatype

.show_factors()

Show any set factors

.set_contrasts()

Encode factor levels using one of the supported coding schemes (contr.treatment by default) or a custom one using numeric codes that sum-to-zero

.show_contrasts()

Show any set contrasts

.anova()

Calculate a Type-III ANOVA table using joint_tests() in R

.summary_anova()

Get a nicely formatted table containing .result_anova

.emmeans()

Compute marginal means/contrasts at specific factor levels

Basics#

In the previous tutorial we estimated a model with a categorical predictor \(cyl\)

from pymer4.models import lm
from pymer4 import load_dataset
import numpy as np
from polars import col

mtcars = load_dataset('mtcars').select('mpg', 'wt', 'cyl')
model = lm('mpg ~ wt * cyl', data=mtcars)

We specified our factor variable and transformed a continuous variable

model.set_factors({'cyl': ['4', '6', '8']})
model.set_transforms({'wt': 'center'})

We saw that factor variables use treatment contrasts (dummy-coding) by default just like R. This means each parameter represents the difference between the first level and the nth level.

For \(cyl\), that would be \(6-4\) and \(8-4\) as \(4\) is the first level of \(cyl\) based on the order we provided to .set_factors().

model.show_contrasts()

{'cyl': 'contr.treatment'}

Using .set_* methods changes the underlying .data attribute of a model while the coresponding .unset_* methods undo those changes

model.data.head()
shape: (5, 4)
mpgwtcylwt_orig
f64f64enumf64
21.0-0.59725"6"2.62
21.0-0.34225"6"2.875
22.8-0.89725"4"2.32
21.4-0.00225"6"3.215
18.70.22275"8"3.44
# wt has been replaced by its original value previously
# stored in wt_orig
model.unset_transforms()
model.data.head()
shape: (5, 3)
mpgwtcyl
f64f64i64
21.02.626
21.02.8756
22.82.324
21.43.2156
18.73.448
# cyl is now back to an integer type
model.unset_factors()
model.data.head()
shape: (5, 4)
mpgwtcylwt_orig
f64f64i64f64
21.0-0.5972562.62
21.0-0.3422562.875
22.8-0.8972542.32
21.4-0.0022563.215
18.70.2227583.44

For factor variables there are several ways to work with .set_contrasts().
The first option is to use any of the default contrasts supported by R

# Mean-center weight again
model.set_transforms({'wt': 'center'})

model.set_factors({'cyl': ['4', '6', '8']})
model.set_contrasts({'cyl': 'contr.sum'})

model.show_contrasts()

{'cyl': 'contr.sum'}

Alternatively we can test a specific hypothesis by using a custom contrast across factor levels with

  • the order based on the level order set when using .set_factors() or alphabetically by default.

  • numerical codes that sum-to-zero which will be automatically converted to the contrast matrix expected by R

If we inspect the data we can see that there seems to be an approximately linear trend across levels of \(cyl\)

import seaborn as sns

grid = sns.FacetGrid(data=model.data, palette='Purples', aspect=2)

grid.map(sns.pointplot, 'cyl', 'mpg', color='gray', order=['4','6','8'], alpha=.25, errorbar=None, markersize=10)

grid.map_dataframe(sns.stripplot, 'cyl', 'mpg', order=['4','6','8'], hue='cyl', palette='Purples', legend=True);
../_images/ea322db105ca23fce4868430752b1c192c679246fc1b32d450df7b90eb24a8d0.png

We can test this by specifying a linear contrast across levels of cyl such that 4 < 6 < 8. The normalize argument simply scales our contrasts for us to match the polynomial used in R:

Now \(cyl1\) below reflects this specific comparison: a linear trend across levels
\(cyl2\) represents an automatically computed orthogonal comparison, in this case a quadratic trend across \(cyl\) levels.

model.set_contrasts({'cyl':[-1, 0, 1]}, normalize=True)
model.fit()
model.params
shape: (6, 2)
termestimate
strf64
"(Intercept)"19.22742
"wt"-3.539856
"cyl1"-3.244839
"cyl2"-0.290422
"wt:cyl1"2.442762
"wt:cyl2"-0.9305

We can inspect the .design_matrix after fitting the model to see how our contrast was converted for R and how the second orthogonal contrast was automatically calculated. The code below just grabs and renames the unique rows to make the display easier to see.

The first column is each level of \(cyl\).
The second column is our linear contrast [-1, 0, 1] converted to R’s format of contrast codes.
The third column is the orthogonal contrast auto-calculated for us, which in this case is a quadratic trend

(
    model.design_matrix
    .select('cyl1', 'cyl2')
    .unique(maintain_order=True)
    .with_columns(model.data['cyl'].unique(maintain_order=True))
    .sort('cyl')
    .select('cyl', 'cyl1', 'cyl2')
    .rename({'cyl': 'level', 'cyl1': 'contrast1', 'cyl2': 'contrast2'})
)
shape: (3, 3)
levelcontrast1contrast2
enumf64f64
"4"-0.7071070.408248
"6"0.0-0.816497
"8"0.7071070.408248

Our custom normalized linear contrast is equivalent to using the default contr.poly scheme, which models differences across factor levels as orthogonal polynomials.

Notice that while the parameter names have changes, the estimates are the same as above:

model.set_contrasts({'cyl': 'contr.poly'})
model.fit()
model.params
shape: (6, 2)
termestimate
strf64
"(Intercept)"19.22742
"wt"-3.539856
"cyl.L"-3.244839
"cyl.Q"-0.290422
"wt:cyl.L"2.442762
"wt:cyl.Q"-0.9305

Using .set_contrasts() has a direct influence on what a model’s parameter estimates represent. This will influence the output of what you see with .summary(), .params, and .result_fit

This is one way to work with categorical predictors, sometimes called “planned-comparisons”. Later in this tutorial we’ll see an alternative, more general approach using marginal effects estimation.

Common statistical tests are regression#

A variety of common statistical tests that involve comparisons between 2 or more groups, are just linear regression models with categorical predictors:

Often called

GLM predictor(s) \(X\)

Example

T-test

1 two-level categorical

y ~ group
e.g. group levels:
['a', 'b']

One-way ANOVA

1 three+ level categorical

y ~ group
e.g. group levels:
['a', 'b', 'c']

ANCOVA

1 three+ level categorical & 1+ continuous

y ~ group + covariate

Moderation

1 three+ level categorical & 1+ continuous & interaction

y ~ group * covariate

Factorial ANOVA

2 or more two+ level categoricals

y ~ group * category * color

2-levels (t-test)#

Let’s start with a simple example using a categorical variable \(Student\) with 2 levels to compare mean differences in \(Balance\)

\[ Balance \sim \beta_0 + \beta_1 * Student\]
Hide code cell source 
df = load_dataset("credit")

grid = sns.catplot(
    x="Student",
    y="Balance",
    kind="box",
    hue="Student",
    order=["Yes", "No"],
    aspect=0.75,
    height=4,
    data=df.to_pandas(),
)
grid.map(
    sns.stripplot, "Student", "Balance", color="black", alpha=0.25, order=["Yes", "No"]
).set_ylabels('Credit Balance');
../_images/1cebad84e251cb532d4d0c980ae359eb49262ddb612754f80e93dab153fbf821.png

After creating the model, we can use .set_factors() treat \(Student\) as a categorical variable. Its levels will be alphabetical by default

model = lm('Balance ~ Student', data=df)

model.set_factors('Student')

model.show_factors()

{'Student': ['No', 'Yes']}

This means that our parameter estimate \(\beta_1\) will reflect \(Student_{Yes} - Student_{No}\)

model.fit(summary=True)
Formula: lm(Balance~Student)
Number of observations: 400
Confidence intervals: parametric
---------------------
R-squared: 0.0671
R-squared-adj: 0.0647
F(1, 398) = 28.622, p = <.001
Log-likelihood: -3005
AIC: 6016 | BIC: 6028
Residual error: 444.626
Estimate SE CI-low CI-high T-stat df p
(Intercept) 480.369 23.434 434.300 526.439 20.499 398 <.001 ***
StudentYes 396.456 74.104 250.771 542.140 5.350 398 <.001 ***
Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1

The t-test on the \(\beta_1\) is equivalent to a two-sample t-test:

Hide code cell source 
# compare to scipy
from scipy.stats import ttest_ind

student = df.filter(col('Student') == 'Yes').select('Balance')
non_student = df.filter(col('Student') == 'No').select('Balance')

results = ttest_ind(student, non_student)
print(f"Scipy T-test:\ndiff = {student.mean().item() - non_student.mean().item():.3f} t= {results.statistic[0]:.2f} p = {results.pvalue[0]:.3f}")

Scipy T-test:
diff = 396.456 t= 5.35 p = 0.000

3-levels (one-way ANOVA)#

With more than 2 levels of a factor you may want to calculate an F-table, i.e. an ANOVA

penguins = load_dataset('penguins')
Hide code cell source 
grid = sns.catplot(
    x="species",
    y="body_mass_g",
    kind="box",
    hue="species",
    order=["Adelie", "Chinstrap", "Gentoo"],
    aspect=1,
    height=4,
    data=penguins.to_pandas(),
)
grid.map(
    sns.stripplot,
    "species",
    "body_mass_g",
    color="black",
    alpha=0.25,
    order=["Adelie", "Chinstrap", "Gentoo"],
).set_axis_labels("Species", "Body Mass (g)").figure.suptitle("Penguin Weights", x=.6, fontsize=14);
../_images/a5642103285ad929fe58bf33e7aa43df24aeece5ea4838ff6cdf05452d3d8b85.png

Let’s say we’re interested in the main-effect of \(Species\) on \(Body Mass\)

We can use just .set_factors() and call the .anova() method

model = lm('body_mass_g ~ species', data=penguins)
model.set_factors('species')
model.anova(summary=True)
ANOVA (Type III tests)
model term df1 df2 F_ratio p_value
species 2 339 343.626 <.001 ***
Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1

Using .anova() is the most convenient way to get omnibus tests on categorical variables (i.e. main effects) especially with more complicated designs. By default it will completely ignore any contrasts you have set (or default contrasts) and use the appropriate ones in the face of both balanced and unbalanced group designs (i.e. orthogonal type-III SS) via emmeans::joint_tests.

This makes it easy to inspect a summary of the hypothesis tests setup by your parameterization without worrying about the correct ANOVA coding-schemes for categorical predictors. Just set your factors with .set_factors() and call .anova()

# Unchanged by ANOVA and determines what we see in .summary()
model.show_contrasts()

{'species': 'contr.treatment'}

# Default treatment contrast with Adelie
# as the reference because alphabetical order
model.params
shape: (3, 2)
termestimate
strf64
"(Intercept)"3700.662252
"speciesChinstrap"32.425984
"speciesGentoo"1375.354009

Marginal Effects Estimation#

Previously we saw how using .set_contrasts() allows us to parameterize our model so we can test specific hypotheses about the levels of a factor variable (e.g. treatment-coding, linear-trend, etc). However, to test a different hypothesis requires re-fitting a new model with new contrasts.

An alternative approach is to compare predictions from a fitted model using predictors set/fixed to particular values. Aggregating and comparing these model-predictions is known as marginal effects estimation. It’s an incredibly powerful general-purpose approach that encompasses traditional techniques like “post-hoc tests” as well as modern approaches like counter-factual estimations. It also makes it easier to reason about more complicated generalized models with interactions, without struggling to potentially mis-interpret parameter estimates alone.

In pymer4 you can compare marginal estimates of each factor levels using .emmeans(), with the order based on how we used .set_factors() and generate arbitrary predictions and comparisons using .empredict()

Using a factor variable as the first argument to .emmeans() provides marginal estimates at each level of the factor while any other model parameters are held at their mean:

model.emmeans('species')
shape: (3, 6)
speciesemmeanSEdflower_CLupper_CL
catf64f64f64f64f64
"Adelie"3700.66225237.619354339.03610.3910843790.93342
"Chinstrap"3733.08823556.059001339.03598.5694063867.607064
"Gentoo"5076.0162641.681876339.04975.9966915176.035829

By using the contrasts argument, you can setup various comparisons of interest with multiple-comparisons correction. Using the 'pairwise' will perform all possible pairwise comparisons. In this case we have only 1 predictor with 3 levels:

model.emmeans('species', contrasts='pairwise', p_adjust='sidak')
shape: (3, 8)
contrastestimateSEdflower_CLupper_CLt_ratiop_value
strf64f64f64f64f64f64f64
"Adelie - Chinstrap"-32.42598467.511684339.0-194.426598129.574631-0.4803020.949888
"Adelie - Gentoo"-1375.35400956.147971339.0-1510.086328-1240.621689-24.4951690.0
"Chinstrap - Gentoo"-1342.92802569.856928339.0-1510.556273-1175.299777-19.2239780.0

We can also can test a specific comparison using numeric codes in the same order as factor levels

model.emmeans('species', contrasts={'adelie_vs_others': [1, -0.5, -0.5]})
shape: (1, 8)
contrastestimateSEdflower_CLupper_CLt_ratiop_value
strf64f64f64f64f64f64f64
"adelie_vs_others"-703.88999651.33433339.0-804.863928-602.916064-13.7118772.3870e-34

Compare this to the approach we used before with .set_contrasts(), which allowed us to parameterize a model to test specific comparisons, thus changing what parameter estimates represent.
Using .emmeans() compares model predictions, making it a more generally usefully approach that allows you to ask many questions from a single estimated model.

Parameterizing our model the same way as our marginal estimation above matches the estimated value:

# With auto-solve for other orthogonal contrasts
# speceis1 represents the same comparison as above
model.set_contrasts({'species': [1, -.5, -.5]})
model.fit()
model.params
shape: (3, 2)
termestimate
strf64
"(Intercept)"4169.922249
"species1"-703.889996
"species2"949.593513

Multiple predictors (factorial ANOVA)#

Hide code cell source 
poker = load_dataset('poker')
grid = sns.catplot(
    x="hand",
    y="balance",
    kind="box",
    hue="skill",
    order=["bad", "neutral", "good"],
    hue_order=['average', 'expert'],
    aspect=1,
    height=4,
    data=poker.to_pandas(),
)
grid.map_dataframe(
    sns.stripplot,
    "hand",
    "balance",
    hue='skill',
    palette='dark:black',
    dodge=True,
    alpha=0.25,
    order=["bad", "neutral", "good"],
    hue_order=['average', 'expert'],
).set_axis_labels("Hand Quality", "Money wond ($)").figure.suptitle("Poker Winnings", x=.5, fontsize=14);
../_images/dc8de98ae131af9c8a7c97a6979f34e5eb19f5fb9590b7e30f188d4f53f3dfaa.png

Let’s say we wanted estimate a 2x3 factorial ANOVA. Like before we can use .set_factors() and then call .anova(). Specifying the factor levels isn’t necessary (alphabetic by default), but makes it easier to explore the marginal estimates later.

poker = load_dataset('poker')

model = lm('balance ~ hand * skill', data=poker)

# Set factors
model.set_factors({'hand': ['bad', 'neutral','good'], 'skill':['average', 'expert']})

# Calculate ANOVA
model.anova(summary=True)
ANOVA (Type III tests)
model term df1 df2 F_ratio p_value
hand 2 294 79.169 <.001 ***
skill 1 294 2.434 0.1198
hand:skill 2 294 7.083 <.001 ***
Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1

More Marginal Effects Estimation#

Main effect of skill#

Like before, we can unpack the F-test results using comparison between marginal means.
The F-test result tells are there’s no main-effect of skill, which we can verify by calculating the difference or contrast between each level of skill while averaging over levels of hand.

Notice the p-value is the same and that \(F=t^2\)

skill_contrast = [1, -1]
model.emmeans('skill', contrasts={'avg_vs_expert': skill_contrast})

R messages: 
NOTE: Results may be misleading due to involvement in interactions
shape: (1, 8)
contrastestimateSEdflower_CLupper_CLt_ratiop_value
strf64f64f64f64f64f64f64
"avg_vs_expert"-0.7243330.464243294.0-1.6379940.189328-1.5602460.119778

We can verify the estimate matches what would have gotten had we calculated the ANOVA cell means ourselves:

average = model.data.filter(col('skill') == 'average').select('balance')
expert = model.data.filter(col('skill') == 'expert').select('balance')
average.mean() - expert.mean()
shape: (1, 1)
balance
f64
-0.724333

Main effect of hand#

The F-test for hand indicates a main effect but it doesn’t tell us why.
In the figure above, it looks like there’s an increasing linear trend as we go from bad -> neutral -> good.

Let’s test that, while averaging over levels of skill

# Specified in the order of the factor levels, which we set above
hand_contrast = [-1, 0, 1]
model.emmeans('hand', contrasts={'linear_trend':hand_contrast})

R messages: 
NOTE: Results may be misleading due to involvement in interactions
shape: (1, 8)
contrastestimateSEdflower_CLupper_CLt_ratiop_value
strf64f64f64f64f64f64f64
"linear_trend"7.08490.568579294.05.9658998.20390112.4607066.6673e-29

We can again verify this estimate matches what we would have gotten through manual calculation:

bad = model.data.filter(col("hand") == "bad").select("balance").mean().item()
neutral = model.data.filter(col("hand") == "neutral").select("balance").mean().item()
good = model.data.filter(col("hand") == "good").select("balance").mean().item()

np.dot([bad, neutral, good], hand_contrast)

np.float64(7.084900000000006)

Moderated by skill#

Both the printed message from the previous contrast, as well as the F-test on the interaction tells us that this trend is going to vary depending upon which level of skill we look at.
We can use by to subset out contrast for each level of skill to see how they differ

model.emmeans("hand", by="skill", contrasts={"linear_trend": hand_contrast})
shape: (2, 9)
contrastskillestimateSEdflower_CLupper_CLt_ratiop_value
catcatf64f64f64f64f64f64f64
"linear_trend""average"9.2110.804093294.07.62849310.79350711.4551482.2922e-25
"linear_trend""expert"4.95880.804093294.03.3762936.5413076.1669512.2960e-9

It does look like the trend is steeper for average players vs expert players. Let’s test that specific interaction by using interaction which allow us to declare a comparison for the by="skill" variable, in the same way that contrasts lets us decare a comparison for the "hand" variable

skill_contrast = [1, -1]

model.emmeans(
    "hand",
    by="skill",
    contrasts={"linear_trend": hand_contrast},
    interaction={"avg_vs_expert": skill_contrast},
)
shape: (1, 8)
contrastestimateSEdflower_CLupper_CLt_ratiop_value
strf64f64f64f64f64f64f64
"avg_vs_expert"4.25221.137159294.02.0141976.4902033.739320.000222

Since skill only has 2 levels, we could have used pairwise instead:

model.emmeans(
    "hand",
    by="skill",
    contrasts={"linear_trend": hand_contrast},
    interaction="pairwise",
)
shape: (1, 8)
contrastestimateSEdflower_CLupper_CLt_ratiop_value
strf64f64f64f64f64f64f64
"linear_trend average - linear_…4.25221.137159294.02.0141976.4902033.739320.000222

This is a little more work to calculate manually but matches our estimate:

Hide code cell source 
bad_a = (
    model.data.filter(col("hand") == "bad", col("skill") == "average")
    .select("balance")
    .mean()
    .item()
)
neutral_a = (
    model.data.filter(col("hand") == "neutral", col("skill") == "average")
    .select("balance")
    .mean()
    .item()
)
good_a = (
    model.data.filter(col("hand") == "good", col("skill") == "average")
    .select("balance")
    .mean()
    .item()
)
bad_e = (
    model.data.filter(col("hand") == "bad", col("skill") == "expert")
    .select("balance")
    .mean()
    .item()
)
neutral_e = (
    model.data.filter(col("hand") == "neutral", col("skill") == "expert")
    .select("balance")
    .mean()
    .item()
)
good_e = (
    model.data.filter(col("hand") == "good", col("skill") == "expert")
    .select("balance")
    .mean()
    .item()
)

# Toggle code above to see calculations
np.dot([bad_a, neutral_a, good_a], hand_contrast) - np.dot(
    [bad_e, neutral_e, good_e], hand_contrast
)

np.float64(4.252199999999995)

Comparison to re-parameterizing model#

Alternatively, we can parameterize our regression model such that our coefficients reflect the comparisons we want, so called “planned comparisons”. Rather than use .anova(), we’ll now want to use .fit() and inspect the \(k-1\) parameters for each \(k\) level categorical predictor.

Let’s change the default treatment-coding to the comparisons we care about:

# Reuse from before
# linear contrast across levels of hand 
# average - expert for skill
model.set_contrasts({'hand': hand_contrast, 'skill': skill_contrast})
model.show_contrasts()

{'hand': array([-1,  0,  1]), 'skill': array([ 1, -1])}

Now the estimates capture what we care about.
Since hand has 3 levels, we need 2 parameters to fully capture it

  • hand1 = linear contrast (bad < neutral < good) averaged over skill

  • hand2 = automatically computed quadratic contrast averaged over skill

Since skill had 2 levels, we just need 1 parameter to capture it

  • skill1 = average - expert difference averaged over hand

The first interaction parameter reflect linear trend difference we calculated above:

  • hand1:skill1 = linear contrast average skill - linear contrast expert skill

  • hand2:skill1 = automatically computed quadratic contrast average skill - quadratic contrast expert skill

model.fit(summary=True)
Formula: lm(balance~hand*skill)
Number of observations: 300
Confidence intervals: parametric
---------------------
R-squared: 0.3731
R-squared-adj: 0.3624
F(5, 294) = 34.988, p = <.001
Log-likelihood: -840
AIC: 1694 | BIC: 1720
Residual error: 4.02
Estimate SE CI-low CI-high T-stat df p
(Intercept) 9.771 0.232 9.315 10.228 42.096 294 <.001 ***
hand1 7.085 0.569 5.966 8.204 12.461 294 <.001 ***
hand2 −0.704 0.402 −1.496 0.087 −1.752 294 0.08083 .
skill1 −0.724 0.464 −1.638 0.189 −1.560 294 0.1198
hand1:skill1 4.252 1.137 2.014 6.490 3.739 294 <.001 ***
hand2:skill1 0.344 0.804 −1.238 1.927 0.428 294 0.6687
Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1

All pairwise comparisons#

Though not encouraged, we can also calculate the pairwise comparisons between every combination of hand and skill

model.emmeans(['hand', 'skill'], contrasts='pairwise', p_adjust='bonferroni')
shape: (15, 8)
contrastestimateSEdflower_CLupper_CLt_ratiop_value
strf64f64f64f64f64f64f64
"bad average - neutral average"-5.25720.804093294.0-7.636816-2.877584-6.5380534.1276e-9
"bad average - good average"-9.2110.804093294.0-11.590616-6.831384-11.4551483.4383e-24
"bad average - bad expert"-2.70980.804093294.0-5.089416-0.330184-3.370010.012781
"bad average - neutral expert"-6.26280.804093294.0-8.642416-3.883184-7.7886551.7546e-12
"bad average - good expert"-7.66860.804093294.0-10.048216-5.288984-9.5369618.7939e-18
……………………
"good average - neutral expert"2.94820.804093294.00.5685845.3278163.6664930.004376
"good average - good expert"1.54240.804093294.0-0.8372163.9220161.9181870.840836
"bad expert - neutral expert"-3.5530.804093294.0-5.932616-1.173384-4.4186450.00021
"bad expert - good expert"-4.95880.804093294.0-7.338416-2.579184-6.1669513.4440e-8
"neutral expert - good expert"-1.40580.804093294.0-3.7854160.973816-1.7483061.0

With continuous predictor(s) (e.g. ANCOVA/moderation)#

Hide code cell source 
grid = sns.lmplot(
    x="flipper_length_mm",
    y="body_mass_g",
    hue="species",
    hue_order=["Adelie", "Chinstrap", "Gentoo"],
    aspect=1,
    height=4,
    data=penguins.to_pandas(),
)
grid.set_axis_labels("Flipper Length", "Body Mass (g)").figure.suptitle("Penguin Weights", x=.6, fontsize=14);
../_images/87bcd955c245e9d0348046c9812f4affdeddd9db050f46519bcda31cc5869a22.png

Models with mixed continuous and factor models work exactly the same.
Here we explore whether there’s an interaction between species and flipper length, i.e. whether a penguin’s species changes the relationship between flipper length and body mass.

model = lm('body_mass_g ~ species * flipper_length_mm', data=penguins)

# Categorical
model.set_factors('species')
model.set_contrasts({'species':'contr.sum'})

# Continuous
model.set_transforms({'flipper_length_mm': 'center'})

# Calculate ANOVA
model.anova(summary=True)
ANOVA (Type III tests)
model term df1 df2 F_ratio p_value
species 2 336 3.934 0.02047 *
flipper_length_mm 1 336 168.292 <.001 ***
species:flipper_length_mm 2 336 5.532 0.004327 **
Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1

While we could try to understand these relationships just based on the estimated parameters, it’s much easier to unpack them using marginal effects estimation

model.params
shape: (6, 2)
termestimate
strf64
"(Intercept)"4052.294709
"species1"8.254155
"species2"-143.169979
"flipper_length_mm"40.675862
"species1:flipper_length_mm"-7.844172
"species2:flipper_length_mm"-6.102468

Even more Marginal Effects Estimation#

\(flipper\_length\_mm\) above is the slope of flipper_length_mm averaged over species since we used sum (deviation) coding

model.emmeans('flipper_length_mm')

R messages: 
NOTE: Results may be misleading due to involvement in interactions
shape: (1, 5)
flipper_length_mm_trendSEdflower_CLupper_CL
f64f64f64f64f64
40.6758623.135487336.034.50820446.84352

\(species1\) is the first parameter of our sum (deviation) coding which represents factor levels as a differences from the grand-mean, when flipper_length_mm is held at its mean.

species_means = model.emmeans('species').select('emmean')
grand_mean = species_means.mean().item()
adelie_mean = species_means[0]
adelie_mean - grand_mean

R messages: 
NOTE: Results may be misleading due to involvement in interactions
shape: (1, 1)
emmean
f64
8.254155

\(species2\) is the difference between the next factor level (Chinstrap) and the grand-mean

species_means[1] - grand_mean
shape: (1, 1)
emmean
f64
-143.169979

With marginal estimation we can also ask questions about our continuous predictor. For example, passing using flipper_length_mm as our first variable and species as our by variable we can see how the slope of \(flipper length\) varies across \(species\)

model.emmeans('flipper_length_mm', by='species')
shape: (3, 6)
speciesflipper_length_mm_trendSEdflower_CLupper_CL
catf64f64f64f64f64
"Adelie"32.831694.627184336.021.72783643.935543
"Chinstrap"34.5733946.348364336.019.33922549.807563
"Gentoo"54.6225025.173874336.042.20675767.038247

Like before, we can use the contrasts argument to specify a test across levels of our by variable, i.e. a test across the different slopes above.

Here we test whether the relationship between flipper_length and weight gets stronger across species: Adelie < Chinstrap < Gentoo, i.e. that the slope increases

model.emmeans("flipper_length_mm", by="species", contrasts={"linear": [-1, 0, 1]}, normalize=True)
shape: (1, 8)
contrastestimateSEdflower_CLupper_CLt_ratiop_value
strf64f64f64f64f64f64f64
"linear"15.4084314.908146336.05.75386525.0629973.1393590.001843

Counter-factual estimates & comparisons#

We can also calculate marginal means at arbitrary values of flipper_length_mm using the at keyword and passing in 1 or more values. Let’s test a hypothetical question: how would species differ in weight if we were to observe a very short flipper_length_mm of 150mm?

We never actually observe this value of flipper_length_mm but can we can use the model to generate a counter-factual comparisons to get an answer:

model.emmeans('species', at={'flipper_length_mm': 150}, contrasts='pairwise')

R messages: 
NOTE: Results may be misleading due to involvement in interactions
shape: (3, 8)
contrastestimateSEdflower_CLupper_CLt_ratiop_value
strf64f64f64f64f64f64f64
"Adelie - Chinstrap"240.103356348.902095336.0-597.1570411077.3637520.6881680.868765
"Adelie - Gentoo"982.82199396.284903336.031.8570661933.7869152.4800890.040319
"Chinstrap - Gentoo"742.718635456.726766336.0-353.2886371838.7259071.6261770.282721

For even more control we can use .empredict() to generate predictions for arbitrary combinations of our predictors. This allows us to ask very specific questions like:

How much would an Adelie penguin weigh if it had a flipper length of 210mm?

.empredict() (and .emmeans()) will automatically handle transformations from .set_transforms() and return values of the predictor on the transformed scale along with predictions and uncertainty

model.empredict({'flipper_length_mm': 210, 'species': 'Adelie'})
shape: (1, 7)
speciesflipper_length_mmpredictionSEdflower_CLupper_CL
catf64f64f64f64f64f64
"Adelie"9.0847954358.81804697.537889336.04166.9562024550.67989

Using the special 'data' value allows us to easily reference all the observed values for a predictor. For example, we can generate unit-level predictions using all observed flipper lengths only for Chinstrap and Gentoo penguins

model.empredict({'flipper_length_mm': 'data', 'species': ['Chinstrap', 'Gentoo']}).head()
shape: (5, 7)
speciesflipper_length_mmpredictionSEdflower_CLupper_CL
catf64f64f64f64f64f64
"Chinstrap"-19.9152053220.588515104.285847336.03015.4531033425.723926
"Gentoo"-19.9152053099.392229190.185526336.02725.2879063473.496551
"Chinstrap"-14.9152053393.45548476.869652336.03242.2490813544.661886
"Gentoo"-14.9152053372.504738164.781292336.03048.3717983696.637678
"Chinstrap"-5.9152053704.61602945.244783336.03615.6173073793.61475
Contents