From 887d2f2be911a739357a0a9add84c4eb46e5704f Mon Sep 17 00:00:00 2001
From: Phillip Alday <me@phillipalday.com>
Date: Fri, 22 Aug 2025 13:31:50 -0500
Subject: [PATCH 1/2] remove dataset and datasets

The dataset functionality is now housed in MixedModelsDatasets. By fully separating out this functionality, we set ourselves up for easier and lighter weight upgrades to which and how we store datasets.
---
 benchmark/benchmarks.jl     |  2 +-
 docs/Project.toml           |  1 +
 docs/src/GaussHermite.md    |  4 ++--
 docs/src/bootstrap.md       |  6 +++---
 docs/src/constructors.md    | 24 ++++++++++++------------
 docs/src/derivatives.md     |  6 +++---
 docs/src/index.md           |  6 +++---
 docs/src/mime.md            |  8 ++++----
 docs/src/optimization.md    | 10 +++++-----
 docs/src/prediction.md      |  5 +++--
 docs/src/rankdeficiency.md  |  3 ++-
 src/MixedModels.jl          |  4 ++--
 test/FactorReTerm.jl        |  3 ++-
 test/bootstrap.jl           |  3 ++-
 test/likelihoodratiotest.jl |  3 ++-
 test/mime.jl                |  3 ++-
 test/misc.jl                |  2 +-
 test/modelcache.jl          |  2 +-
 test/optsummary.jl          |  2 +-
 test/pirls.jl               |  2 +-
 test/predict.jl             |  2 +-
 test/utilities.jl           | 10 ----------
 22 files changed, 54 insertions(+), 57 deletions(-)

diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl
index 4c05d5c0d..27601c62f 100644
--- a/benchmark/benchmarks.jl
+++ b/benchmark/benchmarks.jl
@@ -1,5 +1,5 @@
 using BenchmarkTools, MixedModels
-using MixedModels: dataset
+using MixedModelsDatasets: dataset
 
 const SUITE = BenchmarkGroup()
 
diff --git a/docs/Project.toml b/docs/Project.toml
index 6c7cfad17..17ebe3b0c 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -9,6 +9,7 @@ ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 FreqTables = "da1fdf0e-e0ff-5433-a45f-9bb5ff651cb1"
 Gadfly = "c91e804a-d5a3-530f-b6f0-dfbca275c004"
 MixedModels = "ff71e718-51f3-5ec2-a782-8ffcbfa3c316"
+MixedModelsDatasets = "7e9fb7ac-9f67-43bf-b2c8-96ba0796cbb6"
 ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 StatsAPI = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
diff --git a/docs/src/GaussHermite.md b/docs/src/GaussHermite.md
index 966b742ba..e04ebad4b 100644
--- a/docs/src/GaussHermite.md
+++ b/docs/src/GaussHermite.md
@@ -73,7 +73,7 @@ The definition of `MixedModels.GHnorm` is similar to the `gausshermitenorm` func
 GHnorm
 ```
 ```@example Main
-using MixedModels
+using MixedModels, MixedModelsDatasets
 GHnorm(3)
 ```
 
@@ -100,7 +100,7 @@ Several covariates were recorded including the woman's age (centered at the mean
 The version of the data used here is that used in review of multilevel modeling software conducted by the Center for Multilevel Modelling, currently at University of Bristol (http://www.bristol.ac.uk/cmm/learning/mmsoftware/data-rev.html).
 These data are available as the `:contra` dataset.
 ```@example Main
-contra = DataFrame(MixedModels.dataset(:contra))
+contra = DataFrame(MixedModelsDatasets.dataset(:contra))
 describe(contra)
 ```
 
diff --git a/docs/src/bootstrap.md b/docs/src/bootstrap.md
index dbe51e620..9fd5d1ec0 100644
--- a/docs/src/bootstrap.md
+++ b/docs/src/bootstrap.md
@@ -27,12 +27,12 @@ package for [`R`](https://www.r-project.org) is fit by
 ```@example Main
 using DataFrames
 using Gadfly          # plotting package
-using MixedModels
+using MixedModels, MixedModelsDatasets
 using Random
 ```
 
 ```@example Main
-dyestuff = MixedModels.dataset(:dyestuff)
+dyestuff = MixedModelsDatasets.dataset(:dyestuff)
 m1 = fit(MixedModel, @formula(yield ~ 1 + (1 | batch)), dyestuff)
 ```
 
@@ -88,7 +88,7 @@ However, it is not as straightforward to detect singularity in vector-valued ran
 
 For example, if we bootstrap a model fit to the `sleepstudy` data
 ```@example Main
-sleepstudy = MixedModels.dataset(:sleepstudy)
+sleepstudy = MixedModelsDatasets.dataset(:sleepstudy)
 contrasts = Dict(:subj => Grouping())
 m2 = let f = @formula reaction ~ 1+days+(1+days|subj)
     fit(MixedModel, f, sleepstudy; contrasts)
diff --git a/docs/src/constructors.md b/docs/src/constructors.md
index b6cf4a642..3df61d6c8 100644
--- a/docs/src/constructors.md
+++ b/docs/src/constructors.md
@@ -14,8 +14,8 @@ using DisplayAs
 ```
 
 ```@example Main
-using DataFrames, MixedModels, StatsModels
-dyestuff = MixedModels.dataset(:dyestuff)
+using DataFrames, MixedModels, MixedModelsDatasets, StatsModels
+dyestuff = MixedModelsDatasets.dataset(:dyestuff)
 ```
 
 ```@example Main
@@ -71,7 +71,7 @@ Notice that both are equivalent.
 
 ```@example Main
 using BenchmarkTools
-dyestuff2 = MixedModels.dataset(:dyestuff2)
+dyestuff2 = MixedModelsDatasets.dataset(:dyestuff2)
 @benchmark fit(MixedModel, $fm, $dyestuff2)
 ```
 
@@ -110,7 +110,7 @@ It corresponds to a shift in the intercept for each level of the grouping factor
 The *sleepstudy* data are observations of reaction time, `reaction`, on several subjects, `subj`, after 0 to 9 days of sleep deprivation, `days`.
 A model with random intercepts and random slopes for each subject, allowing for within-subject correlation of the slope and intercept, is fit as
 ```@example Main
-sleepstudy = MixedModels.dataset(:sleepstudy)
+sleepstudy = MixedModelsDatasets.dataset(:sleepstudy)
 fm2 = fit(MixedModel, @formula(reaction ~ 1 + days + (1 + days|subj)), sleepstudy)
 DisplayAs.Text(ans) # hide
 ```
@@ -120,14 +120,14 @@ DisplayAs.Text(ans) # hide
 A model for the *Penicillin* data incorporates random effects for the plate, and for the sample.
 As every sample is used on every plate these two factors are *crossed*.
 ```@example Main
-penicillin = MixedModels.dataset(:penicillin)
+penicillin = MixedModelsDatasets.dataset(:penicillin)
 fm3 = fit(MixedModel, @formula(diameter ~ 1 + (1|plate) + (1|sample)), penicillin)
 DisplayAs.Text(ans) # hide
 ```
 
 In contrast, the `cask` grouping factor is *nested* within the `batch` grouping factor in the *Pastes* data.
 ```@example Main
-pastes = DataFrame(MixedModels.dataset(:pastes))
+pastes = DataFrame(MixedModelsDatasets.dataset(:pastes))
 describe(pastes)
 ```
 This can be expressed using the solidus (the "`/`" character) to separate grouping factors, read "`cask` nested within `batch`":
@@ -147,7 +147,7 @@ In observational studies it is common to encounter *partially crossed* grouping
 For example, the *InstEval* data are course evaluations by students, `s`, of instructors, `d`.
 Additional covariates include the academic department, `dept`, in which the course was given and `service`, whether or not it was a service course.
 ```@example Main
-insteval = MixedModels.dataset(:insteval)
+insteval = MixedModelsDatasets.dataset(:insteval)
 fm5 = fit(MixedModel, @formula(y ~ 1 + service * dept + (1|s) + (1|d)), insteval)
 DisplayAs.Text(ans) # hide
 ```
@@ -222,7 +222,7 @@ To create a GLMM representation, the distribution family for the response, and p
 You can either use `fit(MixedModel, ...)` or `glmm(...)` to fit the model. For instance:
 
 ```@example Main
-verbagg = MixedModels.dataset(:verbagg)
+verbagg = MixedModelsDatasets.dataset(:verbagg)
 verbaggform = @formula(r2 ~ 1 + anger + gender + btype + situ + mode + (1|subj) + (1|item));
 gm1 = fit(MixedModel, verbaggform, verbagg, Bernoulli())
 DisplayAs.Text(ans) # hide
@@ -258,7 +258,7 @@ The optional argument `nAGQ=k` causes evaluation of the deviance function to use
 adaptive Gauss-Hermite quadrature rule.
 This method only applies to models with a single, simple, scalar random-effects term, such as
 ```@example Main
-contraception = MixedModels.dataset(:contra)
+contraception = MixedModelsDatasets.dataset(:contra)
 contraform = @formula(use ~ 1 + age + abs2(age) + livch + urban + (1|dist));
 bernoulli = Bernoulli()
 deviances = Dict{Symbol,Float64}()
@@ -511,7 +511,7 @@ sum(leverage(fm2))
 
 When a model converges to a singular covariance, such as
 ```@example Main
-fm3 = fit(MixedModel, @formula(yield ~ 1+(1|batch)), MixedModels.dataset(:dyestuff2))
+fm3 = fit(MixedModel, @formula(yield ~ 1+(1|batch)), MixedModelsDatasets.dataset(:dyestuff2))
 DisplayAs.Text(ans) # hide
 ```
 the effective degrees of freedom is the lower bound.
@@ -522,7 +522,7 @@ sum(leverage(fm3))
 Models for which the estimates of the variances of the random effects are large relative to the residual variance have effective degrees of freedom close to the upper bound.
 ```@example Main
 fm4 = fit(MixedModel, @formula(diameter ~ 1+(1|plate)+(1|sample)),
-    MixedModels.dataset(:penicillin))
+    MixedModelsDatasets.dataset(:penicillin))
 DisplayAs.Text(ans) # hide
 ```
 ```@example Main
@@ -532,7 +532,7 @@ sum(leverage(fm4))
 Also, a model fit by the REML criterion generally has larger estimates of the variance components and hence a larger effective degrees of freedom.
 ```@example Main
 fm4r = fit(MixedModel, @formula(diameter ~ 1+(1|plate)+(1|sample)),
-    MixedModels.dataset(:penicillin), REML=true)
+    MixedModelsDatasets.dataset(:penicillin), REML=true)
 DisplayAs.Text(ans) # hide
 ```
 ```@example Main
diff --git a/docs/src/derivatives.md b/docs/src/derivatives.md
index 52a290e6e..c4ed7d305 100644
--- a/docs/src/derivatives.md
+++ b/docs/src/derivatives.md
@@ -15,9 +15,9 @@ ForwardDiff.hessian(::LinearMixedModel{T}, ::Vector{T}) where {T}
 ### Exact zero at optimum for trivial models
 
 ```@example Derivatives
-using MixedModels, ForwardDiff
+using MixedModels, MixedModelsDatasets, ForwardDiff
 using DisplayAs # hide
-fm1 = lmm(@formula(yield ~ 1 + (1|batch)), MixedModels.dataset(:dyestuff2))
+fm1 = lmm(@formula(yield ~ 1 + (1|batch)), MixedModelsDatasets.dataset(:dyestuff2))
 DisplayAs.Text(ans) # hide
 ```
 
@@ -32,7 +32,7 @@ ForwardDiff.hessian(fm1)
 ### Approximate zero at optimum for non trivial models
 
 ```@example Derivatives
-fm2 = lmm(@formula(reaction ~ 1 + days + (1+days|subj)), MixedModels.dataset(:sleepstudy))
+fm2 = lmm(@formula(reaction ~ 1 + days + (1+days|subj)), MixedModelsDatasets.dataset(:sleepstudy))
 DisplayAs.Text(ans) # hide
 ```
 
diff --git a/docs/src/index.md b/docs/src/index.md
index fd6eeddb8..ae4f9d6f1 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -16,8 +16,8 @@ You can fit a model using a `lmer`-style model formula using `@formula` and a da
 Here is a short example of how to fit a linear mixed-effects modeling using the `dyestuff` dataset:
 
 ```@example Main
-using DataFrames, MixedModels           # load packages
-dyestuff = MixedModels.dataset(:dyestuff);              # load dataset
+using DataFrames, MixedModels, MixedModelsDatasets  # load packages
+dyestuff = MixedModelsDatasets.dataset(:dyestuff);   # load dataset
 
 lmod = lmm(@formula(yield ~ 1 + (1|batch)), dyestuff)   # fit the model!
 DisplayAs.Text(ans) # hide
@@ -28,7 +28,7 @@ A quick example of generalized linear model using the `verbagg` dataset:
 
 ```@example Main
 using DataFrames, MixedModels               # load packages
-verbagg = MixedModels.dataset(:verbagg);    # load dataset
+verbagg = MixedModelsDatasets.dataset(:verbagg);    # load dataset
 
 frm = @formula(r2 ~ 1 + anger + gender + btype + situ + mode + (1|subj) + (1|item));
 bernmod = glmm(frm, verbagg, Bernoulli())   # fit the model!
diff --git a/docs/src/mime.md b/docs/src/mime.md
index 9dd1d16c8..7a25f42af 100644
--- a/docs/src/mime.md
+++ b/docs/src/mime.md
@@ -12,7 +12,7 @@ Packages like `IJulia` and `Documenter` can often detect the presence of these d
 
 
 ```@example Main
-using MixedModels
+using MixedModels, MixedModelsDatasets
 form = @formula(rt_trunc ~ 1 + spkr * prec * load +
                           (1 + load | item) +
                           (1 + spkr + prec + load | subj))
@@ -21,7 +21,7 @@ contr = Dict(:spkr => EffectsCoding(),
              :load => EffectsCoding(),
              :item => Grouping(),
              :subj => Grouping())
-kbm = fit(MixedModel, form, MixedModels.dataset(:kb07); contrasts=contr)
+kbm = fit(MixedModel, form, MixedModelsDatasets.dataset(:kb07); contrasts=contr)
 ```
 
 Note that the display here is more succinct than the standard REPL display:
@@ -52,8 +52,8 @@ kbm.optsum
 ```
 
 ```@example Main
-m0 = fit(MixedModel, @formula(reaction ~ 1 + (1|subj)), MixedModels.dataset(:sleepstudy))
-m1 = fit(MixedModel, @formula(reaction ~ 1 + days + (1+days|subj)), MixedModels.dataset(:sleepstudy))
+m0 = fit(MixedModel, @formula(reaction ~ 1 + (1|subj)), MixedModelsDatasets.dataset(:sleepstudy))
+m1 = fit(MixedModel, @formula(reaction ~ 1 + days + (1+days|subj)), MixedModelsDatasets.dataset(:sleepstudy))
 MixedModels.likelihoodratiotest(m0,m1)
 ```
 
diff --git a/docs/src/optimization.md b/docs/src/optimization.md
index a0aadb559..8a8541cf4 100644
--- a/docs/src/optimization.md
+++ b/docs/src/optimization.md
@@ -101,10 +101,10 @@ For the simple example
 using DisplayAs
 ```
 ```@example Main
-using BenchmarkTools, DataFrames, MixedModels
+using BenchmarkTools, DataFrames, MixedModels, MixedModelsDatasets
 ```
 ```@example Main
-dyestuff = MixedModels.dataset(:dyestuff)
+dyestuff = MixedModelsDatasets.dataset(:dyestuff)
 fm1 = fit(MixedModel, @formula(yield ~ 1 + (1|batch)), dyestuff)
 DisplayAs.Text(ans) # hide
 ```
@@ -133,7 +133,7 @@ Furthermore, there is only one block in $\Lambda_\theta$.
 
 For a vector-valued random-effects term, as in
 ```@example Main
-sleepstudy = MixedModels.dataset(:sleepstudy)
+sleepstudy = MixedModelsDatasets.dataset(:sleepstudy)
 fm2 = fit(MixedModel, @formula(reaction ~ 1+days+(1+days|subj)), sleepstudy)
 DisplayAs.Text(ans) # hide
 ```
@@ -168,7 +168,7 @@ MixedModels.getθ(t31)
 Random-effects terms with distinct grouping factors generate distinct elements of the `reterms` field of the `LinearMixedModel` object.
 Multiple `ReMat` objects are sorted by decreasing numbers of random effects.
 ```@example Main
-penicillin = MixedModels.dataset(:penicillin)
+penicillin = MixedModelsDatasets.dataset(:penicillin)
 fm4 = fit(MixedModel,
     @formula(diameter ~ 1 + (1|sample) + (1|plate)),
     penicillin)
@@ -314,7 +314,7 @@ Poisson
 A `GeneralizedLinearMixedModel` object is generated from a formula, data frame and distribution family.
 
 ```@example Main
-verbagg = MixedModels.dataset(:verbagg)
+verbagg = MixedModelsDatasets.dataset(:verbagg)
 const vaform = @formula(r2 ~ 1 + anger + gender + btype + situ + (1|subj) + (1|item));
 mdl = GeneralizedLinearMixedModel(vaform, verbagg, Bernoulli());
 typeof(mdl)
diff --git a/docs/src/prediction.md b/docs/src/prediction.md
index 6980b2d31..751193488 100644
--- a/docs/src/prediction.md
+++ b/docs/src/prediction.md
@@ -5,10 +5,11 @@ We recommend the [MixedModelsSim.jl](https://github.com/RePsychLing/MixedModelsS
 ```@example Main
 using DataFrames
 using MixedModels
+using MixedModelsDatasets
 using StatsBase
 using DisplayAs # hide
 # use a DataFrame to make it easier to change things later
-slp = DataFrame(MixedModels.dataset(:sleepstudy))
+slp = DataFrame(MixedModelsDatasets.dataset(:sleepstudy))
 slpm = fit(MixedModel, @formula(reaction ~ 1 + days + (1|subj)), slp)
 DisplayAs.Text(slpm) # hide
 ```
@@ -80,7 +81,7 @@ predict(slpm, slp2; new_re_levels=:population)
 For generalized linear mixed models, there is an additional keyword argument to `predict`: `type` specifies whether the predictions are returned on the scale of the linear predictor (`:linpred`) or on the level of the response `(:response)` (i.e. the level at which the values were originally observed).
 
 ```@example Main
-cbpp = DataFrame(MixedModels.dataset(:cbpp))
+cbpp = DataFrame(MixedModelsDatasets.dataset(:cbpp))
 cbpp.rate = cbpp.incid ./ cbpp.hsz
 gm = fit(MixedModel, @formula(rate ~ 1 + period + (1|herd)), cbpp, Binomial(), wts=float(cbpp.hsz))
 predict(gm, cbpp; type=:response) ≈ fitted(gm)
diff --git a/docs/src/rankdeficiency.md b/docs/src/rankdeficiency.md
index 44b874f2a..9960c975a 100644
--- a/docs/src/rankdeficiency.md
+++ b/docs/src/rankdeficiency.md
@@ -2,6 +2,7 @@
 
 ```@setup Main
 using MixedModels
+using MixedModelsDatasets
 using DisplayAs
 ```
 
@@ -85,7 +86,7 @@ In addition to handling naturally occurring rank deficiency in the random effect
 For example, we can use `fulldummy` to fit both an intercept term and $n$ indicator variables in the random effects for a categorical variable with $n$ levels instead of the usual $n-1$ contrasts.
 
 ```@example Main
-kb07 = MixedModels.dataset(:kb07)
+kb07 = MixedModelsDatasets.dataset(:kb07)
 contrasts = Dict(var => HelmertCoding() for var in (:spkr, :prec, :load))
 fit(MixedModel, @formula(rt_raw ~ spkr * prec * load + (1|subj) + (1+prec|item)), kb07; contrasts=contrasts)
 DisplayAs.Text(ans) # hide
diff --git a/src/MixedModels.jl b/src/MixedModels.jl
index 852be8fdd..111c1e9ba 100644
--- a/src/MixedModels.jl
+++ b/src/MixedModels.jl
@@ -19,7 +19,7 @@ using LinearAlgebra: UpperTriangular, cond, diag, diagind, dot, eigen, isdiag
 using LinearAlgebra: ldiv!, lmul!, logdet, mul!, norm, normalize, normalize!, qr
 using LinearAlgebra: rank, rdiv!, rmul!, svd, tril!
 using Markdown: Markdown
-using MixedModelsDatasets: dataset, datasets
+using MixedModelsDatasets: dataset
 using PooledArrays: PooledArrays, PooledArray
 using NLopt: NLopt
 using PrecompileTools: PrecompileTools, @setup_workload, @compile_workload
@@ -167,7 +167,7 @@ export @formula,
 # TODO: move this to the correct spot in list once we've decided on name
 export savereplicates, restorereplicates
 
-@compat public rePCA, PCA, dataset, datasets
+@compat public rePCA, PCA, opt_params, optimizers
 
 """
     MixedModel
diff --git a/test/FactorReTerm.jl b/test/FactorReTerm.jl
index 8f1993795..516042e3f 100644
--- a/test/FactorReTerm.jl
+++ b/test/FactorReTerm.jl
@@ -6,7 +6,8 @@ using SparseArrays
 using StatsModels
 using Test
 
-using MixedModels: dataset, levels, modelcols, nlevs
+using MixedModels: levels, modelcols, nlevs
+using MixedModelsDatasets: dataset
 
 const LMM = LinearMixedModel
 
diff --git a/test/bootstrap.jl b/test/bootstrap.jl
index 16764e904..777a7c74c 100644
--- a/test/bootstrap.jl
+++ b/test/bootstrap.jl
@@ -9,7 +9,8 @@ using Suppressor
 using Tables
 using Test
 
-using MixedModels: dataset, MixedModelBootstrap
+using MixedModels: MixedModelBootstrap
+using MixedModelsDatasets: dataset
 
 include("modelcache.jl")
 
diff --git a/test/likelihoodratiotest.jl b/test/likelihoodratiotest.jl
index 86b869f77..d24f4d378 100644
--- a/test/likelihoodratiotest.jl
+++ b/test/likelihoodratiotest.jl
@@ -4,7 +4,8 @@ using MixedModels
 using Suppressor
 using Test
 
-using MixedModels: dataset, likelihoodratiotest
+using MixedModels: likelihoodratiotest
+using MixedModelsDatasets: dataset
 using GLM: ProbitLink
 using StatsModels: lrtest, isnested
 
diff --git a/test/mime.jl b/test/mime.jl
index 3d6b33a02..0f380a340 100644
--- a/test/mime.jl
+++ b/test/mime.jl
@@ -2,8 +2,9 @@ using MixedModels
 using Suppressor
 using Test
 
-using MixedModels: dataset, likelihoodratiotest
+using MixedModels: likelihoodratiotest
 using MixedModels: pirls!, setβθ!, setθ!, updateL!
+using MixedModelsDatasets: dataset
 
 include("modelcache.jl")
 
diff --git a/test/misc.jl b/test/misc.jl
index dc2d3b740..5d973514e 100644
--- a/test/misc.jl
+++ b/test/misc.jl
@@ -1,7 +1,7 @@
 using MixedModels
 using Test
 
-using MixedModels: dataset
+using MixedModelsDatasets: dataset
 
 @testset "formula misspecification" begin
     dyestuff = dataset(:dyestuff)
diff --git a/test/modelcache.jl b/test/modelcache.jl
index 259d62bdd..c45a7dff0 100644
--- a/test/modelcache.jl
+++ b/test/modelcache.jl
@@ -1,5 +1,5 @@
 using MixedModels
-using MixedModels: dataset
+using MixedModelsDatasets: dataset
 
 @isdefined(gfms) || const global gfms = Dict(
     :cbpp => [@formula((incid / hsz) ~ 1 + period + (1 | herd))],
diff --git a/test/optsummary.jl b/test/optsummary.jl
index ad6b1a9a0..45b300754 100644
--- a/test/optsummary.jl
+++ b/test/optsummary.jl
@@ -1,5 +1,5 @@
 using MixedModels
-using MixedModels: dataset
+using MixedModelsDatasets: dataset
 using Test
 
 include("modelcache.jl")
diff --git a/test/pirls.jl b/test/pirls.jl
index 5e9ff53d2..aeb16f834 100644
--- a/test/pirls.jl
+++ b/test/pirls.jl
@@ -7,7 +7,7 @@ using Tables
 using Test
 
 using GLM: Link
-using MixedModels: dataset
+using MixedModelsDatasets: dataset
 
 include("modelcache.jl")
 
diff --git a/test/predict.jl b/test/predict.jl
index c34c58529..b3d35ab82 100644
--- a/test/predict.jl
+++ b/test/predict.jl
@@ -8,7 +8,7 @@ using Tables
 using Test
 
 using GLM: Link, linkfun, linkinv
-using MixedModels: dataset
+using MixedModelsDatasets: dataset
 
 include("modelcache.jl")
 
diff --git a/test/utilities.jl b/test/utilities.jl
index 1f5a8458a..3d92ef888 100644
--- a/test/utilities.jl
+++ b/test/utilities.jl
@@ -44,16 +44,6 @@ end
     @test_logs (:warn, r"hide_progress") replicate(string, 1; hide_progress=true)
 end
 
-@testset "datasets" begin
-    @test isa(MixedModels.datasets(), Vector{String})
-    @test length(MixedModels.dataset(:dyestuff)) == 2
-    @test length(MixedModels.dataset("dyestuff")) == 2
-    dyestuff = MixedModels.dataset(:dyestuff)
-    @test keys(dyestuff) == [:batch, :yield]
-    @test length(dyestuff.batch) == 30
-    @test_throws ArgumentError MixedModels.dataset(:foo)
-end
-
 @testset "PCA" begin
     io = IOBuffer()
     pca = models(:kb07)[3].PCA.item

From 0400525000706935729271efb05f4415430ec7ff Mon Sep 17 00:00:00 2001
From: Phillip Alday <me@phillipalday.com>
Date: Fri, 22 Aug 2025 13:32:28 -0500
Subject: [PATCH 2/2] NEWS

---
 NEWS.md | 1 +
 1 file changed, 1 insertion(+)

diff --git a/NEWS.md b/NEWS.md
index 790eabd58..a16bbc05f 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -8,6 +8,7 @@ MixedModels v5.0.0 Release Notes
 - Internal code around the default optimizer has been restructured. In particular, the NLopt backend has been moved to a submodule, which will make it easier to move it to an extension if we promote another backend to the default. [#853]
 - Internal code around optimization in profiling has been restructuring so that fitting done during calls to `profile` respect the `backend` and `optimizer` settings. [#853]
 - The `prfit!` convenience function has been removed. [#853]
+- The `dataset` and `datasets` functions have been removed. They are now housed in `MixedModelsDatasets`.[#854]
 
 MixedModels v4.38.0 Release Notes
 ==============================