From 43bab38760aaaf5bf39a8da1872c54937fb6f5f1 Mon Sep 17 00:00:00 2001 From: Ajinkya Kokandakar Date: Sat, 22 Mar 2025 17:25:53 -0500 Subject: [PATCH 1/5] added StatsAPI.mss(m::LinearMixedModel) --- src/linearmixedmodel.jl | 2 ++ test/fit.jl | 2 ++ 2 files changed, 4 insertions(+) diff --git a/src/linearmixedmodel.jl b/src/linearmixedmodel.jl index 7c7b761b6..01ca9e317 100644 --- a/src/linearmixedmodel.jl +++ b/src/linearmixedmodel.jl @@ -1212,6 +1212,8 @@ function StatsAPI.stderror(m::LinearMixedModel{T}) where {T} return stderror!(similar(pivot(m), T), m) end +StatsAPI.mss(m::LinearMixedModel) = sum(abs.(mean(m.y) - fitted(m))) + """ updateA!(m::LinearMixedModel) diff --git a/test/fit.jl b/test/fit.jl index a057093e5..529c907f9 100644 --- a/test/fit.jl +++ b/test/fit.jl @@ -11,6 +11,8 @@ using Test @test deviance(m1) ≈ deviance(m2) @test isa(lmm(@formula(yield ~ 1 + (1|batch)), MixedModels.dataset(:dyestuff); progress=false, REML = true), LinearMixedModel) + @test_broken mss(m1) ≈ 30780.69 + # example from https://github.com/JuliaStats/MixedModels.jl/issues/194 # copied from tetst/pls.jl data = ( From b15ba2e76a29d59804939076515d06cabacd5872 Mon Sep 17 00:00:00 2001 From: Ajinkya Kokandakar Date: Sat, 29 Mar 2025 22:11:06 -0500 Subject: [PATCH 2/5] fixed MSS formula for linear mixed effects models --- src/linearmixedmodel.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/linearmixedmodel.jl b/src/linearmixedmodel.jl index 01ca9e317..c3507cf6f 100644 --- a/src/linearmixedmodel.jl +++ b/src/linearmixedmodel.jl @@ -1212,7 +1212,7 @@ function StatsAPI.stderror(m::LinearMixedModel{T}) where {T} return stderror!(similar(pivot(m), T), m) end -StatsAPI.mss(m::LinearMixedModel) = sum(abs.(mean(m.y) - fitted(m))) +StatsAPI.mss(m::LinearMixedModel) = sum(abs2.(mean(m.y) .- fitted(m))) """ updateA!(m::LinearMixedModel) From 2e3bb5c8aeb2e24bdf88a7ff9e0fc19aa7d21157 Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Sun, 30 Mar 2025 22:02:06 -0500 Subject: [PATCH 3/5] clean up mss definition a little --- src/MixedModels.jl | 11 ++++++----- src/linearmixedmodel.jl | 9 ++++++++- src/mixedmodel.jl | 2 ++ test/fit.jl | 4 +++- 4 files changed, 19 insertions(+), 7 deletions(-) diff --git a/src/MixedModels.jl b/src/MixedModels.jl index d0ec6df98..5a00b7e2c 100644 --- a/src/MixedModels.jl +++ b/src/MixedModels.jl @@ -30,10 +30,9 @@ using SparseArrays: nonzeros, nzrange, rowvals, sparse using StaticArrays: StaticArrays, SVector using Statistics: Statistics, mean, quantile, std using StatsAPI: StatsAPI, aic, aicc, bic, coef, coefnames, coeftable, confint, deviance -using StatsAPI: dof, dof_residual, fit, fit!, fitted, isfitted, islinear, leverage -using StatsAPI: - loglikelihood, meanresponse, modelmatrix, nobs, pvalue, predict, r2, residuals -using StatsAPI: response, responsename, stderror, vcov, weights +using StatsAPI: dof, dof_residual, fit, fit!, fitted, isfitted, islinear, leverage, mss +using StatsAPI: loglikelihood, meanresponse, modelmatrix, nobs, pvalue, predict, r2 +using StatsAPI: residuals, response, responsename, stderror, vcov, weights using StatsBase: StatsBase, CoefTable, model_response, summarystats using StatsFuns: log2π, normccdf using StatsModels: StatsModels, AbstractContrasts, AbstractTerm, CategoricalTerm @@ -41,6 +40,7 @@ using StatsModels: ConstantTerm, DummyCoding, EffectsCoding, FormulaTerm, Functi using StatsModels: HelmertCoding, HypothesisCoding, InteractionTerm, InterceptTerm using StatsModels: MatrixTerm, SeqDiffCoding, TableRegressionModel, Term using StatsModels: apply_schema, drop_term, formula, lrtest, modelcols, term, @formula +using StatsModels: hasintercept using StructTypes: StructTypes using Tables: Tables, columntable using TypedTables: TypedTables, DictTable, FlexTable, Table @@ -107,6 +107,7 @@ export @formula, fnames, GHnorm, glmm, + hasintercept, isfitted, islinear, issingular, @@ -118,6 +119,7 @@ export @formula, lowerbd, lrtest, meanresponse, + mss, modelmatrix, model_response, nobs, @@ -214,7 +216,6 @@ include("profile/profile.jl") include("nlopt.jl") include("prima.jl") - # aliases with non-unicode function names const settheta! = setθ! const profilesigma = profileσ diff --git a/src/linearmixedmodel.jl b/src/linearmixedmodel.jl index c3507cf6f..b72080af5 100644 --- a/src/linearmixedmodel.jl +++ b/src/linearmixedmodel.jl @@ -1212,7 +1212,14 @@ function StatsAPI.stderror(m::LinearMixedModel{T}) where {T} return stderror!(similar(pivot(m), T), m) end -StatsAPI.mss(m::LinearMixedModel) = sum(abs2.(mean(m.y) .- fitted(m))) +function StatsAPI.mss(m::LinearMixedModel) + hasintercept(m) && return sum(abs2, meanresponse(m) .- fitted(m)) + throw( + ArgumentError( + "Mean sum of squares is defined only for models with an intercept term." + ), + ) +end """ updateA!(m::LinearMixedModel) diff --git a/src/mixedmodel.jl b/src/mixedmodel.jl index 16618be7c..973709665 100644 --- a/src/mixedmodel.jl +++ b/src/mixedmodel.jl @@ -57,6 +57,8 @@ function StatsAPI.dof_residual(m::MixedModel) return nobs(m) - dof(m) end +StatsModels.hasintercept(m::MixedModel) = hasintercept(formula(m)) + """ issingular(m::MixedModel, θ=m.θ; atol::Real=0, rtol::Real=atol>0 ? 0 : √eps) diff --git a/test/fit.jl b/test/fit.jl index 529c907f9..01af5bae3 100644 --- a/test/fit.jl +++ b/test/fit.jl @@ -11,7 +11,9 @@ using Test @test deviance(m1) ≈ deviance(m2) @test isa(lmm(@formula(yield ~ 1 + (1|batch)), MixedModels.dataset(:dyestuff); progress=false, REML = true), LinearMixedModel) - @test_broken mss(m1) ≈ 30780.69 + @test mss(m1) ≈ 30780.69 atol=0.005 + m0 = fit(MixedModel, @formula(yield ~ 0 + (1|batch)), MixedModels.dataset(:dyestuff); progress=false) + @test_throws ArgumentError("Mean sum of squares is defined only for models with an intercept term.") mss(m0) # example from https://github.com/JuliaStats/MixedModels.jl/issues/194 # copied from tetst/pls.jl From eb75ca9897ccca771cec8c23ac073eed0db8ffaf Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Sun, 30 Mar 2025 22:16:42 -0500 Subject: [PATCH 4/5] tweak tolerance --- test/fit.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/fit.jl b/test/fit.jl index 01af5bae3..f9fa0faf5 100644 --- a/test/fit.jl +++ b/test/fit.jl @@ -11,7 +11,7 @@ using Test @test deviance(m1) ≈ deviance(m2) @test isa(lmm(@formula(yield ~ 1 + (1|batch)), MixedModels.dataset(:dyestuff); progress=false, REML = true), LinearMixedModel) - @test mss(m1) ≈ 30780.69 atol=0.005 + @test mss(m1) ≈ 30780.69 rtol=0.005 m0 = fit(MixedModel, @formula(yield ~ 0 + (1|batch)), MixedModels.dataset(:dyestuff); progress=false) @test_throws ArgumentError("Mean sum of squares is defined only for models with an intercept term.") mss(m0) From c7c365f2cccc2c6d4115288b620516b774af2fb2 Mon Sep 17 00:00:00 2001 From: Ajinkya Kokandakar Date: Sun, 30 Mar 2025 23:53:30 -0500 Subject: [PATCH 5/5] tweaked error message for MSS --- src/linearmixedmodel.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/linearmixedmodel.jl b/src/linearmixedmodel.jl index b72080af5..2bb3b10da 100644 --- a/src/linearmixedmodel.jl +++ b/src/linearmixedmodel.jl @@ -1216,7 +1216,7 @@ function StatsAPI.mss(m::LinearMixedModel) hasintercept(m) && return sum(abs2, meanresponse(m) .- fitted(m)) throw( ArgumentError( - "Mean sum of squares is defined only for models with an intercept term." + "Model sum of squares (MSS) is defined only for models with an intercept term." ), ) end