hiv_ex.rmd

---
title: "Bayesian CPM for HIV example"
output:
  html_document:
    toc: no
    toc_depth: 3
    number_sections: false
    code_folding: hide
    theme: paper
---

<!-- HIV Biomarker example -->

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)

rm(list=ls())
libs <- c("rstan", "brms" ,"dplyr", "stringr", "readr", "tidyr", "magrittr", "bayesCPM","bayesplot","cowplot","loo")
invisible(lapply(libs, library, character.only = TRUE))

dir <- file.path("~/Dropbox/njames/school/PhD/orise_ra/bayes_cpm")

# figdir<-file.path(dir,"biostat_sem","fig") # old
figdir<-file.path(dir,"paper","fig") # for paper

# read in functions
#source(file.path(dir,"cpm_functions.r"))

set.seed(88711)

# call this once to distribute MCMC chains across cpu cores:
options(mc.cores=parallel::detectCores())
```

## HIV example using Stan and linear regression

```{r hiv_dat}
datdir <- file.path(dir,"data")
hivdat0 <- read_csv(file.path(datdir,"koethe-biomarker-data.csv"))

# make censored and transformed outcomes, center predictors 
hivdat0 %<>% mutate(il_6_cens=if_else(il_6==0,-1,0),
                    sqrt_il_6=sqrt(il_6),
                    il_1_beta_cens=if_else(il_1_beta==0,-1,0),
                    bmi_c = bmi-mean(bmi),
                    age_c = age-mean(age),
                    cd4_c = cd4-mean(cd4),
                    cd4_100_c= (cd4_c)/100)

# change outcomes to ordered factors
hivdat <- hivdat0 %>% mutate(il_6=ordered(il_6),
                             il_1_beta=ordered(il_1_beta),
                             hscrp=ordered(hscrp), 
                             cd14=ordered(cd14),
                             leptin=ordered(leptin))

# for plots
seq_ln <- 30
newdat1 <- data.frame(bmi_c=seq(-10,15,length=seq_ln),
                     age_c=rep(0,seq_ln),
                     male=rep(1,seq_ln),
                     nonwhite=rep(0,seq_ln),
                     study.linc=rep(1,seq_ln),
                     cd4_100_c=rep(0,seq_ln),
                     smoker=rep(0,seq_ln))

```

```{r}
summary(hivdat0)
hivdat0
ggplot(data=hivdat0)+geom_histogram(aes(x=il_1_beta))
ggplot(data=hivdat0)+geom_histogram(aes(x=il_6))
ggplot(data=hivdat0)+geom_histogram(aes(x=leptin))
```

### with brms

```{r, eval=FALSE}
# censored normal linear regression with outcome transformation for conditional mean
brm_fit_il_6 <- brm(sqrt_il_6|cens(il_6_cens) ~ bmi_c + age_c + male + nonwhite +
                      study.linc + cd4_100_c + smoker,
                    data = hivdat0, family=gaussian)

# two censored quantile regression models for conditional median and 90th quantile
# ! not sure these are correct
brm_fit_il_6_q50 <-brm(bf(sqrt_il_6|cens(il_6_cens) ~ bmi_c + age_c + male + 
                            nonwhite + study.linc + cd4_100_c + smoker,
                          quantile = 0.5), data = hivdat0, 
                       family = asym_laplace())

brm_fit_il_6_q90 <-brm(bf(sqrt_il_6|cens(il_6_cens) ~ bmi_c + age_c + male + 
                            nonwhite + study.linc + cd4_100_c + smoker,
                          quantile = 0.1), data = hivdat0, 
                       family = asym_laplace())


brm_il_6_mn_dat0<-fitted(brm_fit_il_6,newdata=newdat1)
brm_il_6_md_dat0<-fitted(brm_fit_il_6_q50,newdata=newdat1)
brm_il_6_q90_dat0<-fitted(brm_fit_il_6_q90,newdata=newdat1)

brm_il_6_mn_dat<-brm_il_6_mn_dat0 %>% 
  as_tibble() %>% mutate(bmi_c=seq(-10,15,length=seq_ln)) %>%
  mutate(statistic="mean") %>% 
  rename(`il-6`=Estimate)

brm_il_6_md_dat<-brm_il_6_md_dat0 %>% 
#  as_tibble() %>% mutate(bmi=seq(20,45,length=seq_ln)) %>%
  as_tibble() %>% mutate(bmi_c=seq(-10,15,length=seq_ln)) %>%
  mutate(statistic="median") %>%
  rename(`il-6`=Estimate)

brm_il_6_q90_dat<-brm_il_6_q90_dat0 %>% 
 # as_tibble() %>% mutate(bmi=seq(20,45,length=seq_ln)) %>%
  as_tibble() %>% mutate(bmi_c=seq(-10,15,length=seq_ln)) %>%
  mutate(statistic="q90") %>%
  rename(`il-6`=Estimate)

brm_il_6_plt_dat <- rbind(brm_il_6_mn_dat, brm_il_6_md_dat, brm_il_6_q90_dat) 

ggplot(brm_il_6_plt_dat) + geom_line(aes(x=bmi_c,y=`il-6`,color=statistic)) +
geom_ribbon(aes(x=bmi_c,ymin=Q2.5,ymax=Q97.5,fill=statistic),alpha=0.2)+ylab("IL-6")+xlab("BMI (centered)")


#marginal_effects(brm_fit_il_6, "bmi")

brm_il_6_mn_dat %>% as_tibble() %>% mutate(bmi_c=seq(-10,15,length=seq_ln)) %>%
  ggplot() + geom_line(aes(x=bmi_c,y=`il-6`)) +
  geom_ribbon(aes(x=bmi_c,ymin=Q2.5,ymax=Q97.5),alpha=0.2)


plot(brm_fit_il_6)

plot(brm_fit_il_6_q50)

# posterior samples
brm_post<-as.data.frame(brm_fit_il_6)

```


## HIV example using Bayesian CPM model

### IL-6

```{r il_6_fit}
#!! find better inits?
fit_il_6_probit <- bayes_cpm(il_6 ~ bmi_c + age_c + male + 
                               nonwhite + study.linc + cd4_100_c + smoker,
                             data=hivdat, link="probit", seed=77845,
                             iter=4000, warmup=2000, chains=2, init="0", 
                             control = list(adapt_delta = 0.95))

fit_il_6_probit2 <- bayes_cpm(il_6 ~ bmi_c + age_c + male + 
                               nonwhite + study.linc + cd4_100_c + smoker,
                             data=hivdat, link="probit",
                             dir_prior_conc = function(n) 1/(0.8 + 0.35*n),
                             seed=77845,
                             iter=4000, warmup=2000, chains=2, init="0", 
                             control = list(adapt_delta = 0.95))

fit_il_6_logit <- bayes_cpm(il_6 ~ bmi_c + age_c + male + 
                              nonwhite + study.linc + cd4_100_c + smoker,
                            data=hivdat, link="logistic", seed=77845,
                            iter=4000, warmup=2000, chains=2, init="0", 
                            control = list(adapt_delta = 0.95))

fit_il_6_logit2 <- bayes_cpm(il_6 ~ bmi_c + age_c + male + 
                              nonwhite + study.linc + cd4_100_c + smoker,
                            data=hivdat, link="logistic",
                            dir_prior_conc = function(n) 1/(0.8 + 0.35*n),
                            seed=77845,
                            iter=4000, warmup=2000, chains=2, init="0", 
                            control = list(adapt_delta = 0.95))

fit_il_6_loglog <- bayes_cpm(il_6 ~ bmi_c + age_c + male + 
                              nonwhite + study.linc + cd4_100_c + smoker,
                            data=hivdat, link="loglog", seed=77845,
                            iter=4000, warmup=2000, chains=2, init="0", 
                            control = list(adapt_delta = 0.95))

fit_il_6_loglog2 <- bayes_cpm(il_6 ~ bmi_c + age_c + male + 
                              nonwhite + study.linc + cd4_100_c + smoker,
                            data=hivdat, link="loglog",
                            dir_prior_conc = function(n) 1/(0.8 + 0.35*n),
                            seed=77845,
                            iter=4000, warmup=2000, chains=2, init="0", 
                            control = list(adapt_delta = 0.95))

if(0){
fit_il_6_nocnt <- bayes_cpm(il_6 ~ bmi + age + male + nonwhite + study.linc +
                        cd4_100 + smoker, data=hivdat, link="probit", 
                        seed=77845, iter=4000, warmup=2000, chains=2, init="0", 
                        control = list(adapt_delta = 0.95))
}
```

```{r il_6_mcmc_chk}
qplot(summary(fit_il_6_probit$stanfit)$summary[,"Rhat"])
qplot(summary(fit_il_6_probit2$stanfit)$summary[,"Rhat"])
qplot(summary(fit_il_6_logit$stanfit)$summary[,"Rhat"])
qplot(summary(fit_il_6_logit2$stanfit)$summary[,"Rhat"])
qplot(summary(fit_il_6_loglog$stanfit)$summary[,"Rhat"])
qplot(summary(fit_il_6_loglog2$stanfit)$summary[,"Rhat"])

main_pars <- c("b[1]","b[2]","b[3]","b[4]","b[5]","b[6]","b[7]")
main_pars_labs <-c("bmi_c", "age_c", "male", "nonwhite", "study.linc", "cd4_100_c", "smoker")
#main_pars_labs2 <-c("BMI", "Age", "Male", "Nonwhite", "Study.linc", "CD4/100", "Smoker")
main_pars_labs2 <-c("BMI~(per~kg/m^2)", "Age~(per~1~year)", "Male~(vs.~Female)", "Non-white~(vs.~White)", "LiNC~study~(vs.~AIAC~study)", "CD4~(per~100~cells/mm^3)", "Smoker~(vs.~non-smoker)")


summary(fit_il_6_probit$stanfit, pars=main_pars)$summary %>% round(3) %>% set_rownames(main_pars_labs)

summary(fit_il_6_logit$stanfit, pars=main_pars)$summary %>% round(3) %>% set_rownames(main_pars_labs)

rstan::traceplot(fit_il_6_probit$stanfit, pars=main_pars) 
rstan::traceplot(fit_il_6_probit2$stanfit, pars=main_pars) 
rstan::traceplot(fit_il_6_logit$stanfit, pars=main_pars) 
rstan::traceplot(fit_il_6_logit2$stanfit, pars=main_pars) 
rstan::traceplot(fit_il_6_loglog$stanfit, pars=main_pars) 
rstan::traceplot(fit_il_6_loglog2$stanfit, pars=main_pars) 
# https://stackoverflow.com/questions/57411900/changing-the-title-of-individual-plots-within-a-panel-of-traceplots-in-ggplot-st
```


```{r il_6_loo}
# use updated model w/ output compatible with loo & ppc
#https://cloud.r-project.org/web/packages/loo/vignettes/loo2-example.html

ll_il6_probit<-extract_log_lik(fit_il_6_probit$stanfit)
ll_il6_probit2<-extract_log_lik(fit_il_6_probit2$stanfit)
ll_il6_logit<-extract_log_lik(fit_il_6_logit$stanfit)
ll_il6_logit2<-extract_log_lik(fit_il_6_logit2$stanfit)
ll_il6_loglog<-extract_log_lik(fit_il_6_loglog$stanfit)
ll_il6_loglog2<-extract_log_lik(fit_il_6_loglog2$stanfit)

# loglog w/ 1/(0.8 + 0.35*J) is preferred model
loo_comp1<-loo_compare(loo(ll_il6_probit), loo(ll_il6_probit2), loo(ll_il6_logit), loo(ll_il6_logit2),
 loo(ll_il6_loglog), loo(ll_il6_loglog2))

plot(loo(ll_il6_logit))
plot(loo(ll_il6_probit))
```


```{r il_6_pars}
# il_6_posta <- as.matrix(fit_il_6_loglog2$stanfit, pars=main_pars[c(1,2,6)])  %>%
# mcmc_intervals(prob_outer = 0.95) + scale_y_discrete(labels =main_pars_labs2[c(1,2,6)] )
# 
# il_6_postb <- as.matrix(fit_il_6_loglog2$stanfit, pars=main_pars[c(3,4,5,7)])  %>%
# mcmc_intervals(prob_outer = 0.95) + scale_y_discrete(labels = main_pars_labs2[c(3,4,5,7)])
# 
# plot_grid(il_6_posta,il_6_postb,ncol=1)
# 
# ggsave(file.path(figdir,"il_6_post.png"),width=6,height=4)


il_6_posta <- as.matrix(fit_il_6_loglog2$stanfit, pars=main_pars[c(1,2,6)]) %>%
mcmc_intervals(prob_outer = 0.95) + scale_y_discrete(labels =parse(text=main_pars_labs2[c(1,2,6)]))

il_6_postb <- as.matrix(fit_il_6_loglog2$stanfit, pars=main_pars[c(3,4,5,7)])  %>%
mcmc_intervals(prob_outer = 0.95) + scale_y_discrete(labels = parse(text=main_pars_labs2[c(3,4,5,7)]))

pg_il_6_a<-plot_grid(il_6_posta,il_6_postb,ncol=1)

pg_il_6_a

ggsave(file.path(figdir,"il_6_post_v2.png"),width=6,height=4)
ggsave(file.path(figdir,"il_6_post_v2.tiff"),width=6,height=4,bg='white',dpi=600, compression = "lzw")

# stan_plot(fit_il_6_loglog2$stanfit, pars=main_pars) + 
#   scale_y_discrete(breaks=paste0("b[",1:7,"]"),
#         labels=letters[1:7])

# cutpoints
il_6_ct_plt<-plot(fit_il_6_loglog2$stanfit,pars=c(paste0("cutpoints[",1:173,"]")),
                  ci_level = 0.5)

# il_6_ct_plt$data %>% mutate(y_rev=rev(y)) %>% 
#   ggplot(aes(x=y_rev,y=m)) + geom_point()

pg_il_6_b <- il_6_ct_plt$data %>% mutate(y_rev=rev(y)) %>% 
  ggplot(aes(x=y_rev,y=m)) + geom_line(aes(color='median')) +
  geom_ribbon(aes(ymin=l,ymax=h,fill='50% int'))+
  geom_ribbon(aes(ymin=ll,ymax=hh,fill='95% int'))+
  scale_color_manual(labels = c('median'),
                     values = c('blue')) +
  scale_fill_manual(labels= c('50% int','95% int'),
                    values = alpha('blue',c(0.4,0.2))) +
  labs(color='', fill='') +
  ylab(expression(hat(gamma)[j]))+xlab("j")

pg_il_6_b

ggsave(file.path(figdir,"il_6_post_ct.png"),width=6,height=4)
ggsave(file.path(figdir,"il_6_post_ct.tiff"),width=6,height=4,bg='white',dpi=600, compression = "lzw")

plot_grid(pg_il_6_a,pg_il_6_b,nrow=1)
ggsave(file.path(figdir,"il_6_post_est_ct.tiff"),width=8,height=4,bg='white',dpi=600, compression = "lzw")

# +
#   theme(axis.title.x=element_text(size=fctsiz),
#         axis.title.y = element_text(size=fctsiz),
#         axis.text =  element_text(size=atxtsz),
#         strip.text = element_text(size=fctsiz),
#         legend.text = element_text(size=fctsiz))

if(0){
traceplot(fit_il_6_loglog2$stanfit, pars=c(paste0("cutpoints[",1:20,"]")))
traceplot(fit_il_6_loglog2$stanfit, pars=c(paste0("cutpoints[",21:40,"]")))
traceplot(fit_il_6_loglog2$stanfit, pars=c(paste0("cutpoints[",41:60,"]")))
traceplot(fit_il_6_loglog2$stanfit, pars=c(paste0("cutpoints[",61:80,"]")))
traceplot(fit_il_6_loglog2$stanfit, pars=c(paste0("cutpoints[",81:100,"]")))
traceplot(fit_il_6_loglog2$stanfit, pars=c(paste0("cutpoints[",101:120,"]")))
traceplot(fit_il_6_loglog2$stanfit, pars=c(paste0("cutpoints[",121:140,"]")))
traceplot(fit_il_6_loglog2$stanfit, pars=c(paste0("cutpoints[",141:160,"]")))
traceplot(fit_il_6_loglog2$stanfit, pars=c(paste0("cutpoints[",161:173,"]")))
}

# cell probs
if(0){
plot(fit_il_6_loglog2$stanfit,pars=c(paste0("pi[",1:20,"]")))
plot(fit_il_6_loglog2$stanfit,pars=c(paste0("pi[",21:40,"]")))
plot(fit_il_6_loglog2$stanfit,pars=c(paste0("pi[",41:60,"]")))
plot(fit_il_6_loglog2$stanfit,pars=c(paste0("pi[",61:80,"]")))
plot(fit_il_6_loglog2$stanfit,pars=c(paste0("pi[",81:100,"]")))
plot(fit_il_6_loglog2$stanfit,pars=c(paste0("pi[",101:120,"]")))
plot(fit_il_6_loglog2$stanfit,pars=c(paste0("pi[",121:140,"]")))
plot(fit_il_6_loglog2$stanfit,pars=c(paste0("pi[",141:160,"]")))
plot(fit_il_6_loglog2$stanfit,pars=c(paste0("pi[",161:174,"]")))
}

```

```{r il_6_postpred}
dr <- 750
il_6_ppd <- getPostPred(fit_il_6_loglog2,newdata=hivdat[,main_pars_labs],draws=dr)

il_6_ppd_df <- il_6_ppd %>% as_tibble() %>% mutate(samp=1:n()) %>% 
  rename_with(~ gsub("V", "", .x,)) %>% pivot_longer(1:dr,names_to="rep")

il_6_dat_df<-hivdat0 %>% select(il_6) %>% mutate(rep=999,samp=1:n()) %>% rename(value=il_6)

y_reps_il_6 <- rbind(il_6_ppd_df,il_6_dat_df) %>% mutate(rep=as.numeric(rep))

ndr_plt <- 10
y_reps_il_6 %>% filter(rep<=ndr_plt|rep>998) %>% ggplot(aes(x=value, group=rep, col = if_else(rep>998,'y','yrep'))) + geom_density() + 
scale_color_manual(name = '', values = setNames(c('blue','lightblue'), c('y', 'yrep'))) + ylab('') + xlab('IL-6')

ggsave(file.path(figdir,"il_6_postpred.png"),width=6,height=4)
ggsave(file.path(figdir,"il_6_postpred.tiff"),width=6,height=4,dpi=600, compression = "lzw")
```


```{r il_6_ppp}
# posterior predictive p-values (aka Bayesian p-values)

# variance
ppp_var_il_6 <- y_reps_il_6 %>% group_by(rep) %>%
  summarize(ppp_v=var(value),.groups="drop") %>% 
  mutate(true_v=ifelse(rep>998,ppp_v,NA),
         true_v=max(true_v,na.rm=TRUE),
         pp_ind=ppp_v>true_v) %>% 
  filter(rep<998) %>% pull(pp_ind) %>% mean()
  
# skewness
library(e1071)  
ppp_skew_il_6 <- y_reps_il_6 %>% group_by(rep) %>%
  summarize(ppp_v=skewness(value),.groups="drop") %>% 
  mutate(true_v=ifelse(rep>998,ppp_v,NA),
         true_v=max(true_v,na.rm=TRUE),
         pp_ind=ppp_v>true_v) %>% 
  filter(rep<998) %>% pull(pp_ind) %>% mean()

# proportion censored
ppp_cens_il_6 <- y_reps_il_6 %>% group_by(rep) %>%
  summarize(ppp_v=mean(value==0),.groups="drop") %>% 
  mutate(true_v=ifelse(rep>998,ppp_v,NA),
         true_v=max(true_v,na.rm=TRUE),
         pp_ind=ppp_v>true_v) %>% 
  filter(rep<998) %>% pull(pp_ind) %>% mean()

c(ppp_var_il_6,ppp_skew_il_6,ppp_cens_il_6)
```



<!-- estimated transformation -->
```{r il_6_trans}
cp_il_6 <- summary(fit_il_6_loglog2$stanfit, pars=c(paste0("cutpoints[",1:173,"]")))$summary[,"50%"]

# cp_il_6_nocnt <- summary(fit_il_6_nocnt$stanfit, pars=c(paste0("cutpoints[",1:173,"]")))$summary[,"50%"]

il_6_trans_dat <- cbind(fit_il_6_loglog2$standata$truey0[2:174], cp_il_6) %>% as.data.frame() %>% rename(`IL-6`=V1)

# il_6_trans_dat_nocnt <- cbind(fit_il_6_nocnt$standata$truey0[2:174], cp_il_6_nocnt) %>% as.data.frame() %>% rename(`IL-6`=V1)

qplot(hivdat0$il_6,bins=35)+xlab("IL-6") + 
  annotate(geom="text", x=18, y=32, 
           label="3% below detection limit", color="red")
#ggsave(file.path(figdir,"il_6_hist.png"),width=6,height=4)

ggplot(data=il_6_trans_dat, aes(x=`IL-6`,y=cp_il_6, color="est")) + geom_step() + ylab(expression(gamma[j])) + xlim(c(1,10)) + 
  stat_function(fun=log, mapping = aes(color = "log")) + 
  stat_function(fun=sqrt, mapping = aes(color = "sqrt")) +
  scale_color_manual(name = "transformation",
                       values = c("black", "blue", "red"), # Color specification
                       labels = c("Bayes CPM estimate","log", "sqrt"))
  
# take log transform and then is distributed as Gumbel ?
# take sqrt transform and then is distributed as Gumbel ?

ggsave(file.path(figdir,"il_6_trans.png"),width=6,height=4)


ggplot(data=il_6_trans_dat, aes(x=`IL-6`,y=cp_il_6)) + geom_step() + ylab(expression(gamma[j])) + xlim(c(1,15)) + ylim(c(-2.5,4)) #+ 
 # scale_color_manual(name = "transformation",
 #                      values = c("black", "blue", "red"), # Color specification
 #                      labels = c("Bayes CPM estimate","log", "sqrt"))
  
# take log transform and then is distributed as Gumbel ?
# take sqrt transform and then is distributed as Gumbel ?

ggsave(file.path(figdir,"il_6_trans_v2.png"),width=6,height=4)
ggsave(file.path(figdir,"il_6_trans_v2.tiff"),width=6,height=4,dpi=600, compression = "lzw")

# ggplot(data=il_6_trans_dat_nocnt, aes(x=`IL-6`,y=cp_il_6_nocnt, color="est")) + geom_step() + ylab(expression(gamma[j])) + xlim(c(1,10)) + 
#   stat_function(fun=log, mapping = aes(color = "log")) + 
#   stat_function(fun=sqrt, mapping = aes(color = "sqrt")) +
#   scale_color_manual(name = "transformation",
#                        values = c("black", "blue", "red"), # Color specification
#                        labels = c("Bayes CPM estimate","log", "sqrt"))
# 
# ggsave(file.path(figdir,"il_6_trans_nocnt.png"),width=6,height=4)

```

```{r il_6_bmi_dat}
mn_summ_il_6 <- getMean(fit_il_6_loglog2, newdat1)
#mn_vals_il_6 <- getMean(fit_il_6_loglog2, newdat1, summ=FALSE)
med_summ_il_6 <- getQuantile(fit_il_6_loglog2, newdat1, q=0.5)
q90_summ_il_6 <- getQuantile(fit_il_6_loglog2, newdat1, q=0.9)
```

```{r il_6_bmi_plt}
il_6_mn_dat<-mn_summ_il_6 %>% 
  select(med_mn, mn_q2.5, mn_q97.5, bmi_c) %>% 
  mutate(statistic="mean") %>% 
  rename(`il-6`=med_mn,q2.5=mn_q2.5,q97.5=mn_q97.5)

il_6_med_dat<-med_summ_il_6 %>% 
  select(med_qtile, qtile_q2.5, qtile_q97.5, bmi_c) %>%
  mutate(statistic="median") %>%
  rename(`il-6`=med_qtile,q2.5=qtile_q2.5,q97.5=qtile_q97.5)

il_6_q90_dat<-q90_summ_il_6 %>% 
  select(med_qtile, qtile_q2.5, qtile_q97.5, bmi_c) %>%
  mutate(statistic="q90") %>%
  rename(`il-6`=med_qtile,q2.5=qtile_q2.5,q97.5=qtile_q97.5)

il_6_plt_dat <- rbind(il_6_mn_dat,il_6_med_dat,il_6_q90_dat) 

ggplot(il_6_plt_dat) + geom_line(aes(x=bmi_c,y=`il-6`,color=statistic)) +
geom_ribbon(aes(x=bmi_c,ymin=q2.5,ymax=q97.5,fill=statistic),alpha=0.2)+ylab("IL-6")+xlab("difference from mean BMI")

ggsave(file.path(figdir,"il_6_bmi.png"),width=6,height=3.5)
ggsave(file.path(figdir,"il_6_bmi.tiff"),width=6,height=3.5,dpi=600, compression = "lzw")
```

```{r, eval=FALSE}

if (0){
mn_summ_brm_il_6_0 <- fitted(brm_fit_il_6, newdata=newdat1)
mn_summ_brm_il_6<-cbind(newdat1,mn_summ_brm_il_6_0)

mn_summ_il_6 %>% 
  ggplot(aes(x=bmi,y=med_mn)) + geom_line() +
  geom_ribbon(aes(ymin=mn_q2.5,ymax=mn_q97.5),alpha=0.2) +
  geom_line(data=mn_summ_brm_il_6,aes(x=bmi,y=Estimate^2,color="red"))+
  geom_ribbon(data=mn_summ_brm_il_6,aes(x=bmi,ymin=Q2.5^2,ymax=Q97.5^2,fill="red"),alpha=0.2,inherit.aes = FALSE)
}

mn_summ_il_6 %>% 
  ggplot(aes(x=bmi,y=med_mn)) + geom_line() +
  geom_ribbon(aes(ymin=mn_q5,ymax=mn_q95),alpha=0.2)

med_summ_il_6 %>% 
  ggplot(aes(x=bmi,y=med_qtile)) + geom_line() +
  geom_ribbon(aes(ymin=qtile_q5,ymax=qtile_q95),alpha=0.2)
```


### IL-1-beta 


```{r il_1_beta_fit}
fit_il_1_beta_probit <- bayes_cpm(il_1_beta ~ bmi_c + age_c + male + 
                               nonwhite + study.linc + cd4_100_c + smoker,
                             data=hivdat, link="probit", seed=3845,
                             iter=4000, warmup=2000, chains=2, init="0", 
                             control = list(adapt_delta = 0.95))

fit_il_1_beta_probit2 <- bayes_cpm(il_1_beta ~ bmi_c + age_c + male + 
                               nonwhite + study.linc + cd4_100_c + smoker,
                             data=hivdat, link="probit",
                             dir_prior_conc = function(n) 1/(0.8 + 0.35*n),
                             seed=3845,
                             iter=4000, warmup=2000, chains=2, init="0", 
                             control = list(adapt_delta = 0.95))

fit_il_1_beta_logit <- bayes_cpm(il_1_beta ~ bmi_c + age_c + male + 
                              nonwhite + study.linc + cd4_100_c + smoker,
                            data=hivdat, link="logistic", seed=6734,
                            iter=4000, warmup=2000, chains=2, init="0", 
                            control = list(adapt_delta = 0.95))


fit_il_1_beta_logit2 <- bayes_cpm(il_1_beta ~ bmi_c + age_c + male + 
                              nonwhite + study.linc + cd4_100_c + smoker,
                            data=hivdat, link="logistic", 
                            dir_prior_conc = function(n) 1/(0.8 + 0.35*n),
                            seed=6734,
                            iter=4000, warmup=2000, chains=2, init="0", 
                            control = list(adapt_delta = 0.95))


fit_il_1_beta_loglog <- bayes_cpm(il_1_beta ~ bmi_c + age_c + male + 
                              nonwhite + study.linc + cd4_100_c + smoker,
                            data=hivdat, link="loglog", seed=6734,
                            iter=4000, warmup=2000, chains=2, init="0", 
                            control = list(adapt_delta = 0.95))


fit_il_1_beta_loglog2 <- bayes_cpm(il_1_beta ~ bmi_c + age_c + male + 
                              nonwhite + study.linc + cd4_100_c + smoker,
                            data=hivdat, link="loglog", 
                            dir_prior_conc = function(n) 1/(0.8 + 0.35*n),
                            seed=6734,
                            iter=4000, warmup=2000, chains=2, init="0", 
                            control = list(adapt_delta = 0.95))


```


```{r il_1_beta_mcmc_chk}
qplot(summary(fit_il_1_beta_probit$stanfit)$summary[,"Rhat"])
qplot(summary(fit_il_1_beta_probit2$stanfit)$summary[,"Rhat"])
qplot(summary(fit_il_1_beta_logit$stanfit)$summary[,"Rhat"])
qplot(summary(fit_il_1_beta_logit2$stanfit)$summary[,"Rhat"])
qplot(summary(fit_il_1_beta_loglog$stanfit)$summary[,"Rhat"])
qplot(summary(fit_il_1_beta_loglog2$stanfit)$summary[,"Rhat"])

main_pars <- c("b[1]","b[2]","b[3]","b[4]","b[5]","b[6]","b[7]")
main_pars_labs <-c("bmi_c", "age_c", "male", "nonwhite", "study.linc", "cd4_100_c", "smoker")
main_pars_labs2 <-c("BMI~(per~kg/m^2)", "Age~(per~1~year)", "Male~(vs.~Female)", "Non-white~(vs.~White)", "LiNC~study~(vs.~AIAC~study)", "CD4~(per~100~cells/mm^3)", "Smoker~(vs.~non-smoker)")

summary(fit_il_1_beta_probit$stanfit, pars=main_pars)$summary %>% round(3) %>% set_rownames(main_pars_labs)

summary(fit_il_1_beta_logit$stanfit, pars=main_pars)$summary %>% round(3) %>% set_rownames(main_pars_labs)

rstan::traceplot(fit_il_1_beta_probit$stanfit, pars=main_pars) 
rstan::traceplot(fit_il_1_beta_probit2$stanfit, pars=main_pars) 
rstan::traceplot(fit_il_1_beta_logit$stanfit, pars=main_pars)
rstan::traceplot(fit_il_1_beta_logit2$stanfit, pars=main_pars)
rstan::traceplot(fit_il_1_beta_loglog$stanfit, pars=main_pars)
rstan::traceplot(fit_il_1_beta_loglog2$stanfit, pars=main_pars)

```


```{r il_1_beta_loo}
ll_il_1_beta_logit<-extract_log_lik(fit_il_1_beta_logit$stanfit)
ll_il_1_beta_logit2<-extract_log_lik(fit_il_1_beta_logit2$stanfit)
ll_il_1_beta_probit<-extract_log_lik(fit_il_1_beta_probit$stanfit)
ll_il_1_beta_probit2<-extract_log_lik(fit_il_1_beta_probit2$stanfit)
ll_il_1_beta_loglog<-extract_log_lik(fit_il_1_beta_loglog$stanfit)
ll_il_1_beta_loglog2<-extract_log_lik(fit_il_1_beta_loglog2$stanfit)


# loglog2 is preferred model
loo_comp2<-loo_compare(loo(ll_il_1_beta_logit), loo(ll_il_1_beta_logit2), loo(ll_il_1_beta_probit),loo(ll_il_1_beta_probit2),
                       loo(ll_il_1_beta_loglog),loo(ll_il_1_beta_loglog2))

plot(loo(ll_il_1_beta_logit))
plot(loo(ll_il_1_beta_logit2))
plot(loo(ll_il_1_beta_probit))
plot(loo(ll_il_1_beta_loglog))

```


```{r il_1_beta_pars}
# il_1_beta_posta <- as.matrix(fit_il_1_beta_loglog2$stanfit, pars=main_pars[c(1,2,6)]) %>%
# mcmc_intervals(prob_outer = 0.95) + scale_y_discrete(labels = main_pars_labs2[c(1,2,6)])
# 
# il_1_beta_postb <- as.matrix(fit_il_1_beta_loglog2$stanfit, pars=main_pars[c(3,4,5,7)])  %>%
# mcmc_intervals(prob_outer = 0.95) + scale_y_discrete(labels = main_pars_labs2[c(3,4,5,7)])
# 
# plot_grid(il_1_beta_posta,il_1_beta_postb,ncol=1)
# 
# ggsave(file.path(figdir,"il_1_beta_post.png"),width=6,height=4)

il_1_beta_posta <- as.matrix(fit_il_1_beta_loglog2$stanfit, pars=main_pars[c(1,2,6)]) %>%
mcmc_intervals(prob_outer = 0.95) + scale_y_discrete(labels = parse(text=main_pars_labs2[c(1,2,6)]))

il_1_beta_postb <- as.matrix(fit_il_1_beta_loglog2$stanfit, pars=main_pars[c(3,4,5,7)])  %>%
mcmc_intervals(prob_outer = 0.95) + scale_y_discrete(labels = parse(text=main_pars_labs2[c(3,4,5,7)]))

pg_il_1_beta_a <- plot_grid(il_1_beta_posta,il_1_beta_postb,ncol=1)

pg_il_1_beta_a
ggsave(file.path(figdir,"il_1_beta_post_v2.png"),width=6,height=4)
ggsave(file.path(figdir,"il_1_beta_post_v2.tiff"),width=6,height=4,bg='white',dpi=600, compression = "lzw")

# stan_plot(fit_il_1_beta_loglog2$stanfit, pars=main_pars) + 
#   scale_y_discrete(breaks=paste0("b[",1:7,"]"),
#         labels=letters[1:7])

# cutpoints
il_1_beta_ct_plt<-plot(fit_il_1_beta_loglog2$stanfit,pars=c(paste0("cutpoints[",1:113,"]")),
                  ci_level = 0.5)

# il_1_beta_ct_plt$data %>% mutate(y_rev=rev(y)) %>% 
#   ggplot(aes(x=y_rev,y=m)) + geom_point()

pg_il_1_beta_b <- il_1_beta_ct_plt$data %>% mutate(y_rev=rev(y)) %>% 
  ggplot(aes(x=y_rev,y=m)) + geom_line(aes(color='median')) +
  geom_ribbon(aes(ymin=l,ymax=h,fill='50% int'))+
  geom_ribbon(aes(ymin=ll,ymax=hh,fill='95% int'))+
  scale_color_manual(labels = c('median'),
                     values = c('blue')) +
  scale_fill_manual(labels= c('50% int','95% int'),
                    values = alpha('blue',c(0.4,0.2))) +
  labs(color='', fill='') +
  ylab(expression(hat(gamma)[j]))+xlab("j")

pg_il_1_beta_b
ggsave(file.path(figdir,"il_1_beta_post_ct.png"),width=6,height=4)
ggsave(file.path(figdir,"il_1_beta_post_ct.tiff"),width=6,height=4,bg='white',dpi=600, compression = "lzw")

plot_grid(pg_il_1_beta_a,pg_il_1_beta_b,nrow=1)
ggsave(file.path(figdir,"il_1_beta_post_est_ct.tiff"),width=8,height=4,bg='white',dpi=600, compression = "lzw")


if(0){
traceplot(fit_il_1_beta_loglog2$stanfit, pars=c(paste0("cutpoints[",1:20,"]")))
traceplot(fit_il_1_beta_loglog2$stanfit, pars=c(paste0("cutpoints[",21:40,"]")))
traceplot(fit_il_1_beta_loglog2$stanfit, pars=c(paste0("cutpoints[",41:60,"]")))
traceplot(fit_il_1_beta_loglog2$stanfit, pars=c(paste0("cutpoints[",61:80,"]")))
traceplot(fit_il_1_beta_loglog2$stanfit, pars=c(paste0("cutpoints[",81:100,"]")))
traceplot(fit_il_1_beta_loglog2$stanfit, pars=c(paste0("cutpoints[",101:113,"]")))
}

# cell probs
if(0){
plot(fit_il_1_beta_loglog2$stanfit,pars=c(paste0("pi[",1:20,"]")))
plot(fit_il_1_beta_loglog2$stanfit,pars=c(paste0("pi[",21:40,"]")))
plot(fit_il_1_beta_loglog2$stanfit,pars=c(paste0("pi[",41:60,"]")))
plot(fit_il_1_beta_loglog2$stanfit,pars=c(paste0("pi[",61:80,"]")))
plot(fit_il_1_beta_loglog2$stanfit,pars=c(paste0("pi[",81:100,"]")))
plot(fit_il_1_beta_loglog2$stanfit,pars=c(paste0("pi[",101:114,"]")))
}


# summary(fit_il_1_beta, pars=main_pars)$summary %>% round(3) %>% set_rownames(main_pars_labs)
# plot(fit_il_1_beta, pars=main_pars)
# 
# as.matrix(fit_il_1_beta, pars=main_pars)  %>%
# mcmc_intervals() + scale_y_discrete(labels = c("bmi_c", "age_c", "male", "nonwhite",
#                                  "study.linc", "cd4_100_c", "smoker"))
# 
# fit_il_1_betaa<-as.matrix(fit_il_1_beta, pars=main_pars[c(1,2,6)])  %>%
# mcmc_intervals() + scale_y_discrete(labels =main_pars_labs2[c(1,2,6)] )
# 
# fit_il_1_betab<-as.matrix(fit_il_1_beta, pars=main_pars[c(3,4,5,7)])  %>%
# mcmc_intervals() + scale_y_discrete(labels = main_pars_labs2[c(3,4,5,7)])
# plot_grid(fit_il_1_betaa,fit_il_1_betab,ncol=1)
# 
# ggsave(file.path(figdir,"il_1_beta_post.png"),width=6,height=4)
```

```{r il_1_beta_postpred}
dr <- 750
il_1_beta_ppd <- getPostPred(fit_il_1_beta_loglog2,newdata=hivdat[,main_pars_labs],draws=dr)

il_1_beta_ppd_df <- il_1_beta_ppd %>% as_tibble() %>% mutate(samp=1:n()) %>% 
  rename_with(~ gsub("V", "", .x,)) %>% pivot_longer(1:dr,names_to="rep")

il_1_beta_dat_df<-hivdat0 %>% select(il_1_beta) %>% mutate(rep=999,samp=1:n()) %>% rename(value=il_1_beta)

y_reps_il_1_beta <- rbind(il_1_beta_ppd_df,il_1_beta_dat_df) %>% mutate(rep=as.numeric(rep))

ndr_plt <- 10
y_reps_il_1_beta %>% filter(rep<=ndr_plt|rep>998) %>% ggplot(aes(x=value, group=rep, col = if_else(rep>998,'y','yrep'))) + geom_density() +
scale_color_manual(name = '', values = setNames(c('blue','lightblue'), c('y', 'yrep'))) + ylab('') + xlab(expression("IL-1-"*beta))

ggsave(file.path(figdir,"il_1_beta_postpred.png"),width=6,height=4)
ggsave(file.path(figdir,"il_1_beta_postpred.tiff"),width=6,height=4,dpi=600, compression = "lzw")
```

```{r}
# see https://cloud.r-project.org/web/packages/loo/vignettes/loo2-example.html
il_1_beta_ppd2 <- getPostPred(fit_il_1_beta_loglog2,newdata=hivdat[,main_pars_labs],draws=4000)

loo_il_1_beta_loglog2 <- loo(fit_il_1_beta_loglog2, save_psis = TRUE)
ppc_loo_pit_overlay(
  y = as.numeric(as.character(hivdat$il_1_beta)),
  yrep = t(il_1_beta_ppd2),
  lw = weights(loo_il_1_beta_probit2$psis_object)
)
```

```{r il_1_beta_ppp}
# posterior predictive p-values (aka Bayesian p-values)

# variance
ppp_var_il_1_beta <- y_reps_il_1_beta %>% group_by(rep) %>%
  summarize(ppp_v=var(value),.groups="drop") %>% 
  mutate(true_v=ifelse(rep>998,ppp_v,NA),
         true_v=max(true_v,na.rm=TRUE),
         pp_ind=ppp_v>true_v) %>% 
  filter(rep<998) %>% pull(pp_ind) %>% mean()
  
# skewness
library(e1071)  
ppp_skew_il_1_beta <- y_reps_il_1_beta %>% group_by(rep) %>%
  summarize(ppp_v=skewness(value),.groups="drop") %>% 
  mutate(true_v=ifelse(rep>998,ppp_v,NA),
         true_v=max(true_v,na.rm=TRUE),
         pp_ind=ppp_v>true_v) %>% 
  filter(rep<998) %>% pull(pp_ind) %>% mean()

# proportion censored
ppp_cens_il_1_beta <- y_reps_il_1_beta %>% group_by(rep) %>%
  summarize(ppp_v=mean(value==0),.groups="drop") %>% 
  mutate(true_v=ifelse(rep>998,ppp_v,NA),
         true_v=max(true_v,na.rm=TRUE),
         pp_ind=ppp_v>true_v) %>% 
  filter(rep<998) %>% pull(pp_ind) %>% mean()

c(ppp_var_il_1_beta, ppp_skew_il_1_beta, ppp_cens_il_1_beta)
```

```{r il_1_beta_trans}
cp_il_1_beta <- summary(fit_il_1_beta_loglog2$stanfit, pars=c(paste0("cutpoints[",1:113,"]")))$summary[,"mean"]

il_1_beta_trans_dat <- cbind(fit_il_1_beta_loglog2$standata$truey0[2:114], cp_il_1_beta) %>% as.data.frame() %>%
  rename(`IL-1-beta`=V1)

qplot(hivdat0$il_1_beta,bins=35)+xlab(expression(IL-1-beta)) + 
  annotate(geom="text", x=19, y=98, label="39% below detection limit", color="red")
#ggsave(file.path(figdir,"il_1_beta_hist.png"),width=6,height=4)

ggplot(data=il_1_beta_trans_dat,aes(x=`IL-1-beta`,y=cp_il_1_beta,color = "est")) + geom_step() + ylab(expression(gamma[j])) + ylim(c(-0.5,2.5)) + xlim(c(0,3.0)) + xlab(expression(IL-1-beta)) +
  stat_function(fun=log, mapping = aes(color = "log")) + 
  stat_function(fun=sqrt, mapping = aes(color = "sqrt")) +
  scale_color_manual(name = "transformation",
                       values = c("black", "blue", "red"), # Color specification
                       labels = c("Bayes CPM estimate","log", "sqrt"))

ggsave(file.path(figdir,"il_1_beta_trans.png"),width=6,height=4)


ggplot(data=il_1_beta_trans_dat,aes(x=`IL-1-beta`,y=cp_il_1_beta)) + geom_step() + ylab(expression(gamma[j])) + ylim(c(-0.5,2.5)) + xlim(c(0,3.0)) + xlab(expression(IL-1-beta))

ggsave(file.path(figdir,"il_1_beta_trans_v2.png"),width=6,height=4)
ggsave(file.path(figdir,"il_1_beta_trans_v2.tiff"),width=6,height=4,dpi=600, compression = "lzw")
```

```{r il_1_beta_bmi_dat}
mn_summ_il_1_beta <- getMean(fit_il_1_beta_loglog2, newdat1)
med_summ_il_1_beta <- getQuantile(fit_il_1_beta_loglog2, newdat1, q=0.5)
q90_summ_il_1_beta <- getQuantile(fit_il_1_beta_loglog2, newdat1, q=0.9)
```

```{r il_1_beta_bmi_plt}
il_1_beta_mn_dat <- mn_summ_il_1_beta %>% 
  select(med_mn, mn_q2.5, mn_q97.5, bmi_c) %>% 
  mutate(statistic="mean") %>% 
  rename(`il-1-beta`=med_mn,q2.5=mn_q2.5,q97.5=mn_q97.5)

il_1_beta_med_dat <- med_summ_il_1_beta %>% 
  select(med_qtile, qtile_q2.5, qtile_q97.5, bmi_c) %>%
  mutate(statistic="median") %>%
  rename(`il-1-beta`=med_qtile,q2.5=qtile_q2.5,q97.5=qtile_q97.5)

il_1_beta_q90_dat <- q90_summ_il_1_beta %>% 
  select(med_qtile, qtile_q2.5, qtile_q97.5, bmi_c) %>%
  mutate(statistic="q90") %>%
  rename(`il-1-beta`=med_qtile,q2.5=qtile_q2.5,q97.5=qtile_q97.5)

il_1_beta_plt_dat <- rbind(il_1_beta_mn_dat,il_1_beta_med_dat,il_1_beta_q90_dat) 

ggplot(il_1_beta_plt_dat) + geom_line(aes(x=bmi_c,y=`il-1-beta`,color=statistic)) +
geom_ribbon(aes(x=bmi_c,ymin=q2.5,ymax=q97.5,fill=statistic),alpha=0.2)+
  xlab("difference from mean BMI")+ylab(expression(IL-1-beta))

ggsave(file.path(figdir,"il_1_beta_bmi.png"),width=6,height=3.5)
ggsave(file.path(figdir,"il_1_beta_bmi.tiff"),width=6,height=3.5,dpi=600, compression = "lzw")
```