6  第8章 消費者の選択と離散選択モデル

6.1 顕示選好の理論

Kreps (1988) “Notes on the Theory of Choice”に基づいて，単一個人意思決定モデルの基礎的な考え方を説明します。

「選択」(choice)という現象は，ある集合Aがある個人にとっての実際に選択可能な選択肢の集合であったときに，その中から個人によって1つの選択肢a \in Aが選ばれる，という現象をいいます。 ここに，「選ぶ」というのは，第三者に観察できる行為として選ぶ，ということに限って用いられます。 また，問題となっている集合Aの各要素a \in Aを「選択肢」(alternative)と呼びます。

選択現象を問題にしている選択肢aが確率的構造を持っているとき，その選択肢をとくにクジ(lottery)と呼びます。

説明するべき対象を選択肢とします。

X \supseteq A \mapsto a

ここでAは選択肢の集合で，aはその選択肢を選んだ結果となります。Xは整合性を問題にできる選択肢集合のいちばん最大のもの，となります。

X上の，この人の選択関数 C : X \subseteq A \mapsto C(A) \subseteq X が安定的に存在すると仮定します。 ここでC(A)は，Aの中からこの人が選んだ選択肢集合を表します。

Xを考察する選択行動で，選択の対象となりうる潜在的な選択肢を集めた集合を表す。 選択集合とは，あるA \in \mathcal{P}(X)が実際に選択可能な選択肢の集合であったときに，その中から1つの選択肢が選ばれる現象をいいます。 「選ばれる」といういうとき，現に選ばれた1つの選択肢以外にも「それでも同様に良かった」というものがあるとき，それらの「そのどれでも良かった」選択肢を全て集めた集合Cを，Aからの選択結果とみなします。

選択関数を定義すると，

C: \mathcal{P}(X) \to \mathcal{P}(X)が， \forall A \in \mathcal{P}(X); \ C(A) \subseteq A であるとき，X上の選択関数(choice function)である，といいます。 ここで，\mathcal{P}(X)は，Xの部分集合1つ1つを全部集めてできる冪集合(power set)のうち，Xの中から二者択一の結果だけを集めたデータを表します。

6.1.1 選好モデル

二項関係(binary relation)を定義します。 Xの（たまたま同じであることを許す）2つの要素について語られた関係が，X上の二項関係(binary relation)であるとは， X \times Xのすべての要素(x,y)について，xはyに対してその関係にあるか，またはその関係にないか，を必ず一方をいうことができる場合を言います。

つまり，選択肢集合Aに選択肢xとyがあり，その人の選択結果C(\{x,y\}) = ?を見ます。 その人が選択するものは，

\begin{aligned} C(\{x,y\}) &= \phi \quad \text{$x$も$y$も選べたのに選ばなかった。}\\ C(\{x,y\}) &= \{x\} \quad xを選び，yを選ばなかった\\ C(\{x,y\}) &= \{y\} \quad yを選び，xを選ばなかった\\ C(\{x,y\}) &= \{x,y\} \quad xもyも選んだ。 \end{aligned}

6.1.2 選好

X 上の選択関数 C: \mathcal{P}(X) \rightarrow \mathcal{P}(X) が与えられたとき，X の(それがたまたま同じであることを許す)2つの要素についての関係を，各 (x,y) \in X \times X に対して，y \not \in C(\{x,y\}) なら，x は y に対してその関係にあり，y \in C(\{x,y\}) なら x は y に対してその関係にない，というやり方で定められるとき，その関係(それは二項関係である)を，(選択関数 C:\mathcal{P}(X) \rightarrow \mathcal{P}(X) によって生成された) X 上の強選好(strict preference)といいます。

この二項関係を記号 \succ で表わし，x が y に対して強選好の関係にあるとき，x を y よりも厳密に選好する，あるいは x は y よりも厳密に選好されるといって，x \succ y と書き，x が y に対して強選好の関係にないとき，x は y よりも厳密に選好されるわけではない，といって x \not \succ y と書きます。

ある X 上の強選好 \succ は，それが定義されたもとの選択関数 C : \mathcal{P}(X) \rightarrow \mathcal{P}(X)を明示する必要があるときは，\succ _cと書く。

x,yの二者択一の結果，C(\{ x,y \})を， \begin{aligned} & x \not \succ y \Longleftrightarrow y \in C(\{x,y\})\leftarrow \text{$x$は$y$より厳密に選好されるわけではない。}\\ & x \succ y \Longleftrightarrow y \not \in C(\{x,y\})\leftarrow xはyより厳密に選好される。 \end{aligned}

のように書いたものを，この人の(強)選好という。

つぎに、X 上の強選好 \succ が与えられたとき、各A \in \mathcal{P}(X)に対し、 C(A, \succ) = \{ x \in A \mid \nexists y \in A \quad \text{s.t.} \quad y \succ x \} を対応づける（集合体集合の）関数 C(\cdot , \succ) : \mathcal{P}(X) \rightarrow \mathcal{P}(X) を、X 上の選択の選好モデル(preference model)といいます。

まとめると、 C : \mathcal{P}(X) \rightarrow C(A) \in \mathcal{P}(X) のうち、「二者択一の結果を、Xの各ペア(x,y) ごとに、x \succ y か x \not \succ y かを書いて、書いたものをX上の（この人の）強選好」といいます。 選好モデルとは、 C : \mathcal{P}(X) ∋ A \rightarrow C(A) \in \mathcal{P}(X) を \begin{aligned} & C(\cdot, \succ) : \mathcal{P}(X) ∋ A \rightarrow C(A, \succ) \in \mathcal{P}(X)\\ & \text{such that} \quad C(A, \succ) = \{ x \in A \mid \nexists y \in A \quad \text{s.t.} \quad y \succ x \} \end{aligned} となり、この選好\succをもつ人が、実際にAに直面したときに選ぶと予測するということを意味します。 これはC(\cdot , \succ)のデータから予測されます。

C(\cdot , \succ) : \mathcal{P}(X) ∋ A \rightarrow C(A, \succ) \in \mathcal{P}(X)は、C : \mathcal{P}(X) ∋ A \rightarrow C(A) \in \mathcal{P}(X)のうち、|X|^2個のデータで、C : \mathcal{P}(X) ∋ A \rightarrow C(A) \in \mathcal{P}(X)全体の2^{|X|}個の値を予測しているのです。

合理性モデルでは、「もし・・・なら、\forall A \in \mathcal{P}; C(A, \succ) = C(A)となります。

ふりかえると、

• 説明する対象は、選択関数 C : \mathcal{P}(X) \in A \mapsto C(A) \in \mathcal{P}(X) である。
• データを取ってくるパラメータとしての選好\succ (つまり、X上のすべての二者択一の結果)を考える。
• 選択を説明するモデルは選好モデルC(A, \succ) = \{x \in A \mid \nexists y \in A \quad \text{s.t.} \quad y \succ x \}
• 選択理論とは、「C: \mathcal{P}(X) \rightarrow \mathcal{P}(X) が・・・の条件を満たすとき、Cのうち××から作った$C(A, ××) = C(A)」が成り立つような「･･･」と「××」を考えます。

ただ、選好\succとしてどんなものでも認めてしまうと、

\exists \succ, \exists A \quad \text{s.t.} \quad C(A, \succ) = \emptyset となる場合があります。 たとえば、A = X = \{x, y, z\}で、 C(\{x,y\}) = \{x\}、 C(\{y,z\}) = \{y\}、 C(\{x,z\}) = \{z\} かつ、 C(A) \not = \emptyset. 場合、C(\{x,y,z\}) = \emptyset となります。

つまり、

• z \succ x より x \not \in C(\{x,y, z\}, \succ)
• x \succ y より y \not \in C(\{x,y, z\}, \succ)
• y \succ z より z \not \in C(\{x,y, z\}, \succ)

すなわち、C(\{x,y,z\}, \succ) = \emptyset となります。

て、次のように選択関数 C_{\succ} : \mathcal{P}(X) \rightarrow \mathcal{P}(X) を定めます。

では、どのような二者択一行動を示す人なら、その人の選好モデルは完結性があるのかを考えます。 次の2つの条件が満たされるとき、選好\succは完結性(completeness)を持つといいます。

1. 反対称性(asymmetry) : \forall x,y \in X;, 「x \succ y かつ y \succ x」ということはない。つまり二者択一行動がC(\{x,y\}) = \{\emptyset\}ということがない。
2. 負推移性(negative transitivity) : \forall x,y \in X;, 「もし x \not \succ y かつ y \not \succ z」ならば「x \not \succ z」である。 つまり二者択一行動としては、もしy \in C(\{x,y\}) かつ z \in C(\{y,z\}) ならば、z \in C(\{x,z\}) である。

これより、X 上の強選好 \succ は、次のように定義できます。

Important反対称性・負推移性

\nexists (x,y) \in X \times X; \quad \text{s.t} \quad x \succ y, \quad y \succ x であるとき、反対称的(assymmetric)である、あるいは反対称性を満たすといい、 \forall x,y,z \in X; \quad \text{if} \quad x \not \succ y, \quad y \not \succ z, \quad \text{then} \quad x \not \succ z であるとき、負推移的(negatively transitive)である、あるいは負推移性を満たすといいます

Important選好関係

X 上の強選好 \succ は、反対称的かつ負推移的であるとき、X上の選好関係(preference relation)であるといいます。

選好\succが反対称的かつ負推移的であるとき、推移的(acyclic)であるといい、 x_1 \succ x_2 かつ x_2 \succ x_3 かつ \cdots かつ x_{n-1} \succ x_n なら x_1 \not = x_n となります。

6.1.3 Xが有限集合の場合

X 上の強選好 \succ が \forall x \in X; \quad x \not \succ x であるとき、非反射的(irreflexive)である、あるいは非反射性を満たすといいます。

Note推移性

X 上の強選好 \succ は、 \forall x, y, z \in X; \quad \text{if} \quad x \succ y, \ y \succ x, \quad \text{then} \quad x \succ z であるとき、推移的(transitive)である、あるいは推移性を満たすといいます。

証明は省略（背理法により示せます）。

Note非周期性

次に、X 上の強選好 \succ は、

\begin{aligned} &\forall x_1, x_2, \cdots, x_n \in X;\\ &\quad \text{if} \quad x_1 \succ x_2, \ x_2 \succ x_3, \ \dots, \ x_{n-1} \succ x_n, \quad \text{then} \quad x_1 \not = x_n \end{aligned} であるとき、非周期的(acyclic)である、あるいは非周期性を満たすといいます。

証明は省略（帰納法により示せます）。

Important

X 上の強選好 \succ は、X上の選好関係ならば、非反射的かつ推移的かつ非周期的である。

X上の選択関数 C : \mathcal{P}(X) \rightarrow \mathcal{P}(X)は、

\begin{aligned} \forall A \in \mathcal{P}(X); \quad \text{if} \quad A \not = \emptyset \quad \text{then} \quad C(A) \not = \emptyset \end{aligned} であるとき、非空(nonempty)であるといいます。

これらの定理と定義から、次の定理が示せます。

Important

X を有限集合(|X| < \infty)とします。 X 上の強選好\succは、X上の(反対称的と負推移性を満たす)選好関係ならば、A \not = \emptysetのとき、C(\cdot, \succ): \mathcal{P}(X) \rightarrow \mathcal{P}(X)は非空となります。

が選好関係であるとき、Xが有限集合であれば、任意のA \in \mathcal{P}(X)に対し、C(A, \succ) \not = \emptyset となります。

ここまでの説明により、ある個人の選択行動を説明するために必要な条件は、 ある個人が(有限の)選択肢の中から選び取るとき、その人の選好が反対称的かつ負推移的であれば、選択肢集合が非空なら、その人がもつ選択関数は非空である、ということがわかりました。

6.2 顕示選考理論

選択肢集合Xが有限であるとします。 すると、X上の強選好\succがX上の選好関係であるならば、 C(\cdot , \succ) : \mathcal{P}(X) \rightarrow \mathcal{P}(X)は非空まつハウタッカーの公理(Houthakker’s axiom)を満たします。

Importantハウタッカーの公理

X上の選択関数C : \mathcal{P}(X) \rightarrow \mathcal{P}(X)は、 \begin{aligned} & \forall x,y \in X, \quad \forall A,B \in \mathcal{P}(X); \\ & \quad \text{if} \quad x,y \in A \cap B, \quad x \in C(A), \quad y \in C(B)\\ & \quad \text{then} \quad x \in C(B), \quad y \in C(A) \end{aligned} であるとき、ハウタッカーの公理を満たす、といいます。

具体例で説明します。 たとえば、あなたの買い物の選択肢として

• A=\{\text{みかん, リンゴ, バナナ}\}、
• B=\{\text{バナナ, リンゴ, ぶどう}\}

があるとします。

このとき、あなたは

• A の中からリンゴを選び（C(A)=\{\text{リンゴ}\}）、
• B の中からバナナを選ぶ（C(B)=\{\text{バナナ}\}）

とします。

リンゴとバナナはいずれも A\cap B に含まれています。 ハウタッカーの公理は、この状況のもとでは、 選択肢集合が B に変わってもリンゴが選ばれ、 選択肢集合が A に変わってもバナナが選ばれなければならない、 ということを要請します。

つまり、共通して含まれる選択肢については、 集合ごとに選択の優劣が入れ替わるような行動は許されない、 という選択の一貫性を表しています。

Noteセンの\alpha規則性

X上の選択関数C : \mathcal{P}(X) \rightarrow \mathcal{P}(X)は、

\begin{aligned} & \forall x \in X, \quad \forall A,B \in \mathcal{P}(X);\\ & \quad \text{if} \quad x \in B \subset A, \quad x \in C(A),\\ & \quad \text{then} \quad x \in C(B) \end{aligned}

であるとき、センの\alpha規則性(Sen’s property \alpha)を満たす、といいます。

Tipセンの\alpha規則性の具体例

あなたの選択肢の全集合をXとし、 A=\{\text{みかん, リンゴ, バナナ}\}、 B=\{\text{リンゴ, バナナ}\} とします。このとき B \subset A が成り立っています。

いま、あなたが集合 A からリンゴを選んでいる、すなわち C(A)=\{\text{リンゴ}\} であるとします。リンゴは B にも含まれています。

センの \alpha 規則性は、この状況において、選択肢集合が A からその部分集合である B に縮小されても、リンゴは引き続き選ばれなければならない、すなわち C(B)=\{\text{リンゴ}\} であることを要請します。

言い換えると、ある選択肢がより大きな集合の中で選ばれていたのであれば、その選択肢を排除しない範囲で選択肢集合を減らしても、その選択肢が選ばれなくなることは許されない、という性質を表しています。

Noteセンの\beta規則性

X上の選択関数C : \mathcal{P}(X) \rightarrow \mathcal{P}(X)は、 \begin{aligned} &\forall x,y \in X, \quad \forall A,B \in \mathcal{P}(X); \\ &\quad \text{if} \quad x,y \in C(A), \quad A \subset B, \quad y \in C(B), \\ &\quad \text{then} \quad x \in C(B) \end{aligned} であるとき、センの\beta規則性(Sen’s property \beta)を満たす、といいます。

Tipセンの\beta規則性の具体例

選択肢の全集合を X とし、 A=\{\text{リンゴ, バナナ}\}、 B=\{\text{みかん, リンゴ, バナナ}\} とします。このとき A \subset B が成り立っています。

いま、集合 A からはリンゴとバナナの両方を選んでいる、すなわち C(A)=\{\text{リンゴ, バナナ}\} とします。

次に、選択肢が増えた集合 B からはバナナを選んでいる、すなわち C(B)=\{\text{バナナ}\} であるとします。

センの \beta 規則性は、この状況のもとでは、集合 A で同時に選ばれていたリンゴとバナナのうち、集合 B に移ってもバナナが選ばれているのであれば、リンゴもまた集合 B から選ばれなければならない、ということを要請します。

したがって、この規則性を満たす選択関数では、 C(B) は \{\text{バナナ}\} ではなく、少なくともリンゴを含み、たとえば C(B)=\{\text{リンゴ, バナナ}\} でなければならない、ということになります。

この例は、集合を拡大したときに、もともと同時に選ばれていた選択肢の一部だけが残るような振る舞いを排除する性質を具体的に示しています。

これらの公理と規則性から次の定理が示せます。

Important定理

X上の選択関数C : \mathcal{P}(X) \rightarrow \mathcal{P}(X)が非空とします。 このとき、C: \mathcal{P}(X) \rightarrow \mathcal{P}(X) がハウタッカーの公理を満たすための必要十分条件は、それがセンの\alpha規則性と\beta規則性を満たすことです。

具体例で説明してみます。

X=\{\text{みかん, リンゴ, バナナ, ぶどう}\} とします。あなたの好みが

\text{リンゴ} \sim \text{バナナ} \succ \text{みかん} \succ \text{ぶどう}

だとします。 つまりリンゴとバナナは同程度に最も好ましく、次にみかん、最後にぶどう、という順番です。

この好みに基づく選択関数 C を、「与えられた集合の中の“最も好ましいもの”をすべて選ぶ」と定めます。 つまり任意の A\subseteq X に対して

• A にリンゴまたはバナナが含まれていれば、C(A) = A\cap\{\text{リンゴ, バナナ}\}
• そうでなければ（リンゴもバナナも無いなら）、残りの中で最上位のもの（ここではみかん、なければぶどう）を選ぶ

というルールです。たとえば

• C(\{\text{みかん, リンゴ, バナナ}\})=\{\text{リンゴ, バナナ}\},
• C(\{\text{バナナ, ぶどう}\})=\{\text{バナナ}\},
• C(\{\text{みかん, ぶどう}\})=\{\text{みかん}\}

のようになります。

この選択関数C はセンの \alpha 規則性を満たします。理由は、ある集合 A で選ばれた x は A の中で「最上位」なので、x を含む部分集合 BA に選択肢を減らしても、x より上位のものが新たに現れることはなく、x は依然として最上位（最大元）だからです。

また 規則性も満たします。実際、ある集合 A で x,yC(A) ということは、x と y がともに A の中で最上位（この例ではリンゴとバナナのように同率首位）ということです。集合を広げた BA でも、もし y が依然として選ばれているなら（yC(B)）、y より上位の選択肢は B に存在しないので、同率首位である x も排除されず xC(B) となります。

ここから、ハウタッカーの公理を満たしていることを、同じ果物で具体的に確かめます。たとえば

A={}, B={}

とおくと、

C(A)={}, C(B)={}

です。いま x=, y= とすると、確かに x,yAB かつ xC(A)、yC(B) が成り立ちます。そして結論として要求される xC(B)、yC(A) も実際に成り立っています（どちらの集合でもリンゴとバナナが選ばれているため）。

また、X上の選択関数C : \mathcal{P}(X) \rightarrow \mathcal{P}(X) が非空( A\not = \emptysetならC(A) \not = \emptyset )かつハウタッカーの公理を満たすならば、\succ_cはX上の選好関係\succであり(反対称性と負推移性を満たす)であり、しかもA \not = \emptysetならば、

\forall A \in \mathcal{P}(X); \quad C(A, \succ_c) = C(A) が成り立ちます。

pacman::p_load(tidyverse, readxl, modelsummary, mfx, DescTools, mlogit)

The downloaded binary packages are in
    /var/folders/3k/7grsq_p50s341kh1q563qy040000gn/T//RtmprL1cVJ/downloaded_packages

The downloaded binary packages are in
    /var/folders/3k/7grsq_p50s341kh1q563qy040000gn/T//RtmprL1cVJ/downloaded_packages

The downloaded binary packages are in
    /var/folders/3k/7grsq_p50s341kh1q563qy040000gn/T//RtmprL1cVJ/downloaded_packages

# コード8-1
df1 <- data.frame(
    Y = c(0, 0, 0, 0, 0, 1, 1, 1, 1),
    X = c(3.4, 5.22, 7.06, 2.81, 4.11, 10.34, 13.67, 15.99, 9.09)
  )
lpm1 <- lm(Y ~ X, data = df1)
pred_lpm1 <- predict(lpm1)
pred_lpm1

           1            2            3            4            5            6 
-0.006342894  0.173357681  0.355032986 -0.064597476  0.063760077  0.678888967 
           7            8            9 
 1.007681776  1.236750640  0.555468243 

# コード8-2
choice_df <- readxl::read_xlsx("data/choice_data.xlsx")
head(choice_df)

# コード8-3
probit1 <- glm(y1 ~ p1 + a1,
               family = binomial(link = probit),
               data = choice_df)
summary(probit1)

Call:
glm(formula = y1 ~ p1 + a1, family = binomial(link = probit), 
    data = choice_df)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  2.04117    0.22350   9.133  < 2e-16 ***
p1          -0.24720    0.02213 -11.171  < 2e-16 ***
a1           0.62763    0.08656   7.251 4.13e-13 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1384.5  on 999  degrees of freedom
Residual deviance: 1205.5  on 997  degrees of freedom
AIC: 1211.5

Number of Fisher Scoring iterations: 3

# コード8-4
ols1 <- lm(y1 ~ p1 + a1,
           data = choice_df)
logit1 <- glm(y1 ~ p1 + a1,
              family = binomial(link = logit),
              data = choice_df)
models <- list()
models[["Linear probability model"]] <- ols1
models[["Probit model"]] <- probit1
models[["Logit model"]] <- logit1
modelsummary(models,
             title = "モデル比較",
             notes = "Values in [ ] show robust standard errors",
             stars = TRUE,
             statistic = "std.error",
             vcov = "robust",
             gof_map = "nobs")

モデル比較

	Linear probability model	Probit model	Logit model
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001
Values in [ ] show robust standard errors
(Intercept)	1.196***	2.041***	3.332***
	(0.065)	(0.212)	(0.361)
p1	-0.085***	-0.247***	-0.404***
	(0.006)	(0.021)	(0.036)
a1	0.221***	0.628***	1.035***
	(0.029)	(0.087)	(0.145)
Num.Obs.	1000	1000	1000

# 適合度指標において，サンプルサイズ ("nobs") のみ表示するという指示

# コード8-5
# Average Marginal Effects（限界効果の平均）
probit1_ame <- probitmfx(y1 ~ p1 + a1,
                         data = choice_df, atmean = FALSE)
probit1_ame

Call:
probitmfx(formula = y1 ~ p1 + a1, data = choice_df, atmean = FALSE)

Marginal Effects:
        dF/dx  Std. Err.        z     P>|z|    
p1 -0.0847786  0.0060787 -13.9469 < 2.2e-16 ***
a1  0.2188866  0.0289816   7.5526 4.266e-14 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

dF/dx is for discrete change for the following variables:

[1] "a1"

# コード8-7
# Marginal Effects at Mean（平均値における限界効果）
probit1_mem <- probitmfx(y1 ~ p1 + a1,
                         data = choice_df, atmean = TRUE)
probit1_mem

Call:
probitmfx(formula = y1 ~ p1 + a1, data = choice_df, atmean = TRUE)

Marginal Effects:
        dF/dx  Std. Err.        z     P>|z|    
p1 -0.0984235  0.0088074 -11.1750 < 2.2e-16 ***
a1  0.2446609  0.0324366   7.5427 4.602e-14 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

dF/dx is for discrete change for the following variables:

[1] "a1"

# コード8-8
marginal <- c(
  "p1 marginal" = "p1",
  "a1 marginal" = "a1"
)

marg_model <- list(
  "(1) Probit_AME" <- probit1_ame,
  "(2) Probit_MEM" <- probit1_mem
)

marg_sum <- modelsummary(
  marg_model,
  statistic = "std.error",
  coef_map = marginal,
  stars = TRUE,
  shape = term:component ~ model,
  add_rows = data.frame("Marginal effect type",
                        "Average Marginal Effect",
                        "Marginal Effect at Mean"),
  title = "限界効果サマリー",
  gof_map = "nobs"
)
marg_sum

限界効果サマリー

	(1)	(2)
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001
p1	-0.085***	-0.098***
	(0.006)	(0.009)
a1	0.219***	0.245***
	(0.029)	(0.032)
Num.Obs.	1000	1000
Marginal effect type	Average Marginal Effect	Marginal Effect at Mean

# コード8-9
DescTools::PseudoR2(probit1)

 McFadden 
0.1293345 

# コード8-10

# 価格差変数作成
choice_df <- choice_df %>%
  mutate(p_ratio = p1 - p2)

probit2 <- glm(y1 ~ p_ratio + a1 + a2,
               family = binomial(link = probit),
               data = choice_df)
summary(probit2)

Call:
glm(formula = y1 ~ p_ratio + a1 + a2, family = binomial(link = probit), 
    data = choice_df)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -0.54634    0.07756  -7.044 1.87e-12 ***
p_ratio     -0.22180    0.01524 -14.559  < 2e-16 ***
a1           0.66314    0.09162   7.238 4.55e-13 ***
a2          -0.32291    0.09884  -3.267  0.00109 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1384.5  on 999  degrees of freedom
Residual deviance: 1066.9  on 996  degrees of freedom
AIC: 1074.9

Number of Fisher Scoring iterations: 4

# コード8-11
DescTools::PseudoR2(probit2)

McFadden 
0.229384 

# コード8-12
probit2_ame <- probitmfx(y1 ~ p_ratio + a1 + a2,
                         data = choice_df, atmean = FALSE)
probit2_ame

Call:
probitmfx(formula = y1 ~ p_ratio + a1 + a2, data = choice_df, 
    atmean = FALSE)

Marginal Effects:
             dF/dx  Std. Err.        z     P>|z|    
p_ratio -0.0669199  0.0029405 -22.7581 < 2.2e-16 ***
a1       0.2036848  0.0270475   7.5306 5.049e-14 ***
a2      -0.0969891  0.0291848  -3.3233 0.0008897 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

dF/dx is for discrete change for the following variables:

[1] "a1" "a2"

# コード8-13
data(Cracker, package = "mlogit")
head(Cracker, n = 20)

# コード8-14
cracker <- mlogit.data(Cracker, choice = "choice", shape = "wide",
                       varying = c(2:13))
head(cracker, n = 20)

~~~~~~~
 first 20 observations out of 13168 
~~~~~~~
   id choice      alt disp feat price chid    idx
1   1  FALSE  keebler    0    0    88    1 1:bler
2   1   TRUE  nabisco    0    0   120    1 1:isco
3   1  FALSE  private    0    0    71    1 1:vate
4   1  FALSE sunshine    0    0    98    1 1:hine
5   1  FALSE  keebler    0    0   109    2 2:bler
6   1   TRUE  nabisco    0    0    99    2 2:isco
7   1  FALSE  private    0    0    71    2 2:vate
8   1  FALSE sunshine    0    0    99    2 2:hine
9   1  FALSE  keebler    0    0   109    3 3:bler
10  1  FALSE  nabisco    0    0   109    3 3:isco
11  1  FALSE  private    0    0    78    3 3:vate
12  1   TRUE sunshine    1    0    49    3 3:hine
13  1  FALSE  keebler    0    0   109    4 4:bler
14  1   TRUE  nabisco    0    0    89    4 4:isco
15  1  FALSE  private    0    0    78    4 4:vate
16  1  FALSE sunshine    0    0   103    4 4:hine
17  1  FALSE  keebler    0    0   109    5 5:bler
18  1   TRUE  nabisco    0    0   119    5 5:isco
19  1  FALSE  private    0    0    64    5 5:vate
20  1  FALSE sunshine    0    0   109    5 5:hine

~~~ indexes ~~~~
   chid      alt
1     1  keebler
2     1  nabisco
3     1  private
4     1 sunshine
5     2  keebler
6     2  nabisco
7     2  private
8     2 sunshine
9     3  keebler
10    3  nabisco
11    3  private
12    3 sunshine
13    4  keebler
14    4  nabisco
15    4  private
16    4 sunshine
17    5  keebler
18    5  nabisco
19    5  private
20    5 sunshine
indexes:  1, 2 

# コード8-15
ml_cracker <- mlogit(choice ~ price | 1 | disp + feat, probit = FALSE,
                     data = cracker)

summary(ml_cracker)

Call:
mlogit(formula = choice ~ price | 1 | disp + feat, data = cracker, 
    probit = FALSE, method = "nr")

Frequencies of alternatives:choice
 keebler  nabisco  private sunshine 
0.068651 0.544350 0.314399 0.072600 

nr method
5 iterations, 0h:0m:0s 
g'(-H)^-1g = 0.000258 
successive function values within tolerance limits 

Coefficients :
                       Estimate Std. Error  z-value  Pr(>|z|)    
(Intercept):nabisco   2.0011226  0.0831336  24.0712 < 2.2e-16 ***
(Intercept):private   0.3193793  0.1238052   2.5797 0.0098888 ** 
(Intercept):sunshine -0.5435987  0.1139407  -4.7709 1.834e-06 ***
price                -0.0300966  0.0021082 -14.2761 < 2.2e-16 ***
disp:keebler          0.2999316  0.2070369   1.4487 0.1474251    
disp:nabisco          0.1011111  0.0773633   1.3070 0.1912249    
disp:private         -0.2244555  0.1495236  -1.5011 0.1333199    
disp:sunshine         0.4818670  0.1672400   2.8813 0.0039605 ** 
feat:keebler          0.6678542  0.2581732   2.5868 0.0096859 ** 
feat:nabisco          0.6047773  0.1404077   4.3073 1.653e-05 ***
feat:private          0.1726981  0.2004446   0.8616 0.3889213    
feat:sunshine         0.8304895  0.2340938   3.5477 0.0003886 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Log-Likelihood: -3333.8
McFadden R^2:  0.052798 
Likelihood ratio test : chisq = 371.66 (p.value = < 2.22e-16)

# コード8-16
# パラメータの推定値（beta）の抽出
beta <- probit2$coefficients

# パラメータbeta と説明変数の線形結合を作成（beta_0 は定数項なので1 をかける）
est_probit2 <- cbind(1, choice_df$p_ratio, choice_df$a1,
                     choice_df$a2) %*% beta

# pnorm によって標準正規分布の分布関数による計算を実行する。
# 製品1(c1) を選ぶ確率の予測値
pred_probit_c1 <- pnorm(est_probit2)
mean(pred_probit_c1)

[1] 0.4789384

# コード8-17
# ロジットモデルの場合
logit2 <- glm(y1 ~ p_ratio + a1 + a2,
              family = binomial(link = logit),
              data = choice_df)
beta_logit <- logit2$coefficients
est_logit2 <- cbind(1, choice_df$p_ratio, choice_df$a1,
                    choice_df$a2) %*% beta_logit
pred_logit_c1 <- (exp(est_logit2)) / (1 + exp(est_logit2))
mean(pred_logit_c1)

[1] 0.479

# コード8-18
choice_df_v <- choice_df %>%
  mutate(p1_v = p1 - 1,
         p_ratio_v = p1_v - p2)
est_probit_v <- cbind(1, choice_df_v$p_ratio_v, choice_df_v$a1,
                      choice_df_v$a2) %*% beta

# pnorm によって標準正規分布の分布関数による計算を実行する。
# 製品1(c1) を選ぶ確率の予測値
pred_probit_v <- pnorm(est_probit_v)
mean(pred_probit_v)

[1] 0.5459071