Commit 5f68ca50 authored by Chris Jewell's avatar Chris Jewell
Browse files

Implemented uncertainty on observations

1. Added function to read in negative test data
2. Implemented a sampling layer for testing data.
3. Enabled move and occult samplers on removal times.
parent 595ea2eb
......@@ -137,6 +137,35 @@ due to missing values ({100. * (orig_len - line_listing.shape[0])/orig_len}%)"
)
def read_phe_neg_cases(path, date_low, date_high, ltlas=None):
"""Reads negative test cases from PHE negative testing data
:param path: path to Negative testing data
:param date_low: lower date bound
:param date_high: upper date bound
:param returns: a Pandas data frame of LTLAs x dates
"""
dat = pd.read_csv(path, usecols=["LTLA", "earliestspecimendate"])
dat.columns = ["lad19cd", "date"]
dat["date"] = pd.to_datetime(dat["date"], format="%Y-%m-%d")
dat["lad19cd"] = _merge_ltla(dat["lad19cd"])
dat = dat[(date_low <= dat["date"]) & (dat["date"] < date_high)]
counts = dat.groupby(["date", "lad19cd"]).size()
counts.name = "count"
dates = pd.date_range(date_low, date_high, closed="left")
if ltlas is None:
ltlas = counts.index.levels[1]
index = pd.MultiIndex.from_product(
[dates, ltlas], names=["date", "lad19cd"]
)
counts = counts.reindex(index, fill_value=0)
return counts.reset_index().pivot(
index="lad19cd", columns="date", values="count"
)
def read_tier_restriction_data(
tier_restriction_csv, lad19cd_lookup, date_low, date_high
):
......
......@@ -17,7 +17,7 @@ from gemlib.mcmc import MultiScanKernel
from gemlib.mcmc import AdaptiveRandomWalkMetropolis
from gemlib.mcmc import Posterior
from covid.data import read_phe_cases
from covid.data import read_phe_cases, read_phe_neg_cases
from covid.cli_arg_parse import cli_args
import model_spec
......@@ -61,6 +61,14 @@ if __name__ == "__main__":
pillar=config["data"]["pillar"],
).astype(DTYPE)
neg_tests = read_phe_neg_cases(
config["data"]["negative_tests"],
date_low=inference_period[0],
date_high=inference_period[1],
).astype(DTYPE)
covar_data["total_tests"] = cases + neg_tests
# Impute censored events, return cases
events = model_spec.impute_censored_events(cases)
......@@ -81,6 +89,7 @@ if __name__ == "__main__":
start_time = state.shape[1] - cases.shape[1]
initial_state = state[:, start_time, :]
events = events[:, start_time:, :]
observed_events = events[..., 2]
num_metapop = covar_data["N"].shape[0]
########################################################
......@@ -114,6 +123,7 @@ if __name__ == "__main__":
beta1=block1[0],
xi=block1[1:],
seir=events,
y=observed_events,
)
)
......@@ -191,7 +201,8 @@ if __name__ == "__main__":
),
cumulative_event_offset=initial_state,
nmax=config["mcmc"]["occult_nmax"],
t_range=(events.shape[1] - 21, events.shape[1]),
# t_range=(events.shape[1] - 21, events.shape[1]),
t_range=(0, events.shape[1]),
name=name,
),
name=name,
......@@ -207,8 +218,10 @@ if __name__ == "__main__":
kernel_list=[
(0, make_partially_observed_step(0, None, 1, "se_events")),
(0, make_partially_observed_step(1, 0, 2, "ei_events")),
(0, make_partially_observed_step(1, 2, None, "ir_events")),
(0, make_occults_step(None, 0, 1, "se_occults")),
(0, make_occults_step(0, 1, 2, "ei_occults")),
(0, make_occults_step(1, 2, None, "ir_occults")),
],
name="gibbs1",
),
......@@ -248,8 +261,10 @@ if __name__ == "__main__":
res1 = res0[2].inner_results
results_dict["move/S->E"] = get_move_results(res1[0])
results_dict["move/E->I"] = get_move_results(res1[1])
results_dict["occult/S->E"] = get_move_results(res1[2])
results_dict["occult/E->I"] = get_move_results(res1[3])
results_dict["move/I->R"] = get_move_results(res1[2])
results_dict["occult/S->E"] = get_move_results(res1[3])
results_dict["occult/E->I"] = get_move_results(res1[4])
results_dict["occult/I->R"] = get_move_results(res1[5])
return results_dict
......
......@@ -112,7 +112,6 @@ def CovidUK(covariates, initial_state, initial_step, num_steps, priors):
)
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)
......@@ -123,12 +122,6 @@ def CovidUK(covariates, initial_state, initial_step, num_steps, priors):
)
def gamma0():
# return tfd.Gamma(
# concentration=tf.constant(
# priors["gamma"]["concentration"], dtype=DTYPE
# ),
# rate=tf.constant(priors["gamma"]["rate"], dtype=DTYPE),
# )
return tfd.Normal(
loc=tf.constant(0.0, dtype=DTYPE),
scale=tf.constant(100.0, dtype=DTYPE),
......@@ -209,6 +202,15 @@ def CovidUK(covariates, initial_state, initial_step, num_steps, priors):
num_steps=num_steps,
)
def y(seir):
"""Observation noise process"""
p_hat = seir[..., -1] / covariates["total_tests"]
return tfd.Independent(
tfd.Binomial(total_count=covariates["total_tests"], probs=p_hat),
reinterpreted_batch_ndims=2,
)
return tfd.JointDistributionNamed(
dict(
beta1=beta1,
......@@ -219,6 +221,7 @@ def CovidUK(covariates, initial_state, initial_step, num_steps, priors):
gamma0=gamma0,
gamma1=gamma1,
seir=seir,
y=y,
)
)
......
Markdown is supported
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