## Neal Grantham

[email protected]

# Statistical Rethinking homework solutions with Turing.jl

using Distributions  # v0.25
using LinearAlgebra
using ParetoSmooth
using DataFrames     # v1.3
using StatsPlots     # v0.14
using BSplines
using Turing         # v0.19
using CSV            # v0.9

standardize(x) = (x .- mean(x)) ./ std(x);


## Homework 1

1. Suppose the globe tossing data (Chapter 2) had turned out to be 4 water and 11 land. Construct the posterior distribution, using grid approximation. Use the same flat prior as in the book.

w = 4;
l = 11;

p_grid = range(0, 1, length=1000);
prior_grid = pdf.(Uniform(0, 1), p_grid);
like_grid = pdf.(Binomial.(w + l, p_grid), w);
post_grid = like_grid .* prior_grid;
post_grid /= sum(post_grid);

plot(p_grid, post_grid, ylab="Density", xlab="p", legend=false)


2. Now suppose the data are 4 water and 2 land. Compute the posterior again, but this time use a prior that is zero below p = 0.5 and a constant above p = 0.5. This corresponds to prior information that a majority of the Earth’s surface is water.

w = 4;
l = 2;

prior_grid = pdf.(Uniform(0.5, 1.), p_grid);
like_grid = pdf.(Binomial.(w + l, p_grid), w);
post_grid = like_grid .* prior_grid;
post_grid /= sum(post_grid);

plot(p_grid, post_grid, ylab="Density", xlab="p", legend=false)


3. For the posterior distribution from 2, compute 89% percentile and HPDI intervals. Compare the widths of these intervals. Which is wider? Why? If you had only the information in the interval, what might you misunderstand about the shape of the posterior distribution?

@model function globetoss(n, w)
p ~ Uniform(0.5, 1.)
w ~ Binomial(n, p)
end;

post = sample(globetoss(w + l, w), MH(), 10_000);

quantile(post, q=[0.055, 0.945])

## Quantiles
##   parameters      5.5%     94.5%
##       Symbol   Float64   Float64
##
##            p    0.5243    0.8790

hpd(post, alpha=0.11)

## HPD
##   parameters     lower     upper
##       Symbol   Float64   Float64
##
##            p    0.5002    0.8427


The 89% percentile interval is wider than the 89% HPD interval because the 89% HPD interval is, by definition, the narrowest region with 89% of the posterior probability. If we only had information in the interval, we might completely miss the fact that the distribution is truncated on 0.5 to 1.0!

OPTIONAL CHALLENGE. Suppose there is bias in sampling so that Land is more likely than Water to be recorded. Specifically, assume that 1-in-5 (20%) of Water samples are accidentally recorded instead as Land. First, write a generative simulation of this sampling process. Assuming the true proportion of Water is 0.70, what proportion does your simulation tend to produce instead? Second, using a simulated sample of 20 tosses, compute the unbiased posterior distribution of the true proportion of Water.

If the true proportion of Water is 0.7, and the proportion of Water samples that are accurately recorded as Water is 0.8, then we expect the simulation to tend to produce a proportion of Water of 0.7 * 0.8 = 0.56.

BiasedBinomial(n, p) = Binomial(n, 0.8p);

@model function globetoss(n, w)
p ~ Uniform(0, 1)
w ~ BiasedBinomial(n, p)
end;

n = 20;
w = 11;

post = sample(globetoss(n, w), MH(), 10_000)

## Chains MCMC chain (10000×2×1 Array{Float64, 3}):
##
## Iterations        = 1:1:10000
## Number of chains  = 1
## Samples per chain = 10000
## Wall duration     = 0.32 seconds
## Compute duration  = 0.32 seconds
## parameters        = p
## internals         = lp
##
## Summary Statistics
##   parameters      mean       std   naive_se      mcse         ess      rhat    ⋯
##       Symbol   Float64   Float64    Float64   Float64     Float64   Float64    ⋯
##
##            p    0.6804    0.1284     0.0013    0.0022   3381.0880    1.0003    ⋯
##                                                                 1 column omitted
##
## Quantiles
##   parameters      2.5%     25.0%     50.0%     75.0%     97.5%
##       Symbol   Float64   Float64   Float64   Float64   Float64
##
##            p    0.4259    0.5913    0.6838    0.7695    0.9247

density(post)


Notice that the posterior distribution is unbiased. Despite recording 11 Water out of 20 tosses (55% Water), the posterior distribution is concentrated around 0.7, the true proportion of Water.

## Homework 2

1. Construct a linear regression of weight as predicted by height, using the adults (age 18 or greater) from the Howell1 dataset. The heights listed below were recorded in the !Kung census, but weights were not recorded for these individuals. Provide predicted weights and 89% compatibility intervals for each of these individuals. That is, fill in the table below, using model-based predictions.

iheightexpected weight89% interval
1140
2160
3175
howell = CSV.read(
DataFrame,
delim=';'
)

## 544×4 DataFrame
##  Row │ height   weight    age      male
##      │ Float64  Float64   Float64  Int64
## ─────┼───────────────────────────────────
##    1 │ 151.765  47.8256      63.0      1
##    2 │ 139.7    36.4858      63.0      0
##    3 │ 136.525  31.8648      65.0      0
##    4 │ 156.845  53.0419      41.0      1
##    5 │ 145.415  41.2769      51.0      0
##    6 │ 163.83   62.9926      35.0      1
##    7 │ 149.225  38.2435      32.0      0
##    8 │ 168.91   55.48        27.0      1
##   ⋮  │    ⋮        ⋮         ⋮       ⋮
##  538 │ 142.875  34.2462      31.0      0
##  539 │  76.835   8.02291      1.0      1
##  540 │ 145.415  31.1278      17.0      1
##  541 │ 162.56   52.1631      31.0      1
##  542 │ 156.21   54.0625      21.0      0
##  543 │  71.12    8.05126      0.0      1
##  544 │ 158.75   52.5316      68.0      1
##                          529 rows omitted

adult = howell[howell.age .>= 18, :];


@model function adult_stature(w, h, h̄)
α ~ Normal(50, 10)
β ~ LogNormal(0, 1)
σ ~ Exponential(1)
μ = α .+ β .* (h .- h̄)
w ~ MvNormal(μ, σ)
end;

post = sample(model, NUTS(), 10_000);

mean(preds)  # expected weights

## Mean
##   parameters      mean
##       Symbol   Float64
##
##         w[1]   35.8370
##         w[2]   48.3595
##         w[3]   57.7743

quantile(preds, q=[0.055, 0.945])  # 89% intervals

## Quantiles
##   parameters      5.5%     94.5%
##       Symbol   Float64   Float64
##
##         w[1]   29.2081   42.6792
##         w[2]   41.5726   55.0810
##         w[3]   50.9816   64.6083


2. From the Howell1 dataset, consider only the people younger than 13 years old. Estimate the causal association between age and weight. Assume that age influences weight through two paths. First, age influences height, and height influences weight. Second, age directly influences weight through age-related changes in muscle growth and body proportions. All of this implies this causal model (DAG):

Use a linear regression to estimate the total (not just direct) causal effect of each year of growth on weight. Be sure to carefully consider the priors. Try using prior predictive simulation to assess what they imply.

@model function child_stature(w, a)
α ~ Normal(4, 1)
β ~ LogNormal(0, 0.5)
σ ~ Exponential(0.5)
μ = α .+ β .* a
w ~ MvNormal(μ, σ)
end;

child = howell[howell.age .< 13, :];


We assess our priors with prior predictive simulations.

model = child_stature(child.weight, child.age);
prior = sample(model, Prior(), 1_000);
age_seq = collect(0:12);
prior_preds = predict(child_stature(missing, age_seq), prior);
prior_preds = dropdims(prior_preds.value, dims=3)';
plot(repeat(age_seq, size(prior_preds, 2)), vec(prior_preds), ylab="Child Weight (kg)", xlab="Child Age", title="Prior Predictive Simulations", alpha=0.1, seriestype=:scatter, legend=false)


The simulations appear to represent plausible relationships between child age and weight.

What does the data look like?

plot(child.age, child.weight, ylab="Child Weight (kg)", xlab="Child Age", title="Howell's Dobe !Kung Census", seriestype=:scatter, legend=false)


Now we fit a linear regression model to the data and plot the marginal posterior distribution of β, which represents the effect of each year of growth on child weight.

post = sample(model, NUTS(), 10_000)

## Chains MCMC chain (10000×15×1 Array{Float64, 3}):
##
## Iterations        = 1001:1:11000
## Number of chains  = 1
## Samples per chain = 10000
## Wall duration     = 2.73 seconds
## Compute duration  = 2.73 seconds
## parameters        = α, β, σ
## internals         = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size
##
## Summary Statistics
##   parameters      mean       std   naive_se      mcse         ess      rhat    ⋯
##       Symbol   Float64   Float64    Float64   Float64     Float64   Float64    ⋯
##
##            α    7.0735    0.3389     0.0034    0.0047   4914.6959    1.0001    ⋯
##            β    1.3865    0.0528     0.0005    0.0007   5077.9655    1.0000    ⋯
##            σ    2.5251    0.1436     0.0014    0.0019   6465.1057    0.9999    ⋯
##                                                                 1 column omitted
##
## Quantiles
##   parameters      2.5%     25.0%     50.0%     75.0%     97.5%
##       Symbol   Float64   Float64   Float64   Float64   Float64
##
##            α    6.4042    6.8446    7.0713    7.3010    7.7484
##            β    1.2837    1.3518    1.3863    1.4217    1.4927
##            σ    2.2570    2.4268    2.5194    2.6178    2.8275

density(post["β"], ylab="Density", xlab="β", legend=false)


3. Now suppose the causal association between age and weight might be different for boys and girls. Use a single linear regression, with a categorical variable for sex, to estimate the total causal effect of age on weight separately for boys and girls. How do girls and boys differ? Provide one or more posterior contrasts as a summary.

@model function child_stature(w, a, s)
α ~ filldist(Normal(4, 1), 2)
β ~ filldist(LogNormal(0, 0.5), 2)
σ ~ Exponential(0.5)
μ = α[s] .+ β[s] .* a
w ~ MvNormal(μ, σ)
end;

model = child_stature(child.weight, child.age, child.male .+ 1);
post = sample(model, NUTS(), 10_000)

## Chains MCMC chain (10000×17×1 Array{Float64, 3}):
##
## Iterations        = 1001:1:11000
## Number of chains  = 1
## Samples per chain = 10000
## Wall duration     = 5.19 seconds
## Compute duration  = 5.19 seconds
## parameters        = α[1], α[2], β[1], β[2], σ
## internals         = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size
##
## Summary Statistics
##   parameters      mean       std   naive_se      mcse         ess      rhat    ⋯
##       Symbol   Float64   Float64    Float64   Float64     Float64   Float64    ⋯
##
##         α[1]    6.5110    0.4451     0.0045    0.0047   6904.7900    1.0000    ⋯
##         α[2]    7.1157    0.4634     0.0046    0.0053   6351.7104    1.0001    ⋯
##         β[1]    1.3552    0.0705     0.0007    0.0007   7058.9105    0.9999    ⋯
##         β[2]    1.4798    0.0721     0.0007    0.0008   6497.3970    1.0004    ⋯
##            σ    2.4642    0.1460     0.0015    0.0016   7756.6571    1.0000    ⋯
##                                                                 1 column omitted
##
## Quantiles
##   parameters      2.5%     25.0%     50.0%     75.0%     97.5%
##       Symbol   Float64   Float64   Float64   Float64   Float64
##
##         α[1]    5.6317    6.2189    6.5110    6.8026    7.3883
##         α[2]    6.1999    6.7999    7.1251    7.4348    8.0154
##         β[1]    1.2166    1.3083    1.3554    1.4020    1.4965
##         β[2]    1.3403    1.4305    1.4782    1.5287    1.6224
##            σ    2.1978    2.3646    2.4593    2.5545    2.7703

age_seq = range(0, 12, length=50);
w_f = predict(child_stature(missing, age_seq, 1), post);
w_m = predict(child_stature(missing, age_seq, 2), post);
w_diff = dropdims(w_m.value .- w_f.value, dims=3)';

w_diff_mean = vec(mean(w_diff, dims=2));
w_diff_50 = mapslices(x -> quantile(x, [0.25, 0.75]), w_diff, dims=2);
w_diff_70 = mapslices(x -> quantile(x, [0.15, 0.85]), w_diff, dims=2);
w_diff_90 = mapslices(x -> quantile(x, [0.05, 0.95]), w_diff, dims=2);

plot(age_seq, w_diff_mean, ribbon=abs.(w_diff_mean .- w_diff_90), ylab="Male - Female Child Weight Difference (kg)", xlab="Child Age");
plot!(age_seq, w_diff_mean, ribbon=abs.(w_diff_mean .- w_diff_70), legend=false);
plot!(age_seq, w_diff_mean, ribbon=abs.(w_diff_mean .- w_diff_50), legend=false);
plot!(age_seq, fill(0, length(age_seq)), line=:dash, color=:white, legend=false)


4. OPTIONAL CHALLENGE. The data in oxboys are growth records for 26 boys measured over 9 periods. I want you to model their growth. Specifically, model the increments in growth from one period (Occasion in the data table) to the next. Each increment is simply the difference between height in one occasion and height in the previous occasion. Since none of these boys shrunk during the study, all of the growth increments are greater than zero. Estimate the posterior distribution of these increments. Constrain the distribution so it is always positive — it should not be possible for the model to think that boys can shrink from year to year. Finally compute the posterior distribution of the total growth over all 9 occasions.

oxboys = CSV.read(
DataFrame,
delim=';'
);
oxboys = combine(groupby(oxboys, :Subject), :height => diff => :diff);
oxboys = transform(groupby(oxboys, :Subject), :Subject => (x -> 1:length(x)) => :time)

## 208×3 DataFrame
##  Row │ Subject  diff     time
##      │ Int64    Float64  Int64
## ─────┼─────────────────────────
##    1 │       1      2.9      1
##    2 │       1      1.4      2
##    3 │       1      2.3      3
##    4 │       1      0.6      4
##    5 │       1      2.5      5
##    6 │       1      1.5      6
##    7 │       1      1.6      7
##    8 │       1      2.5      8
##   ⋮  │    ⋮        ⋮       ⋮
##  202 │      26      0.8      2
##  203 │      26      1.6      3
##  204 │      26      1.7      4
##  205 │      26      0.5      5
##  206 │      26      2.9      6
##  207 │      26      0.8      7
##  208 │      26      0.5      8
##                193 rows omitted

@model function oxboy_growth(d, t)
μ ~ filldist(Normal(0, 0.2), length(unique(t)))
σ ~ Exponential(0.2)
d ~ MvLogNormal(μ[t], σ)
end;

prior = sample(oxboy_growth(oxboys.diff, oxboys.time), Prior(), 1_000);
prior_preds = predict(oxboy_growth(missing, 1:8), prior);
prior_preds = dropdims(prior_preds.value, dims=3)';
plot(repeat(1:8, size(prior_preds, 2)), vec(prior_preds), ylab="Height increment (cm)", xlab="Time", title="Prior Predictive Simulations", alpha=0.1, seriestype=:scatter, legend=false)


Not great prior predictive simulations. A downside to using the Log-Normal distribution here is that it produces skewed, unrealistically high height increments.

Let’s take a look at the actual data.

plot(oxboys.time, oxboys.diff, ylab="Height increment (cm)", xlab="Time", seriestype=:scatter, legend=false)


model = oxboy_growth(oxboys.diff, oxboys.time);
post = sample(model, MH(), 10_000)

## Chains MCMC chain (10000×10×1 Array{Float64, 3}):
##
## Iterations        = 1:1:10000
## Number of chains  = 1
## Samples per chain = 10000
## Wall duration     = 1.27 seconds
## Compute duration  = 1.27 seconds
## parameters        = μ[1], μ[2], μ[3], μ[4], μ[5], μ[6], μ[7], μ[8], σ
## internals         = lp
##
## Summary Statistics
##   parameters      mean       std   naive_se      mcse       ess      rhat   es ⋯
##       Symbol   Float64   Float64    Float64   Float64   Float64   Float64      ⋯
##
##         μ[1]    0.1066    0.0640     0.0006    0.0060   28.6541    1.1949      ⋯
##         μ[2]    0.0383    0.1057     0.0011    0.0103   24.9660    1.0002      ⋯
##         μ[3]    0.3841    0.1318     0.0013    0.0129   25.4927    1.1567      ⋯
##         μ[4]    0.0904    0.1118     0.0011    0.0108   26.1662    1.0012      ⋯
##         μ[5]   -0.0287    0.1563     0.0016    0.0153   21.6931    2.3444      ⋯
##         μ[6]    0.3885    0.1773     0.0018    0.0176   22.2689    1.7977      ⋯
##         μ[7]    0.1270    0.1957     0.0020    0.0192   21.8922    2.1461      ⋯
##         μ[8]    0.0392    0.1565     0.0016    0.0149   26.0338    1.2843      ⋯
##            σ    0.6287    0.0303     0.0003    0.0028   34.1437    1.0116      ⋯
##                                                                 1 column omitted
##
## Quantiles
##   parameters      2.5%     25.0%     50.0%     75.0%     97.5%
##       Symbol   Float64   Float64   Float64   Float64   Float64
##
##         μ[1]   -0.0525    0.0540    0.1406    0.1406    0.2141
##         μ[2]   -0.1046   -0.0178    0.0321    0.0321    0.3453
##         μ[3]   -0.0298    0.3492    0.3760    0.4851    0.4851
##         μ[4]   -0.1932    0.1032    0.1032    0.1926    0.1926
##         μ[5]   -0.1766   -0.1766   -0.0400    0.0496    0.2723
##         μ[6]   -0.0027    0.2869    0.2869    0.5695    0.5695
##         μ[7]   -0.0537   -0.0537    0.0118    0.3640    0.3982
##         μ[8]   -0.1523   -0.1523    0.1228    0.1228    0.3487
##            σ    0.5276    0.6318    0.6318    0.6446    0.6446

d_sim = predict(oxboy_growth(missing, 1:8), post);
d_sim = dropdims(d_sim.value, dims=3)';
d_mean = vec(mean(d_sim, dims=2));
d_90 = mapslices(x -> quantile(x, [0.05, 0.95]), d_sim, dims=2);

plot(oxboys.time, oxboys.diff, ylab="Height increment (cm)", xlab="Time", seriestype=:scatter, alpha=0.3, legend=false);
plot!(1:8, d_mean, seriestype=:scatter, yerror=(d_mean .- d_90), markersize=8, legend=false)


growth_sim = vec(sum(d_sim, dims=1));
density(growth_sim, ylab="Density", xlab="Total growth (cm)", legend=false)


## Homework 3

1. The first two problems are based on the same data. The data in foxes are 116 foxes from 30 different urban groups in England. These fox groups are like street gangs. Group size (groupsize) varies from 2 to 8 individuals. Each group maintains its own (almost exclusive) urban territory. Some territories are larger than others. The area variable encodes this information. Some territories also have more avgfood than others. And food influences the weight of each fox. Assume this DAG:

where F is avgfood, G is groupsize, A is area, and W is weight.

Use the backdoor criterion and estimate the total causal influence of A on F . What effect would increasing the area of a territory have on the amount of food inside it?

foxes = CSV.read(
DataFrame,
delim=';'
)

## 116×5 DataFrame
##  Row │ group  avgfood  groupsize  area     weight
##      │ Int64  Float64  Int64      Float64  Float64
## ─────┼─────────────────────────────────────────────
##    1 │     1     0.37          2     1.09     5.02
##    2 │     1     0.37          2     1.09     2.84
##    3 │     2     0.53          2     2.05     5.33
##    4 │     2     0.53          2     2.05     6.07
##    5 │     3     0.49          2     2.12     5.85
##    6 │     3     0.49          2     2.12     3.25
##    7 │     4     0.45          2     1.29     4.53
##    8 │     4     0.45          2     1.29     4.09
##   ⋮  │   ⋮       ⋮         ⋮         ⋮        ⋮
##  110 │    29     0.67          4     2.75     5.0
##  111 │    29     0.67          4     2.75     3.81
##  112 │    29     0.67          4     2.75     4.81
##  113 │    29     0.67          4     2.75     3.94
##  114 │    30     0.41          3     1.91     3.16
##  115 │    30     0.41          3     1.91     2.78
##  116 │    30     0.41          3     1.91     3.86
##                                    101 rows omitted

@model function fox_area_on_food_total(f, a)
α ~ Normal(0, 0.2)
β ~ Normal(0, 0.5)
σ ~ Exponential(1)
μ = α .+ β .* a
f ~ MvNormal(μ, σ)
end;

standardize(x) = (x .- mean(x)) ./ std(x);

model = fox_area_on_food_total(
standardize(foxes.avgfood),
standardize(foxes.area)
);
post = sample(model, NUTS(), 10_000)

## Chains MCMC chain (10000×15×1 Array{Float64, 3}):
##
## Iterations        = 1001:1:11000
## Number of chains  = 1
## Samples per chain = 10000
## Wall duration     = 3.27 seconds
## Compute duration  = 3.27 seconds
## parameters        = α, β, σ
## internals         = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size
##
## Summary Statistics
##   parameters      mean       std   naive_se      mcse          ess      rhat   ⋯
##       Symbol   Float64   Float64    Float64   Float64      Float64   Float64   ⋯
##
##            α   -0.0006    0.0430     0.0004    0.0004   14291.7230    0.9999   ⋯
##            β    0.8762    0.0437     0.0004    0.0004   14647.8061    0.9999   ⋯
##            σ    0.4757    0.0323     0.0003    0.0003   14247.4144    1.0000   ⋯
##                                                                 1 column omitted
##
## Quantiles
##   parameters      2.5%     25.0%     50.0%     75.0%     97.5%
##       Symbol   Float64   Float64   Float64   Float64   Float64
##
##            α   -0.0855   -0.0291   -0.0009    0.0278    0.0844
##            β    0.7889    0.8472    0.8763    0.9053    0.9627
##            σ    0.4172    0.4531    0.4739    0.4964    0.5447

area_seq = range(-2, 2, length=50);
f_preds = predict(fox_area_on_food_total(missing, area_seq), post);
f_preds = dropdims(f_preds.value, dims=3)';

f_preds_mean = vec(mean(f_preds, dims=2));
f_preds_50 = mapslices(x -> quantile(x, [0.25, 0.75]), f_preds, dims=2);
f_preds_70 = mapslices(x -> quantile(x, [0.15, 0.85]), f_preds, dims=2);
f_preds_90 = mapslices(x -> quantile(x, [0.05, 0.95]), f_preds, dims=2);

plot(area_seq, f_preds_mean, ribbon=abs.(f_preds_mean .- f_preds_90), ylab="Average Food (Standardized)", xlab="Area (Standardized)", legend=false);
plot!(area_seq, f_preds_mean, ribbon=abs.(f_preds_mean .- f_preds_70), legend=false);
plot!(area_seq, f_preds_mean, ribbon=abs.(f_preds_mean .- f_preds_50), legend=false)


There is a strong positive relationship between area and food — the larger an area, the more food available within it.

2. Now infer both the total and direct causal effects of adding food F to a territory on the weight W of foxes. Which covariates do you need to adjust for in each case? In light of your estimates from this problem and the previous one, what do you think is going on with these foxes? Feel free to speculate — all that matters is that you justify your speculation.

To answer this we will fit two models, one that captures the total causal effect of food on weight, and another that captures the direct causal effect of food on weight.

For the total model, we simply regress weight on food.

@model function fox_food_on_weight_total(w, f)
α ~ Normal(0, 0.2)
β ~ Normal(0, 0.5)
σ ~ Exponential(1)
μ = α .+ β .* f
w ~ MvNormal(μ, σ)
end;

model_total = fox_food_on_weight_total(
standardize(foxes.weight),
standardize(foxes.avgfood)
);
post_total = sample(model_total, NUTS(), 10_000)

## Chains MCMC chain (10000×15×1 Array{Float64, 3}):
##
## Iterations        = 1001:1:11000
## Number of chains  = 1
## Samples per chain = 10000
## Wall duration     = 2.31 seconds
## Compute duration  = 2.31 seconds
## parameters        = α, β, σ
## internals         = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size
##
## Summary Statistics
##   parameters      mean       std   naive_se      mcse          ess      rhat   ⋯
##       Symbol   Float64   Float64    Float64   Float64      Float64   Float64   ⋯
##
##            α   -0.0014    0.0845     0.0008    0.0007   15177.4868    0.9999   ⋯
##            β   -0.0240    0.0924     0.0009    0.0008   15388.6857    0.9999   ⋯
##            σ    1.0088    0.0659     0.0007    0.0006   14648.4465    0.9999   ⋯
##                                                                 1 column omitted
##
## Quantiles
##   parameters      2.5%     25.0%     50.0%     75.0%     97.5%
##       Symbol   Float64   Float64   Float64   Float64   Float64
##
##            α   -0.1666   -0.0576   -0.0005    0.0551    0.1629
##            β   -0.2055   -0.0861   -0.0244    0.0366    0.1573
##            σ    0.8883    0.9630    1.0060    1.0516    1.1458

density(post_total[:β], ylab="Density", xlab="Parameter Value", label="β")


The total effect of food on weight appears negligible. In the total model, β (coefficient of food) is estimated to be approximately centered at zero with an equal spread of both positive and negative values.

Now for the direct model, we regress weight on food and group size, a mediator of food and weight.

@model function fox_food_on_weight_direct(w, f, g)
α ~ Normal(0, 0.2)
β ~ Normal(0, 0.5)
γ ~ Normal(0, 0.5)
σ ~ Exponential(1)
μ = α .+ β .* f .+ γ .* g
w ~ MvNormal(μ, σ)
end;

model_direct = fox_food_on_weight_direct(
standardize(foxes.weight),
standardize(foxes.avgfood),
standardize(foxes.groupsize)
);
post_direct = sample(model_direct, NUTS(), 10_000)

## Chains MCMC chain (10000×16×1 Array{Float64, 3}):
##
## Iterations        = 1001:1:11000
## Number of chains  = 1
## Samples per chain = 10000
## Wall duration     = 4.84 seconds
## Compute duration  = 4.84 seconds
## parameters        = α, β, γ, σ
## internals         = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size
##
## Summary Statistics
##   parameters      mean       std   naive_se      mcse         ess      rhat    ⋯
##       Symbol   Float64   Float64    Float64   Float64     Float64   Float64    ⋯
##
##            α    0.0016    0.0815     0.0008    0.0008   9342.7736    1.0000    ⋯
##            β    0.4727    0.1827     0.0018    0.0022   6287.6286    1.0000    ⋯
##            γ   -0.5682    0.1827     0.0018    0.0022   6440.0546    1.0002    ⋯
##            σ    0.9621    0.0626     0.0006    0.0006   9208.1096    1.0003    ⋯
##                                                                 1 column omitted
##
## Quantiles
##   parameters      2.5%     25.0%     50.0%     75.0%     97.5%
##       Symbol   Float64   Float64   Float64   Float64   Float64
##
##            α   -0.1606   -0.0518    0.0012    0.0546    0.1624
##            β    0.1087    0.3528    0.4728    0.5950    0.8263
##            γ   -0.9222   -0.6931   -0.5714   -0.4460   -0.2011
##            σ    0.8475    0.9178    0.9594    1.0027    1.0924

density([post_direct[:β], post_direct[:γ]], ylab="Density", xlab="Parameter Value", label=["β" "γ"])


Unlike the total effect, the direct effect of food on weight is reliably positive in a model that also contains group size, which is reliably negative of an almost equal magnitude. In the direct model, β (coefficient of food) is estimated to be centered at 0.5 with most of its density between 0 and 1, while γ (coefficient of group size) is estimated to be centered at -0.6 with most of its density between -1 and 0.

Taken together, the model results suggest that the effect of food on a fox’s weight is almost entirely offset by group size. Put another way, more food available in an urban territory leads to an increase in the number of foxes within that territory, preventing any potential increase in a fox’s weight due to available food.

3. Reconsider the Table 2 Fallacy (from Lecture 6), this time with an unobserved confound U that influences both smoking S and stroke Y . Here’s the modified DAG:

First use the backdoor criterion to determine the adjustment set that allows you to estimate the causal effect of_ X on Y , i.e. P(Y ∣ do(X)) . Second explain the proper interpretation of each coefficient implied by the regression model that corresponds to the adjustment set. Which coefficients (slopes) are causal and which are not? There is not need to fit any models. Just think through the implications.

Alright, no code required for this problem. Let’s think through everything logically.

The backdoor criterion says we need to block any paths that connect X and Y that enter X from the backdoor. There are two such paths.

The first path is X ← A → Y . Here, A is a confounder of X and Y . To block the path, we include A in our adjustment set.

The second path is X ← S → Y . S too is a confounder of X and Y so we also include S in our adjustment set to the block path.

Now, with the adjustment set {A, S} we can properly estimate the causal effect of X on Y . Great! So what’s the problem?

The unobserved variable U wreaks havoc on the other coefficients in the model, contributing to a Table 2 Fallacy.

Take the variable A . As we have conditioned on A and S , and S is a collider of A and U , our estimate of the effect of A on Y in our model is biased by U .

Similarly, U is a confounder of S and Y , biasing the estimate of the effect of S on Y as well.

The takeaway here is that even with an adjustment set that allows us to properly estimate the causal effect of X on Y , we are not guaranteed unbiased estimates of other coefficients in the model which may be polluted by unobserved variables we are unable to control for.

4. OPTIONAL CHALLENGE. Write a synthetic data simulation for the causal model shown in Problem 3. Be sure to include the unobserved confound in the simulation. Choose any functional relationships that you like — you don’t have to get the epidemiology correct. You just need to honor the causal structure. Then design a regression model to estimate the influence of X on Y and use it on your synthetic data. How lrge of a sample do you need to reliably estimate P(Y ∣ do(X)) ? Define “reliably” as you like, but justify your definition.

First, we define a function to generate synthetic data. As the DAG in Problem 3 models the risk of stroke, represented by Y , we call the function generate_strokes() (at the risk of sounding like a diabolical supervillain).

function generate_strokes(n; β = 1)
a = randn(n)
u = randn(n)
s = rand.(Normal.(a .+ u))
x = rand.(Normal.(a .+ s))
y = rand.(Normal.(a .+ s .+ u .+ β .* x))
DataFrame(:x => x, :y => y, :a => a, :s => s)
end;


We’ll also define a regression model called model_strokes().

@model function model_strokes(y, x, a, s)
α ~ Normal(0, 1)
β ~ Normal(0, 1)
γ ~ Normal(0, 1)
δ ~ Normal(0, 1)
σ ~ Exponential(1)
μ = α .+ β .* x .+ γ .* a .+ δ .* s
y ~ MvNormal(μ, σ)
end;


And finally a simulation function called simulate_strokes() which generates a dataset given a sample size n, fits a model, samples from the posterior, and makes a density plot of β which it adds to an existing base plot p.

function simulate_strokes!(p, n)
strokes = generate_strokes(n)
model = model_strokes(strokes.y, strokes.x, strokes.a, strokes.s)
post = sample(model, NUTS(), 10_000)
density!(p, post[:β])
end;


Now let’s investigate the accuracy of our estimate of β under three sample sizes: 10, 100, and 1,000.

using Random
Random.seed!(42);

p = plot([1, 1], [0, 12], ylab="Density", xlab="β", xlims=[-1, 3], line=:dash, color=:white, legend=false);
p10 = deepcopy(p);
p100 = deepcopy(p);
p1000 = deepcopy(p);
for _ in 1:5
simulate_strokes!(p10, 10);
simulate_strokes!(p100, 100);
simulate_strokes!(p1000, 1000);
end;

plot!(p10, title="n = 10")


plot!(p100, title="n = 100")


plot!(p1000, title="n = 1000")


Foregoing a formal definition of “reliable” (as the question insists), it appears as though the estimates of β become reliable-ish with a sample size approaching 100.

## Homework 4

1. Revisit the marriage, age, and happiness collider bias example from Chapter 6. Run models m6.9 and m6.10 (pages 178-179). Compare these two models using PSIS and WAIC. Which model is expected to make better predictions, according to these criteria? On the basis of the causal model, how should you interpret the parameter estimates from the model preferred by PSIS and WAIC?

function generate_persons(; n_years=1000, max_age=65, n_births=20, age_of_marriage=18)
# https://github.com/rmcelreath/rethinking/blob/master/R/sim_happiness.R

invlogit(x) = 1 ./ (exp(-x) .+ 1)
age = fill(1, n_births)
married = fill(false, n_births)
happiness = collect(range(-2, 2, length=n_births))

for _ in 1:n_years
age .+= 1
append!(age, fill(1, n_births))
append!(married, fill(false, n_births))
append!(happiness, collect(range(-2, 2, length=n_births)))

for i in eachindex(age, married, happiness)
if age[i] >= age_of_marriage && !married[i]
married[i] = rand(Bernoulli(invlogit(happiness[i] - 4)))
end
end

deaths = findall(age .> max_age)
deleteat!(age, deaths)
deleteat!(married, deaths)
deleteat!(happiness, deaths)
end

DataFrame(
age = age,
married = convert.(Int, married),
happiness = happiness
)
end;

persons = generate_persons();
adults = persons[persons.age .> 17, :];

@model function with_collider(h, a, m)
α ~ filldist(Normal(0, 1), 2)
β ~ Normal(0, 2)
σ ~ Exponential(1)
μ = α[m] .+ β .* a
h .~ Normal.(μ, σ)
end;

model_with_collider = with_collider(
);
post_with_collider = sample(model_with_collider, NUTS(), 10_000)

## Chains MCMC chain (10000×16×1 Array{Float64, 3}):
##
## Iterations        = 1001:1:11000
## Number of chains  = 1
## Samples per chain = 10000
## Wall duration     = 16.23 seconds
## Compute duration  = 16.23 seconds
## parameters        = α[1], α[2], β, σ
## internals         = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size
##
## Summary Statistics
##   parameters      mean       std   naive_se      mcse         ess      rhat    ⋯
##       Symbol   Float64   Float64    Float64   Float64     Float64   Float64    ⋯
##
##         α[1]   -0.2066    0.0641     0.0006    0.0009   5215.8579    1.0000    ⋯
##         α[2]    1.2387    0.0879     0.0009    0.0013   4748.6592    1.0003    ⋯
##            β   -0.7526    0.1150     0.0011    0.0017   4351.4088    0.9999    ⋯
##            σ    1.0109    0.0233     0.0002    0.0003   7107.2851    1.0000    ⋯
##                                                                 1 column omitted
##
## Quantiles
##   parameters      2.5%     25.0%     50.0%     75.0%     97.5%
##       Symbol   Float64   Float64   Float64   Float64   Float64
##
##         α[1]   -0.3321   -0.2498   -0.2072   -0.1634   -0.0791
##         α[2]    1.0675    1.1796    1.2387    1.2964    1.4118
##            β   -0.9782   -0.8299   -0.7523   -0.6766   -0.5229
##            σ    0.9662    0.9952    1.0103    1.0261    1.0586

@model function without_collider(h, a)
α ~ Normal(0, 1)
β ~ Normal(0, 2)
σ ~ Exponential(1)
μ = α .+ β .* a
h .~ Normal.(μ, σ)
end;

model_without_collider = without_collider(
);
post_without_collider = sample(model_without_collider, NUTS(), 10_000)

## Chains MCMC chain (10000×15×1 Array{Float64, 3}):
##
## Iterations        = 1001:1:11000
## Number of chains  = 1
## Samples per chain = 10000
## Wall duration     = 7.74 seconds
## Compute duration  = 7.74 seconds
## parameters        = α, β, σ
## internals         = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size
##
## Summary Statistics
##   parameters      mean       std   naive_se      mcse         ess      rhat    ⋯
##       Symbol   Float64   Float64    Float64   Float64     Float64   Float64    ⋯
##
##            α   -0.0003    0.0772     0.0008    0.0011   4709.7550    0.9999    ⋯
##            β    0.0004    0.1337     0.0013    0.0020   4485.6562    0.9999    ⋯
##            σ    1.2160    0.0280     0.0003    0.0004   6816.9221    0.9999    ⋯
##                                                                 1 column omitted
##
## Quantiles
##   parameters      2.5%     25.0%     50.0%     75.0%     97.5%
##       Symbol   Float64   Float64   Float64   Float64   Float64
##
##            α   -0.1505   -0.0518    0.0004    0.0510    0.1525
##            β   -0.2633   -0.0914    0.0000    0.0887    0.2645
##            σ    1.1630    1.1968    1.2151    1.2347    1.2731

psis_loo(model_with_collider, post_with_collider, source=:mcmc)
# -1360

psis_loo(model_without_collider, post_without_collider, source=:mcmc)
# -1551


2. Reconsider the urban fox analysis from Homework 3. On the basis of PSIS and WAIC scores, which combination of variables best predicts body weight ( W , weight)? How would you interpret the estimates from the best scoring model?

foxes.area = standardize(foxes.area);
foxes.avgfood = standardize(foxes.avgfood);
foxes.groupsize = standardize(foxes.groupsize);
foxes.weight = standardize(foxes.weight);

@model function fox1(w, x)
α ~ Normal(0, 0.2)
β ~ Normal(0, 0.5)
σ ~ Exponential(1)
μ = α .+ β .* x
w .~ Normal.(μ, σ)
end;

@model function fox2(w, x, y)
α ~ Normal(0, 0.2)
β ~ Normal(0, 0.5)
γ ~ Normal(0, 0.5)
σ ~ Exponential(1)
μ = α .+ β .* x .+ γ .* y
w .~ Normal.(μ, σ)
end;

@model function fox3(w, x, y, z)
α ~ Normal(0, 0.2)
β ~ Normal(0, 0.5)
γ ~ Normal(0, 0.5)
δ ~ Normal(0, 0.5)
σ ~ Exponential(1)
μ = α .+ β .* x .+ γ .* y .+ δ .* z
w .~ Normal.(μ, σ)
end;

model_f = fox1(foxes.weight, foxes.avgfood);
model_g = fox1(foxes.weight, foxes.groupsize);
model_fa = fox2(foxes.weight, foxes.avgfood, foxes.area);
model_fg = fox2(foxes.weight, foxes.avgfood, foxes.groupsize);
model_fga = fox3(foxes.weight, foxes.avgfood, foxes.groupsize, foxes.area);

post_f = sample(model_f, NUTS(), 10_000);
post_g = sample(model_g, NUTS(), 10_000);
post_fa = sample(model_fa, NUTS(), 10_000);
post_fg = sample(model_fg, NUTS(), 10_000);
post_fga = sample(model_fga, NUTS(), 10_000);

psis_loo(model_f, post_f, source=:mcmc)      # -166.60
psis_loo(model_g, post_g, source=:mcmc)      # -165.24
psis_loo(model_fa, post_fa, source=:mcmc)    # -167.17
psis_loo(model_fg, post_fg, source=:mcmc)    # -161.78
psis_loo(model_fga, post_fga, source=:mcmc)  # -161.43


3. Build a predictive model of the relationship show on the cover of the book, the relationship between the timing of cherry blossoms and March temperature in the same year. The data are in cherry_blossoms. Consider at least two functions to predict doy with temp. Compare them with PSIS or WAIC.

Suppose March temperatures reach 9 degrees by the year 2050. What does your best model predict for the predictive distribution of the day-in-year the cherry trees will blossom?

blossoms = CSV.read(
DataFrame,
delim=';',
missingstring="NA"
)

## 1215×5 DataFrame
##   Row │ year   doy      temp      temp_upper  temp_lower
##       │ Int64  Int64?   Float64?  Float64?    Float64?
## ──────┼──────────────────────────────────────────────────
##     1 │   801  missing   missing     missing     missing
##     2 │   802  missing   missing     missing     missing
##     3 │   803  missing   missing     missing     missing
##     4 │   804  missing   missing     missing     missing
##     5 │   805  missing   missing     missing     missing
##     6 │   806  missing   missing     missing     missing
##     7 │   807  missing   missing     missing     missing
##     8 │   808  missing   missing     missing     missing
##   ⋮   │   ⋮       ⋮        ⋮          ⋮           ⋮
##  1209 │  2009       95   missing     missing     missing
##  1210 │  2010       95   missing     missing     missing
##  1211 │  2011       99   missing     missing     missing
##  1212 │  2012      101   missing     missing     missing
##  1213 │  2013       93   missing     missing     missing
##  1214 │  2014       94   missing     missing     missing
##  1215 │  2015       93   missing     missing     missing
##                                         1200 rows omitted

dropmissing!(blossoms, [:doy, :temp])

## 787×5 DataFrame
##  Row │ year   doy    temp     temp_upper  temp_lower
##      │ Int64  Int64  Float64  Float64?    Float64?
## ─────┼───────────────────────────────────────────────
##    1 │   851    108     7.38       12.1         2.66
##    2 │   864    100     6.42        8.69        4.14
##    3 │   866    106     6.44        8.11        4.77
##    4 │   889    104     6.83        8.48        5.19
##    5 │   891    109     6.98        8.96        5.0
##    6 │   892    108     7.11        9.11        5.11
##    7 │   894    106     6.98        8.4         5.55
##    8 │   895    104     7.08        8.57        5.58
##   ⋮  │   ⋮      ⋮       ⋮         ⋮           ⋮
##  781 │  1974     99     8.18        8.78        7.58
##  782 │  1975    100     8.18        8.76        7.6
##  783 │  1976     99     8.2         8.77        7.63
##  784 │  1977     93     8.22        8.78        7.66
##  785 │  1978    104     8.2         8.78        7.61
##  786 │  1979     97     8.28        8.83        7.73
##  787 │  1980    102     8.3         8.86        7.74
##                                      772 rows omitted

plot(blossoms.temp, blossoms.doy, ylab="Day in year of first blossom", xlab="March temperature (°C)", seriestype=:scatter, legend=false)


mean_doy = mean(blossoms.doy);
std_doy = std(blossoms.doy);
blossoms.doy = standardize(blossoms.doy);

mean_temp = mean(blossoms.temp);
std_temp = std(blossoms.temp);
blossoms.temp = standardize(blossoms.temp);

@model function blossom_line(d, t)
α ~ Normal(0, 0.2)
β ~ Normal(0, 0.5)
σ ~ Exponential(1)
μ = α .+ β .* t
d .~ Normal.(μ, σ)
end;

model_line = blossom_line(blossoms.doy, blossoms.temp);
post_line = sample(model_line, NUTS(), 10_000)

## Chains MCMC chain (10000×15×1 Array{Float64, 3}):
##
## Iterations        = 1001:1:11000
## Number of chains  = 1
## Samples per chain = 10000
## Wall duration     = 4.06 seconds
## Compute duration  = 4.06 seconds
## parameters        = α, β, σ
## internals         = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size
##
## Summary Statistics
##   parameters      mean       std   naive_se      mcse          ess      rhat   ⋯
##       Symbol   Float64   Float64    Float64   Float64      Float64   Float64   ⋯
##
##            α    0.0002    0.0328     0.0003    0.0002   17025.7063    1.0000   ⋯
##            β   -0.3252    0.0341     0.0003    0.0003   12714.6852    1.0001   ⋯
##            σ    0.9464    0.0237     0.0002    0.0002   15815.9987    0.9999   ⋯
##                                                                 1 column omitted
##
## Quantiles
##   parameters      2.5%     25.0%     50.0%     75.0%     97.5%
##       Symbol   Float64   Float64   Float64   Float64   Float64
##
##            α   -0.0650   -0.0213    0.0002    0.0220    0.0646
##            β   -0.3920   -0.3484   -0.3254   -0.3023   -0.2571
##            σ    0.9012    0.9302    0.9460    0.9620    0.9946

@model function blossom_splines(d, B)
α ~ Normal(0, 0.2)
β ~ filldist(Normal(0, 0.5), size(B, 2))
σ ~ Exponential(1)
μ = α .+ B * β
@. d ~ Normal(μ, σ)
end;

n_knots = 10;
breakpoints = range(extrema(blossoms.temp)..., length=n_knots);
basis = BSplineBasis(4, breakpoints);
B = basismatrix(basis, blossoms.temp);

model_splines = blossom_splines(blossoms.doy, B);
post_splines = sample(model_splines, NUTS(), 10_000)

## Chains MCMC chain (10000×26×1 Array{Float64, 3}):
##
## Iterations        = 1001:1:11000
## Number of chains  = 1
## Samples per chain = 10000
## Wall duration     = 29.45 seconds
## Compute duration  = 29.45 seconds
## parameters        = α, β[1], β[2], β[3], β[4], β[5], β[6], β[7], β[8], β[9], β[10], β[11], β[12], σ
## internals         = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size
##
## Summary Statistics
##   parameters      mean       std   naive_se      mcse          ess      rhat   ⋯
##       Symbol   Float64   Float64    Float64   Float64      Float64   Float64   ⋯
##
##            α   -0.0920    0.1246     0.0012    0.0019    5108.6872    1.0000   ⋯
##         β[1]    0.2760    0.3875     0.0039    0.0033   12795.9630    0.9999   ⋯
##         β[2]    0.7123    0.3345     0.0033    0.0032   10393.2794    0.9999   ⋯
##         β[3]    0.3454    0.2635     0.0026    0.0030    8340.9922    1.0001   ⋯
##         β[4]    0.5412    0.2152     0.0022    0.0025    7627.4806    1.0000   ⋯
##         β[5]    0.1595    0.1895     0.0019    0.0022    7354.2051    1.0002   ⋯
##         β[6]   -0.0165    0.1950     0.0020    0.0023    7583.9913    0.9999   ⋯
##         β[7]   -0.2038    0.2202     0.0022    0.0025    8218.6683    0.9999   ⋯
##         β[8]   -0.3293    0.2631     0.0026    0.0026    9071.6733    1.0000   ⋯
##         β[9]   -0.6002    0.3069     0.0031    0.0033    9888.1501    1.0001   ⋯
##        β[10]   -0.4692    0.3845     0.0038    0.0039   10985.7487    1.0000   ⋯
##        β[11]   -0.4978    0.4036     0.0040    0.0039   13242.1135    1.0000   ⋯
##        β[12]   -0.4522    0.3939     0.0039    0.0034   13846.2617    1.0001   ⋯
##            σ    0.9505    0.0241     0.0002    0.0002   15595.1343    1.0001   ⋯
##                                                                 1 column omitted
##
## Quantiles
##   parameters      2.5%     25.0%     50.0%     75.0%     97.5%
##       Symbol   Float64   Float64   Float64   Float64   Float64
##
##            α   -0.3364   -0.1757   -0.0924   -0.0095    0.1543
##         β[1]   -0.4796    0.0200    0.2770    0.5403    1.0330
##         β[2]    0.0638    0.4851    0.7111    0.9404    1.3618
##         β[3]   -0.1765    0.1700    0.3477    0.5150    0.8650
##         β[4]    0.1226    0.3963    0.5395    0.6861    0.9674
##         β[5]   -0.2101    0.0320    0.1594    0.2902    0.5289
##         β[6]   -0.3900   -0.1477   -0.0192    0.1163    0.3688
##         β[7]   -0.6345   -0.3535   -0.2031   -0.0538    0.2251
##         β[8]   -0.8527   -0.5060   -0.3288   -0.1523    0.1810
##         β[9]   -1.1981   -0.8084   -0.5997   -0.3943    0.0043
##        β[10]   -1.2143   -0.7340   -0.4683   -0.2046    0.2831
##        β[11]   -1.2896   -0.7674   -0.5000   -0.2315    0.3007
##        β[12]   -1.2090   -0.7199   -0.4528   -0.1860    0.3177
##            σ    0.9045    0.9341    0.9502    0.9664    0.9987

psis_loo(model_line, post_line, source=:mcmc);        # -1074.60
psis_loo(model_splines, post_splines, source=:mcmc);  # -1080.86

temp_seq = collect(range(-3, 5, length=100));
μ_post = post_line["α"] .+ post_line["β"] * temp_seq';

plot(std_temp .* blossoms.temp .+ mean_temp, std_doy .* blossoms.doy .+ mean_doy, ylab="Day in year of first blossom", xlab="March temperature (°C)", seriestype=:scatter, legend=false, alpha=0.4);
plot!(std_temp .* temp_seq .+ mean_temp, std_doy .* vec(mean(μ_post, dims=1)) .+ mean_doy, legend=false)


μ_post_at_temp_9C = rand.(Normal.(post_line["α"] .+ post_line["β"] .* (9 .- mean_temp) ./ std_temp, post_line["σ"]));

density(std_doy .* μ_post_at_temp_9C .+ mean_doy, ylab="Density", xlab="Predicted day in year of first bloom if March temp. is 9 °C", legend=false)


## Homework 5

1. The data in grants are outcomes for scientific funding applications for the Netherlands Organization for Scientific Research (NWO) from 2010–2012. (see van der Lee and Ellemers <doi:10.1073/pnas.1510159112>). These data have a similar structure to the UCB admissions data discussed in Chapter 11. Draw a DAG for this sample and then use one or more binomial GLMs to estimate the total causal effect of gender on grant awards.

Let G represent the gender of the grant applicant, D represent the discipline of the grant application, and A represent whether or not the grant was awarded.

invlogit(x) = exp(x) / (1 + exp(x));

DataFrame,
delim=';'
)

## 18×4 DataFrame
##  Row │ discipline           gender   applications  awards
##      │ String31             String1  Int64         Int64
## ─────┼────────────────────────────────────────────────────
##    1 │ Chemical sciences    m                  83      22
##    2 │ Chemical sciences    f                  39      10
##    3 │ Physical sciences    m                 135      26
##    4 │ Physical sciences    f                  39       9
##    5 │ Physics              m                  67      18
##    6 │ Physics              f                   9       2
##    7 │ Humanities           m                 230      33
##    8 │ Humanities           f                 166      32
##   ⋮  │          ⋮              ⋮          ⋮          ⋮
##   12 │ Interdisciplinary    f                  78      17
##   13 │ Earth/life sciences  m                 156      38
##   14 │ Earth/life sciences  f                 126      18
##   15 │ Social sciences      m                 425      65
##   16 │ Social sciences      f                 409      47
##   17 │ Medical sciences     m                 245      46
##   18 │ Medical sciences     f                 260      29
##                                             3 rows omitted

grants.gender = [gender == "f" ? 1 : 2 for gender in grants.gender];
disciplines = unique(grants.discipline);
grants.discipline = [
findfirst(d -> d == discipline, disciplines)
for discipline in grants.discipline
];

@model function grant_gender_on_award_total(a, g, n)
α ~ filldist(Normal(0, 1.5), 2)
p = invlogit.(α[g])
a .~ Binomial.(n, p)
end;

model_total = grant_gender_on_award_total(
grants.awards,
grants.gender,
grants.applications
);
post_total = sample(model_total, NUTS(), 10_000)

## Chains MCMC chain (10000×14×1 Array{Float64, 3}):
##
## Iterations        = 1001:1:11000
## Number of chains  = 1
## Samples per chain = 10000
## Wall duration     = 7.75 seconds
## Compute duration  = 7.75 seconds
## parameters        = α[1], α[2]
## internals         = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size
##
## Summary Statistics
##   parameters      mean       std   naive_se      mcse          ess      rhat   ⋯
##       Symbol   Float64   Float64    Float64   Float64      Float64   Float64   ⋯
##
##         α[1]   -1.7383    0.0823     0.0008    0.0009    9852.4335    0.9999   ⋯
##         α[2]   -1.5330    0.0644     0.0006    0.0005   10554.9990    1.0000   ⋯
##                                                                 1 column omitted
##
## Quantiles
##   parameters      2.5%     25.0%     50.0%     75.0%     97.5%
##       Symbol   Float64   Float64   Float64   Float64   Float64
##
##         α[1]   -1.9014   -1.7942   -1.7368   -1.6832   -1.5792
##         α[2]   -1.6601   -1.5769   -1.5320   -1.4894   -1.4077

gender_adv_total = invlogit.(post_total["α[2]"]) - invlogit.(post_total["α[1]"]);

density(gender_adv_total, ylab="Density", xlab="Gender contrast (probability)", legend=false);
plot!([0, 0], [0, 34], line=:dash, color=:white);
annotate!(-0.005, 32, text("women advantaged", :white, :right, 8));
annotate!(0.005, 32, text("men advantaged", :white, :left, 8))


2. Now estimate the direct causal effect of gender on grant awards. Compute the average direct causal effect of gender, weighting each discipline in proportion to the number of applications in the sample. Refer to the marginal effect example in Lecture 9 for help.

@model function grant_gender_on_award_direct(a, g, n, d)
α ~ filldist(Normal(0, 1.5), length(unique(g)), length(unique(d)))
for i in 1:length(a)
a[i] ~ Binomial(n[i], invlogit(α[g[i], d[i]]))
end
end;

model_direct = grant_gender_on_award_direct(
grants.awards,
grants.gender,
grants.applications,
grants.discipline
);
post_direct = sample(model_direct, NUTS(), 10_000)

## Chains MCMC chain (10000×30×1 Array{Float64, 3}):
##
## Iterations        = 1001:1:11000
## Number of chains  = 1
## Samples per chain = 10000
## Wall duration     = 4.78 seconds
## Compute duration  = 4.78 seconds
## parameters        = α[1,1], α[2,1], α[1,2], α[2,2], α[1,3], α[2,3], α[1,4], α[2,4], α[1,5], α[2,5], α[1,6], α[2,6], α[1,7], α[2,7], α[1,8], α[2,8], α[1,9], α[2,9]
## internals         = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size
##
## Summary Statistics
##   parameters      mean       std   naive_se      mcse          ess      rhat   ⋯
##       Symbol   Float64   Float64    Float64   Float64      Float64   Float64   ⋯
##
##       α[1,1]   -1.0354    0.3565     0.0036    0.0027   17314.5238    0.9999   ⋯
##       α[2,1]   -1.0055    0.2426     0.0024    0.0017   16975.3056    1.0001   ⋯
##       α[1,2]   -1.1640    0.3646     0.0036    0.0030   14896.8231    1.0000   ⋯
##       α[2,2]   -1.4197    0.2212     0.0022    0.0018   15476.4679    0.9999   ⋯
##       α[1,3]   -1.0686    0.6990     0.0070    0.0053   16916.5101    0.9999   ⋯
##       α[2,3]   -0.9901    0.2763     0.0028    0.0023   16650.9139    0.9999   ⋯
##       α[1,4]   -1.4195    0.1938     0.0019    0.0016   18235.8213    0.9999   ⋯
##       α[2,4]   -1.7713    0.1867     0.0019    0.0014   15472.5930    0.9999   ⋯
##       α[1,5]   -1.2981    0.3077     0.0031    0.0022   16814.2082    0.9999   ⋯
##       α[2,5]   -1.6534    0.1953     0.0020    0.0015   16313.0272    0.9999   ⋯
##       α[1,6]   -1.2575    0.2727     0.0027    0.0021   17493.6539    1.0001   ⋯
##       α[2,6]   -1.9976    0.2970     0.0030    0.0024   17122.6507    1.0001   ⋯
##       α[1,7]   -1.7645    0.2512     0.0025    0.0016   19403.7417    0.9999   ⋯
##       α[2,7]   -1.1263    0.1853     0.0019    0.0014   20398.8360    0.9999   ⋯
##       α[1,8]   -2.0284    0.1535     0.0015    0.0014   14443.1747    0.9999   ⋯
##       α[2,8]   -1.7047    0.1320     0.0013    0.0009   17872.9843    0.9999   ⋯
##       α[1,9]   -2.0521    0.1935     0.0019    0.0014   20877.7466    0.9999   ⋯
##       α[2,9]   -1.4575    0.1642     0.0016    0.0012   19516.6949    0.9999   ⋯
##                                                                 1 column omitted
##
## Quantiles
##   parameters      2.5%     25.0%     50.0%     75.0%     97.5%
##       Symbol   Float64   Float64   Float64   Float64   Float64
##
##       α[1,1]   -1.7650   -1.2663   -1.0244   -0.7895   -0.3700
##       α[2,1]   -1.4856   -1.1677   -1.0000   -0.8394   -0.5447
##       α[1,2]   -1.9138   -1.3993   -1.1510   -0.9162   -0.4844
##       α[2,2]   -1.8670   -1.5689   -1.4159   -1.2702   -0.9980
##       α[1,3]   -2.5336   -1.5215   -1.0411   -0.5953    0.2414
##       α[2,3]   -1.5518   -1.1774   -0.9813   -0.8001   -0.4703
##       α[1,4]   -1.8128   -1.5487   -1.4144   -1.2878   -1.0526
##       α[2,4]   -2.1513   -1.8963   -1.7662   -1.6446   -1.4145
##       α[1,5]   -1.9183   -1.4995   -1.2913   -1.0873   -0.7207
##       α[2,5]   -2.0521   -1.7784   -1.6465   -1.5223   -1.2849
##       α[1,6]   -1.8120   -1.4371   -1.2509   -1.0728   -0.7437
##       α[2,6]   -2.6055   -2.1927   -1.9856   -1.7926   -1.4437
##       α[1,7]   -2.2604   -1.9313   -1.7592   -1.5887   -1.2964
##       α[2,7]   -1.4985   -1.2483   -1.1249   -1.0006   -0.7719
##       α[1,8]   -2.3439   -2.1274   -2.0226   -1.9249   -1.7398
##       α[2,8]   -1.9689   -1.7928   -1.7031   -1.6146   -1.4503
##       α[1,9]   -2.4342   -2.1809   -2.0468   -1.9189   -1.6866
##       α[2,9]   -1.7929   -1.5640   -1.4528   -1.3469   -1.1426

gender_adv_direct = [
invlogit.(post_direct["α[2,$i]"]) - invlogit.(post_direct["α[1,$i]"])
for i in 1:length(disciplines)
];

density(gender_adv_direct[:, 1:5], ylab="Density", xlab="Gender contrast (probability)", labels=["Chemical sciences" "Physical sciences" "Physics" "Humanities" "Technical sciences"], legend=:left);
plot!([0, 0], [0, 15], line=:dash, color=:white, label="");
annotate!(-0.025, 13, text("women advantaged", :white, :right, 8));
annotate!(0.025, 13, text("men advantaged", :white, :left, 8))


density(gender_adv_direct[:, 6:9], ylab="Density", xlab="Gender contrast (probability)", labels=["Interdisciplinary" "Earth/life sciences" "Social sciences" "Medical sciences"], legend=:left);
plot!([0, 0], [0, 21], line=:dash, color=:white, label="");
annotate!(-0.025, 19, text("women advantaged", :white, :right, 8));
annotate!(0.025, 19, text("men advantaged", :white, :left, 8))


3. Considering the total effect (problem 1) and direct effect (problem 2) of gender, what causes contribute to the average difference between women and men in award rate in this sample? It is not necessary to say whether or not there is evidence of discrimination. Simply explain how the direct effects you have estimated make sense (or not) of the total effect.

## Homework 6

The theme of this homework is tadpoles. You must keep them alive.

1. Conduct a prior predictive simulation for the Reedfrog model. By this I mean to simulate the prior distribution of tank survival probabilities αj . Start by using this prior:

\\begin{aligned} \\alpha\_j &\\sim \\text{Normal}(\\bar{\\alpha}, \\sigma) \\\\ \\bar{\\alpha} &\\sim \\text{Normal}(0, 1) \\\\ \\sigma &\\sim \\text{Exponential}(1) \\\\ \\end{aligned}

Be sure to transform the αj values to the probability scale for plotting and summary. How does increasing the width of the prior on σ change the prior distribution αj ? You might try Exponential(10) and Exponential(0.1) for example.

# redefine a version of invlogit that's more computationally stable
invlogit(x) = exp(x - log(1 + exp(x)));

σ_dists = [Exponential(0.1), Exponential(1), Exponential(10)];

function sim_priors(σ_dist, n=10_000)
σ = rand(σ_dist, n)
ᾱ = rand(Normal(0, 1), n)
α = rand.(Normal.(ᾱ, σ))
invlogit.(α)
end;

sims = hcat(sim_priors.(σ_dists)...);

density(sims, xlims=[0.025, 0.975], ylab="Density",xlab="Simulated prior survival probability", labels=["Exp(0.1)" "Exp(1)" "Exp(10)"], legend=:top)


2. Revisit the Reedfrog survival data. Start with the varying effects model from the book and the lecture. Then modify it to estimate the causal effects of the treatement variables pred and size, including how size might modify the effect of predation. An easy approach is to estimate the effect for each combination of pred and size. Justify your model with a DAG of this experiment.

tadpoles = CSV.read(
DataFrame,
delim=';'
)

## 48×5 DataFrame
##  Row │ density  pred     size     surv   propsurv
##      │ Int64    String7  String7  Int64  Float64
## ─────┼────────────────────────────────────────────
##    1 │      10  no       big          9  0.9
##    2 │      10  no       big         10  1.0
##    3 │      10  no       big          7  0.7
##    4 │      10  no       big         10  1.0
##    5 │      10  no       small        9  0.9
##    6 │      10  no       small        9  0.9
##    7 │      10  no       small       10  1.0
##    8 │      10  no       small        9  0.9
##   ⋮  │    ⋮        ⋮        ⋮       ⋮       ⋮
##   42 │      35  pred     big         12  0.342857
##   43 │      35  pred     big         13  0.371429
##   44 │      35  pred     big         14  0.4
##   45 │      35  pred     small       22  0.628571
##   46 │      35  pred     small       12  0.342857
##   47 │      35  pred     small       31  0.885714
##   48 │      35  pred     small       17  0.485714
##                                    33 rows omitted

tadpoles.tank = collect(1:size(tadpoles, 1));
tadpoles.pred = [pred == "no" ? 1 : 2 for pred in tadpoles.pred];
tadpoles.size = [size == "small" ? 1 : 2 for size in tadpoles.size];

@model function survival_vareff(surv, density, tank, pred, size)
σ ~ Exponential(1)
ᾱ ~ Normal(0, 1.5)
α ~ filldist(Normal(ᾱ, σ), length(unique(tank)))
β ~ filldist(Normal(0, 1.5), length(unique(pred)), length(unique(size)))
for i in 1:length(surv)
surv[i] ~ Binomial(density[i], invlogit(α[i] + β[pred[i], size[i]]))
end
end;

model_vareff = survival_vareff(
);
post_vareff = sample(model_vareff, NUTS(), 10_000);

βs = ["β[1,1]", "β[1,2]", "β[2,1]", "β[2,2]"];

density([post_vareff[β] for β in βs], ylab="Density", xlab="β", labels=["Small, no predators" "Big, no predators" "Small, predators" "Big, predators"], legend=:topleft);
plot!([0, 0], [0, 0.55], line=:dash, color=:white, label="")


3. Now estimate the causal effect of density on survival. Consider whether pred modifies the effect of density. There are several good ways to include density in your Binomial GLM. You could treat it as a continuous regression variable (possibly standardized). Or you could convert it to an ordered category (with three levels). Compare the σ (tank standard deviation) posterior distribution to σ from your model in Problem 2. How are they different? Why?

@model function survival_with_density(surv, density, tank, pred, size)
σ ~ Exponential(1)
ᾱ ~ Normal(0, 1.5)
α ~ filldist(Normal(ᾱ, σ), length(unique(tank)))
β ~ filldist(Normal(0, 1.5), length(unique(pred)), length(unique(size)))
γ ~ filldist(Normal(0, 1.5), length(unique(pred)))
for i in 1:length(surv)
surv[i] ~ Binomial(density[i], invlogit(α[i] + β[pred[i], size[i]] + γ[pred[i]] * density[i]))
end
end;

model_with_density = survival_with_density(
);
post_with_density = sample(model_with_density, NUTS(), 10_000);

βs = ["β[1,1]", "β[1,2]", "β[2,1]", "β[2,2]"];

density([post_with_density[β] for β in βs], ylab="Density", xlab="β", labels=["Small, no predators" "Big, no predators" "Small, predators" "Big, predators"], legend=:topleft);
plot!([0, 0], [0, 0.6], line=:dash, color=:white, label="")


γs = ["γ[1]", "γ[2]"];

density([post_with_density[γ] for γ in γs], ylab="Density", xlab="γ", labels=["No predators" "Predators"], legend=:topleft);
plot!([0, 0], [0, 30], line=:dash, color=:white, label="")


density([post_vareff["σ"], post_with_density["σ"]], ylab="Density", xlab="σ", labels=["Binomial GLM without density" "Binomial GLM with density"])


## Homework 7

1. The data in bangladesh are 1934 women from the 1989 Bangladesh Fertility Survey. For each woman, we know which district she lived, her number of living_children, her age_centered, whether she lived in an urban center, and finally whether or not she used contraception (use_contraception). In the first problem, I only want you to investigate the proportion of women using contraception in each distirct. Use partial pooling (varying effects). Then compare the varying effects estimates to the raw empirical proportion in each district. Explain the differences between the estimates and the data. Note that district number 54 is absent in the data. This may cause some problems in indexing the parameters. Pay special attention to district number 54’s estimate.

bangladesh = CSV.read(
DataFrame,
delim=';'
);

rename!(
"use.contraception" => :use_contraception,
"living.children" => :living_children,
"age.centered" => :age_centered
)

## 1934×6 DataFrame
##   Row │ woman  district  use_contraception  living_children  age_centered  urb ⋯
##       │ Int64  Int64     Int64              Int64            Float64       Int ⋯
## ──────┼─────────────────────────────────────────────────────────────────────────
##     1 │     1         1                  0                4       18.44        ⋯
##     2 │     2         1                  0                1       -5.5599
##     3 │     3         1                  0                3        1.44
##     4 │     4         1                  0                4        8.44
##     5 │     5         1                  0                1      -13.559       ⋯
##     6 │     6         1                  0                1      -11.56
##     7 │     7         1                  0                4       18.44
##     8 │     8         1                  0                4       -3.5599
##   ⋮   │   ⋮       ⋮              ⋮                 ⋮              ⋮          ⋮ ⋱
##  1928 │  1928        61                  1                3       -9.5599      ⋯
##  1929 │  1929        61                  0                3       -2.5599
##  1930 │  1930        61                  0                4       14.44
##  1931 │  1931        61                  0                3       -4.5599
##  1932 │  1932        61                  0                4       14.44        ⋯
##  1933 │  1933        61                  0                1      -13.56
##  1934 │  1934        61                  0                4       10.44
##                                                   1 column and 1919 rows omitted

bangladesh_districts = combine(
:woman => length => :n_total,
:use_contraception => sum => :n_use_contraception
)

## 60×3 DataFrame
##  Row │ district  n_total  n_use_contraception
##      │ Int64     Int64    Int64
## ─────┼────────────────────────────────────────
##    1 │        1      117                   30
##    2 │        2       20                    7
##    3 │        3        2                    2
##    4 │        4       30                   15
##    5 │        5       39                   14
##    6 │        6       65                   19
##    7 │        7       18                    5
##    8 │        8       37                   14
##   ⋮  │    ⋮         ⋮              ⋮
##   54 │       55        6                    1
##   55 │       56       45                   26
##   56 │       57       27                    5
##   57 │       58       33                   15
##   58 │       59       10                    1
##   59 │       60       32                    7
##   60 │       61       42                    9
##                                45 rows omitted

@model function bangladesh_district_contraception_use(d, n, k)
σ ~ Exponential(1)
ᾱ ~ Normal(0, 1.5)
α ~ filldist(Normal(ᾱ, σ), length(d))
k .~ Binomial.(n, invlogit.(α))
end;

);
NUTS(),
10_000
);

post_p = invlogit.(
);

p_df = DataFrame(
estimate = vec(mean(post_p, dims=1)),
);
sort!(p_df, :empirical)

## 60×3 DataFrame
##  Row │ district  estimate  empirical
##      │ Int64     Float64   Float64
## ─────┼───────────────────────────────
##    1 │       11  0.182979  0.0
##    2 │       49  0.30574   0.0
##    3 │       24  0.244268  0.0714286
##    4 │       10  0.251211  0.0769231
##    5 │       59  0.276014  0.1
##    6 │       55  0.320645  0.166667
##    7 │       27  0.23933   0.181818
##    8 │       57  0.261808  0.185185
##   ⋮  │    ⋮         ⋮          ⋮
##   54 │       37  0.448628  0.538462
##   55 │       42  0.443219  0.545455
##   56 │       16  0.46939   0.55
##   57 │       56  0.52331   0.577778
##   58 │       14  0.597111  0.627119
##   59 │       34  0.567127  0.657143
##   60 │        3  0.444418  1.0
##                       45 rows omitted

plot(1:nrow(p_df), [p_df.empirical p_df.estimate], seriestype=:scatter, labels=["Empirical proportion of contraception use" "Estimated proportion of contraception use"], legend=:topleft)


2. First, draw a DAG that includes all five variables. Use contraception C , Age A , Children K , Urban U , and District D . You don’t have to be an expert on fertility. But do think about which variables can influence which other variables. Second, design an estimation strategy to identify both the total and direct causal effects of living in an urban center on contraceptive use. Those are your estimands. Consider causal relationships among the variables. Then use your DAG to justify an adjustment set that will yield the estimate of the causal effect of urban living on contraceptive use. Do not run a statistical model (yet). I just want you to try to design an analysis. There is no firm right answer. Just apply the backdoor criterion and rules of d-separation (the elemental confounds) correctly to the DAG you design.

3. Now build one or more statistical models to estimate the total and the direct causal effects of urban living on contraceptive use. Again include district as a simple varying effect (as in problem 1) so that each district has its own average contraceptive use. You may also want to stratify the effect of urban living by district. If you do, think carefully about how to do this statistically.

@model function bangladesh_urban_on_contraception_use_total(c, u, d, n_districts)
σ ~ Exponential(1)
ᾱ ~ Normal(0, 1.5)
α ~ filldist(Normal(ᾱ, σ), n_districts)
β ~ Normal(0, 0.5)
for i in 1:length(c)
c[i] ~ Bernoulli(invlogit(α[d[i]] + β * u[i]))
end
end;

);
NUTS(),
5_000
);



@model function bangladesh_urban_on_contraception_use_direct(c, u, d, a, k, n_districts)
σ ~ Exponential(1)
ᾱ ~ Normal(0, 1.5)
α ~ filldist(Normal(ᾱ, σ), n_districts)
β ~ Normal(0, 0.5)
γ ~ Normal(0, 0.5)
δ ~ Dirichlet(length(unique(k)), 2)
for i in 1:length(c)
c[i] ~ Bernoulli(invlogit(α[d[i]] + β * u[i] + γ * a[i] + δ[k[i]]))
end
end;

);
NUTS(),
5_000
);



## Homework 8

This assignment is difficult to complete without support for LKJCholesky in Turing.jl. There is an open issue to address this at github.com/TuringLang/Turing.jl/issues/1629.

invlogit(x) = 1 ./ (exp(-x) .+ 1);

DataFrame
)

## 153×9 DataFrame
##  Row │ dyad_id  A      B      like_AB  like_BA  dislike_AB  dislike_BA  A_name ⋯
##      │ Int64    Int64  Int64  Int64    Int64    Int64       Int64       String ⋯
## ─────┼──────────────────────────────────────────────────────────────────────────
##    1 │       1      1      2        0        3           0           0  ROMUL  ⋯
##    2 │       2      1      3        3        3           0           0  ROMUL
##    3 │       3      1      4        0        0           0           0  ROMUL
##    4 │       4      1      5        0        0           0           0  ROMUL
##    5 │       5      1      6        0        0           0           0  ROMUL  ⋯
##    6 │       6      1      7        0        0           0           0  ROMUL
##    7 │       7      1      8        0        0           0           0  ROMUL
##    8 │       8      1      9        0        0           0           1  ROMUL
##   ⋮  │    ⋮       ⋮      ⋮       ⋮        ⋮         ⋮           ⋮          ⋮   ⋱
##  147 │     147     14     18        0        0           0           0  ALBERT ⋯
##  148 │     148     15     16        0        1           0           0  AMAND
##  149 │     149     15     17        0        0           0           0  AMAND
##  150 │     150     15     18        0        0           0           0  AMAND
##  151 │     151     16     17        0        0           0           0  BASIL  ⋯
##  152 │     152     16     18        0        0           0           0  BASIL
##  153 │     153     17     18        3        3           0           0  ELIAS
##                                                   2 columns and 138 rows omitted

@model function monk_reciprocity(like_ab, like_ba, dyad_id, n_dyads)
σ ~ Exponential(1)
P ~ LKJ(2, 1)
d[i, :] = σ .* cholesky(P).L * randn(2)
end
α ~ Normal(0, 1)
p_ab = invlogit.(α .+ d[:, 1])
p_ba = invlogit.(α .+ d[:, 2])
like_ab .~ Binomial.(3, p_ab)
like_ba .~ Binomial.(3, p_ba)
end;

model_monk_reciprocity = monk_reciprocity(
monks.like_AB,
monks.like_BA,
nrow(monks)
);
@time post_monk_reciprocity = sample(model_monk_reciprocity, NUTS(), 500)

## 200.584295 seconds (2.73 G allocations: 283.401 GiB, 19.40% gc time, 0.01% compilation time)

## Chains MCMC chain (500×18×1 Array{Float64, 3}):
##
## Iterations        = 251:1:750
## Number of chains  = 1
## Samples per chain = 500
## Wall duration     = 187.29 seconds
## Compute duration  = 187.29 seconds
## parameters        = σ, P[1,1], P[2,1], P[1,2], P[2,2], α
## internals         = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size
##
## Summary Statistics
##   parameters      mean       std   naive_se      mcse       ess      rhat   es ⋯
##       Symbol   Float64   Float64    Float64   Float64   Float64   Float64      ⋯
##
##            σ    0.3219    0.0000     0.0000    0.0000    2.0080    0.9980      ⋯
##       P[1,1]    1.0000    0.0000     0.0000    0.0000       NaN       NaN      ⋯
##       P[2,1]   -0.5501    0.0000     0.0000    0.0000       NaN       NaN      ⋯
##       P[1,2]   -0.5501    0.0000     0.0000    0.0000       NaN       NaN      ⋯
##       P[2,2]    1.0000    0.0000     0.0000    0.0000    2.0080    0.9980      ⋯
##            α   -1.8666    0.0000     0.0000    0.0000    2.0080    0.9980      ⋯
##                                                                 1 column omitted
##
## Quantiles
##   parameters      2.5%     25.0%     50.0%     75.0%     97.5%
##       Symbol   Float64   Float64   Float64   Float64   Float64
##
##            σ    0.3219    0.3219    0.3219    0.3219    0.3219
##       P[1,1]    1.0000    1.0000    1.0000    1.0000    1.0000
##       P[2,1]   -0.5501   -0.5501   -0.5501   -0.5501   -0.5501
##       P[1,2]   -0.5501   -0.5501   -0.5501   -0.5501   -0.5501
##       P[2,2]    1.0000    1.0000    1.0000    1.0000    1.0000
##            α   -1.8666   -1.8666   -1.8666   -1.8666   -1.8666


## Homework 9

trips = CSV.read(
DataFrame,
delim=';'
)

## 14364×8 DataFrame
##    Row │ month  day     year   id     age    kg.meat  hours       datatype
##        │ Int64  Int64?  Int64  Int64  Int64  Float64  Float64?    Float64
## ───────┼───────────────────────────────────────────────────────────────────
##      1 │    10       2   1981   3043     67     0.0         6.97       1.0
##      2 │    10       3   1981   3043     67     0.0         9.0        1.0
##      3 │    10       4   1981   3043     67     0.0         1.55       1.0
##      4 │    10       5   1981   3043     67     4.5         8.0        1.0
##      5 │    10       6   1981   3043     67     0.0         3.0        1.0
##      6 │    10       7   1981   3043     67     0.0         7.5        1.0
##      7 │    10       8   1981   3043     67     0.0         9.0        1.0
##      8 │     2       6   1982   3043     68     0.0         7.3        1.0
##    ⋮   │   ⋮      ⋮       ⋮      ⋮      ⋮       ⋮         ⋮          ⋮
##  14358 │     8      26   2006  60121     21     4.96  missing          3.0
##  14359 │     8      27   2006  60121     21     0.0   missing          3.0
##  14360 │     2       2   2007  60121     22     0.0   missing          3.0
##  14361 │     2       3   2007  60121     22     0.0   missing          3.0
##  14362 │     2      21   2007  60121     22     0.0   missing          3.0
##  14363 │     2      22   2007  60121     22     4.87  missing          3.0
##  14364 │     3       5   2007  60121     22     3.15  missing          3.0
##                                                          14349 rows omitted

rename!(trips, "kg.meat" => :kg_meat);
trips.success = convert.(Int, trips.kg_meat .> 0.);

@model function hunting1(s, a)
α ~ Beta(4, 4)
β ~ Exponential(2)
γ ~ Exponential(2)
δ ~ Exponential(0.5)
p = α .* exp.(-γ .* a) .* (1 .- exp.(-β .* a)) .^ δ
s .~ Bernoulli.(p)
end;

model_hunting1 = hunting1(trips.success, trips.age ./ maximum(trips.age));
@time post_hunting1 = sample(model_hunting1, NUTS(), 1_000)

##  70.565193 seconds (77.91 M allocations: 63.124 GiB, 2.56% gc time, 0.03% compilation time)

## Chains MCMC chain (1000×16×1 Array{Float64, 3}):
##
## Iterations        = 501:1:1500
## Number of chains  = 1
## Samples per chain = 1000
## Wall duration     = 58.44 seconds
## Compute duration  = 58.44 seconds
## parameters        = α, β, γ, δ
## internals         = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size
##
## Summary Statistics
##   parameters      mean       std   naive_se      mcse        ess      rhat   e ⋯
##       Symbol   Float64   Float64    Float64   Float64    Float64   Float64     ⋯
##
##            α    0.8954    0.0403     0.0013    0.0019   392.0014    1.0021     ⋯
##            β    7.0077    0.4904     0.0155    0.0233   349.7137    1.0053     ⋯
##            γ    0.7696    0.0588     0.0019    0.0027   416.8462    1.0013     ⋯
##            δ    5.4265    0.8002     0.0253    0.0358   364.8487    1.0034     ⋯
##                                                                 1 column omitted
##
## Quantiles
##   parameters      2.5%     25.0%     50.0%     75.0%     97.5%
##       Symbol   Float64   Float64   Float64   Float64   Float64
##
##            α    0.8054    0.8704    0.9004    0.9244    0.9626
##            β    6.0889    6.6834    6.9955    7.3301    7.9452
##            γ    0.6362    0.7337    0.7785    0.8117    0.8707
##            δ    3.9814    4.8914    5.3698    5.9156    7.2558

age_seq = collect(range(10, 80, length=100) ./ maximum(trips.age));
prob_success = vec(post_hunting1["α"]) .* exp.(-vec(post_hunting1["γ"]) * age_seq') .* (1 .- exp.(-post_hunting1["β"] * age_seq')) .^ post_hunting1["δ"];
mean_prob_success = vec(mean(prob_success, dims=1));

plot(age_seq .* maximum(trips.age), mean_prob_success, ylab="Probability of successful hunt", xlab="Age of hunter", legend=false)


@model function hunting2(s, a, hunter_id, n_hunters)
σ ~ Exponential(1)
τ ~ Exponential(1)
η ~ Exponential(1)
P ~ LKJ(2, η)
v = Matrix(undef, n_hunters, 2)
z ~ filldist(Normal(0, 1), 2, n_hunters)
D = diagm([σ, τ])
for i in 1:n_hunters
v[i, :] = D * cholesky(P).L * vec(z[:, i])
end
α ~ Beta(4, 4)
β ~ Normal(0, 0.5)
γ ~ Normal(0, 0.5)
δ ~ Exponential(0.5)
ε = exp.(β .+ vec(v[1:n_hunters, 1]))
ζ = exp.(γ .+ vec(v[1:n_hunters, 2]))
p = α .* exp.(-ζ[hunter_id] .* a) .* ((1 .- exp.(-ε[hunter_id] .* a)) .^ δ)
s .~ Bernoulli.(p)
end;

# reduce size of dataset for demonstration purposes!
trips2 = trips[1:111, :];

trips2.hunter_id = [findfirst(id .== unique(trips2.id)) for id in trips2.id];

model_hunting2 = hunting2(trips2.success, trips2.age ./ maximum(trips2.age), trips2.hunter_id, maximum(trips2.hunter_id));
post_hunting2 = sample(model_hunting2, NUTS(), 4000);

post_hunting2[["α", "β", "γ", "δ", "σ", "τ"]]

## Chains MCMC chain (4000×6×1 Array{Float64, 3}):
##
## Iterations        = 1001:1:5000
## Number of chains  = 1
## Samples per chain = 4000
## Wall duration     = 63.85 seconds
## Compute duration  = 63.85 seconds
## parameters        = σ, τ, α, β, γ, δ
## internals         =
##
## Summary Statistics
##   parameters      mean       std   naive_se      mcse         ess      rhat    ⋯
##       Symbol   Float64   Float64    Float64   Float64     Float64   Float64    ⋯
##
##            σ    0.8937    0.8794     0.0139    0.0133   3643.9991    0.9998    ⋯
##            τ    0.4947    0.4469     0.0071    0.0112   1647.0174    1.0003    ⋯
##            α    0.5722    0.1218     0.0019    0.0024   3251.8017    0.9998    ⋯
##            β    0.0744    0.4869     0.0077    0.0082   3560.4339    1.0006    ⋯
##            γ   -0.2992    0.3665     0.0058    0.0072   2742.6807    0.9998    ⋯
##            δ    0.3093    0.3012     0.0048    0.0054   3243.0843    0.9998    ⋯
##                                                                 1 column omitted
##
## Quantiles
##   parameters      2.5%     25.0%     50.0%     75.0%     97.5%
##       Symbol   Float64   Float64   Float64   Float64   Float64
##
##            σ    0.0202    0.2604    0.6226    1.2457    3.2864
##            τ    0.0167    0.1760    0.3842    0.6886    1.5920
##            α    0.3540    0.4824    0.5651    0.6600    0.8156
##            β   -0.9074   -0.2477    0.0736    0.3976    1.0232
##            γ   -1.0965   -0.5331   -0.2765   -0.0409    0.3427
##            δ    0.0063    0.0886    0.2202    0.4376    1.0983


Updated January 27, 2022 with Homework 1 and Homework 2 solutions

Updated February 10, 2022 with Homework 3 solutions

Updated February 24, 2022 with Homework 4 solutions

Updated March 10, 2022 with Homework 5 solutions

Updated March 24, 2022 with Homework 6 solutions

Updated April 7, 2022 with Homework 7 solutions

Updated April 21, 2022 with Homework 8 solutions

Updated May 5, 2022 with Homework 9 solutions

February 24, 2022  @nsgrantham