Chapter12.Rmd

---
title: "Chapter12"
author: "Julin N Maloof"
date: "11/13/2019"
output: 
  html_document: 
    keep_md: yes
---

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

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

# Problems

## 11E1
_What is the difference between an ordered categorical variable and an unordered one? Define and then give an example of each._

An unordered categorical variable is a grouping variable where the groups do not have any particular relationship to one another.  Example: US states.

An ordered categorical cariable is a goruping variable where the different groups do have a relationship with one another (can be sorted) but the step size between groups is unknown.  Example: small, medium, large.


## 11E2

_What kind of link function does an ordered logistic regression employ? How does it differ from an ordinary logit link?_

Still a logit, but here the we are using cumulative probabilities.  The probability of obserbing that outcome or any "lesser" outcome.


## 11E3
_When count data are zero-inflated, using a model that ignores zero-inflation will tend to induce which kind of inferential error?_

Will underestimate the rate of the event occuring.

## 11E4
_Over-dispersion is common in count data. Give an example of a natural process that might produce over-dispersed counts. Can you also give an example of a process that might produce under- dispersed counts?_

The rate of leaves falling off of a evergeen tree could be overdispersed if you did not take weather/wind conditions into account.

One would expect underdispersion to occur in a system that has a negative feedback loop (depending on the time delay)

## 11M1

_At a certain university, employees are annually rated from 1 to 4 on their productivity, with 1 being least productive and 4 most productive. In a certain department at this certain university in a certain year, the numbers of employees receiving each rating were (from 1 to 4): 12, 36, 7, 41. Compute the log cumulative odds of each rating._

```{r}
ratings <- c(12,36,7,41)
names(ratings) <- as.character(1:4)
ratings
probs <- ratings / sum(ratings) # probability of each rating
cumprobs <- cumsum(probs)
cumprobs
```


```{r}
cum_odds <- map_dbl(1:4, ~ cumprobs[.]/ (1-cumprobs[.]) )
cum_odds

log(cum_odds)
```

or...
```{r}
cum_log_odds <- map_dbl(1:4, ~ log(cumprobs[.]) - log(1-cumprobs[.]))

cum_log_odds
```

## 11M2
_Make a version of Figure 12.5 for the employee ratings data given just above._
```{r}
dat <- tibble(
  rating = names(probs),
  probability = probs,
  cumulative = cumprobs
)
dat
```

```{r}
dat %>% 
  ggplot(aes(x=rating, xend=rating)) +
  geom_segment(aes(yend=cumulative), y=0, col="gray70", lwd=2, position=position_nudge(x=-0.02)) +
  geom_segment(aes(y=cumulative-probability, yend=cumulative), col="skyblue", lwd=2, position = position_nudge(x=.02)) +
  geom_line(aes(x=rating, y = cumulative), group=1) +
  geom_point(aes(y = cumulative), shape=21, fill="white", size=2) +
  geom_text(aes(y = (cumulative - probability/2), label=rating), color="skyblue", position=position_nudge(x=.1)) +
  ylab("cumulative proportion") +
  theme_bw()
  
```


## 11M3 (optional)

_Can you modify the derivation of the zero-inflated Poisson distribution (ZIPoisson) from the chapter to construct a zero-inflated binomial distribution?_

## 11H1

_In 2014, a paper was published that was entitled “Female hurricanes are deadlier than male hurricanes.” As the title suggests, the paper claimed that hurricanes with female names have caused greater loss of life, and the explanation given is that people unconsciously rate female hurricanes as less dangerous and so are less likely to evacuate._

_Statisticians severely criticized the paper after publication. Here, you’ll explore the complete data used in the paper and consider the hypothesis that hurricanes with female names are deadlier. _

_In this problem, you’ll focus on predicting deaths using femininity of each hurricane’s name._

_Fit and interpret the simplest possible model, a Poisson model of deaths using femininity as a predictor. You can use map or map2stan. Compare the model to an intercept-only Poisson model of deaths. How strong is the association between femininity of name and deaths? Which storms does the model fit (retrodict) well? Which storms does it fit poorly?_

```{r}
library(rethinking)
data(Hurricanes)
?Hurricanes
head(Hurricanes)
```

```{r}
Hurricanes %>%
  ggplot(aes(x=femininity, y=deaths)) +
  geom_point() +
  geom_smooth()

Hurricanes %>%
  ggplot(aes(x=femininity, y=deaths)) +
  geom_point() +
  geom_smooth() +
  scale_x_log10()
```

Doesn't look very promising...

Prior predictions.  deaths can go higher, so should I modify alpha?
```{r}
curve(dlnorm(x, 0, 10), from=0, to=100, n=200)
curve(dlnorm(x, 3, .5), from=0, to=100, n=200, add=TRUE, col="blue")
curve(dlnorm(x, 3, 1), from=0, to=100, n=200, add=TRUE, col="red")
```

But, we are really thinking about average rate of deaths, which is around 20, so seems OK

Intercept only model
```{r}
mh1 <- ulam(alist( 
  deaths ~ dpois(lambda),
  log(lambda) <- alpha,
  alpha ~ dnorm(3, .5)),
  data=Hurricanes,
  chains = 4,
  cores = 4,
  log_lik = TRUE)
```

```{r}
trankplot(mh1)
traceplot(mh1)
```

```{r}
precis(mh1)
```

```{r}
exp(3.03)
mean(Hurricanes$deaths)
```


Model with feminity as predictor.  If I don't standardize femininity, what prior should I use?

```{r}
set.seed(10)
N <- 100
a <- rnorm( N , 3 , .5 )
b <- rnorm( N , 0 , 0.2 )
plot( NULL , xlim=c(1,10) , ylim=c(0,500) )
for ( i in 1:N ) curve( exp( a[i] + b[i]*x ) , add=TRUE , col=col.alpha("black",0.5) )
```

Seems reasonable

```{r}
mh1.2 <- ulam(alist( 
  deaths ~ dpois(lambda),
  log(lambda) <- alpha + b*femininity,
  alpha ~ dnorm(3, .5),
  b ~ dnorm(0, 0.2)),
  data=Hurricanes,
  chains = 4,
  cores = 4,
  log_lik = TRUE)
```

```{r}
traceplot(mh1.2)
trankplot(mh1.2)
pairs(mh1.2)
```

```{r}
precis(mh1.2)
```


```{r}
compare(mh1, mh1.2)
```


While model h1.2 shows strong support for a non-zero beta (effect of feminity), it does not fit the data significantly better than the intercept only model.

## 11H2

_Counts are nearly always over-dispersed relative to Poisson. So fit a gamma-Poisson (aka negative-binomial) model to predict deaths using femininity. Show that the over-dispersed model no longer shows as precise a positive association between femininity and deaths, with an 89% interval that overlaps zero. Can you explain why the association diminished in strength?_

```{r}
mh2 <- ulam(alist( 
  deaths ~ dgampois(lambda, phi),
  log(lambda) <- alpha + b*femininity,
  alpha ~ dnorm(3, .5),
  b ~ dnorm(0, 0.2),
  phi ~ dexp(1)),
  data=Hurricanes,
  chains = 4,
  cores = 4,
  log_lik = TRUE)
```

```{r}
traceplot(mh2)
trankplot(mh2)
pairs(mh2)
```

```{r}
precis(mh2)
```

Now we are much less certain about `b` being positive.  The support interval crosses 0.

Optional: 11H3, 11H4

## 11H6

_The data in data(Fish) are records of visits to a national park. See ?Fish for details. The question of interest is how many fish an average visitor takes per hour, when fishing. The problem is that not everyone tried to fish, so the fish_caught numbers are zero-inflated. As with the monks example in the chapter, there is a process that determines who is fishing (working) and another process that determines fish per hour (manuscripts per day), conditional on fishing (working). We want to model both. Otherwise we’ll end up with an underestimate of rate of fish extraction from the park._

_You will model these data using zero-inflated Poisson GLMs. Predict fish_caught as a function of any of the other variables you think are relevant. One thing you must do, however, is use a proper Poisson offset/exposure in the Poisson portion of the zero-inflated model. Then use the hours variable to construct the offset. This will adjust the model for the differing amount of time individuals spent in the park._

```{r}
data(Fish)
?Fish
head(Fish)
summary(Fish)
```

Livebait could be important as could child, and, of course, persons.  If I am including children then persons should get split into adults, because otherwise children are getting counted twice.

First fit this without ZI:

```{r}
f <- Fish
f$log_hours <- log(f$hours)
f$adults <- f$persons - f$child
mh6.1 <- ulam(alist(
  fish_caught ~ dpois(lambda),
  log(lambda) <- log_hours + a + bLive*livebait + bAdult*adults + bChild*child,
  a ~ dnorm(3, .5),
  c(bLive, bAdult, bChild) ~ dnorm(0, .5)),
  data = f,
  chains = 4,
  cores = 4,
  log_lik = TRUE)
```

```{r}
traceplot(mh6.1)
trankplot(mh6.1)
pairs(mh6.1)
```

```{r}
precis(mh6.1)
```


Now the zero inflated
```{r}
mh6.2 <- ulam(alist(
  fish_caught ~ dzipois(p, lambda),
  logit(p) <- ap,
  log(lambda) <- log_hours + a + bLive*livebait + bAdult*adults + bChild*child,
  ap ~ dnorm(0, 1.5),
  a ~ dnorm(3, .5),
  c(bLive, bAdult, bChild) ~ dnorm(0, .5)),
  data = f,
  chains = 4,
  cores = 4,
  log_lik = TRUE)
```

```{r}
traceplot(mh6.2)
trankplot(mh6.2)
pairs(mh6.2)
```
```{r}
precis(mh6.1)
```


```{r}
precis(mh6.2)
```

```{r}
inv_logit(-1.07) # estimate 25% of parties don't fish
exp(-2.45) # average fish caught per hour at intercept condition.  Hmm this should be 0.  Center predictors?
exp(-2.45+.57) # acerage fish caught per hour by adult
```

Let's look at predictions

```{r}
pred1 <- link(mh6.1)
pred2 <- link(mh6.2)
pred2_fish <- pred2$lambda*pred2$p[,1] # have to multiple fishing success by probability of fishing
pred_obs <- as_tibble(
  cbind(f, 
        pred1=colMeans(pred1),
        pred2=colMeans(pred2_fish),
        low.89.1=apply(pred1,2,HPDI)[1,],
        high.89.1=apply(pred1,2,HPDI)[2,],
        low.89.2=apply(pred2_fish,2,HPDI)[1,],
        high.89.2=apply(pred2_fish,2,HPDI)[2,]))
head(pred_obs)
```

```{r}
pl <- pred_obs %>% 
  select(fish_caught, ends_with("1"), ends_with("2")) %>%
  gather(key="model", value="predicted", pred1, pred2) %>%
  ggplot(aes(x=fish_caught, y=predicted, color=model)) +
  geom_point(alpha=.5) +
  geom_abline(intercept=0, slope=1) +
  coord_fixed() +
  facet_wrap(~model)
pl
```

```{r}
pl + coord_cartesian(xlim=c(0,50))
```

```{r}
f %>% arrange(desc(fish_caught)) %>% head(20)
```


## PDF week 7 problem 1
_1. In the Trolley data—data(Trolley)—we saw how education level (modeled as an ordered category) is associated with responses. Is this association causal? One plausible confound is that education is also associated with age, through a causal process: People are older when they finish school than when they begin it._

_Reconsider the Trolley data in this light. Draw a DAG that represents hypothetical causal relationships among response, education, and age. Which statical model or models do you need to evaluate the causal influence of education on responses?_

```{r}
library(dagitty)
g <- dagitty("dag{
  E -> R;
  A -> R;
  A -> E
}")

coordinates(g) <- list(
  x=c(A=1,R=2,E=3),
  y=c(A=0,R=1,E=0))
plot(g)
```


_Fit these models to the trolley data. What do you conclude about the causal relationships among these three variables?_


```{r}
data(Trolley)
d <- Trolley
levels(d$edu)

edu_levels <- c( 6 , 1 , 8 , 4 , 7 , 2 , 5 , 3 )
d$edu_new <- edu_levels[ d$edu ]
```


```{r}
library(gtools)
set.seed(999)
delta <- rdirichlet( 10 , alpha=rep(2,7) )
```

```{r}
dat <- list(
  R = d$response ,
  action = d$action,
  intention = d$intention,
  contact = d$contact,
  age = scale(d$age),
  E = as.integer( d$edu_new ), # edu_new as an index
  alpha = rep(2.1,7) )           # delta prior

system.time({
  mpdf1 <- ulam(
  alist(
    R ~ ordered_logistic( phi , kappa ),
    phi <- bE*sum( delta_j[1:E] ) + bA*action + bI*intention + bC*contact + bAge*age,
    kappa ~ normal( 0 , 1.5 ),
    c(bA,bI,bC,bE, bAge) ~ normal( 0 , 1 ),
    vector[8]: delta_j <<- append_row( 0 , delta ),
    simplex[7]: delta ~ dirichlet( alpha )
  ),
  data=dat , chains=3 , cores=3, iter = 2000 )
})
```


```{r}
precis(mpdf1, depth = 2, omit="kappa")
```

OK this quite changes things.  Now education has a postive effect whereas age has a negative effect.  This suggests that indeed the apparent negative education effect was driven by Age and we have closed that backdoor.


## PDF week 7 problem  2

_Consider one more variable in the Trolley data: Gender. Suppose that gender might influence education as well as response directly. Draw the DAG now that includes response, education, age, and gender._

```{r}
g <- dagitty("dag{
  G -> R;
  G -> E;
  E -> R;
  A -> R;
  A -> E
}")

coordinates(g) <- list(
  x=c(A=1,R=2,E=3,G=4),
  y=c(A=0,R=1,E=0,G=0))
plot(g)
```

_Using only the DAG, is it possible that the inferences from Problem 1 are confounded by gender? If so, define any additional models you need to infer the causal influence of education on response. What do you conclude?_

I don't think that the inference about age is incorrect, but the one about education could be.

```{r}
d %>% ggplot(aes(x=as.factor(male), y = edu_new)) +
  geom_violin()
```


```{r}
dat <- list(
  R = d$response ,
  action = d$action,
  intention = d$intention,
  contact = d$contact,
  age = scale(d$age),
  male = d$male,
  E = as.integer( d$edu_new ), # edu_new as an index
  alpha = rep(2.1,7) )           # delta prior

system.time({
  mpdf2 <- ulam(
  alist(
    R ~ ordered_logistic( phi , kappa ),
    phi <- bE*sum( delta_j[1:E] ) + bA*action + bI*intention + bC*contact + bAge*age +bMale*male,
    kappa ~ normal( 0 , 1.5 ),
    c(bA,bI,bC,bE, bAge, bMale) ~ normal( 0 , 1 ),
    vector[8]: delta_j <<- append_row( 0 , delta ),
    simplex[7]: delta ~ dirichlet( alpha )
  ),
  data=dat , chains=3 , cores=3, iter = 2000 )
})
```

```{r}
precis(mpdf2, depth = 2, omit = "kappa")
```

Gender plays a major role in the Response.  When both Gender and Age are considered, there is no evidence that education is relevant.  Older people and women are both much more likely to judge a scenario as immoral.

# Book Code


```{r, eval=FALSE}
## R code 12.1
pbar <- 0.5
theta <- 5
curve( dbeta2(x,pbar,theta) , from=0 , to=1 ,
       xlab="probability" , ylab="Density" )
```


```{r, eval=FALSE}
## R code 12.2
library(rethinking)
data(UCBadmit)
d <- UCBadmit
d$gid <- ifelse( d$applicant.gender=="male" , 1L , 2L )
dat <- list( A=d$admit , N=d$applications , gid=d$gid )
m12.1 <- ulam(
  alist(
    A ~ dbetabinom( N , pbar , theta ),
    logit(pbar) <- a[gid],
    a[gid] ~ dnorm( 0 , 1.5 ),
    theta ~ dexp(1)
  ), data=dat , chains=4 )
```


```{r, eval=FALSE}
## R code 12.3
post <- extract.samples( m12.1 )
post$da <- post$a[,1] - post$a[,2]
precis( post , depth=2 )
```


```{r, eval=FALSE}
## R code 12.4
gid <- 2
# draw posterior mean beta distribution
curve( dbeta2(x,mean(logistic(post$a[,gid])),mean(post$theta)) , from=0 , to=1 ,
       ylab="Density" , xlab="probability admit", ylim=c(0,3) , lwd=2 )

# draw 50 beta distributions sampled from posterior
for ( i in 1:50 ) {
  p <- logistic( post$a[i,gid] )
  theta <- post$theta[i]
  curve( dbeta2(x,p,theta) , add=TRUE , col=col.alpha("black",0.2) )
}
mtext( "distribution of female admission rates" )
```


```{r, eval=FALSE}
## R code 12.5
#postcheck( m12.1 )
```


```{r, eval=FALSE}
## R code 12.6
library(rethinking)
data(Kline)
d <- Kline
d$P <- standardize( log(d$population) )
d$contact_id <- ifelse( d$contact=="high" , 2L , 1L )

dat2 <- list(
  T = d$total_tools,
  P = d$population,
  cid = d$contact_id )

m12.3 <- ulam(
  alist(
    T ~ dgampois( lambda , phi ),
    lambda <- exp(a[cid])*P^b[cid] / g,
    a[cid] ~ dnorm(1,1),
    b[cid] ~ dexp(1),
    g ~ dexp(1),
    phi ~ dexp(1)
  ), data=dat2 , chains=4 , log_lik=TRUE )
```


```{r, eval=FALSE}
## R code 12.7
# define parameters
prob_drink <- 0.2 # 20% of days
rate_work <- 1    # average 1 manuscript per day

# sample one year of production
N <- 365

# simulate days monks drink
set.seed(365)
drink <- rbinom( N , 1 , prob_drink )

# simulate manuscripts completed
y <- (1-drink)*rpois( N , rate_work )
```


```{r, eval=FALSE}
## R code 12.8
simplehist( y , xlab="manuscripts completed" , lwd=4 )
zeros_drink <- sum(drink)
zeros_work <- sum(y==0 & drink==0)
zeros_total <- sum(y==0)
lines( c(0,0) , c(zeros_work,zeros_total) , lwd=4 , col=rangi2 )
```


```{r, eval=FALSE}
## R code 12.9
m12.4 <- ulam(
  alist(
    y ~ dzipois( p , lambda ),
    logit(p) <- ap,
    log(lambda) <- al,
    ap ~ dnorm( -1.5 , 1 ),
    al ~ dnorm( 1 , 0.5 )
  ) , data=list(y=as.integer(y)) , chains=4 )
precis( m12.4 )
```


```{r, eval=FALSE}
## R code 12.10
inv_logit(-1.28) # probability drink
exp(0.01)       # rate finish manuscripts, when not drinking

## R code 12.11
m12.4_alt <- ulam(
  alist(
    y|y>0 ~ custom( log1m(p) + poisson_lpmf(y|lambda) ),
    y|y==0 ~ custom( log_mix( p , 0 , poisson_lpmf(0|lambda) ) ),
    logit(p) <- ap,
    log(lambda) <- al,
    ap ~ dnorm(-1.5,1),
    al ~ dnorm(1,0.5)
  ) , data=list(y=as.integer(y)) , chains=4 )
```


```{r, eval=FALSE}
## R code 12.12
library(rethinking)
data(Trolley)
d <- Trolley
```


```{r, eval=FALSE}
## R code 12.13
simplehist( d$response , xlim=c(1,7) , xlab="response" )
```


```{r, eval=FALSE}
## R code 12.14
# discrete proportion of each response value
pr_k <- table( d$response ) / nrow(d)

# cumsum converts to cumulative proportions
cum_pr_k <- cumsum( pr_k )

# plot
plot( 1:7 , cum_pr_k , type="b" , xlab="response" ,
      ylab="cumulative proportion" , ylim=c(0,1) )
```


```{r, eval=FALSE}
## R code 12.15
logit <- function(x) log(x/(1-x)) # convenience function
( lco <- logit( cum_pr_k ) )
```


```{r, eval=FALSE}
## R code 12.16
m12.5 <- ulam(
  alist(
    R ~ dordlogit( 0 , cutpoints ),
    cutpoints ~ dnorm( 0 , 1.5 )
  ) ,
  data=list( R=d$response ), chains=4 , cores=4 )
```


```{r, eval=FALSE}
## R code 12.17
m12.5q <- quap(
  alist(
    response ~ dordlogit( 0 , c(a1,a2,a3,a4,a5,a6) ),
    c(a1,a2,a3,a4,a5,a6) ~ dnorm( 0 , 1.5 )
  ) , data=d ,
  start=list(a1=-2,a2=-1,a3=0,a4=1,a5=2,a6=2.5) )
```


```{r, eval=FALSE}
## R code 12.18
precis( m12.5 , depth=2 )
```


```{r, eval=FALSE}
## R code 12.19
inv_logit(coef(m12.5))
```


```{r, eval=FALSE}
## R code 12.20
( pk <- dordlogit( 1:7 , 0 , coef(m12.5) ) )
```


```{r, eval=FALSE}
## R code 12.21
sum( pk*(1:7) )
```


```{r, eval=FALSE}
## R code 12.22
( pk <- dordlogit( 1:7 , 0 , coef(m12.5)-0.5 ) )
```


```{r, eval=FALSE}
## R code 12.23
sum( pk*(1:7) )
```


```{r, eval=FALSE}
## R code 12.24
dat <- list(
  R = d$response,
  A = d$action,
  I = d$intention,
  C = d$contact )
m12.6 <- ulam(
  alist(
    R ~ dordlogit( phi , cutpoints ),
    phi <- bA*A + bC*C + BI*I ,
    BI <- bI + bIA*A + bIC*C ,
    c(bA,bI,bC,bIA,bIC) ~ dnorm( 0 , 0.5 ),
    cutpoints ~ dnorm( 0 , 1.5 )
  ) , data=dat , chains=4 , cores=4 )
precis( m12.6 )
```


```{r, eval=FALSE}
## R code 12.25
plot( precis(m12.6) , xlim=c(-1.4,0) )
```


```{r, eval=FALSE}
## R code 12.26
plot( NULL , type="n" , xlab="intention" , ylab="probability" ,
      xlim=c(0,1) , ylim=c(0,1) , xaxp=c(0,1,1) , yaxp=c(0,1,2) )

## R code 12.27
kA <- 0     # value for action
kC <- 0     # value for contact
kI <- 0:1   # values of intention to calculate over
pdat <- data.frame(A=kA,C=kC,I=kI)
phi <- link( m12.6 , data=pdat )$phi

## R code 12.28
post <- extract.samples( m12.6 )
for ( s in 1:50 ) {
  pk <- pordlogit( 1:6 , phi[s,] , post$cutpoints[s,] )
  for ( i in 1:6 ) lines( kI , pk[,i] , col=col.alpha("black",0.1) )
}
```


```{r, eval=FALSE}
## R code 12.29
kA <- 0     # value for action
kC <- 1     # value for contact
kI <- 0:1   # values of intention to calculate over
pdat <- data.frame(A=kA,C=kC,I=kI)
s <- sim( m12.6 , data=pdat )
simplehist( s , xlab="response" )
```


```{r, eval=FALSE}
## R code 12.30
library(rethinking)
data(Trolley)
d <- Trolley
levels(d$edu)
```


```{r, eval=FALSE}
## R code 12.31
edu_levels <- c( 6 , 1 , 8 , 4 , 7 , 2 , 5 , 3 )
d$edu_new <- edu_levels[ d$edu ]
```


```{r, eval=FALSE}
## R code 12.32
library(gtools)
set.seed(1805)
delta <- rdirichlet( 10 , alpha=rep(2,7) )
str(delta)
```


```{r, eval=FALSE}
## R code 12.33
h <- 3
plot( NULL , xlim=c(1,7) , ylim=c(0,0.4) , xlab="index" , ylab="probability" )
for ( i in 1:nrow(delta) ) lines( 1:7 , delta[i,] , type="b" ,
                                  pch=ifelse(i==h,16,1) , lwd=ifelse(i==h,4,1.5) ,
                                  col=ifelse(i==h,"black",col.alpha("black",0.7)) )
```

30 minutes!!
```{r, eval=FALSE}
## R code 12.34
dat <- list(
  R = d$response ,
  action = d$action,
  intention = d$intention,
  contact = d$contact,
  E = as.integer( d$edu_new ), # edu_new as an index
  alpha = rep(2.1,7) )           # delta prior

m12.5 <- ulam(
  alist(
    R ~ ordered_logistic( phi , kappa ),
    phi <- bE*sum( delta_j[1:E] ) + bA*action + bI*intention + bC*contact,
    kappa ~ normal( 0 , 1.5 ),
    c(bA,bI,bC,bE) ~ normal( 0 , 1 ),
    vector[8]: delta_j <<- append_row( 0 , delta ),
    simplex[7]: delta ~ dirichlet( alpha )
  ),
  data=dat , chains=3 , cores=3 )
```


```{r, eval=FALSE}
## R code 12.35
precis( m12.5 , depth=2 , omit="kappa" )
```


```{r, eval=FALSE}
## R code 12.36
delta_labels <- c("Elem","MidSch","SHS","HSG","SCol","Bach","Mast","Grad")
pairs( m12.5 , pars="delta" , labels=delta_labels )
```


```{r, eval=FALSE}
## R code 12.37
dat$edu_norm <- normalize( d$edu_new )
m12.6 <- ulam(
  alist(
    R ~ ordered_logistic( mu , cutpoints ),
    mu <- bE*edu_norm + bA*action + bI*intention + bC*contact,
    c(bA,bI,bC,bE) ~ normal( 0 , 1 ),
    cutpoints ~ normal( 0 , 1.5 )
  ), data=dat , chains=3 , cores=3 )
precis( m12.6 )
```