diff --git a/Project.toml b/Project.toml index 7cda4a0..abc03af 100644 --- a/Project.toml +++ b/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" diff --git a/README.md b/README.md index 7af0566..612bcd6 100644 --- a/README.md +++ b/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 diff --git a/docs/Project.toml b/docs/Project.toml index e43e58d..7b18cfb 100644 --- a/docs/Project.toml +++ b/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" diff --git a/docs/make.jl b/docs/make.jl index 1e0371c..5e8d762 100644 --- a/docs/make.jl +++ b/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", ] ) diff --git a/docs/src/application_to_real_data.md b/docs/src/application_to_real_data.md new file mode 100644 index 0000000..6ea7dc2 --- /dev/null +++ b/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. diff --git a/docs/src/example_ordinal_data.md b/docs/src/example_ordinal_data.md index e608c1a..0b46e6b 100644 --- a/docs/src/example_ordinal_data.md +++ b/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 diff --git a/docs/src/getting_started.md b/docs/src/getting_started.md index 5df6f98..cb48c30 100644 --- a/docs/src/getting_started.md +++ b/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,22 +54,22 @@ 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. @@ -75,7 +77,7 @@ 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,7 +117,7 @@ nothing #hide Generate the response matrix `Y`: -```@example +```@example Y = X*B*Z' + E; nothing #hide @@ -123,24 +125,40 @@ 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 @@ -148,7 +166,7 @@ The coefficient estimates can be accessed using `coef(est)`. Predicted values an 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) ``` diff --git a/docs/src/varShrinkage_example.md b/docs/src/varShrinkage_example.md index f6814af..085e5cf 100644 --- a/docs/src/varShrinkage_example.md +++ b/docs/src/varShrinkage_example.md @@ -48,7 +48,9 @@ correlations. ```@example varshrinkage using MatrixLM, LinearAlgebra, Random, Statistics -Random.seed!(1) +using StableRNGs + +rng = StableRNG(2026) # Matrix dimensions n = 100 @@ -57,8 +59,8 @@ p = 5 q = 4 # Row and column predictors -X = randn(n, p) -Z = randn(m, q) +X = randn(rng, n, p) +Z = randn(rng, m, q) # True coefficient matrix B = [ @@ -83,8 +85,8 @@ standard_deviations = sqrt.(variances) rho = 0.6 # Generate correlated, heteroskedastic errors -common_noise = randn(n, 1) -independent_noise = randn(n, m) +common_noise = randn(rng, n, 1) +independent_noise = randn(rng, n, m) E = ( sqrt(rho) .* common_noise .+ @@ -113,6 +115,14 @@ dat = RawData( 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 varshrinkage +show(dat) +``` + ## Model Estimation Without Variance Shrinkage Set the positional Boolean argument to `false` to estimate the residual covariance matrix without shrinkage. @@ -128,6 +138,13 @@ est_no_shrinkage = mlm( 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 varshrinkage +show(est_no_shrinkage) +``` + The fitted object contains the coefficient estimates in `B`, the estimated error covariance matrix in `sigma`, and the coefficient variance estimates in `varB`. ```@example varshrinkage @@ -147,9 +164,10 @@ est_shrinkage = mlm( addZIntercept=false ) -nothing #hide +show(est_shrinkage) ``` + The estimator separately shrinks 1. arctanh-transformed sample correlations toward their common mean, and diff --git a/src/MatrixLM.jl b/src/MatrixLM.jl index e0bb545..65eea81 100644 --- a/src/MatrixLM.jl +++ b/src/MatrixLM.jl @@ -43,10 +43,14 @@ module MatrixLM include("mlm.jl") export Mlm, mlm, t_stat - # # Estimate extraction + # Estimate extraction include("summary.jl") export confint, summary + # Show methods + include("show.jl") + export show + # Predictions and residuals include("predict.jl") export coef, predict, fitted, residuals diff --git a/src/show.jl b/src/show.jl new file mode 100644 index 0000000..2fe1a2d --- /dev/null +++ b/src/show.jl @@ -0,0 +1,87 @@ +function row_preview(mat::AbstractMatrix; max_cols::Int=3) + if size(mat, 1) == 0 || size(mat, 2) == 0 + return "[]" + end + + ncols = min(size(mat, 2), max_cols) + values = [string(mat[1, j]) for j in 1:ncols] + preview = "[" * join(values, ", ") * (size(mat, 2) > ncols ? ", …" : "") * "]" + return preview +end + +""" + show(io::IO, data::RawData) + +Display a readable summary of the matrices and dimensions stored in a +`RawData` object. +""" +function Base.show(io::IO, data::RawData) + println(io, "RawData") + println(io, " Response matrix Y: $(data.n) × $(data.m)") + println(io, " Design matrix X: $(data.n) × $(data.p)") + println(io, " Design matrix Z: $(data.m) × $(data.q)") + println(io, " X includes intercept: $(data.predictors.hasXIntercept)") + println(io, " Z includes intercept: $(data.predictors.hasZIntercept)") + println(io, " Preview of Y first row (first $(min(3, size(get_Y(data), 2))) columns): " * + row_preview(round.(get_Y(data), digits=4))) + println(io, " Preview of X first row (first $(min(3, size(get_X(data), 2))) columns): " * + row_preview(round.(get_X(data), digits=4))) + println(io, " Preview of Z first row (first $(min(3, size(get_Z(data), 2))) columns): " * + row_preview(round.(get_Z(data), digits=4))) +end + +# """ +# show(io::IO, model::Mlm) + +# Display a compact summary of a fitted matrix linear model. +# """ +# function Base.show(io::IO, model::Mlm) +# print( +# io, +# "Mlm($(size(model.B, 1)) × $(size(model.B, 2)) coefficient matrix)", +# ) +# end + +""" + show(io::IO, model::Mlm) + +Display a readable summary of a fitted matrix linear model. +""" +function Base.show(io::IO, model::Mlm) + data = model.data + + println(io, "Matrix linear model fit") + println(io, " Observations: $(data.n)") + println(io, " Responses: $(data.m)") + println(io, " Row predictors: $(data.p)") + println(io, " Column predictors: $(data.q)") + println( + io, + " Coefficient matrix B: " * + "$(size(model.B, 1)) × $(size(model.B, 2))", + ) + println( + io, + " Residual covariance matrix sigma(Σ): " * + "$(size(model.sigma, 1)) × $(size(model.sigma, 2))", + ) + println(io, " Weighted fit: $(model.weights !== nothing)") + println( + io, + " Preview of B first row (first $(min(3, size(model.B, 2))) columns): " * + row_preview(round.(model.B, digits=4)), + ) + println( + io, + " Preview of sigma(Σ) first row (first $(min(3, size(model.sigma, 2))) columns): " * + row_preview(round.(model.sigma, digits=4)), + ) + + if model.targetType === nothing + println(io, " Covariance shrinkage: none") + else + println(io, " Covariance shrinkage target: $(model.targetType)") + println(io, " Shrinkage coefficient: (correlation = $(round(model.lambda.correlation, + digits=4)), variance = $(round(model.lambda.variance, digits=4)))") + end +end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index dd0de4f..be8b24c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -28,4 +28,5 @@ using Test include("contr_test.jl") include("calc_sigma_test.jl") include("summary_test.jl") + include("show_test.jl") end diff --git a/test/show_test.jl b/test/show_test.jl new file mode 100644 index 0000000..c0ba75b --- /dev/null +++ b/test/show_test.jl @@ -0,0 +1,57 @@ + +@testset "row_preview" begin + empty_mat = zeros(0, 0) + single_row_mat = [1.0 2.0 3.0] + wide_mat = reshape(1.0:8.0, (2, 4)) + + @test MatrixLM.row_preview(empty_mat) == "[]" + @test MatrixLM.row_preview(single_row_mat) == "[1.0, 2.0, 3.0]" + @test MatrixLM.row_preview(wide_mat) == "[1.0, 3.0, 5.0, …]" + @test MatrixLM.row_preview(wide_mat; max_cols=2) == "[1.0, 3.0, …]" +end + +@testset "show methods" begin + Y = randn(10, 5) + X = hcat(ones(10), randn(10)) + Z = hcat(ones(5), randn(5)) + + response = Response(Y) + predictors = Predictors(X, Z, true, true) + data = RawData(response, predictors) + + detailed_data = sprint(show, MIME"text/plain"(), data) + + @test occursin("RawData", detailed_data) + @test occursin("Response matrix Y: 10 × 5", detailed_data) + @test occursin("Design matrix X: 10 × 2", detailed_data) + @test occursin("Design matrix Z: 5 × 2", detailed_data) + @test occursin("X includes intercept: true", detailed_data) + @test occursin("Z includes intercept: true", detailed_data) + @test occursin("Preview of Y first row", detailed_data) + @test occursin("Preview of X first row", detailed_data) + @test occursin("Preview of Z first row", detailed_data) + + model = mlm(data) + + detailed_model = sprint(show, MIME"text/plain"(), model) + + @test occursin("Matrix linear model fit", detailed_model) + @test occursin("Observations: 10", detailed_model) + @test occursin("Responses: 5", detailed_model) + @test occursin("Row predictors: 2", detailed_model) + @test occursin("Column predictors: 2", detailed_model) + @test occursin("Coefficient matrix B: 2 × 2", detailed_model) + @test occursin("Residual covariance matrix sigma(Σ): 5 × 5", detailed_model) + @test occursin("Weighted fit: false", detailed_model) + @test occursin("Preview of B first row", detailed_model) + @test occursin("Preview of sigma(Σ) first row", detailed_model) + @test occursin("Covariance shrinkage: none", detailed_model) + + model_shrunk = mlm(data; targetType = "R") + detailed_model_shrunk = sprint(show, MIME"text/plain"(), model_shrunk) + + @test occursin("Covariance shrinkage target: R", detailed_model_shrunk) + @test occursin("Shrinkage coefficient:", detailed_model_shrunk) + @test occursin("correlation =", detailed_model_shrunk) + @test occursin("variance =", detailed_model_shrunk) +end \ No newline at end of file