Exercises chapter 11.Rmd

# Exercises chapter 11

```{r}
library(rethinking)
data(chimpanzees)
d <- chimpanzees

d$treatment <- 1 + d$prosoc_left + 2*d$condition
xtabs( ~ treatment + prosoc_left + condition , d )

m11.1 <- quap( 
  alist(
    pulled_left ~ dbinom( 1 , p ) ,
    logit(p) <- a ,
    a ~ dnorm( 0 , 10 )
  ) , data=d )

set.seed(1999) 
prior <- extract.prior( m11.1 , n=1e4 )

p <- inv_logit( prior$a )
dens( p , adj=0.1 )


m11.2 <- quap(
  alist(
    pulled_left ~ dbinom( 1 , p ) ,
    logit(p) <- a + b[treatment] ,
    a ~ dnorm( 0 , 1.5 ),
    b[treatment] ~ dnorm( 0 , 10 )
  ) , data=d )
set.seed(1999)
prior <- extract.prior( m11.2 , n=1e4 )
p <- sapply( 1:4 , function(k) inv_logit( prior$a + prior$b[,k] ) )
dens( abs( p[,1] - p[,2] ) , adj=0.1 )


m11.3 <- quap( 
  alist(
    pulled_left ~ dbinom( 1 , p ) ,
    logit(p) <- a + b[treatment] ,
    a ~ dnorm( 0 , 1.5 ),
    b[treatment] ~ dnorm( 0 , 0.5 )
  ) , data=d )
set.seed(1999)
prior <- extract.prior( m11.3 , n=1e4 )
p <- sapply( 1:4 , function(k) inv_logit( prior$a + prior$b[,k] ) )
mean( abs( p[,1] - p[,2] ) )


# trimmed data list 
dat_list <- list(
  pulled_left = d$pulled_left,
  actor = d$actor,
  treatment = as.integer(d$treatment) )

m11.4 <- ulam(
  alist(
    pulled_left ~ dbinom( 1 , p ) ,
    logit(p) <- a[actor] + b[treatment] ,
    a[actor] ~ dnorm( 0 , 1.5 ),
    b[treatment] ~ dnorm( 0 , 0.5 )
  ) , data=dat_list , chains=4 , log_lik=TRUE )
precis( m11.4 , depth=2 )

post <- extract.samples(m11.4)
p_left <- inv_logit( post$a )
plot( precis( as.data.frame(p_left) ) , xlim=c(0,1) )

labs <- c("R/N","L/N","R/P","L/P") 
plot( precis( m11.4 , depth=2 , pars="b" ) , labels=labs )

diffs <- list( 
  db13 = post$b[,1] - post$b[,3],
  db24 = post$b[,2] - post$b[,4] )
plot( precis(diffs) )

pl <- by( d$pulled_left , list( d$actor , d$treatment ) , mean )
pl[1,]

plot( NULL , xlim=c(1,28) , ylim=c(0,1) , xlab="" ,
  ylab="proportion left lever" , xaxt="n" , yaxt="n" )
axis( 2 , at=c(0,0.5,1) , labels=c(0,0.5,1) )
abline( h=0.5 , lty=2 )
for ( j in 1:7 ) abline( v=(j-1)*4+4.5 , lwd=0.5 )
for ( j in 1:7 ) text( (j-1)*4+2.5 , 1.1 , concat("actor ",j) , xpd=TRUE )
for ( j in (1:7)[-2] ) {
  lines( (j-1)*4+c(1,3) , pl[j,c(1,3)] , lwd=2 , col=rangi2 )
  lines( (j-1)*4+c(2,4) , pl[j,c(2,4)] , lwd=2 , col=rangi2 )
}
points( 1:28 , t(pl) , pch=16 , col="white" , cex=1.7 )
points( 1:28 , t(pl) , pch=c(1,1,16,16) , col=rangi2 , lwd=2 )
yoff <- 0.01
text( 1 , pl[1,1]-yoff , "R/N" , pos=1 , cex=0.8 )
text( 2 , pl[1,2]+yoff , "L/N" , pos=3 , cex=0.8 )
text( 3 , pl[1,3]-yoff , "R/P" , pos=1 , cex=0.8 )
text( 4 , pl[1,4]+yoff , "L/P" , pos=3 , cex=0.8 )
mtext( "observed proportions\n" )

dat <- list( actor=rep(1:7,each=4) , treatment=rep(1:4,times=7) )
p_post <- link( m11.4 , data=dat )
p_mu <- apply( p_post , 2 , mean )
p_ci <- apply( p_post , 2 , PI )

d$side <- d$prosoc_left + 1 # right 1, left 2
d$cond <- d$condition + 1 # no partner 1, partner 2

dat_list2 <- list(
  pulled_left = d$pulled_left,
  actor = d$actor,
  side = d$side,
  cond = d$cond )

m11.5 <- ulam(
  alist(
    pulled_left ~ dbinom( 1 , p ) ,
    logit(p) <- a[actor] + bs[side] + bc[cond] ,
    a[actor] ~ dnorm( 0 , 1.5 ),
    bs[side] ~ dnorm( 0 , 0.5 ),
    bc[cond] ~ dnorm( 0 , 0.5 )
  ) , data=dat_list2 , chains=4 , log_lik=TRUE )

compare( m11.5 , m11.4 , func=PSIS )

```


```{r}
data(chimpanzees) 
d <- chimpanzees
d$treatment <- 1 + d$prosoc_left + 2*d$condition
d$side <- d$prosoc_left + 1 # right 1, left 2
d$cond <- d$condition + 1 # no partner 1, partner 2

d_aggregated <- aggregate(
  d$pulled_left ,
  list( treatment=d$treatment , actor=d$actor ,
    side=d$side , cond=d$cond ) ,
  sum )
colnames(d_aggregated)[5] <- "left_pulls"

dat <- with( d_aggregated , list(
  left_pulls = left_pulls,
  treatment = treatment,
  actor = actor,
  side = side,
  cond = cond ) )

m11.6 <- ulam(
  alist(
    left_pulls ~ dbinom( 18 , p ) ,
    logit(p) <- a[actor] + b[treatment] ,
    a[actor] ~ dnorm( 0 , 1.5 ) ,
    b[treatment] ~ dnorm( 0 , 0.5 )
  ) , data=dat , chains=4 , log_lik=TRUE )

compare( m11.6 , m11.4 , func=PSIS )

# deviance of aggregated 6-in-9 
-2*dbinom(6,9,0.2,log=TRUE)
# deviance of dis-aggregated
-2*sum(dbern(c(1,1,1,1,1,1,0,0,0),0.2,log=TRUE))

```
```{r}
library(rethinking)
data(UCBadmit)
d <- UCBadmit

dat_list <- list( 
  admit = d$admit,
  applications = d$applications,
  gid = ifelse( d$applicant.gender=="male" , 1 , 2 )
)
m11.7 <- ulam(
  alist(
    admit ~ dbinom( applications , p ) ,
    logit(p) <- a[gid] ,
    a[gid] ~ dnorm( 0 , 1.5 )
  ) , data=dat_list , chains=4 )
precis( m11.7 , depth=2 )

post <- extract.samples(m11.7) 
diff_a <- post$a[,1] - post$a[,2]
diff_p <- inv_logit(post$a[,1]) - inv_logit(post$a[,2])
precis( list( diff_a=diff_a , diff_p=diff_p ) )

postcheck( m11.7 )
# draw lines connecting points from same dept
for ( i in 1:6 ) {
  x <- 1 + 2*(i-1)
  y1 <- d$admit[x]/d$applications[x]
  y2 <- d$admit[x+1]/d$applications[x+1]
  lines( c(x,x+1) , c(y1,y2) , col=rangi2 , lwd=2 )
  text( x+0.5 , (y1+y2)/2 + 0.05 , d$dept[x] , cex=0.8 , col=rangi2 )
}

dat_list$dept_id <- rep(1:6,each=2) 
m11.8 <- ulam(
  alist(
    admit ~ dbinom( applications , p ) ,
    logit(p) <- a[gid] + delta[dept_id] ,
    a[gid] ~ dnorm( 0 , 1.5 ) ,
    delta[dept_id] ~ dnorm( 0 , 1.5 )
  ) , data=dat_list , chains=4 , iter=4000 )
precis( m11.8 , depth=2 )

post <- extract.samples(m11.8) 
diff_a <- post$a[,1] - post$a[,2]
diff_p <- inv_logit(post$a[,1]) - inv_logit(post$a[,2])
precis( list( diff_a=diff_a , diff_p=diff_p ) )

pg <- with( dat_list , sapply( 1:6 , function(k)
applications[dept_id==k]/sum(applications[dept_id==k]) ) )
rownames(pg) <- c("male","female")
colnames(pg) <- unique(d$dept)
round( pg , 2 )

```
11.2 Poisson regression

```{r}
library(rethinking) 
data(Kline)
d <- Kline
d
d$P <- scale( log(d$population) )
d$contact_id <- ifelse( d$contact=="high" , 2 , 1 )

curve( dlnorm( x , 3 , 0.5 ) , from=0 , to=100 , n=200 )

N <- 100 
a <- rnorm( N , 3 , 0.5 )
b <- rnorm( N , 0 , 0.2 )
plot( NULL , xlim=c(-2,2) , ylim=c(0,100) )
for ( i in 1:N ) curve( exp( a[i] + b[i]*x ) , add=TRUE , col=grau() )

x_seq <- seq( from=log(100) , to=log(200000) , length.out=100 )
lambda <- sapply( x_seq , function(x) exp( a + b*x ) )
plot( NULL , xlim=range(x_seq) , ylim=c(0,500) , xlab="log population" ,
ylab="total tools" )
for ( i in 1:N ) lines( x_seq , lambda[i,] , col=grau() , lwd=1.5 )

plot( NULL , xlim=range(exp(x_seq)) , ylim=c(0,500) , xlab="population" ,
ylab="total tools" )
for ( i in 1:N ) lines( exp(x_seq) , lambda[i,] , col=grau() , lwd=1.5 )


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

# intercept only
m11.9 <- ulam(
  alist(
    T ~ dpois( lambda ),
    log(lambda) <- a,
    a ~ dnorm( 3 , 0.5 )
  ), data=dat , chains=4 , log_lik=TRUE )

# interaction model
m11.10 <- ulam(
  alist(
    T ~ dpois( lambda ),
    log(lambda) <- a[cid] + b[cid]*P,
    a[cid] ~ dnorm( 3 , 0.5 ),
    b[cid] ~ dnorm( 0 , 0.2 )
  ), data=dat , chains=4 , log_lik=TRUE )

compare( m11.9 , m11.10 , func=PSIS )


k <- PSIS( m11.10 , pointwise=TRUE )$k 
plot( dat$P , dat$T , xlab="log population (std)" , ylab="total tools" ,
  col=rangi2 , pch=ifelse( dat$cid==1 , 1 , 16 ) , lwd=2 ,
  ylim=c(0,75) , cex=1+normalize(k) )

# set up the horizontal axis values to compute predictions at
ns <- 100
P_seq <- seq( from=-1.4 , to=3 , length.out=ns )

# predictions for cid=1 (low contact)
lambda <- link( m11.10 , data=data.frame( P=P_seq , cid=1 ) )
lmu <- apply( lambda , 2 , mean )
lci <- apply( lambda , 2 , PI )
lines( P_seq , lmu , lty=2 , lwd=1.5 )
shade( lci , P_seq , xpd=TRUE )

# predictions for cid=2 (high contact)
lambda <- link( m11.10 , data=data.frame( P=P_seq , cid=2 ) )
lmu <- apply( lambda , 2 , mean )
lci <- apply( lambda , 2 , PI )
lines( P_seq , lmu , lty=1 , lwd=1.5 )
shade( lci , P_seq , xpd=TRUE )


plot( d$population , d$total_tools , xlab="population" , ylab="total tools" , 
  col=rangi2 , pch=ifelse( dat$cid==1 , 1 , 16 ) , lwd=2 ,
  ylim=c(0,75) , cex=1+normalize(k) )

ns <- 100
P_seq <- seq( from=-5 , to=3 , length.out=ns )
# 1.53 is sd of log(population)
# 9 is mean of log(population)
pop_seq <- exp( P_seq*1.53 + 9 )

lambda <- link( m11.10 , data=data.frame( P=P_seq , cid=1 ) )
lmu <- apply( lambda , 2 , mean )
lci <- apply( lambda , 2 , PI )
lines( pop_seq , lmu , lty=2 , lwd=1.5 )
shade( lci , pop_seq , xpd=TRUE )

lambda <- link( m11.10 , data=data.frame( P=P_seq , cid=2 ) )
lmu <- apply( lambda , 2 , mean )
lci <- apply( lambda , 2 , PI )
lines( pop_seq , lmu , lty=1 , lwd=1.5 )
shade( lci , pop_seq , xpd=TRUE )


num_days <- 30 
y <- rpois( num_days , 1.5 )
num_weeks <- 4
y_new <- rpois( num_weeks , 0.5*7 )
y_all <- c( y , y_new )
exposure <- c( rep(1,30) , rep(7,4) )
monastery <- c( rep(0,30) , rep(1,4) )
d <- data.frame( y=y_all , days=exposure , monastery=monastery )

# compute the offset
d$log_days <- log( d$days )

# fit the model
m11.12 <- quap(
  alist(
    y ~ dpois( lambda ),
    log(lambda) <- log_days + a + b*monastery,
    a ~ dnorm( 0 , 1 ),
    b ~ dnorm( 0 , 1 )
  ), data=d )

post <- extract.samples( m11.12 )
lambda_old <- exp( post$a )
lambda_new <- exp( post$a + post$b )
precis( data.frame( lambda_old , lambda_new ) )

```

```{r}
# simulate career choices among 500 individuals
N <- 500             # number of individuals
income <- c(1,2,5)   # expected income of each career
score <- 0.5*income  # scores for each career, based on income
# next line converts scores to probabilities
p <- softmax(score[1],score[2],score[3])

# now simulate choice
# outcome career holds event type values, not counts
career <- rep(NA,N) # empty vector of choices for each individual
# sample chosen career for each individual
set.seed(34302)
for ( i in 1:N ) career[i] <- sample( 1:3 , size=1 , prob=p )

code_m11.13 <- " 
data{
    int N; // number of individuals
    int K; // number of possible careers
    int career[N]; // outcome
    vector[K] career_income;
}
parameters{
    vector[K-1] a; // intercepts
    real<lower=0> b; // association of income with choice
}
model{
    vector[K] p;
    vector[K] s;
    a ~ normal( 0 , 1 );
    b ~ normal( 0 , 0.5 );
    s[1] = a[1] + b*career_income[1];
    s[2] = a[2] + b*career_income[2];
    s[3] = 0; // pivot
    p = softmax( s );
    career ~ categorical( p );
}
"
dat_list <- list( N=N , K=3 , career=career , career_income=income ) 
m11.13 <- stan( model_code=code_m11.13 , data=dat_list , chains=4, iter = 4000 )
precis( m11.13 , 2 )


post <- extract.samples( m11.13 )

# set up logit scores
s1 <- with( post , a[,1] + b*income[1] )
s2_orig <- with( post , a[,2] + b*income[2] )
s2_new <- with( post , a[,2] + b*income[2]*2 )

# compute probabilities for original and counterfactual
p_orig <- sapply( 1:length(post$b) , function(i)
    softmax( c(s1[i],s2_orig[i],0) ) )
p_new <- sapply( 1:length(post$b) , function(i)
    softmax( c(s1[i],s2_new[i],0) ) )

# summarize
p_diff <- p_new[2,] - p_orig[2,]
precis( p_diff )


N <- 500
# simulate family incomes for each individual
family_income <- runif(N)
# assign a unique coefficient for each type of event
b <- c(-2,0,2)
career <- rep(NA,N) # empty vector of choices for each individual
for ( i in 1:N ) {
    score <- 0.5*(1:3) + b*family_income[i]
    p <- softmax(score[1],score[2],score[3])
    career[i] <- sample( 1:3 , size=1 , prob=p )
}

code_m11.14 <- "
data{
    int N; // number of observations
    int K; // number of outcome values
    int career[N]; // outcome
    real family_income[N];
}
parameters{
    vector[K-1] a; // intercepts
    vector[K-1] b; // coefficients on family income
}
model{
    vector[K] p;
    vector[K] s;
    a ~ normal(0,1.5);
    b ~ normal(0,1);
    for ( i in 1:N ) {
        for ( j in 1:(K-1) ) s[j] = a[j] + b[j]*family_income[i];
        s[K] = 0; // the pivot
        p = softmax( s );
        career[i] ~ categorical( p );
    }
}
"

dat_list <- list( N=N , K=3 , career=career , family_income=family_income )
m11.14 <- stan( model_code=code_m11.14 , data=dat_list , chains=4, iter = 4000 )
precis( m11.14 , 2 )

```


```{r}
library(rethinking)
data(UCBadmit)
d <- UCBadmit

# binomial model of overall admission probability
m_binom <- quap(
    alist(
        admit ~ dbinom(applications,p),
        logit(p) <- a,
        a ~ dnorm( 0 , 1.5 )
    ), data=d )

# Poisson model of overall admission rate and rejection rate
# 'reject' is a reserved word in Stan, cannot use as variable name
dat <- list( admit=d$admit , rej=d$reject )
m_pois <- ulam(
    alist(
        admit ~ dpois(lambda1),
        rej ~ dpois(lambda2),
        log(lambda1) <- a1,
        log(lambda2) <- a2,
        c(a1,a2) ~ dnorm(0,1.5)
    ), data=dat , chains=3 , cores=3 )

inv_logit(coef(m_binom))
k <- coef(m_pois) 
a1 <- k['a1']; a2 <- k['a2']
exp(a1)/(exp(a1)+exp(a2))

```



################ Exercises ################ 
### easy

11E1. If an event has probability 0.35, what are the log-odds of this event?
```{r}
log(0.35/0.65)
```
11E2. If an event has log-odds 3.2, what is the probability of this event?
```{r}
1-(1/exp(3.2))
```
11E3. Suppose that a coefficient in a logistic regression has value 1.7. What does this imply about
the proportional change in odds of the outcome?
```{r}
exp(1.7)
```
11E4. Why do Poisson regressions sometimes require the use of an offset? Provide an example.


### hard

11H1. Use WAIC or PSIS to compare the chimpanzee model that includes a unique intercept for
each actor, m11.4 (page 330), to the simpler models fit in the same section. Interpret the results.

```{r}
library(rethinking)
data(chimpanzees)
d <- chimpanzees
d$treatment <- 1 + d$prosoc_left + 2*d$condition

dat_list <- list(pulled_left = d$pulled_left)

m11.1 <- ulam( 
    alist(
        pulled_left ~ dbinom( 1 , p ) ,
        logit(p) <- a ,
        a ~ dnorm( 0 , 1.5 )
    ) , data=dat_list, chains=4 , cores=4, iter = 4000, log_lik=TRUE  )

prior <- extract.prior( m11.1 , n=1e4 )
p <- inv_logit( prior$a )
dat_list <- list(pulled_left = d$pulled_left, treatment = d$treatment)


m11.3 <- ulam(
    alist(
        pulled_left ~ dbinom( 1 , p ) ,
        logit(p) <- a + b[treatment] ,
        a ~ dnorm( 0 , 1.5 ),
        b[treatment] ~ dnorm( 0 , 0.5 )
    ) , data=dat_list, chains=4 ,cores=4, iter = 4000, log_lik=TRUE  )

# trimmed data list 
dat_list <- list(
pulled_left = d$pulled_left,
actor = d$actor,
treatment = as.integer(d$treatment) )

m11.4 <- ulam(
    alist(
        pulled_left ~ dbinom( 1 , p ) ,
        logit(p) <- a[actor] + b[treatment] ,
        a[actor] ~ dnorm( 0 , 1.5 ),
        b[treatment] ~ dnorm( 0 , 0.5 )
    ) , data=dat_list , chains=4 , log_lik=TRUE )

precis( m11.4 , depth=2 )

compare(m11.1 , m11.3 , m11.4 , func=PSIS )

```
Since different model classes are used, the models are not really comparable in this way.


11.H2

```{r}
library(rethinking)
library(MASS);data(eagles)
dd <- eagles
dd$P <- as.integer(ifelse(dd$P == 'L',2,1))
dd$V <- as.integer(ifelse(dd$V == 'L',2,1))
dd$A <- as.integer(ifelse(dd$A == 'A',2,1))


H2 <- ulam(
    alist(
        y ~ dbinom( n , p ) ,
        logit(p) <- a + bp*P + bv*V + ba*A ,
        a ~ dnorm( 0 , 1.5 ),
        bp ~ dnorm( 0 , 0.5 ),
        bv ~ dnorm( 0 , 0.5 ),
        ba ~ dnorm( 0 , 0.5 )
    ) , data=dd , chains=4 , log_lik=TRUE, iter = 4000 )

H2q <- quap(
    alist(
        y ~ dbinom( n , p ) ,
        logit(p) <- a + bp*P + bv*V + ba*A ,
        a ~ dnorm( 0 , 1.5 ),
        bp ~ dnorm( 0 , 0.5 ),
        bv ~ dnorm( 0 , 0.5 ),
        ba ~ dnorm( 0 , 0.5 )
    ) , data=dd )

precis(H2, depth = 2)
precis(H2q, depth = 2)

```
a) The quap model gives almost identical outcomes as the stan model and therefore seems fine



b)
Body size of the pirate is positively correlated with success, and so is age/maturity, whereas body size of the victim is negatively correlated with success.

```{r}
postcheck(H2)
dd$scenario <- c(1:8)

plot( dd$n , dd$y , xlab="Total attempts" , ylab="Successful attempts" ,
col=rangi2 , pch=dd$scenario , lwd=2 )

# set up the horizontal axis values to compute predictions at
n_seq <- seq( from=1, to=40 , length.out=40 )

# predictions for P=1 small pirate
lambda <- link( H2 , data=data.frame( n=n_seq , P=1, V=2, A=2 ) )
lmu <- apply( lambda , 2 , mean )
lci <- apply( lambda , 2 , PI )
lines( n_seq , lmu , lty=2 , lwd=1.5 )
shade( lci , n_seq , xpd=TRUE )

# predictions for P=2 large pirate
lambda <- link( H2 , data=data.frame( n=n_seq , P=2, V=2, A=2 ) )
lmu <- apply( lambda , 2 , mean )
lci <- apply( lambda , 2 , PI )
lines( n_seq , lmu , lty=1 , lwd=1.5 )
shade( lci , n_seq , xpd=TRUE )



```
c)
```{r}
H2c <- ulam(
    alist(
        y ~ dbinom( n , p ) ,
        logit(p) <- a + bp*P + bv*V + ba*A + bpa*P*A,
        a ~ dnorm( 0 , 1.5 ),
        bp ~ dnorm( 0 , 0.5 ),
        bv ~ dnorm( 0 , 0.5 ),
        ba ~ dnorm( 0 , 0.5 ),
        bpa ~ dnorm( 0 , 0.5 )
    ) , data=dd , chains=4 , log_lik=TRUE, iter = 4000 )

precis(H2c, depth = 2)

compare(H2, H2c)
coeftab(H2, H2c)

```
The interaction does not seem to significantly improve the predictive power of the model.




11H3
(a) Model the relationship between density and percent cover, using a log-link (same as the example
in the book and lecture). Use weakly informative priors of your choosing. Check the quadratic
approximation again, by comparing quap to ulam. Then plot the expected counts and their 89% interval
against percent cover. In which ways does the model do a good job? A bad job?
```{r}
library(rethinking)
data(salamanders)
d <- salamanders

dat = list(SALAMAN = d$SALAMAN,PCTCOVER = d$PCTCOVER)
           
H11.3 <- ulam(
  alist(
    SALAMAN ~ dpois( lambda ),
    log(lambda) <- a + b*PCTCOVER,
    a ~ dnorm( 3 , 0.5 ),
    b ~ dnorm( 0 , 0.2 )
  ), data=dat , chains=4 , cores = 4, iter = 4000, log_lik=TRUE )

           
H11.3q <- quap(
  alist(
    SALAMAN ~ dpois( lambda ),
    log(lambda) <- a + b*PCTCOVER,
    a ~ dnorm( 3 , 0.5 ),
    b ~ dnorm( 0 , 0.2 )
  ), data=dat )

precis(H11.3, depth = 2)
precis(H11.3q, depth = 2)


plot( d$PCTCOVER , d$SALAMAN , xlab="Percent of ground cover" , ylab="Number of Salamanders" ,
col=rangi2 , lwd=2 )

# set up the horizontal axis values to compute predictions at
pct_seq <- seq( from=1, to=100 , length.out=100 )

# predictions for salamander count base on ground cover
lambda <- link( H11.3 , data=data.frame( PCTCOVER=pct_seq ) )
lmu <- apply( lambda , 2 , mean )
lci <- apply( lambda , 2 , PI )
lines( pct_seq , lmu , lty=2 , lwd=1.5 )
shade( lci , pct_seq , xpd=TRUE )

```
The ulam & quap model seem to make the same predictions.
The model does a good job at predicting the salamander counts in the lower ranges of the ground cover, but at the higher percentages there is just too much variation and the prediction seems more certain than is warranted.


(b) Can you improve the model by using the other predictor, FORESTAGE? Try any models you
think useful. Can you explain why FORESTAGE helps or does not help with prediction?
```{r}
dat = list(SALAMAN = d$SALAMAN, PCTCOVER = d$PCTCOVER, FORESTAGE = log(d$FORESTAGE))
dat$FORESTAGE[30] <- 0

H11.3b1 <- ulam(
  alist(
    SALAMAN ~ dpois( lambda ),
    log(lambda) <- a + b*FORESTAGE,
    a ~ dnorm( 3 , 0.5 ),
    b ~ dnorm( 0 , 0.2 )
  ), data=dat , chains=4 , cores = 4, iter = 4000, log_lik=TRUE )

precis(H11.3b1, depth = 2)

dat = list(SALAMAN = d$SALAMAN, PCTCOVER = d$PCTCOVER, FORESTAGE = log(d$FORESTAGE))
dat$FORESTAGE[30] <- 0

H11.3b2 <- ulam(
  alist(
    SALAMAN ~ dpois( lambda ),
    log(lambda) <- a + b*FORESTAGE + c*PCTCOVER,
    a ~ dnorm( 3 , 0.5 ),
    b ~ dnorm( 0 , 0.2 ),
    c ~ dnorm( 0 , 0.2 )
  ), data=dat , chains=4 , cores = 4, iter = 4000, log_lik=TRUE )

precis(H11.3b2, depth = 2)

H11.3b3 <- ulam(
  alist(
    SALAMAN ~ dpois( lambda ),
    log(lambda) <- a + b*FORESTAGE + c*PCTCOVER + bc*FORESTAGE*PCTCOVER,
    a ~ dnorm( 3 , 0.1 ),
    b ~ dnorm( 0 , 0.5 ),
    c ~ dnorm( 0 , 0.5 ),
    bc ~ dnorm( 0 , 0.5 )
  ), data=dat , chains=4 , cores = 4, iter = 4000, log_lik=TRUE )

precis(H11.3b3, depth = 2)

compare(H11.3b1, H11.3b2, H11.3b3)


plot( dat$FORESTAGE , dat$SALAMAN , xlab="Log forest age" , ylab="Number of Salamanders" ,
col=rangi2 , lwd=2 )

# set up the horizontal axis values to compute predictions at
log_seq <- seq( from=0, to=7 , length.out=100 )

# predictions for salamander count base on ground cover
lambda <- link( H11.3b1 , data=data.frame( FORESTAGE=log_seq ) )
lmu <- apply( lambda , 2 , mean )
lci <- apply( lambda , 2 , PI )
lines( log_seq , lmu , lty=2 , lwd=1.5 )
shade( lci , log_seq , xpd=TRUE )
```


11H6. The data in data(Primates301) are 301 primate species and associated measures. In this
problem, you will consider how brain size is associated with social learning. There are three parts.
(a) Model the number of observations of social_learning for each species as a function of the
log brain size. Use a Poisson distribution for the social_learning outcome variable. Interpret the
resulting posterior. 

```{r}
library(rethinking)
data(Primates301)
d <- Primates301

# remove rows with NA in brain or social_learning column
dd <- d[complete.cases(d$social_learning),]
dd <- dd[complete.cases(dd$brain),]

dat <- list(socl = dd$social_learning, logBS = log(dd$brain))

# Check priors
N <- 100 
a <- rnorm( N , 0 ,5 )
b <- rnorm( N , 0 , 1 )
plot( NULL , xlim=c(0,4) , ylim=c(0,100) )
for ( i in 1:N ) curve( exp( a[i] + b[i]*x ) , add=TRUE , col=grau() )

# we want the priors to be quite vague because there can be a large variation and skew in the amount of observed social learning

# Build model
H11.6 <- ulam(
  alist(
    socl ~ dpois( lambda ),
    log(lambda) <- a + bB*logBS ,
    a ~ dnorm( 0 , 5 ),
    bB ~ dnorm( 0 , 1 )
  ), data=dat , chains=4 , cores = 4, iter = 4000, log_lik=TRUE )

precis(H11.6, depth = 2)

# Plot data & posterior predictions
plot( dat$logBS , dat$socl , xlab="log brain size" , ylab="social learning observations" , col=rangi2 , lwd=2 )

# set up the horizontal axis values to compute predictions at
ns <- 100
bs_seq <- seq( from=0 , to=6.5 , length.out=ns )

# predictions 
lambda <- link( H11.6 , data=data.frame( logBS=bs_seq ) )
lmu <- apply( lambda , 2 , mean )
lci <- apply( lambda , 2 , PI )
lines( bs_seq , lmu , lty=2 , lwd=1.5 )
shade( lci , bs_seq , xpd=TRUE )

```
The posterior fits the data quite well however, this is strongly affected by most of the social learning observations being zero. 
This makes it very hard for the model to make good predictions for the species where social learning is actually observed.
At larger brain sizes the model is too confident in its predictions.

(b) Some species are studied much more than others. So the number of reported
instances of social_learning could be a product of research effort. Use the research_effort
variable, specifically its logarithm, as an additional predictor variable. Interpret the coefficient for log
research_effort. How does this model differ from the previous one? 

```{r}
dd <- dd[complete.cases(dd$research_effort),]
dat <- list(socl = dd$social_learning, logBS = log(dd$brain), logRE = log(dd$research_effort))

# Build model
H11.6b <- ulam(
  alist(
    socl ~ dpois( lambda ),
    log(lambda) <- a + bB*logBS + bR*logRE,
    a ~ dnorm( 0 , 5 ),
    bB ~ dnorm( 0 , 1 ),
    bR ~ dnorm( 0 , 1 )
  ), data=dat , chains=4 , cores = 4, iter = 4000, log_lik=TRUE )

precis(H11.6b, depth = 2)
compare(H11.6, H11.6b, func = PSIS)
coeftab(H11.6, H11.6b)


```
There seems to be a strong positive correlation between research effort and the amount of social learning observed. This does not sound strange to me, the more you look for something, the more likely you are to observe it. 
There is still a positive correlation between brain size and observed social learning as well, however, this accounts for a much smaller part of the predicted observed social learning than in the previous, simpler model.

(c) Draw a DAG to represent
how you think the variables social_learning, brain, and research_effort interact. Justify the
DAG with the measured associations in the two models above (and any other models you used).

```{r}
# check whether there is an interaction
H11.6c <- ulam(
  alist(
    socl ~ dpois( lambda ),
    log(lambda) <- a + bB*logBS + bR*logRE + bBR*logBS*logRE,
    a ~ dnorm( 0 , 5 ),
    bB ~ dnorm( 0 , 1 ),
    bR ~ dnorm( 0 , 1 ),
    bBR ~ dnorm( 0 , 1 )
  ), data=dat , chains=6 , cores = 6, iter = 6000, log_lik=TRUE )

traceplot(H11.6c)
precis(H11.6c, depth = 2)
compare(H11.6, H11.6b, H11.6c, func = PSIS)
coeftab(H11.6, H11.6b, H11.6c)

# Plot data brain & research effort
plot( dat$logBS , dat$logRE , xlab="log brain size" , ylab="log research effort" , col=rangi2 , lwd=2 )
abline(lm(dat$logRE~dat$logBS),col='green')

plot( dd$brain , dd$research_effort , xlab="brain size" , ylab="research effort" , col=rangi2 , lwd=2 )
abline(lm(dd$brain~dd$research_effort),col='green')

# Plot posterior predictions as research effort is kept constant
plot( dat$logBS , dat$socl , xlab="log brain size" , ylab="social learning observations" , col=rangi2 , lwd=2 )

# set up the horizontal axis values to compute predictions
ns <- 100
bs_seq <- seq( from=0 , to=6.5 , length.out=ns )

# predictions 
lambda <- link( H11.6b , data=data.frame( logBS=bs_seq, logRE=mean(dat$logRE) ) )
lmu <- apply( lambda , 2 , mean )
lci <- apply( lambda , 2 , PI )
lines( bs_seq , lmu , lty=2 , lwd=1.5 )
shade( lci , bs_seq , xpd=TRUE )

library(dagitty)
g <- dagitty('dag {
BS [pos="0,2"]
RE [pos="0,0"]
SL [pos="1,1"]

BS -> RE -> SL
BS -> SL
}')
plot(g)

```
The model including an interaction between brain size and research effort performs worse than the model without according to the PSIS values, and the value given for the interaction predictor is around zero.
However, looking at the data there seems to be a correlation between brain size and research effort, and when creating a counterfactual plot where research effort is kept constant, brain size itself lose basically all its predictive power.
Therefore, I believe that the social learning observations are driven to a large part by research effort, and that research effort is driven by brain size. Brain size itself also plays a role, but seemingly less than research effort. Predictions are complicated by the fact that there are only a handful of primate species with a brain size above 300, and all of these are heavily studied.