6  Putting it together

#TODO: Finish page.

6.1 Supplementary results

6.1.1 Varying fixed-lag resampling \(L\)

Fixed-lag resampling introduces a potential source of bias. Fundamentally, we are approximating the smoothing posterior \(P(X_t|y_{1:T})\) with \(P(X_t|y_{1:t+L})\), ignoring any direct association between the current hidden state \(X_t\) and data collected more than \(L\) steps in the future. This can be an issue if the process does not “forget” fast enough (for example, if \(R_t\) varies very slowly over time), or if \(L\) is smaller than the maximum¹ generation time.

Compared to hidden-state estimation, the choice of \(L\) is generally less important when estimating fixed parameters using PMMH, as only the filtering distribution \(P(X_t|y_{1:t})\) is required for this task. In the Markovian setting, all past-state resampling could be disabled entirely, but as the renewal model is inherently non-Markovian, past-state resampling is required. However, particle degeneracy is not a concern, so \(L\) could be set to \(t\) (full past-state resampling), but this is computationally expensive. Using fixed-lag resampling thus strikes a balance between these two extremes, provided \(L\) is large enough that hidden-state trajectories are consistent over the maximum generation time.

To test the effect of varying \(L\), we fit the third example model from Section 8.3, with \(L = 5, 10, 20, 40, 80\). For context, the mean generation time is assumed to be 6.5 days, with \(L = 20\) covering 98.8% of the probability mass of the generation time distribution. The similarity of posterior means and 95% credible intervals of the parameters \(\sigma\) and \(\phi\) suggest there is little systematic impact of \(L\) (the variation observed arises from Monte Carlo error) on parameter estimation in this example. Greater variation is seen in estimates of the reproduction number at low \(t\), where there is limited data. Despite this, there is a high level of agreement between estimates for later estimates, despite the known problems with using \(L < 20\).

using Plots; gr()
using Measures, CSV, DataFrames


include("../src/LoadData.jl")
include("../src/PMMH.jl")
include("../src/FitModel.jl")
include("../src/MarginalPosterior.jl")
include("../src/Support.jl")

Y = loadData("NZCOVID")

function ExampleModel3(θ, Y::DataFrame, opts::Dict)
    
    # Extract frequently used options
    T = opts["T"] # Number of time steps
    N = opts["N"] # Number of particles to use
    L = opts["L"] # Length of fixed-lag resampling
    
    # Define the serial interval
    ω = pdf.(Gamma(2.36, 2.74), 1:100) # (Unnormalised) serial interval
    ω = ω/sum(ω) # Normalise the serial interval
    
    # Initialise output matrices
    R = zeros(N, T) # Matrix to store particle values
    I = zeros(N, T) # Local infections
    W = zeros(N, T) # Matrix to store model weights
    
    # Sample from initial distribution
    R[:,1] = rand(Uniform(0, 10), N)
    I[:,1] = rand.(opts["pI0"], N)
    
    # Run the filter
    for tt = 2:T
        
        # Project according to the state-space model
        R[:,tt] = exp.(rand.(Normal.(log.(R[:,tt-1]), θ[1])))
        
        # Calculate force-of-infection and sample infections
        Λ = sum(I[:,tt-1:-1:1]  .* opts["ω"][1:tt-1]', dims=2)
        I[:,tt] = rand.(Poisson.(R[:,tt] .* Λ))
        
        # Weight according to the observation model
        r = 1/θ[2]
        p = 1 ./ (1 .+ θ[2] * I[:,tt])
        W[:,tt] = pdf.(NegativeBinomial.(r, p), Y.local[tt])
        
        # Resample
        inds = wsample(1:N, W[:,tt], N; replace=true)
        R[:, max(tt - L, 1):tt] = R[inds, max(tt - L, 1):tt]
        I[:, max(tt - L, 1):tt] = I[inds, max(tt - L, 1):tt]
        
    end
    
    # Store output as three-dimensional array
    X = zeros(N, T, 2)
    X[:,:,1] = R
    X[:,:,2] = I
    
    return(X, W)
    
end


opts = Dict(

# Bootstrap filter options
"T" => size(Y, 1), # Number of time-steps
"N" => 1000, # Number of particles
"L" => 50, # Fixed-lag resampling length
"ω" => pdf.(Gamma(2.36, 2.74), 1:100), # Serial interval
"pR0" => Uniform(0, 10), # Prior on Rt at t = 0
"pI0" => Dirac(1),
"predictiveValues" => false, # Whether to calculate predictive cases

# PMMH options
"nChains" => 3, # Number of chains
"chunkSize" => 100, # Number of iterations per chunk
"maxChunks" => 50, # Maximum number of chunks
"maxRhat" => 1.02,  # Stopping criterion: maximum Rhat value
"minESS" => 200, # Stopping criterion: minimum effective sample size
"showChunkProgress" => true, # Whether to show progress of each chunk
"propStdDevInit" => sqrt.([0.1, 0.01]), # Initial proposal standard deviation (this is adaptively fit)
"paramPriors" => [Uniform(0, 1), Uniform(0, 1)],
"initialParamSamplers" => [Uniform(0.1, 0.3), Uniform(0.01, 0.03)],
"paramLimits" => [(0, 1), (0, 1)],
"paramNames" => ["σ", "ϕ"],

# Marginal posterior options
"posteriorNumberOfParticles" => 10000, # Number of particles to use for marginal posterior
"posteriorParamSamples" => 100,
"stateNames" => ["Rt", "It"]
)



function runVaryingL()
    
    Lvals = [5, 10, 20, 40, 80]
    
    df_states_all = DataFrame()
    df_params_all = DataFrame()
    
    for L in Lvals
        
        println("L = ", L)
        
        # Update the options with the new L value
        opts_in = deepcopy(opts)
        opts_in["L"] = L
        
        # Run the PMMH algorithm
        (df_states, df_params, _, _) = fitModel(ExampleModel3, Y, opts_in; verbose=true)
        df_states[!,:L] .= L
        df_params[!,:L] .= L
        
        # Store
        df_states_all = vcat(df_states_all, df_states)
        df_params_all = vcat(df_params_all, df_params)
        
    end
    
    # Return the results
    return(df_states_all, df_params_all)
    
end

# # Run all models (this takes a few minutes) and save the results
# (df_states_all, df_params_all) = runVaryingL()
# CSV.write("results/df_states_all.csv", df_states_all)
# CSV.write("results/df_params_all.csv", df_params_all)


function plotVaryingL_animated()
    df_states_all = CSV.read("results/df_states_all.csv", DataFrame)
    df_params_all = CSV.read("results/df_params_all.csv", DataFrame)
    Lvals = sort(unique(df_states_all.L))

    dfsig = df_params_all[df_params_all.param .== "σ", :]
    plot_sigma = plot(xlab="L", ylab="σ", xscale=:log2, xticks=(dfsig.L, string.(dfsig.L)), label=false, ylims=(0, 0.5), xlims=(4, 110))
    plot_sigma = scatter!(dfsig.L, dfsig.m, yerror=(dfsig.m-dfsig.l, dfsig.u-dfsig.m), label=false)

    dfphi = df_params_all[df_params_all.param .== "ϕ", :]
    plot_phi = plot(xlab="L", ylab="ϕ", xscale=:log2, xticks=(dfphi.L, string.(dfphi.L)), label=false, ylims=(0, 0.1), xlims=(4, 110))
    plot_phi = scatter!(dfphi.L, dfphi.m, yerror=(dfphi.m-dfphi.l, dfphi.u-dfphi.m), label=false)

    anim = Animation()
for i in 1:length(Lvals)
    plot_R = plot(xlab="Date", ylab="Reproduction number", label=false, legend=:topright, ylim=(0, 10), color_palette=:Dark2_8)
    for (j, L) in enumerate(Lvals[1:i])
        dfL = df_states_all[df_states_all.L .== L, :]
        dfLR = dfL[dfL.variable .== "Rt", :]
        plot_R = plot!(dfLR.date, mean.(dfLR.mean), ribbon=(dfLR.mean-dfLR.lower, dfLR.upper-dfLR.mean), label="L = $L", color=palette(:Dark2_8)[j])
    end
    p = plot(plot(plot_sigma, plot_phi, layout=(1,2)), plot_R, layout=(2,1), size=(600,400), margins=3mm)
    frame(anim, p)
end
gif(anim, "varying_L.gif", fps=1, show_msg=false)

end

plotVaryingL_animated()

---

1. Many generation time distributions are assumed to have support on \((0, \infty)\). In this case, \(L\) should be chosen such that it is greater than most of the possible generation times, such as the \(99^{th}\) percentile, for example.