イントロダクション
この vignette
は、順序以上の尺度にある離散予測変数の特別な取り扱い方である単調効果についてである（Bürkner
& Charpentier, in
review）。単調（すなわち、応答と単調に増加または減少する関係）を持つものとしてモデルしたい予測変数は、整数値または順序付き因子のいずれかでなければならない。連続予測変数とは異なり、予測変数のカテゴリ（または整数）は、応答変数でのそれらの効果に関して等距離であることは仮定されない。その代わりに、隣接する予測変数カテゴリ（または整数）間の距離がデータから推定され、カテゴリ間で変化することがある。これは、次のようにパラメータ化することで実現される。
1つのパラメータ \(b\)
は、通常の回帰パラメータと同様、効果の方向と大きさに注意する。単調効果が線形モデルで使用される場合、
\(b\)
は、順序予測変数の2つの隣接するカテゴリ間の期待平均差として解釈されうる。追加のパラメータ・ベクトル、
\(\zeta\)
は、連続する予測変数のカテゴリ間の正規化距離を推定し、それによって単調効果の形状が定義される。単一の単調予測変数、
\(x\) については、オブザベーション
\(n\)
の線形予測変数の項は、以下のようになる。
\[\eta_n = b D \sum_{i = 1}^{x_n}
\zeta_i\]
パラメータ \(b\)
は任意の実数値をとることができる。一方、 \(\zeta\)
は単項式であり、以下を満たすことを意味する。 $_i $ と \(\sum_{i = 1}^D \zeta_i = 1\) を満たし、
\(D\) は \(\zeta\) の要素数である。 同様に、 \(D\)
はカテゴリ数(またはデータ中の最大整数)から1を引いた値である。なぜなら、表記を簡単にするためにカテゴリを0から数え始めるからである。
シンプルな単調モデル
単調効果の主な用途は、連続または無秩序なカテゴリ予測因子として誤って扱うことなく、この方法でモデル化することができる順序予測因子である。例えば、心理学では、この種のデータはLikertスケール項目の形で遍在しており、この仮定を検証することなく、便宜上連続であるとして扱われることがよくある。例えば、0から100までの任意の尺度で測定された年収（ドル）と生活満足度の関係に興味があるとする。通常、人々は正確な収入を尋ねられることはない。例えば、「2万円以下」、「2万円以上4万円以下」、「4万円以上10万円以下」、「10万円以上」というように、ある一定の階級に自分をランク付けするよう求められるのである。ここでは、説明のためにいくつかのシミュレーションデータを使用する。
income_options <- c("below_20", "20_to_40", "40_to_100", "greater_100")
income <- factor(sample(income_options, 100, TRUE),
levels = income_options, ordered = TRUE)
mean_ls <- c(30, 60, 70, 75)
ls <- mean_ls[income] + rnorm(100, sd = 7)
dat <- data.frame(income, ls)
次に、income
をモノトニック効果としてモデル化したデータの分析を進める。
fit1 <- brm(ls ~ mo(income), data = dat)
要約メソッドで得られるのは
Family: gaussian
Links: mu = identity; sigma = identity
Formula: ls ~ mo(income)
Data: dat (Number of observations: 100)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 31.78 1.39 29.12 34.56 1.00 2517 2311
moincome 14.52 0.63 13.28 15.76 1.00 2601 2365
Simplex Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
moincome1[1] 0.64 0.04 0.57 0.72 1.00 2728 2527
moincome1[2] 0.19 0.05 0.10 0.29 1.00 3061 1831
moincome1[3] 0.16 0.04 0.08 0.25 1.00 3492 2431
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 6.80 0.50 5.91 7.84 1.00 3484 2642
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
plot(fit1, variable = "simo", regex = TRUE)
plot(conditional_effects(fit1))
income
のシンプレクスパラメータの分布を見ると、plot
法で示されるように、最初の2つのカテゴリ間で最大の差（最小と最大のカテゴリ間の差の約70％）があることがわかる。
ここで、単調モデルと2つの一般的な代替モデルを比較してみよう。(a)
income が連続的であると仮定する。
dat$income_num <- as.numeric(dat$income)
fit2 <- brm(ls ~ income_num, data = dat)
Family: gaussian
Links: mu = identity; sigma = identity
Formula: ls ~ income_num
Data: dat (Number of observations: 100)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 24.06 2.09 19.91 28.17 1.00 3819 3031
income_num 13.77 0.75 12.26 15.25 1.00 3987 2906
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 8.92 0.64 7.74 10.28 1.00 3180 2706
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
または (b) income、順不同の因子であるとする。
contrasts(dat$income) <- contr.treatment(4)
fit3 <- brm(ls ~ income, data = dat)
Family: gaussian
Links: mu = identity; sigma = identity
Formula: ls ~ income
Data: dat (Number of observations: 100)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 31.58 1.37 28.92 34.29 1.00 3167 2589
income2 28.30 1.91 24.50 31.96 1.00 3609 3192
income3 36.62 2.09 32.58 40.87 1.00 3322 3053
income4 43.87 1.88 40.18 47.54 1.00 3521 2960
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 6.82 0.49 5.95 7.84 1.00 3779 3107
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
3つのモデルの適合度は、リーブワンアウトクロスバリデーション（leave-one-out
cross-validation）を用いて簡単に比較することができる。
Output of model 'fit1':
Computed from 4000 by 100 log-likelihood matrix
Estimate SE
elpd_loo -335.7 6.6
p_loo 4.8 0.7
looic 671.5 13.1
------
Monte Carlo SE of elpd_loo is 0.0.
All Pareto k estimates are good (k < 0.5).
See help('pareto-k-diagnostic') for details.
Output of model 'fit2':
Computed from 4000 by 100 log-likelihood matrix
Estimate SE
elpd_loo -361.8 6.6
p_loo 2.6 0.4
looic 723.5 13.1
------
Monte Carlo SE of elpd_loo is 0.0.
All Pareto k estimates are good (k < 0.5).
See help('pareto-k-diagnostic') for details.
Output of model 'fit3':
Computed from 4000 by 100 log-likelihood matrix
Estimate SE
elpd_loo -335.6 6.5
p_loo 4.7 0.7
looic 671.2 13.0
------
Monte Carlo SE of elpd_loo is 0.0.
All Pareto k estimates are good (k < 0.5).
See help('pareto-k-diagnostic') for details.
Model comparisons:
elpd_diff se_diff
fit3 0.0 0.0
fit1 -0.2 0.2
fit2 -26.2 6.3
単調モデルは連続モデルよりも適合度が高い。これは、income
と ls
の関係が非線形であることを考えると、驚くにはあたらない。
この例では、単調モデルと非順序因子モデルがほぼ同じ適合度であるが、これは他のデータセットではそうではないだろう。
事前分布の設定
これまでの単調モデルでは、隣接するカテゴリ間の差はすべて同じ事前分布を持つ、あるいは正しく定式化されている、と暗黙のうちに仮定してきた。以下では、この仮定をどのように変更するかを示したい。シンプレックス・パラメータの正規の事前分布はディリクレ分布で、ベータ分布の多変量一般化である。
これは、すべての有効なシンプレックス（すなわち、, $_i $ と \(\sum_{i = 1}^D \zeta_i = 1\)
）では0でなく、それ以外では0である。ディリクレ事前分布は、 \(\zeta\) と同じ長さの1つのパラメータ \(\alpha\) を持つ。 \(\alpha_i\) が高いほど、 \(\zeta_i\)
の値が高くなるa-priori確率が高い。
データを見る前に、同じ金額の追加的なお金は、一般にお金を持っていない人々にとってより重要であると予想していたとする。これは、
\(\zeta_1\)
のa-priori値（’20以下’と’20_to_40’の差）が高く、したがって、 \(\alpha_1\) の値も高くなる。私たちは、 \(\alpha_1 = 2\) と \(\alpha_2 = \alpha_3 =
1\)、後者はデフォルト値として \(\alpha\)
を選んだ。モデルを適合させるためには、次のように書く。
prior4 <- prior(dirichlet(c(2, 1, 1)), class = "simo", coef = "moincome1")
fit4 <- brm(ls ~ mo(income), data = dat,
prior = prior4, sample_prior = TRUE)
"moincome1" の末尾にある 1
は、初めて単調効果を扱うときには奇妙に見えるかもしれない。しかし、複数の単調変数の相互作用がモデルに含まれている場合、1つの単調項が複数の単調パラメータに関連することがあるので、これは必要である。
Family: gaussian
Links: mu = identity; sigma = identity
Formula: ls ~ mo(income)
Data: dat (Number of observations: 100)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 31.77 1.33 29.17 34.30 1.00 2643 2568
moincome 14.52 0.61 13.33 15.73 1.00 2687 2464
Simplex Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
moincome1[1] 0.65 0.04 0.57 0.72 1.00 3260 2524
moincome1[2] 0.19 0.05 0.09 0.29 1.00 3339 1933
moincome1[3] 0.16 0.04 0.07 0.25 1.00 3380 2512
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 6.81 0.51 5.93 7.90 1.00 3154 2872
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
sample_prior = TRUE
を使って、simo_moincome1
の事前分布からのドローも取得し、可視化できるようにした。
plot(fit4, variable = "prior_simo", regex = TRUE, N = 3)
プロットでわかるように simo_moincome1 [1]
は平均して2倍高い。 simo_moincome1 [2] と
simo_moincome1 [3] に設定した結果、 \(\alpha_1\)。
モノトーン変数の相互作用のモデル化
さらに、参加者に年齢を聞いたとする。
dat$age <- rnorm(100, mean = 40, sd = 10)
私たちは年齢の主効果だけでなく、所得と年齢の交互作用にも興味がある。単調変数との交互作用は、*
演算子を用いて通常の方法で指定することができる。
fit5 <- brm(ls ~ mo(income)*age, data = dat)
Family: gaussian
Links: mu = identity; sigma = identity
Formula: ls ~ mo(income) * age
Data: dat (Number of observations: 100)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 33.15 4.43 24.73 41.84 1.00 1453 2040
age -0.03 0.11 -0.25 0.17 1.00 1407 1888
moincome 12.57 2.16 8.58 17.20 1.01 888 1507
moincome:age 0.05 0.05 -0.07 0.14 1.00 875 1400
Simplex Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
moincome1[1] 0.70 0.08 0.56 0.88 1.00 1336 1528
moincome1[2] 0.15 0.07 0.02 0.29 1.00 1850 1786
moincome1[3] 0.14 0.06 0.02 0.27 1.00 2075 1449
moincome:age1[1] 0.32 0.22 0.01 0.79 1.00 2363 1611
moincome:age1[2] 0.38 0.23 0.02 0.84 1.00 2781 2092
moincome:age1[3] 0.30 0.21 0.01 0.79 1.00 2517 2333
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 6.85 0.52 5.91 7.96 1.00 3377 2665
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
conditional_effects(fit5, "income:age")
単調なグループレベル効果のモデル化
サンプルデータの100人は10都市から集められ、1都市あたり10人であるとする。そこで、データに
city という識別子を追加し、ls
に都市に関連したバリエーションを追加する。
dat$city <- rep(1:10, each = 10)
var_city <- rnorm(10, sd = 10)
dat$ls <- dat$ls + var_city[dat$city]
以下のコードで、切片と income
の効果が都市によって異なると仮定し、マルチレベルモデルを適合させた。
fit6 <- brm(ls ~ mo(income)*age + (mo(income) | city), data = dat)
Family: gaussian
Links: mu = identity; sigma = identity
Formula: ls ~ mo(income) * age + (mo(income) | city)
Data: dat (Number of observations: 100)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Group-Level Effects:
~city (Number of levels: 10)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 13.74 4.32 7.71 24.10 1.00 1593 2007
sd(moincome) 1.06 0.87 0.04 3.33 1.00 1769 2345
cor(Intercept,moincome) 0.15 0.53 -0.87 0.96 1.00 5215 2622
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 36.25 6.58 23.25 49.32 1.01 1223 1527
age -0.05 0.11 -0.28 0.18 1.00 1938 2221
moincome 13.10 2.40 8.76 18.16 1.00 1665 2379
moincome:age 0.04 0.06 -0.08 0.14 1.00 1610 2230
Simplex Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
moincome1[1] 0.69 0.09 0.54 0.87 1.00 2090 2236
moincome1[2] 0.17 0.07 0.02 0.30 1.00 2652 2222
moincome1[3] 0.14 0.06 0.02 0.27 1.00 3165 2019
moincome:age1[1] 0.33 0.23 0.01 0.81 1.00 3757 2202
moincome:age1[2] 0.36 0.23 0.02 0.85 1.00 3748 3014
moincome:age1[3] 0.31 0.22 0.02 0.80 1.00 3746 3119
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 6.71 0.54 5.74 7.87 1.00 4629 2877
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
を見ると、income
の効果は都市間でほとんど差がないことがわかる。今回のデータでは、データ・シミュレーションにおいて、income
が都市間で同じ効果を持つと仮定していることを考えると、これはそれほど驚くことではない。