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 7c7b761b6..2bb3b10da 100644
--- a/src/linearmixedmodel.jl
+++ b/src/linearmixedmodel.jl
@@ -1212,6 +1212,15 @@ function StatsAPI.stderror(m::LinearMixedModel{T}) where {T}
     return stderror!(similar(pivot(m), T), m)
 end
 
+function StatsAPI.mss(m::LinearMixedModel)
+    hasintercept(m) && return sum(abs2, meanresponse(m) .- fitted(m))
+    throw(
+        ArgumentError(
+            "Model sum of squares (MSS) 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 a057093e5..f9fa0faf5 100644
--- a/test/fit.jl
+++ b/test/fit.jl
@@ -11,6 +11,10 @@ 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 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)
+
     # example from https://github.com/JuliaStats/MixedModels.jl/issues/194
     # copied from tetst/pls.jl
     data = (