diff --git a/NEWS.md b/NEWS.md
index 7980035c7..67cbb1e66 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,4 +1,6 @@
 * New kwarg `amalgamate` can be used to disable amalgation of random0effects terms sharing a single grouping variable. Generally, `amalgamate=false` will result in a slower fit, but may improve convergence in some pathological cases. Note that this feature is experimental and changes to it are **not** considered breakin. [#673]
+* More informative error messages when passing a `Distribution` or `Link` type instead of the desired instance. [#698]
+* More informative error message on the intentional decision not to define methods for the coefficient of determination. [#698]
 
 MixedModels v4.16.0 Release Notes
 ==============================
@@ -449,4 +451,5 @@ Package dependencies
 [#681]: https://github.com/JuliaStats/MixedModels.jl/issues/681
 [#682]: https://github.com/JuliaStats/MixedModels.jl/issues/682
 [#694]: https://github.com/JuliaStats/MixedModels.jl/issues/694
+[#698]: https://github.com/JuliaStats/MixedModels.jl/issues/698
 [#703]: https://github.com/JuliaStats/MixedModels.jl/issues/703
diff --git a/src/generalizedlinearmixedmodel.jl b/src/generalizedlinearmixedmodel.jl
index 01d0b94ef..dacefefc6 100644
--- a/src/generalizedlinearmixedmodel.jl
+++ b/src/generalizedlinearmixedmodel.jl
@@ -328,6 +328,26 @@ end
 
 StatsAPI.fitted(m::GeneralizedLinearMixedModel) = m.resp.mu
 
+function GeneralizedLinearMixedModel(
+    f::FormulaTerm,
+    tbl,
+    d::Type,
+    args...;
+    kwargs...,
+)
+    throw(ArgumentError("Expected a Distribution instance (`$d()`), got a type (`$d`)."))
+end
+
+function GeneralizedLinearMixedModel(
+    f::FormulaTerm,
+    tbl,
+    d::Distribution,
+    l::Type;
+    kwargs...
+)
+    throw(ArgumentError("Expected a Link instance (`$l()`), got a type (`$l`)."))
+end
+
 function GeneralizedLinearMixedModel(
     f::FormulaTerm,
     tbl,
diff --git a/src/mixedmodel.jl b/src/mixedmodel.jl
index 12054c07a..ff4e6c7c2 100644
--- a/src/mixedmodel.jl
+++ b/src/mixedmodel.jl
@@ -74,6 +74,28 @@ issingular(m::GeneralizedLinearMixedModel, θ=m.optsum.final) = any(lowerbd(m) .
 # FIXME: better to base this on m.optsum.returnvalue
 StatsAPI.isfitted(m::MixedModel) = m.optsum.feval > 0
 
+function StatsAPI.fit(
+    ::Type{<:MixedModel},
+    f::FormulaTerm,
+    tbl,
+    d::Type,
+    args...;
+    kwargs...
+)
+    throw(ArgumentError("Expected a Distribution instance (`$d()`), got a type (`$d`)."))
+end
+
+function StatsAPI.fit(
+    ::Type{<:MixedModel},
+    f::FormulaTerm,
+    tbl,
+    d::Distribution,
+    l::Type;
+    kwargs...
+)
+    throw(ArgumentError("Expected a Link instance (`$l()`), got a type (`$l`)."))
+end
+
 StatsAPI.meanresponse(m::MixedModel) = mean(m.y)
 
 """
@@ -97,6 +119,24 @@ function retbl(mat, trm)
 )
 end
 
+StatsAPI.adjr2(m::MixedModel) = r2(m)
+
+function StatsAPI.r2(m::MixedModel)
+    @error (
+        """There is no uniquely defined coefficient of determination for mixed models
+         that has all the properties of the corresponding value for classical 
+         linear models. The GLMM FAQ provides more detail:
+         
+         https://bbolker.github.io/mixedmodels-misc/glmmFAQ.html#how-do-i-compute-a-coefficient-of-determination-r2-or-an-analogue-for-glmms
+
+
+         Alternatively, MixedModelsExtras provides a naive implementation, but 
+         the warnings there and in the FAQ should be taken seriously!
+         """
+    )
+    throw(MethodError(r2, (m,)))
+end
+
 """
     raneftables(m::MixedModel; uscale = false)
 
diff --git a/test/pirls.jl b/test/pirls.jl
index 77a70c4a6..b50f0fa1a 100644
--- a/test/pirls.jl
+++ b/test/pirls.jl
@@ -18,6 +18,20 @@ include("modelcache.jl")
     @test MixedModel(f, d, Bernoulli(), ProbitLink()) isa GeneralizedLinearMixedModel
 end
 
+@testset "Type for instance" begin
+    vaform = @formula(r2 ~ 1 + anger + gender + btype + situ + (1|subj) + (1|item))
+    verbagg = dataset(:verbagg)
+    @test_throws ArgumentError fit(MixedModel, vaform, verbagg, Bernoulli, LogitLink)
+    @test_throws ArgumentError fit(MixedModel, vaform, verbagg, Bernoulli(), LogitLink)
+    @test_throws ArgumentError fit(MixedModel, vaform, verbagg, Bernoulli, LogitLink())
+    @test_throws ArgumentError fit(GeneralizedLinearMixedModel, vaform, verbagg, Bernoulli, LogitLink)
+    @test_throws ArgumentError fit(GeneralizedLinearMixedModel, vaform, verbagg, Bernoulli(), LogitLink)
+    @test_throws ArgumentError fit(GeneralizedLinearMixedModel, vaform, verbagg, Bernoulli, LogitLink())
+    @test_throws ArgumentError GeneralizedLinearMixedModel(vaform, verbagg, Bernoulli, LogitLink)
+    @test_throws ArgumentError GeneralizedLinearMixedModel(vaform, verbagg, Bernoulli(), LogitLink)
+    @test_throws ArgumentError GeneralizedLinearMixedModel(vaform, verbagg, Bernoulli, LogitLink())
+end
+
 @testset "contra" begin
     contra = dataset(:contra)
     thin = 5
diff --git a/test/pls.jl b/test/pls.jl
index 19ec323eb..4902555e3 100644
--- a/test/pls.jl
+++ b/test/pls.jl
@@ -1,3 +1,4 @@
+using GLM # bring r2 into scope
 using LinearAlgebra
 using MixedModels
 using PooledArrays
@@ -688,3 +689,17 @@ end
     # but this is surprisingly hard to trigger in a reliable way across platforms
     # just because of the vagaries of floating point.
 end
+
+@testset "methods we don't define" begin
+    m = first(models(:sleepstudy))
+    for f in [r2, adjr2]
+        @test_logs (:error,) begin
+            try
+                f(m)
+            catch
+                # capture the error, do nothing
+            end
+        end
+        @test_throws MethodError @suppress f(m)
+    end
+end