orm_ex.rmd

---
title: "orm() documentation examples using Bayesian CPM"
output:
  html_document:
    toc: no
    toc_depth: 3
    number_sections: false
    code_folding: hide
    theme: paper
---

<!-- 
file based on orm_ex_with_rstanarm.rmd in orise_ra > rstanarm folder 
modified to us more generic Stan code (still based on stan_polr)
in addition to stan_polr
-->

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

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

# CPM functions
#dir<-getwd()
#source(file.path(dir,"cpm_functions.r"))

# for _dt versions of cpm functions
# library(data.table)
# library(dtplyr)


set.seed(4991)

# 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
# if (0){
# ord_mod_file1 <- read_file(file.path(getwd(),"ordinal_model_1.stan"))
# ord_mod1 <- stan_model(model_code = ord_mod_file1)
# saveRDS(ord_mod1, file = file.path(getwd(),"ordinal_model_1.rds"))
# }
# 
# ord_mod1 <- readRDS(file.path(getwd(),"ordinal_model_1.rds"))
```


## orm() example 1

```{r ex1.1, cache=TRUE}
set.seed(1)
n <- 100
y <- round(runif(n), 2)
x1 <- sample(c(-1,0,1), n, TRUE)
x2 <- sample(c(-1,0,1), n, TRUE)
dat1 <- data.frame(y=ordered(y),x1,x2)
rm(n, y, x1, x2)

num_y <- length(levels(dat1$y))
num_par <- 2

# orm
g <- orm(y ~ x1 + x2, data=dat1, eps=1e-5)

# try equivalent mods using other functions/packages
f <- lrm(y ~ x1 + x2, data=dat1, eps=1e-5)
f2 <- clm(factor(y) ~ x1 + x2, data=dat1)
p <- polr(factor(y) ~ x1 + x2, data=dat1)
p_summ <- summary(p)

g
```

```{r ex1.2, cache=TRUE}
mod_data1 <- mkStanDat(dat1, outcome="y", preds = c("x1","x2"), link=1)

fit_ord1 <- bayes_cpm(mod_data1, data=mod_data1, seed=12345, 
                 iter=8000, warmup=2000, chains=4,
                 control = list(adapt_delta = 0.8))

summary(fit_ord1,pars=c("b[1]","b[2]"))$summary %>% round(3)

b_summ<-summary(fit_ord1)$summary

plot(fit_ord1,pars="pi")
plot(fit_ord1,pars=c(paste0("pi[",1:30,"]"))) # approx 1/67+1/67 ??
summary(fit_ord1,pars=c(paste0("pi[",1:30,"]")))$summary # approx 1/67+1/67 or 2/67+1/67??

(1+(1/67))/(67)

(2+(1/67))/(67)

#means
#0.009593060
#0.019069737
#0.029651789
mns<-summary(fit_ord1,pars=c(paste0("pi")))$summary[,'mean']
sum(mns)

cbind(mns*n,table(y))
cbind(table(y)/n,mns)

#use table(y)/n as inits for pi?

#medians
meds<-summary(fit_ord1,pars=c(paste0("pi")))$summary[,'50%']
#c(0.006747834, 0.016088113, 0.026757347)

# b_summ[(num_y+1):(num_y+num_par),1] # coefficients
# b_summ[(num_y+num_par+1):(2*num_y+num_par-1),1] # cutpoints

mod_data1_pis <- table(mod_data1$y)/mod_data1$N
mod_data1_cp_init <- qlogis(cumsum(mod_data1_pis))[1:66]

ini <- list(cutpoints=mod_data1_cp_init )
ini_list <- lapply(1:4, function(x) lapply(ini,jitter))


fit_ord1b <- bayes_cpm(ord_mod1, data=mod_data1, seed=12345, 
                 iter=8000, warmup=2000, chains=4, 
                 init = ini_list,
                 control = list(adapt_delta = 0.8))

fit_ord1b
```

<!--- modify below here to use bayesCPM --->

```{r ex1.3, cache=TRUE}
# compare intercepts from lrm, orm, ordinal (clm), polr, and bayes_cpm
ints<-cbind(coef(f)[1:(num_y-1)],
      coef(g)[1:(num_y-1)],
      coef(f2)[1:(num_y-1)],
      coef(p_summ)[(num_par+1):(num_y+1),1],
      b_summ[(num_y+num_par+1):(2*num_y+num_par-1),6])

# coefs for rstanarm are negative of coefs from lrm, orm 
qplot(coef(g)[1:(num_y-1)], -1*b_summ[(num_y+num_par+1):(2*num_y+num_par-1),6],
      xlim=c(-5.5,5.5), ylim=c(-5.5,5.5)) + 
  xlab("orm intercept coefs") + ylab("-1 * bayes cpm intercept coefs (median)") +
  geom_abline(slope=1,intercept=0)

# compare intercepts for all 5 mods
colnames(ints)<-c("lrm","orm","clm","polr","bayes_cpm")

ints_plt <- ints %>% as_tibble() %>% 
  mutate(row=1:nrow(ints),lrm=-lrm,orm=-orm) %>% 
  tidyr::gather(colnames(ints),key="mod",value="value") 

# all are identical except bayes_cpm in tails
ints_plt %>% ggplot(aes(x=row,y=value, color=mod)) + geom_point()
```

```{r ex1.4, cache=TRUE}
# covariate coefficients
coef(f)[67:68] # lrm
coef(g)[67:68] # orm
coef(f2)[67:68] # ordinal
coef(p) # polr
b_summ[(num_y+1):(num_y+num_par),6] # bayes cpm

# difference between NPMLE and Bayes CPM
sum(abs(coef(p)-b_summ[(num_y+1):(num_y+num_par),6]))

rm(f,g,f2,p,p_summ,b_summ, fit_ord1, mod_data1, ints, ints_plt, dat1, num_par, num_y)
```

## orm() example 2a

```{r ex2a.1}
set.seed(1)
n <- 300
x1 <- c(rep(0,150), rep(1,150))
y <- rnorm(n) + 3*x1
dat2<-data.frame(y=ordered(y),x1)

g <- orm(y ~ x1, data=dat2)

k <- coef(g)
i <- num.intercepts(g)
h <- orm(y ~ x1, family=probit, data=dat2)
ll <- orm(y ~ x1, family=loglog, data=dat2)
cll <- orm(y ~ x1, family=cloglog, data=dat2)
cau <- orm(y ~ x1, family=cauchit, data=dat2)
x <- 1:i
z <- list(logistic=list(x=x, y=coef(g)[1:i]),
          probit  =list(x=x, y=coef(h)[1:i]),
          loglog  =list(x=x, y=coef(ll)[1:i]),
          cloglog =list(x=x, y=coef(cll)[1:i]))
labcurve(z, pl=TRUE, col=1:4, ylab='Intercept')
```

```{r ex2a.2, cache=TRUE}
mod_data3cau <- mod_data3cll <- mod_data3ll <- mod_data3h <- mod_data3 <- mkStanDat(dat2, outcome="y", preds = c("x1"), link=1)

mod_data3h$link <- 2
mod_data3ll$link <- 3 
mod_data3cll$link <- 4
mod_data3cau$link <- 5 
```

```{r ex2a.3, cache=TRUE}
#probit
bg <- sampling(ord_mod1, data=mod_data3, seed=6472, 
                 iter=3000, warmup=2000, chains=4,
                 control = list(adapt_delta = 0.8))
#plot(bg,pars="pi")
#plot(bh,pars=c(paste0("pi[",1:20,"]"))) approx 1/300 + 1/300
#summary(bg,pars=c(paste0("pi[",1:20,"]")))$summary

mns_bg<-summary(bg,pars=c(paste0("pi")))$summary[,'mean']
sum(mns_bg)

cbind(mns_bg*n,table(y))
```

```{r ex2a.4, cache=TRUE}
# logistic
bh <- sampling(ord_mod1, data=mod_data3h, seed=6472, 
                 iter=4000, warmup=2000, chains=4,
                 control = list(adapt_delta = 0.9))
```

```{r ex2a.5, cache=TRUE}
# log-log
bll <- sampling(ord_mod1, data=mod_data3ll, seed=6472, 
                 iter=3000, warmup=2000, chains=3,
                 control = list(adapt_delta = 0.9))
```

```{r ex2a.6.0, cache=TRUE, eval=FALSE}
# complementary log-log
#!! takes a really long time
#!! didn't converge w/  seed=12345 and:
# iter=3500, warmup=1500, chains=2, control = list(adapt_delta = 0.9)
# iter=5500, warmup=3500, chains=2, control = list(adapt_delta = 0.9)
# iter=5500, warmup=3500, chains=2, control = list(adapt_delta = 0.95)
# iter=4000, warmup=3000, chains=3, control = list(adapt_delta = 0.99)
# iter=7000, warmup=6000, chains=2, control = list(adapt_delta = 0.99)
bcll <- sampling(ord_mod1, data=mod_data3cll, seed=6472, 
                 iter=8000, warmup=7000, chains=2,
                 control = list(adapt_delta = 0.95, max_treedepth=15))
```

```{r,eval=FALSE}
traceplot(bcll,pars=c("cutpoints[1]"))
```

```{r ex2a.6.1, eval=FALSE}
# try optimizing?
bcll_opt <- optimizing(ord_mod1, data=mod_data3cll, seed=6472)
# get non-zero return code
```

```{r ex2a.6.2, eval=FALSE}
# try initializing at MLEs or 
#ll
#summary(bcll,pars=c(paste0("cutpoints[",1:299,"]")))$summary[,'mean']

# somehow even slower!!
#init_list0 <- list(cutpoints=-coef(ll)[1:299])
#init_list <- lapply(1:2, function(x) lapply(init_list0 ,jitter))
#bcll_i <- sampling(ord_mod1, data=mod_data3cll, seed=6472, 
#                 iter=4000, warmup=2000, chains=2, init=init_list,
#                 control = list(adapt_delta = 0.95, max_treedepth=15))

# try initialize based on pis
mod_data3cll_pis <- table(mod_data3cll$y)/mod_data3cll$N
mod_data3cll_cp_init <- log(-log1p(-cumsum(mod_data3cll_pis)))[1:299]

ini <- list(cutpoints=mod_data3cll_cp_init)
ini_list <- lapply(1:2, function(x) lapply(ini,jitter))

# still slow and doesn't converge
bcll_i2 <- sampling(ord_mod1, data=mod_data3cll, seed=6472, 
                 iter=4000, warmup=2000, chains=2, init=ini_list,
                 control = list(adapt_delta = 0.9, max_treedepth=10))

```

```{r ex2a.7, cache=TRUE}
bcau <- sampling(ord_mod1, data=mod_data3cau, seed=6472, 
                 iter=4000, warmup=2000, chains=3,
                 control = list(adapt_delta = 0.9))
```

```{r ex2a.8, cache=TRUE}
summary(bg,pars=c("b[1]"))$summary

uylevs<-length(levels(dat2$y))
levseq<-1:(uylevs - 1)

#rownames(summary(bg,pars=paste0("cutpoints[",levseq,"]"))$summary)

get_cutpoints<-function(fit,summ_stat='mean'){
  summary(fit, pars=paste0("cutpoints[",levseq,"]"), probs=c(0.5))$summary[,summ_stat]
}

cpg<-get_cutpoints(bg)
cph<-get_cutpoints(bh)
cpll<-get_cutpoints(bll)
cpcll<-get_cutpoints(bcll)

# cloglog and loglog seem backward from specified links but matches orm output
# this is because orm estimates P(Y>y) rather than P(Y<y) so 'opposite' link function needs to be used for non-symmetric links
z2 <- list(logistic=list(x=levseq, y=-1*cpg),
           probit=list(x=levseq, y=-1*cph),
           loglog=list(x=levseq, y=-1*cpcll),
           cloglog=list(x=levseq, y=-1*cpll)
)

labcurve(z2, pl=TRUE, col=1:4, ylab='Intercept')

```

## orm() example 2b

```{r ex2b.1, cache=TRUE}

# mean in each group
tapply(as.numeric(as.character(dat2$y)), dat2$x1, mean)

# using linear regression
dat2$y_num<-as.numeric(as.character(dat2$y))
coef(ols(y_num~x1,data=dat2))

# make Mean function to calc mean from orm
#what is importance of interceptRef?? why is this a parameter?
m <- Mean(g)
m(w3 <- k[3] + k['x1']*c(0,1), interceptRef=3)
m(w47 <- k[47] + k['x1']*c(0,1), interceptRef=47)
m(w <- k[1] + k['x1']*c(0,1), interceptRef=1)

mh <- Mean(h)
mh( wh <- coef(h)[1] + coef(h)['x1']*c(0,1), interceptRef=1 )
```

```{r ex2b.2, cache=TRUE}
# using CPM
m_bg <- getMean(bg, mod_data3, newdata=data.frame(x1=c(0,1)))
t(m_bg[,"med_mn"])

m_bh <- getMean(bh, mod_data3h, newdata=data.frame(x1=c(0,1)))
t(m_bh[,"med_mn"])
```

```{r ex2b.3, eval=FALSE}
qu <- Quantile(g)
# Compare model estimated and empirical quantiles
cq <- function(y) {
  cat('orm est', qu(.1, w, interceptRef = 1),
      'emp est', tapply(y, x1, quantile, probs=.1), '\n')
  cat('orm est', qu(.5, w, interceptRef = 1), 
      'emp est', tapply(y, x1, quantile, probs=.5), '\n')
  cat('orm est', qu(.9, w, interceptRef = 1),
      'emp est', tapply(y, x1, quantile, probs=.9), '\n')
}
cq(y)

q0.1_bg <-getQuantile(bg, mod_data3, newdata = data.frame(x1=c(0,1)), q=0.1)
q0.5_bg <-getQuantile(bg, mod_data3, newdata = data.frame(x1=c(0,1)), q=0.5)
q0.9_bg <-getQuantile(bg, mod_data3, newdata = data.frame(x1=c(0,1)), q=0.9)

rbind(t(q0.1_bg[,'med_qtile']), t(q0.5_bg[,'med_qtile']), t(q0.9_bg[,'med_qtile']))
```

## orm() example 3

```{r ex3.1, eval=FALSE}
# Try on log-normal model
rm(g,k,m,qu)
g <- orm(exp(y) ~ x1)
g
k <- coef(g)
plot(k[1:i])
m <- Mean(g)
# mean from orm() model
m(w <- k[1] + k['x1']*c(0,1),interceptRef = 1)
# empirical mean
tapply(exp(y), x1, mean)


qu <- Quantile(g)
# quantiles from orm() and empirical
cq(exp(y))
```


```{r ex3.2, cache=TRUE}
# Try on log-normal model
rm(bg, q0.1_bg, q0.5_bg, q0.9_bg)
dat2$exp_y<- ordered(exp(dat2$y_num))
mod_data4 <- mkStanDat(dat2, outcome="exp_y", preds = c("x1"), link=1)

bg <- sampling(ord_mod1, data=mod_data4, seed=99183, 
                 iter=4000, warmup=2000, chains=4,
                 control = list(adapt_delta = 0.8))

m_bg <- getMean(bg, mod_data4, newdata=data.frame(x1=c(0,1)))
t(m_bg[,"med_mn"])

q0.1_bg <-getQuantile(bg, mod_data4, newdata = data.frame(x1=c(0,1)), q=0.1)
q0.5_bg <-getQuantile(bg, mod_data4, newdata = data.frame(x1=c(0,1)), q=0.5)
q0.9_bg <-getQuantile(bg, mod_data4, newdata = data.frame(x1=c(0,1)), q=0.9)

rbind(t(q0.1_bg[,'med_qtile']), t(q0.5_bg[,'med_qtile']), t(q0.9_bg[,'med_qtile']))
```


## orm() example 4a

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


dd <- datadist(x1); options(datadist='dd')
f <- ols(y ~ x1)
g <- orm(y ~ x1, family=probit)
h <- orm(y ~ x1, family=logistic)
w <- orm(y ~ x1, family=cloglog)
mg <- Mean(g); mh <- Mean(h); mw <- Mean(w)
r <- rbind(ols      = Predict(f, conf.int=FALSE),
           probit   = Predict(g, conf.int=FALSE, fun=mg),
           logistic = Predict(h, conf.int=FALSE, fun=mh),
           cloglog  = Predict(w, conf.int=FALSE, fun=mw))
plot(r, groups='.set.')
```


```{r ex4a.2, cache=TRUE}
mod_data4g <- mod_data4h <- mkStanDat(dat3, outcome="y", preds = c("x1"), link=1)
mod_data4g$link <- 2
rm(bg,bh)
```

```{r ex4a.3, cache=TRUE}
bf <- stan_glm(y_num ~ x1, data=dat3, family=gaussian)
```

```{r ex4a.4, cache=TRUE}
bg <- sampling(ord_mod1, data=mod_data4g, seed=6472, 
                 iter=3250, warmup=2000, chains=4,
                 control = list(adapt_delta = 0.8))
```

```{r ex4a.5, cache=TRUE}
bh <- sampling(ord_mod1, data=mod_data4h, seed=6472, 
                 iter=3250, warmup=2000, chains=4,
                 control = list(adapt_delta = 0.8))
```


```{r, eval=FALSE}
# evaluate time to run getCDF vs. getCDF_dt
system.time( cdf_bg <- getCDF(bg, mod_data4g, newdata=data.frame(x1=c(-1,0,1))) )
system.time( cdf_bg_dt <- getCDF_dt(bg, mod_data4g, newdata=data.frame(x1=c(-1,0,1))) )

callCDF <- function() getCDF(bg, mod_data4g, newdata=data.frame(x1=c(-1,0,1)))
callCDF_dt <- function() getCDF_dt(bg, mod_data4g, newdata=data.frame(x1=c(-1,0,1)))

library(microbenchmark)
res <- microbenchmark(callCDF(), callCDF_dt(), times=25)

print(res)
ggplot2::autoplot(res)


cdf_bg_dt %>% filter(ndrow %in% c(1,2,3)) %>% ggplot(aes(group=ndrow)) +
  geom_ribbon(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("") + ylab("Conditional CDF") 

```

```{r, eval=FALSE}
#plots, etc for 11/13/19 seminar
bg_df<-as.data.frame(bg)

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

#head(bg_df[,c("cutpoints[1]","cutpoints[2]","cutpoints[3]","cutpoints[73]","cutpoints[74]","b[1]")])


dir<-getwd()
figdir<-file.path(dir,"biostat_sem","fig")

library(pammtools)

cdf_bg <- getCDF(bg, mod_data4g, newdata=data.frame(x1=c(-1,0,1))) 

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 = -1", "x = 0", "x = 1")) + 
  scale_color_discrete(name = "covariate value", 
                       labels=c("x = -1", "x = 0", "x = 1")) +
  stat_function(fun=pnorm,color="darkgreen",
                linetype=2, alpha=0.4) + 
  stat_function(fun=function(x) pnorm(x,0.9), color="blue",
                linetype=2, alpha=0.4) +
  stat_function(fun=function(x) pnorm(x,-0.9),color="red",
                linetype=2, alpha=0.4)

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


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

#ggplot(mn_dat,aes(x=mn,group=x1,fill=factor(ndrow)))+geom_density(alpha=0.4)

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 = -1", "x = 0", "x = 1")) + 
  xlab("") + ylab("conditional mean density")
ggsave(file.path(figdir,"cond_mn.png"),width=6,height=3)

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


q50_dat<-getQuantile(bg, mod_data4g, newdata=data.frame(x1=c(-1,0,1)),q=0.50,summ=FALSE)
ggplot(q50_dat,aes(x=qtile,fill=factor(x1)))+geom_density(alpha=0.6,color=NA)+ scale_fill_discrete(name = "covariate value", 
                      labels=c("x = -1", "x = 0", "x = 1")) + 
  xlab("") + ylab("conditional median density")
ggsave(file.path(figdir,"cond_md.png"),width=6,height=3)

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

mean(q50_x0_samps > 0 | q50_x0_samps < -0.5)

```

```{r ex4a.6, cache=TRUE}
xseq <- seq(-1.5, 1.6, length=50)
#this takes a while to run
# see if function can be optimized
m_bg <- getMean(bg, mod_data4g, newdata=data.frame(x1=xseq))
m_bh <- getMean(bh, mod_data4h, newdata=data.frame(x1=xseq))

m_bf <- predict(bf, newdata=data.frame(x1=xseq))

plot(xseq,m_bf,"l",col=2)
lines(xseq,t(m_bg[,"med_mn"]),"l",col=3)
```


## orm example 4b

```{r ex4b.1, cache=TRUE}
# Compare predicted 0.8 quantile with quantile regression
qu <- Quantile(g)
qu80 <- function(lp) qu(.8, lp)
f <- Rq(y ~ x1, tau=.8)
r <- rbind(probit   = Predict(g, conf.int=FALSE, fun=qu80),
           quantreg = Predict(f, conf.int=FALSE))
plot(r, groups='.set.')
```


```{r ex4b.2, cache=TRUE}
library(brms)
#Bayesian quantile regression with brms
q_bf <- brm(bf(y_num ~ x1, quantile = 0.8), data = dat3, family = asym_laplace())
q_bf_pred<-fitted(q_bf, dpar="mu", newdata=data.frame(x1=xseq) )

q_bg <- getQuantile(bg, mod_data4g, newdata=data.frame(x1=xseq), q=0.8)

plot(xseq,q_bf_pred[,1],"l",col=2,ylim=c(-1,2.75))
lines(xseq,t(q_bg[,"med_qtile"]),"l",col=3)

```

## orm() example 5

Verify transformation invariance of ordinal regression

```{r ex5.1, cache=TRUE}
ga <- orm(exp(y) ~ x1, family=probit)
qua <- Quantile(ga)
qua80 <- function(lp) log(qua(.8, lp))
ra <- rbind(logprobit = Predict(ga, conf.int=FALSE, fun=qua80),
           probit    = Predict(g,  conf.int=FALSE, fun=qu80))
plot(ra, groups='.set.')
```

```{r ex5.2, cache=TRUE}
dat3$exp_y <- ordered(exp(dat3$y_num))
mod_data5g <- mkStanDat(dat3, outcome="exp_y", preds = c("x1"), link=1)

bga <- sampling(ord_mod1, data=mod_data5g, seed=6472, 
                 iter=3250, warmup=2000, chains=4,
                 control = list(adapt_delta = 0.8))

q_bga <- getQuantile(bga, mod_data5g, newdata=data.frame(x1=xseq), q=0.8)

plot(xseq,log(t(q_bga[,"med_qtile"])),"l",col=2)
lines(xseq,t(q_bg[,"med_qtile"]),"l",col=3)
```



```{r ex5.3, cache=TRUE}
# Try the same with quantile regression  
# Need to transform x1 using spline to match previous model
fa0 <- Rq(exp(y) ~ x1, tau=.8)
fa <- Rq(exp(y) ~ rcs(x1,5), tau=.8)
r <- rbind(qr    = Predict(f, conf.int=FALSE),
           logqr0 = Predict(fa0, conf.int=FALSE, fun=log),
           logqr = Predict(fa, conf.int=FALSE, fun=log))
plot(r, groups='.set.')
```

```{r ex5.4, cache=TRUE}
# Same with Bayesian quantile regression.  again, Need to transform x1 or will be off!
q_bfa0 <- brm(bf(exp(y_num) ~ x1, quantile = 0.8), data = dat3, family = asym_laplace())

q_bfa <- brm(bf(exp(y_num) ~ s(x1,k=5), quantile = 0.8), data = dat3, family = asym_laplace())

q_bfa0_pred<-fitted(q_bfa0, dpar="mu", newdata=data.frame(x1=xseq) )
q_bfa_pred<-fitted(q_bfa, dpar="mu", newdata=data.frame(x1=xseq) )

plot(xseq, q_bf_pred[,1], "l", col=2, ylim=c(-1,2.75))
lines(xseq, log(q_bfa0_pred[,1]), "l", col=3,ylim=c(-1,2.75))
lines(xseq, log(q_bfa_pred[,1]), "l", col=4,ylim=c(-1,2.75))
```