10.2. Ajustando um modelo SIR a uma epidemia de Influenza

using DifferentialEquations
using Random
using Distributions: Normal, MvNormal
using Statistics
using LsqFit
using Optim
using Plots
using Images

Vamos, agora, ajustar os parâmetros \(\beta\) e \(\gamma\) do modelo SIR aos dados de uma epidemia de Influenza em um escola de meninos (internato) no norte da Inglaterra.
Em um total de 763 residentes, em um período de duas semanas, o número de infectados a cada dia é dado pela lista

\[ \text{Infectados a partir do dia três} = [25, 75, 228, 297, 259, 235, 192, 126, 71, 28, 9, 7]. \]

Fonte: Brauer & Castillo-Chavez (2012).
Há uma lista um pouco diferente de infectados, que inclui desde, o dia da primeira ocorrência:

\[ \text{Infectados a partir do dia zero} = [1, 3, 7, 25, 72, 222, 282, 256, 233, 189, 123, 70, 25, 11, 4]. \]

Fonte: Ledder(2013)
Vamos usar, abaixo, os primeiros dados, de Brauer & Castillo-Chavez (2012).

Ajuste

Nesse caso, o modelo nos dá, a partir dos dados iniciais e dos parâmetros, uma solução \((S(t), I(t), R(t))\).
Mas só temos dados dos números de infectados \(I_j\), a cada dia \(j=0, \ldots, 14\).
Assim, a nossa função objetivo só pode ser construída em cima de \(I(t)\) e dos dados \((t_j, I_j)\).
Em algumas situações, as próprias condições iniciais também podem entrar como parâmetros a serem ajustados, mas no caso mais controlado acima, vamos fixar as seguinte condições iniciais:
- \(I(0) = I_j = 1\) residente;
- \(S(0) = 763 - 1 = 762\) residentes; e
- \(R(0) = 0\) residentes.
Definimos a taxa de evolução do sistema \((S,I)\) do modelo SIR:

# Data from Brauer & Castillo-Chavez (2012)
N = 763
dias_amostra = [3:14...]
infectados_amostra = [25, 75, 228, 297, 259, 235, 192, 126, 71, 28, 9, 7]

# Data from Ledder
#dias_amostra = [1:15...]
#infectados_amostra = [1, 3, 7, 25, 72, 222, 282, 256, 233, 189, 123, 70, 25, 11, 4]

nothing

plot(dias_amostra, infectados_amostra, linestyle=:dash)
scatter!(dias_amostra, infectados_amostra, color=1)
plot!(xlabel="dias", ylabel="infectados",
    title="Número de infectados a cada dia na escola no Reino Unido", titlefont=10)

Chute inicial dos parâmetros

Observe que a epidemia cresce e decresce rapidamente, sugerindo uma alta taxa de transmissão e um curto tempo de recuperação.
Por conta disso, vamos chutar, inicialmente, os seguintes valores dos parâmetros:
- \(\beta = 1.5\)
- \(\gamma = 0.5\) (dois dias de recuperação)
- \(\mathcal{R} = \beta/\gamma = 3.0\)

ρ₀ = [1.5, 0.5]
u₀ = [N-1, 1]
tspan = (0.0, 14.0)

(0.0, 14.0)

function dudt_SIR!(du, u, p, t)
    β, γ = p
    S, I = u
    R = N - S - I
    I_nov = β * I * S / N # novos infectados
    du[1] = - I_nov
    du[2] = I_nov - γ * I
end

dudt_SIR! (generic function with 1 method)

model(t, ρ) = solve(ODEProblem(dudt_SIR!, u₀, tspan, ρ), Tsit5())(t)[2]

model (generic function with 1 method)

plot(dias_amostra, infectados_amostra, linestyle=:dash, label=false)
scatter!(dias_amostra, infectados_amostra, color=1, label="Dados")
tempos = 1.0:0.25:14.0
plot!(tempos, model.(tempos, Ref(ρ₀)), label="Modelo β=$(ρ₀[1]), γ=$(ρ₀[2])")
plot!(xlabel="dias", ylabel="infectados", legend=:topright,
    title="Número de infectados a cada dia na escola no Reino Unido", titlefont=10)

Melhorando o chute inicial

Podemos ajustar um pouco mais aumentando a taxa de transmissão:
- \(\beta = 1.8\) (taxa de contágio)
- \(\gamma = 0.5\) (dois dias de recuperação)
- \(\mathcal{R} = \beta/\gamma = 3.6\)

ρ₀ = [1.8, 0.5]

2-element Vector{Float64}:
 1.8
 0.5

plot(dias_amostra, infectados_amostra, linestyle=:dash, label=false)
scatter!(dias_amostra, infectados_amostra, color=1, label="dados")
tempos = 1.0:0.25:14.0
plot!(tempos, model.(tempos, Ref(ρ₀)), label="Modelo β=$(ρ₀[1]), γ=$(ρ₀[2])")
plot!(xlabel="dias", ylabel="infectados", legend=:topright,
    title="Número de infectados a cada dia na escola no Reino Unido", titlefont=10)

Assim parece mais próximo, mas vamos fazer o ajuste usando LsqFit, com esse novo chute inicial.

Ajuste com LsqFit

fit = curve_fit((t,ρ) -> model.(t, Ref(ρ)), dias_amostra, infectados_amostra, ρ₀)

LsqFit.LsqFitResult{Vector{Float64}, Vector{Float64}, Matrix{Float64}, Vect
or{Int64}}([1.6682491395946721, 0.44168556986258944], [11.300158076906051, 
26.713035784328596, -16.76828047337159, -9.563924212963059, 18.957607286142
775, -11.591611941519062, -27.387005722741605, -9.670669111835466, 9.454353
054971762, 26.983336226454668, 28.314904151597133, 18.215012626217053], [10
0.56002432271706 -104.13885253850472; 313.99606007189016 -345.1498581600921
5; … ; -75.41290593639617 -143.58093947773298; -53.65106603063737 -111.8684
258497654], true, Int64[])

ρ = fit.param

2-element Vector{Float64}:
 1.6682491395946721
 0.44168556986258944

plot(dias_amostra, infectados_amostra, linestyle=:dash, label=false)
scatter!(dias_amostra, infectados_amostra, color=1, label="Dados")
tempos = 1.0:0.25:14.0
plot!(tempos, model.(tempos, Ref(ρ)), label="Modelo ajustado")
plot!(xlabel="dias", ylabel="infectados", legend=:topright,
    title="Número de infectados a cada dia na escola no Reino Unido", titlefont=10)

covar = estimate_covar(fit)

2×2 Matrix{Float64}:
 0.000861757  0.000194578
 0.000194578  0.000313972

Outro chute inicial

A escolha de um bom chute inicial é importante, como podemos ver na tentativa de ajuste abaixo, com outro chute inicial.
Acabou achando os mesmos parâmetros, mas com mais dificuldade.
Em outros casos, pode nem funcionar.

curve_fit((t,ρ) -> model.(t, Ref(ρ)), dias_amostra, infectados_amostra, [0.5, 0.2]).param

2-element Vector{Float64}:
 1.6682492617500464
 0.4416855505070483

Ajuste com o Optim

Também podemos usar o Optim.
Ele é mais flexível (mais liberdade na escolha e na construção da função custo).
E costuma funcionar melhor.
Ele só não nos dá, automaticamente, a matrix de covariância.

custo(ρ) = sum(abs2, infectados_amostra - model.(dias_amostra, Ref(ρ)))

custo (generic function with 1 method)

ρ₀ = [1.5, 0.5]
result = optimize(custo, ρ₀)

* Status: success

 * Candidate solution
    Final objective value:     4.502256e+03

 * Found with
    Algorithm:     Nelder-Mead

 * Convergence measures
    √(Σ(yᵢ-ȳ)²)/n ≤ 1.0e-08

 * Work counters
    Seconds run:   0  (vs limit Inf)
    Iterations:    68
    f(x) calls:    145

Optim.minimizer(result)

2-element Vector{Float64}:
 1.6682492852905042
 0.44168581226548664

Propagação de incertezas

Usando a matrix de covariância obtida com o LsqFit.
E usando o método de Monte-Carlo.
Podemos calcular os intervalos de confiança ao longo do tempo.

num_amostras = 200
parameters = rand(MersenneTwister(21001), MvNormal(ρ, covar), num_amostras)

2×200 Matrix{Float64}:
 1.65479  1.64094   1.64311   1.64178   …  1.65732   1.6368    1.64749
 0.46985  0.435592  0.429651  0.403384     0.418683  0.469833  0.456223

scatter(parameters[1,:], parameters[2,:])

Incerteza

plot(dias_amostra, infectados_amostra, linestyle=:dash, label=false)
scatter!(dias_amostra, infectados_amostra, color=1, label="Dados")
tempos = 1.0:0.25:14.0
simulations = fill(0.0, length(tempos), num_amostras)
for n in 1:num_amostras
    ρ̃ = parameters[:,n]
    simulations[:,n] = model.(tempos, Ref(ρ̃))
    plot!(tempos, simulations[:,n], label=false, alpha=0.1, color=2)
end
plot!(tempos, mean(simulations, dims=2), color=6, label="média simulações")
plot!(xlabel="dias", ylabel="infectados", legend=:topright,
    title="Número de infectados a cada dia na escola no Reino Unido", titlefont=10)

Exercícios

Considere o modelo SEIR e ajuste os seus parâmetros aos dados da escola do Reino Unido considerados no texto.
Idem para o modelo SIAR.
Idem para o modelo SEIAR.

(CC BY-NC-ND 4.0) Attribution-NonCommercial-NoDerivatives 4.0 International

Last modified: July 10, 2022. Built with Franklin.jl.