%load_ext watermark
%watermark -v -m -p numpy,pandas,pymc -g
CPython 3.4.2 IPython 4.0.0-dev numpy 1.10.0.dev-fb898ce pandas 0.15.2-50-g4aa0e0a pymc 2.3.4 compiler : GCC 4.2.1 Compatible Apple LLVM 6.0 (clang-600.0.54) system : Darwin release : 14.3.0 machine : x86_64 processor : i386 CPU cores : 4 interpreter: 64bit Git hash :
%matplotlib inline
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
import pymc as pm
from dateutil.parser import parse
from datetime import datetime
data_dir = './data/'
_adult_data = pd.read_csv(data_dir+'adult_data.csv')
_adult_data.head(5)
/usr/local/lib/python3.4/site-packages/pandas/io/parsers.py:1159: DtypeWarning: Columns (40,111,112,132,188,193) have mixed types. Specify dtype option on import or set low_memory=False. data = self._reader.read(nrows)
COUNTY | STATEID | BIRTH | AGE | UNIT | SEX | RACE | ETHNIC | OUTCOME | BACTEREM | ... | regioneast | agegroup | _clinsynd | vaccineserotype | PENR | PENS | LEVR | LEVS | ERYTHR | ERYTHS | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | WILLIAMSON | TNF6092 | 05nov1972 | 37 | 3 | 1 | 1 | 9 | 1 | 0 | ... | 0 | 3 | 3 | 2 | 1 | 0 | 0 | 1 | 1 | 0 |
1 | SHELBY | TND6006 | 27jun1927 | 80 | 3 | 2 | 1 | 2 | 1 | 0 | ... | 0 | 4 | 3 | 5 | 0 | 1 | 0 | 1 | 0 | 1 |
2 | DAVIDSON | TN43378 | 15apr1958 | 45 | 3 | 1 | NaN | 9 | 1 | 1 | ... | 0 | 3 | 1 | 3 | 0 | 1 | 0 | 1 | 0 | 1 |
3 | SHELBY | TNC6014 | 28jun1950 | 56 | 3 | 2 | 1 | 9 | 1 | 0 | ... | 0 | 3 | 9 | NaN | NaN | NaN | NaN | NaN | NaN | NaN |
4 | RUTHERFORD | TNG6394 | 23sep1953 | 58 | 3 | 2 | 4 | 2 | 1 | 0 | ... | 0 | 3 | 3 | 2 | 1 | 0 | 0 | 1 | 1 | 0 |
5 rows × 269 columns
_adult_2013 = pd.read_csv(data_dir+'2013 adult merged CRF and lab info.csv')
Recode RACE variable
_adult_2013['RACE'] = 1
_adult_2013.loc[_adult_2013.BLACK==1, 'RACE'] = 2
_adult_2013.loc[_adult_2013.AMERIND==1, 'RACE'] = 3
_adult_2013.loc[_adult_2013.ASIAN==1, 'RACE'] = 4
_adult_2013.loc[_adult_2013.PACIFIC==1, 'RACE'] = 5
_adult_2013.loc[_adult_2013.UNKRACE==1, 'RACE'] = 9
_adult_data = _adult_data.append(_adult_2013, ignore_index=True)
keep_columns = ['COUNTY','STATEID','RACE','calcage','agegroup','regioneast','CULT','year',
'pneumo7','pneumo13','pneumo23','pneumononPCV','PENR','LEVR','ERYTHR']
adult_data = _adult_data[keep_columns].copy()
def parse_date(d):
if d.find('-')==-1:
return datetime.strptime(d, '%d%b%Y')
return datetime.strptime(d, '%d-%b-%Y')
adult_data['culture'] = adult_data.CULT.apply(parse_date)
#adult_data['culture'] = [datetime.strptime(d, '%d%b%Y') for d in adult_data.CULT]
adult_data['over64'] = (adult_data.agegroup==4).astype(int)
adult_data.agegroup.value_counts()
3 3548 4 1933 dtype: int64
Obtain subset of interest
#data_subset = adult_data[(adult_data.over64==0) & (adult_data.regioneast==1)]
data_subset = adult_data.copy()
covnames = ['over64', 'regioneast']
# covnames = ['over64','regioneast','pneumo7',
# 'pneumo13','pneumo23','pneumononPCV','PENR','LEVR','ERYTHR'
adult_resampled = data_subset.groupby([pd.Grouper(freq='AS', key='culture'), 'over64', 'regioneast'])['CULT'].count()
adult_resampled = adult_resampled.reset_index().rename(columns={'CULT':'count', 'culture':'year'})
adult_resampled['year'] = [d.year-2001 for d in adult_resampled.year]
adult_resampled
year | over64 | regioneast | count | |
---|---|---|---|---|
0 | 0 | 0 | 0 | 207 |
1 | 0 | 0 | 1 | 66 |
2 | 0 | 1 | 0 | 119 |
3 | 0 | 1 | 1 | 60 |
4 | 1 | 0 | 0 | 188 |
5 | 1 | 0 | 1 | 53 |
6 | 1 | 1 | 0 | 88 |
7 | 1 | 1 | 1 | 40 |
8 | 2 | 0 | 0 | 175 |
9 | 2 | 0 | 1 | 60 |
10 | 2 | 1 | 0 | 110 |
11 | 2 | 1 | 1 | 39 |
12 | 3 | 0 | 0 | 161 |
13 | 3 | 0 | 1 | 50 |
14 | 3 | 1 | 0 | 87 |
15 | 3 | 1 | 1 | 39 |
16 | 4 | 0 | 0 | 221 |
17 | 4 | 0 | 1 | 91 |
18 | 4 | 1 | 0 | 90 |
19 | 4 | 1 | 1 | 53 |
20 | 5 | 0 | 0 | 234 |
21 | 5 | 0 | 1 | 79 |
22 | 5 | 1 | 0 | 100 |
23 | 5 | 1 | 1 | 61 |
24 | 6 | 0 | 0 | 219 |
25 | 6 | 0 | 1 | 81 |
26 | 6 | 1 | 0 | 99 |
27 | 6 | 1 | 1 | 55 |
28 | 7 | 0 | 0 | 232 |
29 | 7 | 0 | 1 | 110 |
30 | 7 | 1 | 0 | 117 |
31 | 7 | 1 | 1 | 61 |
32 | 8 | 0 | 0 | 214 |
33 | 8 | 0 | 1 | 93 |
34 | 8 | 1 | 0 | 95 |
35 | 8 | 1 | 1 | 43 |
36 | 9 | 0 | 0 | 198 |
37 | 9 | 0 | 1 | 81 |
38 | 9 | 1 | 0 | 106 |
39 | 9 | 1 | 1 | 48 |
40 | 10 | 0 | 0 | 186 |
41 | 10 | 0 | 1 | 82 |
42 | 10 | 1 | 0 | 104 |
43 | 10 | 1 | 1 | 61 |
44 | 11 | 0 | 0 | 170 |
45 | 11 | 0 | 1 | 71 |
46 | 11 | 1 | 0 | 84 |
47 | 11 | 1 | 1 | 43 |
48 | 12 | 0 | 0 | 151 |
49 | 12 | 0 | 1 | 75 |
50 | 12 | 1 | 0 | 92 |
51 | 12 | 1 | 1 | 39 |
# adult_data = pd.DataFrame(np.c_[[a.year-2001 for a in adult_resampled.index], adult_resampled.values],
# columns=['year', 'count'])
adult_resampled.tail()
year | over64 | regioneast | count | |
---|---|---|---|---|
47 | 11 | 1 | 1 | 43 |
48 | 12 | 0 | 0 | 151 |
49 | 12 | 0 | 1 | 75 |
50 | 12 | 1 | 0 | 92 |
51 | 12 | 1 | 1 | 39 |
denom = pd.read_csv(data_dir+'adult_denominator.csv')
denom['over64'] = (denom._newagegroupadult!='18to65').astype(int)
denom['year'] = denom.year - 1
denom.rename(columns={'_newagegroupadultpop':'n'}, inplace=True)
denom = denom[denom._newagegroupadult!='total']
denom.head()
year | _newagegroupadult | regioneast | n | over64 | |
---|---|---|---|---|---|
0 | 0 | 18to65 | 0 | 1377668.2 | 0 |
1 | 1 | 18to65 | 0 | 1393999.8 | 0 |
2 | 2 | 18to65 | 0 | 1410903.8 | 0 |
3 | 3 | 18to65 | 0 | 1432392.4 | 0 |
4 | 4 | 18to65 | 0 | 1456641.2 | 0 |
Add 2013 census data
denom13_raw = pd.read_csv(data_dir+'2013 adult census data.csv')
denom13 = denom13_raw.drop('COUNTY', axis=1).groupby(['_newagegroupadult', 'year', 'regioneast']).sum().reset_index()
denom13['over64'] = (denom13._newagegroupadult!='18to65').astype(int)
denom13['year'] = denom.year.max() + 1
denom13.rename(columns={'_newagegroupadultpop':'n'}, inplace=True)
denom13 = denom13[denom13._newagegroupadult!='total'].drop('_newagegroupadult', axis=1)
denom13.head()
year | regioneast | n | over64 | |
---|---|---|---|---|
0 | 12 | 0 | 1617941.2 | 0 |
1 | 12 | 1 | 503405.0 | 0 |
2 | 12 | 0 | 287483.0 | 1 |
3 | 12 | 1 | 117198.0 | 1 |
# denom_subset = denom.ix[(denom.over64==0) & (denom.regioneast==1), ['year', 'n']]
denom_subset = denom[['year', 'n', 'regioneast', 'over64']].append(denom13, ignore_index=True)
data_merged = pd.merge(denom_subset, adult_resampled, on=['year', 'regioneast', 'over64'])
data_merged = data_merged[['year', 'n', 'count', 'regioneast', 'over64']]
data_merged['post_2010'] = (data_merged.year > 9).astype(int)
count = data_merged['count'].values
year = data_merged['year'].values
n = data_merged['n'].values
X = data_merged[['regioneast', 'over64']].values
post_2010 = data_merged.post_2010
nyears = data_merged.year.unique().shape[0]
nyears
13
X.shape
(52, 2)
data_merged.groupby(['regioneast', 'over64']).plot(x='year', y='count', ax=plt.subplot(), legend=False);
data_merged.head()
year | n | count | regioneast | over64 | post_2010 | |
---|---|---|---|---|---|---|
0 | 0 | 1377668.2 | 207 | 0 | 0 | 0 |
1 | 1 | 1393999.8 | 188 | 0 | 0 | 0 |
2 | 2 | 1410903.8 | 175 | 0 | 0 | 0 |
3 | 3 | 1432392.4 | 161 | 0 | 0 | 0 |
4 | 4 | 1456641.2 | 221 | 0 | 0 | 0 |
data_merged.shape
(52, 6)
Model year effect as first-order autocorrelated random effect
pcv13 = pm.Normal('pcv13', 0, 0.01, value=0)
mu_alpha = pm.Normal('mu_alpha', 0, 0.01, value=0)
sigma_alpha = pm.Uniform('sigma_alpha', 0, 100, value=10)
rho = pm.Uniform('rho', -1, 1, value=0)
mu = pm.Lambda('mu', lambda mu=mu_alpha, d=pcv13: [mu]*nyears + d*post_2010[:nyears])
off_diag = np.eye(nyears, k=1)
Sigma = pm.Lambda('Sigma', lambda s=sigma_alpha, r=rho: (np.eye(nyears) + off_diag*r + off_diag.T*r)*(s**2))
alpha = pm.MvNormalCov('alpha', mu, Sigma, value=[0]*nyears)
mu_alpha = pm.Normal('mu_alpha', 0, 0.01, value=0)
sigma_alpha = pm.Uniform('sigma_alpha', 0, 100, value=10)
alpha = pm.Normal('alpha', mu_alpha, sigma_alpha**-2, value=[0]*nyears)
Covariate effects
beta = pm.Normal('beta', 0, 0.001, value=[0]*X.shape[1])
Poisson regression expected value
rate = pm.Lambda('rate', lambda alpha=alpha, beta=beta: np.exp(alpha[year] + X.dot(beta)))
rate_west_under65 = pm.Lambda('rate_west_under65', lambda r=rate: r[:12])
rate_west_over64 = pm.Lambda('rate_west_over64', lambda r=rate: r[12:24])
rate_east_under65 = pm.Lambda('rate_east_under65', lambda r=rate: r[24:36])
rate_east_over64 = pm.Lambda('rate_east_over64', lambda r=rate: r[36:])
theta = pm.Lambda('theta', lambda rate=rate: rate * n)
Likelihood function
cases = pm.Poisson('cases', theta, value=count, observed=True)
# cases = pm.Binomial('cases', n.astype(int), rate, value=count, observed=True)
Posterior predictive distribution for all observations (for goodness-of-fit)
cases_pred = pm.Poisson('cases_pred', theta, value=count)
# cases_pred = pm.Binomial('cases_pred', n.astype(int), rate, value=count)
Run MCMC sampler
M = pm.MCMC(locals())
M.sample(50000, 40000)
[-----------------100%-----------------] 50000 of 50000 complete in 45.2 sec
Posterior mean
pm.Matplot.summary_plot(alpha)
Could not calculate Gelman-Rubin statistics. Requires multiple chains of equal length.
pm.Matplot.plot(beta)
Plotting beta_0 Plotting beta_1
pm.Matplot.plot(pcv13)
Plotting pcv13
beta.summary()
beta: Mean SD MC Error 95% HPD interval ------------------------------------------------------------------ 0.187 0.029 0.001 [ 0.132 0.244] 1.141 0.029 0.002 [ 1.085 1.198] Posterior quantiles: 2.5 25 50 75 97.5 |---------------|===============|===============|---------------| 0.128 0.169 0.188 0.207 0.242 1.084 1.122 1.14 1.162 1.198
pm.Matplot.plot(mu_alpha)
Plotting mu_alpha
pm.Matplot.summary_plot(rate)
Could not calculate Gelman-Rubin statistics. Requires multiple chains of equal length.
Goodness of fit plots (posterior predictive checks)
#pm.Matplot.gof_plot(cases_pred, count)
data_subset = adult_data.copy()
data_subset.head()
COUNTY | STATEID | RACE | calcage | agegroup | regioneast | CULT | year | pneumo7 | pneumo13 | pneumo23 | pneumononPCV | PENR | LEVR | ERYTHR | culture | over64 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | WILLIAMSON | TNF6092 | 1 | 37.275837 | 3 | 0 | 14feb2010 | 10 | 0 | 1 | 0 | 0 | 1 | 0 | 1 | 2010-02-14 | 0 |
1 | SHELBY | TND6006 | 1 | 80.528404 | 4 | 0 | 06jan2008 | 8 | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 2008-01-06 | 1 |
2 | DAVIDSON | TN43378 | NaN | 45.889118 | 3 | 0 | 05mar2004 | 4 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 2004-03-05 | 0 |
3 | SHELBY | TNC6014 | 1 | 56.594112 | 3 | 0 | 31jan2007 | 7 | NaN | NaN | NaN | NaN | NaN | NaN | NaN | 2007-01-31 | 0 |
4 | RUTHERFORD | TNG6394 | 4 | 58.028748 | 3 | 0 | 04oct2011 | 11 | 0 | 1 | 0 | 0 | 1 | 0 | 1 | 2011-10-04 | 0 |
data_grouped = data_subset.groupby([pd.Grouper(freq='AS', key='culture'), 'over64', 'regioneast'])['LEVR']
# data_grouped = data_subset.groupby([pd.Grouper(freq='AS', key='culture'), 'over64'])['LEVR']
adult_resampled = pd.DataFrame(dict(count = data_grouped.sum())).reset_index().rename(columns={'culture':'year'})
adult_resampled['year'] = [d.year-2001 for d in adult_resampled.year]
adult_resampled['post_2010'] = (adult_resampled.year > 9).astype(int)
adult_resampled
year | over64 | regioneast | count | post_2010 | |
---|---|---|---|---|---|
0 | 0 | 0 | 0 | 0 | 0 |
1 | 0 | 0 | 1 | 1 | 0 |
2 | 0 | 1 | 0 | 3 | 0 |
3 | 0 | 1 | 1 | 0 | 0 |
4 | 1 | 0 | 0 | 0 | 0 |
5 | 1 | 0 | 1 | 0 | 0 |
6 | 1 | 1 | 0 | 0 | 0 |
7 | 1 | 1 | 1 | 0 | 0 |
8 | 2 | 0 | 0 | 1 | 0 |
9 | 2 | 0 | 1 | 0 | 0 |
10 | 2 | 1 | 0 | 3 | 0 |
11 | 2 | 1 | 1 | 0 | 0 |
12 | 3 | 0 | 0 | 0 | 0 |
13 | 3 | 0 | 1 | 0 | 0 |
14 | 3 | 1 | 0 | 0 | 0 |
15 | 3 | 1 | 1 | 0 | 0 |
16 | 4 | 0 | 0 | 2 | 0 |
17 | 4 | 0 | 1 | 0 | 0 |
18 | 4 | 1 | 0 | 1 | 0 |
19 | 4 | 1 | 1 | 0 | 0 |
20 | 5 | 0 | 0 | 0 | 0 |
21 | 5 | 0 | 1 | 0 | 0 |
22 | 5 | 1 | 0 | 1 | 0 |
23 | 5 | 1 | 1 | 0 | 0 |
24 | 6 | 0 | 0 | 1 | 0 |
25 | 6 | 0 | 1 | 0 | 0 |
26 | 6 | 1 | 0 | 2 | 0 |
27 | 6 | 1 | 1 | 0 | 0 |
28 | 7 | 0 | 0 | 3 | 0 |
29 | 7 | 0 | 1 | 0 | 0 |
30 | 7 | 1 | 0 | 0 | 0 |
31 | 7 | 1 | 1 | 0 | 0 |
32 | 8 | 0 | 0 | 0 | 0 |
33 | 8 | 0 | 1 | 0 | 0 |
34 | 8 | 1 | 0 | 1 | 0 |
35 | 8 | 1 | 1 | 0 | 0 |
36 | 9 | 0 | 0 | 0 | 0 |
37 | 9 | 0 | 1 | 0 | 0 |
38 | 9 | 1 | 0 | 0 | 0 |
39 | 9 | 1 | 1 | 0 | 0 |
40 | 10 | 0 | 0 | 0 | 1 |
41 | 10 | 0 | 1 | 0 | 1 |
42 | 10 | 1 | 0 | 2 | 1 |
43 | 10 | 1 | 1 | 0 | 1 |
44 | 11 | 0 | 0 | 2 | 1 |
45 | 11 | 0 | 1 | 0 | 1 |
46 | 11 | 1 | 0 | 0 | 1 |
47 | 11 | 1 | 1 | 1 | 1 |
48 | 12 | 0 | 0 | 1 | 1 |
49 | 12 | 0 | 1 | 0 | 1 |
50 | 12 | 1 | 0 | 0 | 1 |
51 | 12 | 1 | 1 | 0 | 1 |
data_merged = pd.merge(denom_subset, adult_resampled, on=['year', 'regioneast', 'over64'])
data_merged = data_merged[['year', 'n', 'count', 'regioneast', 'over64']]
data_merged['post_2010'] = (data_merged.year > 9).astype(int)
data_merged.head()
year | n | count | regioneast | over64 | post_2010 | |
---|---|---|---|---|---|---|
0 | 0 | 1377668.2 | 0 | 0 | 0 | 0 |
1 | 1 | 1393999.8 | 0 | 0 | 0 | 0 |
2 | 2 | 1410903.8 | 1 | 0 | 0 | 0 |
3 | 3 | 1432392.4 | 0 | 0 | 0 | 0 |
4 | 4 | 1456641.2 | 2 | 0 | 0 | 0 |
def resample_data(variable, regional=True):
groupby_ind = [pd.Grouper(freq='AS', key='culture'), 'over64']
if regional:
groupby_ind += ['regioneast']
data_grouped = data_subset.groupby(groupby_ind)[variable]
adult_resampled = pd.DataFrame(dict(n = data_grouped.count(),
count = data_grouped.sum())).reset_index().rename(columns={'culture':'year'})
adult_resampled['year'] = [d.year-2001 for d in adult_resampled.year]
adult_resampled['post_2010'] = (adult_resampled.year > 9).astype(int)
return adult_resampled
denom_subset.groupby(['over64', 'year']).n.sum().reset_index()
over64 | year | n | |
---|---|---|---|
0 | 0 | 0 | 1822327.6 |
1 | 0 | 1 | 1843273.0 |
2 | 0 | 2 | 1864774.0 |
3 | 0 | 3 | 1891701.6 |
4 | 0 | 4 | 1921488.6 |
5 | 0 | 5 | 1958459.8 |
6 | 0 | 6 | 1991248.6 |
7 | 0 | 7 | 2014026.8 |
8 | 0 | 8 | 2038777.6 |
9 | 0 | 9 | 2062440.8 |
10 | 0 | 10 | 2086422.2 |
11 | 0 | 11 | 2107241.4 |
12 | 0 | 12 | 2121346.2 |
13 | 1 | 0 | 306863.0 |
14 | 1 | 1 | 307936.0 |
15 | 1 | 2 | 311000.0 |
16 | 1 | 3 | 314250.0 |
17 | 1 | 4 | 320068.0 |
18 | 1 | 5 | 327686.0 |
19 | 1 | 6 | 334378.0 |
20 | 1 | 7 | 344093.0 |
21 | 1 | 8 | 352267.0 |
22 | 1 | 9 | 360882.0 |
23 | 1 | 10 | 369262.0 |
24 | 1 | 11 | 387854.0 |
25 | 1 | 12 | 404681.0 |
def rate_model(variable, regional=True):
groupby_ind = [pd.Grouper(freq='AS', key='culture'), 'over64']
if regional:
groupby_ind += ['regioneast']
data_grouped = data_subset.groupby(groupby_ind)[variable]
adult_resampled = pd.DataFrame(dict(count = data_grouped.sum())).reset_index().rename(columns={'culture':'year'})
adult_resampled['year'] = [d.year-2001 for d in adult_resampled.year]
adult_resampled['post_2010'] = (adult_resampled.year > 9).astype(int)
if regional:
data_merged = pd.merge(denom_subset, adult_resampled, on=['year', 'regioneast', 'over64'])
data_merged = data_merged[['year', 'n', 'count', 'regioneast', 'over64']]
else:
denom_subset_state = denom_subset.groupby(['over64', 'year']).n.sum().reset_index()
data_merged = pd.merge(denom_subset_state, adult_resampled, on=['year', 'over64'])
data_merged = data_merged[['year', 'n', 'count', 'over64']]
data_merged['post_2010'] = (data_merged.year > 9).astype(int)
count = data_merged['count'].values
year = data_merged['year'].values
# Make rates per 1000
n = (data_merged['n']/1000).astype(int).values
if regional:
X = data_merged[['regioneast', 'over64']].values
regioneast = X[:,0].astype(bool)
regionwest = ~regioneast
k = X.shape[1]
over64 = X[:,1].astype(bool)
under65 = ~over64
else:
X = data_merged['over64'].values
k = 1
over64 = X.astype(bool)
under65 = ~over64
nyears = len(set(year))
mu_alpha = pm.Normal('mu_alpha', 0, 0.01, value=0)
sigma_alpha = pm.Uniform('sigma_alpha', 0, 100, value=10)
rho = pm.Uniform('rho', -1, 1, value=0)
# mu = pm.Lambda('mu', lambda mu=mu_alpha, d=pcv13: [mu]*nyears + d*post_2010[:12])
off_diag = np.eye(nyears, k=1)
Sigma = pm.Lambda('Sigma', lambda s=sigma_alpha, r=rho: (np.eye(nyears) + off_diag*r + off_diag.T*r)*(s**2))
alpha = pm.MvNormalCov('alpha', [mu_alpha]*nyears, Sigma, value=[0]*nyears)
mu_alpha = pm.Normal('mu_alpha', 0, 0.01, value=0)
sigma_alpha = pm.Uniform('sigma_alpha', 0, 100, value=10)
alpha = pm.Normal('alpha', mu_alpha, sigma_alpha**-2, value=[0]*nyears)
# Covariate effects
beta = pm.Normal('beta', 0, 0.001, value=[0]*k)
# Noise term
tau = pm.Exponential('tau', 1, value=1)
epsilon = pm.Normal('epsilon', 0, tau, value=np.zeros(len(n)))
# Poisson regression expected value
if regional:
rate = pm.Lambda('rate', lambda alpha=alpha, beta=beta, e=epsilon: pm.invlogit(alpha[year] + X.dot(beta) + e))
else:
rate = pm.Lambda('rate', lambda alpha=alpha, beta=beta, e=epsilon: pm.invlogit(alpha[year] + X*beta + e))
# Calculate region/age-specific rates (per 100K)
if regional:
rate_west_under65 = pm.Lambda('rate_west_under65', lambda r=rate: r[regionwest & under65]*100)
rate_west_over64 = pm.Lambda('rate_west_over64', lambda r=rate: r[regionwest & over64]*100)
rate_east_under65 = pm.Lambda('rate_east_under65', lambda r=rate: r[regioneast & under65]*100)
rate_east_over64 = pm.Lambda('rate_east_over64', lambda r=rate: r[regioneast & over64]*100)
else:
rate_under65 = pm.Lambda('rate_under65', lambda r=rate: r[under65]*100)
rate_over64 = pm.Lambda('rate_over64', lambda r=rate: r[over64]*100)
# Likelihood function
y = pm.Binomial('y', n, rate, value=count, observed=True)
# Posterior predictive distribution for all observations (for goodness-of-fit)
y_pred = pm.Binomial('y_pred', n, rate, value=count)
return(locals())
M_pneumo7 = pm.MCMC(rate_model('pneumo7'))
M_pneumo7.sample(50000, 40000)
[-----------------100%-----------------] 50000 of 50000 complete in 28.9 sec
M_pneumo7.sample(50000, 40000)
[-----------------100%-----------------] 50000 of 50000 complete in 27.2 sec
pm.Matplot.summary_plot(M_pneumo7.alpha)
M_pneumo7.rate_east_under65.summary()
rate_east_under65: Mean SD MC Error 95% HPD interval ------------------------------------------------------------------ 6.563 1.154 0.072 [ 4.385 8.867] 4.349 0.93 0.059 [ 2.529 6.148] 4.456 0.951 0.063 [ 2.744 6.447] 1.657 0.501 0.035 [ 0.758 2.621] 1.984 0.513 0.035 [ 1.102 3.111] 1.141 0.373 0.028 [ 0.539 1.897] 0.992 0.376 0.03 [ 0.349 1.735] 0.792 0.289 0.023 [ 0.334 1.381] 0.651 0.235 0.018 [ 0.237 1.126] 0.331 0.171 0.014 [ 0.071 0.678] 0.397 0.178 0.014 [ 0.101 0.755] 0.383 0.195 0.017 [ 0.099 0.779] 0.399 0.183 0.015 [ 0.11 0.762] Posterior quantiles: 2.5 25 50 75 97.5 |---------------|===============|===============|---------------| 4.503 5.769 6.507 7.275 9.024 2.696 3.687 4.28 4.944 6.341 2.756 3.776 4.395 5.075 6.473 0.841 1.295 1.617 1.962 2.795 1.094 1.621 1.941 2.308 3.104 0.574 0.864 1.083 1.374 1.976 0.461 0.731 0.923 1.18 1.975 0.374 0.577 0.748 0.961 1.483 0.284 0.483 0.623 0.775 1.234 0.088 0.211 0.298 0.415 0.767 0.131 0.265 0.369 0.498 0.815 0.125 0.237 0.344 0.481 0.855 0.142 0.266 0.366 0.493 0.881
M_pneumo7.rate_west_under65.summary()
rate_west_under65: Mean SD MC Error 95% HPD interval ------------------------------------------------------------------ 5.825 0.611 0.027 [ 4.602 7.018] 4.094 0.517 0.025 [ 3.1 5.095] 2.971 0.437 0.024 [ 2.144 3.829] 1.9 0.324 0.018 [ 1.324 2.578] 2.305 0.368 0.019 [ 1.635 3.026] 1.599 0.304 0.018 [ 1.014 2.226] 0.853 0.214 0.015 [ 0.458 1.292] 0.596 0.169 0.012 [ 0.302 0.936] 0.547 0.171 0.013 [ 0.246 0.873] 0.334 0.122 0.009 [ 0.132 0.576] 0.36 0.127 0.01 [ 0.16 0.609] 0.179 0.077 0.006 [ 0.046 0.328] 0.243 0.096 0.008 [ 0.091 0.429] Posterior quantiles: 2.5 25 50 75 97.5 |---------------|===============|===============|---------------| 4.648 5.418 5.804 6.214 7.091 3.165 3.728 4.074 4.415 5.17 2.181 2.663 2.943 3.261 3.879 1.332 1.667 1.876 2.107 2.595 1.666 2.037 2.281 2.545 3.069 1.028 1.388 1.585 1.791 2.249 0.488 0.7 0.832 0.983 1.332 0.324 0.47 0.576 0.7 0.981 0.273 0.423 0.529 0.645 0.93 0.161 0.245 0.313 0.399 0.632 0.182 0.265 0.339 0.433 0.662 0.064 0.12 0.169 0.225 0.36 0.104 0.171 0.228 0.3 0.458
M_pneumo7.rate_east_over64.summary()
rate_east_over64: Mean SD MC Error 95% HPD interval ------------------------------------------------------------------ 20.474 3.786 0.264 [ 13.695 28.37 ] 8.339 2.423 0.189 [ 3.952 13.102] 13.345 2.97 0.205 [ 7.595 19.214] 9.866 2.832 0.225 [ 4.749 15.748] 5.846 2.066 0.165 [ 2.123 9.93 ] 4.076 1.531 0.127 [ 1.456 7.031] 2.907 1.363 0.116 [ 0.869 5.497] 1.652 0.689 0.056 [ 0.531 2.987] 1.586 0.675 0.054 [ 0.467 2.874] 1.152 0.489 0.038 [ 0.317 2.141] 1.109 0.549 0.045 [ 0.269 2.138] 0.899 0.691 0.062 [ 0.128 2.263] 0.92 0.481 0.041 [ 0.138 1.822] Posterior quantiles: 2.5 25 50 75 97.5 |---------------|===============|===============|---------------| 13.876 17.781 20.181 22.804 28.641 4.232 6.571 8.067 9.876 13.522 8.071 11.212 13.189 15.238 19.809 5.157 7.879 9.498 11.535 16.275 2.611 4.305 5.515 7.168 10.599 1.614 3.054 3.923 4.878 7.492 1.133 2.007 2.672 3.483 6.493 0.674 1.127 1.539 2.049 3.26 0.548 1.088 1.502 1.982 3.137 0.411 0.803 1.087 1.4 2.334 0.387 0.736 1.011 1.333 2.574 0.234 0.489 0.698 1.052 2.8 0.2 0.585 0.872 1.179 2.07
M_pneumo7.rate_west_over64.summary()
rate_west_over64: Mean SD MC Error 95% HPD interval ------------------------------------------------------------------ 19.778 2.447 0.123 [ 15.132 24.547] 13.049 2.175 0.128 [ 8.75 17.372] 11.408 2.014 0.123 [ 7.434 15.313] 5.277 1.343 0.092 [ 2.897 7.964] 4.082 1.098 0.074 [ 2.171 6.4 ] 2.688 0.849 0.062 [ 1.224 4.431] 2.375 0.758 0.057 [ 1.063 3.941] 1.584 0.55 0.041 [ 0.649 2.67 ] 1.443 0.603 0.048 [ 0.386 2.566] 1.098 0.511 0.043 [ 0.289 2.086] 0.958 0.446 0.037 [ 0.282 1.829] 0.539 0.242 0.019 [ 0.165 1.039] 0.658 0.307 0.025 [ 0.175 1.234] Posterior quantiles: 2.5 25 50 75 97.5 |---------------|===============|===============|---------------| 15.236 18.04 19.715 21.466 24.676 8.91 11.568 12.955 14.426 17.565 7.749 10.003 11.342 12.668 15.722 3.096 4.295 5.132 6.077 8.285 2.199 3.3 4.011 4.771 6.473 1.305 2.076 2.568 3.209 4.573 1.122 1.829 2.294 2.853 4.063 0.746 1.173 1.509 1.914 2.841 0.508 1.008 1.392 1.764 2.842 0.395 0.742 0.988 1.367 2.323 0.361 0.623 0.872 1.194 2.006 0.201 0.362 0.495 0.664 1.151 0.221 0.418 0.609 0.844 1.37
def make_rates_df(mod):
rate_west_over64 = pd.DataFrame(mod.rate_west_over64.trace())
rate_west_over64['Region'] = 'West'
rate_west_over64['Age'] = 'Over 64'
rate_east_over64 = pd.DataFrame(mod.rate_east_over64.trace())
rate_east_over64['Region'] = 'East'
rate_east_over64['Age'] = 'Over 64'
rate_west_under65 = pd.DataFrame(mod.rate_west_under65.trace())
rate_west_under65['Region'] = 'West'
rate_west_under65['Age'] = 'Under 65'
rate_east_under65 = pd.DataFrame(mod.rate_east_under65.trace())
rate_east_under65['Region'] = 'East'
rate_east_under65['Age'] = 'Under 65'
rates = pd.concat([rate_west_over64, rate_east_over64, rate_west_under65, rate_east_under65],
ignore_index=True)
names = rates.columns.values.copy()
names[:13] += 2001
rates.columns = names
return rates
import seaborn as sb
def make_factorplot(data_frame):
sb.factorplot(row="Age", col="Region", data=data_frame, kind='box',
row_order=['Under 65', 'Over 64'],
col_order=['West', 'East'],
palette="YlGnBu_d", linewidth=0.7, fliersize=0, aspect=1.25).despine(left=True)
pneumo7_rates = make_rates_df(M_pneumo7)
make_factorplot(pneumo7_rates)
# sb.set_context("talk", font_scale=0.8)
sb.set_style("whitegrid")
f, axes = plt.subplots(figsize=(16,4), sharey=True)
sb.boxplot(data=M_pneumo7.rate_east_over64.trace(), linewidth=0.2, fliersize=0, color='Red')
sb.boxplot(data=M_pneumo7.rate_west_over64.trace(), linewidth=0.2, fliersize=0, color='Blue')
# plt.xlabel('Age group'); plt.ylabel('Susceptibles');
import seaborn as sb
# sb.set_context("talk", font_scale=0.8)
sb.set_style("whitegrid")
f, axes = plt.subplots(figsize=(16,4), sharey=True)
sb.boxplot(data=M_pneumo7.rate_east_under65.trace(), linewidth=0.2, fliersize=0, color='Red')
sb.boxplot(data=M_pneumo7.rate_west_under65.trace(), linewidth=0.2, fliersize=0, color='Blue')
<matplotlib.axes._subplots.AxesSubplot at 0x1232ebb00>
pneumo7_east_over64 = pd.DataFrame(M_pneumo7.rate_east_over64.trace())
pneumo7_east_over64['Region'] = 'East'
pneumo7_east_over64['Age'] = 'Over 64'
pneumo7_west_over64 = pd.DataFrame(M_pneumo7.rate_west_over64.trace())
pneumo7_west_over64['Region'] = 'West'
pneumo7_west_over64['Age'] = 'Over 64'
pneumo7_east_under65 = pd.DataFrame(M_pneumo7.rate_east_under65.trace())
pneumo7_east_under65['Region'] = 'East'
pneumo7_east_under65['Age'] = 'Under 65'
pneumo7_west_under65 = pd.DataFrame(M_pneumo7.rate_west_under65.trace())
pneumo7_west_under65['Region'] = 'West'
pneumo7_west_under65['Age'] = 'Under 65'
pneumo7 = pd.concat([pneumo7_east_over64, pneumo7_east_under65,
pneumo7_west_over64, pneumo7_west_under65], ignore_index=True)
pm.Matplot.summary_plot(M_pneumo7.rate_east_over64)
pm.Matplot.summary_plot(M_pneumo7.rate_west_over64)
pm.Matplot.summary_plot(M_pneumo7.rate_east_under65)
pm.Matplot.summary_plot(M_pneumo7.rate_west_under65)
M_pneumo7_tenn = pm.MCMC(rate_model('pneumo7', regional=False))
M_pneumo7_tenn.sample(50000, 40000)
[-----------------100%-----------------] 50000 of 50000 complete in 28.8 sec
pm.Matplot.plot(M_pneumo7_tenn.beta)
Plotting beta_0
pm.Matplot.summary_plot(M_pneumo7_tenn.rate_over64)
Could not calculate Gelman-Rubin statistics. Requires multiple chains of equal length.
pm.Matplot.summary_plot(M_pneumo7_tenn.rate_under65)
Could not calculate Gelman-Rubin statistics. Requires multiple chains of equal length.
M_pneumo7_tenn.rate_under65.summary()
rate_under65: Mean SD MC Error 95% HPD interval ------------------------------------------------------------------ 5.933 0.541 0.026 [ 4.949 7.111] 4.227 0.473 0.026 [ 3.273 5.148] 3.385 0.418 0.023 [ 2.598 4.177] 1.746 0.312 0.02 [ 1.202 2.42 ] 2.263 0.346 0.021 [ 1.601 2.946] 1.504 0.277 0.018 [ 0.981 2.049] 0.914 0.2 0.014 [ 0.545 1.297] 0.64 0.153 0.012 [ 0.36 0.968] 0.59 0.151 0.012 [ 0.335 0.908] 0.321 0.111 0.01 [ 0.114 0.545] 0.361 0.118 0.01 [ 0.138 0.594] 0.238 0.109 0.01 [ 0.061 0.429] 0.27 0.104 0.009 [ 0.112 0.516] Posterior quantiles: 2.5 25 50 75 97.5 |---------------|===============|===============|---------------| 4.912 5.571 5.908 6.287 7.079 3.334 3.9 4.207 4.537 5.23 2.647 3.101 3.362 3.637 4.276 1.184 1.524 1.72 1.942 2.409 1.619 2.009 2.247 2.49 2.983 1.046 1.301 1.469 1.68 2.139 0.578 0.768 0.895 1.043 1.349 0.378 0.53 0.626 0.731 1.0 0.346 0.477 0.573 0.685 0.928 0.128 0.244 0.307 0.386 0.591 0.151 0.277 0.356 0.435 0.617 0.084 0.16 0.227 0.285 0.51 0.119 0.199 0.25 0.323 0.534
M_pneumo7_tenn.rate_over64.summary()
rate_over64: Mean SD MC Error 95% HPD interval ------------------------------------------------------------------ 20.374 2.407 0.138 [ 16.088 25.37 ] 11.172 1.832 0.129 [ 7.553 14.709] 11.758 1.802 0.109 [ 8.446 15.281] 7.285 1.3 0.085 [ 4.806 9.767] 4.23 1.204 0.098 [ 2.138 6.54 ] 3.11 0.89 0.072 [ 1.565 4.826] 2.46 0.837 0.07 [ 0.972 4.204] 1.597 0.505 0.043 [ 0.735 2.592] 1.354 0.528 0.046 [ 0.439 2.304] 1.19 0.47 0.041 [ 0.434 2.116] 0.977 0.484 0.044 [ 0.127 1.794] 0.725 0.346 0.031 [ 0.242 1.444] 0.806 0.424 0.039 [ 0.161 1.651] Posterior quantiles: 2.5 25 50 75 97.5 |---------------|===============|===============|---------------| 15.935 18.667 20.281 22.02 25.264 7.783 9.876 11.125 12.327 15.064 8.601 10.415 11.642 12.962 15.623 4.924 6.351 7.235 8.134 9.995 2.323 3.362 4.097 4.901 7.057 1.654 2.465 3.046 3.649 5.054 1.101 1.876 2.346 2.919 4.449 0.823 1.208 1.533 1.917 2.708 0.537 0.955 1.289 1.691 2.496 0.504 0.83 1.103 1.485 2.285 0.187 0.615 0.958 1.3 1.972 0.307 0.491 0.641 0.872 1.668 0.21 0.469 0.737 1.049 1.813
M_pneumo13 = pm.MCMC(rate_model('pneumo13'))
M_pneumo13.sample(50000, 40000)
# M_pneumo13.sample(50000, 40000)
[-----------------100%-----------------] 50000 of 50000 complete in 29.2 sec
pneumo13_rates = make_rates_df(M_pneumo13)
make_factorplot(pneumo13_rates)
pm.Matplot.summary_plot(M_pneumo13.rate_east_under65)
Could not calculate Gelman-Rubin statistics. Requires multiple chains of equal length.
pm.Matplot.summary_plot(M_pneumo13.rate_east_over64)
Could not calculate Gelman-Rubin statistics. Requires multiple chains of equal length.
pm.Matplot.summary_plot(M_pneumo13.rate_west_under65)
Could not calculate Gelman-Rubin statistics. Requires multiple chains of equal length.
pm.Matplot.summary_plot(M_pneumo13.rate_west_over64)
Could not calculate Gelman-Rubin statistics. Requires multiple chains of equal length.
M_pneumo13.rate_east_under65.summary()
rate_east_under65: Mean SD MC Error 95% HPD interval ------------------------------------------------------------------ 2.455 0.648 0.058 [ 1.313 3.738] 2.408 0.523 0.044 [ 1.327 3.427] 3.026 0.601 0.049 [ 1.877 4.244] 2.627 0.558 0.046 [ 1.556 3.666] 5.38 1.0 0.082 [ 3.513 7.393] 6.16 0.966 0.074 [ 4.287 8.047] 4.557 0.778 0.062 [ 3.106 6.094] 7.371 1.06 0.08 [ 5.283 9.34 ] 7.475 1.098 0.087 [ 5.245 9.642] 5.186 1.002 0.086 [ 3.365 7.298] 4.091 0.751 0.062 [ 2.668 5.488] 2.816 0.576 0.048 [ 1.784 3.943] 1.806 0.506 0.046 [ 0.9 2.836] Posterior quantiles: 2.5 25 50 75 97.5 |---------------|===============|===============|---------------| 1.459 1.999 2.352 2.797 4.109 1.479 2.057 2.356 2.696 3.658 1.983 2.601 2.97 3.405 4.406 1.572 2.218 2.627 3.019 3.695 3.556 4.691 5.315 6.053 7.475 4.43 5.476 6.102 6.762 8.284 3.167 3.997 4.516 5.061 6.188 5.533 6.602 7.298 8.078 9.631 5.395 6.71 7.418 8.189 9.846 3.408 4.513 5.1 5.802 7.359 2.783 3.54 4.036 4.568 5.695 1.824 2.385 2.777 3.227 3.995 0.987 1.457 1.751 2.071 3.025
M_pneumo7.rate_west_under65.summary()
rate_west_under65: Mean SD MC Error 95% HPD interval ------------------------------------------------------------------ 5.825 0.611 0.027 [ 4.602 7.018] 4.094 0.517 0.025 [ 3.1 5.095] 2.971 0.437 0.024 [ 2.144 3.829] 1.9 0.324 0.018 [ 1.324 2.578] 2.305 0.368 0.019 [ 1.635 3.026] 1.599 0.304 0.018 [ 1.014 2.226] 0.853 0.214 0.015 [ 0.458 1.292] 0.596 0.169 0.012 [ 0.302 0.936] 0.547 0.171 0.013 [ 0.246 0.873] 0.334 0.122 0.009 [ 0.132 0.576] 0.36 0.127 0.01 [ 0.16 0.609] 0.179 0.077 0.006 [ 0.046 0.328] 0.243 0.096 0.008 [ 0.091 0.429] Posterior quantiles: 2.5 25 50 75 97.5 |---------------|===============|===============|---------------| 4.648 5.418 5.804 6.214 7.091 3.165 3.728 4.074 4.415 5.17 2.181 2.663 2.943 3.261 3.879 1.332 1.667 1.876 2.107 2.595 1.666 2.037 2.281 2.545 3.069 1.028 1.388 1.585 1.791 2.249 0.488 0.7 0.832 0.983 1.332 0.324 0.47 0.576 0.7 0.981 0.273 0.423 0.529 0.645 0.93 0.161 0.245 0.313 0.399 0.632 0.182 0.265 0.339 0.433 0.662 0.064 0.12 0.169 0.225 0.36 0.104 0.171 0.228 0.3 0.458
M_pneumo7.rate_east_over64.summary()
rate_east_over64: Mean SD MC Error 95% HPD interval ------------------------------------------------------------------ 20.474 3.786 0.264 [ 13.695 28.37 ] 8.339 2.423 0.189 [ 3.952 13.102] 13.345 2.97 0.205 [ 7.595 19.214] 9.866 2.832 0.225 [ 4.749 15.748] 5.846 2.066 0.165 [ 2.123 9.93 ] 4.076 1.531 0.127 [ 1.456 7.031] 2.907 1.363 0.116 [ 0.869 5.497] 1.652 0.689 0.056 [ 0.531 2.987] 1.586 0.675 0.054 [ 0.467 2.874] 1.152 0.489 0.038 [ 0.317 2.141] 1.109 0.549 0.045 [ 0.269 2.138] 0.899 0.691 0.062 [ 0.128 2.263] 0.92 0.481 0.041 [ 0.138 1.822] Posterior quantiles: 2.5 25 50 75 97.5 |---------------|===============|===============|---------------| 13.876 17.781 20.181 22.804 28.641 4.232 6.571 8.067 9.876 13.522 8.071 11.212 13.189 15.238 19.809 5.157 7.879 9.498 11.535 16.275 2.611 4.305 5.515 7.168 10.599 1.614 3.054 3.923 4.878 7.492 1.133 2.007 2.672 3.483 6.493 0.674 1.127 1.539 2.049 3.26 0.548 1.088 1.502 1.982 3.137 0.411 0.803 1.087 1.4 2.334 0.387 0.736 1.011 1.333 2.574 0.234 0.489 0.698 1.052 2.8 0.2 0.585 0.872 1.179 2.07
M_pneumo7.rate_west_over64.summary()
rate_west_over64: Mean SD MC Error 95% HPD interval ------------------------------------------------------------------ 19.778 2.447 0.123 [ 15.132 24.547] 13.049 2.175 0.128 [ 8.75 17.372] 11.408 2.014 0.123 [ 7.434 15.313] 5.277 1.343 0.092 [ 2.897 7.964] 4.082 1.098 0.074 [ 2.171 6.4 ] 2.688 0.849 0.062 [ 1.224 4.431] 2.375 0.758 0.057 [ 1.063 3.941] 1.584 0.55 0.041 [ 0.649 2.67 ] 1.443 0.603 0.048 [ 0.386 2.566] 1.098 0.511 0.043 [ 0.289 2.086] 0.958 0.446 0.037 [ 0.282 1.829] 0.539 0.242 0.019 [ 0.165 1.039] 0.658 0.307 0.025 [ 0.175 1.234] Posterior quantiles: 2.5 25 50 75 97.5 |---------------|===============|===============|---------------| 15.236 18.04 19.715 21.466 24.676 8.91 11.568 12.955 14.426 17.565 7.749 10.003 11.342 12.668 15.722 3.096 4.295 5.132 6.077 8.285 2.199 3.3 4.011 4.771 6.473 1.305 2.076 2.568 3.209 4.573 1.122 1.829 2.294 2.853 4.063 0.746 1.173 1.509 1.914 2.841 0.508 1.008 1.392 1.764 2.842 0.395 0.742 0.988 1.367 2.323 0.361 0.623 0.872 1.194 2.006 0.201 0.362 0.495 0.664 1.151 0.221 0.418 0.609 0.844 1.37
M_pneumo13_tenn = pm.MCMC(rate_model('pneumo13', regional=False))
M_pneumo13_tenn.sample(50000, 40000)
[-----------------100%-----------------] 50000 of 50000 complete in 28.5 sec
pm.Matplot.summary_plot(M_pneumo13_tenn.alpha)
Could not calculate Gelman-Rubin statistics. Requires multiple chains of equal length.
pm.Matplot.plot(M_pneumo13_tenn.beta)
Plotting beta_0
pm.Matplot.summary_plot(M_pneumo13_tenn.rate_over64)
Could not calculate Gelman-Rubin statistics. Requires multiple chains of equal length.
pm.Matplot.summary_plot(M_pneumo13_tenn.rate_under65)
Could not calculate Gelman-Rubin statistics. Requires multiple chains of equal length.
M_pneumo13_tenn.rate_under65.summary()
rate_under65: Mean SD MC Error 95% HPD interval ------------------------------------------------------------------ 2.386 0.333 0.027 [ 1.8 3.073] 2.728 0.368 0.029 [ 2.026 3.427] 2.892 0.36 0.026 [ 2.255 3.663] 2.552 0.333 0.026 [ 1.919 3.212] 4.029 0.433 0.031 [ 3.209 4.895] 5.545 0.482 0.031 [ 4.664 6.547] 4.859 0.474 0.033 [ 3.988 5.846] 6.011 0.55 0.038 [ 4.994 7.134] 5.79 0.539 0.038 [ 4.705 6.797] 4.433 0.41 0.026 [ 3.604 5.186] 3.855 0.39 0.027 [ 3.111 4.636] 2.287 0.302 0.024 [ 1.724 2.948] 1.674 0.272 0.022 [ 1.104 2.198] Posterior quantiles: 2.5 25 50 75 97.5 |---------------|===============|===============|---------------| 1.83 2.149 2.357 2.588 3.119 2.051 2.465 2.72 2.976 3.459 2.242 2.638 2.87 3.125 3.656 1.991 2.324 2.528 2.746 3.307 3.172 3.735 4.013 4.339 4.888 4.638 5.219 5.534 5.848 6.531 3.988 4.526 4.843 5.186 5.846 4.98 5.648 6.0 6.371 7.131 4.744 5.413 5.791 6.141 6.853 3.671 4.167 4.425 4.692 5.281 3.12 3.597 3.838 4.106 4.654 1.691 2.097 2.274 2.471 2.929 1.153 1.475 1.665 1.841 2.251
M_pneumo13_tenn.rate_over64.summary()
rate_over64: Mean SD MC Error 95% HPD interval ------------------------------------------------------------------ 10.905 1.559 0.123 [ 7.996 13.661] 10.226 1.647 0.139 [ 7.18 13.481] 11.451 1.645 0.127 [ 8.516 14.86 ] 10.294 1.35 0.099 [ 7.885 12.909] 12.215 1.645 0.125 [ 9.078 15.588] 12.409 1.743 0.14 [ 9.207 15.959] 15.947 2.028 0.158 [ 12.244 20.392] 13.887 1.784 0.142 [ 10.732 17.451] 9.957 1.42 0.114 [ 7.619 12.878] 12.889 1.615 0.123 [ 9.922 16.332] 9.989 1.471 0.12 [ 7.098 12.874] 6.488 1.177 0.101 [ 4.204 8.86 ] 4.374 0.87 0.074 [ 2.827 6.114] Posterior quantiles: 2.5 25 50 75 97.5 |---------------|===============|===============|---------------| 8.101 9.732 10.918 12.017 13.847 7.408 9.049 10.145 11.247 13.772 8.606 10.252 11.287 12.541 15.002 7.943 9.309 10.209 11.209 13.044 9.336 11.053 12.105 13.228 15.884 9.307 11.166 12.326 13.536 16.113 12.101 14.545 15.764 17.206 20.316 10.859 12.65 13.696 14.953 17.826 7.611 8.902 9.888 10.892 12.877 9.801 11.81 12.819 13.916 16.23 7.269 8.944 9.961 10.959 13.154 4.36 5.696 6.41 7.145 9.066 2.883 3.796 4.267 4.867 6.26
M_pneumo23 = pm.MCMC(rate_model('pneumo23'))
M_pneumo23.sample(50000, 40000)
# M_pneumo23.sample(50000, 40000)
[-----------------100%-----------------] 50000 of 50000 complete in 30.5 sec
pneumo23_rates = make_rates_df(M_pneumo23)
make_factorplot(pneumo23_rates)
pm.Matplot.summary_plot(M_pneumo23.alpha)
Could not calculate Gelman-Rubin statistics. Requires multiple chains of equal length.
pm.Matplot.summary_plot(M_pneumo23.rate_west_over64)
Could not calculate Gelman-Rubin statistics. Requires multiple chains of equal length.
pm.Matplot.summary_plot(M_pneumo23.rate_west_under65)
Could not calculate Gelman-Rubin statistics. Requires multiple chains of equal length.
pm.Matplot.summary_plot(M_pneumo23.rate_east_over64)
Could not calculate Gelman-Rubin statistics. Requires multiple chains of equal length.
pm.Matplot.summary_plot(M_pneumo23.rate_east_under65)
Could not calculate Gelman-Rubin statistics. Requires multiple chains of equal length.
M_pneumo23.rate_east_under65.summary()
rate_east_under65: Mean SD MC Error 95% HPD interval ------------------------------------------------------------------ 2.026 0.491 0.043 [ 1.06 2.993] 2.159 0.456 0.038 [ 1.273 2.969] 2.338 0.539 0.048 [ 1.435 3.475] 3.422 0.807 0.07 [ 2.021 5.046] 5.625 0.944 0.078 [ 3.941 7.502] 4.266 0.746 0.06 [ 2.663 5.665] 6.613 0.99 0.075 [ 4.692 8.519] 4.773 0.859 0.071 [ 3.147 6.31 ] 4.887 0.755 0.056 [ 3.405 6.309] 5.281 0.958 0.079 [ 3.542 7.108] 5.037 0.803 0.061 [ 3.628 6.674] 5.32 0.805 0.059 [ 3.822 6.908] 6.416 0.998 0.079 [ 4.58 8.43] Posterior quantiles: 2.5 25 50 75 97.5 |---------------|===============|===============|---------------| 1.208 1.713 1.952 2.273 3.257 1.307 1.842 2.157 2.47 3.019 1.502 1.936 2.24 2.662 3.556 2.117 2.836 3.332 3.892 5.281 4.002 4.937 5.564 6.218 7.646 2.843 3.754 4.246 4.73 5.868 4.834 5.907 6.54 7.255 8.707 3.236 4.183 4.726 5.358 6.503 3.488 4.38 4.88 5.353 6.451 3.651 4.558 5.226 5.911 7.295 3.638 4.441 4.959 5.576 6.691 3.84 4.722 5.294 5.873 6.948 4.657 5.735 6.303 7.054 8.555
M_pneumo23.rate_west_under65.summary()
rate_west_under65: Mean SD MC Error 95% HPD interval ------------------------------------------------------------------ 2.663 0.358 0.027 [ 2.002 3.377] 2.672 0.382 0.029 [ 1.964 3.442] 3.074 0.435 0.033 [ 2.262 3.983] 3.252 0.423 0.03 [ 2.378 4.012] 5.824 0.601 0.041 [ 4.707 7.091] 4.251 0.47 0.029 [ 3.349 5.159] 3.228 0.423 0.031 [ 2.367 4.052] 3.254 0.441 0.034 [ 2.429 4.122] 3.154 0.391 0.029 [ 2.479 3.972] 3.727 0.473 0.034 [ 2.834 4.704] 3.91 0.485 0.036 [ 2.971 4.849] 4.055 0.48 0.033 [ 3.166 5.002] 3.623 0.409 0.029 [ 2.867 4.394] Posterior quantiles: 2.5 25 50 75 97.5 |---------------|===============|===============|---------------| 2.04 2.405 2.635 2.894 3.44 1.954 2.409 2.665 2.926 3.435 2.221 2.789 3.066 3.349 3.954 2.414 2.967 3.243 3.544 4.061 4.691 5.427 5.787 6.205 7.082 3.377 3.931 4.234 4.55 5.252 2.462 2.934 3.201 3.494 4.158 2.437 2.943 3.243 3.556 4.134 2.46 2.864 3.138 3.419 3.96 2.834 3.416 3.693 4.03 4.704 3.007 3.57 3.917 4.232 4.903 3.214 3.714 4.027 4.358 5.071 2.89 3.321 3.619 3.901 4.423
M_pneumo23.rate_east_over64.summary()
rate_east_over64: Mean SD MC Error 95% HPD interval ------------------------------------------------------------------ 7.238 1.597 0.137 [ 4.181 10.077] 5.606 1.636 0.15 [ 2.762 8.225] 7.961 1.828 0.153 [ 4.595 11.4 ] 7.811 1.772 0.155 [ 4.78 11.639] 9.571 1.715 0.142 [ 6.433 13.075] 11.645 2.431 0.211 [ 7.372 16.558] 13.184 2.636 0.231 [ 8.533 18.527] 12.789 2.429 0.201 [ 8.633 18.003] 11.389 1.988 0.164 [ 7.817 15.089] 11.332 1.776 0.139 [ 7.854 14.791] 13.345 2.61 0.227 [ 8.266 18.13 ] 12.4 2.373 0.201 [ 7.944 16.85 ] 11.735 2.257 0.188 [ 7.323 16.31 ] Posterior quantiles: 2.5 25 50 75 97.5 |---------------|===============|===============|---------------| 4.31 6.032 7.261 8.313 10.379 3.036 4.447 5.524 6.494 9.35 5.106 6.672 7.738 9.005 12.126 5.023 6.553 7.534 8.825 12.095 6.634 8.368 9.449 10.623 13.408 7.769 9.941 11.291 13.008 17.202 8.572 11.25 13.101 14.886 18.635 8.846 11.093 12.592 14.167 18.62 7.959 9.878 11.267 12.846 15.302 8.06 10.096 11.28 12.477 15.138 8.457 11.505 13.157 15.252 18.391 8.239 10.607 12.35 14.034 17.218 7.627 10.123 11.62 13.142 16.672
M_pneumo23.rate_west_over64.summary()
rate_west_over64: Mean SD MC Error 95% HPD interval ------------------------------------------------------------------ 7.122 1.3 0.107 [ 4.683 9.606] 4.772 1.178 0.104 [ 2.763 7.138] 7.327 1.364 0.111 [ 4.804 9.936] 5.852 1.099 0.091 [ 3.854 8.257] 10.114 1.548 0.122 [ 7.206 13.169] 11.651 1.899 0.157 [ 8.267 15.487] 10.352 1.824 0.151 [ 6.848 13.672] 10.543 1.862 0.155 [ 7.231 14.141] 7.974 1.627 0.143 [ 4.965 11.135] 10.822 1.858 0.152 [ 7.387 14.432] 11.853 1.784 0.145 [ 8.588 15.533] 10.8 1.703 0.137 [ 7.228 13.881] 9.893 1.548 0.119 [ 7.054 13.14 ] Posterior quantiles: 2.5 25 50 75 97.5 |---------------|===============|===============|---------------| 4.852 6.173 6.997 7.971 9.996 2.907 3.927 4.624 5.468 7.472 4.861 6.351 7.271 8.294 10.025 3.931 5.113 5.779 6.491 8.374 7.397 8.986 10.015 11.126 13.384 8.344 10.287 11.457 12.965 15.636 7.186 9.067 10.262 11.484 14.13 7.343 9.164 10.421 11.829 14.349 5.325 6.864 7.723 8.857 12.148 7.465 9.484 10.709 12.091 14.542 8.664 10.585 11.754 12.961 15.626 7.57 9.58 10.784 11.88 14.321 7.054 8.802 9.815 10.864 13.14
M_pneumo23_tenn = pm.MCMC(rate_model('pneumo23', regional=False))
M_pneumo23_tenn.sample(50000, 40000)
[-----------------100%-----------------] 50000 of 50000 complete in 25.9 sec
pm.Matplot.summary_plot(M_pneumo23_tenn.rate_over64)
Could not calculate Gelman-Rubin statistics. Requires multiple chains of equal length.
pm.Matplot.summary_plot(M_pneumo23_tenn.rate_under65)
Could not calculate Gelman-Rubin statistics. Requires multiple chains of equal length.
M_pneumo23_tenn.rate_under65.summary()
rate_under65: Mean SD MC Error 95% HPD interval ------------------------------------------------------------------ 2.438 0.335 0.026 [ 1.798 3.085] 2.57 0.347 0.028 [ 1.94 3.245] 2.751 0.326 0.024 [ 2.127 3.398] 3.345 0.407 0.031 [ 2.592 4.166] 5.884 0.502 0.031 [ 4.931 6.88 ] 4.218 0.445 0.032 [ 3.344 5.059] 3.951 0.402 0.027 [ 3.227 4.809] 3.527 0.414 0.029 [ 2.74 4.335] 3.578 0.359 0.024 [ 2.897 4.292] 4.042 0.404 0.028 [ 3.338 4.869] 4.112 0.42 0.031 [ 3.386 4.991] 4.409 0.423 0.027 [ 3.609 5.287] 4.336 0.442 0.031 [ 3.446 5.228] Posterior quantiles: 2.5 25 50 75 97.5 |---------------|===============|===============|---------------| 1.836 2.212 2.426 2.651 3.125 1.947 2.321 2.553 2.822 3.254 2.135 2.515 2.751 2.97 3.411 2.627 3.063 3.315 3.601 4.231 4.95 5.54 5.869 6.217 6.91 3.375 3.92 4.206 4.508 5.105 3.206 3.675 3.92 4.217 4.799 2.75 3.246 3.51 3.8 4.359 2.9 3.328 3.58 3.824 4.296 3.336 3.76 4.022 4.303 4.869 3.389 3.824 4.075 4.364 5.024 3.59 4.141 4.393 4.674 5.273 3.512 4.05 4.301 4.602 5.307
M_pneumo23_tenn.rate_over64.summary()
rate_over64: Mean SD MC Error 95% HPD interval ------------------------------------------------------------------ 7.056 1.427 0.124 [ 4.597 10.023] 5.227 1.041 0.092 [ 3.094 7.17 ] 8.58 1.431 0.117 [ 6.065 11.581] 6.589 1.095 0.089 [ 4.625 8.775] 9.642 1.409 0.111 [ 6.945 12.436] 10.671 1.675 0.135 [ 7.255 14.054] 11.569 1.571 0.12 [ 8.517 14.638] 11.389 1.738 0.144 [ 8.28 14.732] 8.412 1.26 0.099 [ 5.933 10.807] 10.781 1.499 0.119 [ 7.784 13.639] 12.245 1.633 0.127 [ 8.795 15.372] 11.172 1.495 0.11 [ 8.42 14.23] 10.469 1.412 0.11 [ 7.664 13.235] Posterior quantiles: 2.5 25 50 75 97.5 |---------------|===============|===============|---------------| 4.764 6.019 6.852 7.942 10.227 3.231 4.477 5.24 5.926 7.348 5.959 7.504 8.546 9.552 11.478 4.767 5.794 6.49 7.271 9.102 7.065 8.663 9.63 10.537 12.682 7.474 9.604 10.568 11.705 14.394 8.533 10.519 11.555 12.62 14.665 8.391 10.111 11.306 12.521 15.013 6.03 7.539 8.396 9.24 10.936 8.01 9.718 10.743 11.792 13.898 9.088 11.118 12.165 13.307 15.742 8.457 10.161 11.075 12.127 14.29 7.801 9.484 10.473 11.358 13.401
M_pneumononPCV = pm.MCMC(rate_model('pneumononPCV'))
M_pneumononPCV.sample(50000, 40000)
# M_pneumononPCV.sample(50000, 40000)
[-----------------100%-----------------] 50000 of 50000 complete in 31.7 sec
pneumononpcv_rates = make_rates_df(M_pneumononPCV)
make_factorplot(pneumononpcv_rates)
pm.Matplot.summary_plot(M_pneumononPCV.rate_west_over64)
Could not calculate Gelman-Rubin statistics. Requires multiple chains of equal length.
pm.Matplot.summary_plot(M_pneumononPCV.rate_west_under65)
Could not calculate Gelman-Rubin statistics. Requires multiple chains of equal length.
pm.Matplot.summary_plot(M_pneumononPCV.rate_east_over64)
Could not calculate Gelman-Rubin statistics. Requires multiple chains of equal length.
pm.Matplot.summary_plot(M_pneumononPCV.rate_east_under65)
Could not calculate Gelman-Rubin statistics. Requires multiple chains of equal length.
M_pneumononPCV.rate_east_under65.summary()
rate_east_under65: Mean SD MC Error 95% HPD interval ------------------------------------------------------------------ 3.136 0.636 0.056 [ 1.938 4.371] 2.9 0.564 0.05 [ 1.809 3.989] 3.538 0.629 0.055 [ 2.355 4.724] 4.64 1.021 0.094 [ 2.606 6.408] 6.535 1.07 0.092 [ 4.37 8.541] 6.265 0.772 0.062 [ 4.788 7.826] 8.575 1.184 0.099 [ 6.044 10.671] 8.081 1.173 0.1 [ 5.942 10.597] 7.634 1.143 0.1 [ 5.519 9.872] 7.896 1.042 0.085 [ 5.922 9.916] 8.195 1.242 0.106 [ 5.925 10.565] 7.88 1.108 0.092 [ 5.658 9.971] 8.983 1.121 0.089 [ 6.776 11.036] Posterior quantiles: 2.5 25 50 75 97.5 |---------------|===============|===============|---------------| 2.07 2.701 3.057 3.508 4.62 1.839 2.516 2.863 3.267 4.089 2.409 3.099 3.527 3.948 4.824 2.79 3.916 4.575 5.352 6.78 4.689 5.799 6.514 7.148 8.982 4.757 5.741 6.26 6.786 7.818 6.323 7.753 8.572 9.331 11.041 6.087 7.278 7.976 8.755 10.811 5.631 6.819 7.548 8.387 10.12 5.989 7.165 7.836 8.632 10.024 5.999 7.263 8.13 9.098 10.664 5.753 7.101 7.915 8.606 10.092 6.896 8.189 8.94 9.769 11.22
M_pneumononPCV.rate_west_under65.summary()
rate_west_under65: Mean SD MC Error 95% HPD interval ------------------------------------------------------------------ 3.86 0.511 0.041 [ 2.83 4.822] 4.2 0.498 0.041 [ 3.258 5.132] 5.015 0.551 0.043 [ 3.903 6.021] 5.643 0.594 0.046 [ 4.574 6.871] 7.987 0.657 0.046 [ 6.744 9.309] 7.543 0.669 0.049 [ 6.156 8.823] 6.56 0.592 0.042 [ 5.451 7.76 ] 6.439 0.564 0.04 [ 5.413 7.628] 5.734 0.529 0.039 [ 4.694 6.805] 6.479 0.6 0.043 [ 5.287 7.595] 6.246 0.56 0.04 [ 5.195 7.395] 6.675 0.589 0.041 [ 5.509 7.796] 5.868 0.515 0.035 [ 4.919 6.972] Posterior quantiles: 2.5 25 50 75 97.5 |---------------|===============|===============|---------------| 2.896 3.531 3.847 4.199 4.932 3.288 3.851 4.185 4.547 5.178 4.006 4.63 5.001 5.372 6.224 4.628 5.217 5.578 6.02 6.976 6.833 7.532 7.927 8.402 9.428 6.159 7.124 7.552 7.977 8.836 5.478 6.15 6.545 6.933 7.823 5.423 6.043 6.417 6.802 7.658 4.671 5.381 5.729 6.075 6.802 5.269 6.083 6.487 6.906 7.584 5.197 5.858 6.224 6.598 7.426 5.531 6.276 6.659 7.067 7.835 4.772 5.525 5.882 6.212 6.851
M_pneumononPCV.rate_east_over64.summary()
rate_east_over64: Mean SD MC Error 95% HPD interval ------------------------------------------------------------------ 11.389 2.662 0.247 [ 6.039 16.589] 13.413 2.4 0.215 [ 8.388 18.061] 11.775 2.837 0.265 [ 6.775 17.42 ] 17.079 3.494 0.318 [ 10.933 23.805] 21.895 3.621 0.327 [ 15.517 29.513] 30.855 3.864 0.338 [ 23.844 39.233] 26.829 3.815 0.337 [ 20.033 33.841] 26.722 3.498 0.302 [ 20.399 33.821] 21.614 2.854 0.245 [ 16.261 27.248] 20.633 3.465 0.315 [ 13.337 26.773] 29.564 3.559 0.306 [ 22.705 36.875] 24.831 3.277 0.28 [ 18.784 31.394] 18.538 2.77 0.242 [ 13.367 23.897] Posterior quantiles: 2.5 25 50 75 97.5 |---------------|===============|===============|---------------| 6.908 9.631 11.039 12.964 17.777 8.826 11.861 13.278 14.953 18.575 6.834 9.637 11.698 13.62 17.589 11.09 14.29 17.037 19.553 24.227 15.709 19.291 21.503 24.23 29.767 22.333 28.404 30.836 33.498 38.128 19.848 23.93 26.835 29.75 33.684 20.273 24.317 26.588 29.111 33.748 16.82 19.483 21.421 23.398 28.21 14.057 18.148 20.761 22.972 27.534 22.673 27.172 29.487 31.826 36.866 18.647 22.478 24.775 27.075 31.342 13.501 16.585 18.397 20.419 24.184
M_pneumononPCV.rate_west_over64.summary()
rate_west_over64: Mean SD MC Error 95% HPD interval ------------------------------------------------------------------ 12.085 2.097 0.187 [ 8.228 16.473] 11.144 2.108 0.192 [ 6.918 15.107] 17.338 2.181 0.185 [ 13.003 21.258] 17.507 2.275 0.198 [ 13.092 21.762] 20.48 2.608 0.225 [ 15.228 25.851] 20.591 2.502 0.209 [ 15.727 25.255] 20.172 2.42 0.204 [ 15.706 24.941] 25.796 2.537 0.202 [ 21.128 31.012] 22.035 2.306 0.18 [ 17.655 26.498] 22.672 2.298 0.185 [ 18.499 27.253] 23.131 2.49 0.198 [ 18.086 27.993] 18.618 2.217 0.187 [ 14.472 22.733] 22.82 2.15 0.169 [ 18.787 27.008] Posterior quantiles: 2.5 25 50 75 97.5 |---------------|===============|===============|---------------| 8.561 10.608 11.864 13.207 17.137 7.315 9.684 11.05 12.517 15.622 13.115 15.786 17.341 18.933 21.419 13.153 15.918 17.509 19.109 21.86 14.762 18.79 20.516 22.155 25.472 16.206 18.857 20.497 22.212 25.813 16.143 18.443 19.882 21.714 25.597 21.068 24.042 25.638 27.52 30.987 17.755 20.375 21.985 23.589 26.631 18.649 20.975 22.566 24.288 27.44 18.321 21.454 23.114 24.72 28.298 14.484 17.066 18.596 20.209 22.781 18.951 21.325 22.706 24.229 27.308
M_pneumononPCV_tenn = pm.MCMC(rate_model('pneumononPCV', regional=False))
M_pneumononPCV_tenn.sample(50000, 40000)
[-----------------100%-----------------] 50000 of 50000 complete in 28.9 sec
pm.Matplot.summary_plot(M_pneumononPCV_tenn.rate_under65)
Could not calculate Gelman-Rubin statistics. Requires multiple chains of equal length.
pm.Matplot.summary_plot(M_pneumononPCV_tenn.rate_over64)
Could not calculate Gelman-Rubin statistics. Requires multiple chains of equal length.
M_pneumononPCV_tenn.rate_under65.summary()
rate_under65: Mean SD MC Error 95% HPD interval ------------------------------------------------------------------ 3.623 0.413 0.033 [ 2.832 4.447] 3.852 0.454 0.037 [ 2.968 4.774] 4.653 0.491 0.038 [ 3.725 5.672] 5.321 0.504 0.037 [ 4.336 6.344] 7.719 0.585 0.038 [ 6.612 8.846] 6.99 0.566 0.039 [ 5.864 8.074] 7.1 0.586 0.041 [ 5.983 8.244] 6.913 0.575 0.04 [ 5.774 8.006] 6.254 0.526 0.037 [ 5.23 7.256] 6.913 0.591 0.043 [ 5.82 8.093] 6.737 0.538 0.037 [ 5.779 7.878] 6.973 0.51 0.035 [ 5.89 7.93] 6.677 0.552 0.039 [ 5.59 7.732] Posterior quantiles: 2.5 25 50 75 97.5 |---------------|===============|===============|---------------| 2.879 3.336 3.586 3.888 4.505 2.935 3.554 3.838 4.156 4.752 3.761 4.325 4.598 4.949 5.725 4.345 4.989 5.309 5.632 6.356 6.62 7.299 7.715 8.134 8.856 5.9 6.608 6.99 7.366 8.133 6.017 6.687 7.086 7.496 8.3 5.821 6.522 6.891 7.314 8.075 5.298 5.88 6.212 6.595 7.343 5.82 6.497 6.894 7.324 8.093 5.714 6.378 6.731 7.099 7.824 5.987 6.634 6.964 7.282 8.062 5.66 6.289 6.647 7.054 7.827
M_pneumononPCV_tenn.rate_over64.summary()
rate_over64: Mean SD MC Error 95% HPD interval ------------------------------------------------------------------ 12.327 1.85 0.158 [ 8.619 15.819] 11.427 1.777 0.153 [ 8.12 15.149] 16.347 1.681 0.132 [ 13.393 19.786] 17.547 2.121 0.176 [ 13.353 21.534] 20.144 2.095 0.164 [ 16.247 24.312] 24.217 2.249 0.169 [ 19.847 28.538] 21.922 2.362 0.191 [ 17.748 26.604] 25.963 2.395 0.187 [ 21.335 30.509] 21.733 1.962 0.148 [ 17.56 25.654] 22.394 2.099 0.158 [ 18.247 26.492] 25.68 2.156 0.158 [ 21.116 29.561] 19.96 1.952 0.149 [ 15.961 23.681] 21.229 2.131 0.166 [ 17.317 25.308] Posterior quantiles: 2.5 25 50 75 97.5 |---------------|===============|===============|---------------| 8.71 11.112 12.304 13.552 16.061 8.212 10.193 11.35 12.565 15.287 13.446 15.135 16.234 17.386 19.94 13.449 16.056 17.524 18.999 21.712 16.12 18.671 20.074 21.594 24.243 19.954 22.652 24.143 25.762 28.764 18.056 20.195 21.657 23.484 27.033 21.251 24.39 25.985 27.607 30.477 17.819 20.433 21.625 22.97 25.989 18.444 20.951 22.369 23.787 26.813 21.623 24.14 25.637 27.126 30.168 16.306 18.687 19.91 21.115 24.22 17.346 19.669 21.211 22.642 25.442
M_levr = pm.MCMC(rate_model('LEVR'))
M_levr.sample(50000, 40000)
# M_levr.sample(50000, 40000)
[-----------------100%-----------------] 50000 of 50000 complete in 34.2 sec
levr_rates = make_rates_df(M_levr)
make_factorplot(levr_rates)
pm.Matplot.summary_plot(M_levr.rate_west_over64)
Could not calculate Gelman-Rubin statistics. Requires multiple chains of equal length.
pm.Matplot.summary_plot(M_levr.rate_west_under65)
Could not calculate Gelman-Rubin statistics. Requires multiple chains of equal length.
pm.Matplot.summary_plot(M_levr.rate_east_over64)
Could not calculate Gelman-Rubin statistics. Requires multiple chains of equal length.
pm.Matplot.summary_plot(M_levr.rate_east_under65)
Could not calculate Gelman-Rubin statistics. Requires multiple chains of equal length.
M_levr.rate_east_under65.summary()
rate_east_under65: Mean SD MC Error 95% HPD interval ------------------------------------------------------------------ 0.025 0.031 0.002 [ 0. 0.083] 0.008 0.014 0.001 [ 0. 0.03] 0.021 0.031 0.002 [ 0. 0.073] 0.01 0.017 0.001 [ 0. 0.036] 0.016 0.023 0.001 [ 0. 0.052] 0.01 0.013 0.001 [ 0. 0.036] 0.013 0.019 0.001 [ 0. 0.046] 0.013 0.018 0.001 [ 0. 0.042] 0.008 0.012 0.001 [ 0. 0.028] 0.01 0.022 0.002 [ 0. 0.05] 0.011 0.018 0.001 [ 0. 0.035] 0.019 0.034 0.002 [ 0. 0.071] 0.008 0.011 0.001 [ 0. 0.028] Posterior quantiles: 2.5 25 50 75 97.5 |---------------|===============|===============|---------------| 0.001 0.007 0.015 0.032 0.113 0.0 0.001 0.003 0.008 0.044 0.001 0.005 0.011 0.024 0.095 0.0 0.002 0.005 0.011 0.056 0.001 0.004 0.009 0.019 0.076 0.0 0.002 0.005 0.012 0.047 0.0 0.002 0.006 0.015 0.068 0.0 0.003 0.007 0.015 0.057 0.0 0.002 0.004 0.009 0.04 0.0 0.001 0.003 0.009 0.073 0.0 0.002 0.006 0.013 0.051 0.001 0.003 0.008 0.021 0.105 0.0 0.002 0.004 0.01 0.039
M_levr.rate_west_under65.summary()
rate_west_under65: Mean SD MC Error 95% HPD interval ------------------------------------------------------------------ 0.047 0.04 0.003 [ 0.001 0.123] 0.03 0.029 0.002 [ 0. 0.091] 0.069 0.049 0.004 [ 0.006 0.166] 0.032 0.028 0.002 [ 0.001 0.089] 0.099 0.073 0.006 [ 0.006 0.247] 0.036 0.036 0.003 [ 0.001 0.098] 0.062 0.048 0.004 [ 0.005 0.163] 0.122 0.079 0.006 [ 0.015 0.272] 0.033 0.03 0.003 [ 0.001 0.096] 0.027 0.024 0.002 [ 0. 0.076] 0.032 0.034 0.003 [ 0. 0.094] 0.08 0.058 0.004 [ 0.008 0.198] 0.047 0.035 0.002 [ 0.001 0.113] Posterior quantiles: 2.5 25 50 75 97.5 |---------------|===============|===============|---------------| 0.005 0.018 0.034 0.064 0.146 0.001 0.01 0.021 0.042 0.112 0.012 0.032 0.058 0.094 0.2 0.001 0.011 0.024 0.044 0.105 0.015 0.049 0.08 0.131 0.29 0.004 0.012 0.025 0.047 0.124 0.012 0.029 0.049 0.077 0.201 0.025 0.064 0.1 0.162 0.317 0.003 0.012 0.023 0.043 0.117 0.001 0.007 0.022 0.038 0.09 0.002 0.009 0.023 0.042 0.129 0.015 0.041 0.064 0.103 0.24 0.006 0.022 0.039 0.064 0.138
M_levr.rate_east_over64.summary()
rate_east_over64: Mean SD MC Error 95% HPD interval ------------------------------------------------------------------ 0.141 0.188 0.012 [ 0.001 0.504] 0.064 0.095 0.007 [ 0. 0.258] 0.106 0.153 0.009 [ 0. 0.331] 0.051 0.064 0.004 [ 0. 0.171] 0.116 0.194 0.013 [ 0.001 0.412] 0.073 0.095 0.006 [ 0. 0.252] 0.103 0.145 0.01 [ 0. 0.385] 0.063 0.076 0.004 [ 0.001 0.207] 0.07 0.108 0.008 [ 0. 0.266] 0.042 0.06 0.004 [ 0. 0.159] 0.083 0.113 0.007 [ 0. 0.292] 0.176 0.205 0.011 [ 0.001 0.546] 0.057 0.085 0.006 [ 0. 0.207] Posterior quantiles: 2.5 25 50 75 97.5 |---------------|===============|===============|---------------| 0.008 0.038 0.079 0.164 0.689 0.0 0.008 0.03 0.08 0.328 0.004 0.027 0.06 0.129 0.481 0.001 0.012 0.029 0.063 0.221 0.003 0.021 0.056 0.131 0.623 0.002 0.016 0.041 0.091 0.32 0.003 0.021 0.052 0.119 0.527 0.003 0.017 0.038 0.077 0.27 0.001 0.009 0.032 0.084 0.377 0.001 0.006 0.02 0.052 0.212 0.002 0.017 0.045 0.101 0.385 0.008 0.048 0.112 0.227 0.733 0.001 0.01 0.03 0.072 0.269
M_levr.rate_west_over64.summary()
rate_west_over64: Mean SD MC Error 95% HPD interval ------------------------------------------------------------------ 0.889 0.517 0.031 [ 0.133 1.904] 0.195 0.174 0.013 [ 0.004 0.543] 0.919 0.562 0.036 [ 0.109 2.028] 0.177 0.162 0.011 [ 0.007 0.46 ] 0.374 0.303 0.021 [ 0.009 0.969] 0.375 0.275 0.017 [ 0.015 0.896] 0.645 0.417 0.025 [ 0.081 1.492] 0.28 0.241 0.018 [ 0.003 0.764] 0.371 0.275 0.019 [ 0.003 0.89 ] 0.178 0.165 0.012 [ 0.001 0.499] 0.538 0.342 0.022 [ 0.064 1.205] 0.314 0.25 0.018 [ 0.012 0.795] 0.192 0.166 0.013 [ 0.002 0.536] Posterior quantiles: 2.5 25 50 75 97.5 |---------------|===============|===============|---------------| 0.206 0.514 0.781 1.152 2.18 0.009 0.064 0.144 0.279 0.639 0.207 0.518 0.797 1.178 2.366 0.02 0.07 0.132 0.232 0.603 0.043 0.167 0.294 0.485 1.198 0.055 0.181 0.306 0.496 1.073 0.13 0.342 0.546 0.837 1.715 0.02 0.104 0.21 0.384 0.913 0.028 0.175 0.314 0.496 1.076 0.01 0.06 0.128 0.25 0.618 0.12 0.302 0.459 0.692 1.474 0.044 0.137 0.249 0.419 0.964 0.009 0.072 0.149 0.26 0.613
M_levr_tenn = pm.MCMC(rate_model('LEVR', regional=False))
M_levr_tenn.sample(50000, 40000)
[-----------------100%-----------------] 50000 of 50000 complete in 32.9 sec
pm.Matplot.summary_plot(M_levr_tenn.rate_under65)
Could not calculate Gelman-Rubin statistics. Requires multiple chains of equal length.
pm.Matplot.summary_plot(M_levr_tenn.rate_over64)
Could not calculate Gelman-Rubin statistics. Requires multiple chains of equal length.
M_levr_tenn.rate_under65.summary()
rate_under65: Mean SD MC Error 95% HPD interval ------------------------------------------------------------------ 0.048 0.034 0.003 [ 0.005 0.11 ] 0.028 0.021 0.002 [ 0.003 0.074] 0.056 0.037 0.003 [ 0.005 0.13 ] 0.028 0.021 0.002 [ 0.001 0.07 ] 0.072 0.04 0.003 [ 0.011 0.149] 0.035 0.024 0.002 [ 0.002 0.08 ] 0.048 0.034 0.003 [ 0.002 0.117] 0.093 0.058 0.004 [ 0.017 0.211] 0.03 0.019 0.001 [ 0.003 0.069] 0.028 0.023 0.002 [ 0.001 0.076] 0.031 0.023 0.002 [ 0.003 0.076] 0.064 0.038 0.003 [ 0.011 0.136] 0.041 0.026 0.002 [ 0.004 0.093] Posterior quantiles: 2.5 25 50 75 97.5 |---------------|===============|===============|---------------| 0.01 0.026 0.039 0.06 0.133 0.006 0.014 0.021 0.035 0.087 0.011 0.03 0.048 0.072 0.154 0.004 0.013 0.022 0.037 0.083 0.017 0.043 0.064 0.095 0.167 0.006 0.018 0.029 0.047 0.097 0.007 0.022 0.038 0.063 0.137 0.023 0.053 0.08 0.118 0.243 0.005 0.017 0.026 0.04 0.077 0.004 0.013 0.021 0.035 0.091 0.005 0.015 0.024 0.04 0.094 0.016 0.037 0.055 0.082 0.157 0.009 0.021 0.034 0.053 0.11
M_levr_tenn.rate_over64.summary()
rate_over64: Mean SD MC Error 95% HPD interval ------------------------------------------------------------------ 0.622 0.348 0.023 [ 0.105 1.282] 0.187 0.149 0.011 [ 0.016 0.466] 0.616 0.372 0.023 [ 0.075 1.357] 0.161 0.128 0.009 [ 0.006 0.407] 0.34 0.232 0.016 [ 0.021 0.779] 0.267 0.189 0.012 [ 0.016 0.627] 0.41 0.264 0.016 [ 0.027 0.93 ] 0.226 0.17 0.011 [ 0.017 0.563] 0.294 0.204 0.014 [ 0.011 0.701] 0.178 0.144 0.01 [ 0.011 0.466] 0.348 0.231 0.015 [ 0.039 0.784] 0.315 0.205 0.013 [ 0.025 0.712] 0.194 0.146 0.01 [ 0.007 0.483] Posterior quantiles: 2.5 25 50 75 97.5 |---------------|===============|===============|---------------| 0.164 0.366 0.557 0.796 1.482 0.032 0.092 0.147 0.233 0.592 0.15 0.343 0.526 0.799 1.559 0.02 0.07 0.126 0.217 0.483 0.062 0.179 0.288 0.432 0.932 0.052 0.136 0.22 0.346 0.771 0.084 0.228 0.356 0.52 1.106 0.038 0.103 0.181 0.3 0.671 0.054 0.151 0.239 0.379 0.831 0.024 0.079 0.139 0.234 0.581 0.073 0.187 0.294 0.45 0.92 0.05 0.159 0.274 0.425 0.818 0.021 0.097 0.157 0.248 0.611
M_penr = pm.MCMC(rate_model('PENR'))
M_penr.sample(50000, 40000)
# M_penr.sample(50000, 40000)
[-----------------100%-----------------] 50000 of 50000 complete in 31.9 sec
penr_rates = make_rates_df(M_penr)
make_factorplot(penr_rates)
pm.Matplot.summary_plot(M_penr.rate_west_over64)
Could not calculate Gelman-Rubin statistics. Requires multiple chains of equal length.
pm.Matplot.summary_plot(M_penr.rate_west_under65)
Could not calculate Gelman-Rubin statistics. Requires multiple chains of equal length.
pm.Matplot.summary_plot(M_penr.rate_east_over64)
Could not calculate Gelman-Rubin statistics. Requires multiple chains of equal length.
pm.Matplot.summary_plot(M_penr.rate_east_under65)
Could not calculate Gelman-Rubin statistics. Requires multiple chains of equal length.
M_penr.rate_east_under65.summary()
rate_east_under65: Mean SD MC Error 95% HPD interval ------------------------------------------------------------------ 2.683 0.659 0.049 [ 1.482 4.011] 0.564 0.235 0.019 [ 0.139 0.974] 0.758 0.273 0.022 [ 0.281 1.303] 0.795 0.333 0.028 [ 0.224 1.446] 1.284 0.422 0.033 [ 0.575 2.12 ] 0.945 0.389 0.033 [ 0.288 1.752] 1.324 0.392 0.029 [ 0.662 2.111] 1.636 0.464 0.036 [ 0.777 2.513] 1.851 0.484 0.035 [ 0.884 2.756] 1.396 0.409 0.032 [ 0.662 2.173] 1.661 0.451 0.034 [ 0.813 2.511] 0.532 0.19 0.015 [ 0.207 0.92 ] 0.555 0.249 0.021 [ 0.164 1.061] Posterior quantiles: 2.5 25 50 75 97.5 |---------------|===============|===============|---------------| 1.531 2.187 2.663 3.09 4.081 0.207 0.395 0.533 0.697 1.124 0.336 0.571 0.714 0.894 1.453 0.272 0.532 0.771 0.992 1.54 0.63 0.985 1.232 1.521 2.295 0.358 0.665 0.871 1.168 1.877 0.732 1.029 1.264 1.577 2.211 0.855 1.299 1.594 1.931 2.637 1.019 1.511 1.812 2.148 2.936 0.73 1.089 1.35 1.664 2.298 0.914 1.345 1.614 1.931 2.69 0.236 0.394 0.498 0.637 0.987 0.2 0.371 0.516 0.693 1.169
M_penr.rate_west_under65.summary()
rate_west_under65: Mean SD MC Error 95% HPD interval ------------------------------------------------------------------ 1.208 0.25 0.018 [ 0.756 1.703] 0.482 0.145 0.012 [ 0.253 0.794] 0.603 0.176 0.014 [ 0.296 0.956] 0.573 0.166 0.013 [ 0.284 0.906] 0.803 0.189 0.013 [ 0.479 1.186] 0.83 0.197 0.015 [ 0.457 1.208] 0.8 0.179 0.013 [ 0.489 1.147] 1.039 0.22 0.016 [ 0.617 1.477] 1.123 0.234 0.017 [ 0.683 1.578] 1.304 0.265 0.016 [ 0.845 1.862] 1.35 0.268 0.018 [ 0.838 1.894] 0.342 0.119 0.01 [ 0.151 0.58 ] 0.257 0.088 0.007 [ 0.113 0.437] Posterior quantiles: 2.5 25 50 75 97.5 |---------------|===============|===============|---------------| 0.779 1.028 1.187 1.372 1.737 0.256 0.374 0.463 0.568 0.807 0.314 0.477 0.585 0.711 0.989 0.302 0.454 0.562 0.67 0.936 0.501 0.665 0.784 0.916 1.231 0.474 0.689 0.818 0.964 1.239 0.507 0.672 0.783 0.909 1.179 0.664 0.886 1.017 1.164 1.544 0.696 0.958 1.11 1.277 1.593 0.846 1.109 1.289 1.476 1.865 0.878 1.168 1.321 1.521 1.952 0.165 0.251 0.328 0.41 0.611 0.127 0.192 0.242 0.307 0.475
M_penr.rate_east_over64.summary()
rate_east_over64: Mean SD MC Error 95% HPD interval ------------------------------------------------------------------ 9.18 2.409 0.187 [ 4.957 13.942] 2.246 1.053 0.091 [ 0.61 4.201] 3.169 1.308 0.111 [ 0.968 5.76 ] 2.414 0.982 0.085 [ 1.014 4.515] 6.489 1.857 0.139 [ 3.229 10.003] 2.942 1.235 0.107 [ 0.918 5.471] 3.061 1.261 0.11 [ 1.119 5.506] 5.255 1.897 0.165 [ 2.146 9.381] 4.625 1.493 0.121 [ 1.895 7.414] 4.828 1.639 0.136 [ 1.903 7.777] 5.926 1.67 0.129 [ 2.889 8.983] 1.371 0.613 0.053 [ 0.416 2.588] 1.578 0.717 0.062 [ 0.376 2.857] Posterior quantiles: 2.5 25 50 75 97.5 |---------------|===============|===============|---------------| 5.151 7.351 8.936 10.834 14.28 0.762 1.501 2.072 2.794 4.75 1.181 2.253 2.939 3.898 6.288 1.114 1.654 2.211 2.966 4.816 3.536 5.114 6.297 7.7 10.559 1.177 2.058 2.718 3.619 5.923 1.322 2.022 2.889 3.836 6.038 2.291 3.833 5.084 6.449 9.621 2.115 3.561 4.522 5.561 7.823 2.289 3.583 4.646 5.828 8.545 3.212 4.708 5.766 7.021 9.429 0.513 0.901 1.282 1.698 2.806 0.58 1.057 1.452 1.962 3.357
M_penr.rate_west_over64.summary()
rate_west_over64: Mean SD MC Error 95% HPD interval ------------------------------------------------------------------ 7.035 1.579 0.113 [ 4.182 10.259] 1.921 0.664 0.052 [ 0.776 3.208] 2.566 0.847 0.069 [ 1.22 4.386] 2.149 0.68 0.055 [ 0.883 3.483] 2.539 0.804 0.065 [ 1.153 4.109] 2.064 0.731 0.057 [ 0.838 3.372] 3.641 0.959 0.069 [ 2.018 5.574] 3.643 1.026 0.081 [ 1.815 5.649] 2.267 0.837 0.072 [ 0.746 3.894] 2.816 0.783 0.06 [ 1.445 4.389] 3.543 0.994 0.073 [ 1.785 5.527] 1.08 0.4 0.033 [ 0.424 1.895] 1.011 0.403 0.033 [ 0.366 1.781] Posterior quantiles: 2.5 25 50 75 97.5 |---------------|===============|===============|---------------| 4.28 5.932 6.928 8.031 10.529 0.903 1.417 1.821 2.351 3.436 1.307 1.96 2.422 3.024 4.582 1.025 1.663 2.044 2.571 3.687 1.255 1.936 2.45 3.046 4.314 1.006 1.572 1.945 2.419 3.831 2.088 2.914 3.536 4.279 5.729 1.902 2.886 3.593 4.301 5.869 0.93 1.649 2.172 2.755 4.208 1.543 2.241 2.725 3.315 4.528 1.903 2.807 3.436 4.211 5.698 0.485 0.791 1.014 1.296 2.036 0.449 0.721 0.936 1.231 1.998
M_penr_tenn = pm.MCMC(rate_model('PENR', regional=False))
M_penr_tenn.sample(50000, 40000)
[-----------------100%-----------------] 50000 of 50000 complete in 35.6 sec
pm.Matplot.summary_plot(M_penr_tenn.rate_under65)
Could not calculate Gelman-Rubin statistics. Requires multiple chains of equal length.
pm.Matplot.summary_plot(M_penr_tenn.rate_over64)
Could not calculate Gelman-Rubin statistics. Requires multiple chains of equal length.
M_penr_tenn.rate_under65.summary()
rate_under65: Mean SD MC Error 95% HPD interval ------------------------------------------------------------------ 1.484 0.266 0.017 [ 0.96 1.985] 0.483 0.133 0.01 [ 0.25 0.749] 0.614 0.155 0.012 [ 0.326 0.927] 0.618 0.177 0.014 [ 0.284 0.954] 0.943 0.205 0.014 [ 0.571 1.361] 0.897 0.197 0.014 [ 0.53 1.307] 0.971 0.212 0.014 [ 0.537 1.352] 1.187 0.23 0.015 [ 0.764 1.639] 1.359 0.249 0.015 [ 0.901 1.863] 1.344 0.246 0.017 [ 0.895 1.845] 1.432 0.242 0.015 [ 0.968 1.9 ] 0.383 0.11 0.009 [ 0.179 0.591] 0.343 0.127 0.011 [ 0.139 0.598] Posterior quantiles: 2.5 25 50 75 97.5 |---------------|===============|===============|---------------| 1.005 1.291 1.477 1.652 2.053 0.27 0.388 0.465 0.566 0.775 0.357 0.5 0.597 0.704 0.972 0.32 0.495 0.605 0.721 1.022 0.595 0.795 0.919 1.079 1.39 0.561 0.758 0.879 1.012 1.35 0.587 0.822 0.965 1.104 1.409 0.787 1.025 1.167 1.34 1.672 0.921 1.179 1.35 1.516 1.909 0.904 1.177 1.319 1.498 1.871 1.011 1.264 1.412 1.586 1.957 0.194 0.302 0.376 0.45 0.621 0.155 0.245 0.319 0.426 0.629
M_penr_tenn.rate_over64.summary()
rate_over64: Mean SD MC Error 95% HPD interval ------------------------------------------------------------------ 8.196 1.524 0.102 [ 5.435 11.498] 2.067 0.658 0.052 [ 0.8 3.26] 2.989 0.816 0.059 [ 1.5 4.588] 2.364 0.725 0.056 [ 1.1 3.801] 3.571 0.93 0.067 [ 1.844 5.476] 2.177 0.623 0.048 [ 1.04 3.435] 3.53 0.873 0.056 [ 1.908 5.244] 4.028 0.967 0.07 [ 2.395 6.021] 2.437 0.732 0.056 [ 1.217 4.028] 3.389 0.832 0.058 [ 1.879 4.991] 4.204 0.909 0.058 [ 2.38 5.961] 1.218 0.459 0.038 [ 0.43 2.181] 1.238 0.443 0.036 [ 0.471 2.08 ] Posterior quantiles: 2.5 25 50 75 97.5 |---------------|===============|===============|---------------| 5.517 7.137 8.091 9.109 11.67 0.997 1.603 1.984 2.445 3.581 1.585 2.427 2.95 3.493 4.771 1.259 1.851 2.29 2.733 4.13 1.915 2.91 3.474 4.17 5.571 1.152 1.747 2.109 2.547 3.627 2.006 2.911 3.445 4.083 5.395 2.55 3.349 3.885 4.564 6.314 1.218 1.91 2.372 2.9 4.031 1.991 2.8 3.326 3.9 5.234 2.521 3.591 4.155 4.767 6.158 0.442 0.89 1.16 1.498 2.209 0.547 0.913 1.18 1.512 2.235
M_erythr = pm.MCMC(rate_model('ERYTHR'))
M_erythr.sample(50000, 40000)
# M_erythr.sample(50000, 40000)
[-----------------100%-----------------] 50000 of 50000 complete in 41.0 sec
erythr_rates = make_rates_df(M_erythr)
make_factorplot(erythr_rates)
pm.Matplot.summary_plot(M_erythr.rate_west_over64)
Could not calculate Gelman-Rubin statistics. Requires multiple chains of equal length.
pm.Matplot.summary_plot(M_erythr.rate_west_under65)
Could not calculate Gelman-Rubin statistics. Requires multiple chains of equal length.
pm.Matplot.summary_plot(M_erythr.rate_east_over64)
Could not calculate Gelman-Rubin statistics. Requires multiple chains of equal length.
pm.Matplot.summary_plot(M_erythr.rate_east_under65)
Could not calculate Gelman-Rubin statistics. Requires multiple chains of equal length.
M_erythr.rate_east_under65.summary()
rate_east_under65: Mean SD MC Error 95% HPD interval ------------------------------------------------------------------ 3.471 0.767 0.067 [ 2.137 5.076] 2.315 0.545 0.048 [ 1.32 3.384] 2.155 0.581 0.052 [ 1.104 3.321] 1.893 0.436 0.038 [ 1.192 2.803] 3.526 0.726 0.06 [ 2.207 4.954] 3.188 0.589 0.05 [ 2.222 4.498] 3.211 0.718 0.062 [ 1.938 4.742] 4.048 0.855 0.074 [ 2.324 5.66 ] 4.66 0.797 0.065 [ 3.238 6.219] 4.484 0.832 0.07 [ 2.945 6.135] 5.041 0.827 0.065 [ 3.505 6.682] 2.893 0.611 0.053 [ 1.82 4.086] 3.083 0.641 0.055 [ 1.989 4.491] Posterior quantiles: 2.5 25 50 75 97.5 |---------------|===============|===============|---------------| 2.218 2.935 3.351 3.905 5.236 1.391 1.931 2.253 2.652 3.497 1.207 1.728 2.088 2.51 3.499 1.214 1.56 1.848 2.171 2.849 2.333 3.012 3.437 3.939 5.199 2.166 2.752 3.146 3.592 4.457 1.807 2.737 3.201 3.669 4.637 2.474 3.457 4.002 4.616 5.836 3.274 4.064 4.633 5.197 6.314 3.08 3.866 4.397 4.998 6.377 3.549 4.439 5.0 5.59 6.763 1.939 2.425 2.831 3.298 4.245 2.065 2.629 2.987 3.456 4.621
M_erythr.rate_west_under65.summary()
rate_west_under65: Mean SD MC Error 95% HPD interval ------------------------------------------------------------------ 2.316 0.374 0.03 [ 1.562 2.984] 1.913 0.336 0.028 [ 1.315 2.61 ] 2.604 0.392 0.032 [ 1.877 3.38 ] 2.191 0.35 0.029 [ 1.604 2.912] 3.143 0.41 0.03 [ 2.425 3.92 ] 3.865 0.429 0.031 [ 2.932 4.633] 3.166 0.454 0.037 [ 2.273 4.041] 2.69 0.369 0.029 [ 2.033 3.487] 3.123 0.402 0.03 [ 2.371 3.908] 3.369 0.439 0.033 [ 2.521 4.214] 3.23 0.439 0.033 [ 2.412 4.075] 3.055 0.398 0.03 [ 2.235 3.793] 2.749 0.38 0.03 [ 1.997 3.485] Posterior quantiles: 2.5 25 50 75 97.5 |---------------|===============|===============|---------------| 1.589 2.051 2.307 2.583 3.026 1.361 1.674 1.868 2.112 2.677 1.929 2.335 2.574 2.851 3.452 1.634 1.945 2.155 2.39 2.972 2.434 2.84 3.129 3.424 3.947 3.005 3.585 3.859 4.148 4.741 2.306 2.85 3.152 3.464 4.093 1.99 2.432 2.669 2.935 3.456 2.422 2.824 3.11 3.38 3.968 2.573 3.067 3.339 3.642 4.315 2.469 2.908 3.197 3.526 4.16 2.306 2.772 3.036 3.335 3.883 2.042 2.494 2.741 2.987 3.572
M_erythr.rate_east_over64.summary()
rate_east_over64: Mean SD MC Error 95% HPD interval ------------------------------------------------------------------ 13.341 2.575 0.221 [ 8.555 18.297] 6.673 1.787 0.161 [ 3.461 10.128] 11.744 2.888 0.265 [ 7.076 17.908] 6.91 1.589 0.138 [ 3.909 10.047] 14.014 3.074 0.276 [ 7.544 20.147] 10.243 2.124 0.188 [ 6.464 14.547] 14.281 2.856 0.253 [ 9.266 19.662] 13.908 2.471 0.21 [ 9.564 19.319] 12.332 2.154 0.178 [ 8.664 17.256] 14.091 2.55 0.218 [ 9.764 19.21 ] 14.966 2.726 0.237 [ 9.639 19.99 ] 11.039 2.426 0.215 [ 6.278 16.203] 8.701 2.078 0.186 [ 5.229 12.607] Posterior quantiles: 2.5 25 50 75 97.5 |---------------|===============|===============|---------------| 9.054 11.435 13.079 15.027 18.976 3.61 5.371 6.621 7.794 10.428 7.379 9.734 11.156 13.405 18.381 4.283 5.75 6.697 7.932 10.501 7.687 12.028 13.901 16.027 20.406 6.643 8.664 10.066 11.619 14.777 9.617 12.073 14.039 16.26 20.168 9.381 12.118 13.748 15.496 19.226 8.178 10.781 12.227 13.787 16.825 9.869 12.137 13.821 15.832 19.387 9.909 13.115 14.897 16.709 20.594 6.434 9.465 10.811 12.52 16.434 5.485 7.104 8.424 10.022 13.181
M_pneumo7.rate_west_over64.summary()
rate_west_over64: Mean SD MC Error 95% HPD interval ------------------------------------------------------------------ 19.778 2.447 0.123 [ 15.132 24.547] 13.049 2.175 0.128 [ 8.75 17.372] 11.408 2.014 0.123 [ 7.434 15.313] 5.277 1.343 0.092 [ 2.897 7.964] 4.082 1.098 0.074 [ 2.171 6.4 ] 2.688 0.849 0.062 [ 1.224 4.431] 2.375 0.758 0.057 [ 1.063 3.941] 1.584 0.55 0.041 [ 0.649 2.67 ] 1.443 0.603 0.048 [ 0.386 2.566] 1.098 0.511 0.043 [ 0.289 2.086] 0.958 0.446 0.037 [ 0.282 1.829] 0.539 0.242 0.019 [ 0.165 1.039] 0.658 0.307 0.025 [ 0.175 1.234] Posterior quantiles: 2.5 25 50 75 97.5 |---------------|===============|===============|---------------| 15.236 18.04 19.715 21.466 24.676 8.91 11.568 12.955 14.426 17.565 7.749 10.003 11.342 12.668 15.722 3.096 4.295 5.132 6.077 8.285 2.199 3.3 4.011 4.771 6.473 1.305 2.076 2.568 3.209 4.573 1.122 1.829 2.294 2.853 4.063 0.746 1.173 1.509 1.914 2.841 0.508 1.008 1.392 1.764 2.842 0.395 0.742 0.988 1.367 2.323 0.361 0.623 0.872 1.194 2.006 0.201 0.362 0.495 0.664 1.151 0.221 0.418 0.609 0.844 1.37
M_erythr_tenn = pm.MCMC(rate_model('ERYTHR', regional=False))
M_erythr_tenn.sample(50000, 40000)
[-----------------100%-----------------] 50000 of 50000 complete in 29.7 sec
pm.Matplot.plot(M_erythr_tenn.beta)
Plotting beta_0
pm.Matplot.summary_plot(M_erythr_tenn.rate_under65)
Could not calculate Gelman-Rubin statistics. Requires multiple chains of equal length.
pm.Matplot.summary_plot(M_erythr_tenn.rate_over64)
Could not calculate Gelman-Rubin statistics. Requires multiple chains of equal length.
M_erythr_tenn.rate_under65.summary()
rate_under65: Mean SD MC Error 95% HPD interval ------------------------------------------------------------------ 2.362 0.333 0.027 [ 1.716 2.992] 2.014 0.319 0.026 [ 1.348 2.584] 2.393 0.305 0.022 [ 1.8 2.985] 1.921 0.292 0.024 [ 1.325 2.48 ] 3.281 0.348 0.024 [ 2.59 3.973] 3.778 0.401 0.028 [ 3.049 4.574] 3.121 0.382 0.028 [ 2.369 3.874] 2.963 0.368 0.028 [ 2.228 3.654] 3.642 0.416 0.031 [ 2.836 4.417] 3.739 0.427 0.03 [ 2.928 4.612] 3.596 0.409 0.028 [ 2.863 4.442] 3.067 0.356 0.027 [ 2.425 3.79 ] 2.839 0.361 0.027 [ 2.134 3.522] Posterior quantiles: 2.5 25 50 75 97.5 |---------------|===============|===============|---------------| 1.732 2.136 2.344 2.592 3.025 1.392 1.787 2.017 2.231 2.64 1.831 2.184 2.377 2.588 3.035 1.326 1.731 1.913 2.121 2.486 2.597 3.048 3.283 3.509 3.983 3.041 3.498 3.762 4.048 4.573 2.444 2.864 3.09 3.347 3.984 2.278 2.703 2.954 3.193 3.741 2.897 3.33 3.616 3.933 4.506 2.925 3.451 3.722 4.018 4.611 2.868 3.314 3.566 3.842 4.449 2.464 2.815 3.034 3.281 3.869 2.179 2.587 2.825 3.069 3.6
M_erythr_tenn.rate_over64.summary()
rate_over64: Mean SD MC Error 95% HPD interval ------------------------------------------------------------------ 13.854 1.682 0.123 [ 10.836 17.413] 6.4 1.263 0.109 [ 4.305 8.929] 10.151 1.581 0.126 [ 7.47 13.482] 8.333 1.418 0.117 [ 5.73 11.17] 11.296 1.647 0.13 [ 8.366 14.781] 8.961 1.337 0.108 [ 6.432 11.554] 11.67 1.66 0.131 [ 8.479 14.781] 11.46 1.562 0.12 [ 8.445 14.499] 8.646 1.147 0.084 [ 6.66 11.159] 10.901 1.389 0.102 [ 8.335 13.582] 13.165 1.767 0.137 [ 9.794 16.677] 8.799 1.349 0.105 [ 6.371 11.524] 10.295 1.416 0.108 [ 7.627 13.097] Posterior quantiles: 2.5 25 50 75 97.5 |---------------|===============|===============|---------------| 10.773 12.669 13.766 14.933 17.399 4.271 5.424 6.337 7.303 8.903 7.099 9.024 10.105 11.312 13.267 5.849 7.31 8.247 9.263 11.356 8.343 10.149 11.17 12.363 14.772 6.489 8.002 8.889 9.911 11.633 8.574 10.473 11.631 12.821 14.955 8.44 10.415 11.477 12.494 14.495 6.665 7.862 8.521 9.337 11.183 8.447 9.867 10.873 11.811 13.768 9.982 11.924 13.005 14.265 17.011 6.308 7.839 8.723 9.709 11.485 7.748 9.292 10.231 11.263 13.248