Improve documentation, display methods, and real-data examples

Project.toml

@@ -1,8 +1,11 @@
 name = "MatrixLM"
 uuid = "37290134-6146-11e9-0c71-a5c489be1f53"
-version = "1.0.0"
+version = "1.1.0"
 authors = ["Jane Liang", "Gregory Farage", "Chenhao Zhao", "Harsh Vardhan Dubey", "Saunak Sen"]
 
+[workspace]
+projects = ["docs"]
+
 [deps]
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"

README.md

@@ -2,10 +2,10 @@
 
 [![CI](https://github.com/senresearch/MatrixLM.jl/actions/workflows/ci.yml/badge.svg)](https://github.com/senresearch/MatrixLM.jl/actions/workflows/ci.yml)
 [![codecov](https://codecov.io/gh/senresearch/MatrixLM.jl/branch/main/graph/badge.svg?token=uHM6utUQoi)](https://codecov.io/gh/senresearch/MatrixLM.jl)
-[![GNU Affero license](https://img.shields.io/badge/license%20-GNU%20AFFERO%20GENERAL%20PUBLIC%20-green)](https://github.com/senresearch/MatrixLM.jl/blob/main/LICENSE.md)
 [![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://senresearch.github.io/MatrixLM.jl/stable)
 [![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://senresearch.github.io/MatrixLM.jl/dev) 
 [![Pkg Status](https://www.repostatus.org/badges/latest/active.svg)](https://www.repostatus.org/#active)
+[![GNU Affero license](https://img.shields.io/badge/license%20-GNU%20AFFERO%20GENERAL%20PUBLIC%20-green)](https://github.com/senresearch/MatrixLM.jl/blob/main/LICENSE.md)
 
 ## Description
 

docs/Project.toml

@@ -3,7 +3,10 @@ DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 MatrixLM = "37290134-6146-11e9-0c71-a5c489be1f53"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 StatsModels = "3eaba693-59b7-5ba5-a881-562e759f1c8d"
 
 [compat]
-Documenter = "1"
+DataFrames = "1.8.1"
+Documenter = "1.17.0"
+StatsModels = "0.7.7"

docs/make.jl

@@ -1,3 +1,11 @@
+using Pkg
+
+# Build docs in the docs-specific environment so documentation-only
+# dependencies (for example StableRNGs used in examples) are available.
+Pkg.activate(@__DIR__)
+Pkg.develop(PackageSpec(path=joinpath(@__DIR__, "..")))
+Pkg.instantiate()
+
 using MatrixLM
 using Documenter
 
@@ -11,6 +19,7 @@ makedocs(; modules=[MatrixLM], sitename="MatrixLM.jl", pages=[
         "Getting Started" => "getting_started.md",
         "Example: MLM for ordinal predictors" => "example_ordinal_data.md",
         "Example: Variance Shrinkage with MatrixLM.jl" => "varShrinkage_example.md",
+        "Example: Application to real data" => "application_to_real_data.md",
         "Types and Functions" => "functions.md",
     ]
 )

docs/src/application_to_real_data.md

@@ -0,0 +1,46 @@
+# Applications to Real Data
+
+`MatrixLM.jl` has been used to analyze real datasets in metabolomics and high-throughput chemical genetic screening. The repositories below provide code, data-processing workflows, and examples associated with published applications of matrix linear models.
+
+## Metabolomics
+
+The [`mlm-metabolomics-supplement`](https://github.com/senresearch/mlm-metabolomics-supplement) repository contains tools, scripts, and notebooks for preprocessing, analyzing, and visualizing metabolomics data from three studies:
+
+* Statin-Associated Muscle Symptoms (SAMS)
+* Pansteatitis in Mozambique tilapia
+* Chronic Obstructive Pulmonary Disease (COPD)
+
+Each study is organized in a separate directory containing descriptions of the corresponding data-processing and analysis workflow.
+
+These applications demonstrate how matrix linear models can relate individual-level characteristics, such as disease or treatment status, to metabolite measurements while incorporating metabolite-level information, such as chemical composition or metabolite class.
+
+**Reference**
+
+Farage, G., Zhao, C., Choi, H. Y., et al. (2025). Matrix Linear Models for Connecting Metabolite Composition to Individual Characteristics. *Metabolites*, 15(2), 140.
+https://doi.org/10.3390/metabo15020140
+
+## High-Throughput Chemical Genetic Screens
+
+The [`mlm_gs_supplement`](https://github.com/senresearch/mlm_gs_supplement) repository contains code used to reproduce the analyses presented in *Matrix Linear Models for High-Throughput Chemical Genetic Screens*.
+
+The study applies matrix linear models to an *Escherichia coli* chemical genetic screen to examine relationships between mutant strains and experimental conditions. The analyses were primarily performed in Julia, while visualizations were produced in R.
+
+Instructions for obtaining and preparing the data are provided in the supplemental repository.
+
+**Reference**
+
+Liang, J. W., Nichols, R. J., & Sen, Ś. (2019). Matrix Linear Models for High-Throughput Chemical Genetic Screens. *Genetics*, 212(4), 1063–1073.
+https://doi.org/10.1534/genetics.119.302299
+
+## Using These Examples
+
+These repositories contain complete research workflows rather than small introductory examples. Users who are new to `MatrixLM.jl` may wish to begin with the package tutorials before consulting these applications.
+
+The supplemental repositories are useful for understanding how matrix linear models can be incorporated into larger workflows involving:
+
+* data preprocessing;
+* construction of sample- and response-level design matrices;
+* model fitting;
+* extraction of coefficients and p-values;
+* visualization and interpretation of results; and
+* saving analysis outputs.

docs/src/example_ordinal_data.md

@@ -36,7 +36,9 @@ This matrix is generated as a normally distributed matrix
 
 ```@example
 using MatrixLM, DataFrames, Random, Plots, StatsModels, Statistics
-Random.seed!(1)
+using StableRNGs
+
+rng = StableRNG(2026)
 
 # Dimensions of matrices 
 n = 100
@@ -46,10 +48,10 @@ m = 250
 q = 4
 
 # Generate data with 1 ordinal categorical variable.
-dfX = DataFrame(catvar=rand(1:5, n));
+dfX = DataFrame(catvar=rand(rng, 1:5, n));
 levels_catvar = sort(unique(dfX.catvar));
 
-nothing #hide
+first(dfX, 5)
 ```
 
 To derive the predictor design matrix, we employ the `design_matrix()` function and apply contrast coding using the `SeqDiffCoding()` system. This particular coding system is particularly useful for testing hypotheses related to "sequential differences" between the levels of our ordinal predictor.
@@ -62,7 +64,7 @@ X_ctrst = Dict(
 X = design_matrix(@mlmformula(1 + catvar), dfX, X_ctrst);
 p = size(X, 2);
 
-nothing #hide
+X[1:5, :]
 ```
 
 The design matrix `X` has `p = 5` columns defined as:
@@ -76,14 +78,14 @@ We randomly generate a dataframe `Z` that provides information about whether a r
 The `FullDummyCoding` system generates one indicator (1 or 0) column for each level, including the base level. This technique is sometimes called one-hot encoding and is widely used for categorical variables.
 
 ```@example
-dfZ = DataFrame(attribute= rand(["A", "B", "C", "D"], m))
+dfZ = DataFrame(attribute= rand(rng, ["A", "B", "C", "D"], m))
 Z_ctrst = Dict(
              :attribute => StatsModels.FullDummyCoding(),
           )
 
 Z = design_matrix(@mlmformula(0 + attribute), dfZ, Z_ctrst);
 
-nothing #hide
+Z[1:5, :]
 ```
 
 The design matrix `Z` has `q = 4` columns defined as:
@@ -95,7 +97,7 @@ Z_names = MatrixLM.design_matrix_names(@mlmformula(0 + attribute), dfZ, Z_ctrst)
 The error matrix `E` is obtained as follows:
 
 ```@example
-E = randn(n,m).*4;
+E = randn(rng, n,m).*4;
 
 nothing #hide
 ```
@@ -135,6 +137,13 @@ dat = RawData(Response(Y), Predictors(X, Z, true, false));
 nothing #hide
 ```
 
+We can use the function `show()` to respectively display 
+a readable summary of the matrices and dimensions stored in a `RawData` object.
+
+```@example
+show(dat)
+```
+
 ## Model estimation
 
 
@@ -145,6 +154,14 @@ est = mlm(dat; addXIntercept=false, addZIntercept=false); # Model estimation
 nothing # hide
 ```
 
+We can use the function `show()` to respectively display 
+a readable summary of the matrices and dimensions stored in a `Mlm` object.
+
+```@example
+show(est)
+```
+
+
 ## Model predictions and residuals
 
 

docs/src/getting_started.md

@@ -39,9 +39,11 @@ We generate this matrix as a normally distributed matrix
 (all values from N(0,σ=4)), adding a layer of randomness to our simulation.
 
 
-```@example
+```@example 
 using MatrixLM, DataFrames, Random, Plots, StatsModels, Statistics
-Random.seed!(1)
+using StableRNGs
+
+rng = StableRNG(2026)
 
 # Dimensions of matrices 
 n = 100
@@ -52,30 +54,30 @@ q = 10
 
 # Generate data with two categorical variables and 3 numerical variables.
 dfX = hcat(
-    DataFrame(catvar1=string.(rand(0:1, n)), catvar2=rand(["A", "B", "C", "D"], n)), 
-    DataFrame(rand(n,3), ["var3", "var4", "var5"])
+    DataFrame(catvar1=string.(rand(rng, 0:1, n)), catvar2=rand(rng, ["A", "B", "C", "D"], n)), 
+    DataFrame(rand(rng, n,3), ["var3", "var4", "var5"])
     );
 
-nothing #hide
+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.
 
-```@example
+```@example 
 # 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
 
-nothing #hide
+X[1:5, :]
 ```
 
 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](https://juliastats.org/StatsModels.jl/stable/contrasts/#How-to-specify-contrast-coding).
 This provides comprehensive instructions and examples.
 
-```@example
+```@example 
 # Convert dataframe to prediction matrix
 my_ctrst = Dict(
              :catvar1 => DummyCoding(base = "0"),
@@ -84,21 +86,21 @@ my_ctrst = Dict(
 
 X = design_matrix(@mlmformula(catvar1 + catvar2 + var3 + var4 + var5), dfX, my_ctrst);
 
-nothing #hide
+X[1:5, :]
 ```
 
 Randomly generate some data for column covariates `Z` and the error matrix `E`:
 
-```@example
-Z = rand(m,q);
-E = randn(n,m).*4;
+```@example 
+Z = rand(rng, m,q);
+E = randn(rng, n,m).*4;
 
 nothing #hide
 ```
 
 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:
 
-```@example
+```@example 
 # (p,q)
 B = [
     2.0   3.0   4.0   5.0  6.0  7.0  8.0  0.0 0.5 -2.0;
@@ -115,40 +117,56 @@ nothing #hide
 
 Generate the response matrix `Y`:
 
-```@example
+```@example 
 Y = X*B*Z' + E;
 
 nothing #hide
 ```
 
 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.
 
-```@example
+```@example 
 # Construct a RawData object
 dat = RawData(Response(Y), Predictors(X, Z));
 
 nothing #hide
 ```
 
+We can use the function `show()` to respectively display 
+a readable summary of the matrices and dimensions stored in a `RawData` object.
+
+```@example
+show(dat)
+```
+
+
 ## 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. 
 
-```@example
+```@example 
 est = mlm(dat; addXIntercept=false, addZIntercept=false); # Model estimation
 
 nothing #hide
 ```
 
+We can use the function `show()` to respectively display 
+a readable summary of the matrices and dimensions stored in a `Mlm` object.
+
+```@example 
+show(est)
+```
+
+
 ## 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.
 
-```@example
+```@example 
 esti_coef = coef(est); # Get the coefficients of the model
 
 plot(
@@ -164,7 +182,7 @@ plot(
 
 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.
 
-```@example
+```@example 
 preds = predict(est); # Prediction value
 
 plot(
@@ -182,7 +200,7 @@ 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.
 
-```@example
+```@example 
 resids = residuals(est);
 
 plot(
@@ -202,7 +220,7 @@ plot(
 
 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`.
 
-```@example
+```@example 
 tStats = t_stat(est);
 
 nothing #hide
@@ -214,7 +232,7 @@ 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.
 
-```@example
+```@example 
 nPerms = 500
 tStats, pVals = mlm_perms(dat, nPerms, addXIntercept=false, addZIntercept=false);
 
@@ -257,7 +275,7 @@ The `summary` method has three keyword arguments:
 
 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.
 
-```@example
+```@example 
 MatrixLM.summary(est, alpha = 0.05, permutation_test = false)
 ```