Commit 595ea2eb authored by Chris Jewell's avatar Chris Jewell
Browse files

As of analysis 2020-11-30

CHANGES:

1. Re-introduced Tier lockdown covariates.  National lockdown covariate is set to 0 because
   we cannot identify this from the time-varying baseline force of infection.

2. Added inference on variance parameter for the temporal Gaussian process on beta.
parent 548066ad
......@@ -2,9 +2,9 @@
# Enqueues COVID-19 pipelines
CASES_FILE="data/Anonymised Combined Line List 20201123.csv"
DATE_LOW="2020-08-28"
DATE_HIGH="2020-11-20"
CASES_FILE="data/Anonymised Combined Line List 20201130.csv"
DATE_LOW="2020-09-04"
DATE_HIGH="2020-11-27"
TEMPLATE_CONFIG=template_config.yaml
......
......@@ -109,9 +109,10 @@ if __name__ == "__main__":
beta2=block0[0],
gamma0=block0[1],
gamma1=block0[2],
sigma=block0[3],
beta3=tf.concat([block0[4:6], [0.0]], axis=-1),
beta1=block1[0],
beta3=[0.0, 0.0, 0.0], # block1[1:4],
xi=block1[1:], # block1[4:],
xi=block1[1:],
seir=events,
)
)
......@@ -135,8 +136,13 @@ if __name__ == "__main__":
covariance_burnin=200,
),
bijector=tfp.bijectors.Blockwise(
bijectors=[tfp.bijectors.Exp(), tfp.bijectors.Identity()],
block_sizes=[1, 2],
bijectors=[
tfp.bijectors.Exp(),
tfp.bijectors.Identity(),
tfp.bijectors.Exp(),
tfp.bijectors.Identity(),
],
block_sizes=[1, 2, 1, 2],
),
name=name,
)
......@@ -290,9 +296,9 @@ if __name__ == "__main__":
tf.random.set_seed(2)
current_state = [
np.array([0.2, 0.0, 0.0], dtype=DTYPE),
np.array([0.2, 0.0, 0.0, 0.1, 0.0, 0.0], dtype=DTYPE),
np.zeros(
model.model["xi"](0.0).event_shape[-1]
model.model["xi"](0.0, 0.1).event_shape[-1]
# + model.model["beta3"]().event_shape[-1]
+ 1,
dtype=DTYPE,
......@@ -312,6 +318,8 @@ if __name__ == "__main__":
"beta2": (samples[0][:, 0], (NUM_BURST_SAMPLES,)),
"gamma0": (samples[0][:, 1], (NUM_BURST_SAMPLES,)),
"gamma1": (samples[0][:, 2], (NUM_BURST_SAMPLES,)),
"sigma": (samples[0][:, 3], (NUM_BURST_SAMPLES,)),
"beta3": (samples[0][:, 4:], (NUM_BURST_SAMPLES, 2)),
"beta1": (samples[1][:, 0], (NUM_BURST_SAMPLES,)),
"xi": (
samples[1][:, 1:],
......@@ -347,6 +355,8 @@ if __name__ == "__main__":
"beta2": samples[0][:, 0],
"gamma0": samples[0][:, 1],
"gamma1": samples[0][:, 2],
"sigma": samples[0][:, 3],
"beta3": samples[0][:, 4:],
"beta1": samples[1][:, 0],
"xi": samples[1][:, 1:],
"events": samples[2],
......
......@@ -73,6 +73,16 @@ def impute_censored_events(cases):
return tf.stack([se_events, ei_events, ir_events], axis=-1)
def conditional_gp(gp, observations, new_index_points):
param = gp.parameters
param["observation_index_points"] = param["index_points"]
param["observations"] = observations
param["index_points"] = new_index_points
return tfd.GaussianProcessRegressionModel(**param)
def CovidUK(covariates, initial_state, initial_step, num_steps, priors):
def beta1():
return tfd.Normal(
......@@ -95,12 +105,18 @@ def CovidUK(covariates, initial_state, initial_step, num_steps, priors):
sample_shape=covariates["L"].shape[-1],
)
def xi(beta1):
sigma = tf.constant(0.4, dtype=DTYPE)
def sigma():
return tfd.Gamma(
concentration=tf.constant(2.0, dtype=DTYPE),
rate=tf.constant(20.0, dtype=DTYPE),
)
def xi(beta1, sigma):
# sigma = tf.constant(0.1, dtype=DTYPE)
phi = tf.constant(24.0, dtype=DTYPE)
kernel = tfp.math.psd_kernels.MaternThreeHalves(sigma, phi)
idx_pts = tf.cast(tf.range(num_steps // XI_FREQ) * XI_FREQ, dtype=DTYPE)
return tfd.GaussianProcess(
return tfd.GaussianProcessRegressionModel(
kernel,
mean_fn=lambda idx: beta1,
index_points=idx_pts[:, tf.newaxis],
......@@ -132,7 +148,7 @@ def CovidUK(covariates, initial_state, initial_step, num_steps, priors):
gamma1 = tf.convert_to_tensor(gamma1, DTYPE)
L = tf.convert_to_tensor(covariates["L"], DTYPE)
L = L - tf.reduce_mean(L, axis=0)
L = L - tf.reduce_mean(L, axis=(0, 1))
weekday = tf.convert_to_tensor(covariates["weekday"], DTYPE)
weekday = weekday - tf.reduce_mean(weekday, axis=-1)
......@@ -198,6 +214,7 @@ def CovidUK(covariates, initial_state, initial_step, num_steps, priors):
beta1=beta1,
beta2=beta2,
beta3=beta3,
sigma=sigma,
xi=xi,
gamma0=gamma0,
gamma1=gamma1,
......@@ -222,14 +239,14 @@ def next_generation_matrix_fn(covar_data, param):
def fn(t, state):
L = tf.convert_to_tensor(covar_data["L"], DTYPE)
L = L - tf.reduce_mean(L, axis=0)
L = L - tf.reduce_mean(L, axis=(0, 1))
C = tf.convert_to_tensor(covar_data["C"], dtype=DTYPE)
C = tf.linalg.set_diag(C, -tf.reduce_sum(C, axis=-2))
C = tf.linalg.set_diag(C, tf.zeros(C.shape[0], dtype=DTYPE))
Cstar = C + tf.transpose(C)
Cstar = tf.linalg.set_diag(Cstar, -tf.reduce_sum(C, axis=-2))
W = tf.constant(covar_data["W"], dtype=DTYPE)
N = tf.constant(covar_data["N"], dtype=DTYPE)
......@@ -242,7 +259,7 @@ def next_generation_matrix_fn(covar_data, param):
xi = tf.gather(param["xi"], xi_idx)
L_idx = tf.clip_by_value(tf.cast(t, tf.int64), 0, L.shape[0] - 1)
Lt = tf.gather(L, L_idx)
Lt = L[-1] # Last timepoint
xB = tf.linalg.matvec(Lt, param["beta3"])
beta = tf.math.exp(xi + xB)
......
......@@ -25,7 +25,7 @@ GIS_TEMPLATE = "data/UK2019mod_pop.gpkg"
# Reproduction number calculation
def calc_R_it(param, events, init_state, covar_data):
def calc_R_it(param, events, init_state, covar_data, priors):
"""Calculates effective reproduction number for batches of metapopulations
:param theta: a tensor of batched theta parameters [B] + theta.shape
:param xi: a tensor of batched xi parameters [B] + xi.shape
......@@ -36,21 +36,36 @@ def calc_R_it(param, events, init_state, covar_data):
"""
def r_fn(args):
beta1_, beta2_, xi_, gamma0_, events_ = args
beta1_, beta2_, beta_3, sigma_, xi_, gamma0_, events_ = args
t = events_.shape[-2] - 1
state = compute_state(init_state, events_, model_spec.STOICHIOMETRY)
state = tf.gather(state, t, axis=-2) # State on final inference day
model = model_spec.CovidUK(
covariates=covar_data,
initial_state=init_state,
initial_step=0,
num_steps=events_.shape[-2],
priors=priors,
)
xi_pred = model_spec.conditional_gp(
model.model["xi"](beta1_, sigma_),
xi_,
tf.constant(
[events.shape[-2] + model_spec.XI_FREQ], dtype=model_spec.DTYPE
)[:, tf.newaxis],
)
par = dict(
beta1=beta1_,
beta2=beta2_,
beta3=tf.constant(
[0.0, 0.0, 0.0], dtype=model_spec.DTYPE
), # xi_[1:4],
beta3=tf.concat([beta_3, [0.0]], axis=-1),
sigma=sigma_,
gamma0=gamma0_,
xi=xi_,
xi=xi_, # tf.reshape(xi_pred.sample(), [1]),
)
print("xi shape:", par["xi"].shape)
ngm_fn = model_spec.next_generation_matrix_fn(covar_data, par)
ngm = ngm_fn(t, state)
return ngm
......@@ -60,6 +75,8 @@ def calc_R_it(param, events, init_state, covar_data):
elems=(
param["beta1"],
param["beta2"],
param["beta3"],
param["sigma"],
param["xi"],
param["gamma0"],
events,
......@@ -82,12 +99,12 @@ def predicted_incidence(param, init_state, init_step, num_steps, priors):
"""
def sim_fn(args):
beta1_, beta2_, xi_, gamma0_, gamma1_, init_ = args
beta1_, beta2_, beta3_, sigma_, xi_, gamma0_, gamma1_, init_ = args
par = dict(
beta1=beta1_,
beta2=beta2_,
beta3=[0.0, 0.0, 0.0],
beta3=tf.concat([beta3_, [0.0]], axis=-1),
gamma0=gamma0_,
gamma1=gamma1_,
xi=xi_,
......@@ -108,6 +125,8 @@ def predicted_incidence(param, init_state, init_step, num_steps, priors):
elems=(
param["beta1"],
param["beta2"],
param["beta3"],
param["sigma"],
param["xi"],
param["gamma0"],
param["gamma1"],
......@@ -176,6 +195,12 @@ if __name__ == "__main__":
param = dict(
beta1=posterior["samples/beta1"][idx],
beta2=posterior["samples/beta2"][idx],
beta3=posterior["samples/beta3"][
idx,
],
sigma=posterior["samples/sigma"][
idx,
],
xi=posterior["samples/xi"][idx],
gamma0=posterior["samples/gamma0"][idx],
gamma1=posterior["samples/gamma1"][idx],
......@@ -195,7 +220,9 @@ if __name__ == "__main__":
priors=config["mcmc"]["prior"],
)
ngms = calc_R_it(param, events, init_state, covar_data)
ngms = calc_R_it(
param, events, init_state, covar_data, config["mcmc"]["prior"]
)
b, _ = power_iteration(ngms)
rt = rayleigh_quotient(ngms, b)
q = np.arange(0.05, 1.0, 0.05)
......
......@@ -3,7 +3,7 @@
data:
mobility_matrix: data/mergedflows.csv
population_size: data/c2019modagepop.csv
commute_volume: data/201113_OFF_SEN_COVID19_road_traffic_national_table.xlsx
commute_volume: data/201127_OFF_SEN_COVID19_road_traffic_national_table.xlsx
reported_cases: data/Anonymised Combined Line List 20201109.csv
case_date_type: specimen
pillar: both
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment