simple_ex.rmd

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

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

#rm(list=ls())
libs <- c("rstan", "rms", "dplyr", "stringr", "readr", "pammtools", "bayesCPM")
invisible(lapply(libs, library, character.only = TRUE))

dir <- file.path("/home/nathan/Dropbox/njames/school/PhD/orise_ra/bayes_cpm")
# CPM functions
#source(file.path(dir,"cpm_functions.r"))

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

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

```{r}
# compile/read in CPM models
# concentration (alpha) is given as a scalar parameter along with in data
#ord_mod1 <- readRDS(file.path(dir,"ordinal_model_1.rds"))
```

<!-- similar to orm() example 4a -->

```{r ex4a.1, cache=TRUE}
# Compare predicted mean with ols for a continuous x
set.seed(1567)
n <- 100
x1 <- rnorm(n)
# y <- 0.9*x1 + rnorm(n) # old w/ normal error
y <- 0.9*x1 + rlogis(n)
dat3 <- data.frame(y=ordered(y),y_num=y,x1)
```

```{r ex4a.2, cache=TRUE}
# mod_data  <- mkStanDat(dat3, outcome="y", preds = c("x1"), link=2) # probit link
# mod_data  <- mkStanDat(dat3, outcome="y", preds = c("x1"), link=1)
# 
# bg <- bayes_cpm(mod_data, seed=6472, 
#                  iter=3250, warmup=2000, chains=4,
#                  control = list(adapt_delta = 0.8))
#bg <- sampling(ord_mod1, data=mod_data, seed=6472, 
#                 iter=3250, warmup=2000, chains=4,
#                 control = list(adapt_delta = 0.8))

bg <- bayes_cpm(y~x1, data=dat3, link="logistic", seed=6472, 
                 iter=3250, warmup=2000, chains=4,
                 control = list(adapt_delta = 0.8))
```


```{r, eval=FALSE}
#plots, etc for biostat seminar
bg_df <- as.data.frame(bg$stanfit)

#head(bg_df[,c("cutpoints[98]")])
#head(bg_df[,c("b[1]")])

head(bg_df[,c("cutpoints[1]","cutpoints[2]","cutpoints[3]","cutpoints[98]","cutpoints[99]","b[1]")])
```

```{r}
cdf_bg <- getCDF(bg, newdata=data.frame(x1=c(-2,0,2))) 

cdf_bg %>% filter(ndrow %in% c(1,2,3)) %>% ggplot(aes(group=ndrow)) +
  geom_stepribbon(aes(x=yval, ymin=cdf_q5, ymax=cdf_q95, 
                      fill=factor(ndrow)) , alpha=0.4) +
  geom_step(aes(x=yval, y=med_cdf,color=factor(ndrow))) +
  xlab("y") + ylab("Conditional CDF") + 
  scale_fill_discrete(name = "covariate value", 
                      labels=c("x = -2", "x = 0", "x = 2")) + 
  scale_color_discrete(name = "covariate value", 
                       labels=c("x = -2", "x = 0", "x = 2")) +
  stat_function(fun=plogis,color="darkgreen",
                linetype=2, alpha=0.4) + 
  stat_function(fun=function(x) plogis(x,2*0.9), color="blue",
                linetype=2, alpha=0.4) +
  stat_function(fun=function(x) plogis(x,-2*0.9),color="red",
                linetype=2, alpha=0.4)

ggsave(file.path(figdir,"cond_cdf2.png"),width=6,height=3)


mn_dat<-getMean(bg, newdata=data.frame(x1=c(-2,0,2)),summ=FALSE)

ggplot(mn_dat,aes(x=mn,fill=factor(x1)))+geom_density(alpha=0.6,color=NA)+ 
  scale_fill_discrete(name = "covariate value", 
                      labels=c("x = -2", "x = 0", "x = 2")) + 
  xlab("") + ylab("conditional mean density")

ggsave(file.path(figdir,"cond_mn2.png"),width=6,height=3)

mn_dat %>% filter(x1==2) %>% pull(mn) %>% quantile(probs=c(0.025,0.975))

q50_dat <- getQuantile(bg, newdata=data.frame(x1=c(-2,0,2)),q=0.50,summ=FALSE)

ggplot(q50_dat,aes(x=qtile,fill=factor(x1)))+
  geom_density(alpha=0.6, color=NA, adjust=3)+ 
  scale_fill_discrete(name = "covariate value", 
                      labels=c("x = -2", "x = 0", "x = 2")) + 
  xlab("") + ylab("conditional median density")

ggsave(file.path(figdir,"cond_md2.png"),width=6,height=3)

q50_x0_samps<-q50_dat %>% filter(x1==0) %>% pull(qtile)

mean(q50_x0_samps > 0.25 | q50_x0_samps < -0.25)

```