In dchodge/pilot.hcw: What the Package Does (One Line, Title Case)
library(latex2exp)
library(patchwork)
library(here)
require(plyr)
require(data.table)
require(tidyverse)
library(devtools)
library(loo)
# Required for this analysis
require(reshape2)
require(foreach)
require(doParallel)
require(bayesplot)
require(coda)
require(ggplot2)
require(viridis)
require(ggpubr)
library(mvoutlier)
require(foreach)
require(doParallel)
require(Rcpp)
require(extrafont)
require(showtext)
font_add("Avenir", "/Users/davidhodgson/Documents/personal/fonts/avenir_ff/AvenirLTStd-Book.otf")
showtext_auto()
theme_set(theme_light())
library(wesanderson)
Load packages and data
Currently using the serosolver
#install_github("https://github.com/dchodge/serosolver", ref = "custom_ab")
#install("/Users/davidhodgson/Documents/research/flu/serosolver_vac")
#install(here::here("rcppfunchcw"))
devtools::load_all()
require(serosolver)
require(rcppfunchcw)
load(here::here("data", "hanam_data.RDS")) # hcwpre
Plot sensible stuff for Lognormal model
# Log normal
load(file = here::here("data", paste0("modelinfo_cross_", hanam$study_name_short, ".RDS")))
load(file = here::here("outputs", "hanam", "fits", "dis_dep", "full", "LN_post_full.Rdata"))
all_post_test_full_raw <- post_cross$post_raw
all_post_test_full <- post_cross$post_trim
sero_data <- all_models_hanam[[2]]$create_post$titre_data
sero_data_base <- sero_data %>%
group_by(virus_isol_yr, sample_time) %>%
summarise(titre_mean = mean(titre_value, rm.na = TRUE)) %>%
pivot_wider(names_from = sample_time, values_from = titre_mean) %>%
mutate(day4 = (`4` - `0`), day7 = (`7` - `0`), day14 = (`14` - `0`), day21 = (`21` - `0`), day280 = (`280` - `0`) ) %>%
select(virus_isol_yr, base = `0`, day4, day7, day14, day21, day280) %>%
pivot_longer(c(day4, day7, day14, day21, day280), names_to = "time", values_to = "abs") %>%
mutate(time = case_when(time == "day4"~"4", time == "day7"~"7", time == "day14"~"14", time == "day21"~"21", time == "day280"~"280"))
parsfunctionalLN <- all_post_test_full %>% bind_rows %>% select(mu_LN, sigma_LN, alpha_LN)
mu_LN_v <- parsfunctionalLN$mu_LN
sigma_LN_v <- parsfunctionalLN$sigma_LN
alpha_LN_v <- parsfunctionalLN$alpha_LN
functionalLN <- function(x, mu, sigma, alpha) alpha * dlnorm(x, mu, sigma)
dataplot <- sero_data_base %>% filter(virus_isol_yr == 2016) %>% mutate(time = as.numeric(time) / 7)
ggplot() +
xlim(0, 52) + ylim(0, 4) +
geom_point(data = dataplot, aes(x = time, y = abs), size = 3, alpha = 0.8) +
lapply(seq(1000, 2000, 10),
function(x) geom_function(fun = functionalLN,
args = list(mu = mu_LN_v[x], sigma = sigma_LN_v[x], alpha = alpha_LN_v[x]),
color = "red", alpha = 0.5, size = 0.1)
) +
theme_minimal() + labs(x = "Week since vaccination", y = "Titre increase")
# Log normal
load(file = here::here("data", paste0("modelinfo_cross_", hanam$study_name_short, ".RDS")))
load(file = here::here("outputs", "hanam", "fits", "dis_dep", "full", "LN_post_full.Rdata"))
all_post_test_full_raw <- post_cross$post_raw
all_post_test_full <- post_cross$post_trim
model_no <- 2
chain_short <- as.data.frame(all_post_test_full_raw$theta_chain) %>% filter(chain_no == 1)
infection_histories <- all_post_test_full_raw$inf_chain %>% filter(chain_no == 1)
titre_dat <- all_models_hanam[[model_no]]$create_post$titre_data
vaccination_histories <- all_models_hanam[[model_no]]$create_post$vaccination_histories
antigenic_map <- all_models_hanam[[model_no]]$create_post$antigenic_map
strain_isolation_times <- all_models_hanam[[model_no]]$create_post$strain_isolation_times
par_tab <- all_models_hanam[[model_no]]$create_post$par_tab
vaccination_histories <- all_models_hanam[[model_no]]$create_post$vac_history
custom_ab_kin_func <- all_models_hanam[[model_no]]$custom_ab_kin_func
custom_ab_kin_func_1 <- all_models_hanam[[model_no]]$custom_ab_kin_func
custom_antigenic_maps_func <- all_models_hanam[[model_no]]$custom_antigenic_maps_func
expand_titredat <- FALSE
nsamp <- 100
titre_before_infection <- FALSE
measurement_indices_by_time <- NULL
individuals <- unique(titre_dat$individual)
mu_indices <- NULL
plot_infection_histories_long(
chain = chain_short,
infection_histories = infection_histories,
vaccination_histories = vaccination_histories,
titre_dat = titre_dat,
individuals = 1,
antigenic_map = antigenic_map,
strain_isolation_times = strain_isolation_times,
par_tab = par_tab,
custom_ab_kin_func = custom_ab_kin_func,
custom_antigenic_maps_func = custom_antigenic_maps_func,
size_plot = 25)
Plot sensible stuff for SP-model
# Set Point
load(file = here::here("data", paste0("modelinfo_cross_", hanam$study_name_short, ".RDS")))
load(file = here::here("outputs", "hanam", "fits", "dis_dep", "full", "SP_post_full.Rdata"))
all_post_test_full_raw <- post_cross$post_raw
all_post_test_full <- post_cross$post_trim
sero_data <- all_models_hanam[[3]]$create_post$titre_data
sero_data_base <- sero_data %>%
group_by(virus_isol_yr, sample_time) %>%
summarise(titre_mean = mean(titre_value, na.rm = TRUE)) %>%
pivot_wider(names_from = sample_time, values_from = titre_mean) %>%
mutate(day4 = (`4` - `0`), day7 = (`7` - `0`), day14 = (`14` - `0`), day21 = (`21` - `0`), day280 = (`280` - `0`) ) %>%
select(virus_isol_yr, base = `0`, day4, day7, day14, day21, day280) %>%
pivot_longer(c(day4, day7, day14, day21, day280), names_to = "time", values_to = "abs") %>%
mutate(time = case_when(time == "day4"~"4", time == "day7"~"7", time == "day14"~"14", time == "day21"~"21", time == "day280"~"280"))
dataplot <- sero_data_base %>% filter(virus_isol_yr == 2016) %>% mutate(time = as.numeric(time) / 7)
dataplot_2 <- sero_data %>% filter(virus_isol_yr == 2016) %>% select(pid, virus_isol_yr, sample_time, titre) %>%
pivot_wider(names_from = sample_time, values_from = titre) %>%
mutate(day4 = (`4` - `0`), day7 = (`7` - `0`), day14 = (`14` - `0`), day21 = (`21` - `0`), day280 = (`280` - `0`) ) %>%
select(virus_isol_yr, base = `0`, day4, day7, day14, day21, day280) %>%
pivot_longer(c(day4, day7, day14, day21, day280), names_to = "time", values_to = "abs") %>%
mutate(time = case_when(time == "day4"~"4", time == "day7"~"7", time == "day14"~"14", time == "day21"~"21", time == "day280"~"280")) %>%
mutate(time = as.numeric(time) / 7) %>%
group_by(time, abs)
dataplot_2_full <- dataplot_2 %>% summarise(n = n())
dataplot_2_sum <- dataplot_2 %>% ungroup %>% as.data.frame %>% group_by(time) %>% summarise(abs_mean = mean(abs, na.rm = TRUE))
parsfunctionalSP <- all_post_test_full[[1]][[2]] %>% select(th1_SP, th2_SP, th3_SP, alpha_SP)
th1_SP_v <- parsfunctionalSP$th1_SP
th2_SP_v <- parsfunctionalSP$th2_SP
th3_SP_v <- parsfunctionalSP$th3_S
alpha_SP_v <- parsfunctionalSP$alpha_SP
boost_to_set_point <- function(x, th1, th2, th3, alpha) {
alpha * (th3 / (th2 + th3) * (1- exp(-(th2 + th3)* x)) - (th1 - th3) / (th1 - th2 - th3) * (exp(-th1 * x) - exp(- (th2 + th3) * x)))
}
ggplot() +
xlim(0, 52) + ylim(0, 4) +
geom_point(data = dataplot_2_full, aes(x = time, y = abs, size = n), shape = 21, fill = "gray", alpha = 0.5) +
geom_point(data = dataplot_2_sum, aes(x = time, y = abs_mean), shape = 22, size = 5, fill = "black", alpha = 0.5) +
lapply(seq(500, 1000, 5),
function(x) geom_function(fun = boost_to_set_point,
args = list(th1 = th1_SP_v[x], th2 = th2_SP_v[x], th3 = th3_SP_v[x], alpha = alpha_SP_v[x]),
color = "red", alpha = 0.5, size = 0.1)
) +
lapply(seq(500, 1000, 5),
function(x) geom_function(fun = functionalLN,
args = list(mu = mu_LN_v[x], sigma = sigma_LN_v[x], alpha = alpha_LN_v[x]),
color = "blue", alpha = 0.5, size = 0.1)
) +
theme_bw() + labs(x = "Week since vaccination", y = "Titre increase")
ggsave(file = here::here("outputs", "hanam", "analysis", "figs", "figD.pdf"))
Look at cross immunity
melt_antigenic_coords <- function(anti.map.in) { # anti.map.in can be vector or matrix - rows give inf_times, columns give location
# Calculate antigenic distances
if (is.null(dim(anti.map.in))) { # check if input map is one or 2 dimensions
# If 1D antigenic 'line' defined, calculate distances directory from input
(dmatrix <- sapply(anti.map.in, function(x) {
y <- abs(anti.map.in - x)
y
}))
} else { # If 2D antigenic map defined, calculate distances directory from input
(dmatrix <- apply(
anti.map.in, 1,
function(x) {
y <- sqrt(colSums(apply(anti.map.in, 1, function(y) {
(y - x)^2
})))
y
}
))
}
}
# Log normal
load(file = here::here("data", paste0("modelinfo_cross_", hanam$study_name_short, ".RDS")))
load(file = here::here("outputs", "hanam", "fits", "dis_dep", "full", "LN_post_full.Rdata"))
all_post_test_full_raw_ln <- post_cross$post_raw
all_post_test_full_ln <- post_cross$post_trim
load(file = here::here("outputs", "hanam", "fits", "dis_dep", "full", "SP_post_full.Rdata"))
all_post_test_full_raw_sp <- post_cross$post_raw
all_post_test_full_sp <- post_cross$post_trim
antigenic_map <- all_models_hanam[[1]]$create_post$antigenic_map
ln_cross <- all_post_test_full_ln[[1]][[2]] %>% select(mu, sigma1, sigma1_vac)
sp_cross <- all_post_test_full_sp[[1]][[2]] %>% select(mu, sigma1, sigma1_vac)
get_rel_vec <- function(sigma) {
antigenic_map_melted <- c(melt_antigenic_coords(antigenic_map[, c("x_coord", "y_coord")]))
temp <- 1 - antigenic_map_melted * sigma
temp[temp < 0] = 0
tempmat <- temp %>% matrix(nrow = 49)
tempmat[, 49]
}
ln_mu <- ln_cross[, 1]
sp_mu <- sp_cross[, 1]
bind_rows(
data.frame(model = "log normal", mu = ln_mu),
data.frame(model = "set point", mu = sp_mu)
) %>%
ggplot() +
geom_boxplot(aes(x = model, y = mu)) +
labs(x = "Model", y = "Long-term boost due to natural infection")
ln_n_s <- mean(ln_cross[, 2])
ln_v_s <- mean(ln_cross[, 2] * ln_cross[, 3])
sp_n_s <- mean(sp_cross[, 2])
sp_v_s <- mean(sp_cross[, 2] * sp_cross[, 3])
df_sigma_plot_natural <- bind_rows(
data.frame(type = "log normal", time = 1968:2016, cross = get_rel_vec(ln_n_s)),
data.frame(type = "set point", time = 1968:2016, cross = get_rel_vec(sp_n_s))
)
df_sigma_plot_vaccination <- bind_rows(
data.frame(type = "log normal", time = 1968:2016, cross = get_rel_vec(ln_v_s)),
data.frame(type = "set point", time = 1968:2016, cross = get_rel_vec(sp_v_s))
)
p1 <- df_sigma_plot_natural %>%
ggplot() +
geom_bar(aes(x = time, y = cross, fill = type), stat="identity",position = "identity") +
facet_grid(rows = vars(type)) +
labs(x = "Year", y = "Relative boosting", fill = "Model", title = "Long-term cross immunity due to natural infection")
p2 <- df_sigma_plot_vaccination %>%
ggplot() +
geom_bar(aes(x = time, y = cross, fill = type), stat="identity",position = "identity") +
facet_grid(rows = vars(type)) +
labs(x = "Year", y = "Relative boosting", fill = "Model", title = "Cross immunity due to vaccination infection")
require(patchwork)
p1 / p2 + plot_layout(guides = "collect")
ln_post <- all_post_test_full_ln[[1]][[2]]
ln_n_s <- mean(ln_post$sigma1)
ln_v_s <- mean(ln_post$sigma1 * ln_post$sigma1_vac)
functionalLN <- function(x, mu, sigma, alpha) alpha * dlnorm(x, mu, sigma)
mu_B_LN <- ln_post$mu[100]
mu_LN <- ln_post$mu_LN[100]
sigma_LN <- ln_post$sigma_LN[100]
alpha_LN <- ln_post$alpha_LN[100]
sp_post <- all_post_test_full_sp[[1]][[2]]
sp_n_s <- mean(sp_post$sigma1)
sp_v_s <- mean(sp_post$sigma1 * sp_post$sigma1_vac)
boost_to_set_point <- function(x, th1, th2, th3, alpha) {
alpha * (th3 / (th2 + th3) * (1- exp(-(th2 + th3)* x)) - (th1 - th3) / (th1 - th2 - th3) * (exp(-th1 * x) - exp(- (th2 + th3) * x)))
}
mu_B_SP <- sp_post$mu[100]
th1_SP <- sp_post$th1_SP[100]
th2_SP <- sp_post$th2_SP[100]
th3_SP <- sp_post$th3_SP[100]
alpha_SP <- sp_post$alpha_SP[100]
data_plot_all <- list()
for(i in 1:52){
data_plot_all[[i]] <- data.frame(
x = i,
y = 1968:2016,
lognormal = get_rel_vec(ln_v_s) * functionalLN(i, mu_LN, sigma_LN, alpha_LN),
setpoint = get_rel_vec(sp_v_s) * boost_to_set_point(i, th1_SP, th2_SP, th3_SP, alpha_SP)
)
}
df_data_plot_all <- data_plot_all %>% bind_rows
df_data_plot_all_long <- df_data_plot_all %>% pivot_longer(!c(x,y), names_to = "model_type", values_to = "titre_boost")
df_data_plot_all_long %>%
ggplot() +
geom_tile(aes(x = x, y = y, fill = titre_boost), color = "black") +
scale_fill_gradient(low = "white",
high = "red") +
labs(x = "Weeks post vaccination", y = "Virus year", fill = "Titre boost") +
facet_grid(vars(model_type))
ggsave(file = here::here("outputs", "hanam", "analysis", "figs", "figE.pdf"), height = 5, width = 8)
i_vals <- seq(1, 52 * 10, 52)
data_plot_all_nat <- list()
for(i in 1:10){
data_plot_all_nat[[i]] <- data.frame(
x = i,
y = 1968:2016,
lognormal_nat = get_rel_vec(ln_n_s) * mu_B_LN,
sp_nat = get_rel_vec(sp_n_s) * mu_B_SP,
lognormal_vac = get_rel_vec(ln_v_s) * functionalLN(i_vals[i], mu_LN, sigma_LN, alpha_LN),
sp_vac = get_rel_vec(sp_v_s) * boost_to_set_point(i_vals[i], th1_SP, th2_SP, th3_SP, alpha_SP)
)
}
df_data_plot_all_nat <- data_plot_all_nat %>% bind_rows
df_data_plot_all_nat_long <- df_data_plot_all_nat %>% pivot_longer(!c(x,y), names_to = "model_type", values_to = "titre_boost")
df_data_plot_all_nat_long %>%
ggplot() +
geom_tile(aes(x = x, y = y, fill = titre_boost), color = "black") +
scale_fill_gradient(low = "white",
high = "blue") +
labs(x = "Years post vaccination", y = "Virus year", fill = "Titre boost") +
facet_grid(vars(model_type))
ggsave(file = here::here("outputs", "hanam", "analysis", "figs", "figF.pdf"), height = 5, width = 8)
model_no <- 3
chain_short <- as.data.frame(all_post_test_full_raw$theta_chain) %>% filter(chain_no == 1)
infection_histories <- all_post_test_full_raw$inf_chain %>% filter(chain_no == 1)
titre_dat <- all_models_hanam[[model_no]]$create_post$titre_data
vaccination_histories <- all_models_hanam[[model_no]]$create_post$vaccination_histories
antigenic_map <- all_models_hanam[[model_no]]$create_post$antigenic_map
strain_isolation_times <- all_models_hanam[[model_no]]$create_post$strain_isolation_times
par_tab <- all_models_hanam[[model_no]]$create_post$par_tab
vaccination_histories <- all_models_hanam[[model_no]]$create_post$vac_history
custom_ab_kin_func <- all_models_hanam[[model_no]]$custom_ab_kin_func
custom_ab_kin_func_1 <- all_models_hanam[[model_no]]$custom_ab_kin_func
custom_antigenic_maps_func <- all_models_hanam[[model_no]]$custom_antigenic_maps_func
expand_titredat <- FALSE
nsamp <- 100
titre_before_infection <- FALSE
measurement_indices_by_time <- NULL
individuals <- unique(titre_dat$individual)
mu_indices <- NULL
plot_infection_histories_long(
chain = chain_short,
infection_histories = infection_histories,
vaccination_histories = vaccination_histories,
titre_dat = titre_dat,
individuals = 1,
antigenic_map = antigenic_map,
strain_isolation_times = strain_isolation_times,
par_tab = par_tab,
custom_ab_kin_func = custom_ab_kin_func,
custom_antigenic_maps_func = custom_antigenic_maps_func,
size_plot = 25)
# Log normal
load(file = here::here("data", paste0("modelinfo_cross_", hanam$study_name_short, ".RDS")))
load(file = here::here("outputs", "hanam", "fits", "dis_dep", "full", "LN_post_full.Rdata"))
all_post_test_full_raw_ln <- post_cross$post_raw
all_post_test_full_ln <- post_cross$post_trim
load(file = here::here("outputs", "hanam", "fits", "dis_dep", "full", "SP_post_full.Rdata"))
all_post_test_full_raw_sp <- post_cross$post_raw
all_post_test_full_sp <- post_cross$post_trim
inf_chain <- all_post_test_full_raw_sp$inf_chain
sample_nos <- inf_chain %>% pull(sampno) %>% unique %>% length
chain <- inf_chain %>% pull(chain_no) %>% max
sample_df <- inf_chain %>%
filter(sampno > 100000) %>%
group_by(year = j, sampno) %>%
dplyr::summarise(x = sum(x / (chain * 101))) %>% as.data.frame %>%
filter(year < 49)
sample_df %>%
ggplot() +
geom_boxplot(aes(x = as.character(year + 1967), y = x)) +
scale_x_discrete(guide = guide_axis(check.overlap = TRUE)) +
labs(x = "Year", y = "Attack rate", title = "Estimated attack rate")
load(file = here::here("data", paste0("modelinfo_full_", hcwpre$study_name_short, ".RDS"))) #all_models_hcw_pre_cross
load(file = here::here("outputs", "hcw_pre", "fits", "dis_dep", "full", "trimmed_post_full.Rdata"))
pos <- 5
all_post_test_full <- post_full$post_raw
chain_full <- as.data.frame(all_post_test_full[[pos]]$theta_chain) %>% filter(chain_no == 1)
infection_histories <- all_post_test_full[[pos]]$inf_chain %>% filter(chain_no == 1)
titre_dat <- all_models_hcw_pre_full[[pos]]$create_post$titre_data
vaccination_histories <- all_models_hcw_pre_full[[pos]]$create_post$vaccination_histories
antigenic_map <- all_models_hcw_pre_full[[pos]]$create_post$antigenic_map
strain_isolation_times <- all_models_hcw_pre_full[[pos]]$create_post$strain_isolation_times
par_tab <- all_models_hcw_pre_full[[pos]]$create_post$par_tab
vaccination_histories <- all_models_hcw_pre_full[[pos]]$create_post$vac_history
custom_ab_kin_func <- all_models_hcw_pre_full[[pos]]$custom_ab_kin_func
custom_ab_kin_func_1 <- all_models_hcw_pre_full[[pos]]$custom_ab_kin_func
custom_antigenic_maps_func <- all_models_hcw_pre_full[[pos]]$custom_antigenic_maps_func
expand_titredat <- FALSE
nsamp <- 100
titre_before_infection <- FALSE
measurement_indices_by_time <- NULL
individuals <- unique(titre_dat$individual)
mu_indices <- NULL
chain_full <- chain_full %>% filter(sampno > 300000)
y <- get_titre_predictions(
chain = chain_full,
infection_histories = infection_histories,
vaccination_histories = vaccination_histories,
titre_dat = titre_dat,
individuals = 1:49,
antigenic_map = antigenic_map,
strain_isolation_times = strain_isolation_times,
par_tab = par_tab,
custom_ab_kin_func = custom_ab_kin_func,
custom_antigenic_maps_func = custom_antigenic_maps_func,
nsamp = 100
)
plot_infection_histories_long(
chain = chain_full,
infection_histories = infection_histories,
vaccination_histories = vaccination_histories,
titre_dat = titre_dat,
individuals = 21,
antigenic_map = antigenic_map,
strain_isolation_times = strain_isolation_times,
par_tab = par_tab,
custom_ab_kin_func = custom_ab_kin_func,
custom_antigenic_maps_func = custom_antigenic_maps_func
)
chain_full$mu_vac %>% mean
chain_full$mu_short_vac %>% mean
sum_data <- titre_dat %>% group_by(sample_time, titre) %>% summarise(n = n())
weighted_sum_data <- sum_data %>% summarise(w_mean = weighted.mean(titre, n, na.rm = TRUE))
#predictions
titre_pred <- y$predicted_observations
titre_pred <- titre_pred %>% filter(virus == 24192)
titre_pred %>% select(virus, samples, titre, median)
sum_data <- titre_pred %>% group_by(sample_time, titre) %>% summarise(n = n())
weighted_sum_data <- sum_data %>% summarise(w_mean = weighted.mean(titre, n, na.rm = TRUE))
sum_pred <- titre_pred %>% group_by(sample_time, median) %>% summarise(n = n())
weighted_sum_pred <- sum_pred %>% summarise(w_mean = weighted.mean(median, n, na.rm = TRUE))
Post plots
all_post_test_full_trim <- post_full$post_trim
postplots_long <- all_post_test_full_trim[[5]] %>% bind_rows
dec <- 1:12 %>% map(~
(postplots_long[["mu_short_vac"]] * (1.0 - ((1 / 12 + postplots_long[["wane_vac"]]) * .x )))
)
df_dec <- dec %>% setNames(as.character(c(1:12))) %>%
as.data.frame %>% pivot_longer(c(1:12), values_to = "titre", names_to = "time") %>%
mutate(time = str_extract(time, "[0-9,.]+")) %>%
mutate(time = factor(time, levels = as.character(c(1:12))))
df_dec <- df_dec[df_dec$titre > 0, ]
sero_data <- hcwpre$sero_data
sum_data <- sero_data %>% group_by(sample_time, titre_value) %>% summarise(n = n())
weighted_sum_data <- sum_data %>% summarise(w_mean = weighted.mean(titre_value, n, na.rm = TRUE)) %>%
mutate(sample_time = sample_time / 30)
ggplot() +
geom_point(data = sum_data, aes(x = sample_time, y = titre_value, size = n), fill = "black", alpha = 0.5) +
geom_line(data = weighted_sum_data, aes(x = sample_time, y = w_mean), color = "red", size = 1) +
geom_point(data = weighted_sum_data, aes(x = sample_time, y = w_mean), shape = 21, fill = "red", size = 5) +
labs(x = "Time post vaccination", y = "Titre value", size = "Number of samples")
df_dec %>%
ggplot() +
geom_boxplot(aes(x = time, y = titre)) +
labs(x = "Months post vaccination") +
geom_point(data = weighted_sum_data, aes(x = sample_time, y = w_mean), shape = 8, color = "red", size = 5) +
geom_hline(yintercept = postplots_long[["mu"]] * postplots_long[["mu_vac"]] %>% mean)
antigenic_map_melted <- c(melt_antigenic_coords(antigenic_map[, c("x_coord", "y_coord")]))
chain$sigma1 %>% mean
s <- 0.01
lapply(antigenic_map_melted, function(x) max(1 - x * s, 0)) %>% unlist
dchodge/pilot.hcw documentation built on Feb. 2, 2022, 11:29 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
library(latex2exp) library(patchwork) library(here) require(plyr) require(data.table) require(tidyverse) library(devtools) library(loo) # Required for this analysis require(reshape2) require(foreach) require(doParallel) require(bayesplot) require(coda) require(ggplot2) require(viridis) require(ggpubr) library(mvoutlier) require(foreach) require(doParallel) require(Rcpp) require(extrafont) require(showtext) font_add("Avenir", "/Users/davidhodgson/Documents/personal/fonts/avenir_ff/AvenirLTStd-Book.otf") showtext_auto() theme_set(theme_light()) library(wesanderson)
Load packages and data
Currently using the serosolver
#install_github("https://github.com/dchodge/serosolver", ref = "custom_ab") #install("/Users/davidhodgson/Documents/research/flu/serosolver_vac") #install(here::here("rcppfunchcw")) devtools::load_all() require(serosolver) require(rcppfunchcw) load(here::here("data", "hanam_data.RDS")) # hcwpre
Plot sensible stuff for Lognormal model
# Log normal load(file = here::here("data", paste0("modelinfo_cross_", hanam$study_name_short, ".RDS"))) load(file = here::here("outputs", "hanam", "fits", "dis_dep", "full", "LN_post_full.Rdata")) all_post_test_full_raw <- post_cross$post_raw all_post_test_full <- post_cross$post_trim sero_data <- all_models_hanam[[2]]$create_post$titre_data sero_data_base <- sero_data %>% group_by(virus_isol_yr, sample_time) %>% summarise(titre_mean = mean(titre_value, rm.na = TRUE)) %>% pivot_wider(names_from = sample_time, values_from = titre_mean) %>% mutate(day4 = (`4` - `0`), day7 = (`7` - `0`), day14 = (`14` - `0`), day21 = (`21` - `0`), day280 = (`280` - `0`) ) %>% select(virus_isol_yr, base = `0`, day4, day7, day14, day21, day280) %>% pivot_longer(c(day4, day7, day14, day21, day280), names_to = "time", values_to = "abs") %>% mutate(time = case_when(time == "day4"~"4", time == "day7"~"7", time == "day14"~"14", time == "day21"~"21", time == "day280"~"280")) parsfunctionalLN <- all_post_test_full %>% bind_rows %>% select(mu_LN, sigma_LN, alpha_LN) mu_LN_v <- parsfunctionalLN$mu_LN sigma_LN_v <- parsfunctionalLN$sigma_LN alpha_LN_v <- parsfunctionalLN$alpha_LN functionalLN <- function(x, mu, sigma, alpha) alpha * dlnorm(x, mu, sigma) dataplot <- sero_data_base %>% filter(virus_isol_yr == 2016) %>% mutate(time = as.numeric(time) / 7) ggplot() + xlim(0, 52) + ylim(0, 4) + geom_point(data = dataplot, aes(x = time, y = abs), size = 3, alpha = 0.8) + lapply(seq(1000, 2000, 10), function(x) geom_function(fun = functionalLN, args = list(mu = mu_LN_v[x], sigma = sigma_LN_v[x], alpha = alpha_LN_v[x]), color = "red", alpha = 0.5, size = 0.1) ) + theme_minimal() + labs(x = "Week since vaccination", y = "Titre increase") # Log normal load(file = here::here("data", paste0("modelinfo_cross_", hanam$study_name_short, ".RDS"))) load(file = here::here("outputs", "hanam", "fits", "dis_dep", "full", "LN_post_full.Rdata")) all_post_test_full_raw <- post_cross$post_raw all_post_test_full <- post_cross$post_trim model_no <- 2 chain_short <- as.data.frame(all_post_test_full_raw$theta_chain) %>% filter(chain_no == 1) infection_histories <- all_post_test_full_raw$inf_chain %>% filter(chain_no == 1) titre_dat <- all_models_hanam[[model_no]]$create_post$titre_data vaccination_histories <- all_models_hanam[[model_no]]$create_post$vaccination_histories antigenic_map <- all_models_hanam[[model_no]]$create_post$antigenic_map strain_isolation_times <- all_models_hanam[[model_no]]$create_post$strain_isolation_times par_tab <- all_models_hanam[[model_no]]$create_post$par_tab vaccination_histories <- all_models_hanam[[model_no]]$create_post$vac_history custom_ab_kin_func <- all_models_hanam[[model_no]]$custom_ab_kin_func custom_ab_kin_func_1 <- all_models_hanam[[model_no]]$custom_ab_kin_func custom_antigenic_maps_func <- all_models_hanam[[model_no]]$custom_antigenic_maps_func expand_titredat <- FALSE nsamp <- 100 titre_before_infection <- FALSE measurement_indices_by_time <- NULL individuals <- unique(titre_dat$individual) mu_indices <- NULL plot_infection_histories_long( chain = chain_short, infection_histories = infection_histories, vaccination_histories = vaccination_histories, titre_dat = titre_dat, individuals = 1, antigenic_map = antigenic_map, strain_isolation_times = strain_isolation_times, par_tab = par_tab, custom_ab_kin_func = custom_ab_kin_func, custom_antigenic_maps_func = custom_antigenic_maps_func, size_plot = 25)
Plot sensible stuff for SP-model
# Set Point load(file = here::here("data", paste0("modelinfo_cross_", hanam$study_name_short, ".RDS"))) load(file = here::here("outputs", "hanam", "fits", "dis_dep", "full", "SP_post_full.Rdata")) all_post_test_full_raw <- post_cross$post_raw all_post_test_full <- post_cross$post_trim sero_data <- all_models_hanam[[3]]$create_post$titre_data sero_data_base <- sero_data %>% group_by(virus_isol_yr, sample_time) %>% summarise(titre_mean = mean(titre_value, na.rm = TRUE)) %>% pivot_wider(names_from = sample_time, values_from = titre_mean) %>% mutate(day4 = (`4` - `0`), day7 = (`7` - `0`), day14 = (`14` - `0`), day21 = (`21` - `0`), day280 = (`280` - `0`) ) %>% select(virus_isol_yr, base = `0`, day4, day7, day14, day21, day280) %>% pivot_longer(c(day4, day7, day14, day21, day280), names_to = "time", values_to = "abs") %>% mutate(time = case_when(time == "day4"~"4", time == "day7"~"7", time == "day14"~"14", time == "day21"~"21", time == "day280"~"280")) dataplot <- sero_data_base %>% filter(virus_isol_yr == 2016) %>% mutate(time = as.numeric(time) / 7) dataplot_2 <- sero_data %>% filter(virus_isol_yr == 2016) %>% select(pid, virus_isol_yr, sample_time, titre) %>% pivot_wider(names_from = sample_time, values_from = titre) %>% mutate(day4 = (`4` - `0`), day7 = (`7` - `0`), day14 = (`14` - `0`), day21 = (`21` - `0`), day280 = (`280` - `0`) ) %>% select(virus_isol_yr, base = `0`, day4, day7, day14, day21, day280) %>% pivot_longer(c(day4, day7, day14, day21, day280), names_to = "time", values_to = "abs") %>% mutate(time = case_when(time == "day4"~"4", time == "day7"~"7", time == "day14"~"14", time == "day21"~"21", time == "day280"~"280")) %>% mutate(time = as.numeric(time) / 7) %>% group_by(time, abs) dataplot_2_full <- dataplot_2 %>% summarise(n = n()) dataplot_2_sum <- dataplot_2 %>% ungroup %>% as.data.frame %>% group_by(time) %>% summarise(abs_mean = mean(abs, na.rm = TRUE)) parsfunctionalSP <- all_post_test_full[[1]][[2]] %>% select(th1_SP, th2_SP, th3_SP, alpha_SP) th1_SP_v <- parsfunctionalSP$th1_SP th2_SP_v <- parsfunctionalSP$th2_SP th3_SP_v <- parsfunctionalSP$th3_S alpha_SP_v <- parsfunctionalSP$alpha_SP boost_to_set_point <- function(x, th1, th2, th3, alpha) { alpha * (th3 / (th2 + th3) * (1- exp(-(th2 + th3)* x)) - (th1 - th3) / (th1 - th2 - th3) * (exp(-th1 * x) - exp(- (th2 + th3) * x))) } ggplot() + xlim(0, 52) + ylim(0, 4) + geom_point(data = dataplot_2_full, aes(x = time, y = abs, size = n), shape = 21, fill = "gray", alpha = 0.5) + geom_point(data = dataplot_2_sum, aes(x = time, y = abs_mean), shape = 22, size = 5, fill = "black", alpha = 0.5) + lapply(seq(500, 1000, 5), function(x) geom_function(fun = boost_to_set_point, args = list(th1 = th1_SP_v[x], th2 = th2_SP_v[x], th3 = th3_SP_v[x], alpha = alpha_SP_v[x]), color = "red", alpha = 0.5, size = 0.1) ) + lapply(seq(500, 1000, 5), function(x) geom_function(fun = functionalLN, args = list(mu = mu_LN_v[x], sigma = sigma_LN_v[x], alpha = alpha_LN_v[x]), color = "blue", alpha = 0.5, size = 0.1) ) + theme_bw() + labs(x = "Week since vaccination", y = "Titre increase") ggsave(file = here::here("outputs", "hanam", "analysis", "figs", "figD.pdf"))
Look at cross immunity
melt_antigenic_coords <- function(anti.map.in) { # anti.map.in can be vector or matrix - rows give inf_times, columns give location # Calculate antigenic distances if (is.null(dim(anti.map.in))) { # check if input map is one or 2 dimensions # If 1D antigenic 'line' defined, calculate distances directory from input (dmatrix <- sapply(anti.map.in, function(x) { y <- abs(anti.map.in - x) y })) } else { # If 2D antigenic map defined, calculate distances directory from input (dmatrix <- apply( anti.map.in, 1, function(x) { y <- sqrt(colSums(apply(anti.map.in, 1, function(y) { (y - x)^2 }))) y } )) } } # Log normal load(file = here::here("data", paste0("modelinfo_cross_", hanam$study_name_short, ".RDS"))) load(file = here::here("outputs", "hanam", "fits", "dis_dep", "full", "LN_post_full.Rdata")) all_post_test_full_raw_ln <- post_cross$post_raw all_post_test_full_ln <- post_cross$post_trim load(file = here::here("outputs", "hanam", "fits", "dis_dep", "full", "SP_post_full.Rdata")) all_post_test_full_raw_sp <- post_cross$post_raw all_post_test_full_sp <- post_cross$post_trim antigenic_map <- all_models_hanam[[1]]$create_post$antigenic_map ln_cross <- all_post_test_full_ln[[1]][[2]] %>% select(mu, sigma1, sigma1_vac) sp_cross <- all_post_test_full_sp[[1]][[2]] %>% select(mu, sigma1, sigma1_vac) get_rel_vec <- function(sigma) { antigenic_map_melted <- c(melt_antigenic_coords(antigenic_map[, c("x_coord", "y_coord")])) temp <- 1 - antigenic_map_melted * sigma temp[temp < 0] = 0 tempmat <- temp %>% matrix(nrow = 49) tempmat[, 49] } ln_mu <- ln_cross[, 1] sp_mu <- sp_cross[, 1] bind_rows( data.frame(model = "log normal", mu = ln_mu), data.frame(model = "set point", mu = sp_mu) ) %>% ggplot() + geom_boxplot(aes(x = model, y = mu)) + labs(x = "Model", y = "Long-term boost due to natural infection") ln_n_s <- mean(ln_cross[, 2]) ln_v_s <- mean(ln_cross[, 2] * ln_cross[, 3]) sp_n_s <- mean(sp_cross[, 2]) sp_v_s <- mean(sp_cross[, 2] * sp_cross[, 3]) df_sigma_plot_natural <- bind_rows( data.frame(type = "log normal", time = 1968:2016, cross = get_rel_vec(ln_n_s)), data.frame(type = "set point", time = 1968:2016, cross = get_rel_vec(sp_n_s)) ) df_sigma_plot_vaccination <- bind_rows( data.frame(type = "log normal", time = 1968:2016, cross = get_rel_vec(ln_v_s)), data.frame(type = "set point", time = 1968:2016, cross = get_rel_vec(sp_v_s)) ) p1 <- df_sigma_plot_natural %>% ggplot() + geom_bar(aes(x = time, y = cross, fill = type), stat="identity",position = "identity") + facet_grid(rows = vars(type)) + labs(x = "Year", y = "Relative boosting", fill = "Model", title = "Long-term cross immunity due to natural infection") p2 <- df_sigma_plot_vaccination %>% ggplot() + geom_bar(aes(x = time, y = cross, fill = type), stat="identity",position = "identity") + facet_grid(rows = vars(type)) + labs(x = "Year", y = "Relative boosting", fill = "Model", title = "Cross immunity due to vaccination infection") require(patchwork) p1 / p2 + plot_layout(guides = "collect") ln_post <- all_post_test_full_ln[[1]][[2]] ln_n_s <- mean(ln_post$sigma1) ln_v_s <- mean(ln_post$sigma1 * ln_post$sigma1_vac) functionalLN <- function(x, mu, sigma, alpha) alpha * dlnorm(x, mu, sigma) mu_B_LN <- ln_post$mu[100] mu_LN <- ln_post$mu_LN[100] sigma_LN <- ln_post$sigma_LN[100] alpha_LN <- ln_post$alpha_LN[100] sp_post <- all_post_test_full_sp[[1]][[2]] sp_n_s <- mean(sp_post$sigma1) sp_v_s <- mean(sp_post$sigma1 * sp_post$sigma1_vac) boost_to_set_point <- function(x, th1, th2, th3, alpha) { alpha * (th3 / (th2 + th3) * (1- exp(-(th2 + th3)* x)) - (th1 - th3) / (th1 - th2 - th3) * (exp(-th1 * x) - exp(- (th2 + th3) * x))) } mu_B_SP <- sp_post$mu[100] th1_SP <- sp_post$th1_SP[100] th2_SP <- sp_post$th2_SP[100] th3_SP <- sp_post$th3_SP[100] alpha_SP <- sp_post$alpha_SP[100] data_plot_all <- list() for(i in 1:52){ data_plot_all[[i]] <- data.frame( x = i, y = 1968:2016, lognormal = get_rel_vec(ln_v_s) * functionalLN(i, mu_LN, sigma_LN, alpha_LN), setpoint = get_rel_vec(sp_v_s) * boost_to_set_point(i, th1_SP, th2_SP, th3_SP, alpha_SP) ) } df_data_plot_all <- data_plot_all %>% bind_rows df_data_plot_all_long <- df_data_plot_all %>% pivot_longer(!c(x,y), names_to = "model_type", values_to = "titre_boost") df_data_plot_all_long %>% ggplot() + geom_tile(aes(x = x, y = y, fill = titre_boost), color = "black") + scale_fill_gradient(low = "white", high = "red") + labs(x = "Weeks post vaccination", y = "Virus year", fill = "Titre boost") + facet_grid(vars(model_type)) ggsave(file = here::here("outputs", "hanam", "analysis", "figs", "figE.pdf"), height = 5, width = 8) i_vals <- seq(1, 52 * 10, 52) data_plot_all_nat <- list() for(i in 1:10){ data_plot_all_nat[[i]] <- data.frame( x = i, y = 1968:2016, lognormal_nat = get_rel_vec(ln_n_s) * mu_B_LN, sp_nat = get_rel_vec(sp_n_s) * mu_B_SP, lognormal_vac = get_rel_vec(ln_v_s) * functionalLN(i_vals[i], mu_LN, sigma_LN, alpha_LN), sp_vac = get_rel_vec(sp_v_s) * boost_to_set_point(i_vals[i], th1_SP, th2_SP, th3_SP, alpha_SP) ) } df_data_plot_all_nat <- data_plot_all_nat %>% bind_rows df_data_plot_all_nat_long <- df_data_plot_all_nat %>% pivot_longer(!c(x,y), names_to = "model_type", values_to = "titre_boost") df_data_plot_all_nat_long %>% ggplot() + geom_tile(aes(x = x, y = y, fill = titre_boost), color = "black") + scale_fill_gradient(low = "white", high = "blue") + labs(x = "Years post vaccination", y = "Virus year", fill = "Titre boost") + facet_grid(vars(model_type)) ggsave(file = here::here("outputs", "hanam", "analysis", "figs", "figF.pdf"), height = 5, width = 8)
model_no <- 3 chain_short <- as.data.frame(all_post_test_full_raw$theta_chain) %>% filter(chain_no == 1) infection_histories <- all_post_test_full_raw$inf_chain %>% filter(chain_no == 1) titre_dat <- all_models_hanam[[model_no]]$create_post$titre_data vaccination_histories <- all_models_hanam[[model_no]]$create_post$vaccination_histories antigenic_map <- all_models_hanam[[model_no]]$create_post$antigenic_map strain_isolation_times <- all_models_hanam[[model_no]]$create_post$strain_isolation_times par_tab <- all_models_hanam[[model_no]]$create_post$par_tab vaccination_histories <- all_models_hanam[[model_no]]$create_post$vac_history custom_ab_kin_func <- all_models_hanam[[model_no]]$custom_ab_kin_func custom_ab_kin_func_1 <- all_models_hanam[[model_no]]$custom_ab_kin_func custom_antigenic_maps_func <- all_models_hanam[[model_no]]$custom_antigenic_maps_func expand_titredat <- FALSE nsamp <- 100 titre_before_infection <- FALSE measurement_indices_by_time <- NULL individuals <- unique(titre_dat$individual) mu_indices <- NULL plot_infection_histories_long( chain = chain_short, infection_histories = infection_histories, vaccination_histories = vaccination_histories, titre_dat = titre_dat, individuals = 1, antigenic_map = antigenic_map, strain_isolation_times = strain_isolation_times, par_tab = par_tab, custom_ab_kin_func = custom_ab_kin_func, custom_antigenic_maps_func = custom_antigenic_maps_func, size_plot = 25)
# Log normal load(file = here::here("data", paste0("modelinfo_cross_", hanam$study_name_short, ".RDS"))) load(file = here::here("outputs", "hanam", "fits", "dis_dep", "full", "LN_post_full.Rdata")) all_post_test_full_raw_ln <- post_cross$post_raw all_post_test_full_ln <- post_cross$post_trim load(file = here::here("outputs", "hanam", "fits", "dis_dep", "full", "SP_post_full.Rdata")) all_post_test_full_raw_sp <- post_cross$post_raw all_post_test_full_sp <- post_cross$post_trim inf_chain <- all_post_test_full_raw_sp$inf_chain sample_nos <- inf_chain %>% pull(sampno) %>% unique %>% length chain <- inf_chain %>% pull(chain_no) %>% max sample_df <- inf_chain %>% filter(sampno > 100000) %>% group_by(year = j, sampno) %>% dplyr::summarise(x = sum(x / (chain * 101))) %>% as.data.frame %>% filter(year < 49) sample_df %>% ggplot() + geom_boxplot(aes(x = as.character(year + 1967), y = x)) + scale_x_discrete(guide = guide_axis(check.overlap = TRUE)) + labs(x = "Year", y = "Attack rate", title = "Estimated attack rate")
load(file = here::here("data", paste0("modelinfo_full_", hcwpre$study_name_short, ".RDS"))) #all_models_hcw_pre_cross load(file = here::here("outputs", "hcw_pre", "fits", "dis_dep", "full", "trimmed_post_full.Rdata")) pos <- 5 all_post_test_full <- post_full$post_raw chain_full <- as.data.frame(all_post_test_full[[pos]]$theta_chain) %>% filter(chain_no == 1) infection_histories <- all_post_test_full[[pos]]$inf_chain %>% filter(chain_no == 1) titre_dat <- all_models_hcw_pre_full[[pos]]$create_post$titre_data vaccination_histories <- all_models_hcw_pre_full[[pos]]$create_post$vaccination_histories antigenic_map <- all_models_hcw_pre_full[[pos]]$create_post$antigenic_map strain_isolation_times <- all_models_hcw_pre_full[[pos]]$create_post$strain_isolation_times par_tab <- all_models_hcw_pre_full[[pos]]$create_post$par_tab vaccination_histories <- all_models_hcw_pre_full[[pos]]$create_post$vac_history custom_ab_kin_func <- all_models_hcw_pre_full[[pos]]$custom_ab_kin_func custom_ab_kin_func_1 <- all_models_hcw_pre_full[[pos]]$custom_ab_kin_func custom_antigenic_maps_func <- all_models_hcw_pre_full[[pos]]$custom_antigenic_maps_func expand_titredat <- FALSE nsamp <- 100 titre_before_infection <- FALSE measurement_indices_by_time <- NULL individuals <- unique(titre_dat$individual) mu_indices <- NULL chain_full <- chain_full %>% filter(sampno > 300000) y <- get_titre_predictions( chain = chain_full, infection_histories = infection_histories, vaccination_histories = vaccination_histories, titre_dat = titre_dat, individuals = 1:49, antigenic_map = antigenic_map, strain_isolation_times = strain_isolation_times, par_tab = par_tab, custom_ab_kin_func = custom_ab_kin_func, custom_antigenic_maps_func = custom_antigenic_maps_func, nsamp = 100 ) plot_infection_histories_long( chain = chain_full, infection_histories = infection_histories, vaccination_histories = vaccination_histories, titre_dat = titre_dat, individuals = 21, antigenic_map = antigenic_map, strain_isolation_times = strain_isolation_times, par_tab = par_tab, custom_ab_kin_func = custom_ab_kin_func, custom_antigenic_maps_func = custom_antigenic_maps_func ) chain_full$mu_vac %>% mean chain_full$mu_short_vac %>% mean
sum_data <- titre_dat %>% group_by(sample_time, titre) %>% summarise(n = n()) weighted_sum_data <- sum_data %>% summarise(w_mean = weighted.mean(titre, n, na.rm = TRUE)) #predictions titre_pred <- y$predicted_observations titre_pred <- titre_pred %>% filter(virus == 24192) titre_pred %>% select(virus, samples, titre, median) sum_data <- titre_pred %>% group_by(sample_time, titre) %>% summarise(n = n()) weighted_sum_data <- sum_data %>% summarise(w_mean = weighted.mean(titre, n, na.rm = TRUE)) sum_pred <- titre_pred %>% group_by(sample_time, median) %>% summarise(n = n()) weighted_sum_pred <- sum_pred %>% summarise(w_mean = weighted.mean(median, n, na.rm = TRUE))
Post plots
all_post_test_full_trim <- post_full$post_trim postplots_long <- all_post_test_full_trim[[5]] %>% bind_rows dec <- 1:12 %>% map(~ (postplots_long[["mu_short_vac"]] * (1.0 - ((1 / 12 + postplots_long[["wane_vac"]]) * .x ))) ) df_dec <- dec %>% setNames(as.character(c(1:12))) %>% as.data.frame %>% pivot_longer(c(1:12), values_to = "titre", names_to = "time") %>% mutate(time = str_extract(time, "[0-9,.]+")) %>% mutate(time = factor(time, levels = as.character(c(1:12)))) df_dec <- df_dec[df_dec$titre > 0, ] sero_data <- hcwpre$sero_data sum_data <- sero_data %>% group_by(sample_time, titre_value) %>% summarise(n = n()) weighted_sum_data <- sum_data %>% summarise(w_mean = weighted.mean(titre_value, n, na.rm = TRUE)) %>% mutate(sample_time = sample_time / 30) ggplot() + geom_point(data = sum_data, aes(x = sample_time, y = titre_value, size = n), fill = "black", alpha = 0.5) + geom_line(data = weighted_sum_data, aes(x = sample_time, y = w_mean), color = "red", size = 1) + geom_point(data = weighted_sum_data, aes(x = sample_time, y = w_mean), shape = 21, fill = "red", size = 5) + labs(x = "Time post vaccination", y = "Titre value", size = "Number of samples") df_dec %>% ggplot() + geom_boxplot(aes(x = time, y = titre)) + labs(x = "Months post vaccination") + geom_point(data = weighted_sum_data, aes(x = sample_time, y = w_mean), shape = 8, color = "red", size = 5) + geom_hline(yintercept = postplots_long[["mu"]] * postplots_long[["mu_vac"]] %>% mean)
antigenic_map_melted <- c(melt_antigenic_coords(antigenic_map[, c("x_coord", "y_coord")])) chain$sigma1 %>% mean s <- 0.01 lapply(antigenic_map_melted, function(x) max(1 - x * s, 0)) %>% unlist
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.