Epidemiological Models for  COVID-19. Application to Spain prediction

Guillermo Sánchez, Juan M. Rodríguez Díaz, Jesús L. Fidalgo .

http://diarium.usal.es/guillermo/mathematica/modelizacion-covid-19/

Updated :

DateObject[]

Introducción

Este documento sirve de base para la realización de predicciones realizadas por el DOE USAL como contribución a http://matematicas.uclm.es/cemat/covid19/

SetDirectory@NotebookDirectory[];

Nosotros utilizamos Mathematica

Recursos sobre COVID-29 con Mathematica: https://community.wolfram.com/groups/-/m/t/1872608

COVID en el Mundo

The Johns Hopkins University data are available from the Wolfram Data Repository, and are described at https://www.wolframcloud.com/obj/resourcesystem/published/DataRepository/resources/Epidemic-Data-for-Novel-Coronavirus-COVID-19

covid19Data=ResourceData["Epidemic Data for Novel Coronavirus COVID-19"];

The function fitDataWDRprovides convenient access to the data for modeling

?fitDataWDR

ListLogPlot[fitDataWDR[LinguisticAssistant,"1 March 2020"],PlotLabel->"Spain",FrameLabel->{"time (d)","# individuals"},PlotLegends->{"confirmed"-("recovered"+"dead"),"recovered","dead"}]

ListLogPlot[fitDataWDR[LinguisticAssistant,"1 March 2020"],PlotLabel->"Italy",FrameLabel->{"time (d)","# individuals"},PlotLegends->{"confirmed"-("recovered"+"dead"),"recovered","dead"}]

ListLogPlot[fitDataWDR[Entity["Country", "France"],"1 March 2020"],PlotLabel->"France",FrameLabel->{"time (d)","# individuals"},PlotLegends->{"confirmed"-("recovered"+"dead"),"recovered","dead"}]

ListLogPlot[fitDataWDR[Entity["Country", "Germany"],"1 March 2020"],PlotLabel->"Germany",FrameLabel->{"time (d)","# individuals"},PlotLegends->{"confirmed"-("recovered"+"dead"),"recovered","dead"}]

Importación de datos de España por comunidades

Castilla y León descarga desde la JCyL

La información de CyL mas detallada se pueden encontrar en : https://analisis.datosabiertos.jcyl.es/explore/?sort=modified&refine.keyword=COVID-19

Importamos los datos en csv pero usan como separador “;” es necesario convertirlo en “, “

Situación hospitalaria

datosHospitalesCL=SemanticImport["https://analisis.datosabiertos.jcyl.es/explore/dataset/situacion-de-hospitalizados-por-coronavirus-en-castilla-y-leon/download/?format=csv&timezone=Europe/Madrid&lang=es&use_labels_for_header=true&csv_separator=,"]

Sumamos los totales por día y mostramos la información que nos interesa

datosHospitalesCL1=GroupBy[datosHospitalesCL,First,Total];

datosHospitalesCL2=datosHospitalesCL1[All,4;;7]

{hosplanta,hospuci, hospfall,  hospaltas} =Map[datosHospitalesCL2[All, #]&,{ "hospitalizados_planta","hospitalizados_uci","fallecimientos","altas"}];

DateListStepPlot[{hosplanta,hospuci, hospfall,  hospaltas},PlotLegends->{"hospplanta","uci","fallecidos","altas"},ScalingFunctions->"Log"]

Otra forma alternativa

hopistalizadosCyL = Import["https://analisis.datosabiertos.jcyl.es/explore/dataset/situacion-de-hospitalizados-por-coronavirus-en-castilla-y-leon/download/?format=csv&timezone=Europe/Madrid&lang=es&use_labels_for_header=true&csv_separator=%3B", "Table", "FieldSeparators" -> ";", "CurrencyTokens" -> None];

hopistalizadosCyL=Import["https://analisis.datosabiertos.jcyl.es/explore/dataset/situacion-de-hospitalizados-por-coronavirus-en-castilla-y-leon/download/?format=csv&timezone=Europe/Madrid&lang=es&use_labels_for_header=true&csv_separator=,", "csv"];

Nos interesa la siguiente información

cabeceraCL=First[hopistalizadosCyL][[4;;7]]

dataCL=Normal[Take[GroupBy[Rest[hopistalizadosCyL], First, Total],All,{ 4, 7}]];

fechasCL1=dataCL[[All,1]];

datosCL1A = fechasCL1/.dataCL;

{hospplantaCL, hospuuciCL, altasCL, fallCL}=Transpose[datosCL1A ];

TableForm[datosCL1A ,TableHeadings->{fechasCL1,cabeceraCL}]

Situación epidemiológica por provincias

datosCL2=Import["https://analisis.datosabiertos.jcyl.es/explore/dataset/situacion-epidemiologica-coronavirus-en-castilla-y-leon/download/?format=csv&timezone=Europe/Madrid&lang=es&use_labels_for_header=true&csv_separator=,"];

cabeceraCL2= First[datosCL2][[3;;6]]

datosCL2A=Drop[Normal[Take[GroupBy[ Rest[datosCL2], First, Total],All,{ 3, 6}]],5];

fechasCL2=datosCL2A[[All,1]];

datosCL2B = fechasCL2/.datosCL2A;

TableForm[datosCL2B ,TableHeadings->{fechasCL2,cabeceraCL2}]

{casosCL, nuevosCL, altasCL1, fallCL1}=Transpose[datosCL2B ];

Se presentan todos los datos de CL, es necesario empezar el 18 de marzo pues es el primer día del que se dispone datos para todos los campos,  los campos altas y fallecimientos se muestran por duplicado aunque tienen distintos origines, se ve que los datos coinciden en consecuencia los fallecimientos son todos hospitalarios lo mismo que las altas. En el hospitalizados están sumados los de planta y los de UCI

TableForm[Transpose[{casosCL, nuevosCL, hospplantaCL+hospuuciCL, hospuuciCL,altasCL,  altasCL1,fallCL,  fallCL1}],TableHeadings->{fechasCL2,{"casos","nuevos","Hospitalizados","UCI","altas\n hospitalarias","altas","fallecimientos\n hospitalarios","fallecimientos"}}]

Datos para algunas comunidades, descargados desde el Instituto Carlos III

Descarga

Importamos los datos (desde el Instituto Carlos III y eliminamos la última fila pues corresponde a comentarios

datoscomunidades=Take[Drop[Import["https://covid19.isciii.es/resources/serie_historica_acumulados.csv"],-2], All, 7];

El fichero descargado contiene las siguientes notas:
El objetivo de los datos que se publican en esta web es saber el número de casos acumulados a la fecha y que por tanto no se puede deducir que la diferencia entre un día y el anterior es el número de casos nuevos ya que esos casos pueden haber sido recuperados de fechas anteriores. Cualquier inferencia que se haga sobre las diferencias de un día para otro deben hacerse con precaución y son únicamente la responsabilidad el autor. Los datos de estas comunidades son datos de prevalencia (personas ingresadas a fecha de hoy). No reflejan el total de personas que han sido hospitalizadas o ingresadas en UCI  a lo largo del periodo de notificación(CL-CM-MD-VC-MC)

El Ministerio de Sanidad, en sus informes diarios, distingue: “Casos que han precisado hospitalización” y “Casos que han ingresado en UCI” e incluye una nota en la que indica que (referido a  Comunidades de Galicia, Castilla-La Mancha, Castilla y León, Comunidad Valenciana y la Comunidad de Madrid) “Los datos de estas comunidades son datos de prevalencia (personas ingresadas a fecha de hoy). No reflejan el total de personas que han sido hospitalizadas o ingresadas en UCI a lo largo del periodo de notificación”.

Extraemos la primera fila (o sublista) donde se muestran los nombres de los campos, que la usaremos mas adelante

cabecera=First[datoscomunidades]

Nos quedamos con los datos a partir del día 15 de marzo y agrupamos por CC AA

agrupadosCCAA=Normal[GroupBy[Drop[datoscomunidades,457], First]];

idemccaa=agrupadosCCAA[[All,1]]

datosCCAA=idemccaa/.agrupadosCCAA;

Extraemos los datos de cada CCAA

{datosAN,datosAR,datosAS,datosIB,datosCN,datosCB,datosCM,datosCL,datosCT,datosCE,datosVC,datosEX,datosGA,datosMD,datosME,datosMC,datosNC,datosPV,datosRI}=datosCCAA;

{ca, fecha, casosCM, hospCM, uciCM, falCM, recCM}=Transpose[datosMD]

{ca, fechaCL, casosCL, hospCL, uciCL, falCL, recCL}=Transpose[datosCL];

{ca, fecha, casosAN, hospAN, uciAN, falAN, recAN}=Transpose[datosAN];

rango=Range[15,Length[casosCM]+15-1]

Convertimos la fecha a formato ISO

date=Map[DateList[{#,{"Day","Month","Year"}}]&, fecha]

Madrid

{infectadosCM, hospitalizadosCM, muertosCM, recuperadosCM}=Map[Transpose[{rango,#}]&,{casosCM, hospCM, falCM, recCM}];

TableForm[datosMD,TableHeadings->{None,cabecera}]

ListLogPlot[{casosCM, hospCM, uciCM, falCM, recCM},PlotLegends->cabecera[[3;;7]], PlotRange->All,Joined->True]

ListPlot[Differences[casosCM], Filling->Axis]

Castilla León

{infectadosCL, hospitalizadosCL, muertosCL, recuperadosCL}=Map[Transpose[{rango,#}]&,{casosCL, hospCL, falCL, recCL}];

TableForm[datosCL,TableHeadings->{None,cabecera}]

Esta información nos debe permitir ajustar parámetros del modelo

ListLogPlot[{casosCL, hospCL, uciCL, falCL},PlotLegends->cabecera[[3;;7]], PlotRange->All,Joined->True]

DateListPlot[Transpose[{date, casosCL}]]

ListPlot[Differences[casosCL], Filling->Axis]

Andalucía (el criterio es distinto a las comunidades anteriores, aqui los datos de hospitalizados y UCI con los acumulados)

{infectadosAN, hospitalizadosAN, muertosAN, recuperadosAN}=Map[Transpose[{rango,#}]&,{casosAN, hospAN, falAN, recAN}];

TableForm[datosAN,TableHeadings->{None,cabecera}]

ListLogPlot[{casosAN, hospAN, uciAN, falAN, recAN},PlotLegends->cabecera[[3;;7]], PlotRange->All,Joined->True]

ListPlot[Differences[casosAN], Filling->Axis]

Modelos  compartimentales en epidemiología

Los modelos compartimentales tienen numerosas aplicaciones. En epidemiología podemos modelar un proceso epidemiológico por compartimentos donde cada compartimento representa un estado de salud, Las conexiones entre compartimentos indican la dirección y la tasa de cambio de un estado de salud a otro.

Notación

t : tiempo en días, desde el inicio de la epidemia.

S  	Individuos susceptibles de contagiarse.

I 	Individuos infectados.

ℰ 	(del inglés exposed), Individuos que portan la enfermedad pero aun no muestran síntomas (Incubación)

R 	Individuos recuperados

β 	Tasa de contagios (probabilidad de que una persona enferme al estar en contacto con un infectado); 1/β  es el tiempo promedio en que cada infectado trasmite la enfermedad a otra. (Ej.: β = 1/2 significa que cada infectado trasmite la enfermedad a una persona cada dos días)

λ[t] es la fuerza de la infección (force of infection) que es dependiente del tiempo. Tiene en cuenta que durante el periodo que un individuo infecta, su “capacidad de infección varia”.

γ 	Tasa de recuperación; 1/γ es  la duración media de la infección;   (Ej.: γ = 1/14  significa que un infectado se recupera en 14 días)

ζ	Tasa de incubación. 1/ ζ es la duración media de la incubación. (Ej.: ζ  = 1/7  significa que el periodo d e incubación es de 7 días)

𝒩 - Población total.

No usados por ahora

μ - Tasa promedio de defunciones (probabilidad de que un individuo infectado muera debido a la enfermedad).

f - Tasa promedio de pérdida de inmunidad en individuos recobrados.

1/ε Tiempo promedio de incubación.

B - Tasa promedio de nacimientos.

M - Infantes con inmunidad pasiva.

δ - Tiempo promedio de inmunidad temporal.

Como ejemplo de modelo complejo puede verse: https://community.wolfram.com/groups/-/m/t/1897377

Modelo 1: SIR

Un modelo muy simple  es el SIR. En su versión mas simple está formado por tres compartimentos que representan tres situaciones: (S) susceptible (de ser infectado), (I) infectado y (R) recuperados.

Esquemáticamente lo podemos representar como se muestra:

({{𝒮(Overscript[→, β λ[t]])ℐ, "infection"}, {ℐ(Overscript[→, γ])ℛ, "recovery"}})

Una versión muy simple de este modelo consiste en sustituir S, I, and R  por las variables  s, i, y r que varian entre 0 y 1 (que representa toda la población).

Inicialmente suponemos que casi toda la población es susceptible de contagio (es decir s(0)=casi 1), ira disminuyendo con el tiempo al extenderse la infección. La población infectada sera en t= 0 casi 0, (i(0) casi 0), igualmente al inicio no habrá nadie recuperado (r(0)=0) . En cualquier instante debe cumplirse que s + i + r = 1. Se supone que cada individuo infectado transmite su enfermedad a b individuos nuevos cada día, mientras mayor sea b mas rápidamente se expandirá  la enfermedad,
k representa la fracción de la población infectada que se recupera cada día (Por ejemplo, si la enfermedad dura 7 días, entonces k es 1/7).

Clear[t,s,i,r, b, k]

La primera ecuación describe cómo la población susceptible cambia con el tiempo.  Debido a que comenzamos con una población susceptible de 1 que disminuye con el tiempo, podemos suponer que la tasa de cambio será negativa.  Además, la tasa de cambio es proporcional al número de contactos b (más contactos disminuye la población susceptible más rápido). Además, la tasa de cambio es proporcional a la población que en cada momento es susceptible de ser infectada (pero no con la población recuperada):

equationS=
s'[t]==-b s[t] i[t];

La segunda ecuación describe cómo la población infectada cambia con el tiempo.  El lado derecho tiene dos términos:   El primero corresponde al lado derecho de la ecuación anterior cambiada de signo que representa el aumento de la población infectada (si lo pensamos como un compartimento es flujo del compartimento s que entre en el i) , el segundo término  representa la disminución (signo -) de casos infectados (personas que se recuperan):

equationI=
i'[t]==b s[t]i[t]-k i[t];

La última ecuación representa las personas (la fracción del total) que se recuperan (en términos compartimentales es el  flujo del compartimento i al compartimento s)

equationR=
r'[t]==k i[t];

Vamos a resolver la ecuación diferencial  para b = 1/2 (cada infectado trasmite la enfermedad a una persona cada dos días)  y k = 1/14 (esto es: suponemos que un infectado se recupera en 14 dias y durante ese tiempo puede contagiar la enfermedad). Ademas supondremos s[0]==0.9999, i[0]==0.0001,r[0]==0.0000.  Al emplear un método numérico  NDSolve  para resolverla hemos de especificar el tiempo que analizados. Suponemos t = 100.

solution = Block[{b=0.5,k=1/14.0} , NDSolve[{ equationS,equationI,equationR,s[0]==0.9999,i[0]==0.0001,r[0]==0.0000},{s, r, i},{t,100}]]

Ahora podemos representar la solución como sigue:

{sols,solr,soli}=First[ {s, r, i}/.solution]

Plot[{sols[t],soli[t],solr[t]},{t,0,100},PlotRange->{0,1.01},PlotStyle->{Blue,Red,Green}]

Ahora vamos construir un modelo que  la solución en función de b y k.

sirmodel1 = {s'[t]==-b s[t] i[t], i'[t]==b s[t]i[t]-k i[t],r'[t]==k i[t]};

solucion1=ParametricNDSolve[Join[sirmodel1,{s[0]==0.9999,i[0]==0.0001,r[0]==0.0000}],{s, i, r},{t,0,100},{ b, k}]

{sols1,solr1,soli1}= {s, r, i}/.solucion1

Plot[{sols1[0.5, 1/14][t], solr1[0.5, 1/14][t], soli1[0.5, 1/14][t]}, {t, 0, 100},PlotStyle->{Blue,Red,Green}]

Otra forma

solucion2=ParametricNDSolveValue[Join[sirmodel1,{s[0]==0.9999,i[0]==0.0001,r[0]==0.0000}],{s,  r, i},{t,0,100},{ b, k}];

{sols2,solr2,soli2}=solucion2[0.5, 1/14]

Plot[{sols2[t], solr2[t], soli2[t]}, {t, 0, 100},PlotStyle->{Blue,Red,Green}]

Alternativamente podemos representar el total acumulado

Plot[{
sols2[t],
sols2[t]+soli2[t],
sols2[t]+ soli2[t]+ solr2[t]
},{t,0,100},
PlotRange->{0,1.01},
Filling->{1->Axis,2->{1},3->{2}},
PlotStyle->{Blue,Red,Green}
]

Regresión

Supongamos que no conocemos b y k y lo queremos estimar a partir de los datos del número de personas infectadas

sample= Table[{t, soli1[0.5+RandomReal[{-0.1,0.1}], 1/14+RandomReal[{-0.01,0.01}]][t]}, {t, 10, 100, 5}]

ListPlot[sample]

fit=FindFit[sample, soli1[ b,k][t],{{b, 0.5},{k,0.01}},t]

Plot[ soli1[b,k][t]/.fit, {t, 10, 100}, Epilog->{Point@sample}]

Modelo 2: SEIR

Basado en: Epidemiological Models for Influenza and COVID - 19—part 1. Robert B. Nachbar. Revised: 17-Mar-2020

El modelo SEIR considera una nueva clase de individuos  ℰ (del inglés exposed), es decir, aquellos que portan la enfermedad pero que al hallarse en su periodo de incubación no muestran síntomas y pueden o no estar en condición de infectar a otros (el  un individuo  nunca obtiene inmunidad). También incluye a los individuos recuperados. Cada grupo es mutuamente excluyente y la suma de todos es la población total. No incluye los fallecimientos (por que no existan o por que su proporción sea muy pequeña comparada con los otros grupos)

({{𝒮(Overscript[→, β λ[t]])ℰ, "infection"}, {ℰ(Overscript[→, ζ])ℐ, "incubation"}, {ℐ(Overscript[→, γ])ℛ, "recovery"}})

El modelo

SEIRmodel={ℰ'[t]==-ζ ℰ[t]+β 𝒮[t] λ[t],ℐ'[t]==ζ ℰ[t]-γ ℐ[t],ℛ'[t]==γ ℐ[t],𝒮'[t]==-β 𝒮[t] λ[t]};

forceOfInfectionSEIR={λ[t_]:>ℐ[t]};

{susceptibleSEIR,infectedSEIR,exposedSEIR,recoveredSEIR}={𝒮,ℐ,ℰ,ℛ}/.ParametricNDSolve[Join[SEIRmodel/.forceOfInfectionSEIR,{𝒮[0]==𝒩-I0,ℰ[0]==0,ℐ[0]==I0,ℛ[0]==0}],{𝒮,ℐ,ℰ,ℛ},{t,0,100},{𝒩,I0,β,ζ,γ}];

EJEMPLO

Este ejemplo corresponde a un colegio internado inglés (Martcheva, pp. 144-145 [MM][HJR]) de 578 estudiantes.

An epidemic of influenza began on 15 January 1978. One hundred and sixty-six boys (29 per cent) were treated in the sick bays by bed rest and aspirin. Certainly this is not the total number with influenza, as the older boys had their own rooms in which some remained, treating themselves and avoiding detection and supervision as a result of the dislocated curriculum.

influenzaData1=

Keys@influenzaData1

influenzaData1["NumberNotInClass"]

{firstDate1,lastDate1}={(DateObject["Day: ""Sun 22 Jan 1978", DateObject[{1978,1,22},"Day","Gregorian",-5.`]]), (DateObject["Day: ""Sat 4 Feb 1978", DateObject[{1978,2,4},"Day","Gregorian",-5.`]])};

Vamos a ajustar el modelo SEIR utilizando los datos del número de estudiantes que permanecen en cama, que forman el grupo de los infectados

fitData1={{0,0},{1,6},{2,26},{3,74},{4,223},{5,294},{6,258},{7,237},{8,191},{9,124},{10,68},{11,27}};

Otros datos del modelo son:  
El tamaño de la población(𝒩):

𝒩 = 763;

El número de infectados en t = 0.

I0 = 1;

Statistical fit

fit1=NonlinearModelFit[fitData1,infectedSEIR[763,1,β,ζ,γ][t],{{β,0.012},{ζ,0.61},{γ,0.43}},t]

fit1["EstimatedVariance"]

fit1["ParameterTable"]

fitWithDataPlot[fit1,{firstDate1,lastDate1}]

residualsPlot[fit1]

Manual fit

DynamicModule[{localData,localPop,localSusceptible,localExposed,localInfected,localRecovered,beta,zeta,gamma,tMax=50},
localPop=populationSize1;
localData=fitData1;
{localSusceptible,localExposed,localInfected,localRecovered}={susceptibleSEIR,exposedSEIR,infectedSEIR,recoveredSEIR};
Manipulate[
beta=10.^logBeta;
zeta=10.^logZeta;
gamma=10.^logGamma;
Show[
Plot[{
localSusceptible[localPop,1,beta,zeta,gamma][t],
localExposed[localPop,1,beta,zeta,gamma][t],
localInfected[localPop,1,beta,zeta,gamma][t],
localRecovered[localPop,1,beta,zeta,gamma][t]
},{t,0,tmax},
PlotRange->{All,{0,ymax}},
PlotRangePadding->{Automatic,{Scaled[0.03],Scaled[0.03]}},
PlotLegends->{𝒮,ℰ,ℐ,ℛ},
FrameLabel->{"time","# people"},
GridLines->{None,{0,influenzaData1["NumberNotInClass"]}},
PlotStyle->{
ColorData["Crayola"]["Green"],
ColorData["Crayola"]["Sunglow"],
ColorData["Crayola"]["Scarlet"],
ColorData["Crayola"]["NavyBlue"]
}
],
ListPlot[localData,PlotStyle->ColorData["Crayola"]["Scarlet"]]
],

{{logBeta,Log10@0.012,Dynamic@StringForm["β = `` (infection)",TraditionalForm[10.^logBeta]]},Log10@0.0001,Log10@1.0},
{{logZeta,Log10@0.61,Dynamic@StringForm["ζ = `` (incubation)",TraditionalForm[10.^logZeta]]},Log10@(1/5),Log10@(1/(*0.125*)0.01)},
{{logGamma,Log10@0.43,Dynamic@StringForm["γ = `` (recovery)",TraditionalForm[10.^logGamma]]},Log10@(1/5),Log10@(1/(*1*)0.1)},
Delimiter,
{{tmax,15,"(Subscript[t, max])"},5,tMax},
{{ymax,800,"(Subscript[y, max])"},50,800},

TrackedSymbols:>{logBeta,logZeta,logGamma,tmax,ymax},
SaveDefinitions->True
]
]

Modelo 3: SIRM ({β,γ, μ , k}, con (β/(1+(e^(μ (t^k))))),  Ntotal=10^6;I0=10;)

http : //www.imuva.uva.es/covid19/COVID19 - AnalisisComparadoChinaCyL.html

AL modelo SIR se le añade una función de mitigación  β(t)  >0 que modele la forma en que la tasa de contagio va disminuyendo en el tiempo,

Una opción es (UVA toma x0=0)

Clear[𝒮,ℐ,ℛ, β,γ, μ,k]

β[t_]=(β/(1+(E^(μ (t^k)))));

SIRMmodel={𝒮'[t]==-β[t] 𝒮[t] ℐ[t]/Ntotal,ℐ'[t]==β[t] 𝒮[t] ℐ[t]/Ntotal-γ ℐ[t],ℛ'[t]==γ ℐ[t]};

Ntotal=10^6;I0=10;

{susceptibleSIRM,infectedSIRM,recoveredSIRM}={𝒮,ℐ,ℛ}/.ParametricNDSolve[Join[SIRMmodel,{𝒮[0]== Ntotal,ℐ[0]==I0,ℛ[0]==0}],{𝒮,ℐ,ℛ},{t,0,50},{β,γ, μ , k}];

Modelo 4: SIRV ({β,γ,ν}, Ntotal=10^6;I0=10;)

AL modelo SIR se le añade un termino v que tiene en cuenta la inmunizacion de parte de la población. Puede ser de forma natural o por el empleo de vacunas, que no es el caso.

Clear[𝒮,ℐ,ℛ, β,γ, ν]

Ntotal=10^6;I0=10;

SIRVmodel={𝒮'[t]==-β 𝒮[t] ℐ[t]/Ntotal- ν 𝒮[t],ℐ'[t]==β 𝒮[t] ℐ[t]/Ntotal-γ ℐ[t],ℛ'[t]==γ ℐ[t]+ ν 𝒮[t]};

{susceptibleSIRV,infectedSIRV,recoveredSIRV}={𝒮,ℐ,ℛ}/.ParametricNDSolve[Join[SIRVmodel,{𝒮[0]== Ntotal,ℐ[0]==I0,ℛ[0]==0}],{𝒮,ℐ,ℛ},{t,0,50},{β,γ,ν}];

Plot[infectedSIRV[0.1, 0.1, 0.1][t], {t, 1, 50}]

Modelo 5:  SEIRD Model with standard incidence

SEIRD: SEIR model with dead compartment

An effective population size 𝒩 with standard incidence will be used

This assumes that individuals in the population cannot come in contact with all other individuals in the population; this is in contrast to models for influenza and SARS [MM]

𝒩 may be adjusted for individual locations; initially we'll use 100000 individuals

odesSEIRD={
𝒮'[t]==-β 𝒮[t] λ[t],
ℰ'[t]==-ζ ℰ[t]+β 𝒮[t] λ[t],
ℐ'[t]==ζ ℰ[t]-γ ℐ[t]-μ ℐ[t],
ℛ'[t]==γ ℐ[t],
𝒟'[t]==μ ℐ[t]
};
forceOfInfectionSEIRD={λ[t_]:>((ℐ[t])/(𝒩-𝒟[t]))};
varsSEIRD=Cases[odesSEIRD,(v_)'[t]:>v,∞];

{susceptibleSEIRD,exposedSEIRD,infectedSEIRD,recoveredSEIRD,deadSEIRD}={𝒮,ℰ,ℐ,ℛ,𝒟}/.ParametricNDSolve[Join[odesSEIRD/.forceOfInfectionSEIRD,{𝒮[0]==𝒩-I0,ℰ[0]==0,ℐ[0]==I0,ℛ[0]==0,𝒟[0]==0}],varsSEIRD,{t,0,500},{𝒩,I0,β,ζ,γ,μ}];

Modelo 6:  SEIsaRD Model

SEIsaRD: SEIR model with symptomatic and asymptomatic infectious compartments

Standard incidence

The duration of infection for symptomatic and asymptomatic infectious individuals is the same

We have two different recovered compartments for book keeping

death due to disease occurs only for symptomatic infectious individuals

These are the transitions:

Grid[modelGrid={{𝒮(Overscript[→, β λ[t]])ℰ ,"infection"},{ℰ(Overscript[→, (Subscript[ρ, s]) ζ])(Subscript[ℐ, "s"]) ,"incubation, symptomatic"},{ℰ(Overscript[→, (1-(Subscript[ρ, s])) ζ])(Subscript[ℐ, "a"]) ,"incubation, asymptomatic"},{(Subscript[ℐ, "s"])(Overscript[→, γ])(Subscript[ℛ, "s"]) ,"recovery, symptomatic"},{(Subscript[ℐ, "a"])(Overscript[→, γ])(Subscript[ℛ, "a"]) ,"recovery, asymptomatic"},{(Subscript[ℐ, "s"])(Overscript[→, μ])𝒟 ,"death, symptomatic"}},Alignment->{{Center,Left}},Spacings->3,BaseStyle->{18}]
closeCellGroup[]

model=First/@modelGrid;
CompartmentalModelGraph[model,GraphLayout->{"LayeredDigraphEmbedding","RootVertex"->𝒮,"Orientation"->Left},PerformanceGoal->"Quality",BaseStyle->{18},ImageSize->Large]
closeCellGroup[]

odesSEIsaRD={
𝒮'[t]==-β 𝒮[t] λ[t],
ℰ'[t]==β 𝒮[t] λ[t]-ρ ζ ℰ[t]-(1-ρ)ζ ℰ[t],
(Subscript[ℐ, "s"]^′)[t]==ρ ζ ℰ[t]-γ (Subscript[ℐ, "s"])[t]-μ (Subscript[ℐ, "s"])[t],
(Subscript[ℐ, "a"]^′)[t]==(1-ρ)ζ ℰ[t]-γ (Subscript[ℐ, "a"])[t],
(Subscript[ℛ, "s"]^′)[t]==γ (Subscript[ℐ, "s"])[t],
(Subscript[ℛ, "a"]^′)[t]==γ (Subscript[ℐ, "a"])[t],
𝒟'[t]==μ (Subscript[ℐ, "a"])[t]
};
forceOfInfectionSEIsaRD={λ[t_]:>(((Subscript[ℐ, "s"])[t]+(Subscript[ℐ, "a"])[t])/(𝒩-𝒟[t]))};
varsSEIsaRD=Cases[odesSEIsaRD,(v_)'[t]:>v,∞]

{susceptibleSEIsaRD,exposedSEIsaRD,infectedSymptomaticSEIsaRD,infectedAsymptomaticSEIsaRD,recoveredSymptomaticSEIsaRD,recoveredAsymptomaticSEIsaRD,deadSEIsaRD}={𝒮,ℰ,(Subscript[ℐ, "s"]),(Subscript[ℐ, "a"]),(Subscript[ℛ, "s"]),(Subscript[ℛ, "a"]),𝒟}/.ParametricNDSolve[Join[odesSEIsaRD/.forceOfInfectionSEIsaRD,{𝒮[0]==𝒩-I0,ℰ[0]==0,(Subscript[ℐ, "s"])[0]==ρ I0,(Subscript[ℐ, "a"])[0]==(1-ρ)I0,(Subscript[ℛ, "s"])[0]==0,(Subscript[ℛ, "a"])[0]==0,𝒟[0]==0}],varsSEIsaRD,{t,0,500},{𝒩,I0,β,ζ,ρ,γ,μ}];

Ajustes de parámetros y simulación

Necesitamos construir un modelo  para cada comunidad autónoma que permita predecir los siguientes estados :

a. Número de ingresos en UCI
b. Número de enfermos hospitalizados[hospitalizados]
c. Número de fallecimientos[fallecidos]
d. Número de nuevos casos[nuevos]
e. Número de casos confirmados[confirmados]

Utilizamos la información disponible que la hemos descargado al inicio de https://covid19.isciii.es

cabecera={"CCAA Codigo ISO", "Fecha", "Casos", "Hospitalizados", "UCI", "Fallecidos", "Recuperados"};

t : momento de inicio de la pandemia a efectos del estudio. 1 de marzo

Compartimentos

S  	Individuos susceptibles de contagiarse.

Toda la población -  Infectados acumulados en el momento del estudio

ℰ 	(del inglés exposed), Individuos que portan la enfermedad pero aun no muestran síntomas Incubación

I 	Individuos infectados. La cantidad  real (apdo d) son desconocido, lo probable es que corresponda a una proporción de los confirmados (e). La duda es si es simplemente un factor constante de proporcionalidad o no. En principio consideramos (e): Corresponde a los  k*( e) donde  k es desconocido. ? Ademas de los infectados podemos distinguir: infectados en Sintomáticos y Asintomáticos, ademas de los sintomáticos podemos diferenciar los hospitalizados y los no hospitalizados y dentro de los hospitalizados podemos a su vez distinguir los que acaban en UCI y los que no. Si lo anterior lo convertimos en proporciones el modelo se simplifica mucho.

R 	Individuos recuperados (se dispone de la información).

D 	Muertos (disponemos los datos, probablemente es la mas fiable).

Se trataria de construir un modelo que tenga como compartimentos al menos : S, I, ℰ, R, D. Cabe la posibilidad de dividir los infectados en Sintomáticos y Asintomáticos

Sugiero explorar los modelos incluidos en : Epidemiological Models for Influenza and COVID19fullversion.nb. 
https://community.wolfram.com/groups/-/m/t/1897377

https://community.wolfram.com/groups/-/m/t/1895675

Ajustes SIR (pruebas)

No tiene ninguna validez pretenden mostrar el método que vamos a seguir. los ajustes pueden ser buenos pero carecer de sentido los valores de los parámetros

ClearAll[b,k, fit1]

datosAN1=Transpose[{Range[15, 13+Length[casosAN]] ,Differences[casosAN]/1000}]

fit1=NonlinearModelFit[datosAN1, soli1[ b,k][t],{{b, 0.5},{k,0.01}},t]

Plot[ fit1[t], {t, 10, 35},AxesLabel->{"t (days)","Casos/10^6"}, Epilog->{Point@datosAN1}]

fit1["EstimatedVariance"]

fit1["ParameterTable"]

fitWithDataPlot[fit1,{DateObject[{2020, 3, 20}],DateObject[{2020, 4, 30}]}]

residualsPlot[fit1]

Ajustes SIRM

Inicialmente parece que la única información disponible fiable es los hospitalizados, los datos de alta y los muertos, y sobre esto hay dudas sobre la metodología. En definitiva nos quedamos con los hospitalizados y los recuperados.

({{𝒮(Overscript[→, β λ[t]])ℐ, "infection"}, {ℐ(Overscript[→, γ])ℛ, "recovery"}})

Es necesario llegar a unos valores creibles para los parámetros

β 	Tasa de contagios (probabilidad de que una persona enferme al estar en contacto con un infectado); 1/β  es el tiempo promedio en que cada infectado trasmite la enfermedad a otra. (Ej.: β = 1/2 significa que cada infectado trasmite la enfermedad a una persona cada dos días). Deberla ser un valor < 0.5.

λ[t] es la fuerza de la infección (force of infection) que es dependiente del tiempo. Tiene en cuenta que durante el periodo que un individuo infecta, su “capacidad de infección varia”.

γ 	Tasa de recuperación; 1/γ es  la duración media de la infección;   (Ej.: γ = 1/14  significa que un infectado se recupera en 14 días). <0.1

ζ	Tasa de incubación. 1/ ζ es la duración media de la incubación. (Ej.: ζ  = 1/7  significa que el periodo d e incubación es de 7 días)

Se añade  un parámetro de ‘traslación’ (o retardo) y μ, k que corresponden a los  coeficientes de (β/(1+(e^(μ (t^k))))) la  para conseguir que la curva de hospitalizados encaje mejor con los datos:

Clear[𝒮,ℐ,ℛ, β,γ, μ,k]

infectedSIRM

ListPlot[{hospitalizadosCM,recuperadosCM}]

Manipulate[Plot[{infectedSIRM[β,γ, μ, k][t-c],recoveredSIRM[β,γ, μ, k][t-c]}, {t, Max[1,c], 50}, PlotLegends->{"Hospitalizados","Recuperados"}], {{β, 0.68}, 0, 1},{{γ, 0.073}, 0, 1}, {{μ, 0.00098}, 0, 1}, {{k,2}, 0.0001,3} ,{{c,-10.5}, -20, 20}]

Initialization

SetDirectory@NotebookDirectory[];
Get["CompartmentalModeling.wl"]

General

closeCellGroup[] := (SelectionMove[EvaluationCell[], Next,Cell]; 
FrontEndExecute[FrontEndToken["SelectionCloseUnselectedCells"]])

$PlotTheme="Scientific"; (* ColorData[108] *)

Scan[SetOptions[#1,{Axes->False,Frame->True,PlotRange->All}]&,{Plot,LogPlot,LogLinearPlot,LogLogPlot,ListPlot,ListLogPlot,ListLogLinearPlot,ListLogLogPlot,ListLinePlot,ListStepPlot,DateListPlot,DateListStepPlot,DateListLogPlot,Histogram,DateHistogram}
];

Scan[SetOptions[#1,{PlotMarkers->Automatic}]&,{ListPlot,ListLogPlot,ListLogLinearPlot,ListLogLogPlot,DateListPlot,DateListLogPlot}
];

labeledGrid[grid_, columnLabels_, rowLabels_, opts:OptionsPattern[]] := 
	Grid[MapThread[Prepend, {Prepend[grid, columnLabels], Prepend[rowLabels, ""]}], 
		opts,
		Dividers -> {{False, True, False}, {False, True, False}}, 
		Alignment -> {{{Center}}, {{Top}}}
		]

tfn[data_] := TableForm[data, TableDepth -> 2, TableHeadings -> Automatic]
tfc[data_] := TableForm[Rest[data], TableDepth -> 2, TableHeadings -> {None, First[data]}]

Fitting data

(* time is measured in Days *)
toModelTime[t0_DateObject, date_DateObject] := 
	QuantityMagnitude[DateDifference[t0, date]]
toDataDate[t0_DateObject, time:(_Integer|_Real)] := 
	DatePlus[t0, time]

ClearAll[fitDataWDR]
fitDataWDR::usage = "fitDataWDR[locale, startDate] gives data for infected, recovered, and dead \
relative to the given date for the onset of infection for the given locale, either an administrative \
division or a country. For reference, cases were first reported on 31 Dec 2019."; 

fitDataWDR::locnf = "Locale `` not found."; 

Options[fitDataWDR] = {"makeCorrectionForHubei" -> False, "dateRange" -> All}; 

fitDataWDR[locale:Entity["AdministrativeDivision", _] | Entity["Country", _], date0_, opts:OptionsPattern[]] := 
	Module[{covid19Data, confirmed, recovered, dead, correctionQ, dateRange}, 
		covid19Data = ResourceData["Epidemic Data for Novel Coronavirus COVID-19"]; 
		
		Switch[locale, 
			Entity["AdministrativeDivision", _], 
				If[ !MemberQ[covid19Data[All, "AdministrativeDivision"], locale], 
					Message[fitDataWDR::locnf, locale]; 
					Return[$Failed]
					], 
			Entity["Country", _], 
				If[ !MemberQ[covid19Data[All, "Country"], locale], 
					Message[fitDataWDR::locnf, locale]; 
					Return[$Failed]
					]
			]; 
		
		correctionQ = OptionValue["makeCorrectionForHubei"]; 
		If[MatchQ[correctionQ, Except[True | False]], 
			Message[fitDataWDR::opttf, "makeCorrectionForHubei", correctionQ]; 
			Return[$Failed]
			]; 
		
		dateRange = OptionValue["dateRange"]; 
		dateRange = Switch[dateRange, 
			All, 
				covid19Data[[1,"ConfirmedCases"]] /@ {"FirstDate", "LastDate"}, 
			{_, _}, 
				dateRange = (Quiet[DateObject[#1, "Instant"]] & ) /@ 
					(dateRange //. {{All, b_} :> {covid19Data[[1,"ConfirmedCases"]]["FirstDate"], b}, 
						{a_, All} :> {a, covid19Data[[1,"ConfirmedCases"]]["LastDate"]}}); 
				If[MatchQ[dateRange, {DateObject[{__}, _String, _String, _Real], DateObject[{__}, _String, _String, _Real]}], 
					dateRange, 
					Message[fitDataWDR::optvg, "dateRange", OptionValue["dateRange"], "All or a list of two dates"]; 
					Return[$Failed]
					], 
			_, 
				Message[fitDataWDR::optvg, "dateRange", OptionValue["dateRange"], "All or a list of two dates"]; 
				Return[$Failed]
			]; 
		
		Switch[locale, 
			Entity["AdministrativeDivision", _], 
				confirmed = covid19Data[SelectFirst[#AdministrativeDivision == locale & ], "ConfirmedCases"]; 
				recovered = covid19Data[SelectFirst[#AdministrativeDivision == locale & ], "RecoveredCases"]; 
				dead = covid19Data[SelectFirst[#AdministrativeDivision == locale & ], "Deaths"], 
			Entity["Country", _], 
				confirmed = Plus @@ Normal[covid19Data[Select[#Country == locale & ], "ConfirmedCases"]]; 
				recovered = Plus @@ Normal[covid19Data[Select[#Country == locale & ], "RecoveredCases"]]; 
				dead = Plus @@ Normal[covid19Data[Select[#Country == locale & ], "Deaths"]]
			]; 
		
		confirmed = If[correctionQ && 
				MatchQ[locale, Entity["AdministrativeDivision", {"Hubei", "China"}] | Entity["Country", "China"]], 
			TimeSeries[Apply[{#1, If[#1 >= DateObject["13 Feb 2020", "Instant"], #2 - 13332, #2]} & , Normal[confirmed], {1}]], 
			confirmed
			]; 
		
		confirmed = TimeSeriesWindow[confirmed, dateRange]; 
		recovered = TimeSeriesWindow[recovered, dateRange]; 
		dead = TimeSeriesWindow[dead, dateRange]; 
		
		Map[{QuantityMagnitude[DateDifference[date0, First[#1], "Day"]], Last[#1]} & , 
			Normal /@ {confirmed - (recovered + dead), recovered, dead}, {2}]
		]

(*
fitDataWDR[country:Entity["Country",_],date0_,opts:OptionsPattern[]]/;MemberQ[covid19Data[All,"Country"],country]:=
Module[{confirmed,recovered, dead,correctionQ},
correctionQ=OptionValue["makeCorrectionForHubei"];
If[MatchQ[correctionQ,Except[True|False]],
Message[fitDataWDR::opttf,"makeCorrectionForHubei",correctionQ];
Return[$Failed]
];

confirmed=Plus@@Normal@covid19Data[Select[#Country==country&],"ConfirmedCases"];
recovered=Plus@@Normal@covid19Data[Select[#Country==country&],"RecoveredCases"];
dead=Plus@@Normal@covid19Data[Select[#Country==country&],"Deaths"];

confirmed=If[correctionQ && country==Entity["Country","China"],
{#1,If[#1≥DateObject["13 Feb 2020","Instant"],#2-13332,#2]}&@@@Normal@confirmed//TimeSeries,
(* else *)
confirmed
];

Map[{QuantityMagnitude@DateDifference[date0,First@#,"Day"],Last@#}&,Normal/@{confirmed-(recovered+dead),recovered,dead},{2}]
]*)

Fitting error

Exp10[x_] := 10^x

toLog[x_] := Log[x]
fromLog[y_] := Exp[y]
toLogit[x_] := Log[x/(1-x)]
fromLogit[y_]:=Exp[y]/(1 + Exp[y])

Plot[toLogit[x],{x,0,1}]

toLog[x]//fromLog

toLogit[x]//fromLogit//Simplify[#,Assumptions->x∈Reals]&

ClearAll[sumSquaredError]
sumSquaredError::usage = "sumSquaredError[data, model, parameters] gives the sum of \
squared errors for the model with the given parameters for the given data. The \
data and model can be lists of the same length."; 
sumSquaredError::len = "The length of the data list `` is not the same as the length \
of the model list ``."; 
sumSquaredError[fitData:{{_?NumberQ, _?NumberQ}..}, model:Except[{__}], params_] := 
	sumSquaredError[{fitData}, {model}, params]
sumSquaredError[fitDataList:{{{_?NumberQ, _?NumberQ}..}..}, modelList:{__}, 
		params:{(_ -> _?NumberQ)..}] /; Length[fitDataList] == Length[modelList] || 
		Message[sumSquaredError::len, Length[fitDataList], Length[modelList]] := 
	Module[{time, observed, calculated, sse}, 
		sse = MapThread[
			Function[{data, model}, 
				time = First /@ data; 
				observed = Last /@ data; 
				calculated = Quiet[Check[model /@ time /. params, $Failed]]; 
				(calculated - observed)^2
				], 
			{fitDataList, modelList}
			] // Total[#, 2]&; 
		If[FreeQ[sse, $Failed], 
			sse, 
		  (* else *)
			Sqrt[$MaxMachineNumber]
			]
		]

SEIRD

ClearAll[fitErrorSEIRD]
fitErrorSEIRD[fitData_, 
		pop_, I0_, 
		(logBeta_)?NumberQ, (logZeta_)?NumberQ, (logGamma_)?NumberQ, (logMu_)?NumberQ] := 
	Module[{params = Thread[{β, ζ, γ, μ} -> Exp[{logBeta, logZeta, logGamma, logMu}]]}, 
		sumSquaredError[
			fitData, 
			{infectedSEIRD[pop, I0, β, ζ, γ, μ], 
				recoveredSEIRD[pop, I0, β, ζ, γ, μ], 
				deadSEIRD[pop, I0, β, ζ, γ, μ]}, 
			params
			]
		]

ClearAll[plotFitSEIRD]
plotFitSEIRD[fitData_, {pop_, I0_, beta_, zeta_, gamma_, mu_}, opts:OptionsPattern[{Plot, ListPlot}]] := 
	Module[{tmax = 100, ymax = 40000, scale = 50000/pop}, 
		Show[
			Plot[{
				scale*susceptibleSEIRD[pop, I0, beta, zeta, gamma, mu][t], 
				exposedSEIRD[pop, I0, beta, zeta, gamma, mu][t], 
				infectedSEIRD[pop, I0, beta, zeta, gamma, mu][t], 
				recoveredSEIRD[pop, I0, beta, zeta, gamma, mu][t], 
				deadSEIRD[pop, I0, beta, zeta, gamma, mu][t]} // Evaluate, 
				{t, 0, tmax}, 
				PlotRangePadding -> {Automatic, {Scaled[0.03], Scaled[0.03]}}, 
				PlotLegends -> {N[scale]*𝒮, ℰ, ℐ, ℛ, 𝒟}, 
				PlotLabel -> StringForm["𝒩 = ``, SSE = ``", pop, 
					fitErrorSEIRD[fitData, pop, I0, 
						Log[beta], Log[zeta], Log[gamma], Log[mu]]], 
				FrameLabel -> {"time", "# people"}, 
				PlotStyle -> {
					ColorData["Crayola"]["Green"], 
					ColorData["Crayola"]["Sunglow"], 
					ColorData["Crayola"]["Scarlet"], 
					ColorData["Crayola"]["NavyBlue"], 
					ColorData["Crayola"]["OuterSpace"]
					},
				Evaluate@FilterRules[{opts}, Options[Plot]], 
				PlotRange -> {{0, tmax}, {0, ymax}}, 
				GridLines -> {None, None}
				], 
			ListPlot[fitData, 
				PlotStyle -> {
					ColorData["Crayola"]["Scarlet"], 
					ColorData["Crayola"]["NavyBlue"], 
					ColorData["Crayola"]["OuterSpace"]
					},
				Evaluate@FilterRules[{opts}, Options[ListPlot]]
				]
			]
		]

SEIQRD

ClearAll[fitErrorSEIQRD]
fitErrorSEIQRD[fitData_, 
		pop_, I0_, 
		(logBeta_)?NumberQ, (logZeta_)?NumberQ, (logGamma_)?NumberQ, (logDelta_)?NumberQ, 
		(logMu_)?NumberQ] := 
	Module[{params = Thread[{β, ζ, γ, δ, μ} -> Exp[{logBeta, logZeta, logGamma, logDelta, logMu}]]}, 
		sumSquaredError[
			fitData, 
			{quarantinedSEIQRD[pop, I0, β, ζ, γ, δ, μ], 
				recoveredSEIQRD[pop, I0, β, ζ, γ, δ, μ], 
				deadSEIQRD[pop, I0, β, ζ, γ, δ, μ]}, 
			params
			]
		]

ClearAll[plotFitSEIQRD]
plotFitSEIQRD[fitData_, {pop_, I0_, beta_, zeta_, gamma_, delta_, mu_}, opts:OptionsPattern[{Plot, ListPlot}]] := 
	Module[{tmax = 100, ymax = 40000, scale = 50000/pop}, 
		Show[
			Plot[{
				scale*susceptibleSEIQRD[pop, I0, beta, zeta, gamma, delta, mu][t], 
				exposedSEIQRD[pop, I0, beta, zeta, gamma, delta, mu][t], 
				infectedSEIQRD[pop, I0, beta, zeta, gamma, delta, mu][t], 
				quarantinedSEIQRD[pop, I0, beta, zeta, gamma, delta, mu][t], 
				recoveredSEIQRD[pop, I0, beta, zeta, gamma, delta, mu][t], 
				deadSEIQRD[pop, I0, beta, zeta, gamma, delta, mu][t]} // Evaluate, 
				{t, 0, tmax}, 
				PlotRangePadding -> {Automatic, {Scaled[0.03], Scaled[0.03]}}, 
				PlotLegends -> {N[scale]*𝒮, ℰ, ℐ, 𝒬, ℛ, 𝒟}, 
				PlotLabel -> StringForm["𝒩 = ``, SSE = ``", pop,
					fitErrorSEIQRD[fitData, pop, I0, 
						Log[beta], Log[zeta], Log[gamma], Log[delta], Log[mu]]], 
				FrameLabel -> {"time", "# people"}, 
				PlotStyle -> {
					ColorData["Crayola"]["Green"], 
					ColorData["Crayola"]["Sunglow"], 
					ColorData["Crayola"]["Scarlet"], 
					ColorData["Crayola"]["Fuchsia"], 
					ColorData["Crayola"]["NavyBlue"], 
					ColorData["Crayola"]["OuterSpace"]
					},
				Evaluate@FilterRules[{opts}, Options[Plot]], 
				PlotRange -> {{0, tmax}, {0, ymax}}, 
				GridLines -> {None, None}
				], 
			ListPlot[fitData, 
				PlotStyle -> {
					ColorData["Crayola"]["Fuchsia"], 
					ColorData["Crayola"]["NavyBlue"], 
					ColorData["Crayola"]["OuterSpace"]
					},
				Evaluate@FilterRules[{opts}, Options[Plot]]
				]
			]
		]

SEIsaRD

ClearAll[fitErrorSEIsaRD]
fitErrorSEIsaRD[fitData_, 
		pop_, I0_, 
		(logBeta_)?NumberQ, (logZeta_)?NumberQ, (logitRho_)?NumberQ, (logGamma_)?NumberQ, 
		(logMu_)?NumberQ] := 
	Module[{params = MapThread[#1 -> #2[#3]&, 
			{{β, ζ, ρ, γ, μ}, 
			{fromLog, fromLog, fromLogit, fromLog, fromLog}, 
			{logBeta, logZeta, logitRho, logGamma, logMu}}]}, 
		sumSquaredError[
			fitData, 
			{infectedSymptomaticSEIsaRD[pop, I0, β, ζ, ρ, γ, μ], 
				recoveredSymptomaticSEIsaRD[pop, I0, β, ρ, ζ, γ, μ], 
				deadSEIsaRD[pop, I0, β, ζ, ρ, γ, μ]}, 
			params
			]
		]

ClearAll[plotFitSEIsaRD]
plotFitSEIsaRD[fitData_, {pop_, I0_, beta_, zeta_, rho_, gamma_, mu_}, opts:OptionsPattern[{Plot, ListPlot}]] := 
	Module[{tmax = 100, ymax = 40000, scale = 50000/pop}, 
		Show[
			Plot[{
				scale*susceptibleSEIsaRD[pop, I0, beta, zeta, rho, gamma, mu][t], 
				exposedSEIsaRD[pop, I0, beta, zeta, rho, gamma, mu][t], 
				infectedSymptomaticSEIsaRD[pop, I0, beta, zeta, rho, gamma, mu][t], 
				infectedAsymptomaticSEIsaRD[pop, I0, beta, zeta, rho, gamma, mu][t], 
				recoveredSymptomaticSEIsaRD[pop, I0, beta, zeta, rho, gamma, mu][t], 
				recoveredAsymptomaticSEIsaRD[pop, I0, beta, zeta, rho, gamma, mu][t], 
				deadSEIsaRD[pop, I0, beta, zeta, rho, gamma, mu][t]} // Evaluate, 
				{t, 0, tmax}, 
				PlotRangePadding -> {Automatic, {Scaled[0.03], Scaled[0.03]}}, 
				PlotLegends -> {N[scale]*𝒮, ℰ, (Subscript[ℐ, "s"]), (Subscript[ℐ, "a"]), (Subscript[ℛ, "s"]), (Subscript[ℛ, "a"]), 𝒟}, 
				PlotLabel -> StringForm["𝒩 = ``, SSE = ``", pop,
					fitErrorSEIsaRD[fitData, pop, I0, 
						toLog[beta], toLog[zeta], toLogit[rho], toLog[gamma], toLog[mu]]], 
				FrameLabel -> {"time", "# people"}, 
				PlotStyle -> {
					ColorData["Crayola"]["Green"], 
					ColorData["Crayola"]["Sunglow"], 
					ColorData["Crayola"]["Scarlet"], 
					Directive[ColorData["Crayola"]["Scarlet"], Dashed], 
					ColorData["Crayola"]["NavyBlue"], 
					Directive[ColorData["Crayola"]["NavyBlue"], Dashed], 
					ColorData["Crayola"]["OuterSpace"]
					},
				Evaluate@FilterRules[{opts}, Options[Plot]], 
				PlotRange -> {{0, tmax}, {0, ymax}}, 
				GridLines -> {None, None}
				], 
			ListPlot[fitData, 
				PlotStyle -> {
					ColorData["Crayola"]["Scarlet"], 
					ColorData["Crayola"]["NavyBlue"], 
					ColorData["Crayola"]["OuterSpace"]
					},
				Evaluate@FilterRules[{opts}, Options[Plot]]
				]
			]
		]

Fit visualization

Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"];

(* time is measured in Days *)
toModelTime[t0_DateObject, date_DateObject] := 
	QuantityMagnitude[DateDifference[t0, date]]
toDataDate[t0_DateObject, time:(_Integer|_Real)] := 
	DatePlus[t0, time]

ClearAll[fitWithDataPlot]
fitWithDataPlot[fit_FittedModel,{dateMin_DateObject,dateMax_DateObject}]:=Module[{plotData,bands,cd=ColorData[108]},
plotData={toDataDate[dateMin,#1],#2}&@@@fit["Data"];
bands[x_]=Quiet[fit["SinglePredictionBands",ConfidenceLevel->0.95]/.t->x];
Show[
DateListPlot[plotData,
Joined->False,
PlotRange->All
],
Plot[{fit[toModelTime[dateMin,FromAbsoluteTime@t]],bands[toModelTime[dateMin,FromAbsoluteTime@t]]},{t,AbsoluteTime@dateMin,AbsoluteTime@dateMax},
PlotRange->All,
PlotStyle->{cd[2],None},
Filling->{2->{1}},
FillingStyle->{Opacity[0.5,Lighter@cd[2]]}
],
PlotRange->All,
PlotRangePadding->{Automatic,{Scaled[0.03],Scaled[0.1]}},
FrameLabel->{"time (d)","number infected"},
PlotLabel->StringForm["Estimated variance = ``",fit["EstimatedVariance"]],
ImageSize->360
]//Labeled[#,Column[{PointLegend[{cd[1]},{"observed"}],LineLegend[{cd[2]},{"calculated"}],SwatchLegend[{Opacity[0.5,Lighter@cd[2]]},{"95% prediction band"}]},Left],Right]&
]

ClearAll[residualsPlot]
residualsPlot[fit_FittedModel]:=Module[{width=4 72},
GraphicsGrid[{{
ListPlot[Thread[{First/@fit["Data"],fit["FitResiduals"]}],
FrameLabel->{"time (d)","FitResiduals"},
Filling->Axis,
ImageSize->width
],
ListPlot[Thread[{fit["PredictedResponse"],fit["FitResiduals"]}],
FrameLabel->{"PredictedResponse","FitResiduals"},
Filling->Axis,
ImageSize->width
]
}},
Spacings->Scaled[0.05]
]
]

ClearAll[modelSensitivityPlot]
modelSensitivityPlot[modelIn:(_ParametricFunction[__]),fit_FittedModel,{{tMin_?NumberQ,tMax_?NumberQ},{yMin_?NumberQ,yMax_?NumberQ}},scale:{__?NumberQ}]:=Module[{model,sensitivities,cd=ColorData[108],params,t},
params=fit["BestFitParameters"];
model=modelIn[t]/.params;
sensitivities=MapThread[model+(#2{-1,1}D[modelIn,#1][t]/.params)&,{Keys@params,scale}];
TabView@MapThread[Function[{p,bands,sf},p->(Show[
ListPlot[fit["Data"],PlotRange->All],
Plot[{model,bands}//Evaluate,{t,tMin,tMax},
PlotRange->All,
PlotStyle->{cd[2],Lighter@cd[3],Lighter@cd[4]},
Filling->{2->{{1},Opacity[0.5,Lighter@cd[3]]},3->{{1},Opacity[0.5,Lighter@cd[4]]}}
],
PlotRange->{All,{yMin,yMax}},
PlotRangePadding->{Automatic,{Scaled[0.03],Scaled[0.1]}},
FrameLabel->{"time (d)","number infected"},
PlotLabel->"Parameter sensitivity",
Epilog->{Inset[StringForm["scale factor = ``",TraditionalForm[sf]],Scaled[{0.05,0.95}],{-1,1}]}
]//Labeled[#,Column[{PointLegend[{cd[1]},{"observed"}],LineLegend[{cd[2]},{"calculated"}],SwatchLegend[{Opacity[0.5,Lighter@cd[3]]},{"negative sensitivity"}],SwatchLegend[{Opacity[0.5,Lighter@cd[4]]},{"positive sensitivity"}]},Left],Right]&)],{Keys@params,sensitivities,scale}]
]

ClearAll[epidemicSizePlot]
epidemicSizePlot[susceptibleIn:(_ParametricFunction[__]),fit_FittedModel,populationSize_?Positive,actualSize_?Positive,{{tMin_?NumberQ,tMax_?NumberQ},{yMin_?NumberQ,yMax_?NumberQ}}]:=
Module[{susceptible,params,t,domain,size},
params=fit["BestFitParameters"];
susceptible=susceptibleIn/.params;
domain=First@InterpolatingFunctionDomain[susceptible];
size=(populationSize-susceptible[t])/populationSize/.t->Last@domain;
Plot[populationSize-susceptible[t]/.fit["BestFitParameters"]//Evaluate,
{t,tMin,tMax},
PlotRange->{All,{yMin,yMax}},
PlotRangePadding->{Automatic,{Scaled[0.03],Scaled[0.03]}},
PlotLabel->StringForm["Epidemic Size = ``", Quantity[100size,"Percent"]],
FrameLabel->{"time","epidemic size"},
GridLines->{None,{0,actualSize,populationSize}},
PlotStyle->{
ColorData["Crayola"]["Scarlet"]
}
]
]