1-123-T-Mma-PanamaCoronaVirus - Teacher Version

Brian Winkel, Director SIMIODE, Cornwall NY USA
*** DRAFT ***

Source:  Ministerio de Salud de la Repùblica de Panamà.  2020. Casos de Coronavirus COVID-19 en Panamà  
http://minsa.gob.pa/coronavirus-covid19. Accessed 22 March 2020. 

We study the early spread of the Corona virus (COVID-19) in Panama. We are NOT using the traditional epidemiology SIR model with three populations Susceptible, Infected, and Removed. Rather here we are going for a one population model consisting of those who are infected or are cases with the virus.

Gender specific lockdown went into effect 1 April 2020, day 15. Women permitted out on MWF, Men on TThS, No one on Sunday.  Two hour tim elimit based on ID number, e.g., if last digit ends in 7 then hours permitted are 7 AM  to 9 AM, gender permitting.

Data Input

We present the official existing cases data for early March 2020 with Day 1 in our data set being 9 March 2020 and each increment of 1 being another day. Thus the first entry in the data pairs is the day and the second entry is the number of official existing cases of COVID-19 in Panama.  Day 25 is 2 April 2020, the day after the quarantine set in.

Here is a plot of the data.

data = {{1,1},{2,8},{3,14},{4,27},{5,36},{6,43},{7,55},{8,69},{9,86},{10,109},{11,137},{12,200},{13,245},{14,313},{15,345},{16,443},{17,558},{18,674},{19,786},{20,901},{21,989},{22,1075},{23,1181},{24,1317},{25,1475},{26,1673},{27,1801},{28,1988},{29,2100},{30,2249},{31,2528},{32,2752},{33, 2974},{34,3234},{35,3400},{36,3472},{37,3574},{38,3751},{39,4016},{40,4210},{41,4273},{42,4467},{43, 4658},{44,4821},{45,4992},{46,5166}};

qPlot =  ListPlot[{{{25,1475}}},PlotRange->{{0,data[[Length[data],1]]+1},{0,data[[Length[data],2]]*1.1}},AxesLabel->{"Day","Official Cases"},PlotStyle->{{Red,PointSize[.03]}}]

We will presume the total national population of Panama (based on 2017 data) to be 4,099,000 people.

K= 4099000;  (***2017 population of Panama ***)

Here is a plot of the data.

dataPlot = ListPlot[data,PlotRange->{{0,data[[Length[data],1]]+1},{0,data[[Length[data],2]]*1.1}},AxesLabel->{"Day","Official Cases"}]

dataPlot =Show[dataPlot,qPlot]

tend = data[[Length[data],1]]

CONTINUOUS MODELING

Model 1 Analysis - Continuous Model

Remove[r,A,B,CC,K,KK,SSE1,SSE2,SSE3]

ys[t_,r_,KK_] = y[t]/.DSolve[{y'[t]==r y[t](KK - y[t])/KK,y[0]==1},y[t],t][[1]]

SSE[r_,KK_] = Sum[(data[[i,2]] - ys[data[[i,1]],r,KK])^2,{i,1, Length[data]}];

solrKK=FindMinimum[SSE[r,KK],{r,.35},{KK,2800}]

SSE1Continuous = solrKK[[1]]

K1 = KK/.solrKK[[2]]

ysol1[t_] = ys[t,r,KK]/.solrKK[[2]]

Limit[ysol1[t],t->Infinity]

yopt1Plot = Plot[ysol1[t],{t,0,data[[Length[data],1]]+1},PlotRange->{{0,data[[Length[data],1]]+1},{0,data[[Length[data],2]]*1.1}},AxesLabel->{"Day","Official Cases"},PlotStyle->Green]

Show[yopt1Plot,dataPlot]

ContourPlot[SSE[r,KK],{r,0.2,.4},{KK,1500,6000},ContourShading->False,Contours->100]

Residual1Continuous = ListPlot[Table[{data[[i,1]],data[[i,2]]- ysol1[data[[i,1]]]},{i,1, Length[data]}],PlotStyle->Green,AxesLabel->{"Day", "Residual"}]

Model 2 Analysis - Continuous Model

Remove[r,A,B,CC,K,KK,SSE1,SSE2,SSE3]

ys[t_,A_,B_,KK_] = y[t]/.DSolve[{y'[t]==(A y[t]+B)(KK - y[t])/KK,y[0]==1},y[t],t][[1]]

SSE[A_,B_,KK_] = Sum[(data[[i,2]] - ys[data[[i,1]],A,B,KK])^2,{i,1, Length[data]}];

solABKK=FindMinimum[SSE[A,B,KK],{A,.35},{B,.1},{KK,2800}]

SSE2Continuous = solABKK[[1]]

K2 = KK/.solABKK[[2]]

ysol2[t_] = ys[t,A,B,KK]/.solABKK[[2]]

Limit[ysol2[t],t->Infinity]

Plot[ysol2[t],{t,0,365},PlotRange->{{0,180},{0,5000}},AxesLabel->{"Day","Official Cases"},PlotStyle->Purple]

yopt2Plot = Plot[ysol2[t],{t,0,data[[Length[data],1]]+1},PlotRange->{{0,data[[Length[data],1]]+1},{0,data[[Length[data],2]]*1.1}},AxesLabel->{"Day","Official Cases"},PlotStyle->Red]

Show[yopt2Plot,dataPlot]

Residual2Continuous = ListPlot[Table[{data[[i,1]],data[[i,2]]- ysol2[data[[i,1]]]},{i,1, Length[data]}],PlotStyle->Red,AxesLabel->{"Day", "Residual"}]

Model 3 Analysis - Continuous Model

Remove[r,A,B,CC,K,KK,SSE1,SSE2,SSE3]

solm=ParametricNDSolve[{y'[t]==(A y[t]^2+B y[t] + CC)(KK - y[t])/KK,y[0]==1},y,{t,0, data[[Length[data],1]]+1},{A,B,CC,KK}]

y[A,B,CC,KK][t]/.solm

Manipulate[Show[{dataPlot,Plot[y[A,B,CC,KK][t]/.solm,{t,0,100}, PlotRange->{{0,tend},{0,20000}}]}],{A,-1,1},{B,0,10},{CC,0,100},{KK,1,10000}]

Try (from Manipulate) A = -0.000001, B = .13, CC = 9.4, and KK = 6860.

SSE[A_,B_,CC_,KK_] = Sum[(data[[i,2]] - y[A,B,CC,KK][data[[i,1]]]/.solm )^2,{i,1,Length[data]}];

I am having trouble conditioning the FindMiniumum command to get success in Mathematica with this many parameters (5).

Try (from Manipulate) A = -0.000001, B = .13, CC = 9.4, and KK = 6860.

solABKK0=FindMinimum[SSE[A,B,CC,KK],{A,-0.000001},{B,.13},{CC,9.4},{KK,6860},MaxIterations->1000,WorkingPrecision->24,PrecisionGoal->10,AccuracyGoal->∞]

K3 = KK/.solABKK0[[2]]

SSE3Continuous=SSE[A,B,CC,KK]/.solABKK0[[2]]

ysol3[t_] = y[t]/.NDSolve[{y'[t]==(A y[t]^2+B y[t] + CC)(KK - y[t])/KK/.solABKK0[[2]],y[0]==1},y[t],{t,0, data[[Length[data],1]]+1}][[1]]

yopt3Plot = Plot[ysol3[t],{t,0,data[[Length[data],1]]+1},PlotRange->{{0,data[[Length[data],1]]+1},{0,data[[Length[data],2]]*1.1}},AxesLabel->{"Day","Official Cases"}]

Show[yopt3Plot,dataPlot]

Residual3Continuous = ListPlot[Table[{data[[i,1]],data[[i,2]]- ysol3[data[[i,1]]]},{i,1, Length[data]}],PlotStyle->Blue,AxesLabel->{"Day", "Residual"}]

CONTINUOUS VERSION
We now conduct an AIC test on this - Small AIC is better model, even if negative.

We will offer a way of determining which model is best. It is called the Akaike Information Criterion (AIC). A readable reference to AIC can be found in  Ledder, G. 2013. Mathematics for the Life Sciences. New York: Springer starting on page 138.

The Akaike Information Criterion, where k is number of parameters to fit, n is the number of observations, and RSS is the minimum residual sum of square errors (or just sum of square errors) is

                                              AIC = 2 (1+k)+n Log[    

NB: when SSE is small then the Log function cold be negative. 

Smaller (even negative) values of AIC represent a higher level of statistical support for the corresponding model. A difference of 2 in the AIC’s is needed to suggest a significant difference. 

The smallest AIC means best model - in terms of balancing ACCURACY (RSS) and COMPLEXITY (k).  

So after we complete all our model analyses (1) - (2) we will see if one of the models is best.  

We offer the parameters for and sum of square errors for the best parameters in Models (1), (2), and (#) above.

SSE1Continuous

SSE2Continuous

SSE3Continuous

So  Model (3) would appear to give best fit in terms of smallest sum of square errors, but let us apply AIC to see if one model is significantly better than all the others.

m = Length[data] (*** Number of data points. ***)

AIC[n_,RSS_,k_] = n Log[RSS/n] + 2 (k+1) (*** n data points, RSS = SSE, and k = number of parameters. ***)

We now compute the AIC for Models (1) , (2),  and (3) as AIC1, AIC2,  and AIC3 respectively.

AIC1 = AIC[m,SSE1Continuous,2]//Chop

AIC2 = AIC[m,SSE2Continuous,3]//Chop

AIC3= AIC[m,SSE3Continuous,4]//Chop

We see Model (2) and (3) in the Continuous Model are close, but since AIC3 is less by two or more we are likely to go with Model (3) in these analyses.  In  Model (3) we note that A is so small on the order of -0.0000171611 it effectively reduces Model (3) to Model (2),  making Models (2) and the best of the three as AIC1 is very high.

These values of long term cases or “carrying capacitcy” for the firuis in the population in Panama do vary widely with K3 being highest, although modek K2 and K3 are close to the same in other regards. We keep watching parameters in a second notebook, 1-123-T-Mma-PanamaCoronaVirus-Continuous-Predictive Potential-TeacherVersion .

{K1,K2,K3}

Conclusion - Continuous Model

Model (1) information in green, Model (2) information in  red, and Model (3) information in blue. Model (2) and Model (3) are essentially the same, visually at least and according to AIC above. Again, in Model (3) we note that A is so small (e.g., -0.0000191947) it effectively reduces Model (3) to Model (2), effectively making Models (2) and (3) the best of the three as AIC1 is very high, albeit K2 and K3 differ greatly.

Show[dataPlot,yopt1Plot,yopt2Plot,yopt3Plot,PlotLabel->"Model (1)-Green, Model (2)-Red, Model (3)-Blue"]

We plot the residuals for each model for comparisons, but we look at those for Model (2) first.

Residual2Continuous = ListPlot[Table[{data[[i,1]],data[[i,2]]- ysol2[data[[i,1]]]},{i,1, Length[data]}],PlotStyle->Red,AxesLabel->{"Day", "Residual"}]

The plot of our residuals should be normal, i.e. as many values on each side of 0 with a normal shape distribution. Here we see a histogram of the data for Models (2) and (3). This while there is pattern in plots of the residuals in each case - which is not a good sign of a good model.

Residual2HistData = Table[data[[i,2]]- ysol2[data[[i,1]]],{i,1, Length[data]}]

Histogram[Residual2HistData,7]

The plot of our residuals should be normal, i.e. as many values on each side of 0 with a normal shape distribution. Here we see a histogram of the data.

Residual3Continuous

Residual3HistData = Table[data[[i,2]]- ysol3[data[[i,1]]],{i,1, Length[data]}]

Histogram[Residual3HistData,7]

Show[Residual1Continuous,Residual2Continuous,Residual3Continuous]

There are tests for residuals to determine randomness and hence good fit. However, we will stop here and leave it to respective teachers to pursue this.