diff --git a/Project.toml b/Project.toml
index bee758e..1f118ce 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "DynamicHMC"
 uuid = "bbc10e6e-7c05-544b-b16e-64fede858acb"
 authors = ["Tamas K. Papp <tkpapp@gmail.com>"]
-version = "3.4.1"
+version = "3.4.2"
 
 [deps]
 ArgCheck = "dce04be8-c92d-5529-be00-80e4d2c0e197"
@@ -16,12 +16,16 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 TensorCast = "02d47bb6-7ce6-556a-be16-bb1710789e2b"
 
+[extras]
+MCMCDiagnosticTools = "be115224-59cd-429b-ad48-344e309966f0"
+
 [compat]
 ArgCheck = "1, 2"
 DocStringExtensions = "0.8, 0.9"
 FillArrays = "0.13"
 LogDensityProblems = "1, 2"
 LogExpFunctions = "0.3"
+MCMCDiagnosticTools = "0.2"
 Parameters = "0.11, 0.12"
 ProgressMeter = "1"
 TensorCast = "0.4"
diff --git a/src/mcmc.jl b/src/mcmc.jl
index 621a146..0c71169 100644
--- a/src/mcmc.jl
+++ b/src/mcmc.jl
@@ -561,12 +561,15 @@ $(SIGNATURES)
 
 Given a vector of `results`, each containing a property `posterior_matrix` (eg obtained from
 [`mcmc_with_warmup`](@ref) with the same sample length), return a lazy view as an array
-indexed by `[draw_index, parameter_index, chain_index]`.
+indexed by `[draw_index, chain_index, parameter_index]`.
 
 This is useful as an input for eg `MCMCDiagnosticTools.ess_rhat`.
+
+!!! note
+    The ordering is not compatible with MCMCDiagnostictools version < 0.2.
 """
 function stack_posterior_matrices(results)
-    @cast _[i, j, k]:= results[k].posterior_matrix[j, i]
+    @cast _[i, k, j]:= results[k].posterior_matrix[j, i]
 end
 
 """
diff --git a/test/test_mcmc.jl b/test/test_mcmc.jl
index b274dd8..f640f31 100644
--- a/test/test_mcmc.jl
+++ b/test/test_mcmc.jl
@@ -74,7 +74,7 @@ end
     D, N, K = 5, 100, 7
     ℓ = multivariate_normal(ones(5))
     results = fill(mcmc_with_warmup(RNG, ℓ, N; reporter = NoProgressReport()), K)
-    @test size(stack_posterior_matrices(results)) == (N, D, K)
+    @test size(stack_posterior_matrices(results)) == (N, K, D)
     @test size(pool_posterior_matrices(results)) == (D, N * K)
 end