SMC Procedure

This example illustrates how to estimate the time-varying effective reproduction number of COVID-19 using the stochastic susceptible-infected-recovered (SSIR) model introduced in the section Stochastic Susceptible-Infected-Recovered Model with a Time-Varying Transmission Rate. The daily confirmed number of COVID-19 cases in Pennsylvania reported by the Center for Systems Science and Engineering at Johns Hopkins University (JHU CSSE) (Dong, Du, and Gardner 2020) is used.

The following DATA step inputs the number of new COVID-19 cases confirmed in Pennsylvania into the data set penn_state_data:

data penn_state_data;
   input date: anydtdte. cases newcase;
   format date DATE.;
   label date="Date" cases="cumulative cases" newcase="new case";
datalines;
22-Mar-20  509  110
23-Mar-20  698  189
24-Mar-20  946  248
25-Mar-20  1260  314
26-Mar-20  1795  535
27-Mar-20  2345  550

   ... more lines ...   

25-Sep-20  159014  833
26-Sep-20  159933  919
27-Sep-20  160744  811
28-Sep-20  161594  850
29-Sep-20  163180  1586
30-Sep-20  164099  919
1-Oct-20   165044  945
2-Oct-20   166238  1194
;

The following statements plot the number of new cases confirmed in Pennsylvania from March 22 to September 28, 2020. The output is shown in Output 27.3.1.

proc sgplot data=penn_state_data;
   series x=date y=newcase /
      markers markerattrs=(color=red size=5 symbol=circlefilled)
      lineattrs=(color=red);
   yaxis label='New Case';
   xaxis label='Date';
run;

The data start from March 22, 2020, because it is the first day when Pennsylvania had more than 100 new cases and more than 500 cumulative cases. The transition and measurement distributions for the SSIR model are exactly the same as the ones in the section Stochastic Susceptible-Infected-Recovered Model with a Time-Varying Transmission Rate. However, the initialization distribution must be adjusted because the first day of the data (March 22, 2020) is not the day when the first case was reported in Pennsylvania.

Consider the following initialization distribution:

Here all the denotations are defined in the same way as in the section Stochastic Susceptible-Infected-Recovered Model with a Time-Varying Transmission Rate. The cumulative number of reported cases until March 22, 2020, is 509, and the number of new cases that day is 110. Thus, it is reasonable to set initialization distributions for and as earlier. In the initialization distribution for , can be seen as the expected value of the number of undocumented infectious individuals before that day; 110 is subtracted because the 110 cases are confirmed that day.

The following statements add a logical time ID ("t") and split the data set such that the first 181 observations are used as the training set and the last 14 observations are used as the test set:

data penn_state_data;
   set penn_state_data;
   t = _N_;
RUN;

data penn_state_train;
   set penn_state_data;
   if (t le 181) then output;
run;

data penn_state_test;
   set penn_state_data;
   if (t > 181) then output;
run;

The following DATA step loads the training set into a data table in your session that is associated with the mylib libref. The DATA step assumes that your libref is named mylib, but you can substitute any appropriately defined libref.

data mylib.train;
   set penn_state_train;
run;

The following statements use PROC SMC to define the SSIR model and then use the PMCMC method to estimate the three parameters for the simulated data set. Moreover, the population of Pennsylvania is set as a macro variable. The median of the posterior parameter samples will be used as the parameters in the smoothing and forecasting process.

%let N = 12791692; *Pennsylvania population;
proc smc data=mylib.train seed = 1234 nthreads = 32;
   id time = date;
   var newcase;
   statevar S IU ID logBeta RC R_EFF xi;
   smooth nparticle=24000 alg=sir parm=median
      method = fixedlag(20) out(pct=2.5, 97.5)=mylib.sm;
   forecast nparticle=24000 alg=sir lead=14 parm=median
      out(pct=2.5, 97.5)=mylib.fc;
   parm D = 9.3, u1=-1.31, u2=0, p=0.25, rho=0.8, sigma_var=0.16
      /lower=(5.75 -2 -2 0.1 0.35 0.05) upper=(14 2 2 0.5 0.95 2);

   prior D ~ gamma(325.5, 0.0285714286);
   prior sigma_var ~ igamma(11,1);
   prior rho ~ beta(4,1);
   prior p ~ uniform(0.1, 0.5);
   prior u1 ~ normal(-1.31, 0.3);
   prior u2 ~ normal(0,1);

   initstate ID ~ uniform(450, 550);
   initstate IU ~ uniform(0.98*110/p-110, 3*110/p-110);
   initstate RC ~ degenerate(0);
   initstate S ~ degenerate(&N-ID-IU);
   initstate logBeta ~ uniform(log(2.75*1/D), log(3.43*1/D));
   initstate R_EFF ~ degenerate(exp(logBeta)*S*D/&N);
   initstate xi ~ degenerate(exp(logBeta)*S*(IU + ID)/&N);

   state S ~ degenerate(S.l1 - xi.l1);
   state RC ~ degenerate(RC.l1 + 1/D*(IU.l1+ID.l1));
   state IU ~ degenerate((1-1/D)*IU.l1 + xi.l1 - newcase.l1);
   state ID ~ degenerate((1-1/D)*ID.l1 + newcase.l1);
   state logBeta ~ normal(u1 + u2*_time + rho * (logBeta.l1-u1-u2*(_time-1)),
      sqrt(sigma_var*(1-rho**2)));
   state R_EFF ~ degenerate(exp(logBeta)*S*D/&N);
   state xi ~ degenerate(exp(logBeta)*S*(IU + ID)/&N);

   model newcase ~ normal(p*IU, 2*sqrt(p*(1-p)*ifn(IU>1, IU, 1)));

   learn nparticle=8000 algorithm=sir method=PMCMC(nsample=5000
      nchain=8 nbi=4000 thin=4 seed=1
      outpost=mylib.out samset samsum
      sampler=rwm(ntu=200)
      statistics=all diagnostics=all);
run;

The summary tables of the model information are shown in Figure 29.

The summary tables of the posterior samples are shown in Figure 30.

The correlation and covariance matrices of the parameter samples that are generated by the PMCMC method are shown in Figure 31.

The following statements load the smoothing result for the effective reproduction number (VariableID = 6) into the Work library and then sort that data set by the logical time ID:

data smo_s6; set mylib.sm(where=(VariableID=6)); run;
proc sort data=smo_s6; by _t; run;

The following statements create a relevant data table, combine_REFF, in the Work library to record the smoothing estimates and actual values of the effective reproduction number. Then they draw a plot to compare the smoothing estimates to the actual values at each time step. The output is shown in Figure 32, where the shadow band in the plot indicates the range between the 2.5th and 97.5th percentiles of the smoothing estimates and the horizontal reference line denotes the line .

data combine_REFF;
   set smo_s6 (keep=P_1 P_2 Q2 _t);
   set penn_state_train (keep=date);
run;

proc sgplot data=combine_REFF;
   band x=date upper=P_2 lower=P_1 /transparency=.3;
   series x=date y=Q2 /legendlabel='Smoothing Median'
      lineattrs=(color=blue pattern=shortdash);
   yaxis label='R_EFF, Smoothing';
   xaxis label='Date';
   refline 1 / axis=y label='1' lineattrs=(thickness=3);
run;

The following statements load the forecasting result for each state variable into the Work library and then sort each data set by the time ID:

data fc_8; set mylib.fc(where=(VariableID=8)); run;
proc sort data=fc_8; by _t; run;

The following statements create a relevant data table, combine_Y_forecast, in the Work library to record the forecasting estimates and actual values of new cases. Then they draw a plot to compare the forecasting estimates to the actual values at each time step. The plot is shown in Figure 33, where the shadow band indicates the range between the 2.5th and 97.5th percentiles of the forecasting estimates.

data combine_Y_forecast;
   set fc_8 (keep=P_1 P_2 Q2);
   set penn_state_test (keep=date t newcase);
run;

proc sgplot data=combine_Y_forecast;
   band x=date upper=P_2 lower=P_1 /transparency=.3;
   series x=date y=newcase /legendlabel='Actual Value'
      markers markerattrs=(color=red size=5 symbol=circlefilled)
      lineattrs=(color=red);
   series x=date y=Q2 /legendlabel='Forecasting Median'
      markers markerattrs=(color=blue size=5 symbol=trianglefilled)
      lineattrs=(color=blue pattern=shortdash);
   yaxis label='New Case, Forecasting';
   xaxis label='Date';
run;