brms_practice.Rmd

---
title: "brms_practice"
author: "Julin Maloof"
date: "8/16/2020"
output: 
  html_document: 
    keep_md: yes
---

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

```{r}
library(rethinking)
library(brms)
library(tidyverse)
```

## Assignment

Revisit the following homework problems and try to fit the with brms.  Make your first attempt without looking at the rethinking to brms translation, but if you get stuck definitely look!  Compare the coefficients or predictions that you obtain with brms and those with quap or ulam.

* 4H1, 4H2 (you probably need the function `posterior_predict()`)
* From chapter 8 I assigned a tomato problem from my data "Use the tomato.csv data set and evaluate whether hypocotyl length ("hyp") is affected by shade ("trt"), species ("species") and their interaction."
* From chapter 9: 8M1 (remember that the problem numbers were offset it is actually called 9M1 in the Nov 24 PDF)

## 4H1

### original answer

_The weights listed below were recorded in the !Kung census, but heights were not recorded for these individuals. Provide predicted heights and 89% intervals (either HPDI or PI) for each of these individuals. That is, fill in the table below, using model-based predictions._

```{r}
weights <- c(46.95,43.72,64.78,32.59,54.63)
```

These all seem to be adult weights, so I am going to use the simple regression model on adults

```{r}
#4.42
data(Howell1)
d <- Howell1
d2 <- d[ d$age >= 18 , ]

# define the average weight, x-bar
xbar <- mean(d2$weight)

d2 <- c(d2, xbar=xbar)
str(d2)

# fit model
m4.3 <- ulam(
    alist(
        height ~ dnorm( mu , sigma ) ,
        mu <- a + b*( weight - xbar ) ,
        a ~ dnorm( 178 , 20 ) ,
        b ~ dlnorm( 0 , 1 ) ,
        sigma ~ dexp(1)
    ) , chains = 4, cores = 4, refresh=0,
    data=d2 )

```


```{r}
sim.heights <- sim(fit = m4.3,
                   data = data.frame(weight=weights),
                   n=1e4)
dim(sim.heights)
head(sim.heights)
```

```{r}
resultsulam <- tibble(
  individual=1:5,
  weight=weights,
  predicted.height=apply(sim.heights,2,mean)
  )
resultsulam <- cbind(resultsulam,t(apply(sim.heights, 2, HPDI)))
knitr::kable(resultsulam,digits=2)
```

### brms

```{r}
get_prior(height ~ ( weight - xbar ), data = d2)
```


```{r}
d2$weight2 <-d2$weight-d2$xbar
m4.3brms <- brm(
  height ~ weight2  ,
  prior = c(set_prior("normal(178,20)", class="Intercept"),
            set_prior("lognormal(0, 1)", class="b", lb = 0),
            set_prior("exponential(1)", class="sigma")),
    data=d2, refresh=0 )
```
```{r}
precis(m4.3)
```


```{r}
summary(m4.3brms, prob=0.89)
```

now make predictions...
```{r}
brmspost <- posterior_predict(m4.3brms, newdata = data.frame(weight2=(weights-xbar)))

brmsresults <- tibble(
  individual=1:5,
  weight=weights,
  predicted.height=apply(brmspost,2,mean)
  )
brmsresults <- cbind(brmsresults,t(apply(brmspost, 2, HPDI)))
knitr::kable(brmsresults,digits=2)
```

alternative
```{r}
predict(m4.3brms, 
        newdata = data.frame(weight2=(weights-xbar)),
        probs = c(0.055, 0.945)) %>% 
  knitr::kable(digits = 2)
```

compare to ulam
```{r}
knitr::kable(resultsulam,digits=2)
```
## 4H2

### original

## 4H2. 
_Select out all the rows in the Howell1 data with ages below 18 years of age. If you do it right, you should end up with a new data frame with 192 rows in it._

```{r}
d3 <- Howell1 %>% as_tibble() %>% filter(age < 18)
head(d3)
```


_(a) Fit a linear regression to these data, using quap. Present and interpret the estimates. For every 10 units of increase in weight, how much taller does the model predict a child gets?_

first let's take a look
```{r}
qplot(weight, height, data = d3)
```

so not linear here, but based on problem I think we fit a simple regression model

```{r}
# define the average weight, x-bar
xbar <- mean(d3$weight)

d3 <- c(d3, xbar=xbar)

# fit model
m4H2 <- ulam(
    alist(
        height ~ dnorm( mu , sigma ) ,
        mu <- a + b*( weight - xbar ) ,
        a ~ dnorm( 100 , 20 ) ,
        b ~ dlnorm( 0 , 1 ) ,
        sigma ~ dexp( 1 )
    ) ,
    chains = 4,
    cores = 4,
    refresh = 0,
    data=d3 )
```

```{r}
precis(m4H2)
```

For every 10 units increase in weight, the model predicts an increase of 27.2 height

### brms

```{r}
# define the average weight, x-bar
d3$weight2 <- d3$weight - d3$xbar

# fit model
m4H2brms <- brm(
        height ~  weight2,
  prior = c(set_prior("normal(100,20)", class="Intercept"),
            set_prior("lognormal(0, 1)", class="b", lb = 0),
            set_prior("exponential(1)", class="sigma")),
    data=d3, 
  refresh=0 )

```

```{r}
summary(m4H2brms, prob=0.89)
```


_(b) Plot the raw data, with height on the vertical axis and weight on the horizontal axis. Super- impose the MAP regression line and 89% HPDI for the mean. Also superimpose the 89% HPDI for predicted heights._

### original
```{r}
post <- extract.samples(m4H2)
weight.seq <- seq.int(from=min(d3$weight)-1, to=max(d3$weight)+1, by=1)
mu <- link(m4H2, data=data.frame(weight=weight.seq))
mu.mean <- apply( mu , 2 , mean )
mu.HPDI <- apply( mu , 2 , HPDI , prob=0.89 )
sim <- sim(m4H2, data=data.frame(weight=weight.seq))
sim.HPDI <- apply(sim, 2, HPDI, prob=0.89)
```

plot it
```{r}
plot.data <- as_tibble(cbind(weight.seq, mu.mean, t(mu.HPDI), t(sim.HPDI))) %>%
  rename(weight=weight.seq, mu.89l=`|0.89`, mu.89u=`0.89|`, sim.89l=V5, sim.89u=V6)
head(plot.data)
```

```{r}
obs.plot.data <- tibble(weight=d3$weight, height=d3$height)
plot.data %>% 
  ggplot(aes(x=weight)) +
  geom_line(aes(y=mu.mean)) +
  geom_ribbon(aes(ymin=sim.89l, ymax=sim.89u), fill="gray80") +
  geom_ribbon(aes(ymin=mu.89l, ymax=mu.89u), fill="gray60") +
  geom_point(aes(y=height), data=obs.plot.data,color="blue", alpha=.3) +
  ylab("height")
```

### brms

```{r}
mu.brms <- posterior_epred(m4H2brms, newdata = data.frame(weight2=weight.seq-xbar))
mu.mean.brms <- apply( mu.brms , 2 , mean )
mu.HPDI.brms <- apply( mu.brms , 2 , HPDI , prob=0.89 )
pred <- predict(m4H2brms, newdata=data.frame(weight2=weight.seq-xbar), probs = c(0.055, 0.945))
```

plot it
```{r}
plot.data <- as_tibble(cbind(weight.seq, mu.mean.brms, t(mu.HPDI.brms), pred[,c("Q5.5", "Q94.5")])) %>%
  rename(weight=weight.seq, mu.89l=`|0.89`, mu.89u=`0.89|`)
head(plot.data)
```

```{r}
plot.data %>% 
  ggplot(aes(x=weight)) +
  geom_ribbon(aes(ymin=Q5.5, ymax=Q94.5), fill="gray80") +
  geom_ribbon(aes(ymin=mu.89l, ymax=mu.89u), fill="gray60") +
  geom_line(aes(y=mu.mean.brms)) +
  geom_point(aes(y=height), data=obs.plot.data,color="blue", alpha=.3) +
  ylab("height")
```
## Tomato
_Use the tomato.csv (attached) data set and evaluate whether hypocotyl length ("hyp") is affected by shade ("trt"), species ("species") and their interaction._

### original

```{r}
d <- read_csv("Tomato.csv") %>%
  select(hyp, trt, species) %>%
  na.omit()
head(d)
```

Make a plot
```{r}
d %>%
  group_by(species,trt) %>%
  summarize(mean=mean(hyp), 
            sem=sd(hyp)/sqrt(n()), 
            ymax=mean+sem,
            ymin=mean-sem) %>%
  ggplot(aes(x=species, y=mean, ymax=ymax, ymin=ymin, fill=trt)) +
  geom_col(position = "dodge") +
  geom_errorbar(position = position_dodge(width=.9), width=.5)
```

make indices for the factors:
```{r}
d <- d %>%
  mutate(species_i = as.integer(as.factor(species)),
         trt_i = as.integer(as.factor(trt))-1L)
```


fit non interaction model
```{r}
d2 <- d %>% select(hyp, species_i, trt_i)
m1 <- ulam(flist = alist(
  hyp ~ dnorm(mu, sigma),
  mu <- a[species_i] + b*trt_i, # one beta coefficient
  a[species_i] ~ dnorm(25, 5),
  b ~ dnorm(0, 5),
  sigma ~ dexp(1)),
  log_lik = TRUE,
  data=d2, chains = 4, cores = 4, refresh = 0)
```


```{r}
precis(m1, depth=2)
```

now the interaction model:
```{r}
d2 <- d %>% select(hyp, species_i, trt_i)
m2 <- ulam(flist = alist(
  hyp ~ dnorm(mu, sigma),
  mu <- a[species_i] + b_int[species_i]*trt_i, # a beta coefficent for each species
  a[species_i] ~ dnorm(25, 5),
  b_int[species_i] ~ dnorm(0, 1),
  sigma ~ dexp(1)),
  data=d2, chains = 4, cores = 4,
  log_lik = TRUE,
  refresh=0)
```

```{r}
precis(m2, depth = 2)
```
```{r}
compare(m1, m2)
```

non interaction model preferred

```{r}
plot(coeftab(m2))
```


### brms

non-interaction
```{r}
get_prior(hyp~-1 + species + trt, data=d)
```


```{r}
m1brms <- brm(hyp ~ -1 + species + trt,
              prior = c(set_prior("normal(25, 5)",class  = "b"),
                        set_prior("normal(0, 5)", class="b", coef="trtL"),
                        set_prior("exponential(1)", class="sigma")),
              data = d,
              refresh = 0)
m1brms <- add_criterion(m1brms, "waic")

```

make sure prior as what we think:

```{r}
prior_summary(m1brms)
```

yes

```{r}
precis(m1, depth = 2)
```


```{r}
summary(m1brms)
```

```{r}
get_prior(hyp ~ -1 + species*trt, data = d)
```
Hmm, this is still fitting a main "trtL" term and I am not sure how to get rid of that.  So I think I will do the default parameterization...

```{r}
get_prior(hyp ~ species*trt, data = d)
```

```{r}
m2brms <- brm(hyp ~ species * trt,
              prior = c(set_prior("normal(25, 5)", class="Intercept"),
                        set_prior("normal(0, 5)",class  = "b"),
                        set_prior("exponential(1)", class="sigma")),
              data = d,
              refresh = 0)
m2brms <- add_criterion(m2brms, "waic")
```

```{r}
summary(m2brms, prob=.89)
```

```{r}
print(loo_compare(m1brms, m2brms, criterion="waic"), simplify=FALSE)
```
better try normal paramteriazation for m1

```{r}
m1brmsalt <- brm(hyp ~ species + trt,
              prior = c(set_prior("normal(25, 5)", class="Intercept"),
                        set_prior("normal(0, 5)",class  = "b"),
                        set_prior("exponential(1)", class="sigma")),
              data = d,
              refresh = 0)
m1brmsalt <- add_criterion(m1brmsalt, criterion = "waic")
```

```{r}
print(loo_compare(m1brms, m1brmsalt, m2brms, criterion = "waic"), simplify=FALSE)
```


## 9M1