diff --git a/.github/workflows/documenter.yml b/.github/workflows/documenter.yml
index 483fc4b8f..972c81db3 100644
--- a/.github/workflows/documenter.yml
+++ b/.github/workflows/documenter.yml
@@ -24,7 +24,7 @@ jobs:
       - uses: actions/checkout@v4
       - uses: julia-actions/setup-julia@v2
         with:
-          version: 1.8
+          version: 1.9
       - uses: julia-actions/cache@v2
         with:
             cache-compiled: "true"
diff --git a/.github/workflows/minimum.yml b/.github/workflows/minimum.yml
index f7e37a65f..4731ce9db 100644
--- a/.github/workflows/minimum.yml
+++ b/.github/workflows/minimum.yml
@@ -20,7 +20,7 @@ jobs:
     runs-on: ${{ matrix.os }}
     strategy:
       matrix:
-        julia-version: [1.8]
+        julia-version: [1.9]
         julia-arch: [x64]
         os: [ubuntu-22.04, macos-11, windows-2019]
     steps:
diff --git a/Project.toml b/Project.toml
index 225f9aa54..0f8f670a9 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.24.0"
+version = "4.25.0"
 
 [deps]
 Arrow = "69666777-d1a9-59fb-9406-91d4454c9d45"
@@ -28,6 +28,14 @@ StatsModels = "3eaba693-59b7-5ba5-a881-562e759f1c8d"
 StructTypes = "856f2bd8-1eba-4b0a-8007-ebc267875bd4"
 Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
 TypedTables = "9d95f2ec-7b3d-5a63-8d20-e2491e220bb9"
+Vcov = "ec2bfdc2-55df-4fc9-b9ae-4958c2cf2486"
+
+[weakdeps]
+GroupedArrays = "6407cd72-fade-4a84-8a1e-56e431fc1533"
+Vcov = "ec2bfdc2-55df-4fc9-b9ae-4958c2cf2486"
+
+[extensions]
+MixedModelsVcovExt = ["Vcov", "GroupedArrays"]
 
 [compat]
 Aqua = "0.8"
@@ -38,6 +46,7 @@ DataFrames = "1"
 Distributions = "0.21, 0.22, 0.23, 0.24, 0.25"
 ExplicitImports = "1.3"
 GLM = "1.8.2"
+GroupedArrays = "0.3"
 InteractiveUtils = "1"
 JSON3 = "1"
 LinearAlgebra = "1"
@@ -61,7 +70,8 @@ Suppressor = "0.2"
 Tables = "1"
 Test = "1"
 TypedTables = "1"
-julia = "1.8"
+Vcov = "0.7, 0.8"
+julia = "1.9"
 
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
@@ -71,6 +81,7 @@ InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+Vcov = "ec2bfdc2-55df-4fc9-b9ae-4958c2cf2486"
 
 [targets]
-test = ["Aqua", "DataFrames", "ExplicitImports", "InteractiveUtils", "StableRNGs", "Suppressor", "Test"]
+test = ["Aqua", "DataFrames", "ExplicitImports", "InteractiveUtils", "StableRNGs", "Suppressor", "Test", "Vcov"]
diff --git a/ext/MixedModelsVcovExt.jl b/ext/MixedModelsVcovExt.jl
new file mode 100644
index 000000000..a9a99abda
--- /dev/null
+++ b/ext/MixedModelsVcovExt.jl
@@ -0,0 +1,30 @@
+module MixedModelsVcovExt
+
+using MixedModels
+using MixedModels: StatsAPI
+using MixedModels.StatsModels: hasintercept
+
+using Vcov: Vcov, CovarianceEstimator, VcovData
+
+function Vcov.VcovData(m::MixedModel)
+    cm = crossmodelmatrix(m)
+    mm = modelmatrix(m)
+    # from just the fixed effects
+    ddf = nobs(m) - rank(m) - hasintercept(formula(m)) # dof_residual(m) # nobs(m) - round(Int, sum(leverage(m)))
+    resid = response(m) - mm * coef(m) # residuals(m)
+    @static if pkgversion(Vcov) >= v"0.8"
+        # Vcov 0.8 stores the inverse of crossmodelmatrix
+        return VcovData(mm, cm, inv(cm), resid, ddf)
+    else
+        return VcovData(mm, cm, resid, ddf)
+    end
+
+end
+
+# to avoid method ambiguities, we define this separately for each
+# relevant CovarianceEstimator
+StatsAPI.vcov(m::MixedModel, c::Vcov.SimpleCovariance) = vcov(VcovData(m), c)
+StatsAPI.vcov(m::MixedModel, c::Vcov.RobustCovariance) = vcov(VcovData(m), c)
+# StatsAPI.vcov(m::MixedModel, c::Vcov.ClusterCovariance) = vcov(VcovData(m), c)
+
+end # module
diff --git a/src/MixedModels.jl b/src/MixedModels.jl
index d03de81cd..a97a962bc 100644
--- a/src/MixedModels.jl
+++ b/src/MixedModels.jl
@@ -28,7 +28,7 @@ using SparseArrays: SparseArrays, SparseMatrixCSC, SparseVector, dropzeros!, nnz
 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: StatsAPI, aic, aicc, bic, coef, coefnames, coeftable, confint, crossmodelmatrix, deviance
 using StatsAPI: dof, dof_residual, fit, fit!, fitted, isfitted, islinear, leverage
 using StatsAPI: loglikelihood, meanresponse, modelmatrix, nobs, predict, r2, residuals
 using StatsAPI: response, responsename, stderror, vcov, weights
@@ -89,6 +89,7 @@ export @formula,
     condVar,
     condVartables,
     confint,
+    crossmodelmatrix,
     dataset,
     datasets,
     deviance,