Overview

In this section, we demonstrate how to use MatrixLM.jl with a simple example using simulated data.

Matrix Linear Models (MLMs) provide a simple yet robust multivariate framework for encoding relationships and groupings in high-throughput data. MLM's flexibility allows it to encode both categorical and continuous relationships, thereby enhancing the detection of associations between responses and predictors.

Within the scope of the matrix linear model framework, the model is articulated as follows:

\[Y = XBZ^T+E\]

Where

• $Y_{n \times m}$ is the response matrix
• $X_{n \times p}$ is the matrix for main predictors,
• $Z_{m \times q}$ denotes the response attributes matrix based on supervised knowledge,
• $E_{n \times m}$ is the error term,
• $B_{p \times q}$ is the matrix for main and interaction effects.

Data Generation

For simplicity, we assume that both the responses and the predictors are presented as dataframes.

Our dataset consists of a dataframe dfX, which includes 5 predictors (but will require p=7 coefficients, as we will see later). Among these predictors, two are categorical, and three are numerical, spread across n = 100 samples. We then consider a response dataframe Y composed of m = 250 responses. To simulate the Y data, we need to generate the matrices Z,B, and E.

The matrix Z provides information about the response population, which corresponds to the Y's columns $y_{i \in [1, 250]}$. The dimensions of this matrix are set at 250x10.

Given this setup, the coefficient matrix B is designed to have dimensions of 7×10. This matches the number of columns in X (the matrix coding the 5 predictors in the data frame) with the number of information categories in Z. Finally, we construct the noise matrix E containing the error terms. We generate this matrix as a normally distributed matrix (all values from N(0,σ=4)), adding a layer of randomness to our simulation.

using MatrixLM, DataFrames, Random, Plots, StatsModels, Statistics
using StableRNGs

rng = StableRNG(2026)

# Dimensions of matrices
n = 100
m = 250

# Number of groupings designed in the Z matrix
q = 10

# Generate data with two categorical variables and 3 numerical variables.
dfX = hcat(
    DataFrame(catvar1=string.(rand(rng, 0:1, n)), catvar2=rand(rng, ["A", "B", "C", "D"], n)),
    DataFrame(rand(rng, n,3), ["var3", "var4", "var5"])
    );

first(dfX, 5)

Let's use the function design_matrix() to get the predictor model matrix based on the formula expression, including all the variable terms.

# Convert dataframe to prediction matrix
X = design_matrix(@mlmformula(catvar1 + catvar2 + var3 + var4 + var5), dfX)

p = size(X)[2]; # 5 columns in dfX, but p=7 coefs bc catvar2 has 4 levels

X[1:5, :]

5×7 Matrix{Float64}:
 0.0  0.0  0.0  0.0  0.654691   0.84795    0.753778
 1.0  0.0  0.0  0.0  0.0219391  0.438705   0.612088
 0.0  1.0  0.0  0.0  0.771963   0.377601   0.77834
 1.0  0.0  0.0  1.0  0.382814   0.0603695  0.597588
 1.0  0.0  1.0  0.0  0.339513   0.272342   0.812023

We also have the option to specify contrast coding in our model. For a detailed understanding of how to implement contrast coding, please refer to the documentation for contrast coding with StatsModels.jl. This provides comprehensive instructions and examples.

# Convert dataframe to prediction matrix
my_ctrst = Dict(
             :catvar1 => DummyCoding(base = "0"),
             :catvar2 => DummyCoding(base = "A")
           )

X = design_matrix(@mlmformula(catvar1 + catvar2 + var3 + var4 + var5), dfX, my_ctrst);

X[1:5, :]

5×7 Matrix{Float64}:
 0.0  0.0  0.0  0.0  0.654691   0.84795    0.753778
 1.0  0.0  0.0  0.0  0.0219391  0.438705   0.612088
 0.0  1.0  0.0  0.0  0.771963   0.377601   0.77834
 1.0  0.0  0.0  1.0  0.382814   0.0603695  0.597588
 1.0  0.0  1.0  0.0  0.339513   0.272342   0.812023

Randomly generate some data for column covariates Z and the error matrix E:

Z = rand(rng, m,q);
E = randn(rng, n,m).*4;

Next, we will structure the coefficient matrix B following a specific pattern. This approach will facilitate a more effective visualization of the results in the subsequent steps:

# (p,q)
B = [
    2.0   3.0   4.0   5.0  6.0  7.0  8.0  0.0 0.5 -2.0;
    0.01  0.02  0.01  0.09 0.18 0.03 0.14 0.0 0.5 -2.0;
    -1.0  -0.5  0.02  0.49 1.1  2.0  5.0  0.0 0.5 -2.0;
    -0.01 0.02  -0.01 3.0  3.0  7.0  0.14 0.0 0.5 -2.0;
    0.0   0.0   0.0   0.0  3.0  3.0  3.0  3.0 0.5 -2.0;
    3.0   3.0   3.0   3.0  0.08 0.03 0.0  0.0 0.5 -2.0;
    0.01  0.0   3.0   3.0  3.0  3.0 0.04  0.0 0.5 -2.0;
];

Generate the response matrix Y:

Y = X*B*Z' + E;

Now, construct the RawData object consisting of the response variable Y and row/column predictors X and Z. All three matrices must be passed in as 2-dimensional arrays.

# Construct a RawData object
dat = RawData(Response(Y), Predictors(X, Z));

We can use the function show() to respectively display a readable summary of the matrices and dimensions stored in a RawData object.

show(dat)

RawData
  Response matrix Y: 100 × 250
  Design matrix X: 100 × 7
  Design matrix Z: 250 × 10
  X includes intercept: false
  Z includes intercept: false
  Preview of Y first row (first 3 columns): [14.4917, 11.2353, 16.684, …]
  Preview of X first row (first 3 columns): [0.0, 0.0, 0.0, …]
  Preview of Z first row (first 3 columns): [0.7882, 0.0619, 0.0065, …]

Model estimation

Least-squares estimates for matrix linear models can be obtained by running mlm. An object of type Mlm will be returned, with variables for the coefficient estimates (B), the coefficient variance estimates (varB), and the estimated variance of the errors (sigma). By default, mlm estimates both row and column main effects (X and Z intercepts), but this behavior can be suppressed by setting addXIntercept=false and/or addZIntercept=false. Column weights for Y and the target type for variance shrinkage^[1] can be optionally supplied to weights and targetType, respectively.

est = mlm(dat; addXIntercept=false, addZIntercept=false); # Model estimation

We can use the function show() to respectively display a readable summary of the matrices and dimensions stored in a Mlm object.

show(est)

Matrix linear model fit
  Observations: 100
  Responses: 250
  Row predictors: 7
  Column predictors: 10
  Coefficient matrix B: 7 × 10
  Residual covariance matrix sigma(Σ): 250 × 250
  Weighted fit: false
  Preview of B first row (first 3 columns): [2.2743, 3.1679, 3.9293, …]
  Preview of sigma(Σ) first row (first 3 columns): [18.6124, 1.7587, -0.3541, …]
  Covariance shrinkage: none

Model predictions and residuals

The coefficient estimates can be accessed using coef(est). Predicted values and residuals can be obtained by calling predict() and residuals(). By default, both of these functions use the same data used to fit the model. However, a new Predictors object can be passed into predict() as the newPredictors argument, and a new RawData object can be passed into residuals() as the newData argument. For convenience, fitted(est) will return the fitted values by calling predict with the default arguments.

To compare the estimated coefficients with the original matrix B, we will visualize the matrices using heatmaps. This graphical representation allows us to readily see differences and similarities between the two.

esti_coef = coef(est); # Get the coefficients of the model

plot(
    heatmap(B[end:-1:1, :],
            size = (800, 300)),
    heatmap(esti_coef[end:-1:1, :],
            size = (800, 300),
            clims = (-2, 8)),
    title = ["\$ \\mathbf{B}\$" "\$ \\mathbf{\\hat{B}}\$"]
)

Let's employ the same visualization method to compare the predicted values with the original Y response matrix. This allows us to gauge the accuracy of our model predictions.

preds = predict(est); # Prediction value

plot(
    heatmap(Y[end:-1:1, :],
            size = (800, 300)),
    heatmap(preds.Y[end:-1:1, :],
            size = (800, 300),
            # clims = (-2, 8)
            ),
    title = ["\$ \\mathbf{Y}\$" "\$ \\mathbf{\\hat{Y}}\$"]
)

The residuals() function, available in MatrixLM.jl, computes residuals for each observation, helping us evaluate the discrepancy between the model's predictions and the observed data.

resids = residuals(est);

plot(
    heatmap(resids[end:-1:1, :],
            size = (800, 300)),
    histogram(
        (reshape(resids,250*100,1)),
            grid  = false,
            label = "",
            size = (800, 300)),
    title = ["Residuals" "Distribution of the residuals"]
)

T-statistics and permutation test

The t-statistics for an Mlm object (defined as est.B ./ sqrt.(est.varB)) can be obtained by running t_stat. By default, t_stat does not return the corresponding t-statistics for any main effects that were estimated by mlm, but they will be returned if isMainEff=true.

tStats = t_stat(est);

Permutation p-values for the t-statistics can be computed by the mlm_perms function. mlm_perms calls the more general function perm_pvals and will run the permutations in parallel when possible. The illustrative example below runs only 500 permutations; a different number can be specified as the second argument (by modifying the nPerms value). By default, the function used to permute Y is shuffle_rows, which shuffles the rows for Y. Alternative functions for permuting Y, such as shuffle_cols, can be passed into the argument permFun. mlm_perms calls mlm and t_stat, so the user is free to specify keyword arguments for mlm or t_stat; by default, mlm_perms will call both functions using their default behavior.

nPerms = 500
tStats, pVals = mlm_perms(dat, nPerms, addXIntercept=false, addZIntercept=false);

plot(
    heatmap(tStats[end:-1:1, :],
            c = :bluesreds,
            clims = (-40, 40),
            size = (800, 300)),
    heatmap(-log.(pVals[end:-1:1, :]),
            grid = false,
            size = (800, 300)),
    title = ["T Statistics" "- Log(P-values)"]
)

White cells represent p-values of 0, not missing. These p-values should be reported as < 1/500 when using 500 permutations here, and more permutations would give more precision about how small these p-values are.

julia> pVals7×10 Matrix{Float64}:
 0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.532  0.018  0.0
 0.674  0.914  0.962  0.88   0.76   1.0    0.316  0.812  0.802  0.0
 0.0    0.038  0.998  0.468  0.0    0.0    0.0    0.9    0.682  0.0
 0.718  0.942  0.982  0.0    0.0    0.0    0.11   0.908  0.946  0.0
 0.866  0.944  0.982  0.848  0.0    0.0    0.0    0.0    0.912  0.0
 0.0    0.0    0.0    0.0    0.964  0.992  0.98   0.988  0.544  0.0
 0.94   0.934  0.0    0.0    0.0    0.0    0.934  0.722  0.072  0.0

Summary

Call the summary on a fitted MatrixLM object to inspect the estimates, standard errors, confidence intervals (CIs), and the associated row and column labels.

The method returns a dataframe with one row per coefficient and the following columns:

• row_term, col_term: index labels for the row-side (X) and column-side (Z) associated with each coefficient entry (column-major order).
• coef: coefficient estimate for each X–Z pair.
• std_error, t_stat, p_value: respectively, standard error, T-statistics, and p-values for testing.
• ci_lower, ci_upper: two-sided confidence interval bounds at the specified significance level.

The summary method has three keyword arguments:

• alpha: two-sided significance level for CIs and margins of error (default 0.05, 95% CI).
• permutation_test: when true, compute permutation-based t-stats/p-values; when false, use t-distribution.
• nPerms: number of permutations to run when permutation_test is true.

This summary table can be used to identify which row/column covariate combinations are significant and to estimate their effect sizes along with their uncertainty.

MatrixLM.summary(est, alpha = 0.05, permutation_test = false)

References