m8.1.jl · StatisticalRethinking.jl

m8.1stan

m8.1stan is the first model in the Statistical Rethinking book (pp. 249) using Stan. Here we will use Turing's NUTS support, which is currently (2018) the original NUTS by Hoffman & Gelman, http://www.stat.columbia.edu/~gelman/research/published/nuts.pdf and not the one that's in Stan 2.18.2, i.e., Appendix A.5 in: https://arxiv.org/abs/1701.02434

The StatisticalRethinking pkg uses, e.g., Turing, CSV, DataFrames

using StatisticalRethinking

d = CSV.read(joinpath(dirname(Base.pathof(StatisticalRethinking)), "..", "data",
    "rugged.csv"), delim=';')
size(d) # Should be 234x51

(234, 51)

Apply log() to each element in rgdppc_2000 column and add it as a new column

d = hcat(d, map(log, d[Symbol("rgdppc_2000")]))
rename!(d, :x1 => :log_gdp) # Rename our col x1 => log_gdp

234 rows × 52 columns

	isocode	isonum	country	rugged	rugged_popw	rugged_slope	rugged_lsd	rugged_pc	land_area	lat	lon	soil	desert	tropical	dist_coast	near_coast	gemstones	rgdppc_2000	rgdppc_1950_m	rgdppc_1975_m	rgdppc_2000_m	rgdppc_1950_2000_m	q_rule_law	cont_africa	cont_asia	cont_europe	cont_oceania	cont_north_america	cont_south_america	legor_gbr	legor_fra	legor_soc	legor_deu	legor_sca	colony_esp	colony_gbr	colony_fra	colony_prt	colony_oeu	africa_region_n	africa_region_s	africa_region_w	africa_region_e	africa_region_c	slave_exports	dist_slavemkt_atlantic	dist_slavemkt_indian	dist_slavemkt_saharan	dist_slavemkt_redsea	pop_1400	european_descent	log_gdp
	String⍰	Int64⍰	String⍰	Float64⍰	Float64⍰	Float64⍰	Float64⍰	Float64⍰	Float64⍰	Float64⍰	Float64⍰	Float64⍰	Float64⍰	Float64⍰	Float64⍰	Float64⍰	Int64⍰	Float64⍰	Float64⍰	Float64⍰	Float64⍰	Float64⍰	Float64⍰	Int64⍰	Int64⍰	Int64⍰	Int64⍰	Int64⍰	Int64⍰	Int64⍰	Int64⍰	Int64⍰	Int64⍰	Int64⍰	Int64⍰	Int64⍰	Int64⍰	Int64⍰	Int64⍰	Int64⍰	Int64⍰	Int64⍰	Int64⍰	Int64⍰	Float64⍰	Float64⍰	Float64⍰	Float64⍰	Float64⍰	Int64⍰	Float64⍰	Float64⍰
1	ABW	533	Aruba	0.462	0.38	1.226	0.144	0.0	18.0	12.508	-69.97	21.324	0.0	100.0	0.001	100.0	0	missing	missing	missing	missing	missing	missing	0	0	0	0	1	0	0	1	0	0	0	0	0	0	0	0	0	0	0	0	0	0.0	missing	missing	missing	missing	614	missing	missing
2	AFG	4	Afghanistan	2.518	1.469	7.414	0.72	39.004	65209.0	33.833	66.026	27.849	4.356	0.0	0.922	0.0	0	missing	644.756	720.633	565.231	679.791	-1.687	0	1	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	0	0	0	0	0.0	missing	missing	missing	missing	1870829	0.0	missing
3	AGO	24	Angola	0.858	0.714	2.274	0.228	4.906	124670.0	-12.299	17.551	26.676	0.425	44.346	0.428	13.1587	47756	1794.73	1051.82	1073.04	765.215	1106.76	-1.567	1	0	0	0	0	0	0	1	0	0	0	0	0	0	1	0	0	0	0	0	1	3.61e6	5.669	6.981	4.926	3.872	1223208	2.0	7.49261
4	AIA	660	Anguilla	0.013	0.01	0.026	0.006	0.0	9.0	18.231	-63.064	100.0	0.0	100.0	0.0	100.0	0	missing	missing	missing	missing	missing	missing	0	0	0	0	1	0	missing	missing	missing	missing	missing	0	0	0	0	0	0	0	0	0	0	0.0	missing	missing	missing	missing	missing	missing	missing
5	ALB	8	Albania	3.427	1.597	10.451	1.006	62.133	2740.0	41.143	20.07	68.088	0.0	0.0	0.048	94.6919	0	3703.11	1001.34	2289.47	2741.42	1931.78	-0.82	0	0	1	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	0	0	0	0.0	missing	missing	missing	missing	200000	100.0	8.21693
6	AND	20	Andorra	5.717	6.722	17.774	1.616	99.064	47.0	42.551	1.576	0.0	0.0	0.0	0.134	0.0	0	missing	missing	missing	missing	missing	1.515	0	0	1	0	0	0	0	1	0	0	0	0	0	0	0	0	0	0	0	0	0	0.0	missing	missing	missing	missing	missing	missing	missing
7	ANT	530	Netherlands Antilles	0.255	0.039	0.68	0.08	0.0	80.0	12.725	-68.157	24.595	0.0	74.555	0.001	100.0	0	missing	missing	missing	missing	missing	missing	0	0	0	0	1	0	0	1	0	0	0	0	0	0	0	0	0	0	0	0	0	0.0	missing	missing	missing	missing	1802	missing	missing
8	ARE	784	United Arab Emirates	0.769	0.316	2.112	0.191	6.142	8360.0	23.913	54.331	0.0	77.28	0.0	0.065	75.7464	0	20604.5	15797.6	25465.0	17567.9	20120.0	0.913	0	1	0	0	0	0	1	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0.0	missing	missing	missing	missing	19200	0.0	9.93326
9	ARG	32	Argentina	0.775	0.22	2.268	0.226	9.407	273669.0	-35.396	-65.17	35.678	0.0	0.0	0.352	13.0167	0	12173.7	4986.73	8122.5	8543.56	6926.81	0.033	0	0	0	0	0	1	0	1	0	0	0	1	0	0	0	0	0	0	0	0	0	0.0	missing	missing	missing	missing	276632	89.889	9.40703
10	ARM	51	Armenia	2.688	0.934	8.178	0.799	50.556	2820.0	40.294	44.938	30.148	0.0	0.0	0.348	0.0	0	2421.99	missing	missing	4565.03	missing	-0.453	0	1	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	0	0	0	0.0	missing	missing	missing	missing	105743	0.5	7.79234
11	ASM	16	American Samoa	2.647	1.63	7.032	0.792	44.303	20.0	-14.2	-170.388	100.0	0.0	100.0	0.0	100.0	0	missing	missing	missing	missing	missing	missing	0	0	0	1	0	0	1	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0.0	missing	missing	missing	missing	missing	missing	missing
12	ATG	28	Antigua and Barbuda	0.006	0.003	0.012	0.003	0.0	44.0	17.271	-61.8	100.0	0.0	100.0	0.001	100.0	0	10022.0	missing	missing	missing	missing	0.99	0	0	0	0	1	0	1	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0.0	missing	missing	missing	missing	747	missing	9.21254
13	AUS	36	Australia	0.143	0.183	0.405	0.045	0.685	768230.0	-25.733	134.487	14.248	10.889	14.68	0.336	20.9399	264154	25417.4	7411.58	13169.8	21605.3	13184.2	1.773	0	0	0	1	0	0	1	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0.0	missing	missing	missing	missing	200000	89.954	10.1432
14	AUT	40	Austria	3.513	1.152	11.095	1.008	54.307	8245.0	47.589	14.14	55.098	0.0	0.0	0.242	2.25634	0	28987.8	3706.07	11646.4	20691.4	11601.7	1.853	0	0	1	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	0	0	0.0	missing	missing	missing	missing	1250000	98.4	10.2746
15	AZE	31	Azerbaijan	1.672	0.534	5.08	0.49	27.713	8260.5	40.288	47.528	60.957	0.0	0.0	0.584	0.0	0	2570.94	missing	missing	2538.14	missing	-1.007	0	1	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	0	0	0	0.0	missing	missing	missing	missing	226250	4.0	7.85203
16	BDI	108	Burundi	1.78	1.586	4.721	0.5	27.519	2568.0	-3.365	29.887	24.049	0.0	75.354	0.906	0.0	0	621.652	360.142	534.858	495.977	530.763	-1.25	1	0	0	0	0	0	0	1	0	0	0	0	0	0	0	1	0	0	0	1	0	86.997	10.626	2.57	3.719	2.215	239693	0.0	6.43238
17	BEL	56	Belgium	0.388	0.261	1.239	0.109	0.178	3023.0	50.642	4.661	48.967	0.0	0.0	0.113	45.9637	0	27303.0	5462.2	12440.8	20656.5	12223.3	1.413	0	0	1	0	0	0	0	1	0	0	0	0	0	0	0	0	0	0	0	0	0	0.0	missing	missing	missing	missing	773466	98.0	10.2148
18	BEN	204	Benin	0.141	0.099	0.377	0.045	0.014	11062.0	9.65	2.339	74.973	0.0	97.489	0.367	10.2806	0	959.222	1083.63	969.318	1282.9	1082.47	-0.23	1	0	0	0	0	0	0	1	0	0	0	0	0	1	0	0	0	0	1	0	0	457000.0	5.121	9.234	2.835	3.902	348867	0.0	6.86612
19	BFA	854	Burkina Faso	0.236	0.214	0.638	0.066	0.051	27360.0	12.274	-1.747	56.909	0.0	49.82	0.735	0.0	0	998.416	474.313	669.515	921.09	719.511	-0.52	1	0	0	0	0	0	0	1	0	0	0	0	0	1	0	0	0	0	1	0	0	167000.0	4.775	9.299	2.764	4.239	692976	0.0	6.90617
20	BGD	50	Bangladesh	0.186	0.065	0.502	0.053	1.817	13017.0	23.848	90.27	41.883	0.0	75.038	0.071	71.3844	0	1479.09	539.544	528.607	861.7	599.122	-0.807	0	1	0	0	0	0	1	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0.0	missing	missing	missing	missing	9103705	0.0	7.29918
21	BGR	100	Bulgaria	1.479	0.711	4.498	0.427	20.196	11063.0	42.765	25.239	82.101	0.0	0.0	0.149	31.4395	0	5979.17	1651.03	5830.81	5349.92	4567.5	-0.16	0	0	1	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	0	0	0	0.0	missing	missing	missing	missing	894427	98.85	8.69604
22	BHR	48	Bahrain	0.231	0.163	0.627	0.063	0.0	71.0	26.025	50.565	0.0	0.0	0.0	0.002	100.0	0	15928.1	2104.45	3922.26	5058.84	3776.9	0.457	0	1	0	0	0	0	1	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0.0	missing	missing	missing	missing	31897	0.0	9.67584
23	BHS	44	Bahamas	0.055	0.1	0.144	0.017	0.0	1001.0	24.255	-76.61	13.924	0.0	100.0	0.003	100.0	0	16977.2	missing	missing	missing	missing	1.063	0	0	0	0	1	0	1	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0.0	missing	missing	missing	missing	1276	missing	9.73963
24	BIH	70	Bosnia and Herzegovina	2.311	1.288	7.075	0.665	40.253	5120.0	44.175	17.784	76.888	0.0	0.0	0.117	42.7623	0	5294.98	missing	missing	5572.28	missing	-0.44	0	0	1	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	0	0	0	0.0	missing	missing	missing	missing	366341	100.0	8.57451
⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮

Now we need to drop every row where rgdppc2000 == missing When this (https://github.com/JuliaData/DataFrames.jl/pull/1546) hits DataFrame it'll be conceptually easier: i.e., completecases!(d, :rgdppc2000)

notisnan(e) = !ismissing(e)
dd = d[map(notisnan, d[:rgdppc_2000]), :]
size(dd) # should equal 170 x 52

@model m8_1stan(y, x₁, x₂) = begin
    σ ~ Truncated(Cauchy(0, 2), 0, Inf)
    βR ~ Normal(0, 10)
    βA ~ Normal(0, 10)
    βAR ~ Normal(0, 10)
    α ~ Normal(0, 100)

    for i ∈ 1:length(y)
        y[i] ~ Normal(α + βR * x₁[i] + βA * x₂[i] + βAR * x₁[i] * x₂[i], σ)
    end
end

m8_1stan (generic function with 4 methods)

Test to see that the model is sane. Use 2000 for now, as in the book. Need to set the same stepsize and adapt_delta as in Stan...

posterior = sample(m8_1stan(dd[:,:log_gdp], dd[:,:rugged], dd[:,:cont_africa]),
    Turing.NUTS(2000, 1000, 0.95))
describe(posterior)

Output reg. params of interest: Mean SD Naive SE MCSE ESS α 9.2140454953 0.416410339 0.00931121825 0.0303436655 188.324543 βA -1.9414588557 0.373885658 0.00836033746 0.0583949856 40.994586 βR -0.1987645549 0.158902372 0.00355316505 0.0128657961 152.541295 σ 0.9722532977 0.440031013 0.00983939257 0.0203736871 466.473854 βAR 0.3951414223 0.187780491 0.00419889943 0.0276680621 46.062071

Here's the map2stan output: Mean StdDev lower 0.89 upper 0.89 n_eff Rhat a 9.24 0.14 9.03 9.47 291 1 bR -0.21 0.08 -0.32 -0.07 306 1 bA -1.97 0.23 -2.31 -1.58 351 1 bAR 0.40 0.13 0.20 0.63 350 1 sigma 0.95 0.05 0.86 1.03 566 1

This page was generated using Literate.jl.