Exercises chapter 7.Rmd

Exercises chapter 7

```{r}
library(rethinking)
sppnames <- c( "afarensis","africanus","habilis","boisei",
               "rudolfensis","ergaster","sapiens")
brainvolcc <- c( 438 , 452 , 612, 521, 752, 871, 1350 )
masskg <- c( 37.0 , 35.5 , 34.5 , 41.5 , 55.5 , 61.0 , 53.5 )
d <- data.frame( species=sppnames , brain=brainvolcc , mass=masskg )

d$mass_std <- (d$mass - mean(d$mass))/sd(d$mass) 
d$brain_std <- d$brain / max(d$brain)

m7.1 <- quap(
  alist(
    brain_std ~ dnorm( mu , exp(log_sigma) ),
    mu <- a + b*mass_std,
    a ~ dnorm( 0.5 , 1 ),
    b ~ dnorm( 0 , 10 ),
    log_sigma ~ dnorm( 0 , 1 )
  ), data=d )

set.seed(12) 
s <- sim( m7.1 )
r <- apply(s,2,mean) - d$brain_std
resid_var <- var2(r)
outcome_var <- var2( d$brain_std )
1 - resid_var/outcome_var

R2_is_bad <- function( quap_fit ) {
s <- sim( quap_fit , refresh=0 )
r <- apply(s,2,mean) - d$brain_std
1 - var2(r)/var2(d$brain_std)
}


m7.2 <- quap(
  alist(
    brain_std ~ dnorm( mu , exp(log_sigma) ),
    mu <- a + b[1]*mass_std + b[2]*mass_std^2,
    a ~ dnorm( 0.5 , 1 ),
    b ~ dnorm( 0 , 10 ),
    log_sigma ~ dnorm( 0 , 1 )
  ), data=d , start=list(b=rep(0,2)) )

m7.3 <- quap(
  alist(
    brain_std ~ dnorm( mu , exp(log_sigma) ),
    mu <- a + b[1]*mass_std + b[2]*mass_std^2 +
    b[3]*mass_std^3,
    a ~ dnorm( 0.5 , 1 ),
    b ~ dnorm( 0 , 10 ),
    log_sigma ~ dnorm( 0 , 1 )
  ), data=d , start=list(b=rep(0,3)) )

m7.4 <- quap(
  alist(
    brain_std ~ dnorm( mu , exp(log_sigma) ),
    mu <- a + b[1]*mass_std + b[2]*mass_std^2 +
    b[3]*mass_std^3 + b[4]*mass_std^4,
    a ~ dnorm( 0.5 , 1 ),
    b ~ dnorm( 0 , 10 ),
    log_sigma ~ dnorm( 0 , 1 )
  ), data=d , start=list(b=rep(0,4)) )

m7.5 <- quap(
  alist(
    brain_std ~ dnorm( mu , exp(log_sigma) ),
    mu <- a + b[1]*mass_std + b[2]*mass_std^2 +
    b[3]*mass_std^3 + b[4]*mass_std^4 +
    b[5]*mass_std^5,
    a ~ dnorm( 0.5 , 1 ),
    b ~ dnorm( 0 , 10 ),
    log_sigma ~ dnorm( 0 , 1 )
  ), data=d , start=list(b=rep(0,5)) )

m7.6 <- quap( 
  alist(
    brain_std ~ dnorm( mu , 0.001 ),
    mu <- a + b[1]*mass_std + b[2]*mass_std^2 +
    b[3]*mass_std^3 + b[4]*mass_std^4 +
    b[5]*mass_std^5 + b[6]*mass_std^6,
    a ~ dnorm( 0.5 , 1 ),
    b ~ dnorm( 0 , 10 )
  ), data=d , start=list(b=rep(0,6)) )


post <- extract.samples(m7.1) 
mass_seq <- seq( from=min(d$mass_std) , to=max(d$mass_std) , length.out=100 )
l <- link( m7.1 , data=list( mass_std=mass_seq ) )
mu <- apply( l , 2 , mean )
ci <- apply( l , 2 , PI )
plot( brain_std ~ mass_std , data=d )
lines( mass_seq , mu )
shade( ci , mass_seq )

```

## 7.2

```{r}
p <- c( 0.3 , 0.7 )
-sum( p*log(p) )

set.seed(1)
lppd( m7.1 , n=1e4 )

set.seed(1) 
sapply( list(m7.1,m7.2,m7.3,m7.4,m7.5,m7.6) , function(m) sum(lppd(m)) )

```

7.5

```{r}
library(rethinking)
data(WaffleDivorce)
d <- WaffleDivorce
d$A <- standardize( d$MedianAgeMarriage )
d$D <- standardize( d$Divorce )
d$M <- standardize( d$Marriage )

m5.1 <- quap(
  alist(
    D ~ dnorm( mu , sigma ) ,
    mu <- a + bA * A ,
    a ~ dnorm( 0 , 0.2 ) ,
    bA ~ dnorm( 0 , 0.5 ) ,
    sigma ~ dexp( 1 )
  ) , data = d )

m5.2 <- quap(
  alist(
    D ~ dnorm( mu , sigma ) ,
    mu <- a + bM * M ,
    a ~ dnorm( 0 , 0.2 ) ,
    bM ~ dnorm( 0 , 0.5 ) ,
    sigma ~ dexp( 1 )
  ) , data = d )

m5.3 <- quap(
  alist(
    D ~ dnorm( mu , sigma ) ,
    mu <- a + bM*M + bA*A ,
    a ~ dnorm( 0 , 0.2 ) ,
    bM ~ dnorm( 0 , 0.5 ) ,
    bA ~ dnorm( 0 , 0.5 ) ,
    sigma ~ dexp( 1 )
  ) , data = d )

set.seed(24071847) 
compare( m5.1 , m5.2 , m5.3 , func=PSIS )

set.seed(24071847)
PSIS_m5.3 <- PSIS(m5.3,pointwise=TRUE)
set.seed(24071847)
WAIC_m5.3 <- WAIC(m5.3,pointwise=TRUE)
plot( PSIS_m5.3$k , WAIC_m5.3$penalty , xlab="PSIS Pareto k" ,
ylab="WAIC penalty" , col=rangi2 , lwd=2 )


m5.3t <- quap(
  alist(
    D ~ dstudent( 2 , mu , sigma ) ,
    mu <- a + bM*M + bA*A ,
    a ~ dnorm( 0 , 0.2 ) ,
    bM ~ dnorm( 0 , 0.5 ) ,
    bA ~ dnorm( 0 , 0.5 ) ,
    sigma ~ dexp( 1 )
  ) , data = d )

compare(  m5.3 , m5.3t , func=PSIS )

```



######################### Exercises ######################### 
### Easy

7E1. State the three motivating criteria that define information entropy. Try to express each in your
own words.

1) The measure of uncertainty should be continuous.
- To have a useful measure of uncertainty, you need to be able to adjust for new information. If it wasn't continuous it wouldn't be flexible enough

2) The measure of uncertainty should increase as the number of possible events increases.
- More possible outcomes makes it harder to just guess the right outcome --> more uncertainty/entropy

3) The measure of uncertainty should be additive.
- The uncertainty of one event, and the uncertainty of another event, should be added together if you want to know the uncertainty of both events.


7E2. Suppose a coin is weighted such that, when it is tossed and lands on a table, it comes up heads
70% of the time. What is the entropy of this coin?
```{r}
p <- c( 0.3 , 0.7 )
-sum( p*log(p) )
```


7E3. Suppose a four-sided die is loaded such that, when tossed onto a table, it shows “1” 20%, “2”
25%, “3” 25%, and “4” 30% of the time. What is the entropy of this die?
```{r}
p <- c( 0.2 , 0.25, 0.25, 0.3 )
-sum( p*log(p) )
```

7E4. Suppose another four-sided die is loaded such that it never shows “4”. The other three sides
show equally often. What is the entropy of this die?
```{r}
p <- c( 1/3 , 1/3, 1/3 )
-sum( p*log(p) )
```
### Medium


### Hard

7H1. In 2007, The Wall Street Journal published an editorial (“We’re Number One, Alas”) with a graph of 
corporate tax rates in 29 countries plotted against tax revenue. A badly fit curve was drawn in 
(reconstructed at right), seemingly by hand, to make the argument that the relationship between tax 
rate and tax revenue increases and then declines, such that higher tax rates can actually produce 
less tax revenue. I want you to actually fit a curve to these data, found in data(Laffer).
Consider models that use tax rate to predict tax revenue. Compare, using WAIC or PSIS, a
straight-line model to any curved models you like. What do you conclude about the relationship 
between tax rate and tax revenue?

```{r}
data(Laffer)
d <- Laffer
d$Rate = standardize(d$tax_rate)
d$Revenue = standardize(d$tax_revenue)

plot(d$tax_rate, d$tax_revenue)

straight_M <- quap(
  alist(
    Revenue ~ dnorm( mu , sigma ) ,
    mu <- a + bR*Rate ,
    a ~ dnorm( 0 , 0.2 ) ,
    bR ~ dnorm( 0 , 0.5 ) ,
    sigma ~ dexp( 1 )
  ) , data = d )

curve_M <- quap(
  alist(
    Revenue ~ dnorm( mu , sigma ) ,
    mu <- a + bR*Rate + bR2*Rate^2,
    a ~ dnorm( 0 , 0.2 ) ,
    bR ~ dnorm( 0 , 0.5 ) ,
    bR2 ~ dnorm( 0 , 0.5 ) ,
    sigma ~ dexp( 1 )
  ) , data = d )

compare(straight_M, curve_M, func=WAIC)

mu <- link( curve_M , data=list(Rate=seq( from=min(d$Rate)-0.15 , to=max(d$Rate)+0.15 , length.out=30 )))
mu_mean <- apply(mu,2,mean)
mu_PI <- apply(mu,2,PI)
plot( Revenue ~ Rate , data=d )
lines( seq( from=min(d$Rate)-0.15 , to=max(d$Rate)+0.15 , length.out=30 ) , mu_mean , lwd=2 )
shade( mu_PI , seq( from=min(d$Rate)-0.15 , to=max(d$Rate)+0.15 , length.out=30 ) )

mu <- link( straight_M , data=list(Rate=seq( from=min(d$Rate)-0.15 , to=max(d$Rate)+0.15 , length.out=30 )))
mu_mean <- apply(mu,2,mean)
mu_PI <- apply(mu,2,PI)
plot( Revenue ~ Rate , data=d )
lines( seq( from=min(d$Rate)-0.15 , to=max(d$Rate)+0.15 , length.out=30 ) , mu_mean , lwd=2 )
shade( mu_PI , seq( from=min(d$Rate)-0.15 , to=max(d$Rate)+0.15 , length.out=30 ) )

```
It is not clear whether the relationship between tax rate and tax revenue is linear or curved.



7H2. In the Laffer data, there is one country with a high tax revenue that is an outlier. Use PSIS
and WAIC to measure the importance of this outlier in the models you fit in the previous problem.
Then use robust regression with a Student’s t distribution to revisit the curve fitting problem. How
much does a curved relationship depend upon the outlier point?

```{r}
set.seed(1)
PSIS_straight_M <- PSIS(straight_M,pointwise=TRUE)
set.seed(1)
WAIC_straight_M <- WAIC(straight_M,pointwise=TRUE)
plot( PSIS_straight_M$k , WAIC_straight_M$penalty , xlab="PSIS Pareto k" ,
ylab="WAIC penalty" , col=rangi2 , lwd=2 )
max(PSIS_straight_M$k)
max(WAIC_straight_M$penalty)

set.seed(1)
PSIS_curve_M <- PSIS(curve_M,pointwise=TRUE)
set.seed(1)
WAIC_curve_M <- WAIC(curve_M,pointwise=TRUE)
plot( PSIS_curve_M$k , WAIC_curve_M$penalty , xlab="PSIS Pareto k" ,
ylab="WAIC penalty" , col=rangi2 , lwd=2 )
max(PSIS_curve_M$k)
max(WAIC_curve_M$penalty)


curve_Mt <- quap(
  alist(
    Revenue ~ dstudent( 2, mu , sigma ) ,
    mu <- a + bR*Rate + bR2*Rate^2,
    a ~ dnorm( 0 , 0.2 ) ,
    bR ~ dnorm( 0 , 0.5 ) ,
    bR2 ~ dnorm( 0 , 0.5 ) ,
    sigma ~ dexp( 1 )
  ) , data = d )

compare(straight_M, curve_M, curve_Mt, func=WAIC)


mu <- link( curve_Mt , data=list(Rate=seq( from=min(d$Rate)-0.15 , to=max(d$Rate)+0.15 , length.out=30 )))
mu_mean <- apply(mu,2,mean)
mu_PI <- apply(mu,2,PI)
plot( Revenue ~ Rate , data=d )
lines( seq( from=min(d$Rate)-0.15 , to=max(d$Rate)+0.15 , length.out=30 ) , mu_mean , lwd=2 )
shade( mu_PI , seq( from=min(d$Rate)-0.15 , to=max(d$Rate)+0.15 , length.out=30 ) )
```


7H4. Recall the marriage, age, and happiness collider bias example from Chapter 6. Run models
m6.9 and m6.10 again (page 178). Compare these two models using WAIC (or PSIS, they will produce
identical results). Which model is expected to make better predictions? Which model provides the
correct causal inference about the influence of age on happiness? Can you explain why the answers
to these two questions disagree?




```{r}
d <- sim_happiness( seed=1977 , N_years=1000 )
d2 <- d[ d$age>17 , ] # only adults
d2$A <- ( d2$age - 18 ) / ( 65 - 18 )
d2$mid <- d2$married + 1

m6.9 <- quap(
  alist(
    happiness ~ dnorm( mu , sigma ),
    mu <- a[mid] + bA*A,
    a[mid] ~ dnorm( 0 , 1 ),
    bA ~ dnorm( 0 , 2 ),
    sigma ~ dexp(1)
  ) , data=d2 )


m6.10 <- quap( 
  alist(
    happiness ~ dnorm( mu , sigma ),
    mu <- a + bA*A,
    a ~ dnorm( 0 , 1 ),
    bA ~ dnorm( 0 , 2 ),
    sigma ~ dexp(1)
  ) , data=d2 )


compare(m6.9, m6.10, func=WAIC)

```
Model m6.9 is expected to make better predictions
However, model m6.10 provides the correct causal inference about the influence of age on happiness
The two answer disagree because of a collider in the model. Therefore, as long as marriage status is unknown, model m6.9 will make better predictions.


7H5. Revisit the urban fox data, data(foxes), from the previous chapter’s practice problems. Use
WAIC or PSIS based model comparison on five different models, each using weight as the outcome,
and containing these sets of predictor variables:
  (1) avgfood + groupsize + area
  (2) avgfood + groupsize
  (3) groupsize + area
  (4) avgfood
  (5) area
  
Can you explain the relative differences in WAIC scores, using the fox DAG from the previous chapter?
Be sure to pay attention to the standard error of the score differences (dSE).

```{r}
library(rethinking)
library(dagitty)
fox_dag <- dagitty( "dag {
area -> avgfood
weight <- avgfood -> groupsize
groupsize -> weight
}")

coordinates( fox_dag ) <- list( x=c(area=0.2,avgfood=0,weight=0.2,groupsize=0.4) ,
y=c(area=0,avgfood=0.2,groupsize=0.2,weight=0.4) )
drawdag( fox_dag )

#load data & standardize
data(foxes)
d <- foxes
d$A <- standardize( d$area )
d$F <- standardize( d$avgfood )
d$S <- standardize( d$groupsize )
d$W <- standardize( d$weight )

M1 <- quap(
  alist(
    W ~ dnorm( mu , sigma ) ,
    mu <- a + bF*F + bS*S + bA*A,
    a ~ dnorm( 0 , 0.2 ) ,
    bF ~ dnorm( 0 , 0.5 ) ,
    bS ~ dnorm( 0 , 0.5 ) ,
    bA ~ dnorm( 0 , 0.5 ) ,
    sigma ~ dexp( 1 )
  ) , data=d )

M2 <- quap(
  alist(
    ## F -> W <- S 
    W ~ dnorm( mu , sigma ) ,
    mu <- a + bF*F + bS*S ,
    a ~ dnorm( 0 , 0.2 ) ,
    bF ~ dnorm( 0 , 0.5 ) ,
    bS ~ dnorm( 0 , 0.5 ) ,
    sigma ~ dexp( 1 )
  ) , data=d )

M3 <- quap(
  alist(
    ## F -> W <- S 
    W ~ dnorm( mu , sigma ) ,
    mu <- a + bS*S + bA*A,
    a ~ dnorm( 0 , 0.2 ) ,
    bS ~ dnorm( 0 , 0.5 ) ,
    bA ~ dnorm( 0 , 0.5 ) ,
    sigma ~ dexp( 1 )
  ) , data=d )

M4 <- quap(
  alist(
    ## F -> W <- S 
    W ~ dnorm( mu , sigma ) ,
    mu <- a + bF*F,
    a ~ dnorm( 0 , 0.2 ) ,
    bF ~ dnorm( 0 , 0.5 ) ,
    sigma ~ dexp( 1 )
  ) , data=d )

M5 <- quap(
  alist(
    ## F -> W <- S 
    W ~ dnorm( mu , sigma ) ,
    mu <- a + bA*A,
    a ~ dnorm( 0 , 0.2 ) ,
    bA ~ dnorm( 0 , 0.5 ) ,
    sigma ~ dexp( 1 )
  ) , data=d )

compare(M1,M2,M3,M4,M5, func = WAIC)
plot(compare(M1,M2,M3,M4,M5, func = WAIC))
plot(coeftab(M1,M2,M3,M4,M5))

```
All the models containing group size as a predictor perform similar to each other. However, model 4 and 5 seem to be worse at predicting fox weight and have also higher dSE scores indicating more uncertainty. Since groupsize is a negative countereffect to the positive predictors of area and food, models that don't take group size into account make much worse predictions.