diff --git a/Project.toml b/Project.toml
index 572ea1443..25575504f 100644
--- a/Project.toml
+++ b/Project.toml
@@ -18,7 +18,6 @@ Optimisers = "3bd65402-5787-11e9-1adc-39752487f4e2"
 ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
-SimpleUnPack = "ce78b400-467f-4804-87d8-8f486da07d0a"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 
 [weakdeps]
@@ -54,7 +53,6 @@ ProgressMeter = "1.6"
 Random = "1"
 Requires = "1.0"
 ReverseDiff = "1"
-SimpleUnPack = "1.1.0"
 StatsBase = "0.32, 0.33, 0.34"
 Zygote = "0.6"
 julia = "1.7"
@@ -62,12 +60,8 @@ julia = "1.7"
 [extras]
 Bijectors = "76274a88-744f-5084-9051-94815aaf08c4"
 Enzyme = "7da242da-08ed-463a-9acd-ee780be4f1d9"
-ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
-Mooncake = "da2b9cff-9c12-43a0-ae48-6db2b0edb7d6"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
-ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
-Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [targets]
 test = ["Pkg", "Test"]
diff --git a/README.md b/README.md
index 9bd000c75..f734d7714 100644
--- a/README.md
+++ b/README.md
@@ -26,7 +26,6 @@ a `LogDensityProblem` can be implemented as
 
 ```julia
 using LogDensityProblems
-using SimpleUnPack
 
 struct NormalLogNormal{MX,SX,MY,SY}
     μ_x::MX
diff --git a/docs/Project.toml b/docs/Project.toml
index 8dc25a3b1..df1705185 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -10,7 +10,6 @@ LogDensityProblems = "6fdf6af0-433a-55f7-b3ed-c6c6e0b8df7c"
 Optimisers = "3bd65402-5787-11e9-1adc-39752487f4e2"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 QuasiMonteCarlo = "8a4e6c94-4038-4cdc-81c3-7e6ffdb2a71b"
-SimpleUnPack = "ce78b400-467f-4804-87d8-8f486da07d0a"
 StatsFuns = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
 
 [compat]
@@ -25,6 +24,5 @@ LogDensityProblems = "2.1.1"
 Optimisers = "0.3"
 Plots = "1"
 QuasiMonteCarlo = "0.3"
-SimpleUnPack = "1"
 StatsFuns = "1"
 julia = "1.6"
diff --git a/docs/src/examples.md b/docs/src/examples.md
index 15b8907a8..70f8e0ed3 100644
--- a/docs/src/examples.md
+++ b/docs/src/examples.md
@@ -15,7 +15,6 @@ Using the `LogDensityProblems` interface, we the model can be defined as follows
 
 ```@example elboexample
 using LogDensityProblems
-using SimpleUnPack
 
 struct NormalLogNormal{MX,SX,MY,SY}
     μ_x::MX
@@ -25,7 +24,7 @@ struct NormalLogNormal{MX,SX,MY,SY}
 end
 
 function LogDensityProblems.logdensity(model::NormalLogNormal, θ)
-    @unpack μ_x, σ_x, μ_y, Σ_y = model
+    (; μ_x, σ_x, μ_y, Σ_y) = model
     return logpdf(LogNormal(μ_x, σ_x), θ[1]) + logpdf(MvNormal(μ_y, Σ_y), θ[2:end])
 end
 
@@ -59,7 +58,7 @@ Thus, we will use [Bijectors](https://github.com/TuringLang/Bijectors.jl) to mat
 using Bijectors
 
 function Bijectors.bijector(model::NormalLogNormal)
-    @unpack μ_x, σ_x, μ_y, Σ_y = model
+    (; μ_x, σ_x, μ_y, Σ_y) = model
     return Bijectors.Stacked(
         Bijectors.bijector.([LogNormal(μ_x, σ_x), MvNormal(μ_y, Σ_y)]),
         [1:1, 2:(1 + length(μ_y))],
diff --git a/src/AdvancedVI.jl b/src/AdvancedVI.jl
index 1d0c4f502..8433f927e 100644
--- a/src/AdvancedVI.jl
+++ b/src/AdvancedVI.jl
@@ -1,7 +1,6 @@
 
 module AdvancedVI
 
-using SimpleUnPack: @unpack, @pack!
 using Accessors
 
 using Random
diff --git a/src/families/location_scale.jl b/src/families/location_scale.jl
index 22af4b4a9..7885671c7 100644
--- a/src/families/location_scale.jl
+++ b/src/families/location_scale.jl
@@ -56,7 +56,7 @@ function (re::RestructureMeanField)(flat::AbstractVector)
 end
 
 function Optimisers.destructure(q::MvLocationScale{<:Diagonal,D,L,E}) where {D,L,E}
-    @unpack location, scale, dist = q
+    (; location, scale, dist) = q
     flat = vcat(location, diag(scale))
     return flat, RestructureMeanField(q)
 end
@@ -69,19 +69,19 @@ Base.size(q::MvLocationScale) = size(q.location)
 Base.eltype(::Type{<:MvLocationScale{S,D,L,E}}) where {S,D,L,E} = eltype(D)
 
 function StatsBase.entropy(q::MvLocationScale)
-    @unpack location, scale, dist = q
+    (; location, scale, dist) = q
     n_dims = length(location)
     # `convert` is necessary because `entropy` is not type stable upstream
     return n_dims * convert(eltype(location), entropy(dist)) + logdet(scale)
 end
 
 function Distributions.logpdf(q::MvLocationScale, z::AbstractVector{<:Real})
-    @unpack location, scale, dist = q
+    (; location, scale, dist) = q
     return sum(Base.Fix1(logpdf, dist), scale \ (z - location)) - logdet(scale)
 end
 
 function Distributions.rand(q::MvLocationScale)
-    @unpack location, scale, dist = q
+    (; location, scale, dist) = q
     n_dims = length(location)
     return scale * rand(dist, n_dims) + location
 end
@@ -89,7 +89,7 @@ end
 function Distributions.rand(
     rng::AbstractRNG, q::MvLocationScale{S,D,L}, num_samples::Int
 ) where {S,D,L}
-    @unpack location, scale, dist = q
+    (; location, scale, dist) = q
     n_dims = length(location)
     return scale * rand(rng, dist, n_dims, num_samples) .+ location
 end
@@ -98,7 +98,7 @@ end
 function Distributions.rand(
     rng::AbstractRNG, q::MvLocationScale{<:Diagonal,D,L}, num_samples::Int
 ) where {L,D}
-    @unpack location, scale, dist = q
+    (; location, scale, dist) = q
     n_dims = length(location)
     scale_diag = diag(scale)
     return scale_diag .* rand(rng, dist, n_dims, num_samples) .+ location
@@ -107,14 +107,14 @@ end
 function Distributions._rand!(
     rng::AbstractRNG, q::MvLocationScale, x::AbstractVecOrMat{<:Real}
 )
-    @unpack location, scale, dist = q
+    (; location, scale, dist) = q
     rand!(rng, dist, x)
     x[:] = scale * x
     return x .+= location
 end
 
 function Distributions.mean(q::MvLocationScale)
-    @unpack location, scale = q
+    (; location, scale) = q
     return location + scale * Fill(mean(q.dist), length(location))
 end
 
diff --git a/src/families/location_scale_low_rank.jl b/src/families/location_scale_low_rank.jl
index e2044142f..a0849b656 100644
--- a/src/families/location_scale_low_rank.jl
+++ b/src/families/location_scale_low_rank.jl
@@ -52,7 +52,7 @@ Base.size(q::MvLocationScaleLowRank) = size(q.location)
 Base.eltype(::Type{<:MvLocationScaleLowRank{L,SD,SF,D,E}}) where {L,SD,SF,D,E} = eltype(L)
 
 function StatsBase.entropy(q::MvLocationScaleLowRank)
-    @unpack location, scale_diag, scale_factors, dist = q
+    (; location, scale_diag, scale_factors, dist) = q
     n_dims = length(location)
     scale_diag2 = scale_diag .* scale_diag
     UtDinvU = Hermitian(scale_factors' * (scale_factors ./ scale_diag2))
@@ -63,7 +63,7 @@ end
 function Distributions.logpdf(
     q::MvLocationScaleLowRank, z::AbstractVector{<:Real}; non_differntiable::Bool=false
 )
-    @unpack location, scale_diag, scale_factors, dist = q
+    (; location, scale_diag, scale_factors, dist) = q
     μ_base = mean(dist)
     n_dims = length(location)
 
@@ -86,7 +86,7 @@ function Distributions.logpdf(
 end
 
 function Distributions.rand(q::MvLocationScaleLowRank)
-    @unpack location, scale_diag, scale_factors, dist = q
+    (; location, scale_diag, scale_factors, dist) = q
     n_dims = length(location)
     n_factors = size(scale_factors, 2)
     u_diag = rand(dist, n_dims)
@@ -97,7 +97,7 @@ end
 function Distributions.rand(
     rng::AbstractRNG, q::MvLocationScaleLowRank{S,D,L}, num_samples::Int
 ) where {S,D,L}
-    @unpack location, scale_diag, scale_factors, dist = q
+    (; location, scale_diag, scale_factors, dist) = q
     n_dims = length(location)
     n_factors = size(scale_factors, 2)
     u_diag = rand(rng, dist, n_dims, num_samples)
@@ -108,7 +108,7 @@ end
 function Distributions._rand!(
     rng::AbstractRNG, q::MvLocationScaleLowRank, x::AbstractVecOrMat{<:Real}
 )
-    @unpack location, scale_diag, scale_factors, dist = q
+    (; location, scale_diag, scale_factors, dist) = q
 
     rand!(rng, dist, x)
     x[:] = scale_diag .* x
@@ -120,7 +120,7 @@ function Distributions._rand!(
 end
 
 function Distributions.mean(q::MvLocationScaleLowRank)
-    @unpack location, scale_diag, scale_factors = q
+    (; location, scale_diag, scale_factors) = q
     μ = mean(q.dist)
     return location +
            scale_diag .* Fill(μ, length(scale_diag)) +
@@ -128,14 +128,14 @@ function Distributions.mean(q::MvLocationScaleLowRank)
 end
 
 function Distributions.var(q::MvLocationScaleLowRank)
-    @unpack scale_diag, scale_factors = q
+    (; scale_diag, scale_factors) = q
     σ2 = var(q.dist)
     return σ2 *
            (scale_diag .* scale_diag + sum(scale_factors .* scale_factors; dims=2)[:, 1])
 end
 
 function Distributions.cov(q::MvLocationScaleLowRank)
-    @unpack scale_diag, scale_factors = q
+    (; scale_diag, scale_factors) = q
     σ2 = var(q.dist)
     return σ2 * (Diagonal(scale_diag .* scale_diag) + scale_factors * scale_factors')
 end
diff --git a/src/objectives/elbo/repgradelbo.jl b/src/objectives/elbo/repgradelbo.jl
index b8bf63fa8..5b2fd8281 100644
--- a/src/objectives/elbo/repgradelbo.jl
+++ b/src/objectives/elbo/repgradelbo.jl
@@ -86,7 +86,7 @@ function estimate_objective(obj::RepGradELBO, q, prob; n_samples::Int=obj.n_samp
 end
 
 function estimate_repgradelbo_ad_forward(params′, aux)
-    @unpack rng, obj, problem, adtype, restructure, q_stop = aux
+    (; rng, obj, problem, adtype, restructure, q_stop) = aux
     q = restructure_ad_forward(adtype, restructure, params′)
     samples, entropy = reparam_with_entropy(rng, q, q_stop, obj.n_samples, obj.entropy)
     energy = estimate_energy_with_samples(problem, samples)
diff --git a/src/objectives/elbo/scoregradelbo.jl b/src/objectives/elbo/scoregradelbo.jl
index 053c6b3f4..8b9a91c62 100644
--- a/src/objectives/elbo/scoregradelbo.jl
+++ b/src/objectives/elbo/scoregradelbo.jl
@@ -94,7 +94,7 @@ function estimate_objective(obj::ScoreGradELBO, q, prob; n_samples::Int=obj.n_sa
 end
 
 function estimate_scoregradelbo_ad_forward(params′, aux)
-    @unpack rng, obj, problem, adtype, restructure, q_stop = aux
+    (; rng, obj, problem, adtype, restructure, q_stop) = aux
     baseline = compute_control_variate_baseline(
         obj.baseline_history, obj.baseline_window_size
     )
diff --git a/test/Project.toml b/test/Project.toml
index bbf9c4c61..84245c44f 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -15,7 +15,6 @@ PDMats = "90014a1f-27ba-587c-ab20-58faa44d9150"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
-SimpleUnPack = "ce78b400-467f-4804-87d8-8f486da07d0a"
 StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
@@ -39,7 +38,6 @@ Optimisers = "0.2.16, 0.3"
 PDMats = "0.11.7"
 Random = "1"
 ReverseDiff = "1.15.1"
-SimpleUnPack = "1.1.0"
 StableRNGs = "1.0.0"
 Statistics = "1"
 StatsBase = "0.34"
diff --git a/test/inference/repgradelbo_locationscale.jl b/test/inference/repgradelbo_locationscale.jl
index d1f0d7e41..3b78d10be 100644
--- a/test/inference/repgradelbo_locationscale.jl
+++ b/test/inference/repgradelbo_locationscale.jl
@@ -32,7 +32,7 @@ end
         rng = StableRNG(seed)
 
         modelstats = modelconstr(rng, realtype)
-        @unpack model, μ_true, L_true, n_dims, strong_convexity, is_meanfield = modelstats
+        (; model, μ_true, L_true, n_dims, strong_convexity, is_meanfield) = modelstats
 
         T = 1000
         η = 1e-3
diff --git a/test/inference/repgradelbo_locationscale_bijectors.jl b/test/inference/repgradelbo_locationscale_bijectors.jl
index e2a69d62e..1833e2e8e 100644
--- a/test/inference/repgradelbo_locationscale_bijectors.jl
+++ b/test/inference/repgradelbo_locationscale_bijectors.jl
@@ -32,7 +32,7 @@ end
         rng = StableRNG(seed)
 
         modelstats = modelconstr(rng, realtype)
-        @unpack model, μ_true, L_true, n_dims, strong_convexity, is_meanfield = modelstats
+        (; model, μ_true, L_true, n_dims, strong_convexity, is_meanfield) = modelstats
 
         T = 1000
         η = 1e-3
diff --git a/test/inference/scoregradelbo_distributionsad.jl b/test/inference/scoregradelbo_distributionsad.jl
index 9a621b402..e3d78c469 100644
--- a/test/inference/scoregradelbo_distributionsad.jl
+++ b/test/inference/scoregradelbo_distributionsad.jl
@@ -29,7 +29,7 @@ end
         rng = StableRNG(seed)
 
         modelstats = modelconstr(rng, realtype)
-        @unpack model, μ_true, L_true, n_dims, strong_convexity, is_meanfield = modelstats
+        (; model, μ_true, L_true, n_dims, strong_convexity, is_meanfield) = modelstats
 
         T = 1000
         η = 1e-5
diff --git a/test/inference/scoregradelbo_locationscale.jl b/test/inference/scoregradelbo_locationscale.jl
index 753999dee..f035cffe5 100644
--- a/test/inference/scoregradelbo_locationscale.jl
+++ b/test/inference/scoregradelbo_locationscale.jl
@@ -32,7 +32,7 @@ end
         rng = StableRNG(seed)
 
         modelstats = modelconstr(rng, realtype)
-        @unpack model, μ_true, L_true, n_dims, strong_convexity, is_meanfield = modelstats
+        (; model, μ_true, L_true, n_dims, strong_convexity, is_meanfield) = modelstats
 
         T = 1000
         η = 1e-5
diff --git a/test/inference/scoregradelbo_locationscale_bijectors.jl b/test/inference/scoregradelbo_locationscale_bijectors.jl
index bee8234ab..7621f3a76 100644
--- a/test/inference/scoregradelbo_locationscale_bijectors.jl
+++ b/test/inference/scoregradelbo_locationscale_bijectors.jl
@@ -30,7 +30,7 @@ end
         rng = StableRNG(seed)
 
         modelstats = modelconstr(rng, realtype)
-        @unpack model, μ_true, L_true, n_dims, strong_convexity, is_meanfield = modelstats
+        (; model, μ_true, L_true, n_dims, strong_convexity, is_meanfield) = modelstats
 
         T = 1000
         η = 1e-5
diff --git a/test/interface/repgradelbo.jl b/test/interface/repgradelbo.jl
index be835e203..b5bc9d5ec 100644
--- a/test/interface/repgradelbo.jl
+++ b/test/interface/repgradelbo.jl
@@ -7,7 +7,7 @@ using Test
 
     modelstats = normal_meanfield(rng, Float64)
 
-    @unpack model, μ_true, L_true, n_dims, is_meanfield = modelstats
+    (; model, μ_true, L_true, n_dims, is_meanfield) = modelstats
 
     q0 = TuringDiagMvNormal(zeros(Float64, n_dims), ones(Float64, n_dims))
 
@@ -32,7 +32,7 @@ end
     rng = StableRNG(seed)
 
     modelstats = normal_meanfield(rng, Float64)
-    @unpack model, μ_true, L_true, n_dims, is_meanfield = modelstats
+    (; model, μ_true, L_true, n_dims, is_meanfield) = modelstats
 
     ad_backends = [
         ADTypes.AutoForwardDiff(), ADTypes.AutoReverseDiff(), ADTypes.AutoZygote()
diff --git a/test/interface/scoregradelbo.jl b/test/interface/scoregradelbo.jl
index 8a6ebb14f..8b0a34280 100644
--- a/test/interface/scoregradelbo.jl
+++ b/test/interface/scoregradelbo.jl
@@ -7,7 +7,7 @@ using Test
 
     modelstats = normal_meanfield(rng, Float64)
 
-    @unpack model, μ_true, L_true, n_dims, is_meanfield = modelstats
+    (; model, μ_true, L_true, n_dims, is_meanfield) = modelstats
 
     q0 = TuringDiagMvNormal(zeros(Float64, n_dims), ones(Float64, n_dims))
 
diff --git a/test/models/normal.jl b/test/models/normal.jl
index 9fc6ae38a..5826547df 100644
--- a/test/models/normal.jl
+++ b/test/models/normal.jl
@@ -5,7 +5,7 @@ struct TestNormal{M,S}
 end
 
 function LogDensityProblems.logdensity(model::TestNormal, θ)
-    @unpack μ, Σ = model
+    (; μ, Σ) = model
     return logpdf(MvNormal(μ, Σ), θ)
 end
 
diff --git a/test/models/normallognormal.jl b/test/models/normallognormal.jl
index 176aab2f5..00949bc1b 100644
--- a/test/models/normallognormal.jl
+++ b/test/models/normallognormal.jl
@@ -7,7 +7,7 @@ struct NormalLogNormal{MX,SX,MY,SY}
 end
 
 function LogDensityProblems.logdensity(model::NormalLogNormal, θ)
-    @unpack μ_x, σ_x, μ_y, Σ_y = model
+    (; μ_x, σ_x, μ_y, Σ_y) = model
     return logpdf(LogNormal(μ_x, σ_x), θ[1]) + logpdf(MvNormal(μ_y, Σ_y), θ[2:end])
 end
 
@@ -20,7 +20,7 @@ function LogDensityProblems.capabilities(::Type{<:NormalLogNormal})
 end
 
 function Bijectors.bijector(model::NormalLogNormal)
-    @unpack μ_x, σ_x, μ_y, Σ_y = model
+    (; μ_x, σ_x, μ_y, Σ_y) = model
     return Bijectors.Stacked(
         Bijectors.bijector.([LogNormal(μ_x, σ_x), MvNormal(μ_y, Σ_y)]),
         [1:1, 2:(1 + length(μ_y))],
diff --git a/test/runtests.jl b/test/runtests.jl
index 7c0e3129e..ceff9feb4 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -13,7 +13,6 @@ using Optimisers
 using PDMats
 using Pkg
 using Random, StableRNGs
-using SimpleUnPack: @unpack
 using Statistics
 using StatsBase