diff --git a/NEWS.md b/NEWS.md
index af3463200..e8b325165 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,3 +1,7 @@
+MixedModels v4.16.0 Release Notes
+==============================
+* Support for check tolerances in deserialization. [#703]
+
 MixedModels v4.15.0 Release Notes
 ==============================
 * Support for different optimization criteria during the bootstrap. [#694]
@@ -442,3 +446,4 @@ 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
+[#703]: https://github.com/JuliaStats/MixedModels.jl/issues/703
diff --git a/Project.toml b/Project.toml
index 1078a09f3..498deacd3 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "MixedModels"
 uuid = "ff71e718-51f3-5ec2-a782-8ffcbfa3c316"
 author = ["Phillip Alday <me@phillipalday.com>", "Douglas Bates <dmbates@gmail.com>", "Jose Bayoan Santiago Calderon <jbs3hp@virginia.edu>"]
-version = "4.15.0"
+version = "4.16.0"
 
 [deps]
 Arrow = "69666777-d1a9-59fb-9406-91d4454c9d45"
diff --git a/src/serialization.jl b/src/serialization.jl
index 9fb2b25e7..47cfade90 100644
--- a/src/serialization.jl
+++ b/src/serialization.jl
@@ -1,10 +1,13 @@
 """
-    restoreoptsum!(m::LinearMixedModel, io::IO)
-    restoreoptsum!(m::LinearMixedModel, filename)
+    restoreoptsum!(m::LinearMixedModel, io::IO; atol::Real=0, rtol::Real=atol>0 ? 0 : √eps)
+    restoreoptsum!(m::LinearMixedModel, filename; atol::Real=0, rtol::Real=atol>0 ? 0 : √eps)
 
 Read, check, and restore the `optsum` field from a JSON stream or filename.
 """
-function restoreoptsum!(m::LinearMixedModel{T}, io::IO) where {T}
+function restoreoptsum!(
+    m::LinearMixedModel{T}, io::IO; atol::Real=zero(T),
+    rtol::Real=atol > 0 ? zero(T) : √eps(T),
+) where {T}
     dict = JSON3.read(io)
     ops = m.optsum
     okay =
@@ -22,7 +25,9 @@ function restoreoptsum!(m::LinearMixedModel{T}, io::IO) where {T}
     copyto!(ops.initial, dict.initial)
     copyto!(ops.final, dict.final)
     for (v, f) in (:initial => :finitial, :final => :fmin)
-        if !isapprox(objective(updateL!(setθ!(m, getfield(ops, v)))), getfield(ops, f))
+        if !isapprox(
+            objective(updateL!(setθ!(m, getfield(ops, v)))), getfield(ops, f); rtol, atol
+        )
             throw(ArgumentError("model m at $v does not give stored $f"))
         end
     end
@@ -40,9 +45,9 @@ function restoreoptsum!(m::LinearMixedModel{T}, io::IO) where {T}
     return m
 end
 
-function restoreoptsum!(m::LinearMixedModel{T}, filename) where {T}
+function restoreoptsum!(m::LinearMixedModel{T}, filename; kwargs...) where {T}
     open(filename, "r") do io
-        restoreoptsum!(m, io)
+        restoreoptsum!(m, io; kwargs...)
     end
 end
 
diff --git a/test/pls.jl b/test/pls.jl
index 1e1edc2e4..2b2ab3928 100644
--- a/test/pls.jl
+++ b/test/pls.jl
@@ -167,7 +167,7 @@ end
         sigma0row = only(filter(r -> r.p == :σ && iszero(r.ζ), dspr01.tbl))
         @test sigma0row.σ ≈ dspr01.m.σ
         @test sigma0row.β1 ≈ only(dspr01.m.β)
-        @test sigma0row.θ1 ≈ only(dspr01.m.θ)        
+        @test sigma0row.θ1 ≈ only(dspr01.m.θ)
     end
 end
 
@@ -500,14 +500,23 @@ end
 
     @testset "optsumJSON" begin
         fm = last(models(:sleepstudy))
-            # using a IOBuffer for saving JSON
+        # using a IOBuffer for saving JSON
         saveoptsum(seekstart(io), fm)
         m = LinearMixedModel(fm.formula, MixedModels.dataset(:sleepstudy))
         restoreoptsum!(m, seekstart(io))
         @test loglikelihood(fm) ≈ loglikelihood(m)
         @test bic(fm) ≈ bic(m)
         @test coef(fm) ≈ coef(m)
-            # using a temporary file for saving JSON
+
+        fm_mod = deepcopy(fm)
+        fm_mod.optsum.fmin += 1
+        saveoptsum(seekstart(io), fm_mod)
+        @test_throws(ArgumentError("model m at final does not give stored fmin"),
+                     restoreoptsum!(m, seekstart(io)))
+        restoreoptsum!(m, seekstart(io); atol=1)
+        @test m.optsum.fmin - fm.optsum.fmin ≈ 1
+
+        # using a temporary file for saving JSON
         fnm = first(mktemp())
         saveoptsum(fnm, fm)
         m = LinearMixedModel(fm.formula, MixedModels.dataset(:sleepstudy))
@@ -516,6 +525,7 @@ end
         @test bic(fm) ≈ bic(m)
         @test coef(fm) ≈ coef(m)
     end
+
     @testset "profile" begin
         pr = profile(last(models(:sleepstudy)))
         tbl = pr.tbl