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