The SMC Procedure

Although the stochastic susceptible-infected-recovered (SSIR) model with the constant transmission rate, as shown in the section Stochastic Susceptible-Infected-Recovered Model with a Constant Transmission Rate, is simple and easy to understand, it is not realistic because the transition rate is often time-varying as a result of the intervention, seasonality, and other complex reasons. This example illustrates how to define an SSIR model with a time-varying transmission rate in the SMC procedure and then estimate the parameters using the PMCMC method.

Consider the following SSIR model within the framework of the state space model:

Here , a Poisson distributed instance, denotes the number of new infections at time ; denotes the number of susceptible individuals at time t; is the number of new observed infections at time t; is the cumulative number of recovered infections until time t; is the number of documented active infectious individuals at time t; is the number of undocumented active infectious individuals at time t; is the transmission rate at time t; and is the effective reproduction number at time t. Moreover, N denotes the population; p is the diagnosis rate; D denotes the average time to removal after infection; and denotes the recovery rate. The natural logarithm of the transmission rate follows an AR(1) model with coefficients and , correlation parameter , and variance .

In this SSIR model, there are six parameters: D, , , p, , and ; and there are seven state variables: S, , , , , , and . Actually, the model can be represented without because no state variables or dependent variables depend on it. However, is worth investigating: indicates containment of the disease. Moreover, the equation of can be included in the equations of other state variables. For example, can be written as follows:

The equation can be written in the same way. The reason that is still kept as a state variable is because the estimated number or interval of new infections is also useful for public health organizations in a pandemic.

A scenario with a basic reproduction number that is a decreasing function of the time is considered in the data generating process. The function is defined as

where denotes the basic reproduction number at time t, denotes the decreasing rate, and denotes the long run value to which converges when t is large. Note that the preceding function form is unknown to the modeler; the AR(1) process for the natural logarithm of the transmission rate is used to approximate the true unknown function form of .

The following statements simulate 82 observations in which the first 75 observations are used as a training data set and the remaining 7 observations are used as a test data set:

*macro variables;
%let N=20000000; * 20 million people ;
%let ID_0 = 100;
%let IU_0 = 800;
%let RC_0 = 0;
%let D = 8.7;
%let R_0 = 3.1;
%let eta = 0.045;
%let p = 0.25;
data sir_data;
   t = 1;
   S = &N - &ID_0 - &IU_0;
   RC = &RC_0;
   IU = &IU_0;
   ID = &ID_0;
   R = &R_0;
   beta = R*1/&D;
   logbeta = log(beta);
   R_eff = R*S/&N;
   Y = ID;
   xi = beta*S*(IU + ID)/&N;
   sum = S + IU + ID + RC;
   output;
   call streaminit(1);
   do t = 2 to 82;
      Slag = S;
      S = S - xi;
      RC = RC + 1/&D*(IU + ID);
      IU = (1-1/&D)*IU + xi - Y;
      ID = (1-1/&D)*ID + Y;
      Y = RAND('BINOmial', &p, int(IU));
      R = &R_0*exp(-&eta*t) + 1-exp(-&eta*t);
      beta = R*1/&D;
      logbeta = log(beta);
      R_eff = R*S/&N;
      xi = beta*S*(IU + ID)/&N;
      sum = S + IU + ID + RC;
   output;
   end;
run;

The following DATA steps load the first 75 observations of the data set sir_data into the data set sir_train and then load that data set into a data table in your CAS session that is associated with the mycas CAS engine libref. The DATA step assumes that your CAS engine libref is named mycas, but you can substitute any appropriately defined CAS engine libref. The last 7 observations are loaded in to the data set sir_test.

data sir_train;
   set sir_data;
   if (t le 75) then output;
run;

data sir_test;
   set sir_data;
   if (t > 75) then output;
run;

data mycas.sir_train;
   set sir_train;
run;

The following statements use PROC SMC to define the preceding SSIR model and then use the PMCMC method to estimate the three parameters for the simulated data set. The posterior medians will be used as the parameters in the smoothing and forecasting process.

proc smc data=mycas.sir_train seed=1234 nthreads=32;
   id time = t;
   var Y;
   statevar S IU ID logBeta RC R_EFF xi;
   smooth nparticle=20480 alg=sir parm=median
      method = fixedlag(20) out(pct=2.5, 97.5)=mycas.sm;
   forecast nparticle=20480 alg=sir lead=7 parm=median
      out(pct=2.5, 97.5)=mycas.fc;
   parm D = 9.3, u1=-1.31, u2=0, p=0.3, rho=0.8, sigma_var=0.16;

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

   initstate ID ~ degenerate(100);
   initstate IU ~ uniform(0.5*100/p, 1.5*100/p);
   initstate RC ~ degenerate(0);
   initstate S ~ degenerate(&N-ID-IU);
   initstate logBeta ~ uniform(log(2.75*1/D), log(3.3*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 - Y.l1);
   state ID ~ degenerate((1-1/D)*ID.l1 + Y.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 Y ~ normal(p*IU, ifn(IU>1, 2*sqrt(p*(1-p)*IU), 2*sqrt(p*(1-p))));

   learn nparticle=1600 algorithm=sir method=PMCMC(nsample=20000
      nchain=8 nbi=5000 thin=5 seed=1
      outpost=mycas.out samset samsum
      sampler=rwm(ntu=200)
      statistics=all diagnostics=all);
run;

The summary tables of the model information are shown in Output 23.2.1.

The summary tables of the posterior samples are shown in Output 23.2.2.

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

The following statements load the smoothing estimates for and into the Work library and then sort each data set by the time ID:

data smo_s4; set mycas.sm(where=(VariableID=4 and TimeID>=1)); run;
proc sort data=smo_s4; by timeId; run;

data smo_s6; set mycas.sm(where=(VariableID=6 and TimeID>=1)); run;
proc sort data=smo_s6; by timeId; run;

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

data combine_logbeta;
   set smo_s4 (keep=P_1 P_2 Q2);
   set sir_train (keep=t logbeta);
run;

proc sgplot data=combine_logbeta;
   band x=t upper=P_2 lower=P_1 /transparency=.3;
   series x=t y=logbeta /legendlabel='Actual Value' lineattrs=(color=red);
   series x=t y=Q2 /legendlabel='Smoothing Median'
      lineattrs=(color=blue pattern=shortdash);
   yaxis label='logbeta, Smoothing';
   xaxis integer label='Day';
run;

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

data combine_REFF;
   set smo_s6 (keep=P_1 P_2 Q2);
   set sir_train (keep=t R_EFF);
run;

proc sgplot data=combine_REFF;
   band x=t upper=P_2 lower=P_1 /transparency=.3;
   series x=t y=R_EFF /legendlabel='Actual Value' lineattrs=(color=red);
   series x=t y=Q2 /legendlabel='Smoothing Median'
      lineattrs=(color=blue pattern=shortdash);
   yaxis label='R_EFF, Smoothing';
   xaxis integer label='Day';
run;

The following statements load the forecasting estimates for , , and Y into the Work library and then sort each data set by the time ID:

data fc_6; set mycas.fc(where=(VariableID=6)); run;
proc sort data=fc_6; by _t; run;

data fc_7; set mycas.fc(where=(VariableID=7)); run;
proc sort data=fc_7; by _t; run;

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

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

data combine_REFF_forecast;
   set fc_6 (keep=P_1 P_2 Q2);
   set sir_test (keep=t R_EFF);
run;

proc sgplot data=combine_REFF_forecast;
   band x=t upper=P_2 lower=P_1 /transparency=.3;
   scatter x=t y=R_EFF /legendlabel='Actual Value'
      markerattrs=(color=red size=6 symbol=circlefilled);
   series x=t y=Q2 /legendlabel='Forecasting Median'
      markers markerattrs=(color=blue size=6 symbol=triangle)
      lineattrs=(color=blue pattern=solid);
   yaxis label='R_EFF, Forecasting';
   xaxis integer label='Day';
run;

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

data combine_xi_forecast;
   set fc_7 (keep=P_1 P_2 Q2);
   set sir_test (keep=t xi);
run;

proc sgplot data=combine_xi_forecast;
   band x=t upper=P_2 lower=P_1 /transparency=.3;
   scatter x=t y=xi /legendlabel='Actual Value'
      markerattrs=(color=red size=6 symbol=circlefilled);
   series x=t y=Q2 /legendlabel='Forecasting Median'
      markers markerattrs=(color=blue size=6 symbol=triangle)
      lineattrs=(color=blue pattern=solid);
   yaxis label='xi, Forecasting';
   xaxis integer label='Day';
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 Y. Then they draw a plot to compare the forecasting estimates to the actual values at each time step. The plot is shown in Output 23.2.8, 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 sir_test (keep=t Y);
run;

proc sgplot data=combine_Y_forecast;
   band x=t upper=P_2 lower=P_1 /transparency=.3;
   scatter x=t y=Y /legendlabel='Actual Value'
      markerattrs=(color=red size=6 symbol=circlefilled);
   series x=t y=Q2 /legendlabel='Forecasting Median'
      markers markerattrs=(color=blue size=6 symbol=triangle)
      lineattrs=(color=blue pattern=solid);
   yaxis label='Y, Forecasting';
   xaxis integer label='Day';
run;

As shown in this example, by introducing the AR(1) process for the transition rate, this simple and powerful model successfully captures the dynamics in the epidemic spread, even when the true DGP follows a different stochastic process. This model is applied to the real-data example in the section Estimating the Effective Reproduction Number of COVID-19 Cases in Pennsylvania.