Compile bi-weekly confirmed and unconfirmed data by Sao Paulo district

In [28]:

# All confirmed cases, by districtconfirmed_data=lab_subset[lab_subset.CONCLUSION=='CONFIRMED']confirmed_counts=(confirmed_data.groupby(['ONSET','AGE_GROUP']).size().unstack().reindex(dates_index).fillna(0).sum())all_confirmed_cases=(confirmed_counts.reindex_axis(measles_data['AGE_GROUP'].unique()).fillna(0).values.astype(int))

# All clinical cases, by districtclinical_counts=(clinic_subset.groupby(['ONSET','AGE_GROUP']).size().unstack().reindex(dates_index).fillna(0).sum())all_clinical_cases=(clinical_counts.reindex_axis(measles_data['AGE_GROUP'].unique()).fillna(0).values.astype(int))

We will extend a simple SIR disease model, to account for confirmation status, which will be fit using MCMC.

This model fits the series of 2-week infection totals for each age group $a$ as a set of Poisson random variables:

$$Pr(I_{a}(t) | \lambda_a(t)) = \text{Poisson}(\lambda_a(t)) $$

Where the age-specific outbreak intensity at time $t$ is modeled as:

$$\lambda_a(t) = S_a(t-1) \frac{I(t-1)\mathbf{B}}{N_a} $$

where $S_a(t-1)$ is the number of susceptibles in age group $a$ in the previous time period, $I(t-1)$ an age-specific vector of the number of infected individuals in the previous time period, $\mathbf{B}$ a matrix of transmission coefficients (both within- and between-ages), and $N_a$ an estimate of the population of age-$a$ people in Sao Paulo.

The matrix $B$ was constructed from a scalar transmission parameter $\beta$, which was given a vague half-Cauchy prior (scale=25). This was used to represent within-age-group transmission, and hence placed on the diagonal of a square transmission matrix of size $A$. Off-diagonal elements, representing transmission between age groups were scaled by a decay parameter $\delta$ which was used to scale the transmission to adjacent groups according to:

$$\beta \delta^{|a-b|}$$

where a and b are indices of two age group. The resulting transmission matrix is parameterized as follows:

The basic reproductive number $R_0$ was calculated as the largest real-valued eigenvalue of the matrix $\mathbf{B}$. To impose a mild constraint on $R_0$, we applied a Gaussian prior distribution whose 1st and 99th quantiles are 8 and 24, respectively, a reasonable range for a measles outbreak:

In [35]:

frompymcimportMCMC,Matplot,AdaptiveMetropolis,MAP,Slicerfrompymcimport(Uniform,DiscreteUniform,Beta,Binomial,Normal,CompletedDirichlet,Pareto,Poisson,NegativeBinomial,negative_binomial_like,poisson_like,Lognormal,Exponential,binomial_like,TruncatedNormal,Binomial,Gamma,HalfCauchy,normal_like,MvNormalCov,Bernoulli,Uninformative,Multinomial,rmultinomial,rbinomial,Dirichlet,multinomial_like,uniform_like)frompymcimport(Lambda,observed,invlogit,deterministic,potential,stochastic,logit)defmeasles_model(obs_date,confirmation=True,migrant=False,constrain_R=True):''' Confirmation sub-model '''ifconfirmation:# Specify priors on age-specific meansage_classes=np.unique(age_index)μ=Normal("μ",mu=0,tau=0.0001,value=[0]*len(age_classes))σ=HalfCauchy('σ',0,25,value=1)var=σ**2ρ=Uniform('ρ',-1,1,value=0)# Build variance-covariance matrix with first-order correlation # among age classes@deterministicdefΣ(var=var,cor=ρ):I=np.eye(len(age_classes))*varE=np.diag(np.ones(len(age_classes)-1),k=-1)*var*correturnI+E+E.T# Age-specific probabilities of confirmation as multivariate normal # random variablesβ_age=MvNormalCov("β_age",mu=μ,C=Σ,value=[1]*len(age_classes))p_age=Lambda('p_age',lambdab=β_age:invlogit(b))@deterministic(trace=False)defp_confirm(b=β_age):returninvlogit(b[age_index])# Confirmation likelihoodlab_confirmed=Bernoulli('lab_confirmed',p=p_confirm,value=confirmed,observed=True)''' Truncate data at observation period '''obs_index=clinical_counts_2w.index<=obs_dateconfirmed_obs_t=confirmed_counts_2w[obs_index].values.astype(int)clinical_obs_t=clinical_counts_2w[obs_index].values.astype(int)n_periods,n_age_groups=confirmed_obs_t.shape# Index for observation date, used to index out values of interest # from the model.t_obs=obs_index.sum()-1ifconfirmation:@stochastic(dtype=int)defclinical_cases(value=(clinical_obs_t*0.5).astype(int),n=clinical_obs_t,p=p_age):# Binomial confirmation processreturnnp.sum([binomial_like(xi,ni,p)forxi,niinzip(value,n)])I=Lambda('I',lambdaclinical=clinical_cases:clinical+confirmed_obs_t.astype(int))assertI.value.shape==(t_obs+1,n_age_groups)age_dist_init=np.sum(I.value,0)/float(I.value.sum())else:I=confirmed_obs_t+clinical_obs_tassertI.shape==(t_obs+1,n_age_groups)age_dist_init=np.sum(I,0)/float(I.sum())# Calcuate age distribution from observed distribution of infecteds to date_age_dist=Dirichlet('_age_dist',np.ones(n_age_groups),value=age_dist_init[:-1]/age_dist_init.sum())age_dist=CompletedDirichlet('age_dist',_age_dist)@potentialdefage_dist_like(p=age_dist,I=I):returnmultinomial_like(I.sum(0),I.sum(),p)''' Disease transmission model '''# Transmission parameterβ=HalfCauchy('β',0,25,value=5)#[1]*n_age_groups) decay=Beta('decay',1,5,value=0.8)@deterministicdefB(b=β,d=decay):b=np.ones(n_age_groups)*bB=b*np.eye(n_age_groups)foriinrange(1,n_age_groups):B+=np.diag(np.ones(n_age_groups-i)*b[i:]*d**i,k=-i)B+=np.diag(np.ones(n_age_groups-i)*b[:-i]*d**i,k=i)returnB# Downsample annual series to observed age groupsdownsample=lambdax:np.array([x[s].mean()forsinage_slices])@deterministicdefR0(B=B):evs=np.linalg.eigvals(B)returnmax(evs[np.isreal(evs)])ifconstrain_R:@potentialdefconstrain_R0(R0=R0):# Weakly-informative prior to constrain R0 to be within the # typical measles rangereturnnormal_like(R0,16,3.4**-2)A=Lambda('A',lambdaR0=R0:75./(R0-1))lt_sum=downsample(np.tril(FOI_mat).sum(0)[::-1])natural_susc=Lambda('natural_susc',lambdaA=A:np.exp((-1/A)*lt_sum))@deterministicdefp_μ(natural_susc=natural_susc):returndownsample(sia_susc)*downsample(vacc_susc)*natural_suscifTrue:# Following Stan manual chapter 16.2λ_p=Pareto('λ_p',1.5,0.1,value=0.5)a=Lambda('a',lambdamu=p_μ,lam=λ_p:mu*lam,trace=False)b=Lambda('b',lambdamu=p_μ,lam=λ_p:(1-mu)*lam,trace=False)p_susceptible=Beta('p_susceptible',a,b,value=p_μ.value)else:p_σ=HalfCauchy('p_σ',0,5,value=1)m=Lambda('m',lambdap=p_μ:logit(p))θ_p=Normal('theta_p',m,p_σ)p_susceptible=Lambda('p_susceptible',lambdaθ_p=θ_p:invlogit(θ_p))# Estimated total initial susceptiblesS_0=Binomial('S_0',n=N_age.values.astype(int),p=p_susceptible)''' Model of migrant influx of susceptibles '''ifmigrant:# Data augmentation for migrant susceptiblesimaginary_migrants=1000000N_migrant=DiscreteUniform('N_migrant',0,imaginary_migrants,value=100000)μ_age=Uniform('μ_age',15,35,value=25)σ_age=Uniform('σ_age',1,10,value=5)M_age=Normal('M_age',μ_age,σ_age**-2,size=imaginary_migrants,trace=False)@deterministicdefM_0(M=M_age,N=N_migrant):# Take first N augmented susceptiblesM_real=M[:N]# Drop into age groupsM_group=pd.cut(M_real,[0,5,10,15,20,25,30,35,40,100],right=False)returnM_group.value_counts().valuesp_migrant=Lambda('p_migrant',lambdaM_0=M_0,S_0=S_0:M_0/(M_0+S_0))I_migrant=[Binomial('I_migrant_%i'%i,I[i],p_migrant)foriinrange(t_obs+1)]I_local=Lambda('I_local',lambdaI=I,I_m=I_migrant:np.array([Ii-ImiforIi,Imiinzip(I,I_m)]))S=Lambda('S',lambdaI=I,S_0=S_0,M_0=M_0:S_0+M_0-I.cumsum(0))S_local=Lambda('S_local',lambdaI=I_local,S_0=S_0:S_0-I.cumsum(0))else:# Remaining susceptibles at each 2-week periodS=Lambda('S',lambdaI=I,S_0=S_0:S_0-I.cumsum(axis=0))# Check shapeassertS.value.shape==(t_obs+1.,n_age_groups)# Susceptibles at time t, by ageS_age=Lambda('S_age',lambdaS=S:S[-1].astype(int))# Force of infection@deterministicdefλ(B=B,I=I,S=S):returnS*(I.dot(B)/N_age.values)# Check shapeassertλ.value.shape==(t_obs+1,n_age_groups)# FOI in observation periodλ_t=Lambda('λ_t',lambdalam=λ:lam[-1])# Effective reproductive numberR_t=Lambda('R_t',lambdaS=S,R0=R0:S.sum(1)*R0/N_age.sum())R_t_local=Lambda('R_t_local',lambdaS=S_local,R0=R0:S.sum(1)*R0/N_age.sum())# Poisson likelihood for observed cases@potentialdefnew_cases(I=I,lam=λ):returnpoisson_like(I[1:],lam[:-1])''' Vaccination targets '''@deterministicdefvacc_5(S=S_age):# Vaccination of 5 and underp=[0.95]+[0]*(n_age_groups-1)returnrbinomial(S,p)# Proportion of susceptibles vaccinatedpct_5=Lambda('pct_5',lambdaV=vacc_5,S=S_age:V.sum()/S.sum())@deterministicdefvacc_15(S=S_age):# Vaccination of 15 and underp=[0.95]*3+[0]*(n_age_groups-3)returnrbinomial(S,p)# Proportion of susceptibles vaccinatedpct_15=Lambda('pct_15',lambdaV=vacc_15,S=S_age:V.sum()/S.sum())@deterministicdefvacc_30(S=S_age):# Vaccination of 30 and underp=[0.95]*6+[0]*(n_age_groups-6)returnrbinomial(S,p)# Proportion of 30 and under susceptibles vaccinatedpct_30=Lambda('pct_30',lambdaV=vacc_30,S=S_age:V.sum()/S.sum())@deterministicdefvacc_adult(S=S_age):# Vaccination of adults under 30 (and young kids)p=[0.95,0,0,0,0.95,0.95]+[0]*(n_age_groups-6)returnrbinomial(S,p)# Proportion of adults under 30 (and young kids)pct_adult=Lambda('pct_adult',lambdaV=vacc_adult,S=S_age:V.sum()/S.sum())returnlocals()

/Users/fonnescj/anaconda3/lib/python3.5/site-packages/numpy/core/fromnumeric.py:225: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
return reshape(newshape, order=order)

/Users/fonnescj/anaconda3/lib/python3.5/site-packages/numpy/core/fromnumeric.py:225: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
return reshape(newshape, order=order)

/Users/fonnescj/anaconda3/lib/python3.5/site-packages/numpy/core/fromnumeric.py:225: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
return reshape(newshape, order=order)

/Users/fonnescj/anaconda3/lib/python3.5/site-packages/numpy/core/fromnumeric.py:225: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
return reshape(newshape, order=order)

/Users/fonnescj/anaconda3/lib/python3.5/site-packages/numpy/core/fromnumeric.py:225: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
return reshape(newshape, order=order)

/Users/fonnescj/anaconda3/lib/python3.5/site-packages/numpy/core/fromnumeric.py:225: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
return reshape(newshape, order=order)