原文鏈接：http://tecdat.cn/?p=23759
原文出處：拓端數據部落公眾號
簡介
兩階段最小二乘法（2SLS）回歸擬合的線性模型是一種常用的工具變量估計方法。
本文的主要內容是將各種標准的回歸診斷擴展到2SLS。
2SLS估計的回顧
我們需要2SLS回歸的一些基本結果來開發診斷方法，因此我們在此簡單回顧一下該方法。2SLS回歸是由Basmann(1957)和Theil(引自Theil 1971)在20世紀50年代獨立發明的，他們采取了略微不同但又相當的方法，都在下面描述，以得出2SLS估計器。
我們想估計線性模型y=Xβ+ε，其中y是因變量的n×1觀察向量，X是回歸因子的n×p矩陣，通常初始列1s為回歸常數。 β是一個p×1的回歸系數向量，需要根據數據進行估計，ε是一個n×1的誤差向量，假定其分布為Nn(0,σ2In)，其中Nn是多變量正態分布，0是一個n×1的零向量，In是n階單位矩陣。假設X中的一些（也許是全部）回歸因子是內生的，即它們被認為不獨立於ε的意義。因此，β的普通最小二乘法（OLS）估計值bOLS=(X⊤X)-1X⊤y通常是有偏的，而且不一致。
現在假設我們有另一組獨立於ϵ的q工具變量（IVs）Z，其中q≥p。如果q=p，我們可以直接應用IV來估計β，但如果q>p，我們有更多的IV，而不是我們需要的。簡單地拋棄IVs是低效的，2SLS回歸是一個通過合理的方式將IVs的數量減少到p的程序。
2SLS的第一階段通過多元普通最小二乘法對模型矩陣X中的所有回歸變量進行回歸，得到q×p的回歸系數矩陣B=（Z⊤Z）-1Z⊤X，以及擬合值Xˆ=ZB。B的列相當於X的每一列對Z的單獨最小二乘回歸產生的系數。如果X的某些列是外生的，那么這些列也會出現在Z中，因此，XˆX^中與外生調節器有關的列只是復制了X的相應列。
由於XˆX^的列是Z的列的線性組合，它們（漸進地）與ε不相關，使它們成為估計回歸方程的合適IV。這個IV步驟是Theil方法中2SLS的第二個階段。
作為一種替代方法，我們可以通過對XˆX^進行OLS回歸來獲得完全相同的β的估計值b2SLS，產生b2SLS=(Xˆ⊤Xˆ)Xˆ⊤y。這就是巴斯曼的方法，也是 "2SLS "這個名字的由來。
無論我們把第二階段看成是IV估計還是OLS回歸，我們都可以把這兩個階段合並成一個公式。
這就是sem包中的tsls()函數(Fox, Nie, and Byrnes 2020)所做的，但是從開發回歸診斷的角度來看，通過兩個不同的OLS回歸來計算2SLS估計值是有利的。
對2Sls回歸異常-數據診斷
就我們所知，用2SLS擬合的回歸模型的診斷是一個相對被忽視的話題，但Belsley, Kuh和Welsch(1980, 266-68)簡要地討論了這個問題。刪除診斷法直接評估每個案例對擬合回歸模型的影響，方法是刪除案例，重新擬合模型，並注意到回歸系數或其他回歸輸出，如殘差標准差，如何變化。
對於有影響的數據，總是可以通過粗暴的計算來獲得案例刪除診斷，即用每個案例依次刪除來重新擬合模型，但這種方法效率低下，因此在大樣本中沒有吸引力。對於某些類別的統計模型，如廣義線性模型（如Pregibon 1981），對個案刪除診斷的計算要求較低的近似值是可用的，而對於線性模型，有效的 "更新 "公式是可用的（如Belsley, Kuh, and Welsch 1980所描述的），允許精確計算個案刪除診斷的。
事實證明，正如Belsley、Kuh和Welsch所指出的，Phillips（1977年，公式15和16）給出了2SLS回歸的精確更新公式，允許有效地計算個案選擇統計。
其中，b2SLS-i是去除第ii種情況后的2SLS回歸系數向量，以及
這里，yi是第i個案例的因變量值，x⊤ixi⊤是模型矩陣X的第i行，z⊤izi⊤是工具變量模型矩陣Z的第i行。
Belsley, Kuh和Welsch特別研究了（用我們的符號）dfbetai=b2SLS-b2SLS-i的值。他們還討論了殘差標准差s-i的刪除值。
然后，Belsley、Kuh和Welsch計算它們對擬合值(和回歸系數)影響的綜合度量dffits為
其中（如前）x⊤ixi⊤是模型矩陣X的第i行，XˆX^是第二階段回歸變量的模型矩陣。
讓
代表將y轉換為擬合值的n×n矩陣，yˆ=H∗y。在OLS回歸中，類似的量是hat矩陣H=X（X⊤X）-1X⊤。Belsley, Kuh和Welsch指出，H∗與H不同，它不是一個正交投影矩陣，將y正交地投影到X的列所跨越的子空間上。特別是，盡管H∗和H一樣，是等值的（H∗=H∗H∗），並且trace(H∗)=ptrace(H∗)=p，但H∗和H不同，是不對稱的，因此它的對角線元素不能被當作杠桿的總結性措施，也就是說，不能被當作hat值。
Belsley, Kuh和Welsch建議簡單地使用第二階段回歸的hat值。這些是H2=Xˆ(Xˆ⊤Xˆ)-1Xˆ⊤的對角線條目hi=hii。我們在下面討論一些替代方案。
除了hatvalues、dfbeta、s-i和dfits之外，還計算cook距離Di，這基本上是dfits的一個稍有不同的比例版本，它使用總體殘差標准差s來代替刪除的標准差s-i。
因為它們具有相等的方差，並且在正態線性模型下近似於t分布，所以 studentized殘差對於檢測異常值和解決正態分布誤差的假設非常有用。studentized殘差與OLS回歸相類似，定義為
其中ei=yi-x⊤ib 2SLS是第i種情況的因變量殘差。
如前所述，Belsley, Kuh, and Welsch (1980)建議使用第二階段回歸的 hatvalues。這是一個合理的選擇，但是它有可能遺漏那些在第一階段有高杠桿率但在第二階段回歸中沒有的案例。讓h(1)i代表第一階段的hatvalues，h(2)i代表第二階段的hatvalues。如果模型包括一個截距，兩組hatvalues都以1/n和1為界，但第一階段的平均hatvalues是q/n，而第二階段的平均hatvalues是p/n。為了使兩個階段的hatvalues具有可比性，我們將每個hatvalues除以其平均值，h(1∗)i=h(1)iq/n；h(2∗)i=h(2)ip/n。然后我們可以把兩階段的hatvalue定義為每種情況下兩者中較大的一個，hi=(p/n)×max(h(1∗)i,h(2∗)i)，或者定義為它們的幾何平均。
異常數據診斷
標准的R回歸模型通用方法，包括anova()（用於模型比較），predicted()用於計算預測值，model.matrix()（用於模型或第一或第二階段的回歸），print()，residuals()（有幾種），summary()，update()，和vcov()。
例子
數據在Kmenta（1986年）中用來說明（通過2SLS和其他方法）對線性聯立方程計量經濟學模型的估計。這些數據代表了經濟從1922年到1941年的年度時間序列，有以下變量。
Q,人均食品消費P，食品價格與一般消費價格的比率D, 可支配收入F, 前一年農民收到的價格與一般消費價格的比率A, 年為單位的時間
該數據集很小，我們可以對其進行檢查。
估計以下兩個方程式模型，第一個方程式代表需求，第二個代表供應。
變量D、F和A被視為外生變量，當然常數回歸因子（一列1）也是如此，而兩個結構方程中的P是內生解釋變量。由於有四個工具變量可用，第一個結構方程有三個系數，是過度識別的，而第二個結構方程有四個系數，是剛剛識別的。
外生變量的數值是真實的，而內生變量的數值是由Kmenta根據模型生成（即模擬）的，參數的假設值如下。
解決內生變量P和Q的結構方程，可以得到模型的簡化形式
Kmenta獨立地從N(0,1)中抽出20個δ1和δ2的值，然后設定ν1=2δ1和
結構方程估計如下（比較Kmenta 1986, 686）。
默認情況下，summary()會輸出2SLS回歸的三個 "診斷 "測試的結果。這些測試不是本文的重點，所以我們只對它們進行簡單的評論。
-
一個好的工具變量與一個或多個解釋變量高度相關，同時與誤差保持不相關。如果一個內生的回歸者與工具變量只有微弱的關系，那么它的系數將被不精確地估計。在弱工具的診斷測試中，我們希望有一個大的測試統計量和小的p值，Kmenta模型中的兩個回歸方程就是如此。
-
應用於2SLS回歸中，Wu-Hausman檢驗是對內生性的一種檢驗。如果所有的回歸者都是外生的，那么OLS和2SLS的估計都是一致的，並且OLS的估計更有效，但是如果一個或多個回歸者是內生的，那么OLS的估計就不一致了。大的檢驗統計量和小的p值，就像在這個例子中一樣，表明OLS估計器是不一致的，因此，2SLS估計器是首選。
-
Sargan檢驗是對過度識別的檢驗。也就是說，在一個過度識別的回歸方程中，如Kmenta的需求方程中，工具變量比要估計的系數多，工具變量有可能提供關於系數值的沖突信息。因此，大的檢驗統計量和小的Sargan檢驗的pp值表明，該模型被錯誤地指定了。在這個例子中，盡管我們知道（通過數據的構建方式）需求方程是正確的，但我們還是偶然得到了一個適度小的pp值0.084。Sargan檢驗不適用於剛剛確定的回歸方程，其工具變量和系數的數量相等，如Kmenta的供給方程。
lm "類對象的幾個方法與產生的對象正常工作。例如，對象的plot()方法調用了相應的 "lm "方法並產生了可解釋的圖，這里是Kmenta模型中需求方程的2SLS擬合。
-
par(mfrow=c(2, 2))
-
plot(deq)
然而，在這種情況下，我們更喜歡描述的這些診斷圖的版本。
數據表現良好。例如，在第一個結構方程中，學生化殘差的QQ圖和hatvalues、學生化殘差和庫克cook距離的 "影響圖 "都是不明顯的，除了幾個高杠桿但在一起的案例。
qqPlot
influence
影響圖中的圓圈面積與Cook's D成正比，水平線畫在學生化殘差標度的0和±2處（rstudent=2處的水平線不在圖中），垂直線在2×h¯和3×h¯處。
為了產生一個更有趣的例子，我們將把高杠桿的第20種情況（即1941年）的QQ值從Q20=106.232改為Q20=95，這個值完全在數據中QQ的范圍內，但與其他數據不一致。
然后重復對第一個結構方程的2SLS擬合，並將結果與未被破壞的數據進行比較，發現回歸系數有很大變化。
compareCoefs(deq, deq1)
有問題的第20個案例（1941年）通過異常數據回歸診斷法清楚地顯示出來。
qqPlot(deq1)
outlier
influence Plot
avPlots(deq1)
去掉第20種情況，產生的估計系數接近於未被破壞的數據的系數。
compareCoefs(deq, deq1, deq1.20)
估計系數的標准誤差比原來大，因為我們現在有19個而不是20個案例，也因為解釋變量的變化減少了。
發現hatvaues的三種定義是否對這個例子產生了實際的影響，是有一定的意義的。三種hatvalues的散點圖矩陣表明，它們都產生了類似的結果。
最后，讓我們驗證一下刪除診斷的計算結果是否正確。
非線性診斷法
Cook(1993)和Cook and Croos-Dabrera(1998)系統地探討了成分、殘差圖作為非線性診斷的理論屬性。按照這些作者的說法，並把重點放在解釋變量x1上，讓我們假設響應y與x1的部分關系可能是非線性的，由部分回歸函數f(x1)表示，而y與其他xxs的部分關系是線性的，因此，數據的准確模型是。
我們不知道f()，所以改用工作模型來擬合
在我們的例子中，通過2SLS回歸，得到估計的回歸系數a′,b′1,b′2,...,b′k。Cook和Croos-Dabrera的工作表明，只要回歸估計是一致的，XXS是線性相關的。部分殘差b′1x1+e可以被繪制出來，並對x1x1進行平滑處理，以顯示f()的估計。其中e=y-(a′+b′1x1+b′2x2+⋯b′kxk)是因變量殘差。在實踐中，如果x1和其他xxs之間有很強的非線性關系，或者y與另一個與x1相關的x有非線性關系，那么分量加殘差圖就會被分解為f()的准確表示。
Fox和Weisberg（2018）將成分加殘差圖擴展到更復雜的回歸模型，例如可以包括交互作用，將偏殘差添加到預測變量效應圖中。這些圖也可以應用於由2SLS回歸擬合的線性模型。
診斷非線性：一個例子
我們再一次轉向Kmenta的數據和模型的需求方程來說明成分殘差圖，數據再一次表現良好。為一個加法回歸方程中的所有數字解釋變量構建了分量殘差圖。比如說。
crPlots(deq, smooth=list(span=1))
我們在圖中為局部加權回歸loess 平滑（Cleveland, Grosse, and Shyu 1992）設置了一個較大的跨度，因為在數據集中只有n=20個案例。跨度的默認值是2/3。在每個面板中，紅線給出的loess 平滑度與藍線給出的最小二乘線緊密匹配，藍線代表的是解釋變量方向的擬合回歸面，左邊是P，右邊是D。因此，兩種偏關系似乎都是線性的。
CERES圖(Cook 1993)，是成分加殘差圖的一個版本，它使用平滑器而不是線性回歸，因此對預測因子之間的非線性關系更加穩定。
ceresPlots(deq, smooth=list(span=1))
在當前的例子中，這是一個加性模型，我們得到的圖形與之前的基本相同，只是y軸的縮放比例不同。
plot(predictorEffects)
預測效應圖中的藍色陰影區域代表擬合的部分回歸線周圍95%的置信度包絡。
然而，假設我們對數據擬合了錯誤的模型。
-
deq2 <- update(deq, . ~ I((P - 85)^4/10^5) + D)
-
crPlots(deq2, smooth=list(span=1))
因為max(P)/min(P)=113.49/86.50=1.3的比率不比1大多少，所以我們在把變量提高到4次方之前，從P中減去一個比min(P)略小的數字，以引起擬合部分回歸曲線中的非線性。變換后的P的成分加殘差圖清楚地反映了由此產生的缺乏擬合，而D的圖仍然是合理的線性。
帶有部分殘差的預測器效應圖顯示了對同一情況的不同看法，它將P而不是轉換后的P放在橫軸上，並揭示了擬合的非線性部分回歸函數未能捕獲數據的線性模式。
-
plot(predictorEffects(deq2, residuals=TRUE),
-
partial.residuals=list(span=1))
藍線代表擬合模型，紅線代表平滑后的偏殘差；兩條線之間的差異表明缺乏擬合。
非恆定誤差方差
標准的最小二乘法非恆定方差（"異方差"）診斷法可以直接延伸到2SLS回歸中。例如，我們可以繪制殘差與擬合值的對比圖，以發現前者的變異性隨着后者的水平而變化（通常是增加）的趨勢。對於的模型中的需求方程。
-
plot(fitted(deq), rstudent(deq))
-
abline(h=0)
似乎沒有問題。
Fox(2016)提出了這個圖形的一個變體，改編了Tukey的擴散水平圖(Tukey 1977)，繪制了絕對螺距殘差的對數與擬合值的對數，假設后者是正數。如果該圖的擬合線的斜率為b，則方差穩定化的功率轉換由yλ=y1-b給出。因此，如果b>0，建議的變換是順着Tukey的冪和根的階梯，例如，λ=1-b=1/2代表平方根變換，λ=1-b=0代表對數變換，以此類推。對於模型，我們有
這表明分布有隨水平提高的輕微趨勢。λ=-2.45的轉換似乎很強，直到我們注意到QQ的值離0很遠，而且最大和最小的值Qmax/Qmin=106.23/92.42=1.15的比率接近1，所以Q-2.45幾乎是QQ的線性轉換，也就是說，實際上根本沒有轉換。
由Breusch和Pagan（1979）提出的最小二乘回歸中的非恆定誤差方差的普通分數測試，是基於模型的
其中函數g()是未指定的，變量z1,...,zs是誤差方差的預測因子。在最常見的應用中，由Cook和Weisberg（1983）獨立提出，有一個zz，即回歸的擬合值yˆ，盡管使用初級回歸中的回歸者x作為zs也很常見。測試是通過將標准化殘差的平方e2i/σˆ2回歸到zs上實現的，其中σˆ2=∑e2i/n。然后，在誤差方差不變的無效假設下，該輔助回歸的回歸平方和除以2的漸近分布為χ2s。
Breusch-Pagan/Cook-Weisberg檢驗很容易適應2SLS回歸。對於需求方程。
在這里，第一個檢驗是針對擬合值的，第二個更一般的檢驗是針對需求方程中的解釋變量的；這兩個檢驗的p值都很大，表明沒有什么證據反對恆定方差的假說。
2SLS回歸中對非恆定方差的補救方法與最小二乘回歸中的補救方法相似。
-
我們已經提出，如果誤差方差隨着響應水平的提高（或降低），並且因變量是正的，那么我們就可以通過對因變量進行冪變換來穩定誤差方差。
-
另外，如果我們知道誤差的方差與一個常數成正比，那么我們可以對2SLS估計使用反方差加權。
-
最后，我們可以在2SLS中使用系數協方差矩陣的估計（或自舉法：例如，見Davison和Hinkley 1997）來修正非恆定誤差方差的標准誤差，就像Huber（1967）和White（1980；也見Long和Ervin 2000）提出的最小二乘法回歸那樣。
對於例子，非恆定誤差方差的證據是輕微的，標准誤差與傳統的2SLS標准誤差相似，甚至比它略小。
我們將修改數據以反映非恆定誤差方差，像最初那樣從還原形式方程中重新生成數據，將內生變量P和Q表示為外生變量D、F和A的函數，以及還原形式誤差ν1和ν2。
-
-
-
w <- 3*(EQ - min(EQ) + 0.1)/(max(EQ) - min(EQ))
-
v1 <- v1*w # 非恆定方差
-
Q <- EQ + v1
-
P <- EP + v2
-
將采樣的簡化形式誤差v1與Q的期望值作圖，顯示出明顯的異質性模式。
plot(EQ, v1))
然后將需求方程與新的數據集重新匹配，我們得到
而非恆定的誤差方差在診斷中得到了明顯的反映；例如：
ncvTest(deq2)
從建議的QQ的冪變換的極值λ=-23，反映了（正如我們之前指出的）max(Q)/min(Q)並沒有比1大多少。
在我們的例子中，標准誤差與傳統的標准誤差沒有太大區別。
如前所述，bootstrapping 提供了一種替代標准誤差的方法，作為對非恆定誤差方差的修正，實現了個案再抽樣的bootstrapping ，並返回一個適合與boot包中的函數一起使用的 "boot "類對象（Davison和Hinkley 1997；Canty和Ripley 2020）。比如說。
b.deq2 <- Boot(deq2)
在這個例子中，bootstrap的標准誤差比傳統的標准誤差大。
Boostrap置信區間也可以從Boot()返回的對象中計算出來，默認情況下報告BCa（偏差校正，加速）區間（。
confint(b.deq2)
加權的2SLS回歸
假設我們修改了回歸模型y=Xβ+ε，那么現在Nn(0,σ2W-1)其中W=diag{wi}是一個已知反方差權重的n×n對角矩陣；即V(εi)=σ2/wi。如前所述，X的某些列可能與誤差ε相關，但我們有足夠的工具變量Z，它們與誤差無關。
那么加權的2SLS估計是
或者，我們可以將2SLS的兩個階段視為加權最小二乘法（WLS）問題，在每個階段都要最小化加權殘差平方和。
Phillips的2SLS回歸的更新公式也可以針對加權情況進行修改，但是更簡單的方法（是將加權2SLS問題轉換成非加權問題，通過使用W1/2=diag{wi--√}，W的Cholesky平方根將數據轉換為恆定方差。W的平方根特別簡單，因為W是對角線的。然后在Phillips的更新公式中，我們用y∗=W1/2y代替y，用X∗=W1/2X代替X，用Z∗=W1/2Z代替Z。
對於修改后的數據，我們知道需求方程的誤差方差與變量w成反比。這反映了數據的構建方式。因此，加權的2SLS估計被計算為
將求和殘差與擬合值作圖，並測試非恆定誤差方差，並不表明有異方差問題，但有一個相對較大的求和殘差，約為-3，與其他數值相比有些突出。
plot(fitted(deqw), rstudent(deqw))
Bonferroni離群值測試表明，學生化殘差並不是異常大，模型是正確的。
outlierTest
共線性關系診斷
除了異常數據診斷外，Belsley, Kuh和Welsch（1980）還簡要地將他們的共線性關系診斷方法擴展到2SLS回歸中。我們認為，這種將共線性關系同化為數值不穩定性的方法是有缺陷的，因為它考慮到了 "與截距的共線性關系"。也就是說，數值遠離0的回歸者與常數回歸者的乘積之和很大，產生了截距的大標准誤差，只是反映了截距將擬合的回歸面遠遠超出了數據范圍的事實。
Fox和Monette(1992)描述了一種基於廣義方差膨脹因子的最小二乘法擬合的線性模型中串聯性診斷的替代方法。廣義方差膨脹因子采用了系數的估計協方差矩陣，一般適用於有線性預測因子的模型，包括由2SLS估計的線性模型。
例如，對於模型中的需求方程。
sqrt(vif(deq))
取VIF的平方根將它們放在系數標准誤差的刻度上。也就是說，P和D的系數的標准誤差比估計的系數不相關時要大23%。像這里一樣，模型中的每個項只有一個系數時，廣義和普通方差膨脹因子是一致的。P和D的VIFs相等是兩個回歸變量（超越回歸常數）的情況下所特有的。
邊際/條件圖是由car軟件包中的mcPlots()函數生成的，它將添加變量圖疊加到相應的回歸因子邊際散點圖上。這些圖使我們能夠直觀地看到由於共線性關系造成的每個系數估計精度的降低，共線性關系降低了回歸變量相對於其邊際變化的條件變化。 例如，對於需求方程。
Plots(deq)
每個面板中的藍色點代表邊際散點圖，紅色點代表（部分）添加變量圖，箭頭顯示兩組點之間的關系。
結語
仔細的回歸分析需要有效地觀察數據的方法。許多潛在的問題可以通過在擬合回歸模型之前檢查數據來解決，減少（如果不是消除）擬合后診斷的必要性。毫無疑問，采用2SLS的謹慎的數據分析者一直都是這樣做的。然而，擁有允許人們對2SLS擬合的回歸模型進行批評的方法，至少在某些情況下會建議對模型進行改進，或者對數據進行修正。
參考文獻
Basmann, R. L. 1957. “A Generalized Classical Method of Linear Estimation of Coefficients in a Structural Equation.” Econometrica 25: 77–83. https://doi.org/10.2307/1907743.
Belsley, D. A., E. Kuh, and R. E. Welsch. 1980. Regression Diagnostics: Identifying Influential Data and Sources of Collinearity. New York: John Wiley & Sons. Regression Diagnostics | Wiley Series in Probability and Statistics.
Breusch, T. S., and A. R. Pagan. 1979. “A Simple Test for Heteroscedasticity and Random Coefficient Variation.” Econometrica 47: 1287–94. https://doi.org/10.2307/1911963.
最受歡迎的見解
3.matlab中的偏最小二乘回歸（PLSR）和主成分回歸（PCR）
5.R語言回歸中的Hosmer-Lemeshow擬合優度檢驗