For multiple response linear models, it would be very nice to have support from equatiomatic. Here's a simple example with two response variables:
mod_eq <- lm(cbind(read, math) ~ income + educ,
data= heplots::NLSY)
extract_eq(mod_eq)
This gives:
$$
\operatorname{cbind(read,\ math)} = \alpha + \beta_{1}(\operatorname{income}) + \beta_{2}(\operatorname{educ}) +
\alpha + \beta_{4}(\operatorname{income}) + \beta_{5}(\operatorname{educ}) + \epsilon
$$
which renders as

If the parameters are recognized as 2 x 1 vectors, this could simply be given as
$$
\operatorname{(read,\ math)} = \boldsymbol{\alpha} +
\boldsymbol{\beta}_{1}(\operatorname{income}) +
\boldsymbol{\beta}_{2}(\operatorname{educ}) +
\boldsymbol{\epsilon}
$$

Note that heplots does supply a tidy.mlm method:
> tidy(mod_eq)
# A tibble: 6 × 6
response term estimate std.error statistic p.value
<chr> <chr> <dbl> <dbl> <dbl> <dbl>
1 read (Intercept) 15.9 4.26 3.73 0.000241
2 read income 0.0147 0.0327 0.450 0.653
3 read educ 0.947 0.361 2.62 0.00927
4 math (Intercept) 8.76 4.33 2.02 0.0443
5 math income 0.0919 0.0332 2.77 0.00613
6 math educ 1.27 0.367 3.45 0.000652
So one could sort this by term to extract the coefficients.
For multiple response linear models, it would be very nice to have support from
equatiomatic. Here's a simple example with two response variables:This gives:
which renders as
If the parameters are recognized as 2 x 1 vectors, this could simply be given as
Note that
heplotsdoes supply atidy.mlmmethod:So one could sort this by
termto extract the coefficients.