米国における州別離婚率

参考資料

Richard McElreath, Statistical Rethinking: A Bayesian Course with Examples in R and Stan, CRC Press/Taylor & Francis Group, 2016. Chapter14

データの概要

米国における州別の結婚率と離婚率のデータ。サンプル数は州の数と同じ50。今回使用する変数は次の7つ。

• Location: 州の名前
• Population: 州の人口
• MedianAgeMarriage: 結婚年齢の中央値
• Marriage: 結婚率
• Marriage.SE: 結婚率の標準誤差
• Divorce: 離婚率
• Divorce.SE: 離婚率の標準誤差

結婚率と離婚率には標準誤差の情報も与えられている。

library(rethinking)

## Loading required package: rstan

## Loading required package: ggplot2

## Loading required package: StanHeaders

## rstan (Version 2.17.3, GitRev: 2e1f913d3ca3)

## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)

## Loading required package: parallel

## rethinking (Version 1.59)

library(tidyverse)

## ─ Attaching packages ──────────────────── tidyverse 1.2.1 ─

## ✔ tibble  1.4.2     ✔ purrr   0.2.4
## ✔ tidyr   0.8.0     ✔ dplyr   0.7.4
## ✔ readr   1.1.1     ✔ stringr 1.3.0
## ✔ tibble  1.4.2     ✔ forcats 0.3.0

## ─ Conflicts ───────────────────── tidyverse_conflicts() ─
## ✖ tidyr::extract() masks rstan::extract()
## ✖ dplyr::filter()  masks stats::filter()
## ✖ dplyr::lag()     masks stats::lag()
## ✖ purrr::map()     masks rethinking::map()

library(rstan)
library(ggmcmc)
library(bayesplot)

## This is bayesplot version 1.5.0

## - Plotting theme set to bayesplot::theme_default()

## - Online documentation at mc-stan.org/bayesplot

data(WaffleDivorce)
d <- WaffleDivorce
d <- d %>%
  select(Location, Population, MedianAgeMarriage, Marriage, Marriage.SE, Divorce, Divorce.SE)
head(d)

分析の目的

結婚年齢(MedianAgeMarriage)・結婚率(Marriage)が離婚率(Divorce)とどのように関係しているかを明らかにする。このとき、結婚率と離婚率の標準誤差の情報も活用する。

モデル式

離婚率の観測値\(D_{OBS}\)は、平均\(d_{est}\)、標準偏差\(D_{SE}\)の正規分布から生成されると考える。同様に、結婚率の観測値\(M_{OBS}\)は、平均\(m_{est}\)、標準偏差\(M_{SE}\)の正規分布から生成されると考える。\(d_{est}\)が州の真の離婚率であり、その平均値\(\mu\)は結婚年齢\(A\)と結婚率\(m_{est}\)の線形結合で表される。

\[ \begin{align} d_{est}[i] &\sim {\rm Normal}(\mu[i], \sigma) &i&=1,2,...,N\\ \mu[i] &\sim \alpha + \beta_A A[i] + \beta_M m_{est}[i] &i&=1,2,...,N\\ D_{obs}[i] &\sim {\rm Normal}(d_{est}[i], D_{SE}[i]) &i&=1,2,...,N\\ M_{obs}[i] &\sim {\rm Normal}(M_{est}[i], M_{SE}[i]) &i&=1,2,...,N\\ \alpha &\sim {\rm Normal}(0, 10) \\ \beta_A &\sim {\rm Normal}(0, 10) \\ \beta_M &\sim {\rm Normal}(0, 10) \\ \sigma &\sim {\rm Cauchy}(0, 2.5) \end{align} \]

パラメータ推定実行

data <- list(N=nrow(d), D_OBS=d$Divorce, D_SE=d$Divorce.SE,
             M_OBS=d$Marriage, M_SE=d$Marriage.SE, A=d$MedianAgeMarriage)
fit <- stan(file="WaffleDivorce.stan", data=data, iter=500, chains=4)

## In file included from file15546365c6b22.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:14:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core/matrix_vari.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat/fun/Eigen_NumTraits.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Core:531:
## /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
##     #pragma clang diagnostic pop
##                              ^
## In file included from file15546365c6b22.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:14:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core/matrix_vari.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat/fun/Eigen_NumTraits.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Dense:2:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/LU:47:
## /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
##     #pragma clang diagnostic pop
##                              ^
## In file included from file15546365c6b22.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:14:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core/matrix_vari.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat/fun/Eigen_NumTraits.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Dense:3:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Cholesky:12:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Jacobi:29:
## /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
##     #pragma clang diagnostic pop
##                              ^
## In file included from file15546365c6b22.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:14:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core/matrix_vari.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat/fun/Eigen_NumTraits.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Dense:3:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Cholesky:43:
## /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
##     #pragma clang diagnostic pop
##                              ^
## In file included from file15546365c6b22.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:14:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core/matrix_vari.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat/fun/Eigen_NumTraits.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Dense:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/QR:17:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Householder:27:
## /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
##     #pragma clang diagnostic pop
##                              ^
## In file included from file15546365c6b22.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:14:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core/matrix_vari.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat/fun/Eigen_NumTraits.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Dense:5:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/SVD:48:
## /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
##     #pragma clang diagnostic pop
##                              ^
## In file included from file15546365c6b22.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:14:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core/matrix_vari.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat/fun/Eigen_NumTraits.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Dense:6:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Geometry:58:
## /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
##     #pragma clang diagnostic pop
##                              ^
## In file included from file15546365c6b22.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:14:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core/matrix_vari.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat/fun/Eigen_NumTraits.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Dense:7:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Eigenvalues:58:
## /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
##     #pragma clang diagnostic pop
##                              ^
## In file included from file15546365c6b22.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:36:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core/operator_unary_plus.hpp:7:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/scal/fun/constants.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/BH/include/boost/math/constants/constants.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/BH/include/boost/math/tools/convert_from_string.hpp:15:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/BH/include/boost/lexical_cast.hpp:32:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/BH/include/boost/lexical_cast/try_lexical_convert.hpp:42:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/BH/include/boost/lexical_cast/detail/converter_lexical.hpp:52:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/BH/include/boost/container/container_fwd.hpp:61:
## /Library/Frameworks/R.framework/Versions/3.5/Resources/library/BH/include/boost/container/detail/std_fwd.hpp:27:1: warning: inline namespaces are a C++11 feature [-Wc++11-inline-namespace]
## BOOST_MOVE_STD_NS_BEG
## ^
## /Library/Frameworks/R.framework/Versions/3.5/Resources/library/BH/include/boost/move/detail/std_ns_begin.hpp:18:34: note: expanded from macro 'BOOST_MOVE_STD_NS_BEG'
##    #define BOOST_MOVE_STD_NS_BEG _LIBCPP_BEGIN_NAMESPACE_STD
##                                  ^
## /Library/Developer/CommandLineTools/usr/include/c++/v1/__config:439:52: note: expanded from macro '_LIBCPP_BEGIN_NAMESPACE_STD'
## #define _LIBCPP_BEGIN_NAMESPACE_STD namespace std {inline namespace _LIBCPP_NAMESPACE {
##                                                    ^
## In file included from file15546365c6b22.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:12:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat.hpp:83:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/csr_extract_u.hpp:6:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Sparse:26:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/SparseCore:66:
## /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
##     #pragma clang diagnostic pop
##                              ^
## In file included from file15546365c6b22.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:12:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat.hpp:83:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/csr_extract_u.hpp:6:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Sparse:27:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/OrderingMethods:71:
## /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
##     #pragma clang diagnostic pop
##                              ^
## In file included from file15546365c6b22.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:12:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat.hpp:83:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/csr_extract_u.hpp:6:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Sparse:29:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/SparseCholesky:43:
## /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
##     #pragma clang diagnostic pop
##                              ^
## In file included from file15546365c6b22.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:12:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat.hpp:83:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/csr_extract_u.hpp:6:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Sparse:32:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/SparseQR:35:
## /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
##     #pragma clang diagnostic pop
##                              ^
## In file included from file15546365c6b22.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math.hpp:4:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:12:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat.hpp:83:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/csr_extract_u.hpp:6:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/Sparse:33:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/IterativeLinearSolvers:46:
## /Library/Frameworks/R.framework/Versions/3.5/Resources/library/RcppEigen/include/Eigen/src/Core/util/ReenableStupidWarnings.h:10:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
##     #pragma clang diagnostic pop
##                              ^
## In file included from file15546365c6b22.cpp:636:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/rstan/include/rstan/rstaninc.hpp:3:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/rstan/include/rstan/stan_fit.hpp:36:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/src/stan/services/optimize/bfgs.hpp:11:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/src/stan/optimization/bfgs.hpp:9:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/src/stan/optimization/lbfgs_update.hpp:6:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/BH/include/boost/circular_buffer.hpp:54:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/BH/include/boost/circular_buffer/details.hpp:20:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/BH/include/boost/move/move.hpp:30:
## In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/BH/include/boost/move/iterator.hpp:27:
## /Library/Frameworks/R.framework/Versions/3.5/Resources/library/BH/include/boost/move/detail/iterator_traits.hpp:29:1: warning: inline namespaces are a C++11 feature [-Wc++11-inline-namespace]
## BOOST_MOVE_STD_NS_BEG
## ^
## /Library/Frameworks/R.framework/Versions/3.5/Resources/library/BH/include/boost/move/detail/std_ns_begin.hpp:18:34: note: expanded from macro 'BOOST_MOVE_STD_NS_BEG'
##    #define BOOST_MOVE_STD_NS_BEG _LIBCPP_BEGIN_NAMESPACE_STD
##                                  ^
## /Library/Developer/CommandLineTools/usr/include/c++/v1/__config:439:52: note: expanded from macro '_LIBCPP_BEGIN_NAMESPACE_STD'
## #define _LIBCPP_BEGIN_NAMESPACE_STD namespace std {inline namespace _LIBCPP_NAMESPACE {
##                                                    ^
## 15 warnings generated.
## 
## SAMPLING FOR MODEL 'WaffleDivorce' NOW (CHAIN 1).
## Rejecting initial value:
##   Error evaluating the log probability at the initial value.
## Exception: validate transformed params: mu[1] is -21.5682, but must be greater than or equal to 0  (in 'model1554647037e5a_WaffleDivorce' at line 20)
## 
## Rejecting initial value:
##   Error evaluating the log probability at the initial value.
## Exception: validate transformed params: mu[2] is -15.0303, but must be greater than or equal to 0  (in 'model1554647037e5a_WaffleDivorce' at line 20)
## 
## Rejecting initial value:
##   Error evaluating the log probability at the initial value.
## Exception: validate transformed params: mu[2] is 110.961, but must be less than or equal to 100  (in 'model1554647037e5a_WaffleDivorce' at line 20)
## 
## 
## Gradient evaluation took 6e-05 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0.6 seconds.
## Adjust your expectations accordingly!
## 
## 
## Iteration:   1 / 500 [  0%]  (Warmup)
## Iteration:  50 / 500 [ 10%]  (Warmup)
## Iteration: 100 / 500 [ 20%]  (Warmup)
## Iteration: 150 / 500 [ 30%]  (Warmup)
## Iteration: 200 / 500 [ 40%]  (Warmup)
## Iteration: 250 / 500 [ 50%]  (Warmup)
## Iteration: 251 / 500 [ 50%]  (Sampling)
## Iteration: 300 / 500 [ 60%]  (Sampling)
## Iteration: 350 / 500 [ 70%]  (Sampling)
## Iteration: 400 / 500 [ 80%]  (Sampling)
## Iteration: 450 / 500 [ 90%]  (Sampling)
## Iteration: 500 / 500 [100%]  (Sampling)
## 
##  Elapsed Time: 0.482626 seconds (Warm-up)
##                0.528095 seconds (Sampling)
##                1.01072 seconds (Total)
## 
## 
## SAMPLING FOR MODEL 'WaffleDivorce' NOW (CHAIN 2).
## Rejecting initial value:
##   Error evaluating the log probability at the initial value.
## Exception: validate transformed params: mu[1] is -122.056, but must be greater than or equal to 0  (in 'model1554647037e5a_WaffleDivorce' at line 20)
## 
## Rejecting initial value:
##   Error evaluating the log probability at the initial value.
## Exception: validate transformed params: mu[2] is -13.4502, but must be greater than or equal to 0  (in 'model1554647037e5a_WaffleDivorce' at line 20)
## 
## Rejecting initial value:
##   Error evaluating the log probability at the initial value.
## Exception: validate transformed params: mu[2] is 106.469, but must be less than or equal to 100  (in 'model1554647037e5a_WaffleDivorce' at line 20)
## 
## Rejecting initial value:
##   Error evaluating the log probability at the initial value.
## Exception: validate transformed params: mu[1] is -96.6625, but must be greater than or equal to 0  (in 'model1554647037e5a_WaffleDivorce' at line 20)
## 
## Rejecting initial value:
##   Error evaluating the log probability at the initial value.
## Exception: validate transformed params: mu[2] is -8.0192, but must be greater than or equal to 0  (in 'model1554647037e5a_WaffleDivorce' at line 20)
## 
## Rejecting initial value:
##   Error evaluating the log probability at the initial value.
## Exception: validate transformed params: mu[1] is 140.711, but must be less than or equal to 100  (in 'model1554647037e5a_WaffleDivorce' at line 20)
## 
## Rejecting initial value:
##   Error evaluating the log probability at the initial value.
## Exception: validate transformed params: mu[1] is 145.739, but must be less than or equal to 100  (in 'model1554647037e5a_WaffleDivorce' at line 20)
## 
## Rejecting initial value:
##   Error evaluating the log probability at the initial value.
## Exception: validate transformed params: mu[1] is -83.6408, but must be greater than or equal to 0  (in 'model1554647037e5a_WaffleDivorce' at line 20)
## 
## Rejecting initial value:
##   Error evaluating the log probability at the initial value.
## Exception: validate transformed params: mu[1] is -47.1632, but must be greater than or equal to 0  (in 'model1554647037e5a_WaffleDivorce' at line 20)
## 
## Rejecting initial value:
##   Error evaluating the log probability at the initial value.
## Exception: validate transformed params: mu[2] is -1.42394, but must be greater than or equal to 0  (in 'model1554647037e5a_WaffleDivorce' at line 20)
## 
## Rejecting initial value:
##   Error evaluating the log probability at the initial value.
## Exception: validate transformed params: mu[30] is -3.68492, but must be greater than or equal to 0  (in 'model1554647037e5a_WaffleDivorce' at line 20)
## 
## Rejecting initial value:
##   Error evaluating the log probability at the initial value.
## Exception: validate transformed params: mu[1] is -87.2772, but must be greater than or equal to 0  (in 'model1554647037e5a_WaffleDivorce' at line 20)
## 
## Rejecting initial value:
##   Error evaluating the log probability at the initial value.
## Exception: validate transformed params: mu[1] is -155.365, but must be greater than or equal to 0  (in 'model1554647037e5a_WaffleDivorce' at line 20)
## 
## Rejecting initial value:
##   Error evaluating the log probability at the initial value.
## Exception: validate transformed params: mu[24] is -6.96464, but must be greater than or equal to 0  (in 'model1554647037e5a_WaffleDivorce' at line 20)
## 
## Rejecting initial value:
##   Error evaluating the log probability at the initial value.
## Exception: validate transformed params: mu[1] is -14.1656, but must be greater than or equal to 0  (in 'model1554647037e5a_WaffleDivorce' at line 20)
## 
## Rejecting initial value:
##   Error evaluating the log probability at the initial value.
## Exception: validate transformed params: mu[4] is -11.1857, but must be greater than or equal to 0  (in 'model1554647037e5a_WaffleDivorce' at line 20)
## 
## Rejecting initial value:
##   Error evaluating the log probability at the initial value.
## Exception: validate transformed params: mu[2] is -9.62549, but must be greater than or equal to 0  (in 'model1554647037e5a_WaffleDivorce' at line 20)
## 
## Rejecting initial value:
##   Error evaluating the log probability at the initial value.
## Exception: validate transformed params: mu[1] is 166.421, but must be less than or equal to 100  (in 'model1554647037e5a_WaffleDivorce' at line 20)
## 
## Rejecting initial value:
##   Error evaluating the log probability at the initial value.
## Exception: validate transformed params: mu[1] is -47.7702, but must be greater than or equal to 0  (in 'model1554647037e5a_WaffleDivorce' at line 20)
## 
## 
## Gradient evaluation took 2.8e-05 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0.28 seconds.
## Adjust your expectations accordingly!
## 
## 
## Iteration:   1 / 500 [  0%]  (Warmup)
## Iteration:  50 / 500 [ 10%]  (Warmup)
## Iteration: 100 / 500 [ 20%]  (Warmup)
## Iteration: 150 / 500 [ 30%]  (Warmup)
## Iteration: 200 / 500 [ 40%]  (Warmup)
## Iteration: 250 / 500 [ 50%]  (Warmup)
## Iteration: 251 / 500 [ 50%]  (Sampling)
## Iteration: 300 / 500 [ 60%]  (Sampling)
## Iteration: 350 / 500 [ 70%]  (Sampling)
## Iteration: 400 / 500 [ 80%]  (Sampling)
## Iteration: 450 / 500 [ 90%]  (Sampling)
## Iteration: 500 / 500 [100%]  (Sampling)
## 
##  Elapsed Time: 0.512381 seconds (Warm-up)
##                1.29647 seconds (Sampling)
##                1.80885 seconds (Total)
## 
## 
## SAMPLING FOR MODEL 'WaffleDivorce' NOW (CHAIN 3).
## 
## Gradient evaluation took 2.8e-05 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0.28 seconds.
## Adjust your expectations accordingly!
## 
## 
## Iteration:   1 / 500 [  0%]  (Warmup)
## Iteration:  50 / 500 [ 10%]  (Warmup)
## Iteration: 100 / 500 [ 20%]  (Warmup)
## Iteration: 150 / 500 [ 30%]  (Warmup)
## Iteration: 200 / 500 [ 40%]  (Warmup)
## Iteration: 250 / 500 [ 50%]  (Warmup)
## Iteration: 251 / 500 [ 50%]  (Sampling)
## Iteration: 300 / 500 [ 60%]  (Sampling)
## Iteration: 350 / 500 [ 70%]  (Sampling)
## Iteration: 400 / 500 [ 80%]  (Sampling)
## Iteration: 450 / 500 [ 90%]  (Sampling)
## Iteration: 500 / 500 [100%]  (Sampling)
## 
##  Elapsed Time: 0.527645 seconds (Warm-up)
##                0.660888 seconds (Sampling)
##                1.18853 seconds (Total)
## 
## 
## SAMPLING FOR MODEL 'WaffleDivorce' NOW (CHAIN 4).
## Rejecting initial value:
##   Error evaluating the log probability at the initial value.
## Exception: validate transformed params: mu[1] is -1.56846, but must be greater than or equal to 0  (in 'model1554647037e5a_WaffleDivorce' at line 20)
## 
## Rejecting initial value:
##   Error evaluating the log probability at the initial value.
## Exception: validate transformed params: mu[3] is 113.362, but must be less than or equal to 100  (in 'model1554647037e5a_WaffleDivorce' at line 20)
## 
## Rejecting initial value:
##   Error evaluating the log probability at the initial value.
## Exception: validate transformed params: mu[1] is 158.884, but must be less than or equal to 100  (in 'model1554647037e5a_WaffleDivorce' at line 20)
## 
## Rejecting initial value:
##   Error evaluating the log probability at the initial value.
## Exception: validate transformed params: mu[6] is 101.965, but must be less than or equal to 100  (in 'model1554647037e5a_WaffleDivorce' at line 20)
## 
## Rejecting initial value:
##   Error evaluating the log probability at the initial value.
## Exception: validate transformed params: mu[1] is -109.658, but must be greater than or equal to 0  (in 'model1554647037e5a_WaffleDivorce' at line 20)
## 
## Rejecting initial value:
##   Error evaluating the log probability at the initial value.
## Exception: validate transformed params: mu[4] is -6.99572, but must be greater than or equal to 0  (in 'model1554647037e5a_WaffleDivorce' at line 20)
## 
## 
## Gradient evaluation took 2.9e-05 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0.29 seconds.
## Adjust your expectations accordingly!
## 
## 
## Iteration:   1 / 500 [  0%]  (Warmup)
## Iteration:  50 / 500 [ 10%]  (Warmup)
## Iteration: 100 / 500 [ 20%]  (Warmup)
## Iteration: 150 / 500 [ 30%]  (Warmup)
## Iteration: 200 / 500 [ 40%]  (Warmup)
## Iteration: 250 / 500 [ 50%]  (Warmup)
## Iteration: 251 / 500 [ 50%]  (Sampling)
## Iteration: 300 / 500 [ 60%]  (Sampling)
## Iteration: 350 / 500 [ 70%]  (Sampling)
## Iteration: 400 / 500 [ 80%]  (Sampling)
## Iteration: 450 / 500 [ 90%]  (Sampling)
## Iteration: 500 / 500 [100%]  (Sampling)
## 
##  Elapsed Time: 0.730123 seconds (Warm-up)
##                1.3269 seconds (Sampling)
##                2.05702 seconds (Total)

推定結果

posterior <- as.matrix(fit)
ms <- rstan::extract(fit)
ms_df <- as.data.frame(ms)

(1)\(\beta_A\)が負であることから、結婚年齢が低い州ほど離婚率が高いことが分かる。また、\(\beta_M\)が正であることから、結婚率が高いほど離婚率も高いことが分かる。また、結婚率よりも結婚年齢の方が離婚率に対する影響が大きい。

plot_title <- ggtitle("Posterior distributions of coefficients", "with medians and 80% intervals")
mcmc_areas(posterior, pars = c("beta_A", "beta_M"), prob = 0.8) +
  plot_title

(2)推定された離婚率は、離婚率の標準誤差が大きい州ほど観測された離婚率との差が大きくなっている。

d_est <- data.frame(t(apply(ms$d_est, 2, quantile, c(0.025, 0.5, 0.975))))
fig_df <- data.frame(d_est=d_est$X50., D_OBS=data$D_OBS, D_SE=data$D_SE)

ggplot(data=fig_df, aes(x=D_SE, y=d_est-D_OBS)) +
  geom_point() +
  xlab("Divorce observed standard error") +
  ylab("Divorce estimated - divorce observed") +
  geom_hline(yintercept = 0, linetype = "dashed")