diff --git a/Project.toml b/Project.toml
index 9e02e7138..e2f8252e9 100644
--- a/Project.toml
+++ b/Project.toml
@@ -33,6 +33,7 @@ TypedTables = "9d95f2ec-7b3d-5a63-8d20-e2491e220bb9"
 Aqua = "0.8"
 Arrow = "1, 2"
 BSplineKit = "0.14, 0.15, 0.16, 0.17"
+Chairmarks = "1"
 DataAPI = "1"
 DataFrames = "1"
 Distributions = "0.21, 0.22, 0.23, 0.24, 0.25"
@@ -50,6 +51,7 @@ ProgressMeter = "1.7"
 Random = "1"
 SparseArrays = "1"
 StableRNGs = "0.1, 1"
+StandardizedPredictors = "1"
 StaticArrays = "0.11, 0.12, 1"
 Statistics = "1"
 StatsAPI = "1.5"
@@ -65,12 +67,23 @@ julia = "1.8"
 
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
+Chairmarks = "0ca39b1e-fe0b-4e98-acfc-b1656634c4de"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
+StandardizedPredictors = "5064a6a7-f8c2-40e2-8bdc-797ec6f1ae18"
 Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Aqua", "DataFrames", "ExplicitImports", "InteractiveUtils", "StableRNGs", "Suppressor", "Test"]
+test = [
+    "Aqua",
+    "Chairmarks",
+    "DataFrames",
+    "ExplicitImports",
+    "InteractiveUtils",
+    "StableRNGs",
+    "StandardizedPredictors",
+    "Suppressor",
+    "Test"]
diff --git a/bench/runbenchmarks.jl b/bench/runbenchmarks.jl
new file mode 100644
index 000000000..bb4b4b4ab
--- /dev/null
+++ b/bench/runbenchmarks.jl
@@ -0,0 +1,95 @@
+using Chairmarks, MixedModels, StandardizedPredictors
+using MixedModels: dataset, FormulaTerm, Table
+
+@isdefined(contrasts) || const contrasts = Dict{Symbol, Any}()
+
+contrasts[:F] = HelmertCoding()       # mrk17_exp1
+contrasts[:P] = HelmertCoding()       # mrk17_exp1
+contrasts[:Q] = HelmertCoding()       # mrk17_exp1
+contrasts[:lQ] = HelmertCoding()      # mrk17_exp1
+contrasts[:lT] = HelmertCoding()      # mrk17_exp1
+contrasts[:ch] = HelmertCoding()      # contra
+contrasts[:load] = HelmertCoding()    # kb07
+contrasts[:prec] = HelmertCoding()    # kb07
+contrasts[:service] = HelmertCoding() # insteval
+contrasts[:spkr] = HelmertCoding()    # kb07
+contrasts[:height] = Center()         # grouseticks
+contrasts[:gender] = HelmertCoding()  # verbagg
+contrasts[:btype] = EffectsCoding()   # verbagg
+contrasts[:situ] = HelmertCoding()    # verbagg
+contrasts[:mode] = HelmertCoding()    # verbagg
+
+tbl = Table(
+    dsnm = [
+        :dyestuff2, :dyestuff, :pastes, :pastes, :machines, :penicillin,
+        :sleepstudy, :sleepstudy, :sleepstudy, :sleepstudy, :kb07, :kb07, 
+        :mrk17_exp1, :kb07, :insteval, :insteval, :mrk17_exp1, :d3, :ml1m,
+    ],
+    secs = append!(
+        fill(0.1f0, 12),
+        fill(1.0f0, 2),
+        fill(5.0f0, 2),
+        fill(25.0f0, 4)
+    ),
+    frm = FormulaTerm[
+        @formula(yield ~ 1 + (1|batch)),
+        @formula(yield ~ 1 + (1|batch)),
+        @formula(strength ~ 1 + (1 | batch & cask)),
+        @formula(strength ~ 1 + (1 | batch / cask)),
+        @formula(score ~ 1 + (1 | Worker) + (1 | Machine)),
+        @formula(diameter ~ 1 + (1 | plate) + (1 | sample)),
+        @formula(reaction ~ 1 + days + (1 | subj)),
+        @formula(reaction ~ 1 + days + zerocorr(1 + days | subj)),
+        @formula(reaction ~ 1 + days + (1 | subj) + (0 + days | subj)),
+        @formula(reaction ~ 1 + days + (1 + days | subj)),
+        @formula(log(rt_trunc) ~ 1 + spkr + prec + load + (1 | subj) + (1 | item)),
+        @formula(log(rt_trunc) ~ 1 + spkr * prec * load + (1 | subj) + (1 + prec | item)),
+        @formula(1000 / rt ~ 1 + F * P * Q * lQ * lT + (1 | item) + (1 | subj)),
+        @formula(
+            log(rt_trunc) ~
+                1 + spkr * prec * load + (1 + spkr + prec + load | subj) +
+                (1 + spkr + prec + load | item)
+        ),
+        @formula(y ~ 1 + service * dept + (1 | s) + (1 | d)),
+        @formula(y ~ 1 + service + (1 | s) + (1 | d) + (1 | dept)),
+        @formula(
+            1000 / rt ~
+                1 + F * P * Q * lQ * lT + (1 + P + Q + lQ + lT | item) +
+                (1 + F + P + Q + lQ + lT | subj)
+        ),
+        @formula(y ~ 1 + u + (1 + u | g) + (1 + u | h) + (1 + u | i)),
+        @formula(y ~ 1 + (1 | g) + (1 | h)),
+    ]
+)
+
+# linmark(f, d, t) = @b fit(MixedModel, f, dataset(d); contrasts, progress=false) seconds=t
+
+@track (@b (first(tbl.frm), dataset(first(tbl.dsnm))) fit(MixedModel, first(_), last(_); progress=false)).time
+@track (@b (tbl.frm[2], dataset(tbl.dsnm[2])) fit(MixedModel, first(_), last(_); progress=false)).time
+@track (@b (tbl.frm[3], dataset(tbl.dsnm[3])) fit(MixedModel, first(_), last(_); progress=false)).time
+@track (@b (tbl.frm[4], dataset(tbl.dsnm[4])) fit(MixedModel, first(_), last(_); progress=false)).time
+@track (@b (tbl.frm[5], dataset(tbl.dsnm[5])) fit(MixedModel, first(_), last(_); progress=false)).time
+@track (@b (tbl.frm[6], dataset(tbl.dsnm[6])) fit(MixedModel, first(_), last(_); progress=false)).time
+@track (@b (tbl.frm[7], dataset(tbl.dsnm[7])) fit(MixedModel, first(_), last(_); progress=false)).time
+
+# function runbmrk(tbl)
+#     return Table([(; bmk=linmark(f, d, t), dsnm = d, frm=f) for (d, t, f) in tbl])
+# end
+
+# gltbl = Table(
+#     dsnm = [:contra, :contra, :verbagg, :grouseticks],
+#     secs = [2.0, 2.0, 15.0, 15.0],
+#     dist = [Bernoulli(), Bernoulli(), Bernoulli(), Poisson()],
+#     frm = FormulaTerm[
+#         @formula(use ~ 1+age+abs2(age)+urban+livch+(1|urban&dist)),
+#         @formula(use ~ 1+age+abs2(age)+urban+(≠("0"))(livch)+(1+urban|dist)),
+#         @formula(r2 ~ 1+anger+gender+btype+situ+(1|subj)+(1|item)),
+#         @formula(ticks ~ 1+year+height+(1|index)+(1|brood)+(1|location)),
+#     ]
+# )
+
+# glmark(f, d, r, t; init_from_lmm=()) = @b fit(MixedModel, f, dataset(d), r; init_from_lmm, contrasts, progress=false) seconds=t
+
+# function runglbmk(tbl; init_from_lmm=())
+#     return Table((; bmk=glmark(f, d, r, t; init_from_lmm), dsnm=d, dist=r, frm=f) for (d, t, r, f) in tbl)
+# end