隨機係數異質性多項選擇模型(Mixed-Multinominal Discrete Choice Models )
梁雪富 (Alfo Liang)
2022-03-27
- 1 簡介
- 2 方法論
- 3 Ex1: 婦女勞動力參與實證研究
- 4 Ex2: 博士論文發表量的性別差異: Sex Differences in Science: Poisson model
- 5 Health Demand: Ordered probit model
- 6 Panel Data: Poisson Model with Random Parameters
- 7 面板數據隨機效應的隨機參數模型(Panel data with random effects and random parameters models)
- 8 具有觀察到的異質性的隨機參數模型
- 9 條件平均數
- 10 補償變量 (compensenting variations)
1 簡介
離散選擇模型(Discrete choice models )或定性選擇模型(qualitative choice models ) 旨在描述、解釋和預測兩個或多個離散選擇之間的選擇，例如是否購買汽車，在供暖系統之間進行選擇，或在其他應用程序中選擇不同職業（Train 2009）。 這些模型的一個主要特點是，每個人的選擇可以在效用最大化行為的假設下推導出來。本質上，在給定個人的一組屬性、替代方案的屬性以及 一個隨機項，決策者選擇具有最高滿意度效用的方案，此隨機項捕捉影響效用但不包括在模型內的因素。 例如，我從備選方案 \(j (j = 1,..., J )\) 獲得的效用。 效用可以寫成 \(U_{ij} = V_{ij} + \epsilon_{ij}\),其中 \(V_{ij}\) 是效用的確定性部分，取決於可觀察的特徵(observable characteristics ), 而 \(\epsilon_{ij}\) 是隨機部分(random terms)。 如果所有的 \(U_{ij} > U_{ik}, \, k \not= j\)， 則個人 i 選擇備選方案 j。 一旦指定了 \(\epsilon_{ij}\) 的分佈和可觀察輸出決策的性質，就可以使用概率模型來估計行為過程的參數和選擇某個替代方案的概率。
這種稱為隨機效用模型 (random utility model RUM) 的公式一直是推導條件 logit 模型 (McFadden 1974) 的標準方法。 然而，當觀察到的輸出是二進制(binary)、有序(ordered)或計數(count)數據時，它也可以用來描述解釋變量的關係和決策。對於這些類型的因變量(dependent variables)，可以應用傳統的二進制 logit/probit 模型、有序 probit/logit 和計數數據的 Poisson 模型。 有關這些模型的一般概述，請參見 Long (1997) 和 Winkelmann 和 Boes (2006)。
這些方法的一個重要建模缺點是假設樣本中所有個體的偏好相同, 因此其對決策影響固定, 其估計出的係數為固定值且唯一，但是每個個體本質上是異質的( intrinsically individual heterogeneity )，因此同質性(homogeneity)假定是不現實的。例如，考慮以發表文章數量衡量的幼兒數量對科學家生產力的影響。由於有更多的孩子可能意味著更少的工作時間，我們預估平均而言，孩子的數量與出版物的水平之間存在負相關。這種負面關係有可能是全體性的，但它忽略了一個事實，即對於一些科學家來說，擁有更多孩子可能意味著擁有更有條理的生活方式，從而提高他們的生產力（例如，參見 Krapf 等人，2014 年）。因此，負係數可能隱藏顯著的個體異質性，從而導致對具有正係數的個體而言的誤導性推斷。同樣，我們可以假設所有科學家的影響都是負面的，但不利影響的幅度可能因樣本而異, 有的影響大有的影響小, 在這種情況下，存在個體異質性, 而非固定係數可以解釋的。最後，我們可能還會發現估計係數為零（不顯著）。在這種情況下，我們會得出結論，擁有年幼的孩子對生產力並不重要。儘管如此，這可能是由於樣本中科學家之間的異質性抵消了正面和負面影響(一半有正影響,一半有負影響, 導致整體估計係數為零)。在這些情況下，允許個體異質估計性的模型可能更合適。
2 方法論
2.1 具有隨機參數的模型
為了發展和激發隨機參數模型背後的想法，請考慮以下潛在過程
\[\begin{equation} y^* = {\bf x}_{it}^T \beta_{i} + \epsilon_{it} , \\ i=1,...,n; \, t=1,...,T \\ \beta_i ∼ g(\beta_i|\theta) \end{equation}\]
\(y^*\) 是個體 i 在時期 t 的潛在（未觀察到）過程因變數(dependent variables)，\(x_{it}\) 是一組解釋變數(covariates), \(\epsilon_{it}\) 隨機誤差項( random errors), \(g(β_i|\theta)\) 為 \(\theta\) 已知下的機率分配．注意一旦觀測到 \(y_{it}\) 的性質和隨機誤差項PDF 已知 ,則件概率密度函數(PDF) 的潛在過程 \(f(y^*|x_{it}, \beta_i)\) 可以得出． 例如 如果 \(y_{it}\) 是二進制且 $_{it} $ 服從正態分佈，則潛過程成為傳統的 probit 模型； 如果 \(y_{it}\) 是一個有序分類變量並且 \(\epsilon_{it}\) 是邏輯分佈(logistical distribution)的那麼傳統的有序 logit 模型就出現了。
二進制模型(binary model)、有序和泊松(Poisson)模型的 PDF 分別是
\[\begin{equation} f(y_i^* | x_{it},\beta_i)= F(x_{it}^T \beta_i)^{y_{it}} (1-F(x_{it}^T \beta_i)^{1-y_{it}} \end{equation}\]
有序模型 (ordered model) \[\begin{equation} f(y_i^* | x_{it},\beta_i)= \prod_{j=1}^{J} [F(k_j-x_{it}^T \beta_i) -F(k_{j-1}-x_{it}^T \beta_i))]^{y_{it}} \\ J =1, \dots, j-1, \, k_0=- \infty; \, k_J= + \infty \end{equation}\]
泊松(Poisson)模型
\[\begin{equation} \frac{1}{y_{it} !}\exp[-\exp( x_{it}^T \beta_i)] \exp( x_{it}^T \beta_i)^{y_{it}} \end{equation}\]
對於二元和有序模型，F(·) 表示誤差項的累積分佈函數 (CDF)，其中 Normal PDF \(F(\epsilon) = \Phi (\epsilon)\) 表示 probit， Logistic PDF: $F( ) = () $ 表示 logit:
\[ \Lambda (\epsilon) \equiv \frac{ \exp(\epsilon)}{1+\exp( \epsilon)} \]
在等式 (1) 給出的結構模型中，我們允許向量係數 \(\beta_i\) 對於群體中的每個個體都不同。 換句話說，對潛在因變量的邊際效應(marginal effect)是個體特定的。 然而，我們不知道這些參數如何因人而異。 我們所知道的是，它們根據總體 PDF \(g( \beta_i| \theta)\) 而變化，其中 \(\theta\) 表示分佈的動差(moments)，例如必須估計的均值和方差。 一旦給定了 \(g(\beta_i|\theta)\) 和 \(\epsilon\) 的分佈，就會出現一個完全參數化的模型。 在第 4 節中，我們將詳細介紹分佈。
為簡單起見，假設係數向量是獨立正態分佈的，因此對於 \(\beta_i\) 中的第 k 個元素， \(\beta_k ∼ N(\beta_k,\sigma_k^2)\)。 請注意，每個係數可以寫為 \(\beta_{ki} = \beta_k +\sigma_k \omega_i\)，其中 \(\omega_i ∼ N(0,1)\)，或以向量形式寫為 \(\beta_i = \beta+L \omega_i\)，其中 L 是包含標準偏差參數 \(\sigma_k\) 的對角矩陣。 標準偏差參數 \(\sigma_k\) 捕獲了有關每個單獨屬性的單獨異質性的所有信息。 如果 \(\sigma_k=0\)，則模型簡化為固定參數模型，但如果它確實顯著，則表明 \(x_{itk}\) 和 \(y_{it}\) 之間的關係是異質的，僅關注集中趨勢 \(\beta_k\) 將掩蓋有用信息 . 值得注意的是，隨機效應模型是一種特殊情況，其中只有常數是隨機的（參見第 5.3 節）。
2.2 擴展：觀察到的異質性和相關性
隨機參數模型的一個重要且直接的擴展是允許相關係數。 在這種情況下，L 不再是對角矩陣, L 是一個下三角矩陣（也稱為 Cholesky 矩陣），它產生隨機參數的協方差矩陣(variance-covariance matrix) \(LL^⊤ = \Sigma\)。
此外，觀察到的異質性也可以通過允許參數異質性受到部分可觀察到的變量(observed covariants; \(x_i\))影響。 這種類型的模型也稱為分層模型（hierarchical model Greene and Hensher 2010a）。 形式上，參數向量可以寫為
\[ \beta_i = \beta + \Pi s_i + L \omega_i\]
其中$ $ 是參數矩陣，\(s_i \in \bf X\) 是不隨時間變化的某些解釋變量向量，\(\omega ∼N(0,I)\)。 那麼，參數的均值為 \(E(\beta_i)=\beta+ \Pi s_i+LE(\omega)=\beta+ \Pi s_i\)，其協方差為
\[\text{VAR}(\beta_i) = E(L \omega ( \omega L)^⊤) = LE(\omega \omega ^⊤)L^⊤ = LIL^T=LL^T=\Sigma\]
作為上述公式的說明，假設我們對建模一些潛在過程感興趣，該過程由兩個變量解釋。 與每個觀察到的變量 \(\beta_{i1}\) 和 \(\beta_{i2}\) 相關的參數不是固定的，而是它們因個體而異並且是相關的。 此外，假設隨機參數的均值取決於一些觀察到的個體社會經濟變量影響, 如：\(S_i、B_i\) 和 \(C_i\)。 那麼，參數可以寫成：
\[\beta_{1,i} = \beta_1 + \Pi_{1,1} S_i + \Pi_{1,2}B_i + \Pi_{1,3}C_i + s_{11} \omega_{1,i}\]
\[\beta_{2,i} = \beta_2 + \Pi_{2,1} S_i + \Pi_{2,2}B_i + \Pi_{2,3}C_i + s_{21} \omega_{1,i} + s_{22} \omega_{2,i}\]
由於隨機參數的均值是觀測變量的函數，因此本規範放寬了關於均值個體間同質性的假設。 未觀察到的部分解釋了所有其他無法為觀察到的變量捕獲的個體特定因素。 請注意，隨機參數的方差-協方差矩陣為：
\[ \Sigma =LL^T= \begin{pmatrix} s_{11} & 0 \\ s_{21} & s_{22} \end{pmatrix} \begin{pmatrix} s_{11} & s_{21} \\ 0 & s_{22} \end{pmatrix} = \begin{pmatrix} s_{11}^2 & s_{11} s_{21} \\ s_{21}s_{11} & s_{21}^2+ s_{22}^2 \end{pmatrix} \]
條件均值向量為 \(E(\beta_i|s_i) = \beta+\Pi s_i\)，因 \(s_i\) 因個體而異, 因此不同個體有不同的條件均值 \(\beta_i\) 係數.
2.3 估計
在本節中，我將解釋 如何估計具有隨機參數的模型的模擬最大似然 (simulated maximum likelihood SML) 程序。 有關 SML 的更完整的詳細說明方法，請參見 Gourieroux 和 Monfort (1997)； Lee（1992）； Hajivassiliou 和 Ruud (1994) 或 Train (2009)。
設 \({\bf y_i} = {y_{i1}, y_{i2}, . . . , y_{iTi} }\) 是個體 i 做出的選擇序列(sequence of choices)。 假設個體在時間上是獨立的，在給定 \(\beta_i\) 下 \(y_i\)的聯合 PDF 可以寫成
\[ P(y_i|X_i,\beta_i )=\Pi_{t=1}^{T_i} f(y_{it}^*|x_{it} ,\beta_i) \]
其中 \(f(y_{it}^* |x_{it}, \beta_i)\) 見上節。 由於 \(\beta_i\) 是未觀察到的(unobserved)，因此我們需要對其聯合密度進行積分。 無條件 PDF 將是條件概率 在所有可能的 \(\beta\) 值上評估的加權平均值，這取決於 \(\beta_i\) 的分佈：
\[ P_i(\theta) =\int_{\beta_i} P(y_i|X_i, \beta_i)g(\beta_i)d \beta_i \]
式中的概率沒有封閉解(closed solution)，同時在大部分情況下 無法積分，因此很難執行最大似然 (ML) 估計。 但是，如果我們改為使用 $ P_i()$ 的良好近似 \(\tilde P_i(\theta )\) 來形成似然函數，則 ML 估計仍然是可能的。
但是，我們如何獲得 \(\tilde P_i(\theta )\)？ 通過蒙特卡洛積分(Monte Caro Integratal)可以得到一個很好的近似值。 此過程提供了確定性數值積分的替代方法。 在這裡，我們可以使用從分佈 \(g(\beta_i)\) 中隨機抽取來“模擬”積分。對於給定的參數 \(\theta\) 值， \(\beta_i\) 的值是從其分佈中抽取的。 使用 \(\beta_i\) 的採樣，計算公式 5 中的 \(P_i(\theta)\)。 多次重複採樣此過程，採樣的平均值是朗好得近似概率模擬。 形式上，個人 i 的模擬概率是
\[\tilde P_i(\theta) = \frac{1}{R} \Sigma_{i=1}^R \tilde P_i(y_i| {\bf X}, \beta_{ir} ) \]
其中 \(\tilde P_i(y_i| {\bf X}, \beta_{ir}\) 是在第 r 次 \(\beta_i\) 採樣時評估的個體 i 在時期 t 的模擬概率，R 是抽籤的總數。 給定對 i 的獨立性，SML 估計量是最大化的值 \(\theta\)：
\[ \hat \theta_{SML} \equiv \text{arg max}_{\theta \in \Theta} \Sigma_{i=1}^N \log \tilde P_i(\theta).\]
Lee (1992) 和 Hajivassiliou and Ruud (1994) 表明，在正則條件下(regular conditions)，SML 估計量是一致且漸近正態的(consistent and asymptotically normal)。 當採樣次數 R 的上升速度快於觀察數的平方根時(\(\sqrt N\))，估計量漸近等效於最大似然估計量 (asymptotically equivalent to the maximum likelihood estimator)。
3 Ex1: 婦女勞動力參與實證研究
Mroz（1987）從 1976 年收入動態 panel data 研究中提取的數據。 樣本包括 753 名 30 至 60 歲的已婚白人女性。資料變數介紹 lfp：如果妻子是有償勞動力，則為 1； 否則為 0，
k618：6 至 18 歲兒童的數量，
age：妻子的年齡，以歲為單位，
wc :1 如果妻子上過大學； 否則為 0，
hc1：如果丈夫上過大學； 否則為 0，
lwg：妻子估計對數工資率，
inc：不包括妻子工資的家庭收入，
linc：不包括妻子工資的家庭對數收入記錄，
3.1 Probit model
## 載入需要的套件：Formula## 載入需要的套件：maxLik## 載入需要的套件：miscTools##
## Please cite the 'maxLik' package as:
## Henningsen, Arne and Toomet, Ott (2011). maxLik: A package for maximum likelihood estimation in R. Computational Statistics 26(3), 443-458. DOI 10.1007/s00180-010-0217-1.
##
## If you have questions, suggestions, or comments regarding the 'maxLik' package, please use a forum or 'tracker' at maxLik's R-Forge site:
## https://r-forge.r-project.org/projects/maxlik/前面六位受訪婦女的資料 ．
| lfp | k5 | k618 | age | wc | hc | lwg | inc | linc |
|---|---|---|---|---|---|---|---|---|
| 0 | 0 | 3 | 39 | 0 | 0 | 0.8532125 | 28.363 | 3.345086 |
| 0 | 0 | 0 | 60 | 0 | 0 | 1.2249736 | 24.984 | 3.218236 |
| 0 | 0 | 0 | 43 | 0 | 0 | 0.8881401 | 9.952 | 2.297774 |
| 0 | 2 | 3 | 31 | 0 | 0 | 1.1580402 | 10.000 | 2.302585 |
| 0 | 0 | 2 | 40 | 1 | 1 | 1.0828638 | 28.200 | 3.339322 |
| 0 | 0 | 2 | 36 | 0 | 0 | 0.8895015 | 5.330 | 1.673351 |
### Probit model
work.probit <- Rchoice(lfp ~ k5 + k618 + age + wc + hc + lwg + inc,
data = Workmroz, family = binomial('probit'))
work.probit.sum<-summary(work.probit)
knitr::kable(round(work.probit.sum$CoefTable,3),format = 'simple',caption=' Probit Model:Labor Force Participation')| Estimate | Std. Error | z-value | Pr(>|z|) | |
|---|---|---|---|---|
| constant | 1.918 | 0.381 | 5.040 | 0.000 |
| k5 | -0.875 | 0.114 | -7.703 | 0.000 |
| k618 | -0.039 | 0.040 | -0.953 | 0.340 |
| age | -0.038 | 0.008 | -4.971 | 0.000 |
| wc | 0.488 | 0.135 | 3.604 | 0.000 |
| hc | 0.057 | 0.124 | 0.461 | 0.645 |
| lwg | 0.366 | 0.088 | 4.165 | 0.000 |
| inc | -0.021 | 0.005 | -4.297 | 0.000 |
由估計的結果可知, 使婦女勞動參與率下降的因素有:婦女本身年齡越高,小孩數越多小孩, 有五歲小孩以及婦女只有高中學歷以下(但不顯著), 以及不包括妻子工資的家庭收入越高,都會使得 婦女勞動參與率下降． 造成婦女勞動參與率提高的因素有婦女工資率越高, 以及婦女有大學畢業學歷及丈夫有大學畢業學歷(但不顯著).
3.2 Logit model
work.logit <- Rchoice(lfp ~ k5 + k618 + age + wc + hc + lwg + inc,
data = Workmroz, family = binomial('logit'))
work.logit.sum<-summary(work.logit)
knitr::kable(round(work.logit.sum$CoefTable,3),format = 'simple',caption=' Logit Model:Labor Force Participation')| Estimate | Std. Error | z-value | Pr(>|z|) | |
|---|---|---|---|---|
| constant | 3.182 | 0.644 | 4.938 | 0.000 |
| k5 | -1.463 | 0.197 | -7.426 | 0.000 |
| k618 | -0.065 | 0.068 | -0.950 | 0.342 |
| age | -0.063 | 0.013 | -4.918 | 0.000 |
| wc | 0.807 | 0.230 | 3.510 | 0.000 |
| hc | 0.112 | 0.206 | 0.542 | 0.588 |
| lwg | 0.605 | 0.151 | 4.009 | 0.000 |
| inc | -0.034 | 0.008 | -4.196 | 0.000 |
兩個模型估計結果差不多. probit model 估計出的係數大約是logit model 0.6倍;也就是 logit model 估計出的係數大約是 probit model 的 1.7倍.
4 Ex2: 博士論文發表量的性別差異: Sex Differences in Science: Poisson model
博士生畢業以後論文發表量的性別差異會受到哪些因素的影響?資料來源相關期刊相關論文
Long, J. S. (1990). The origins of sex differences in science. Social Forces, 68(4), 1297-1316.
Long, J. S. (1997). Regression models for categorical and limited dependent variables (Vol. 7). Sage.
Long, J. S., & Freese, J. (2006). Regression models for categorical and limited dependent variables using Stata. Stata Press, College Station, TX.
Xie, Y. and Shauman K. A. (1998). Sex Differences in Research Productivity: New Evidence about an Old Puzzle.
American Sociological Review, Vol. 63, No. 6 (Dec., 1998), pp. 847-870. Stable URL: http://www.jstor.org/stable/2657505
- Long, J. S. (1990) 結論如下:
科學社會學 (sociology of science) 清楚地確定了科學生產力和地位的性別差異 (sex differences)。 Long(1990) 研究了導致女性科學家在完成博士培訓後生產力降低的過程。
首先發現與導師的合作被發現是影響女性科學家生產力的最重要因素。
其次發現對於女性來說，生育年幼的孩子會大大減少合作的機會。因此，幼兒的存在對女科學家在研究生學習期間的生產力產生了不利的間接影響。男性不存在這種效果。
除了合作過程的差異外，在影響生產力的資源水平和資源轉化為生產力的機制上，也發現了許多對女性不利和對男性有利的細微差異。小劣勢的集中進一步解釋了職業生涯開始時生產力的性別差異。由於已經發現早期的優勢和劣勢是累積的，因此本文為了解在職業生涯中出現的科學生產力和職位的性別差異提供了必不可少的第一步。
- Xie and Shauman(1998) 結論如下:
許多研究發現，女性科學家發表論文的速度低於男性科學家。到目前為止，對這種一致模式的解釋還沒有出現，研究生產力的性別差異仍然是一個謎。
該報告基於對 1969 年、1973 年、1988 年和 1993 年的四次大型、具有全國代表性的高等教育教師橫斷面調查的數據進行了系統和詳細的分析。
Xie and Shauman(1998)的研究產生了兩個主要發現。首先，研究生產力的性別差異在研究期間有所下降,也就是說女性研究生產力有逐年提高的趨勢，女性研究生產力與男性研究生產力的比例從 1960 年代後期的約 60% 增加到 1980 年代後期和 1990 年代初期的 75% 到 80%。
其次，觀察到的研究生產力的性別差異大部分可歸因於個人特徵、結構立場和婚姻狀況的性別差異。
這些結果表明，研究生產力的性別差異源於結構位置的性別差異，因此對女性在科學中地位的長期提高作出反應。
我們重做此篇論文結果如下:
data("Articles")
knitr::kable(head(Articles),'simple',caption='Doctoral Publications Data')| art | fem | mar | kid5 | phd | ment |
|---|---|---|---|---|---|
| 0 | 0 | 1 | 0 | 2.52 | 7 |
| 0 | 1 | 0 | 0 | 2.05 | 6 |
| 0 | 1 | 0 | 0 | 3.75 | 6 |
| 0 | 0 | 1 | 1 | 1.18 | 3 |
| 0 | 1 | 0 | 0 | 3.75 | 26 |
| 0 | 1 | 1 | 2 | 3.59 | 2 |
Long（1990, 97）的研究數據分析了科學家的出版物水平。 具有以下 6 個變量的 915 個觀察值的數據：
art:博士最後三年的文章，
fem: 1 如果是女科學家； 否則為 0，
mar: 1 如果已婚； 否則為 0，
kid5: 5 歲或以下的兒童人數，
phd: 就讀博士學位學校系所聲望，
ment:近三年導師的文章數量.
4.1 Poisson model with fixed coefficients
art.poisson.fixed <- Rchoice(art ~ fem + mar + kid5 + phd + ment, data = Articles, family = poisson)
art.poisson.fix.sum=summary(art.poisson.fixed)
knitr::kable(round(art.poisson.fix.sum$CoefTable,3),format = 'simple',caption=' Poisson model with fixed coefficients:Doctoral Publications')| Estimate | Std. Error | z-value | Pr(>|z|) | |
|---|---|---|---|---|
| constant | 0.305 | 0.103 | 2.958 | 0.003 |
| fem | -0.225 | 0.055 | -4.112 | 0.000 |
| mar | 0.155 | 0.061 | 2.529 | 0.011 |
| kid5 | -0.185 | 0.040 | -4.607 | 0.000 |
| phd | 0.013 | 0.026 | 0.486 | 0.627 |
| ment | 0.026 | 0.002 | 12.733 | 0.000 |
估計結果顯示: 對數似然函數 (Log Likelihood= -1651) 是使用 Newton-Raphson 算法,在第 7 次迭代到達收斂。 在解釋變數方面，我們可以說，作為一名女性科學家，在其他條件不變的情況下(centurus-paribus)，預期文章數量會減少 0.8 倍（= exp(-0.225)）。 或作為一名女科學家，在其他條件不變的情況下，預期文章數量會減少 20% (= 100 [exp(−0.225) − 1])。 博士系所的聲望對生產力並不重要。
在其他條件不變的情況下， 科學家每增加一位 5 歲以下的兒童,則預期其論文發表數量會減少約 17%((exp(-0.185)-1)*100)。
在其他條件不變的情況下，指導教授每增加一篇論文發表,此科學家預期論文數量會增加約 2.64%((exp(0.026)-1)*100)。 因此指導教授的研究能力(發表論文數量) 深深地影響女性科學家的發表文章的數量以及研究能力.
在其他條件不變的情況下，已婚科學家發表論文的篇數比未婚科學家的篇數多出 約16.8 %(exp(0.155)-1)*100= 16.8%.
4.2 Poisson model with random coefficients
art.poisson.ran <- Rchoice(art ~ fem + mar + kid5 + phd + ment,
data = Articles, family = poisson,
ranp = c(kid5 = "n", phd = "n", ment = "n", mar = "n"))
art.poisson.ran.sum=summary(art.poisson.ran)
knitr::kable(round(art.poisson.ran.sum$CoefTable,3),format = 'simple',
caption='Poisson model with random coefficients:Doctoral Publications')| Estimate | Std. Error | z-value | Pr(>|z|) | |
|---|---|---|---|---|
| constant | 0.193 | 0.132 | 1.461 | 0.144 |
| fem | -0.216 | 0.071 | -3.020 | 0.003 |
| mean.mar | 0.071 | 0.085 | 0.840 | 0.401 |
| mean.kid5 | -0.164 | 0.057 | -2.883 | 0.004 |
| mean.phd | -0.007 | 0.038 | -0.188 | 0.851 |
| mean.ment | 0.031 | 0.004 | 7.980 | 0.000 |
| sd.mar | 0.406 | 0.084 | 4.862 | 0.000 |
| sd.kid5 | 0.129 | 0.088 | 1.467 | 0.142 |
| sd.phd | 0.140 | 0.019 | 7.305 | 0.000 |
| sd.ment | 0.016 | 0.004 | 3.940 | 0.000 |
但是在隨機參數模型係數的估計下,我們發現以上的結論不成立. 因為受到異質性的影響, 平均數的估計不顯著, 但是係數標準車很顯著, 代表有相當的異質性,所以有一半的人是有正影響,有一半的人負影響,使得平均效果為零.
下表為兩種不同的方法估計出來的係數差異. 隨機參數模型係數估計值明顯小於 固定參數模型係數估計值.
knitr::kable(round(data.frame(rbind(fixed=art.poisson.fix.sum$coefficients[1:6],
random=art.poisson.ran.sum$coefficients[1:6])),3),
'pipe',caption='poisson model with fixed and random paraneters')| constant | fem | mar | kid5 | phd | ment | |
|---|---|---|---|---|---|---|
| fixed | 0.305 | -0.225 | 0.155 | -0.185 | 0.013 | 0.026 |
| random | 0.193 | -0.216 | 0.071 | -0.164 | -0.007 | 0.031 |
5 Health Demand: Ordered probit model
Riphahn, R. T., Wambach, A., & Million, A. (2003). Incentive effects in the demand for health care: a bivariate panel count data estimation. Journal of applied econometrics, 18(4), 387-405.
此文在三個方面對有關醫療保健需求的文獻做出了貢獻。
1.: 本文是第一篇在計數資料(count data)中應用雙變量隨機效應 probit 模型(bivariate random effects probit model)，以允許有效估計這種帶有面板數據(panel data)的模型。
2.:本文提供了一種創新的逆選擇(adverse selection)測試，並證實了 高風險個人更有可能獲得補充附加保險(supplemental add-on insurance)。
3.: 本文估計結果表明，根據道德危機理論 (moral hazard):道德危機，即被保險人的醫療需求會超過 未投保人的醫療需求，本文觀察到個體經營者和有小孩的母親門診次數比一般病患頻率要低得多.
資料簡介: 德國醫療保健數據，不平衡面板數據(unblanced panel data)。具有以下 27 個變量的 27326 個觀測值的數據框,研究時間分別為: 1984 1985 1986 1988 1987 1991 1994.
hsat：健康滿意度，0（低），…，10（高）
newhsat:健康滿意度(hsat)分成 4 級(ordered-choice), (0-2) = 0, (3-5)=1, (6-8)=2, (9)=3, (10)=4. 此為因變數(dependent variable),4 級 ordered probe model .
hhkids：家中有16 歲以下兒童 = 1； 否則 = 0
hhinc：以德國馬克計的家庭名目月淨收入(nominal income).
## Ordered probit model
data("Health")
knitr::kable(head(Health),format='simple',caption='Health')| id | female | year | age | hsat | handdum | handper | hhinc | hhkids | educ | married | haupts | reals | fachhs | abitur | univ | working | bluec | whitec | self | beamt | docvis | hospvis | public | addon | hsat2 | newhsat |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | 0 | 1984 | 54 | 8 | 0 | 0 | 3050.000 | 0 | 15 | 1 | 0 | 0 | 1 | 0 | 0 | 1 | 0 | 1 | 0 | 0 | 1 | 0 | 1 | 0 | 8 | 2 |
| 1 | 0 | 1985 | 55 | 8 | 0 | 0 | 4510.050 | 0 | 15 | 1 | 0 | 0 | 1 | 0 | 0 | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 1 | 0 | 8 | 2 |
| 1 | 0 | 1986 | 56 | 7 | 0 | 0 | 3500.000 | 0 | 15 | 1 | 0 | 0 | 1 | 0 | 0 | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 1 | 0 | 7 | 2 |
| 2 | 1 | 1984 | 44 | 7 | 0 | 0 | 3050.000 | 0 | 9 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 7 | 2 |
| 2 | 1 | 1985 | 45 | 8 | 0 | 0 | 3182.779 | 0 | 9 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 1 | 0 | 8 | 2 |
| 2 | 1 | 1986 | 46 | 7 | 0 | 0 | 3500.000 | 0 | 9 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 2 | 0 | 1 | 0 | 7 | 2 |
估計 1988 年的結果
health.probit.1988<- Rchoice(newhsat ~ age + educ + hhinc + married + hhkids,
data = Health, family = ordinal('probit'),subset = year == 1988)
summary(health.probit.1988)##
## Model: ordinal
## Model estimated on: 日 3 27 09時49分33秒 2022
##
## Call:
## Rchoice(formula = newhsat ~ age + educ + hhinc + married + hhkids,
## data = Health, subset = year == 1988, family = ordinal("probit"),
## method = "bfgs")
##
##
## Frequencies of categories:
## y
## 0 1 2 3 4
## 0.05130 0.24827 0.49654 0.11153 0.09235
## The estimation took: 0h:0m:1s
##
## Coefficients:
## Estimate Std. Error z-value Pr(>|z|)
## kappa.1 1.148e+00 3.160e-02 36.335 < 2e-16 ***
## kappa.2 2.548e+00 3.720e-02 68.496 < 2e-16 ***
## kappa.3 3.056e+00 4.073e-02 75.047 < 2e-16 ***
## constant 1.979e+00 1.196e-01 16.552 < 2e-16 ***
## age -1.808e-02 1.629e-03 -11.097 < 2e-16 ***
## educ 3.556e-02 7.141e-03 4.980 6.36e-07 ***
## hhinc 2.587e-05 1.039e-05 2.489 0.0128 *
## married -3.099e-02 4.203e-02 -0.737 0.4609
## hhkids 6.064e-02 3.824e-02 1.586 0.1128
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Optimization of log-likelihood by BFGS maximization
## Log Likelihood: -5753
## Number of observations: 4483
## Number of iterations: 96
## Exit of MLE: successful convergence1984,1994 健康需求比較
health.probit.1984<- Rchoice(newhsat ~ age + educ + hhinc + married + hhkids,
data = Health, family = ordinal('probit'),
subset = year == 1984)
(sum84<-summary(health.probit.1984))##
## Model: ordinal
## Model estimated on: 日 3 27 09時49分34秒 2022
##
## Call:
## Rchoice(formula = newhsat ~ age + educ + hhinc + married + hhkids,
## data = Health, subset = year == 1984, family = ordinal("probit"),
## method = "bfgs")
##
##
## Frequencies of categories:
## y
## 0 1 2 3 4
## 0.06557 0.24522 0.39520 0.10145 0.19257
## The estimation took: 0h:0m:1s
##
## Coefficients:
## Estimate Std. Error z-value Pr(>|z|)
## kappa.1 1.073e+00 3.147e-02 34.088 < 2e-16 ***
## kappa.2 2.169e+00 3.659e-02 59.275 < 2e-16 ***
## kappa.3 2.509e+00 3.840e-02 65.343 < 2e-16 ***
## constant 1.926e+00 1.348e-01 14.285 < 2e-16 ***
## age -2.300e-02 1.758e-03 -13.083 < 2e-16 ***
## educ 4.134e-02 8.118e-03 5.093 3.53e-07 ***
## hhinc 5.239e-05 1.219e-05 4.297 1.73e-05 ***
## married 3.118e-02 4.603e-02 0.677 0.4981
## hhkids 8.800e-02 4.003e-02 2.199 0.0279 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Optimization of log-likelihood by BFGS maximization
## Log Likelihood: -5404
## Number of observations: 3874
## Number of iterations: 95
## Exit of MLE: successful convergencehealth.probit.1994<- Rchoice(newhsat ~ age + educ + hhinc + married + hhkids,
data = Health, family = ordinal('probit'),
subset = year == 1994)
(sum94<-summary(health.probit.1994))##
## Model: ordinal
## Model estimated on: 日 3 27 09時49分35秒 2022
##
## Call:
## Rchoice(formula = newhsat ~ age + educ + hhinc + married + hhkids,
## data = Health, subset = year == 1994, family = ordinal("probit"),
## method = "bfgs")
##
##
## Frequencies of categories:
## y
## 0 1 2 3 4
## 0.04679 0.25022 0.51436 0.11312 0.07551
## The estimation took: 0h:0m:1s
##
## Coefficients:
## Estimate Std. Error z-value Pr(>|z|)
## kappa.1 1.185e+00 3.756e-02 31.561 < 2e-16 ***
## kappa.2 2.654e+00 4.420e-02 60.046 < 2e-16 ***
## kappa.3 3.222e+00 4.927e-02 65.399 < 2e-16 ***
## constant 2.063e+00 1.330e-01 15.508 < 2e-16 ***
## age -1.968e-02 1.894e-03 -10.390 < 2e-16 ***
## educ 2.494e-02 8.145e-03 3.062 0.0022 **
## hhinc 4.571e-05 9.142e-06 5.000 5.73e-07 ***
## married 2.596e-02 4.888e-02 0.531 0.5953
## hhkids 1.293e-02 4.455e-02 0.290 0.7716
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Optimization of log-likelihood by BFGS maximization
## Log Likelihood: -4200
## Number of observations: 3377
## Number of iterations: 136
## Exit of MLE: successful convergenceknitr::kable(round(rbind('1984'=sum84$coefficients,'1994'=sum94$coefficients ),3),'simple',
caption='Health Demand: 1984 and 1994 years')| kappa.1 | kappa.2 | kappa.3 | constant | age | educ | hhinc | married | hhkids | |
|---|---|---|---|---|---|---|---|---|---|
| 1984 | 1.073 | 2.169 | 2.509 | 1.926 | -0.023 | 0.041 | 0 | 0.031 | 0.088 |
| 1994 | 1.185 | 2.654 | 3.222 | 2.063 | -0.020 | 0.025 | 0 | 0.026 | 0.013 |
6 Panel Data: Poisson Model with Random Parameters
6.1 具有橫截面數據的隨機參數模型 (Random parameters models with cross-sectional data)
本節，將展示如何使用橫截面數據(cross-sectional data)來估計這些類型的計數數據模型。 繼續之前的泊松模型，現在我將假設 Kid5、phd 和 ment 的影響不是固定的，而是在整個人群中是異質的。 具體來說，我將假設這些變量的係數是獨立的正態分佈，也就是說，它們之間不存在相關性。 特別是，
\[ \beta_{kid5,ir} = \beta_{kid5} + \sigma_{kid5} \omega_{kid5,ir} \]
\[ \beta_{phd,ir} = \beta_{phd} + \sigma_{phd} \omega_{phd,ir} \]
\[ \beta_{ment,ir} = \beta_{ment} + \sigma_{ment} \omega_{ment,ir} \]
其中 \(\omega_{k,ir} ∼ N(0,1)\)。
6.2 Poisson Model with Random Parameters
art.poisson.ran <- Rchoice(art ~ fem + mar + kid5 + phd + ment,
data = Articles, family = poisson,
ranp = c(kid5 = "n", phd = "n", ment = "n"))
summary(art.poisson.ran) ##
## Model: poisson
## Model estimated on: 日 3 27 09時49分52秒 2022
##
## Call:
## Rchoice(formula = art ~ fem + mar + kid5 + phd + ment, data = Articles,
## family = poisson, ranp = c(kid5 = "n", phd = "n", ment = "n"),
## method = "bfgs", iterlim = 2000)
##
## The estimation took: 0h:0m:17s
##
## Coefficients:
## Estimate Std. Error z-value Pr(>|z|)
## constant 0.225583 0.132500 1.703 0.08866 .
## fem -0.218498 0.070558 -3.097 0.00196 **
## mar 0.156431 0.079121 1.977 0.04803 *
## mean.kid5 -0.197775 0.063472 -3.116 0.00183 **
## mean.phd -0.029942 0.037217 -0.805 0.42109
## mean.ment 0.031110 0.003814 8.158 4.44e-16 ***
## sd.kid5 0.285310 0.089104 3.202 0.00136 **
## sd.phd 0.165405 0.016585 9.973 < 2e-16 ***
## sd.ment 0.015876 0.003535 4.491 7.11e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Optimization of log-likelihood by BFGS maximization
## Log Likelihood: -1574
## Number of observations: 915
## Number of iterations: 81
## Exit of MLE: successful convergence
## Simulation based on 40 Halton draws估計結果顯示: 對數似然函數 (Log Likelihood= -1574, 高於 固定係數模型的 -1651) 是使用 Newton-Raphson 算法在 81 次迭代後收斂。 結果還表明，係數的標準差非常顯著，表明參數確實在總體中有所不同。 由於參數是正態分佈的，我們也可以說：
kid5 平均數係數估計為 -0.19 標準差為 0.28 二者估計值都極為顯著. 因為 mean.kid5/sd.kod5=-0.7, 因此 \(N(x \le -0.7) \approx 0.24\), 也就是約有 24% 的人對 Kid5 的係數為正。 換句話說，對於大約 76% 的博士生來說，孩子不到 6 歲會降低他們的工作效率。
因此 在其他條件不變的情況下， 科學家每增加一位 5 歲以下的兒童,則預期其論文發表數量會減少約 17.3%((exp(-0.19)-1)*100)。
phd 平均數係數估計為 -0.03 (不顯著) 標準差為 0.165(極為顯著).因此博士畢業學校聲望有一半是有正影響,有一半是有負影響造成平均來講是沒有影響.
ment 平均數係數估計為 0.03 (極為顯著) 標準差為 0.016(極為顯著).因此 在其他條件不變的情況下，科學家論文執導教授發表論文的篇數越多則其指導科學家的論文篇數會多出約3 %(exp(0.03)-1)*100.
由以上隨機係數模型種種結果發現,這些結果都是使用固定參數泊松回歸模型不可能進行觀察得到的。
假設現在我們要測試是否 \(H_0 = \sigma_{kid5} = \sigma_{phd} = \sigma_{ment}} = 0\) 。這可以做到
library("lmtest")## 載入需要的套件：zoo##
## 載入套件：'zoo'## 下列物件被遮斷自 'package:base':
##
## as.Date, as.Date.numericknitr::kable(lrtest(art.poisson.fixed, art.poisson.ran),'pipe')| #Df | LogLik | Df | Chisq | Pr(>Chisq) |
|---|---|---|---|---|
| 6 | -1651.056 | NA | NA | NA |
| 9 | -1574.166 | 3 | 153.7807 | 0 |
由概似比檢定(Likelihood ratio test)可以知道隨機係數模型優於固定係數模型.
6.3 參數不同分配設定
我們假定指定參數不同分配設定, 以找出較佳模型, 假設 kid5 ~ uniform distribution, phd ~ t distribution, ment ~ censored normal distribution.
art.poisson.ran2 <- update(art.poisson.ran, ranp = c(kid5 = "u", phd = "t" , ment = "cn"))
summary(art.poisson.ran2)##
## Model: poisson
## Model estimated on: 日 3 27 09時50分18秒 2022
##
## Call:
## Rchoice(formula = art ~ fem + mar + kid5 + phd + ment, data = Articles,
## family = poisson, ranp = c(kid5 = "u", phd = "t", ment = "cn"),
## method = "bfgs", iterlim = 2000)
##
## The estimation took: 0h:0m:25s
##
## Coefficients:
## Estimate Std. Error z-value Pr(>|z|)
## constant 0.253728 0.131734 1.926 0.054096 .
## fem -0.216431 0.068964 -3.138 0.001699 **
## mar 0.148878 0.076371 1.949 0.051248 .
## mean.kid5 -0.221396 0.064356 -3.440 0.000581 ***
## mean.phd -0.104283 0.042852 -2.434 0.014952 *
## mean.ment 0.028434 0.004889 5.816 6.02e-09 ***
## sd.kid5 0.507644 0.123447 4.112 3.92e-05 ***
## sd.phd 0.226169 0.025087 9.016 < 2e-16 ***
## sd.ment 0.022482 0.004485 5.013 5.36e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Optimization of log-likelihood by BFGS maximization
## Log Likelihood: -1576
## Number of observations: 915
## Number of iterations: 72
## Exit of MLE: successful convergence
## Simulation based on 40 Halton draws將兩者不同設定模型就互相比較
library("memisc")
(Ttable <- mtable("Model 1"= art.poisson.ran, "Model 2" = art.poisson.ran2,
summary.stats = c("N", "Log-likelihood", "BIC", "AIC")))##
## Calls:
## Model 1: Rchoice(formula = art ~ fem + mar + kid5 + phd + ment, data = Articles,
## family = poisson, ranp = c(kid5 = "n", phd = "n", ment = "n"),
## method = "bfgs", iterlim = 2000)
## Model 2: Rchoice(formula = art ~ fem + mar + kid5 + phd + ment, data = Articles,
## family = poisson, ranp = c(kid5 = "u", phd = "t", ment = "cn"),
## method = "bfgs", iterlim = 2000)
##
## ==============================================
## Model 1 Model 2
## ----------------------------------------------
## constant 0.226 0.254
## (0.132) (0.132)
## fem -0.218** -0.216**
## (0.071) (0.069)
## mar 0.156* 0.149
## (0.079) (0.076)
## mean.kid5 -0.198** -0.221***
## (0.063) (0.064)
## mean.phd -0.030 -0.104*
## (0.037) (0.043)
## mean.ment 0.031*** 0.028***
## (0.004) (0.005)
## sd.kid5 0.285** 0.508***
## (0.089) (0.123)
## sd.phd 0.165*** 0.226***
## (0.017) (0.025)
## sd.ment 0.016*** 0.022***
## (0.004) (0.004)
## ----------------------------------------------
## N 915 915
## Log-likelihood -1574.166 -1575.816
## BIC 3209.702 3213.003
## AIC 3166.332 3169.632
## ==============================================
## Significance: *** = p < 0.001;
## ** = p < 0.01; * = p < 0.05 knitr::kable(lrtest(art.poisson.ran , art.poisson.ran2 ),'pipe')| #Df | LogLik | Df | Chisq | Pr(>Chisq) |
|---|---|---|---|---|
| 9 | -1574.166 | NA | NA | NA |
| 9 | -1575.816 | 0 | 3.300534 | 0 |
由概似比檢定, AIC, BIC 皆可以發現常態分配假設(model 1)優於其他分配設定(model 2). 因此 係數 kid5, phd , ment 皆以 normal distribution設定較佳.
7 面板數據隨機效應的隨機參數模型(Panel data with random effects and random parameters models)
假設第 i 個觀測值 \(y_i\)其聯合密度 \(f(y_i^*|X_i, \alpha_i, \beta)\)，隨機效應的密度為 \(\alpha_i \sim g(\alpha_i|\theta)\)，其中 $g(_i|) $ 不依賴於可觀測值。 因此，第 i 個觀測值聯合密度為
\[f(y_i^*|X_i, \beta) =\int_{\alpha_i} f(y_i^*|X_i, \alpha_i, \beta) g(\alpha_i|\theta) d\alpha_i \]
該模型可以使用高斯-赫米特正交 (Gauss- Hermite quadrature )來近似概率的積分。
7.1 An binary probit model with random effects and random parameters
Vella, F. and M. Verbeek (1998) “Whose wages do unions raise ? A dynamic model of unionism and wage”, Journal of Applied Econometrics, 13, 163–183.
工會和工資率:從 1980 年到 1987 年對 545 個人的年度觀資料, 觀察次數：4360,國家 : 美國.
expe:經驗，計算為年齡 - 6 - 受教育年限.
rural:個人住在鄉村嗎？
school:受教育年限.
union:工資是通過集體談判確定的嗎?
wage: 工資以美元計的小時工資.
假設 \(\alpha_i \sim N( 0, \sigma^2_\alpha)\);
\(\beta_{wage}\) ~ log-normal distribution.
\[\beta_{\text{wage},ir} = exp(\beta_{wage} +\sigma_k \omega_{k,ir}) \] \(\omega_{k,ir} ∼ N(0,1)\)
library(pglm)## 載入需要的套件：plm##
## 載入套件：'pglm'## 下列物件被遮斷自 'package:Rchoice':
##
## ordinaldata('UnionWage', package = 'pglm')
UnionWage$union1 <- as.numeric(UnionWage$union)-1
union.ran <- Rchoice(union1 ~ exper + school + rural + wage,
data = UnionWage[1:2000, ],
family = binomial('probit'),
ranp = c(constant = "n", wage = "ln"),
R = 10,
panel = TRUE,
index = "id",
print.init = TRUE)##
## Starting Values:
## exper school ruralyes mean.constant mean.wage
## -0.03966858 -0.14012457 0.25144225 0.16690100 -0.55766923
## sd.constant sd.wage
## 0.10000000 0.10000000summary(union.ran)##
## Model: binomial
## Model estimated on: 日 3 27 09時52分01秒 2022
##
## Call:
## Rchoice(formula = union1 ~ exper + school + rural + wage, data = UnionWage[1:2000,
## ], family = binomial("probit"), ranp = c(constant = "n",
## wage = "ln"), R = 10, panel = TRUE, index = "id", print.init = TRUE,
## method = "bfgs", iterlim = 2000)
##
##
## Frequencies of categories:
## y
## 0 1
## 0.7605 0.2395
## The estimation took: 0h:0m:55s
##
## Coefficients:
## Estimate Std. Error z-value Pr(>|z|)
## exper -0.08845 0.01921 -4.604 4.14e-06 ***
## school -0.24919 0.06624 -3.762 0.000168 ***
## ruralyes 0.21317 0.20098 1.061 0.288829
## mean.constant 1.68815 0.79492 2.124 0.033697 *
## mean.wage -2.17001 0.40645 -5.339 9.35e-08 ***
## sd.constant 1.47650 0.14746 10.013 < 2e-16 ***
## sd.wage 2.30410 0.31141 7.399 1.37e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Optimization of log-likelihood by BFGS maximization
## Log Likelihood: -744.4
## Number of observations: 2000
## Number of iterations: 81
## Exit of MLE: successful convergence
## Simulation based on 10 Halton drawssd.constant and mean.constant 顯著不等於零, 因此有隨機效應效果.
pnorm(coef(union.ran)["mean.wage"]/coef(union.ran)["sd.wage"])## mean.wage
## 0.17314717 % 負效應, 83 % 有正效應.
8 具有觀察到的異質性的隨機參數模型
在本節中，我將說明如何估計具有觀察到的異質性的泊松隨機參數模型(random parameter model with observed heterogeneity) 。 在下面的例子中，我假設在 phd 和 ment 的係數中不僅存在未觀察到的異質性，而且在係數的平均值中也存在觀察到的異質性。 具體來說，假設係數因人而異：
\[\beta_{kid5,ir} = \beta_{kid5} + \sigma_{kid5}\omega_{kid5,ir}\]
\[\beta_{phd,ir} = \beta_{phd} + \pi_{phd,fem} fem+ \sigma_{phd}\omega_{phd,ir}\]
\[\beta_{ment,ir} = \beta_{ment} + \pi_{ment,fem} fem +\pi_{ment,phd} phd + \sigma_{ment}\omega_{ment,ir}\]
上面的公式意味著，ment 的係數（或對潛在生產力的邊際效應）因 fem 性別和 博士畢業學校的聲望(phd)而異; phd 的係數因 fem 性別而異 。
## Hierarchical Poisson Model
art.poissonH.ran <- Rchoice(art ~ fem + mar + kid5 + phd + ment | fem + phd,
data = Articles,
ranp = c(kid5 = "n", phd = "n", ment = "n"),
mvar = list(phd = c("fem"), ment = c("fem", "phd")),
family = poisson,
R = 10)
summary(art.poissonH.ran)##
## Model: poisson
## Model estimated on: 日 3 27 09時52分43秒 2022
##
## Call:
## Rchoice(formula = art ~ fem + mar + kid5 + phd + ment | fem +
## phd, data = Articles, family = poisson, ranp = c(kid5 = "n",
## phd = "n", ment = "n"), R = 10, mvar = list(phd = c("fem"),
## ment = c("fem", "phd")), method = "bfgs", iterlim = 2000)
##
## The estimation took: 0h:0m:42s
##
## Coefficients:
## Estimate Std. Error z-value Pr(>|z|)
## constant 0.222646 0.181920 1.224 0.22100
## fem -0.572068 0.217337 -2.632 0.00848 **
## mar 0.176966 0.075165 2.354 0.01855 *
## mean.kid5 -0.259958 0.061981 -4.194 2.74e-05 ***
## mean.phd -0.019339 0.056808 -0.340 0.73354
## mean.ment 0.048094 0.009960 4.829 1.37e-06 ***
## phd.fem 0.133917 0.068995 1.941 0.05226 .
## ment.fem -0.006846 0.007134 -0.960 0.33721
## ment.phd -0.004857 0.002902 -1.674 0.09415 .
## sd.kid5 0.431139 0.085080 5.067 4.03e-07 ***
## sd.phd 0.133125 0.015333 8.682 < 2e-16 ***
## sd.ment 0.014867 0.003050 4.874 1.09e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Optimization of log-likelihood by BFGS maximization
## Log Likelihood: -1577
## Number of observations: 915
## Number of iterations: 99
## Exit of MLE: successful convergence
## Simulation based on 10 Halton draws估計的參數表明，性別僅對博士畢業學校的聲望(phd)平均係數有重要顯著正影響, 性別對ment(導師文章數)的平均係數不重要。 此外，博士畢業學校的聲望降低了ment(導師文章數)的影響。 我們可以使用likelihood-test 測試交互變量是否聯合顯著：
lrtest(art.poissonc.ran, art.poissonH.ran)## Likelihood ratio test
##
## Model 1: art ~ fem + mar + kid5 + phd + ment
## Model 2: art ~ fem + mar + kid5 + phd + ment | fem + phd
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 12 -1570.8
## 2 12 -1576.9 0 12.242 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 lrtest(art.poisson.ran, art.poissonH.ran)## Likelihood ratio test
##
## Model 1: art ~ fem + mar + kid5 + phd + ment
## Model 2: art ~ fem + mar + kid5 + phd + ment | fem + phd
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 9 -1574.2
## 2 12 -1576.9 3 5.4385 0.1424 lrtest(art.poisson.ran, art.poissonc.ran)## Likelihood ratio test
##
## Model 1: art ~ fem + mar + kid5 + phd + ment
## Model 2: art ~ fem + mar + kid5 + phd + ment
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 9 -1574.2
## 2 12 -1570.8 3 6.8037 0.07842 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1art.poissonc.ran (相關性隨機參數模型) 顯著優於 art.poissonH.ran (異質性隨機參數模型) 而 art.poissonc.ran 略優於(10% 顯著性) art.poisson.ran(隨機參數模型).
9 條件平均數
模型參數的估計提供了參數向量(parameter vectors)的非條件估計(unconditional estimate)，但我們可以推導出個體特定的條件估計（見Train 2009；Greene 2012）。 以特定人的資料為條件的隨機參數分佈條件平均值的估計是：
\[\hat{\bar \beta_i} = E(\beta_i \mid \text{data}_i) = \sum_{r=1}^R \left ( \frac{ \hat P(y \vert X ,\beta_i )} {\sum_{r=1}^R \hat P(y \vert X ,\beta_i )} \right ) \hat{\beta_{ir}} \]
\[ \hat \beta_{ir} = \hat \beta_i+ \hat \Pi s_i + \hat L \omega_{ir}\].
注意，ˋ這不是 \(\beta_i\) 的實際估計值，而是隨機參數分佈的條件平均值的估計值（Greene et al. 2014）。 我們還可以透過估計以下方式估計該分佈的變異數：
\[\hat{\bar{\beta_i}}^2 = E(\beta_i^2 \mid \text{data}_i) = \sum_{r=1}^R \left ( \frac{ \hat P(y \vert X ,\beta_i )} {\sum_{r=1}^R \hat P(y \vert X ,\beta_i )} \right ) \hat \beta_{ir}^2 \]
標準差
\[ \sqrt{ E(\beta_i ^2\mid \text{data}_i)}\]
透過條件平均值和條件方差的估計，我們可以計算類似於置信區間(credbile interval)的區間的極限，即平均正負兩個估計標準差。 這將構建一個區間，其中包含至少95%的 \(\beta_i\) 條件分佈（Greene 2012）。
劃製 \(\beta_{wage}\) 直方圖和核密度(kernel density)。
## Plot the distribution of the conditional mean for wage
plot(union.ran, par = "wage", type = "density")
## Plot the conditional mean for the first 20 individuals
plot(union.ran, par = "wage", ind = TRUE, id = 1:20, col = "blue",
main=expression(paste('The first 20 ', E( beta[i]),' credible interval ')))
10 補償變量 (compensenting variations)
令 \(x_{il}\) 表示 \(x_i\) 中的第 l 個元素，\(\beta_l\) 表示相應的參數，而\(x_{im}\) 表示 \(x_i\) 中的第 m 個元素，\(\beta_m\) 表示相應的參數。 現在考慮同時改變 \(x_{il}\) 和 \(x_{im}\)，使得 \(y_i^∗ = 0\)(效用不變)。這需要
\[ 0 = \beta_{il} \Delta x_{il} + \beta_{im} \Delta x_{im} \Rightarrow \frac{\delta x_{il}}{\delta x_{im}} = -\frac{\beta_{im}} { \beta_{il}}\]
由於 \(\beta_{il}, \beta_{im}\) 是隨機的，因此係數的比率是隨機的，並且其分佈遵循係數的聯合分佈。
## Plot the compensating variation
plot(union.ran, par = "wage", wrt = "school", effect = "cv",
main='Compensating variation for wage to school')
10.1 補償變化(compensating variation)檢驗:
補償變化檢驗(非線性假設檢定)，即 phd/ment(減少一篇指導教授的論文,需要增加多少的系所聲望, 才能夠維持博士畢業生的發表論文篇數不變) , 是否為零
\[ H_0：\frac{\text{phd}}{\text{ment}} = 0\]
library("car")## 載入需要的套件：carData##
## 載入套件：'car'## 下列物件被遮斷自 'package:memisc':
##
## recode knitr::kable(deltaMethod(art.poisson.fixed, "phd/ment"),format='simple',
caption='Compensating Variation Test:Fixed Parameters')| Estimate | SE | 2.5 % | 97.5 % | |
|---|---|---|---|---|
| phd/ment | 0.5020048 | 1.043031 | -1.542298 | 2.546307 |
固定參數模型: 補償變量估計值大概是0.5, 標準差1.04, 在 95% 信心區間介於-1.5到2.5之間,因此不拒絕虛無假設為零的狀況.
因此,減少一篇指導教授的論文,不需要增加系所聲望, 也能夠維持博士畢業生的發表論文篇數.
knitr::kable(deltaMethod(art.poisson.ran, "mean.phd/mean.ment"),format='simple',
caption='Compensating Variation Test: Random Parameters ')| Estimate | SE | 2.5 % | 97.5 % | |
|---|---|---|---|---|
| mean.phd/mean.ment | -0.962443 | 1.176899 | -3.269122 | 1.344236 |
隨機參數模型: 補償變量估計值大概是 -0.23, 標準差1.22, 在 95% 信心區間介於-2.6到2.15之間,因此不拒絕虛無假設為零的狀況.
因此,減少一篇指導教授的論文,不需要增加系所聲望, 也能夠維持博士畢業生的發表論文篇數.