Tracking Extra Quantities

Often, there are quantities in models that we might be interested in viewing the values of, but which are not random variables in the model that are explicitly drawn from a distribution.

As a motivating example, the most natural parameterization for a model might not be the most computationally feasible. Consider the following (efficiently reparametrized) implementation of Neal’s funnel (Neal, 2003):

using Turing
setprogress!(false)

@model function Neal()
    # Raw draws
    y_raw ~ Normal(0, 1)
    x_raw ~ arraydist([Normal(0, 1) for i in 1:9])

    # Transform:
    y = 3 * y_raw
    x = exp.(y ./ 2) .* x_raw
    return nothing
end

Precompiling Turing...
    805.4 ms  ? OptimizationBase
   1404.9 ms  ? Optimization
   2163.4 ms  ? OptimizationOptimJL
Info Given Turing was explicitly requested, output will be shown live 
WARNING: redefinition of constant OptimizationBase.OPTIMIZER_MISSING_ERROR_MESSAGE. This may fail, cause incorrect answers, or produce other errors.
WARNING: Method definition (::Type{OptimizationBase.OptimizerMissingError})(Any) in module OptimizationBase at /home/runner/.julia/packages/OptimizationBase/sfIfa/src/solve.jl:23 overwritten at /home/runner/.julia/packages/OptimizationBase/sfIfa/src/solve.jl:177.
ERROR: Method overwriting is not permitted during Module precompilation. Use `__precompile__(false)` to opt-out of precompilation.
   5713.1 ms  ? Turing
   5918.9 ms  ? Turing → TuringOptimExt
WARNING: redefinition of constant OptimizationBase.OPTIMIZER_MISSING_ERROR_MESSAGE. This may fail, cause incorrect answers, or produce other errors.
WARNING: Method definition (::Type{OptimizationBase.OptimizerMissingError})(Any) in module OptimizationBase at /home/runner/.julia/packages/OptimizationBase/sfIfa/src/solve.jl:23 overwritten at /home/runner/.julia/packages/OptimizationBase/sfIfa/src/solve.jl:177.
ERROR: Method overwriting is not permitted during Module precompilation. Use `__precompile__(false)` to opt-out of precompilation.
Precompiling Optimization...
    830.0 ms  ? OptimizationBase
Info Given Optimization was explicitly requested, output will be shown live 
WARNING: redefinition of constant OptimizationBase.OPTIMIZER_MISSING_ERROR_MESSAGE. This may fail, cause incorrect answers, or produce other errors.
WARNING: Method definition (::Type{OptimizationBase.OptimizerMissingError})(Any) in module OptimizationBase at /home/runner/.julia/packages/OptimizationBase/sfIfa/src/solve.jl:23 overwritten at /home/runner/.julia/packages/OptimizationBase/sfIfa/src/solve.jl:177.
ERROR: Method overwriting is not permitted during Module precompilation. Use `__precompile__(false)` to opt-out of precompilation.
   1442.2 ms  ? Optimization
WARNING: redefinition of constant OptimizationBase.OPTIMIZER_MISSING_ERROR_MESSAGE. This may fail, cause incorrect answers, or produce other errors.
WARNING: Method definition (::Type{OptimizationBase.OptimizerMissingError})(Any) in module OptimizationBase at /home/runner/.julia/packages/OptimizationBase/sfIfa/src/solve.jl:23 overwritten at /home/runner/.julia/packages/OptimizationBase/sfIfa/src/solve.jl:177.
ERROR: Method overwriting is not permitted during Module precompilation. Use `__precompile__(false)` to opt-out of precompilation.
Precompiling OptimizationBase...
Info Given OptimizationBase was explicitly requested, output will be shown live 
WARNING: redefinition of constant OptimizationBase.OPTIMIZER_MISSING_ERROR_MESSAGE. This may fail, cause incorrect answers, or produce other errors.
WARNING: Method definition (::Type{OptimizationBase.OptimizerMissingError})(Any) in module OptimizationBase at /home/runner/.julia/packages/OptimizationBase/sfIfa/src/solve.jl:23 overwritten at /home/runner/.julia/packages/OptimizationBase/sfIfa/src/solve.jl:177.
ERROR: Method overwriting is not permitted during Module precompilation. Use `__precompile__(false)` to opt-out of precompilation.
    804.0 ms  ? OptimizationBase
WARNING: redefinition of constant OptimizationBase.OPTIMIZER_MISSING_ERROR_MESSAGE. This may fail, cause incorrect answers, or produce other errors.
WARNING: Method definition (::Type{OptimizationBase.OptimizerMissingError})(Any) in module OptimizationBase at /home/runner/.julia/packages/OptimizationBase/sfIfa/src/solve.jl:23 overwritten at /home/runner/.julia/packages/OptimizationBase/sfIfa/src/solve.jl:177.
ERROR: Method overwriting is not permitted during Module precompilation. Use `__precompile__(false)` to opt-out of precompilation.
WARNING: redefinition of constant OptimizationBase.OPTIMIZER_MISSING_ERROR_MESSAGE. This may fail, cause incorrect answers, or produce other errors.
WARNING: redefinition of constant OptimizationBase.OPTIMIZER_MISSING_ERROR_MESSAGE. This may fail, cause incorrect answers, or produce other errors.
┌ Warning: Replacing docs for `CommonSolve.solve :: Tuple{SciMLBase.OptimizationProblem, Any, Vararg{Any}}` in module `OptimizationBase`
└ @ Base.Docs docs/Docs.jl:243
WARNING: redefinition of constant OptimizationBase.OPTIMIZER_MISSING_ERROR_MESSAGE. This may fail, cause incorrect answers, or produce other errors.
┌ Warning: Replacing docs for `CommonSolve.init :: Tuple{SciMLBase.OptimizationProblem, Any, Vararg{Any}}` in module `OptimizationBase`
└ @ Base.Docs docs/Docs.jl:243
┌ Warning: Replacing docs for `CommonSolve.solve! :: Tuple{SciMLBase.AbstractOptimizationCache}` in module `OptimizationBase`
└ @ Base.Docs docs/Docs.jl:243
Precompiling OptimizationOptimJL...
    837.4 ms  ? OptimizationBase
   1005.2 ms  ? Optimization
Info Given OptimizationOptimJL was explicitly requested, output will be shown live 
┌ Warning: Module Optimization with build ID ffffffff-ffff-ffff-93ac-6c312b9a4a8d is missing from the cache.
│ This may mean Optimization [7f7a1694-90dd-40f0-9382-eb1efda571ba] does not support precompilation but is imported by a module that does.
└ @ Base loading.jl:2541
   1147.7 ms  ? OptimizationOptimJL
┌ Warning: Module Optimization with build ID ffffffff-ffff-ffff-93ac-6c312b9a4a8d is missing from the cache.
│ This may mean Optimization [7f7a1694-90dd-40f0-9382-eb1efda571ba] does not support precompilation but is imported by a module that does.
└ @ Base loading.jl:2541
Precompiling TuringOptimExt...
    814.6 ms  ? OptimizationBase
   1018.6 ms  ? Optimization
   1142.4 ms  ? OptimizationOptimJL
   3863.3 ms  ? Turing
Info Given TuringOptimExt was explicitly requested, output will be shown live 
┌ Warning: Module Turing with build ID ffffffff-ffff-ffff-9066-e5a16071bc0d is missing from the cache.
│ This may mean Turing [fce5fe82-541a-59a6-adf8-730c64b5f9a0] does not support precompilation but is imported by a module that does.
└ @ Base loading.jl:2541
    644.3 ms  ? Turing → TuringOptimExt
┌ Warning: Module Turing with build ID ffffffff-ffff-ffff-9066-e5a16071bc0d is missing from the cache.
│ This may mean Turing [fce5fe82-541a-59a6-adf8-730c64b5f9a0] does not support precompilation but is imported by a module that does.
└ @ Base loading.jl:2541
[ Info: [Turing]: progress logging is disabled globally

Neal (generic function with 2 methods)

In this case, the random variables exposed in the chain (x_raw, y_raw) are not in a helpful form — what we’re after are the deterministically transformed variables x and y.

There are two ways to track these extra quantities in Turing.jl.

Using := (during inference)

The first way is to use the := operator, which behaves exactly like = except that the values of the variables on its left-hand side are automatically added to the chain returned by the sampler. For example:

@model function Neal_coloneq()
    # Raw draws
    y_raw ~ Normal(0, 1)
    x_raw ~ arraydist([Normal(0, 1) for i in 1:9])

    # Transform:
    y := 3 * y_raw
    x := exp.(y ./ 2) .* x_raw
end

sample(Neal_coloneq(), NUTS(), 1000)

┌ Info: Found initial step size
└   ϵ = 1.6

Chains MCMC chain (1000×34×1 Array{Float64, 3}):

Iterations        = 501:1:1500
Number of chains  = 1
Samples per chain = 1000
Wall duration     = 8.19 seconds
Compute duration  = 8.19 seconds
parameters        = y_raw, x_raw[1], x_raw[2], x_raw[3], x_raw[4], x_raw[5], x_raw[6], x_raw[7], x_raw[8], x_raw[9], y, x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9]
internals         = n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size, lp, logprior, loglikelihood

Use `describe(chains)` for summary statistics and quantiles.

Using returned (post-inference)

Alternatively, one can specify the extra quantities as part of the model function’s return statement:

@model function Neal_return()
    # Raw draws
    y_raw ~ Normal(0, 1)
    x_raw ~ arraydist([Normal(0, 1) for i in 1:9])

    # Transform and return as a NamedTuple
    y = 3 * y_raw
    x = exp.(y ./ 2) .* x_raw
    return (x=x, y=y)
end

chain = sample(Neal_return(), NUTS(), 1000)

┌ Info: Found initial step size
└   ϵ = 1.6

Chains MCMC chain (1000×24×1 Array{Float64, 3}):

Iterations        = 501:1:1500
Number of chains  = 1
Samples per chain = 1000
Wall duration     = 1.58 seconds
Compute duration  = 1.58 seconds
parameters        = y_raw, x_raw[1], x_raw[2], x_raw[3], x_raw[4], x_raw[5], x_raw[6], x_raw[7], x_raw[8], x_raw[9]
internals         = n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size, lp, logprior, loglikelihood

Use `describe(chains)` for summary statistics and quantiles.

The sampled chain does not contain x and y, but we can extract the values using the returned function. Calling this function outputs an array:

nts = returned(Neal_return(), chain)

1000×1 Matrix{@NamedTuple{x::Vector{Float64}, y::Float64}}:
 (x = [-0.9761088708359967, -5.204348682724804, 1.4432508265039343, 4.450981614018321, 3.115042983609613, 0.11050499839486816, 0.08596398068289371, 0.3815477315059375, 3.367166729885056], y = 1.703717673028828)
 (x = [0.2303913942136318, -0.25390470540982873, -0.02188270047385421, 0.058660630925560314, 0.5700968837983188, -0.0673035055164295, -0.2102750884269525, -0.5792667540631488, -0.721136046699057], y = -1.1693138929223799)
 (x = [-0.011440632777481325, 0.24534511724383914, 0.04362444018133397, -0.059772664618638074, -0.8091595880119986, 0.2328895681876888, 0.40659182625989726, 0.7736475841459312, 0.6492136629825728], y = -1.146722906157675)
 (x = [0.38816179260737654, 0.46666337707399463, -0.1404805229756647, 0.1504967381913083, -0.3533373863810297, 0.3395836354341499, 0.6300498132681434, 0.19448967484757604, -0.35002888311814795], y = -2.0399443770286436)
 (x = [-2.280544371376041, -3.3622760780356655, 0.3839503044520262, -1.0903555752210654, 2.6817853335229582, -3.5645948713088687, -4.384934562199315, -0.9575978534726586, 2.194664065204745], y = 2.0011886361368387)
 (x = [-1.716170072413519, 3.225451488676409, 3.523941794519998, -1.5769962175693115, 0.9725889947001997, 1.866950809967169, -0.7755537016022253, -2.7651685222809106, -0.5876156147039192], y = 0.8016527285734487)
 (x = [4.951726975785703, -5.074250410368461, -8.353036418359693, 4.184782204304398, -3.410005746264753, 4.752427587292146, -3.147234881753012, 1.6160360741727366, 2.065128743947606], y = 2.3389046346452647)
 (x = [0.17171025163792042, -0.2671852099775187, 0.26312521334586914, -0.19819835653191142, -0.04596335693415075, -0.032746616485977084, -0.12718039812927334, -0.15509933582709978, 0.9143414977282287], y = -2.8633532274072353)
 (x = [0.09515615548952386, -0.053587772348902656, -0.009046222433939424, -0.07574350978133126, -0.041624234219202624, 0.04733008977343868, 0.02752247025310246, -0.06266521695442072, 0.11942061051325717], y = -6.404979426204374)
 (x = [-42.52586622251074, 24.42540557854438, 2.4321232918634452, 29.867706706100492, 10.59007211877523, -18.565333716403437, -8.629106943636614, 8.751427241094486, -54.137294793569986], y = 5.497940661949807)
 ⋮
 (x = [-0.658513825460643, 0.18648828622741537, 0.45553354106317345, -0.6829940850962382, 0.2756310571808607, 1.7132988164784684, 0.2687856861579226, 3.7833540548464484, -0.744198881480148], y = 1.4407582865378221)
 (x = [0.5771418168945011, -0.2759519912628746, -0.0884634325290273, 0.06763678563935, -0.38921931698737194, -0.5006581809068897, 0.2870662123598147, -1.5450104395572921, -0.0353206183710423], y = -1.1486140004817464)
 (x = [-0.04783853889044349, 0.08081243910253795, 0.10735991661352981, -0.03540817669491662, -0.18907057606681746, -0.24150431263091923, 0.31767090372609524, -0.6289179613194281, -0.31201001091935987], y = -2.7730816499171143)
 (x = [-0.08904798418371511, 1.3335440704536925, 1.082274198596173, -0.7188168982527747, 1.1073526618740746, -1.0324345708850426, 1.354590124373184, -0.44446293674611626, 1.0004050325929654], y = -0.1277910324624677)
 (x = [-0.4457758762868777, 0.4652897629082582, 0.3895758610043819, -1.6302850167763125, -0.056660523068051485, -0.8332591120399631, 1.0704440733933138, -1.4399945832021908, 0.07924690294718294], y = -0.6308637671529016)
 (x = [0.07606254617549019, -0.15396983299152167, -0.04502362722747536, -0.2770954782704271, -0.6565418706454927, 0.35789547724181925, -0.03610506858078951, -1.5311368739555167, 0.1523995022759511], y = -1.5875070710627912)
 (x = [0.7955479542761971, 0.7589860873679242, 0.20530415639855876, -0.6549365995051861, 1.342461139600413, -1.1885983801443227, -0.38314485433083584, -1.2203270422536605, -0.12257515412717515], y = 0.0938446652327099)
 (x = [0.29726475537324326, -0.18564281697962415, -0.2560309857915554, 0.33632837421531653, 0.31337422462753717, 0.07360844231236632, -0.27105360432372483, -0.4742416858452971, 0.5185499667926708], y = -1.8042745431453024)
 (x = [-0.09491375005922846, 0.16644096682876838, 0.23169409060925408, -0.2535191688507173, -0.17404304915297975, -0.06951730895210571, -0.006761339943084361, 0.17915551902384122, -0.20379097925151912], y = -2.5773356722575507)

where each element of which is a NamedTuple, as specified in the return statement of the model.

nts[1]

(x = [-0.9761088708359967, -5.204348682724804, 1.4432508265039343, 4.450981614018321, 3.115042983609613, 0.11050499839486816, 0.08596398068289371, 0.3815477315059375, 3.367166729885056], y = 1.703717673028828)

Which to use?

There are some pros and cons of using returned, as opposed to :=.

Firstly, returned is more flexible, as it allows you to track any type of object; := only works with variables that can be inserted into an MCMCChains.Chains object. (Notice that x is a vector, and in the first case where we used :=, reconstructing the vector value of x can also be rather annoying as the chain stores each individual element of x separately.)

A drawback is that naively using returned can lead to unnecessary computation during inference. This is because during the sampling process, the return values are also calculated (since they are part of the model function), but then thrown away. So, if the extra quantities are expensive to compute, this can be a problem.

To avoid this, you will essentially have to create two different models, one for inference and one for post-inference. The simplest way of doing this is to add a parameter to the model argument:

@model function Neal_coloneq_optional(track::Bool)
    # Raw draws
    y_raw ~ Normal(0, 1)
    x_raw ~ arraydist([Normal(0, 1) for i in 1:9])

    if track
        y = 3 * y_raw
        x = exp.(y ./ 2) .* x_raw
        return (x=x, y=y)
    else
        return nothing
    end
end

chain = sample(Neal_coloneq_optional(false), NUTS(), 1000)

┌ Info: Found initial step size
└   ϵ = 1.6

Chains MCMC chain (1000×24×1 Array{Float64, 3}):

Iterations        = 501:1:1500
Number of chains  = 1
Samples per chain = 1000
Wall duration     = 1.65 seconds
Compute duration  = 1.65 seconds
parameters        = y_raw, x_raw[1], x_raw[2], x_raw[3], x_raw[4], x_raw[5], x_raw[6], x_raw[7], x_raw[8], x_raw[9]
internals         = n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size, lp, logprior, loglikelihood

Use `describe(chains)` for summary statistics and quantiles.

The above ensures that x and y are not calculated during inference, but allows us to still use returned to extract them:

returned(Neal_coloneq_optional(true), chain)

1000×1 Matrix{@NamedTuple{x::Vector{Float64}, y::Float64}}:
 (x = [0.6380532198200358, -0.20001903510085747, 0.10544936957548776, 0.6465341921439931, -1.2666990764229111, -0.70500234588546, 0.2810544590554409, -0.5205286463688934, 1.3778375211597105], y = 0.06064093367317169)
 (x = [0.9916709557653861, -2.0728275194689574, -1.1031888717336367, 1.8871737844714132, 0.6523573991039372, -0.21124604969436558, -1.0266370433770664, 0.016367843633954816, -0.9520028661320892], y = 0.8808378286634303)
 (x = [0.9916709557653861, -2.0728275194689574, -1.1031888717336367, 1.8871737844714132, 0.6523573991039372, -0.21124604969436558, -1.0266370433770664, 0.016367843633954816, -0.9520028661320892], y = 0.8808378286634303)
 (x = [0.9916709557653861, -2.0728275194689574, -1.1031888717336367, 1.8871737844714132, 0.6523573991039372, -0.21124604969436558, -1.0266370433770664, 0.016367843633954816, -0.9520028661320892], y = 0.8808378286634303)
 (x = [-0.14654345646913047, 0.02988128864425104, 0.057369516730195164, 0.11313802617462844, -0.008816097312658143, 0.15938754777753752, 0.0004759963443701986, 0.12685932236706945, -0.1882154922201721], y = -4.003330270538967)
 (x = [-1.9112261161877502, -1.0501328665382703, -1.5147343746267197, 1.5501953259001175, 0.31502420272015474, -0.27688168746956393, -0.5775016267306105, 1.5780017910818112, -0.5059757921097513], y = -0.14607204656749184)
 (x = [-0.14008265057089722, -0.1870769772786126, 0.1787393768978083, -0.08435655678800534, -0.06618009866024141, -0.012560611594410465, 0.09604795723751479, 0.08839032716413815, 0.04983034199015561], y = -4.466868084894504)
 (x = [-0.04360199331402114, -0.08590225688749864, 0.017737118973307337, -0.021349240640300342, -0.009883763485161021, -0.02631945812834671, 0.104193768913333, 0.023243873075089322, 0.0017859952391018152], y = -6.7941560883655425)
 (x = [-11.600432163736196, -16.711168624526607, 20.121761087003026, -4.526158406524454, -17.224166308807494, 33.51184352755993, -63.92154268607287, 24.645633571808595, -0.2771454343287743], y = 5.9744566458499975)
 (x = [-0.07860711530692682, -0.8797116123316353, 2.2245231255440476, 0.09971209590884506, -1.2176222831727987, 0.22495328065276052, -0.4393731670862836, 1.3173166868347441, 0.9571430218682557], y = -0.5738517104174914)
 ⋮
 (x = [2.5565322922144254, 1.0042552013810164, 1.1131742175575936, 2.283210171973585, 0.25257222661031314, 0.8175175974751011, -0.41753732004234617, -1.5707151967275355, -0.2489048385532152], y = 0.3123042577399735)
 (x = [-0.8729960829111867, 0.3063000470056581, 0.3534239552381912, -0.254011276048432, -0.3648047691881745, -0.4521346127336582, -0.21467950196092594, 0.2473166829462221, -0.7015477642212055], y = -1.3676739223033147)
 (x = [1.869529439323088, 0.12302693236713168, 0.6267622385806317, 0.551612493568384, -1.2959538395346135, -0.8424337695017085, 2.466580889451787, 3.0584099651256964, -1.7862789781422144], y = 1.685599527542134)
 (x = [-0.5754507164140175, -0.30300316985307896, -0.8245605378713297, 1.1995739149458753, -0.5917707780671618, 0.6192795933747774, -0.6338943024865548, 0.324066673713136, 0.3301890858739846], y = 0.441028313405728)
 (x = [0.0220312788987611, 0.06818716535730764, -0.005477977554137843, -0.22584554963584105, 0.032413663150197676, 0.015189059667802514, -0.03086990090079641, -0.2354487613957679, -0.24632324914493986], y = -2.9951376953070374)
 (x = [0.5086623616413938, 0.25727410444069065, 0.0734628839788139, -0.2505155626464654, -0.09117779844709314, 0.46605224461499, -0.07701899312364045, 0.27825584090389344, 0.11461896293724595], y = -2.5412058405605884)
 (x = [-2.6807116540895777, -2.3580932623133672, 1.7556805711864372, 1.1759932073044164, 0.3736563299446332, -2.261835692338751, 0.5064599310627613, -1.8407745391710881, -0.8932276306470595], y = 1.0409739607623134)
 (x = [1.9491345682844738, -0.3928061208870163, -1.7566428998236463, -0.5016175777033682, -0.5097953117571379, 1.9252693980403146, -0.4298463412264657, 1.3940725839537542, 0.8765806400015815], y = 0.7010118467645461)
 (x = [-1.0479137175416477, 0.02609327627528047, -0.2546327425389955, 0.6201820320429676, 0.2624754223757027, -1.0562807644831405, 0.29553514504322986, -1.4269115230495113, -0.210642080094762], y = -0.08915617007831611)

Another equivalent option is to use a submodel:

@model function Neal()
    y_raw ~ Normal(0, 1)
    x_raw ~ arraydist([Normal(0, 1) for i in 1:9])
    return (x_raw=x_raw, y_raw=y_raw)
end

chain = sample(Neal(), NUTS(), 1000)

@model function Neal_with_extras()
    neal ~ to_submodel(Neal(), false)
    y = 3 * neal.y_raw
    x = exp.(y ./ 2) .* neal.x_raw
    return (x=x, y=y)
end

returned(Neal_with_extras(), chain)

┌ Info: Found initial step size
└   ϵ = 1.6

1000×1 Matrix{@NamedTuple{x::Vector{Float64}, y::Float64}}:
 (x = [7.758972772076274, 8.39101400750792, -2.1843300677542765, -7.124105158879408, -1.310960605391347, -37.12158517601075, -6.250630041420653, 6.364497329411868, 15.489323487891296], y = 5.352134036500144)
 (x = [0.07015531987619926, -0.016846181350943075, 0.054969421117734396, 0.0560856596663242, -0.010014206363082282, 0.17962330252239342, 0.02439570400701714, -0.031135823573792368, -0.0709847021852034], y = -5.348537685608913)
 (x = [-0.24679515283668452, -0.7625842932375881, 1.7348760568488484, 0.2735120979623507, 0.27065868419288613, 0.8932697404338782, -0.3597728550025933, -0.9543678589792086, 1.193909561768782], y = 0.6282289453250869)
 (x = [0.10866219194545504, -1.126382010635536, 0.48828231098492986, 0.11085594150516442, -1.6365886839223487, 0.29952184732077386, -0.8966444521523798, -1.4478962999715463, 1.0300772442308967], y = 0.6003661237934048)
 (x = [-0.0851284470648633, 0.05626309618294626, 0.021310879799102716, 0.0363046412541881, 1.2219752236820987, -0.403047900073925, 0.21839915160777695, 0.7856795308370973, 0.15545586235765643], y = -0.6099640349329358)
 (x = [0.6256445589766079, -1.7085159767719875, -0.19649227468550715, -0.47747899628294427, 1.3784664696711741, -1.286107338676066, 1.412840126718721, -1.4995803630856896, -0.4811574663028566], y = -0.08322796569921487)
 (x = [-0.072336920804889, 1.7436681027389926, 0.5933442475813333, 0.674009341244435, -1.1399078750180842, 0.7901605994988302, -1.1135550547543163, 0.3065339134155228, 0.48849642383206276], y = -0.45932070096337574)
 (x = [-0.8933981489396019, 1.8442085660454073, 0.30962676452449106, 0.22382720363003036, 0.7538398841200197, 0.4037216687430575, -0.21704835158297933, 0.46458058082877834, 1.4648335668709234], y = -0.5098086675573463)
 (x = [-1.3927440408286693, 1.802130826628073, -1.4708655922662885, -0.21398167124729167, 0.002329617315958729, 0.14391114385960993, 0.1492225610126328, 1.3177499597168223, 0.5297838830102867], y = 0.3454026145396995)
 (x = [0.4280037315253775, -0.3479113607418051, -1.3298152275839343, -0.6004002711533453, 0.6545926508383214, -0.5081667196753212, 0.4883709235561456, -0.408354260745982, -0.5728340816052414], y = -0.8027141739589381)
 ⋮
 (x = [3.8998158455789453, -8.764290676089578, -7.310044011530945, -6.206942927028825, -2.2174244076160945, 1.1682507406201346, 6.801646764611466, 1.4733620425139795, -2.1573776931241566], y = 2.6408922088272524)
 (x = [-0.24258979404466735, 0.5527245285592063, 0.6364548800288444, 0.3815198277582032, -0.029040639911697145, -0.14470389207203127, -0.2888797914371354, -0.08973315011596933, 0.25410614132973824], y = -2.656381402174953)
 (x = [0.40672520931479617, -8.702823565791903, -8.802954912683365, -4.089370898632748, -0.9943674474606897, 3.7107773567021285, 1.2434384884220855, 1.4940006734882925, -0.44812117585803274], y = 2.708713906296444)
 (x = [1.519683305270694, 2.298802178988467, -1.4859064046757948, 0.6028788177716148, 1.5734589449727139, -2.9893092639773995, -1.0179210422227523, 1.0279675487489874, -0.9681558913668762], y = 1.6634559370487105)
 (x = [-0.9998465703853113, 1.2433106841779027, -0.6893899090599035, 0.4627641736677376, -0.14846810402177266, -0.2746582147701319, 0.10986320777874371, 0.10259317850116467, 0.47381913118690827], y = -0.5128606770801295)
 (x = [0.998119530875437, -2.906568435983433, 0.7282911906584546, -0.6710806737125761, -0.32867031558246806, 0.44798736352269386, 1.256198373557169, -0.2790012049397262, 0.06250997724931638], y = 0.6121349930816002)
 (x = [-18.282790465430754, 26.584549796206346, -4.256658173303513, 8.031649718077828, 4.749175684946737, 1.625554791802499, -40.02888938569024, 0.2985953212189829, -7.86050481141781], y = 5.749280362092253)
 (x = [-0.9676580739999089, -0.03356510361435681, -0.6195407749282014, 0.2499487442963725, 1.3843229235240526, -1.4643889264210967, -0.30019796545708105, 0.37307387687312427, -0.37390057290126505], y = -0.46457872895205443)
 (x = [-0.950081070723225, 18.400237426492954, 16.95324374584172, 3.7713649016918613, 12.155887858918208, -8.071960799104202, -3.944292663539916, 18.136577136869104, -6.23743594385871], y = 5.124523664944405)

Note that for the returned call to work, the Neal_with_extras() model must have the same variable names as stored in chain. This means the submodel Neal() must not be prefixed, i.e. to_submodel() must be passed a second parameter false.