文 / Chuck Huber, Director of Statistical Outreach at StataCorp
導讀
在之前的文章中，作者演示了如何利用蒙特卡羅模擬計算功率的t檢驗的示例，以及如何把你的模擬集成到Stata的power命令中。本文，作者將向您展示如何為線性和邏輯回歸模型執行這些任務。線性回歸和邏輯回歸程序的策略和整體結構類似于t檢驗示例。不同點在于數據模擬和用于檢驗零假設的模型。
在模擬回歸模型時，選擇現實的回歸參數是很具有挑戰性的。有時，pilot data 或historical data可以提供思路，但通常我們會考慮一系列我們認為有意義的參數值。本文將通過基于全國健康和營養檢查調查 (NHANES) 數據做為示例?？梢酝ㄟ^鍵入webuse nhanes2來下載這些數據的一個版本。
線性回歸示例
假設您正在計劃一項systolic blood pressure (SBP) 研究，并且您認為年齡和性別之間存在相互作用。NHANES 數據集包括變量bpsystol (SBP)、age和sex。下面，我們擬合了一個線性回歸模型，其中包括一個age與sex的交互項，并且所有參數估計的p值都等于 0.000。這并不奇怪，因為該數據集包含 10,351個觀測值。當其他一切保持不變時，隨著樣本量變大，p值變小。
也許您沒有資源為您的研究收集 10,351名參與者的樣本，但如果希望有 80%的功效來檢測0.35的交互參數。您的樣本需要多大？
讓我們首先根據 NHANES 模型的參數估計創建單個偽隨機數據集。我們通過清除Stata的內存開始下面的代碼塊。接下來，我們將隨機種子設置為 15，以便我們可以重現我們的結果并將觀察次數設置為 100。
代碼塊的第四行生成一個名為age的變量 ，其中包括從區間[18,65]上的均勻分布中提取的整數。
第五行使用概率等于0.5的 Bernoulli distribution 生成一個名為female的指示變量?；叵胍幌?，一次試驗的二項式分布等價于 Bernoulli distribution 。
第六行生成了 age和 female交互作用的變量。
第七行生成一個變量e，即回歸模型的誤差項。生成了誤差由均值為 0且標準差為 20的正態分布。值 20基于從NHANES回歸模型估計的root MSE。
代碼塊的末行基于我們的模擬變量和來自NHANES回歸模型的參數估計的線性組合生成變量sbp。
以下是使用regress擬合我們的模擬數據的線性模型的結果。參數估計與我們的輸入參數有些不同，因為我只生成了一個相對較小的數據集。我們可以通過增加樣本大小、抽取大量樣本或兩者兼而有之的方式來減少這種差異。
交互作用項的p值等于 0.420，這在 0.05水平上不具有統計顯著性。顯然，我們需要更大的樣本量。
我們可以使用回歸模型中的p值來檢驗交互項為零的原假設。這在本例中是可行的，因為我們只測試一個參數。但是，如果我們的交互包含一個分類變量，例如race，我們將不得不同時測試多個參數。有時我們希望同時測試多個變量。
Likelihood-ratio tests可以檢驗多種假設，包括同時檢驗多個參數。我將在本示例中向您展示如何使用Likelihood-ratio tests，因為它將涉及到您在研究中可能遇到的其他假設。如果您不熟悉它們，您可以在Stata Base Reference Manual中閱讀有關likelihood-ratio tests的更多信息。
下面的代碼塊顯示了用于計算likelihood-ratio tests的五個步驟中的四個。我們將檢驗交互項的系數為零的原假設。首行擬合包含交互項的“full” 回歸模型。第二行將完整模型的估計值存儲在內存中?！癴ull”這個名字是任意的。我們可以給這個模型的結果命名任何我們喜歡的名字。第三行擬合了省略交互項的 “reduced” 回歸模型。第四行將簡化模型的結果存儲在內存中。
第五步使用lrtest計算完整模型與簡化模型的 likelihood-ratio test。該檢驗產生0.4089的p值，接近上述回歸輸出中報告的 Wald 檢驗。我們不能拒絕交互參數為零的原假設。
您可以鍵入return list以查看存儲在標量r(p)中的p值。您可以使用r(p)來定義reject，就像我們在t測試程序中所做的那樣。
模擬數據和檢驗回歸模型的零假設比t檢驗稍微復雜一些。但是編寫一個程序來自動化這個過程幾乎與t測試示例相同。讓我們考慮下面的代碼塊，它定義了程序simregress。
前三行以capture program、program和version開頭，與我們的t測試程序基本相同。
該程序的語法部分與t測試程序的語法部分類似，但輸入參數的名稱顯然不同。已經包含了樣本大小、alpha 級別和基本回歸參數的輸入參數。我們沒有為模型中的每個可能的參數都包含一個輸入參數，但如果你愿意，也可以這樣做。例如，在我們的程序中將變量age的范圍“hard coded”為18到 65。但是，如果您愿意，也可以包含age上限和下限的輸入參數。我還發現包含描述參數名稱的注釋很有幫助，這樣就不會產生歧義。
下一段代碼嵌入在一個 “quietly”的塊中。像set obs、generate和regress這樣的命令將輸出發送到結果窗口和日志文件（如果你打開了）。將這些命令放在一個quietly塊中會抑制該輸出。
我們已經編寫了創建隨機數據和檢驗零假設的命令。因此，我們可以將該代碼復制到quietly塊中，并用語法定義的相應本地宏替換任何輸入參數。例如，我已將set obs 100更改為set obs 'n'，以便觀察的數量將由語法指定的輸入參數設置。我還為輸入參數指定了與模型中模擬變量相同的名稱。所以'age'*age是syntax定義的輸入參數'age'和模擬生成的變量age的乘積。
likelihood-ratio test的p值存儲在標量r(p)中，我們的程序返回標量reject ，與在我們的t測試程序中完全一樣。
下面，我使用simulate運行simregress 100次并summarized變量reject。結果表明，給定 100名參與者的樣本和關于模型的其他假設，我們將有16%的功效來檢測 0.35的交互參數。
接下來，我們編寫一個名為power_cmd_simregress的程序，以便我們可以將simregress集成到 Stata 的power命令中。power_cmd_simregress的結構與上一篇文章中的power_cmd_ttest相同。首先，我們定義語法和輸入參數并指定它們的默認值。然后，我們運行模擬并總結變量reject。之后，我們返回結果。
我們還編寫一個名為power_cmd_simregress_init的程序?；叵胍幌律弦黄恼?，該程序將允許我們為一系列輸入參數值運行power simregress，包括雙引號中列出的參數。
現在，我們已經準備好使用power simregress了！下面的輸出顯示了當interaction參數等于 0.2到0.4時的模擬功效，增量為0.05，對于大小為 400、500、600和700的樣本。
圖 1 以圖形方式顯示了結果。
圖 1：回歸模型中交互項的估計功效
表格和圖表向我們展示了可以產生 80%功效的幾種參數組合。700名參與者的樣本將使我們有大約 80%的功效來檢測0.30的interaction參數。600名參與者的樣本將使我們有大約80%的能力來檢測 0.33的interaction參數。500名參與者的樣本將使我們有大約80%的功效來檢測大約 0.37的interaction參數。400名參與者的樣本將使我們有大約 80%的能力來檢測 0.40的interaction參數。我們對樣本大小的選擇基于我們想要檢測的interaction參數的大小。
此示例側重于具有兩個協變量的回歸模型中的交互項。但是你可以修改這個例子來模擬你能想象到的幾乎任何類型的回歸模型的能力。在規劃模擬時，我建議執行以下步驟：
寫下感興趣的回歸模型，包括所有參數。
指定協變量的詳細信息，例如age 范圍或females比例。
找到或考慮模型中參數的合理值。
假設替代假設，模擬單個數據集，并擬合模型。
編寫一個程序來創建數據集、擬合模型并使用simulate 來測試程序。
編寫一個程序命名為power_cmd_mymethod，它允許運行您的power模擬。
編寫一個名為power_cmd_ mymethod_init的程序，以便您可以將 numlists 用于所有參數。
讓我們嘗試使用這種方法進行邏輯回歸模型。
邏輯回歸示例
在此示例中，讓我們假設您正計劃對hypertension ( highbp )進行研究。hypertension是二元的，因此我們將使用邏輯回歸來擬合模型并使用優勢比作為效應大小。
第 1 步：寫下模型
模擬功率的首步是寫下模型。logit(highbp)=β0+β1(age)+β2(sex)+β3(age×sex)
我們需要創建變量highbp、age、sex和交互項age×sex。我們還需要指定合理的參數值β0,β1, β2， 和 β3.
步驟 2：指定協變量的詳細信息
接下來，我們需要考慮模型中的協變量。什么樣的年齡值對我們的研究是合理的？我們對老年人感興趣嗎？年輕的成年人？假設我們對 18到 65歲之間的成年人感興趣。年齡分布在區間 [18,65]內是否可能是均勻的，或者我們是否期望在年齡中間出現駝峰狀分布范圍？我們還需要考慮研究中男性和女性的比例。我們是否可能對 50%的男性和 50%的女性進行抽樣？這些是我們在規劃功率計算時需要問自己的問題。
假設我們對 18到65歲之間的成年人感興趣，并且我們相信年齡是均勻分布的。我們還假設樣本將是 50%的女性。一旦我們為age和sex創建變量，就很容易計算交互項age×sex。
步驟 3：為參數指定合理的值
接下來，我們需要考慮模型中參數的合理值。我們可以根據文獻綜述、試點研究結果或公開數據來選擇參數值。
我選擇再次使用 NHANES 數據，因為它包括變量hypertension ( highbp )、age和sex。
輸出模型中每個變量的優勢比估計值。優勢比是指數參數估計(that is,
），所以我們可以指定優勢比的自然對數
作為我們功率模擬中的參數。例如，上面輸出中年齡的優勢比的估計值為 1.04，因此我們可以指定
我們也可以指定
和
第 4 步：Simulate a dataset assuming the alternative hypothesis, and fit the model
接下來，我們根據我們對替代假設下的模型的假設創建一個模擬數據集。下面的代碼塊與我們用來為線性回歸模型創建數據的代碼幾乎相同，但有兩個重要的區別。首先，我們使用generate xb創建參數和模擬變量的線性組合。參數表示為使用NHANES數據估計的優勢比的自然對數。其次，我們使用rlogistic(m,s)從變量xb創建二進制因變量highbp。
然后，我們可以將邏輯回歸模型擬合到我們的模擬數據中。
第 5 步：編寫程序以創建數據集、擬合模型并使用模擬測試程序
接下來，讓我們編寫一個程序，在 alternative hypothesis下創建數據集，擬合邏輯回歸模型，檢驗零假設，并使用simulate運行程序的多次迭代。
下面的代碼塊包含名為simlogit的程序的語法。語法命令中的默認參數值是我們使用 NHANES數據估計的優勢比。我們使用lrtest來檢驗零假設，即age×sex的優勢比等于 1。
然后，我們使用simulate默認參數值運行simlogit100 次。
simulate將假設檢驗的結果保存到名為reject的變量中。reject的平均值是假設樣本量為 500 人時，我們估計檢測age×sex交互項的優勢比為 1.03的能力。
第 6 步：編寫一個名為power_cmd_simlogit的程序
如果我們只對一組特定的假設感興趣，我們可以停止我們的快速模擬。但是編寫一個名為power_cmd_simlogit的附加程序很容易，它允許我們使用 Stata 的power命令為一系列樣本大小創建表格和圖形。
第 7 步：編寫一個名為 power_cmd_simlogit_init的程序
編寫一個名為power_cmd_simlogit_init的程序也很容易，它允許我們為模型中參數的一系列值模擬功率。
使用power simlogit
現在，我們可以使用power simlogit來模擬各種假設的功率。下面的示例模擬了一系列樣本大小和效果大小的功效。樣本大小從400到 1000人不等，以 200為增量。以及age×sex交互項的優勢比范圍從 1.02到 1.05，增量為 0.01。
圖 2：邏輯回歸模型中交互項的估計功效
上面的表格和圖表表明，使用樣本大小和效應大小的四種組合可實現 80%的功效。鑒于我們的假設， 我們將至少有 80%的功效來檢測 600、800和 1000的樣本量的優勢比為 1.04。我們將有80%的功效來檢測1.05的優勢比，樣本量為 400人。
在這篇文章中，我們展示了如何在線性和邏輯回歸模型中模擬交互項的統計功效。您可以根據自己的目的修改上面的示例。
Stata軟件訂購：
如需訂購Stata V17全新版軟件，請聯系Stata中國授權經銷商及合作伙伴北京友萬信息科技有限公司（www.uone-tech.cn）。我司擁有強大的售后服務團隊，聚合國內一線Stata行業專家為客戶提供優質的技術支持服務，并幫助中國用戶建立完善的軟件服務體系。手機/微信：18610597626 郵箱：crystal@uone-tech.cn。
專注分享商業數據分析、金融數據分析、應用統計分析、知識圖譜、機器學習、計量經濟、人工智能、網絡爬蟲、自動化報告與可重復研究等熱門技術內容。定向培養Stata、Python、R語言數據人才，助力產學研政企商協同發展，為中國大數據產業蓄能。合作熱線：010-56548231 郵箱：info@uone-tech.cn。
