+
+
+
+
+Code 14.47
+ +
+
+
+
+
+In [57]:
+
+
+
+
+
+
+d = DataFrame(CSV.File("data/Primates301.csv", missingstring="NA"))
+describe(d)
+
+
+
+
+
+
+
+
+
+
+
+Out[57]:
+
+
+
+
+
+
+
+
+
+
+Code 14.48
+ +
+
+
+
+
+In [62]:
+
+
+
+
+
+
+dstan = d[completecases(d, ["group_size", "body", "brain"]), :]
+spp_obs = dstan.name;
+
+
+
+
+
+Code 14.49
+ +
+
+
+
+
+In [59]:
+
+
+
+
+
+
+dat_list = (
+ N_spp = nrow(dstan),
+ M = standardize(ZScoreTransform, log.(dstan.body)),
+ B = standardize(ZScoreTransform, log.(dstan.brain)),
+ G = standardize(ZScoreTransform, log.(dstan.group_size)),
+)
+
+@model function m14_9(N_spp, M, B, G)
+ σ_sq ~ Exponential()
+ bM ~ Normal(0, 0.5)
+ bG ~ Normal(0, 0.5)
+ a ~ Normal()
+ μ = @. a + bM*M + bG * G
+ B ~ MvNormal(μ, σ_sq)
+end
+
+Random.seed!(1)
+m14_9_ch = sample(m14_9(dat_list...), NUTS(), 1000)
+m14_9_df = DataFrame(m14_9_ch)
+precis(m14_9_df)
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+Code 14.50
+ +
+
+
+
+
+In [61]:
+
+
+
+
+
+
+V = DataFrame(CSV.File("data/Primates301_vcov_matrix.csv"))
+Dmat = DataFrame(CSV.File("data/Primates301_distance_matrix.csv"))
+
+# Drop index columns
+select!(V, Not(:Column1))
+select!(Dmat, Not(:Column1));
+
+p = scatter(Matrix(Dmat), Matrix(V), c=:black, ms=2)
+display("image/png", p)
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+Code 14.51
+ +
+
+
+
+
+In [63]:
+
+
+
+
+
+
+V_inds = [
+ findfirst(x -> x == n, names(V))
+ for n in spp_obs
+];
+
+V_dat = Matrix(V[V_inds, V_inds])
+R = V_dat ./ maximum(V_dat);
+
+@model function m14_10(N_spp, M, B, G, R)
+ σ_sq ~ Exponential()
+ bM ~ Normal(0, 0.5)
+ bG ~ Normal(0, 0.5)
+ a ~ Normal()
+ μ = @. a + bM*M + bG * G
+ Σ = R * σ_sq
+ B ~ MvNormal(μ, Σ)
+end
+
+Random.seed!(1)
+m14_10_ch = sample(m14_10(dat_list..., R), NUTS(), 1000)
+m14_10_df = DataFrame(m14_10_ch)
+precis(m14_10_df)
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+Code 14.52
+ +
+
+
+
+
+In [64]:
+
+
+
+
+
+
+D_inds = [
+ findfirst(x -> x == n, names(Dmat))
+ for n in spp_obs
+];
+
+D_dat = Matrix(Dmat[D_inds, D_inds])
+D_dat ./= maximum(D_dat);
+
+@model function m14_11(N_spp, M, B, G, D_dat)
+ bM ~ Normal(0, 0.5)
+ bG ~ Normal(0, 0.5)
+ a ~ Normal()
+ μ = @. a + bM*M + bG * G
+
+ η² ~ TruncatedNormal(1, 0.25, 1, Inf)
+ ρ² ~ TruncatedNormal(3, 0.25, 3, Inf)
+
+ Σ = η² * exp.(-ρ² * D_dat) + LinearAlgebra.I * (0.01 + η²)
+ B ~ MvNormal(μ, Σ)
+end
+
+Random.seed!(1)
+m14_11_ch = sample(m14_11(dat_list..., D_dat), NUTS(500, 0.65, init_ϵ=0.4), 4000)
+m14_11_df = DataFrame(m14_11_ch)
+precis(m14_11_df)
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+Code 14.53
+ +
+
+
+
+In [65]:
+
+
+
+
+
+
+plot(xlim=(0, maximum(D_dat)), ylim=(0, 1.5), xlab="phylogenetic distance", ylab="covariance")
+
+# posterior
+for r in first(eachrow(m14_11_df), 30)
+ plot!(x -> η² * exp(-ρ²*x), c=:black, alpha=0.5)
+end
+
+Random.seed!(1)
+η = rand(Normal(1, 0.25), 1000)
+ρ = rand(Normal(3, 0.25), 1000)
+d_seq = range(0, 1, length=50)
+
+K = [
+ [
+ η² * exp(-ρ²*d)
+ for d ∈ d_seq
+ ]
+ for (η², ρ²) ∈ zip(η, ρ)
+]
+K = hcat(K...)'
+
+K_μ = mean.(eachcol(K))
+K_pi = PI.(eachcol(K))
+K_pi = vcat(K_pi'...)
+
+plot!(d_seq, [K_μ K_μ], fillrange=K_pi, fillalpha=0.2, c=:black)
+annotate!([
+ (0.5, 0.5, text("prior", 12)),
+ (0.2, 0.2, text("posterior", 12))
+])
+
+
+
+
+
+
+
+
+
+
+Out[65]:
+
+
+
+
+
+