updated analysis
Will King
3 years ago
parent
a7e1f71e12
commit
98fc18e200
@ -0,0 +1,783 @@
|
|||||||
|
---
|
||||||
|
title: "The Effects of market conditions on recruitment and completion of clinical trials"
|
||||||
|
author: "Will King"
|
||||||
|
format: html
|
||||||
|
editor: source
|
||||||
|
---
|
||||||
|
|
||||||
|
|
||||||
|
# Setup
|
||||||
|
|
||||||
|
```{r}
|
||||||
|
library(bayesplot)
|
||||||
|
available_mcmc(pattern = "_nuts_")
|
||||||
|
library(ggplot2)
|
||||||
|
library(patchwork)
|
||||||
|
library(tidyverse)
|
||||||
|
library(rstan)
|
||||||
|
library(tidyr)
|
||||||
|
library(ghibli)
|
||||||
|
#Resources: https://github.com/stan-dev/rstan/wiki/RStan-Getting-Started
|
||||||
|
|
||||||
|
#save unchanged models instead of recompiling
|
||||||
|
rstan_options(auto_write = TRUE)
|
||||||
|
#allow for multithreaded sampling
|
||||||
|
options(mc.cores = parallel::detectCores())
|
||||||
|
|
||||||
|
#test installation, shouldn't get any errors
|
||||||
|
#example(stan_model, package = "rstan", run.dontrun = TRUE)
|
||||||
|
```
|
||||||
|
|
||||||
|
```{r}
|
||||||
|
################ Pull data from database ######################
|
||||||
|
library(RPostgreSQL)
|
||||||
|
|
||||||
|
driver <- dbDriver("PostgreSQL")
|
||||||
|
|
||||||
|
get_data <- function(driver) {
|
||||||
|
|
||||||
|
con <- dbConnect(
|
||||||
|
driver,
|
||||||
|
user='root',
|
||||||
|
password='root',
|
||||||
|
dbname='aact_db',
|
||||||
|
host='will-office'
|
||||||
|
)
|
||||||
|
on.exit(dbDisconnect(con))
|
||||||
|
|
||||||
|
query <- dbSendQuery(
|
||||||
|
con,
|
||||||
|
# "select * from formatted_data_with_planned_enrollment;"
|
||||||
|
"
|
||||||
|
select
|
||||||
|
fdqpe.nct_id
|
||||||
|
--,fdqpe.start_date
|
||||||
|
--,fdqpe.current_enrollment
|
||||||
|
--,fdqpe.enrollment_category
|
||||||
|
,fdqpe.current_status
|
||||||
|
,fdqpe.earliest_date_observed
|
||||||
|
,fdqpe.elapsed_duration
|
||||||
|
,fdqpe.n_brands as identical_brands
|
||||||
|
,ntbtu.brand_name_count
|
||||||
|
,fdqpe.category_id
|
||||||
|
,fdqpe.final_status
|
||||||
|
,fdqpe.h_sdi_val
|
||||||
|
--,fdqpe.h_sdi_u95
|
||||||
|
--,fdqpe.h_sdi_l95
|
||||||
|
,fdqpe.hm_sdi_val
|
||||||
|
--,fdqpe.hm_sdi_u95
|
||||||
|
--,fdqpe.hm_sdi_l95
|
||||||
|
,fdqpe.m_sdi_val
|
||||||
|
--,fdqpe.m_sdi_u95
|
||||||
|
--,fdqpe.m_sdi_l95
|
||||||
|
,fdqpe.lm_sdi_val
|
||||||
|
--,fdqpe.lm_sdi_u95
|
||||||
|
--,fdqpe.lm_sdi_l95
|
||||||
|
,fdqpe.l_sdi_val
|
||||||
|
--,fdqpe.l_sdi_u95
|
||||||
|
--,fdqpe.l_sdi_l95
|
||||||
|
from formatted_data_with_planned_enrollment fdqpe
|
||||||
|
join \"Formularies\".nct_to_brands_through_uspdc ntbtu
|
||||||
|
on fdqpe.nct_id = ntbtu.nct_id
|
||||||
|
order by fdqpe.nct_id, fdqpe.earliest_date_observed
|
||||||
|
;
|
||||||
|
"
|
||||||
|
)
|
||||||
|
df <- fetch(query, n = -1)
|
||||||
|
df <- na.omit(df)
|
||||||
|
|
||||||
|
query2 <-dbSendQuery(con,"select count(*) from \"DiseaseBurden\".icd10_categories ic where \"level\"=1;")
|
||||||
|
n_categories <- fetch(query2, n = -1)
|
||||||
|
|
||||||
|
return(list(data=df,ncat=n_categories))
|
||||||
|
}
|
||||||
|
|
||||||
|
d <- get_data(driver)
|
||||||
|
df <- d$data
|
||||||
|
n_categories <- d$ncat
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
################ Format Data ###########################
|
||||||
|
|
||||||
|
data_formatter <- function(df) {
|
||||||
|
categories <- df["category_id"]
|
||||||
|
|
||||||
|
x <- df["elapsed_duration"]
|
||||||
|
x["identical_brands"] <- asinh(df$identical_brands)
|
||||||
|
x["brand_name_counts"] <- asinh(df$brand_name_count)
|
||||||
|
x["h_sdi_val"] <- asinh(df$h_sdi_val)
|
||||||
|
x["hm_sdi_val"] <- asinh(df$hm_sdi_val)
|
||||||
|
x["m_sdi_val"] <- asinh(df$m_sdi_val)
|
||||||
|
x["lm_sdi_val"] <- asinh(df$lm_sdi_val)
|
||||||
|
x["l_sdi_val"] <- asinh(df$l_sdi_val)
|
||||||
|
|
||||||
|
|
||||||
|
#interaction terms for competitors
|
||||||
|
x["ib*elapsed"] <- x["elapsed_duration"]*x["identical_brands"]
|
||||||
|
x["bnc*elapsed"] <- x["elapsed_duration"] * x["brand_name_counts"]
|
||||||
|
|
||||||
|
y <- ifelse(df["final_status"]=="Terminated",1,0)
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
return(list(x=x,y=y))
|
||||||
|
}
|
||||||
|
|
||||||
|
train <- data_formatter(df)
|
||||||
|
|
||||||
|
categories <- df$category_id
|
||||||
|
|
||||||
|
x <- train$x
|
||||||
|
y <- train$y
|
||||||
|
```
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
### Intervention: Adding a single competitor
|
||||||
|
```{r}
|
||||||
|
inherited_cols <- c(
|
||||||
|
"elapsed_duration"
|
||||||
|
,"identical_brands"
|
||||||
|
,"brand_name_counts"
|
||||||
|
,"h_sdi_val"
|
||||||
|
,"hm_sdi_val"
|
||||||
|
,"m_sdi_val"
|
||||||
|
,"lm_sdi_val"
|
||||||
|
,"l_sdi_val"
|
||||||
|
,"status_NYR"
|
||||||
|
,"status_EBI"
|
||||||
|
,"status_ANR"
|
||||||
|
,"ib*elapsed"
|
||||||
|
,"bnc*elapsed"
|
||||||
|
)
|
||||||
|
```
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
#### Generics
|
||||||
|
|
||||||
|
```{r}
|
||||||
|
#generics intervention
|
||||||
|
brand_intervention_ib <- x[inherited_cols]
|
||||||
|
brand_intervention_ib["identical_brands"] <- asinh(sinh(x$identical_brands)+1) #add a single generic brand
|
||||||
|
brand_intervention_ib["ib*elapsed"] <- brand_intervention_ib$identical_brands * brand_intervention_ib$elapsed_duration
|
||||||
|
```
|
||||||
|
|
||||||
|
```{r}
|
||||||
|
counterfact_marketing_ib <- list(
|
||||||
|
D = ncol(x),#
|
||||||
|
N = nrow(x),
|
||||||
|
L = n_categories$count,
|
||||||
|
y = as.vector(y),
|
||||||
|
ll = as.vector(categories),
|
||||||
|
x = as.matrix(x),
|
||||||
|
mu_mean = 0,
|
||||||
|
mu_stdev = 0.05,
|
||||||
|
sigma_shape = 4,
|
||||||
|
sigma_rate = 20,
|
||||||
|
Nx = nrow(x),
|
||||||
|
llx = as.vector(categories),
|
||||||
|
counterfact_x_tilde = as.matrix(brand_intervention_ib),
|
||||||
|
counterfact_x = as.matrix(x)
|
||||||
|
)
|
||||||
|
```
|
||||||
|
|
||||||
|
```{r}
|
||||||
|
generated_bi <- gqs(
|
||||||
|
fit@stanmodel,
|
||||||
|
data=counterfact_marketing_ib,
|
||||||
|
draws=as.matrix(fit),
|
||||||
|
seed=11021585
|
||||||
|
)
|
||||||
|
```
|
||||||
|
|
||||||
|
### USP DC
|
||||||
|
|
||||||
|
|
||||||
|
```{r}
|
||||||
|
#formulary intervention
|
||||||
|
brand_intervention_bnc <- x[c(inherited_cols,"identical_brands","ib*elapsed")]
|
||||||
|
brand_intervention_bnc["brand_name_counts"] <- asinh(sinh(x$brand_name_counts)+1) #add a single formulary competitor brand
|
||||||
|
brand_intervention_bnc["bnc*elapsed"] <- brand_intervention_bnc$brand_name_counts * brand_intervention_bnc$elapsed_duration
|
||||||
|
```
|
||||||
|
|
||||||
|
```{r}
|
||||||
|
counterfact_marketing_bnc <- list(
|
||||||
|
D = ncol(x),#
|
||||||
|
N = nrow(x),
|
||||||
|
L = n_categories$count,
|
||||||
|
y = as.vector(y),
|
||||||
|
ll = as.vector(categories),
|
||||||
|
x = as.matrix(x),
|
||||||
|
mu_mean = 0,
|
||||||
|
mu_stdev = 0.05,
|
||||||
|
sigma_shape = 4,
|
||||||
|
sigma_rate = 20,
|
||||||
|
Nx = nrow(x),
|
||||||
|
llx = as.vector(categories),
|
||||||
|
counterfact_x_tilde = as.matrix(brand_intervention_bnc),
|
||||||
|
counterfact_x = as.matrix(x)
|
||||||
|
)
|
||||||
|
```
|
||||||
|
|
||||||
|
|
||||||
|
```{r}
|
||||||
|
generated_bnc <- gqs(
|
||||||
|
fit@stanmodel,
|
||||||
|
data=counterfact_marketing_bnc,
|
||||||
|
draws=as.matrix(fit),
|
||||||
|
seed=11021585
|
||||||
|
)
|
||||||
|
```
|
||||||
|
|
||||||
|
|
||||||
|
# Fit Model
|
||||||
|
|
||||||
|
|
||||||
|
```{r}
|
||||||
|
################################# FIT MODEL #########################################
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
fit <- stan(
|
||||||
|
file='Hierarchal_Logistic.stan',
|
||||||
|
data = counterfact_marketing_ib,
|
||||||
|
chains = 4,
|
||||||
|
iter = 5000,
|
||||||
|
seed = 11021585
|
||||||
|
)
|
||||||
|
```
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
```{r}
|
||||||
|
hist(as.vector(extract(generated, pars="p_prior")$p_prior))
|
||||||
|
hist(as.vector(extract(generated, pars="mu_prior")$mu_prior), )
|
||||||
|
hist(as.vector(extract(generated, pars="sigma_prior")$sigma_prior))
|
||||||
|
```
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
```{r}
|
||||||
|
check_hmc_diagnostics(fit)
|
||||||
|
hist(as.vector(extract(generated, pars="p_predicted")$p_predicted))
|
||||||
|
```
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
# Diagnostics
|
||||||
|
|
||||||
|
```{r}
|
||||||
|
#trace plots
|
||||||
|
plot(fit, pars=c("mu"), plotfun="trace")
|
||||||
|
|
||||||
|
|
||||||
|
for (i in 1:4) {
|
||||||
|
print(
|
||||||
|
mcmc_rank_overlay(
|
||||||
|
fit,
|
||||||
|
pars=c(
|
||||||
|
paste0("mu[",4*i-3,"]"),
|
||||||
|
paste0("mu[",4*i-2,"]"),
|
||||||
|
paste0("mu[",4*i-1,"]"),
|
||||||
|
paste0("mu[",4*i,"]")
|
||||||
|
),
|
||||||
|
n_bins=100
|
||||||
|
)+ legend_move("top") +
|
||||||
|
scale_colour_ghibli_d("KikiMedium")
|
||||||
|
)
|
||||||
|
}
|
||||||
|
```
|
||||||
|
|
||||||
|
```{r}
|
||||||
|
plot(fit, pars=c("sigma"), plotfun="trace")
|
||||||
|
|
||||||
|
for (i in 1:4) {
|
||||||
|
print(
|
||||||
|
mcmc_rank_overlay(
|
||||||
|
fit,
|
||||||
|
pars=c(
|
||||||
|
paste0("sigma[",4*i-3,"]"),
|
||||||
|
paste0("sigma[",4*i-2,"]"),
|
||||||
|
paste0("sigma[",4*i-1,"]"),
|
||||||
|
paste0("sigma[",4*i,"]")
|
||||||
|
),
|
||||||
|
n_bins=100
|
||||||
|
)+ legend_move("top") +
|
||||||
|
scale_colour_ghibli_d("KikiMedium")
|
||||||
|
)
|
||||||
|
}
|
||||||
|
```
|
||||||
|
|
||||||
|
```{r}
|
||||||
|
#other diagnostics
|
||||||
|
logpost <- log_posterior(fit)
|
||||||
|
nuts_prmts <- nuts_params(fit)
|
||||||
|
posterior <- as.array(fit)
|
||||||
|
|
||||||
|
```
|
||||||
|
|
||||||
|
```{r}
|
||||||
|
color_scheme_set("darkgray")
|
||||||
|
div_style <- parcoord_style_np(div_color = "green", div_size = 0.05, div_alpha = 0.4)
|
||||||
|
mcmc_parcoord(posterior, regex_pars = "mu", np=nuts_prmts, np_style = div_style, alpha = 0.05)
|
||||||
|
```
|
||||||
|
|
||||||
|
```{r}
|
||||||
|
for (i in 1:4) {
|
||||||
|
mus = sapply(3:0, function(j) paste0("mu[",4*i-j ,"]"))
|
||||||
|
print(
|
||||||
|
mcmc_pairs(
|
||||||
|
posterior,
|
||||||
|
np = nuts_prmts,
|
||||||
|
pars=c(
|
||||||
|
mus,
|
||||||
|
"lp__"
|
||||||
|
),
|
||||||
|
off_diag_args = list(size = 0.75)
|
||||||
|
)
|
||||||
|
)
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
```
|
||||||
|
|
||||||
|
```{r}
|
||||||
|
mcmc_parcoord(posterior,regex_pars = "sigma", np=nuts_prmts, alpha=0.05)
|
||||||
|
```
|
||||||
|
|
||||||
|
```{r}
|
||||||
|
|
||||||
|
for (i in 1:4) {
|
||||||
|
params = sapply(3:0, function(j) paste0("sigma[",4*i-j ,"]"))
|
||||||
|
print(
|
||||||
|
mcmc_pairs(
|
||||||
|
posterior,
|
||||||
|
np = nuts_prmts,
|
||||||
|
pars=c(
|
||||||
|
params,
|
||||||
|
"lp__"
|
||||||
|
),
|
||||||
|
off_diag_args = list(size = 0.75)
|
||||||
|
)
|
||||||
|
)
|
||||||
|
}
|
||||||
|
```
|
||||||
|
|
||||||
|
|
||||||
|
```{r}
|
||||||
|
for (k in 1:22) {
|
||||||
|
for (i in 1:4) {
|
||||||
|
params = sapply(3:0, function(j) paste0("beta[",k,",",4*i-j ,"]"))
|
||||||
|
print(
|
||||||
|
mcmc_pairs(
|
||||||
|
posterior,
|
||||||
|
np = nuts_prmts,
|
||||||
|
pars=c(
|
||||||
|
params,
|
||||||
|
"lp__"
|
||||||
|
),
|
||||||
|
off_diag_args = list(size = 0.75)
|
||||||
|
)
|
||||||
|
)
|
||||||
|
}}
|
||||||
|
```
|
||||||
|
|
||||||
|
|
||||||
|
# Results
|
||||||
|
|
||||||
|
|
||||||
|
```{r}
|
||||||
|
################################# ANALYZE #####################################
|
||||||
|
print(fit)
|
||||||
|
```
|
||||||
|
|
||||||
|
## Result Plots
|
||||||
|
|
||||||
|
|
||||||
|
Note the regular large difference in variance.
|
||||||
|
I would guess those are the beta[1:22,2] values.
|
||||||
|
I wonder if a lot of the variance is due to the 2 values that are sitting out.
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
```{r}
|
||||||
|
beta_list <- list(
|
||||||
|
groups = c(
|
||||||
|
`1`="Infections & Parasites",
|
||||||
|
`2`="Neoplasms",
|
||||||
|
`3`="Blood & Immune system",
|
||||||
|
`4`="Endocrine, Nutritional, and Metabolic",
|
||||||
|
`5`="Mental & Behavioral",
|
||||||
|
`6`="Nervous System",
|
||||||
|
`7`="Eye and Adnexa",
|
||||||
|
`8`="Ear and Mastoid",
|
||||||
|
`9`="Circulatory",
|
||||||
|
`10`="Respiratory",
|
||||||
|
`11`="Digestive",
|
||||||
|
`12`="Skin & Subcutaneaous tissue",
|
||||||
|
`13`="Musculoskeletal",
|
||||||
|
`14`="Genitourinary",
|
||||||
|
`15`="Pregancy, Childbirth, & Puerperium",
|
||||||
|
`16`="Perinatal Period",
|
||||||
|
`17`="Congential",
|
||||||
|
`18`="Symptoms, Signs etc.",
|
||||||
|
`19`="Injury etc.",
|
||||||
|
`20`="External Causes",
|
||||||
|
`21`="Contact with Healthcare",
|
||||||
|
`22`="Special Purposes"
|
||||||
|
),
|
||||||
|
parameters = c(
|
||||||
|
`1`="Elapsed Duration",
|
||||||
|
# brands
|
||||||
|
`2`="asinh(Generic Brands)",
|
||||||
|
`3`="asinh(Competitors USPDC)",
|
||||||
|
# population
|
||||||
|
`4`="asinh(High SDI)",
|
||||||
|
`5`="asinh(High-Medium SDI)",
|
||||||
|
`6`="asinh(Medium SDI)",
|
||||||
|
`7`="asinh(Low-Medium SDI)",
|
||||||
|
`8`="asinh(Low SDI)",
|
||||||
|
#Status
|
||||||
|
`9`="status_NYR",
|
||||||
|
`10`="status_EBI",
|
||||||
|
`11`="status_ANR",
|
||||||
|
#interactions for brands
|
||||||
|
`12`="ib*elapsed",
|
||||||
|
`13`="bnc*elapsed",
|
||||||
|
# interactions for status
|
||||||
|
`14`="sNYR*elapsed",
|
||||||
|
`15`="sEBI*elapsed",
|
||||||
|
`16`="sANR*elapsed"
|
||||||
|
|
||||||
|
)
|
||||||
|
)
|
||||||
|
|
||||||
|
get_parameters <- function(stem,class_list) {
|
||||||
|
#get categories and lengths
|
||||||
|
named <- names(class_list)
|
||||||
|
lengths <- sapply(named, (function (x) length(class_list[[x]])))
|
||||||
|
|
||||||
|
#describe the grid needed
|
||||||
|
iter_list <- sapply(named, (function (x) 1:lengths[x]))
|
||||||
|
|
||||||
|
#generate the list of parameters
|
||||||
|
pardf <- generate_parameter_df(stem, iter_list)
|
||||||
|
|
||||||
|
#add columns with appropriate human-readable names
|
||||||
|
for (name in named) {
|
||||||
|
pardf[paste(name,"_hr",sep="")] <- as.factor(
|
||||||
|
sapply(pardf[name], (function (i) class_list[[name]][i]))
|
||||||
|
)
|
||||||
|
}
|
||||||
|
|
||||||
|
return(pardf)
|
||||||
|
}
|
||||||
|
|
||||||
|
generate_parameter_df <- function(stem, iter_list) {
|
||||||
|
grid <- expand.grid(iter_list)
|
||||||
|
grid["param_name"] <- grid %>% unite(x,colnames(grid),sep=",")
|
||||||
|
grid["param_name"] <- paste(stem,"[",grid$param_name,"]",sep="")
|
||||||
|
return(grid)
|
||||||
|
}
|
||||||
|
|
||||||
|
group_mcmc_areas <- function(
|
||||||
|
stem,# = "beta"
|
||||||
|
class_list,# = beta_list
|
||||||
|
stanfit,# = fit
|
||||||
|
group_id,# = 2
|
||||||
|
rename=TRUE
|
||||||
|
) {
|
||||||
|
#get all parameter names
|
||||||
|
params <- get_parameters(stem,class_list)
|
||||||
|
#filter down to parameters of interest
|
||||||
|
params <- filter(params,groups == group_id)
|
||||||
|
#Get dataframe with only the rows of interest
|
||||||
|
filtdata <- as.data.frame(stanfit)[params$param_name]
|
||||||
|
#rename columns
|
||||||
|
if (rename) dimnames(filtdata)[[2]] <- params$parameters_hr
|
||||||
|
#get group name for title
|
||||||
|
group_name <- class_list$groups[group_id]
|
||||||
|
#create area plot with appropriate title
|
||||||
|
mcmc_areas(filtdata,prob = 0.8, prob_outer = 0.95) +
|
||||||
|
ggtitle(paste("Parameter distributions for ICD-10 class:",group_name))
|
||||||
|
}
|
||||||
|
|
||||||
|
parameter_mcmc_areas <- function(
|
||||||
|
stem,# = "beta"
|
||||||
|
class_list,# = beta_list
|
||||||
|
stanfit,# = fit
|
||||||
|
parameter_id,# = 2
|
||||||
|
rename=TRUE
|
||||||
|
) {
|
||||||
|
#get all parameter names
|
||||||
|
params <- get_parameters(stem,class_list)
|
||||||
|
#filter down to parameters of interest
|
||||||
|
params <- filter(params,parameters == parameter_id)
|
||||||
|
#Get dataframe with only the rows of interest
|
||||||
|
filtdata <- as.data.frame(stanfit)[params$param_name]
|
||||||
|
#rename columns
|
||||||
|
if (rename) dimnames(filtdata)[[2]] <- params$groups_hr
|
||||||
|
#get group name for title
|
||||||
|
parameter_name <- class_list$parameters[parameter_id]
|
||||||
|
#create area plot with appropriate title
|
||||||
|
mcmc_areas(filtdata,prob = 0.8, prob_outer = 0.95) +
|
||||||
|
ggtitle(parameter_name,"Parameter Distribution")
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
```
|
||||||
|
|
||||||
|
```{r}
|
||||||
|
#mcmc_intervals(fit, pars=get_parameters("beta",beta_list)$param_name)
|
||||||
|
```
|
||||||
|
|
||||||
|
### Investigating parameter distributions
|
||||||
|
|
||||||
|
```{r}
|
||||||
|
#g1 <- group_mcmc_areas("beta",beta_list,fit,2)
|
||||||
|
#g2 <- group_mcmc_areas("beta",beta_list,fit,2)
|
||||||
|
#g3 <- group_mcmc_areas("beta",beta_list,fit,2)
|
||||||
|
#g4 <- group_mcmc_areas("beta",beta_list,fit,2)
|
||||||
|
#g5 <- group_mcmc_areas("beta",beta_list,fit,2)
|
||||||
|
#g6 <- group_mcmc_areas("beta",beta_list,fit,2)
|
||||||
|
#g7 <- group_mcmc_areas("beta",beta_list,fit,2)
|
||||||
|
#g8 <- group_mcmc_areas("beta",beta_list,fit,2)
|
||||||
|
#g9 <- group_mcmc_areas("beta",beta_list,fit,2)
|
||||||
|
#g10 <- group_mcmc_areas("beta",beta_list,fit,2)
|
||||||
|
#g11 <- group_mcmc_areas("beta",beta_list,fit,2)
|
||||||
|
#g12 <- group_mcmc_areas("beta",beta_list,fit,2)
|
||||||
|
#g13 <- group_mcmc_areas("beta",beta_list,fit,2)
|
||||||
|
#g14 <- group_mcmc_areas("beta",beta_list,fit,2)
|
||||||
|
#g15 <- group_mcmc_areas("beta",beta_list,fit,2)
|
||||||
|
#g16 <- group_mcmc_areas("beta",beta_list,fit,2)
|
||||||
|
#g17 <- group_mcmc_areas("beta",beta_list,fit,2)
|
||||||
|
#g18 <- group_mcmc_areas("beta",beta_list,fit,2)
|
||||||
|
#g19 <- group_mcmc_areas("beta",beta_list,fit,2)
|
||||||
|
#g20 <- group_mcmc_areas("beta",beta_list,fit,2)
|
||||||
|
#g21 <- group_mcmc_areas("beta",beta_list,fit,2)
|
||||||
|
#g22 <- group_mcmc_areas("beta",beta_list,fit,2)
|
||||||
|
|
||||||
|
|
||||||
|
p1 <- parameter_mcmc_areas("beta",beta_list,fit,1)
|
||||||
|
p2 <- parameter_mcmc_areas("beta",beta_list,fit,2)
|
||||||
|
p3 <- parameter_mcmc_areas("beta",beta_list,fit,3)
|
||||||
|
#p4 <- parameter_mcmc_areas("beta",beta_list,fit,4)
|
||||||
|
#p5 <- parameter_mcmc_areas("beta",beta_list,fit,5)
|
||||||
|
#p6 <- parameter_mcmc_areas("beta",beta_list,fit,6)
|
||||||
|
#p7 <- parameter_mcmc_areas("beta",beta_list,fit,7)
|
||||||
|
#p8 <- parameter_mcmc_areas("beta",beta_list,fit,8)
|
||||||
|
p9 <- parameter_mcmc_areas("beta",beta_list,fit,9)
|
||||||
|
p10 <- parameter_mcmc_areas("beta",beta_list,fit,10)
|
||||||
|
p11 <- parameter_mcmc_areas("beta",beta_list,fit,11)
|
||||||
|
p12 <- parameter_mcmc_areas("beta",beta_list,fit,12)
|
||||||
|
p13 <- parameter_mcmc_areas("beta",beta_list,fit,13)
|
||||||
|
p14 <- parameter_mcmc_areas("beta",beta_list,fit,14)
|
||||||
|
p15 <- parameter_mcmc_areas("beta",beta_list,fit,15)
|
||||||
|
p16 <- parameter_mcmc_areas("beta",beta_list,fit,16)
|
||||||
|
```
|
||||||
|
|
||||||
|
Note these have 95% outer CI and 80% inner (shaded)
|
||||||
|
|
||||||
|
|
||||||
|
1) "Elapsed Duration",
|
||||||
|
2) "asinh(Generic Brands)",
|
||||||
|
3) "asinh(Competitors USPDC)",
|
||||||
|
4) "asinh(High SDI)",
|
||||||
|
5) "asinh(High-Medium SDI)",
|
||||||
|
6) "asinh(Medium SDI)",
|
||||||
|
7) "asinh(Low-Medium SDI)",
|
||||||
|
8) "asinh(Low SDI)",
|
||||||
|
9) "status_NYR",
|
||||||
|
10) "status_EBI",
|
||||||
|
11) "status_ANR",
|
||||||
|
12) "ib*elapsed",
|
||||||
|
13) "bnc*elapsed",
|
||||||
|
14) "sNYR*elapsed",
|
||||||
|
15) "sEBI*elapsed",
|
||||||
|
16) "sANR*elapsed"
|
||||||
|
|
||||||
|
of interest
|
||||||
|
- p1 + p2
|
||||||
|
- p3 + p2
|
||||||
|
- p2 + p12
|
||||||
|
- p3 + p13
|
||||||
|
- p9 + p14
|
||||||
|
- p10 + p15
|
||||||
|
- p11 + p16
|
||||||
|
|
||||||
|
```{r}
|
||||||
|
p1 + p2
|
||||||
|
```
|
||||||
|
|
||||||
|
|
||||||
|
```{r}
|
||||||
|
p2 + p3
|
||||||
|
```
|
||||||
|
|
||||||
|
```{r}
|
||||||
|
p2 + p12
|
||||||
|
```
|
||||||
|
|
||||||
|
```{r}
|
||||||
|
p3 + p13
|
||||||
|
```
|
||||||
|
|
||||||
|
```{r}
|
||||||
|
p9 + p14
|
||||||
|
```
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
```{r}
|
||||||
|
p10 + p15
|
||||||
|
```
|
||||||
|
|
||||||
|
|
||||||
|
```{r}
|
||||||
|
p11 + p16
|
||||||
|
```
|
||||||
|
|
||||||
|
|
||||||
|
# Posterior Prediction
|
||||||
|
|
||||||
|
|
||||||
|
```{r}
|
||||||
|
#TODO: Convert to ggplot, stabilize y axis
|
||||||
|
hist(as.vector(extract(generated, pars="p_predicted_default")$p_predicted_default))
|
||||||
|
hist(as.vector(extract(generated, pars="p_predicted_intervention")$p_predicted_intervention))
|
||||||
|
```
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
## Distribution of Predicted Differences
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
### Intervention: Adding a single competitor
|
||||||
|
|
||||||
|
|
||||||
|
#### Generics
|
||||||
|
|
||||||
|
|
||||||
|
```{r}
|
||||||
|
#TODO: Convert to ggplot, stabilize y axis
|
||||||
|
hist(as.vector(extract(generated_ib, pars="p_predicted_default")$p_predicted_default), bins=100)
|
||||||
|
hist(as.vector(extract(generated_ib, pars="p_predicted_intervention")$p_predicted_intervention), bins=100)
|
||||||
|
hist(as.vector(extract(generated_ib, pars="predicted_difference")$predicted_difference), bins=100)
|
||||||
|
```
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
```{r}
|
||||||
|
|
||||||
|
|
||||||
|
pddf_ib <- data.frame(extract(generated_ib, pars="predicted_difference")$predicted_difference) |>
|
||||||
|
pivot_longer(X1:X1343)
|
||||||
|
|
||||||
|
#TODO: Fix Category names
|
||||||
|
pddf_ib["entry_idx"] <- as.numeric(gsub("\\D","",pddf_ib$name))
|
||||||
|
pddf_ib["category"] <- sapply(pddf_ib$entry_idx, function(i) df$category_id[i])
|
||||||
|
pddf_ib["category_name"] <- sapply(pddf_ib$category, function(i) beta_list$groups[i])
|
||||||
|
```
|
||||||
|
|
||||||
|
|
||||||
|
```{r}
|
||||||
|
|
||||||
|
ggplot(pddf_ib, aes(x=value,)) +
|
||||||
|
geom_histogram(bins=100) +
|
||||||
|
labs(
|
||||||
|
title = "Distribution of predicted differences"
|
||||||
|
,x = "Difference in probability due to intervention"
|
||||||
|
,y = "Predicted counts"
|
||||||
|
) +
|
||||||
|
#xlim(-0.3,0.1) +
|
||||||
|
geom_vline(aes(xintercept = 0), color = "skyblue", linetype="dashed")
|
||||||
|
|
||||||
|
ggplot(pddf_ib, aes(x=value,)) +
|
||||||
|
geom_histogram(bins=100) +
|
||||||
|
facet_wrap(
|
||||||
|
~factor(
|
||||||
|
category_name,
|
||||||
|
levels=beta_list$groups
|
||||||
|
)
|
||||||
|
, labeller = label_wrap_gen(multi_line = TRUE)
|
||||||
|
, ncol=5) +
|
||||||
|
labs(
|
||||||
|
title = "Distribution of predicted differences",
|
||||||
|
subtitle = "By group"
|
||||||
|
,x = "Difference in probability due to intervention"
|
||||||
|
,y = "Predicted counts"
|
||||||
|
) +
|
||||||
|
#xlim(-0.25,0.1) +
|
||||||
|
geom_vline(aes(xintercept = 0), color = "skyblue", linetype="dashed") +
|
||||||
|
theme(strip.text.x = element_text(size = 8))
|
||||||
|
|
||||||
|
```
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
#### USP DC
|
||||||
|
```{r}
|
||||||
|
#TODO: Convert to ggplot, stabilize y axis
|
||||||
|
hist(as.vector(extract(generated_bnc, pars="p_predicted_default")$p_predicted_default), bins=100)
|
||||||
|
hist(as.vector(extract(generated_bnc, pars="p_predicted_intervention")$p_predicted_intervention), bins=100)
|
||||||
|
hist(as.vector(extract(generated_bnc, pars="predicted_difference")$predicted_difference), bins=100)
|
||||||
|
```
|
||||||
|
|
||||||
|
|
||||||
|
```{r}
|
||||||
|
pddf_bnc <- data.frame(extract(generated_bnc, pars="predicted_difference")$predicted_difference) |>
|
||||||
|
pivot_longer(X1:X1343)
|
||||||
|
|
||||||
|
#Add Category names
|
||||||
|
pddf_bnc["entry_idx"] <- as.numeric(gsub("\\D","",pddf_bnc$name))
|
||||||
|
pddf_bnc["category"] <- sapply(pddf_bnc$entry_idx, function(i) df$category_id[i])
|
||||||
|
pddf_bnc["category_name"] <- sapply(pddf_bnc$category, function(i) beta_list$groups[i])
|
||||||
|
|
||||||
|
#add snapshot date
|
||||||
|
pddf_bnc["snapshot_date"] <- sapply(pddf_bnc$entry_idx, function(i) as_date(df$earliest_date_observed[i]))
|
||||||
|
```
|
||||||
|
|
||||||
|
|
||||||
|
```{r}
|
||||||
|
|
||||||
|
ggplot(pddf_bnc, aes(x=value,)) +
|
||||||
|
geom_histogram(bins=100) +
|
||||||
|
labs(
|
||||||
|
title = "Distribution of predicted differences"
|
||||||
|
,x = "Difference in probability due to intervention"
|
||||||
|
,y = "Predicted counts"
|
||||||
|
) +
|
||||||
|
#xlim(-0.3,0.1) +
|
||||||
|
geom_vline(aes(xintercept = 0), color = "skyblue", linetype="dashed")
|
||||||
|
|
||||||
|
ggplot(pddf_bnc, aes(x=value,)) +
|
||||||
|
geom_histogram(bins=100) +
|
||||||
|
facet_wrap(
|
||||||
|
~factor(
|
||||||
|
category_name,
|
||||||
|
levels=beta_list$groups
|
||||||
|
)
|
||||||
|
, labeller = label_wrap_gen(multi_line = TRUE)
|
||||||
|
, ncol=5) +
|
||||||
|
labs(
|
||||||
|
title = "Distribution of predicted differences",
|
||||||
|
subtitle = "By group"
|
||||||
|
,x = "Difference in probability due to intervention"
|
||||||
|
,y = "Predicted counts"
|
||||||
|
) +
|
||||||
|
#xlim(-0.25,0.1) +
|
||||||
|
geom_vline(aes(xintercept = 0), color = "skyblue", linetype="dashed") +
|
||||||
|
theme(strip.text.x = element_text(size = 8))
|
||||||
|
|
||||||
|
```
|
||||||
|
|
||||||
Loading…
Reference in New Issue