diff --git a/Project.toml b/Project.toml
index 301e485f7..fafeec2ba 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "DynamicPPL"
 uuid = "366bfd00-2699-11ea-058f-f148b4cae6d8"
-version = "0.31.0"
+version = "0.31.1"
 
 [deps]
 ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
diff --git a/docs/src/api.md b/docs/src/api.md
index d7531af0f..7448812bf 100644
--- a/docs/src/api.md
+++ b/docs/src/api.md
@@ -14,12 +14,6 @@ These statements are rewritten by `@model` as calls of [internal functions](@ref
 @model
 ```
 
-One can nest models and call another model inside the model function with [`@submodel`](@ref).
-
-```@docs
-@submodel
-```
-
 ### Type
 
 A [`Model`](@ref) can be created by calling the model function, as defined by [`@model`](@ref).
@@ -110,6 +104,34 @@ Similarly, we can [`unfix`](@ref) variables, i.e. return them to their original
 unfix
 ```
 
+## Models within models
+
+One can include models and call another model inside the model function with `left ~ to_submodel(model)`.
+
+```@docs
+to_submodel
+```
+
+Note that a `[to_submodel](@ref)` is only sampleable; one cannot compute `logpdf` for its realizations.
+
+In the past, one would instead embed sub-models using [`@submodel`](@ref), which has been deprecated since the introduction of [`to_submodel(model)`](@ref)
+
+```@docs
+@submodel
+```
+
+In the context of including models within models, it's also useful to prefix the variables in sub-models to avoid variable names clashing:
+
+```@docs
+prefix
+```
+
+Under the hood, [`to_submodel`](@ref) makes use of the following method to indicate that the model it's wrapping is a model over its return-values rather than something else
+
+```@docs
+returned(::Model)
+```
+
 ## Utilities
 
 It is possible to manually increase (or decrease) the accumulated log density from within a model function.
@@ -118,10 +140,10 @@ It is possible to manually increase (or decrease) the accumulated log density fr
 @addlogprob!
 ```
 
-Return values of the model function for a collection of samples can be obtained with [`generated_quantities`](@ref).
+Return values of the model function for a collection of samples can be obtained with [`returned(model, chain)`](@ref).
 
 ```@docs
-generated_quantities
+returned(::DynamicPPL.Model, ::NamedTuple)
 ```
 
 For a chain of samples, one can compute the pointwise log-likelihoods of each observed random variable with [`pointwise_loglikelihoods`](@ref). Similarly, the log-densities of the priors using
diff --git a/ext/DynamicPPLMCMCChainsExt.jl b/ext/DynamicPPLMCMCChainsExt.jl
index 8a2679d09..52200171d 100644
--- a/ext/DynamicPPLMCMCChainsExt.jl
+++ b/ext/DynamicPPLMCMCChainsExt.jl
@@ -43,7 +43,7 @@ function DynamicPPL.varnames(c::MCMCChains.Chains)
 end
 
 """
-    generated_quantities(model::Model, chain::MCMCChains.Chains)
+    returned(model::Model, chain::MCMCChains.Chains)
 
 Execute `model` for each of the samples in `chain` and return an array of the values
 returned by the `model` for each sample.
@@ -63,12 +63,12 @@ m = demo(data)
 chain = sample(m, alg, n)
 # To inspect the `interesting_quantity(θ, x)` where `θ` is replaced by samples
 # from the posterior/`chain`:
-generated_quantities(m, chain) # <= results in a `Vector` of returned values
+returned(m, chain) # <= results in a `Vector` of returned values
                                #    from `interesting_quantity(θ, x)`
 ```
 ## Concrete (and simple)
 ```julia
-julia> using DynamicPPL, Turing
+julia> using Turing
 
 julia> @model function demo(xs)
            s ~ InverseGamma(2, 3)
@@ -87,7 +87,7 @@ julia> model = demo(randn(10));
 
 julia> chain = sample(model, MH(), 10);
 
-julia> generated_quantities(model, chain)
+julia> returned(model, chain)
 10×1 Array{Tuple{Float64},2}:
  (2.1964758025119338,)
  (2.1964758025119338,)
@@ -101,9 +101,7 @@ julia> generated_quantities(model, chain)
  (-0.16489786710222099,)
 ```
 """
-function DynamicPPL.generated_quantities(
-    model::DynamicPPL.Model, chain_full::MCMCChains.Chains
-)
+function DynamicPPL.returned(model::DynamicPPL.Model, chain_full::MCMCChains.Chains)
     chain = MCMCChains.get_sections(chain_full, :parameters)
     varinfo = DynamicPPL.VarInfo(model)
     iters = Iterators.product(1:size(chain, 1), 1:size(chain, 3))
diff --git a/src/DynamicPPL.jl b/src/DynamicPPL.jl
index a5d178125..d4e13d456 100644
--- a/src/DynamicPPL.jl
+++ b/src/DynamicPPL.jl
@@ -86,7 +86,6 @@ export AbstractVarInfo,
     Model,
     getmissings,
     getargnames,
-    generated_quantities,
     extract_priors,
     values_as_in_model,
     # Samplers
@@ -122,6 +121,9 @@ export AbstractVarInfo,
     decondition,
     fix,
     unfix,
+    prefix,
+    returned,
+    to_submodel,
     # Convenience macros
     @addlogprob!,
     @submodel,
@@ -130,7 +132,8 @@ export AbstractVarInfo,
     check_model_and_trace,
     # Deprecated.
     @logprob_str,
-    @prob_str
+    @prob_str,
+    generated_quantities
 
 # Reexport
 using Distributions: loglikelihood
@@ -196,6 +199,8 @@ include("values_as_in_model.jl")
 include("debug_utils.jl")
 using .DebugUtils
 
+include("deprecated.jl")
+
 if !isdefined(Base, :get_extension)
     using Requires
 end
diff --git a/src/compiler.jl b/src/compiler.jl
index 90220cbf5..c67da6f95 100644
--- a/src/compiler.jl
+++ b/src/compiler.jl
@@ -178,6 +178,11 @@ function check_tilde_rhs(@nospecialize(x))
 end
 check_tilde_rhs(x::Distribution) = x
 check_tilde_rhs(x::AbstractArray{<:Distribution}) = x
+check_tilde_rhs(x::ReturnedModelWrapper) = x
+function check_tilde_rhs(x::Sampleable{<:Any,AutoPrefix}) where {AutoPrefix}
+    model = check_tilde_rhs(x.model)
+    return Sampleable{typeof(model),AutoPrefix}(model)
+end
 
 """
     unwrap_right_vn(right, vn)
diff --git a/src/context_implementations.jl b/src/context_implementations.jl
index 489c64c57..462012676 100644
--- a/src/context_implementations.jl
+++ b/src/context_implementations.jl
@@ -103,8 +103,17 @@ By default, calls `tilde_assume(context, right, vn, vi)` and accumulates the log
 probability of `vi` with the returned value.
 """
 function tilde_assume!!(context, right, vn, vi)
-    value, logp, vi = tilde_assume(context, right, vn, vi)
-    return value, acclogp_assume!!(context, vi, logp)
+    return if is_rhs_model(right)
+        # Prefix the variables using the `vn`.
+        rand_like!!(
+            right,
+            should_auto_prefix(right) ? PrefixContext{Symbol(vn)}(context) : context,
+            vi,
+        )
+    else
+        value, logp, vi = tilde_assume(context, right, vn, vi)
+        value, acclogp_assume!!(context, vi, logp)
+    end
 end
 
 # observe
@@ -159,6 +168,11 @@ Falls back to `tilde_observe!!(context, right, left, vi)` ignoring the informati
 and indices; if needed, these can be accessed through this function, though.
 """
 function tilde_observe!!(context, right, left, vname, vi)
+    is_rhs_model(right) && throw(
+        ArgumentError(
+            "`~` with a model on the right-hand side of an observe statement is not supported",
+        ),
+    )
     return tilde_observe!!(context, right, left, vi)
 end
 
@@ -172,6 +186,11 @@ By default, calls `tilde_observe(context, right, left, vi)` and accumulates the
 probability of `vi` with the returned value.
 """
 function tilde_observe!!(context, right, left, vi)
+    is_rhs_model(right) && throw(
+        ArgumentError(
+            "`~` with a model on the right-hand side of an observe statement is not supported",
+        ),
+    )
     logp, vi = tilde_observe(context, right, left, vi)
     return left, acclogp_observe!!(context, vi, logp)
 end
@@ -321,8 +340,13 @@ model inputs), accumulate the log probability, and return the sampled value and
 Falls back to `dot_tilde_assume(context, right, left, vn, vi)`.
 """
 function dot_tilde_assume!!(context, right, left, vn, vi)
+    is_rhs_model(right) && throw(
+        ArgumentError(
+            "`.~` with a model on the right-hand side is not supported; please use `~`"
+        ),
+    )
     value, logp, vi = dot_tilde_assume(context, right, left, vn, vi)
-    return value, acclogp_assume!!(context, vi, logp), vi
+    return value, acclogp_assume!!(context, vi, logp)
 end
 
 # `dot_assume`
@@ -573,6 +597,11 @@ Falls back to `dot_tilde_observe!!(context, right, left, vi)` ignoring the infor
 name and indices; if needed, these can be accessed through this function, though.
 """
 function dot_tilde_observe!!(context, right, left, vn, vi)
+    is_rhs_model(right) && throw(
+        ArgumentError(
+            "`~` with a model on the right-hand side of an observe statement is not supported",
+        ),
+    )
     return dot_tilde_observe!!(context, right, left, vi)
 end
 
@@ -585,6 +614,11 @@ probability, and return the observed value and updated `vi`.
 Falls back to `dot_tilde_observe(context, right, left, vi)`.
 """
 function dot_tilde_observe!!(context, right, left, vi)
+    is_rhs_model(right) && throw(
+        ArgumentError(
+            "`~` with a model on the right-hand side of an observe statement is not supported",
+        ),
+    )
     logp, vi = dot_tilde_observe(context, right, left, vi)
     return left, acclogp_observe!!(context, vi, logp)
 end
diff --git a/src/contexts.jl b/src/contexts.jl
index 5da4208b5..9eb3d5ccb 100644
--- a/src/contexts.jl
+++ b/src/contexts.jl
@@ -281,6 +281,34 @@ function prefix(::PrefixContext{Prefix}, vn::VarName{Sym}) where {Prefix,Sym}
     end
 end
 
+"""
+    prefix(model::Model, x)
+
+Return `model` but with all random variables prefixed by `x`.
+
+If `x` is known at compile-time, use `Val{x}()` to avoid runtime overheads for prefixing.
+
+# Examples
+
+```jldoctest
+julia> using DynamicPPL: prefix
+
+julia> @model demo() = x ~ Dirac(1)
+demo (generic function with 2 methods)
+
+julia> rand(prefix(demo(), :my_prefix))
+(var"my_prefix.x" = 1,)
+
+julia> # One can also use `Val` to avoid runtime overheads.
+       rand(prefix(demo(), Val(:my_prefix)))
+(var"my_prefix.x" = 1,)
+```
+"""
+prefix(model::Model, x) = contextualize(model, PrefixContext{Symbol(x)}(model.context))
+function prefix(model::Model, ::Val{x}) where {x}
+    return contextualize(model, PrefixContext{Symbol(x)}(model.context))
+end
+
 struct ConditionContext{Values,Ctx<:AbstractContext} <: AbstractContext
     values::Values
     context::Ctx
diff --git a/src/deprecated.jl b/src/deprecated.jl
new file mode 100644
index 000000000..0bcaae9b7
--- /dev/null
+++ b/src/deprecated.jl
@@ -0,0 +1 @@
+@deprecate generated_quantities(model, params) returned(model, params)
diff --git a/src/model.jl b/src/model.jl
index 2a1a6db88..c0489b089 100644
--- a/src/model.jl
+++ b/src/model.jl
@@ -223,15 +223,16 @@ true
 ## Nested models
 
 `condition` of course also supports the use of nested models through
-the use of [`@submodel`](@ref).
+the use of [`to_submodel`](@ref).
 
 ```jldoctest condition
 julia> @model demo_inner() = m ~ Normal()
 demo_inner (generic function with 2 methods)
 
 julia> @model function demo_outer()
-           @submodel m = demo_inner()
-           return m
+           # By default, `to_submodel` prefixes the variables using the left-hand side of `~`.
+           inner ~ to_submodel(demo_inner())
+           return inner
        end
 demo_outer (generic function with 2 methods)
 
@@ -240,63 +241,28 @@ julia> model = demo_outer();
 julia> model() ≠ 1.0
 true
 
-julia> conditioned_model = model | (m = 1.0, );
+julia> # To condition the variable inside `demo_inner` we need to refer to it as `inner.m`.
+       conditioned_model = model | (var"inner.m" = 1.0, );
 
 julia> conditioned_model()
 1.0
-```
-
-But one needs to be careful when prefixing variables in the nested models:
-
-```jldoctest condition
-julia> @model function demo_outer_prefix()
-           @submodel prefix="inner" m = demo_inner()
-           return m
-       end
-demo_outer_prefix (generic function with 2 methods)
-
-julia> # (×) This doesn't work now!
-       conditioned_model = demo_outer_prefix() | (m = 1.0, );
-
-julia> conditioned_model() == 1.0
-false
 
-julia> # (✓) `m` in `demo_inner` is referred to as `inner.m` internally, so we do:
-       conditioned_model = demo_outer_prefix() | (var"inner.m" = 1.0, );
+julia> # However, it's not possible to condition `inner` directly.
+       conditioned_model_fail = model | (inner = 1.0, );
 
-julia> conditioned_model()
-1.0
-
-julia> # Note that the above `var"..."` is just standard Julia syntax:
-       keys((var"inner.m" = 1.0, ))
-(Symbol("inner.m"),)
+julia> conditioned_model_fail()
+ERROR: ArgumentError: `~` with a model on the right-hand side of an observe statement is not supported
+[...]
 ```
 
 And similarly when using `Dict`:
 
 ```jldoctest condition
-julia> conditioned_model_dict = demo_outer_prefix() | (@varname(var"inner.m") => 1.0);
+julia> conditioned_model_dict = model | (@varname(var"inner.m") => 1.0);
 
 julia> conditioned_model_dict()
 1.0
 ```
-
-The difference is maybe more obvious once we look at how these different
-in their trace/`VarInfo`:
-
-```jldoctest condition
-julia> keys(VarInfo(demo_outer()))
-1-element Vector{VarName{:m, typeof(identity)}}:
- m
-
-julia> keys(VarInfo(demo_outer_prefix()))
-1-element Vector{VarName{Symbol("inner.m"), typeof(identity)}}:
- inner.m
-```
-
-From this we can tell what the correct way to condition `m` within `demo_inner`
-is in the two different models.
-
 """
 AbstractPPL.condition(model::Model; values...) = condition(model, NamedTuple(values))
 function AbstractPPL.condition(model::Model, value, values...)
@@ -578,15 +544,15 @@ true
 ## Nested models
 
 `fix` of course also supports the use of nested models through
-the use of [`@submodel`](@ref).
+the use of [`to_submodel`](@ref), similar to [`condition`](@ref).
 
 ```jldoctest fix
 julia> @model demo_inner() = m ~ Normal()
 demo_inner (generic function with 2 methods)
 
 julia> @model function demo_outer()
-           @submodel m = demo_inner()
-           return m
+           inner ~ to_submodel(demo_inner())
+           return inner
        end
 demo_outer (generic function with 2 methods)
 
@@ -595,63 +561,36 @@ julia> model = demo_outer();
 julia> model() ≠ 1.0
 true
 
-julia> fixed_model = model | (m = 1.0, );
+julia> fixed_model = fix(model, var"inner.m" = 1.0, );
 
 julia> fixed_model()
 1.0
 ```
 
-But one needs to be careful when prefixing variables in the nested models:
-
-```jldoctest fix
-julia> @model function demo_outer_prefix()
-           @submodel prefix="inner" m = demo_inner()
-           return m
-       end
-demo_outer_prefix (generic function with 2 methods)
-
-julia> # (×) This doesn't work now!
-       fixed_model = demo_outer_prefix() | (m = 1.0, );
-
-julia> fixed_model() == 1.0
-false
+However, unlike [`condition`](@ref), `fix` can also be used to fix the
+return-value of the submodel:
 
-julia> # (✓) `m` in `demo_inner` is referred to as `inner.m` internally, so we do:
-       fixed_model = demo_outer_prefix() | (var"inner.m" = 1.0, );
+```julia
+julia> fixed_model = fix(model, inner = 2.0,);
 
 julia> fixed_model()
-1.0
-
-julia> # Note that the above `var"..."` is just standard Julia syntax:
-       keys((var"inner.m" = 1.0, ))
-(Symbol("inner.m"),)
+2.0
 ```
 
 And similarly when using `Dict`:
 
 ```jldoctest fix
-julia> fixed_model_dict = demo_outer_prefix() | (@varname(var"inner.m") => 1.0);
+julia> fixed_model_dict = fix(model, @varname(var"inner.m") => 1.0);
 
 julia> fixed_model_dict()
 1.0
-```
 
-The difference is maybe more obvious once we look at how these different
-in their trace/`VarInfo`:
+julia> fixed_model_dict = fix(model, @varname(inner) => 2.0);
 
-```jldoctest fix
-julia> keys(VarInfo(demo_outer()))
-1-element Vector{VarName{:m, typeof(identity)}}:
- m
-
-julia> keys(VarInfo(demo_outer_prefix()))
-1-element Vector{VarName{Symbol("inner.m"), typeof(identity)}}:
- inner.m
+julia> fixed_model_dict()
+2.0
 ```
 
-From this we can tell what the correct way to fix `m` within `demo_inner`
-is in the two different models.
-
 ## Difference from `condition`
 
 A very similar functionality is also provided by [`condition`](@ref) which,
@@ -1051,7 +990,9 @@ function Base.rand(rng::Random.AbstractRNG, ::Type{T}, model::Model) where {T}
         evaluate!!(
             model,
             SimpleVarInfo{Float64}(OrderedDict()),
-            SamplingContext(rng, SampleFromPrior(), model.context),
+            # NOTE: Use `leafcontext` here so we a) avoid overriding the leaf context of `model`,
+            # and b) avoid double-stacking the parent contexts.
+            SamplingContext(rng, SampleFromPrior(), leafcontext(model.context)),
         ),
     )
     return values_as(x, T)
@@ -1204,9 +1145,9 @@ function Distributions.loglikelihood(model::Model, chain::AbstractMCMC.AbstractC
 end
 
 """
-    generated_quantities(model::Model, parameters::NamedTuple)
-    generated_quantities(model::Model, values, keys)
-    generated_quantities(model::Model, values, keys)
+    returned(model::Model, parameters::NamedTuple)
+    returned(model::Model, values, keys)
+    returned(model::Model, values, keys)
 
 Execute `model` with variables `keys` set to `values` and return the values returned by the `model`.
 
@@ -1231,18 +1172,257 @@ julia> model = demo(randn(10));
 
 julia> parameters = (; s = 1.0, m_shifted=10.0);
 
-julia> generated_quantities(model, parameters)
+julia> returned(model, parameters)
 (0.0,)
 
-julia> generated_quantities(model, values(parameters), keys(parameters))
+julia> returned(model, values(parameters), keys(parameters))
 (0.0,)
 ```
 """
-function generated_quantities(model::Model, parameters::NamedTuple)
+function returned(model::Model, parameters::NamedTuple)
     fixed_model = fix(model, parameters)
     return fixed_model()
 end
 
-function generated_quantities(model::Model, values, keys)
-    return generated_quantities(model, NamedTuple{keys}(values))
+function returned(model::Model, values, keys)
+    return returned(model, NamedTuple{keys}(values))
 end
+
+"""
+    is_rhs_model(x)
+
+Return `true` if `x` is a model or model wrapper, and `false` otherwise.
+"""
+is_rhs_model(x) = false
+
+"""
+    Distributional
+
+Abstract type for type indicating that something is "distributional".
+"""
+abstract type Distributional end
+
+"""
+    should_auto_prefix(distributional)
+
+Return `true` if the `distributional` should use automatic prefixing, and `false` otherwise.
+"""
+function should_auto_prefix end
+
+"""
+    is_rhs_model(x)
+
+Return `true` if the `distributional` is a model, and `false` otherwise.
+"""
+function is_rhs_model end
+
+"""
+    Sampleable{M} <: Distributional
+
+A wrapper around a model indicating it is sampleable.
+"""
+struct Sampleable{M,AutoPrefix} <: Distributional
+    model::M
+end
+
+should_auto_prefix(::Sampleable{<:Any,AutoPrefix}) where {AutoPrefix} = AutoPrefix
+is_rhs_model(x::Sampleable) = is_rhs_model(x.model)
+
+# TODO: Export this if it end up having a purpose beyond `to_submodel`.
+"""
+    to_sampleable(model[, auto_prefix])
+
+Return a wrapper around `model` indicating it is sampleable.
+
+# Arguments
+- `model::Model`: the model to wrap.
+- `auto_prefix::Bool`: whether to prefix the variables in the model. Default: `true`.
+"""
+to_sampleable(model, auto_prefix::Bool=true) = Sampleable{typeof(model),auto_prefix}(model)
+
+"""
+    rand_like!!(model_wrap, context, varinfo)
+
+Returns a tuple with the first element being the realization and the second the updated varinfo.
+
+# Arguments
+- `model_wrap::ReturnedModelWrapper`: the wrapper of the model to use.
+- `context::AbstractContext`: the context to use for evaluation.
+- `varinfo::AbstractVarInfo`: the varinfo to use for evaluation.
+    """
+function rand_like!!(
+    model_wrap::Sampleable, context::AbstractContext, varinfo::AbstractVarInfo
+)
+    return rand_like!!(model_wrap.model, context, varinfo)
+end
+
+"""
+    ReturnedModelWrapper
+
+A wrapper around a model indicating it is a model over its return values.
+
+This should rarely be constructed explicitly; see [`returned(model)`](@ref) instead.
+"""
+struct ReturnedModelWrapper{M<:Model}
+    model::M
+end
+
+is_rhs_model(::ReturnedModelWrapper) = true
+
+function rand_like!!(
+    model_wrap::ReturnedModelWrapper, context::AbstractContext, varinfo::AbstractVarInfo
+)
+    # Return's the value and the (possibly mutated) varinfo.
+    return _evaluate!!(model_wrap.model, varinfo, context)
+end
+
+"""
+    returned(model)
+
+Return a `model` wrapper indicating that it is a model over its return-values.
+"""
+returned(model::Model) = ReturnedModelWrapper(model)
+
+"""
+    to_submodel(model::Model[, auto_prefix::Bool])
+
+Return a model wrapper indicating that it is a sampleable model over the return-values.
+
+This is mainly meant to be used on the right-hand side of a `~` operator to indicate that
+the model can be sampled from but not necessarily evaluated for its log density.
+
+!!! warning
+    Note that some other operations that one typically associate with expressions of the form
+    `left ~ right` such as [`condition`](@ref), will also not work with `to_submodel`.
+
+!!! warning
+    To avoid variable names clashing between models, it is recommend leave argument `auto_prefix` equal to `true`.
+    If one does not use automatic prefixing, then it's recommended to use [`prefix(::Model, input)`](@ref) explicitly.
+
+# Arguments
+- `model::Model`: the model to wrap.
+- `auto_prefix::Bool`: whether to automatically prefix the variables in the model using the left-hand
+    side of the `~` statement. Default: `true`.
+
+# Examples
+
+## Simple example
+```jldoctest submodel-to_submodel; setup=:(using Distributions)
+julia> @model function demo1(x)
+           x ~ Normal()
+           return 1 + abs(x)
+       end;
+
+julia> @model function demo2(x, y)
+            a ~ to_submodel(demo1(x))
+            return y ~ Uniform(0, a)
+       end;
+```
+
+When we sample from the model `demo2(missing, 0.4)` random variable `x` will be sampled:
+```jldoctest submodel-to_submodel
+julia> vi = VarInfo(demo2(missing, 0.4));
+
+julia> @varname(var\"a.x\") in keys(vi)
+true
+```
+
+The variable `a` is not tracked. However, it will be assigned the return value of `demo1`,
+and can be used in subsequent lines of the model, as shown above.
+```jldoctest submodel-to_submodel
+julia> @varname(a) in keys(vi)
+false
+```
+
+We can check that the log joint probability of the model accumulated in `vi` is correct:
+
+```jldoctest submodel-to_submodel
+julia> x = vi[@varname(var\"a.x\")];
+
+julia> getlogp(vi) ≈ logpdf(Normal(), x) + logpdf(Uniform(0, 1 + abs(x)), 0.4)
+true
+```
+
+## Without automatic prefixing
+As mentioned earlier, by default, the `auto_prefix` argument specifies whether to automatically
+prefix the variables in the submodel. If `auto_prefix=false`, then the variables in the submodel
+will not be prefixed.
+```jldoctest submodel-to_submodel-prefix; setup=:(using Distributions)
+julia> @model function demo1(x)
+           x ~ Normal()
+           return 1 + abs(x)
+       end;
+
+julia> @model function demo2_no_prefix(x, z)
+            a ~ to_submodel(demo1(x), false)
+            return z ~ Uniform(-a, 1)
+       end;
+
+julia> vi = VarInfo(demo2_no_prefix(missing, 0.4));
+
+julia> @varname(x) in keys(vi)  # here we just use `x` instead of `a.x`
+true
+```
+However, not using prefixing is generally not recommended as it can lead to variable name clashes
+unless one is careful. For example, if we're re-using the same model twice in a model, not using prefixing
+will lead to variable name clashes: However, one can manually prefix using the [`prefix(::Model, input)`](@ref):
+```jldoctest submodel-to_submodel-prefix
+julia> @model function demo2(x, y, z)
+            a ~ to_submodel(prefix(demo1(x), :sub1), false)
+            b ~ to_submodel(prefix(demo1(y), :sub2), false)
+            return z ~ Uniform(-a, b)
+       end;
+
+julia> vi = VarInfo(demo2(missing, missing, 0.4));
+
+julia> @varname(var"sub1.x") in keys(vi)
+true
+
+julia> @varname(var"sub2.x") in keys(vi)
+true
+```
+
+Variables `a` and `b` are not tracked, but are assigned the return values of the respective
+calls to `demo1`:
+```jldoctest submodel-to_submodel-prefix
+julia> @varname(a) in keys(vi)
+false
+
+julia> @varname(b) in keys(vi)
+false
+```
+
+We can check that the log joint probability of the model accumulated in `vi` is correct:
+
+```jldoctest submodel-to_submodel-prefix
+julia> sub1_x = vi[@varname(var"sub1.x")];
+
+julia> sub2_x = vi[@varname(var"sub2.x")];
+
+julia> logprior = logpdf(Normal(), sub1_x) + logpdf(Normal(), sub2_x);
+
+julia> loglikelihood = logpdf(Uniform(-1 - abs(sub1_x), 1 + abs(sub2_x)), 0.4);
+
+julia> getlogp(vi) ≈ logprior + loglikelihood
+true
+```
+
+## Usage as likelihood is illegal
+
+Note that it is illegal to use a `to_submodel` model as a likelihood in another model:
+
+```jldoctest submodel-to_submodel-illegal; setup=:(using Distributions)
+julia> @model inner() = x ~ Normal()
+inner (generic function with 2 methods)
+
+julia> @model illegal_likelihood() = a ~ to_submodel(inner())
+illegal_likelihood (generic function with 2 methods)
+
+julia> model = illegal_likelihood() | (a = 1.0,);
+
+julia> model()
+ERROR: ArgumentError: `~` with a model on the right-hand side of an observe statement is not supported
+[...]
+"""
+to_submodel(model::Model, auto_prefix::Bool=true) =
+    to_sampleable(returned(model), auto_prefix)
diff --git a/src/submodel_macro.jl b/src/submodel_macro.jl
index 050bf31fc..e5a8e0617 100644
--- a/src/submodel_macro.jl
+++ b/src/submodel_macro.jl
@@ -4,6 +4,10 @@
 
 Run a Turing `model` nested inside of a Turing model.
 
+!!! warning
+    This is deprecated and will be removed in a future release.
+    Use `left ~ to_submodel(model)` instead (see [`to_submodel`](@ref)).
+
 # Examples
 
 ```jldoctest submodel; setup=:(using Distributions)
@@ -21,6 +25,9 @@ julia> @model function demo2(x, y)
 When we sample from the model `demo2(missing, 0.4)` random variable `x` will be sampled:
 ```jldoctest submodel
 julia> vi = VarInfo(demo2(missing, 0.4));
+┌ Warning: `@submodel model` and `@submodel prefix=... model` are deprecated; see `to_submodel` for the up-to-date syntax.
+│   caller = ip:0x0
+└ @ Core :-1
 
 julia> @varname(x) in keys(vi)
 true
@@ -62,6 +69,10 @@ Valid expressions for `prefix=...` are:
 The prefix makes it possible to run the same Turing model multiple times while
 keeping track of all random variables correctly.
 
+!!! warning
+    This is deprecated and will be removed in a future release.
+    Use `left ~ to_submodel(model)` instead (see [`to_submodel(model)`](@ref)).
+
 # Examples
 ## Example models
 ```jldoctest submodelprefix; setup=:(using Distributions)
@@ -81,6 +92,9 @@ When we sample from the model `demo2(missing, missing, 0.4)` random variables `s
 `sub2.x` will be sampled:
 ```jldoctest submodelprefix
 julia> vi = VarInfo(demo2(missing, missing, 0.4));
+┌ Warning: `@submodel model` and `@submodel prefix=... model` are deprecated; see `to_submodel` for the up-to-date syntax.
+│   caller = ip:0x0
+└ @ Core :-1
 
 julia> @varname(var"sub1.x") in keys(vi)
 true
@@ -124,6 +138,9 @@ julia> # When `prefix` is unspecified, no prefix is used.
 submodel_noprefix (generic function with 2 methods)
 
 julia> @varname(x) in keys(VarInfo(submodel_noprefix()))
+┌ Warning: `@submodel model` and `@submodel prefix=... model` are deprecated; see `to_submodel` for the up-to-date syntax.
+│   caller = ip:0x0
+└ @ Core :-1
 true
 
 julia> # Explicitely don't use any prefix.
@@ -131,6 +148,9 @@ julia> # Explicitely don't use any prefix.
 submodel_prefix_false (generic function with 2 methods)
 
 julia> @varname(x) in keys(VarInfo(submodel_prefix_false()))
+┌ Warning: `@submodel model` and `@submodel prefix=... model` are deprecated; see `to_submodel` for the up-to-date syntax.
+│   caller = ip:0x0
+└ @ Core :-1
 true
 
 julia> # Automatically determined from `a`.
@@ -138,6 +158,9 @@ julia> # Automatically determined from `a`.
 submodel_prefix_true (generic function with 2 methods)
 
 julia> @varname(var"a.x") in keys(VarInfo(submodel_prefix_true()))
+┌ Warning: `@submodel model` and `@submodel prefix=... model` are deprecated; see `to_submodel` for the up-to-date syntax.
+│   caller = ip:0x0
+└ @ Core :-1
 true
 
 julia> # Using a static string.
@@ -145,6 +168,9 @@ julia> # Using a static string.
 submodel_prefix_string (generic function with 2 methods)
 
 julia> @varname(var"my prefix.x") in keys(VarInfo(submodel_prefix_string()))
+┌ Warning: `@submodel model` and `@submodel prefix=... model` are deprecated; see `to_submodel` for the up-to-date syntax.
+│   caller = ip:0x0
+└ @ Core :-1
 true
 
 julia> # Using string interpolation.
@@ -152,6 +178,9 @@ julia> # Using string interpolation.
 submodel_prefix_interpolation (generic function with 2 methods)
 
 julia> @varname(var"inner.x") in keys(VarInfo(submodel_prefix_interpolation()))
+┌ Warning: `@submodel model` and `@submodel prefix=... model` are deprecated; see `to_submodel` for the up-to-date syntax.
+│   caller = ip:0x0
+└ @ Core :-1
 true
 
 julia> # Or using some arbitrary expression.
@@ -159,6 +188,9 @@ julia> # Or using some arbitrary expression.
 submodel_prefix_expr (generic function with 2 methods)
 
 julia> @varname(var"3.x") in keys(VarInfo(submodel_prefix_expr()))
+┌ Warning: `@submodel model` and `@submodel prefix=... model` are deprecated; see `to_submodel` for the up-to-date syntax.
+│   caller = ip:0x0
+└ @ Core :-1
 true
 
 julia> # (×) Automatic prefixing without a left-hand side expression does not work!
@@ -207,6 +239,8 @@ function prefix_submodel_context(prefix::Bool, ctx)
     return ctx
 end
 
+const SUBMODEL_DEPWARN_MSG = "`@submodel model` and `@submodel prefix=... model` are deprecated; see `to_submodel` for the up-to-date syntax."
+
 function submodel(prefix_expr, expr, ctx=esc(:__context__))
     prefix_left, prefix = getargs_assignment(prefix_expr)
     if prefix_left !== :prefix
@@ -225,6 +259,9 @@ function submodel(prefix_expr, expr, ctx=esc(:__context__))
     return if args_assign === nothing
         ctx = prefix_submodel_context(prefix, ctx)
         quote
+            # Raise deprecation warning to let user know that we recommend using `left ~ to_submodel(model)`.
+            $(Base.depwarn)(SUBMODEL_DEPWARN_MSG, Symbol("@submodel"))
+
             $retval, $(esc(:__varinfo__)) = $(_evaluate!!)(
                 $(esc(expr)), $(esc(:__varinfo__)), $(ctx)
             )
@@ -241,6 +278,9 @@ function submodel(prefix_expr, expr, ctx=esc(:__context__))
             )
         end
         quote
+            # Raise deprecation warning to let user know that we recommend using `left ~ to_submodel(model)`.
+            $(Base.depwarn)(SUBMODEL_DEPWARN_MSG, Symbol("@submodel"))
+
             $retval, $(esc(:__varinfo__)) = $(_evaluate!!)(
                 $(esc(R)), $(esc(:__varinfo__)), $(ctx)
             )
diff --git a/test/compiler.jl b/test/compiler.jl
index f2d7e5852..977c1156c 100644
--- a/test/compiler.jl
+++ b/test/compiler.jl
@@ -382,6 +382,27 @@ module Issue537 end
         @test demo2()() == 42
     end
 
+    @testset "@submodel is deprecated" begin
+        @model inner() = x ~ Normal()
+        @model outer() = @submodel x = inner()
+        @test_logs(
+            (
+                :warn,
+                "`@submodel model` and `@submodel prefix=... model` are deprecated; see `to_submodel` for the up-to-date syntax.",
+            ),
+            outer()()
+        )
+
+        @model outer_with_prefix() = @submodel prefix = "sub" x = inner()
+        @test_logs(
+            (
+                :warn,
+                "`@submodel model` and `@submodel prefix=... model` are deprecated; see `to_submodel` for the up-to-date syntax.",
+            ),
+            outer_with_prefix()()
+        )
+    end
+
     @testset "submodel" begin
         # No prefix, 1 level.
         @model function demo1(x)
@@ -469,7 +490,7 @@ module Issue537 end
             num_steps = length(y[1])
             num_obs = length(y)
             @inbounds for i in 1:num_obs
-                @submodel prefix = "ar1_$i" x = AR1(num_steps, α, μ, σ)
+                x ~ to_submodel(prefix(AR1(num_steps, α, μ, σ), "ar1_$i"), false)
                 y[i] ~ MvNormal(x, 0.01 * I)
             end
         end
diff --git a/test/debug_utils.jl b/test/debug_utils.jl
index 9ca1bc1ba..dfa46affc 100644
--- a/test/debug_utils.jl
+++ b/test/debug_utils.jl
@@ -45,14 +45,16 @@
         @testset "submodel" begin
             @model ModelInner() = x ~ Normal()
             @model function ModelOuterBroken()
-                @submodel z = ModelInner()
+                # Without automatic prefixing => `x` s used twice.
+                z ~ to_submodel(ModelInner(), false)
                 return x ~ Normal()
             end
             model = ModelOuterBroken()
             @test_throws ErrorException check_model(model; error_on_failure=true)
 
             @model function ModelOuterWorking()
-                @submodel prefix = true z = ModelInner()
+                # With automatic prefixing => `x` is not duplicated.
+                z ~ to_submodel(ModelInner())
                 x ~ Normal()
                 return z
             end
diff --git a/test/ext/DynamicPPLMCMCChainsExt.jl b/test/ext/DynamicPPLMCMCChainsExt.jl
index c19bf6f2d..e117c4bbc 100644
--- a/test/ext/DynamicPPLMCMCChainsExt.jl
+++ b/test/ext/DynamicPPLMCMCChainsExt.jl
@@ -3,7 +3,7 @@
     model = demo()
 
     chain = MCMCChains.Chains(randn(1000, 2, 1), [:x, :y], Dict(:internals => [:y]))
-    chain_generated = @test_nowarn generated_quantities(model, chain)
+    chain_generated = @test_nowarn returned(model, chain)
     @test size(chain_generated) == (1000, 1)
     @test mean(chain_generated) ≈ 0 atol = 0.1
 end