使用蒙特卡羅模擬計(jì)算功率，*1部分: 基礎(chǔ)知識

時間：2019-12-20作者：北京天演融智軟件有限公司瀏覽：182

功率和樣本大小的計(jì)算是科學(xué)研究計(jì)劃的重要組成部分?？梢允褂肧tata的power命令來計(jì)算許多常用統(tǒng)計(jì)測試的功率和樣本大小需求。但對于較復(fù)雜的模型是沒有簡單公式的，如多層/縱向模型和結(jié)構(gòu)方程模型(SEMs)。蒙特卡羅模擬是計(jì)算復(fù)雜模型的功率和樣本大小要求的一種方法，Stata提供了執(zhí)行此操作所需要的所有工具。甚至可以將模擬集成到Stata的power命令中，可以輕松地為一系列參數(shù)值創(chuàng)建自定義表格和圖形。




比如，下面的自定義程序power simmixed模擬縱向模型的功率，假設(shè)參與者的數(shù)量（級別2）為100到500，每次增加100，每個參與者有5到6個觀察值(級別1)。power simmixed還創(chuàng)建了一個表(未顯示)和下圖來顯示模擬的結(jié)果。




我的同事和我寫了一系列文章展示如何做到這一點(diǎn)。在今天的文章中，將介紹使用模擬計(jì)算功率和樣本大小要求所需的基本工具。在*二篇文章中，將展示如何將模擬集成到Stata的power命令中。 然后，我們將展示線性回歸，邏輯回歸，多級/縱向模型和結(jié)構(gòu)方程模型的具體示例。

基本思路

統(tǒng)計(jì)功率是當(dāng)零假設(shè)為假時拒絕零假設(shè)的概率。功率的計(jì)算是基于一組假設(shè)，例如樣本大小，alpha水平和特定的替代假設(shè)。例如，我們可能希望計(jì)算t檢驗(yàn)的功率，假設(shè)零假設(shè)的樣本均值為70，替代假設(shè)為75，樣本大小為100，alpha水平為0.05。

使用蒙特卡羅模擬計(jì)算功率的基本步驟是

1. 生成假設(shè)替代假設(shè)為真(例如，均值=75)的數(shù)據(jù)集。

2. 使用數(shù)據(jù)集測試零假設(shè)（例如，測試均值= 70）。

3. 保存測試結(jié)果（例如，“拒絕”或“未拒絕”）。

4. 多次重復(fù)步驟1-3（通常為1,000或更多）。

零假設(shè)被拒絕的比例是我們對統(tǒng)計(jì)能力的估算。在上面的例子中，我們可能在1000次迭代中觀察到834次“拒絕”，這使我們的估算功率為0.834或83.4％。

要執(zhí)行這些步驟，需要熟悉Stata的一些編程工具。以下是本文中介紹的主題列表。 如果您熟悉其中一些主題，可以單擊下面的鏈接跳到不熟悉的主題。

主題列表

標(biāo)量和局部宏

創(chuàng)建偽隨機(jī)數(shù)據(jù)集

存儲模型輸出

如何使用程序創(chuàng)建簡單程序

如何使用程序創(chuàng)建有用的程序

如何使用模擬多次運(yùn)行程序

用power onemean來檢查結(jié)果

標(biāo)量和局部宏

標(biāo)量和本地宏是模擬的重要工具，因?yàn)樗鼈冊试S將數(shù)字臨時存儲在內(nèi)存中。例如，可以通過輸入將數(shù)字1存儲到名為i的標(biāo)量中




稍后通過輸入來引用此標(biāo)量







如果對標(biāo)量和變量使用相同的名稱，Stata將會變得混亂?？梢酝ㄟ^為標(biāo)量使用一的名稱來避免這種混淆，或者可以使用scalar()函數(shù)來引用標(biāo)量。




還可以使用本地宏存儲數(shù)字。 例如，可以通過輸入將數(shù)字1存儲到名為i的本地宏。







然后，可以通過在其**引用左單引號（通常位于“1”鍵左側(cè)的鍵盤上）并使用右單引號（通常位于Enter鍵左側(cè)）來引用本地宏。）。




本地宏通常用于定義模擬的輸入?yún)?shù)，許多Stata命令的結(jié)果存儲為標(biāo)量。

創(chuàng)建偽隨機(jī)數(shù)據(jù)集

您還需要生成隨機(jī)數(shù)來進(jìn)行模擬。在這里，將展示一些常用的隨機(jī)數(shù)函數(shù)，并參考Stata Functions Reference Manual以獲取完整列表。首先清除Stata的內(nèi)存，將隨機(jī)數(shù)種子設(shè)置為15，然后使用set obs告訴Stata我們想要創(chuàng)建一個包含200個觀察值的數(shù)據(jù)集。




可以使用runiform()函數(shù)在區(qū)間(0,1)內(nèi)生成均勻分布的隨機(jī)數(shù)。




可以使用類似的runiformint(a,b)函數(shù)在區(qū)間[a,b]上生成均勻分布的隨機(jī)整數(shù)變量。例如，可以輸入以下命令，在區(qū)間[18,65]內(nèi)生成一個人的年齡的隨機(jī)值:




可以使用rbinomial（n，p）函數(shù)生成二項(xiàng)式(n,p)隨機(jī)變量，其中n是試驗(yàn)次數(shù)，p是成功概率。例如，可以通過輸入以下命令為女性生成隨機(jī)指示器：







還可以使用rnormal(m,s)函數(shù)從正常密度生成隨機(jī)值，其中均值等于m，標(biāo)準(zhǔn)差等于s。 例如，可以使用以下命令生成重量和高度的變量。在這里，*了72千克的平均值和15的標(biāo)準(zhǔn)偏差，以及170厘米的平均值和10的標(biāo)準(zhǔn)偏差。




在這個例子中，我們獨(dú)立地生成了weight和height。但是像weight和height這樣的變量很可能是相關(guān)的，可以使用drawnorm來生成相關(guān)的變量。

在下面的示例中，平均值存儲在矩陣m中，標(biāo)準(zhǔn)偏差存儲在矩陣s中，矩陣C中存儲變量之間的相關(guān)性。然后，可以將這些矩陣作為參數(shù)包含在drawingorm選項(xiàng)中以創(chuàng)建相關(guān)變量height和weight.。







上述均值和標(biāo)準(zhǔn)差的估算值與我們的輸入?yún)?shù)相似，下面的估算相關(guān)系數(shù)為0.5049，接近我們在上面的相關(guān)矩陣中*的0.5的值。




存儲模型輸出

許多Stata命令將其結(jié)果存儲為標(biāo)量，宏和矩陣。運(yùn)行命令后，可以通過輸入return list查看存儲結(jié)果的列表。還可以在估算命令（如regress）后輸入ereturn list。

在下面的例子中，我使用ttest命令來測試平均權(quán)重等于70的零假設(shè)。







輸入return list顯示Stata存儲在內(nèi)存中的標(biāo)量列表。




上述均值和標(biāo)準(zhǔn)差的估算值與我們的輸入?yún)?shù)相似，下面的估算相關(guān)系數(shù)為0.5049，接近我們在上面的相關(guān)矩陣中*的0.5的值。







存儲輸出模型

許多Stata命令將其結(jié)果存儲為標(biāo)量，宏和矩陣。運(yùn)行命令后，可以通過輸入return list查看存儲結(jié)果的列表。還可以在估算命令（如regress）后鍵入ereturn list。

在下面的示例中，我使用ttest命令來測試平均權(quán)重等于70的零假設(shè)。




輸入return list會顯示Stata存儲在內(nèi)存中的標(biāo)量列表。







可以將這些標(biāo)量中的任意一個存儲到另一個標(biāo)量中。例如，可以通過輸入將存儲在r(p)中的雙面p值存儲到名為pvalue的標(biāo)量中




如何使用程序創(chuàng)建一個簡單的程序

您將運(yùn)行生成隨機(jī)數(shù)據(jù)的命令，并在運(yùn)行模擬時多次測試零假設(shè)。還需要為模擬定義輸入?yún)?shù)并返回假設(shè)檢驗(yàn)的結(jié)果。有效地完成這些任務(wù)的一種方法是使用program定義自己的Stata程序。下面的代碼塊定義了一個名為myprogram的程序，它接受輸入?yún)?shù)n()，顯示n的值，并返回n的值。




讓我們逐行考慮這個代碼塊。

**行是capture program drop myprogram。定義程序后，必須先使用program drop將程序從內(nèi)存中刪除，然后才能修改和重新定義程序。因?yàn)楹芸赡茉谕瓿芍皶啻涡薷拇顺绦?，所以在重新定義程序之前，它將節(jié)省運(yùn)行program drop的時間。**次運(yùn)行此代碼塊時不會定義程序，所以program drop將返回一個錯誤。在program drop之前輸入capture將捕獲此錯誤并允許代碼塊繼續(xù)運(yùn)行。

*二行是程序program myprogram, rclass。 開始定義程序myprogram。 定義程序后，可以輸入myprogram，Stata將在program和end之間運(yùn)行所有命令。選項(xiàng)rclass告訴Stata我們想要使用return返回程序中的值。

*三行，15.1版本告訴Stata希望使用Stata 15.1中編寫的功能來運(yùn)行程序。 可以在Stata Programming Reference Manual中了解有關(guān)Stata版本控制的更多信息。

*四行是syntax, n(integer)。 定義了程序的語法。該程序要求用戶輸入逗號，后跟n()的整數(shù)。n()的值存儲在名為n的本地宏中，可以將此本地宏稱為“n”。

*五行顯示本地宏n的值。

*六行是return scalar N = `n’。此行指示程序以名為N的標(biāo)量形式返回n的值。
最后一行是end。這告訴Stata你已經(jīng)完成了程序myprogram的定義。

讓我們運(yùn)行myprogram，看看它做了什么。




為n()輸入值50，Stata顯示結(jié)果n = 50。
可以輸入return list，并看到myprogram以標(biāo)量r（N）的形式返回n的值。




非常有用！

如何使用程序創(chuàng)建有用的程序

讓我們定義一個名為simttest的程序，它根據(jù)輸入?yún)?shù)生成一個隨機(jī)數(shù)據(jù)集，測試零假設(shè)，并返回我們的假設(shè)檢驗(yàn)結(jié)果。 下面的代碼塊僅使用注釋來定義程序。




現(xiàn)在，讓我們在下面的代碼塊中填寫詳細(xì)信息。




syntax的定義包括一個必需參數(shù)，n()和四個可選參數(shù)，alpha(),m0(),ma(),和sd()。 可選參數(shù)包含在方括號中。所有輸入?yún)?shù)都包含對輸入值的限制。例如，n()必須是整數(shù)，可選參數(shù)必須是實(shí)數(shù)?？蛇x參數(shù)還包括一個默認(rèn)值，如果用戶在輸入程序名稱時未*輸入?yún)?shù)，則使用該值。例如，除非用戶*不同的值，否則將為m0分配值0。

注意，在此示例中，n()是樣本大小，alpha()是alpha級別，m0()是假設(shè)零假設(shè)的平均值，ma()是假設(shè)備選假設(shè)的平均值，sd()是標(biāo)準(zhǔn)偏差。

然后drawnorm使用輸入?yún)?shù)從均值為`ma’、標(biāo)準(zhǔn)差為`sd’的正態(tài)分布中生成一個觀測值為`n’的樣本。接下來，ttest測試樣本均值等于`m0′的原假設(shè)。

可以通過輸入return scalar reject =（r（p）<`alpha'）來返回假設(shè)檢驗(yàn)的結(jié)果。 回想一下ttest將雙邊p值存儲在標(biāo)量r(p)中。 當(dāng)r(p)小于`alpha’*的alpha級別時，我們的程序使用標(biāo)量拒絕返回值1，否則返回0。

現(xiàn)在，可以輸入simttest和輸入?yún)?shù)來運(yùn)行模擬。







通過輸入return list，可以看到模擬的結(jié)果。







如何使用模擬多次運(yùn)行程序

simttest程序?qū)?zhí)行模擬的一次迭代所需的所有操作。 接下來，需要一種方法來多次運(yùn)行simttest并收集結(jié)果。下面的代碼塊顯示了如何使用simulate來完成這兩項(xiàng)任務(wù)。




參數(shù)reject=r(reject)告訴simulate將r(reject)中返回的結(jié)果保存到名為reject的變量中。 選項(xiàng)reps(100)指示模擬運(yùn)行程序100次。 選項(xiàng)seed(12345)設(shè)置隨機(jī)數(shù)種子，以便我們的結(jié)果可以重現(xiàn)。

冒號之后是simttest以及我們模擬的輸入?yún)?shù)。 某些模擬需要很長時間才能運(yùn)行，并且simulate會在結(jié)果窗口中顯示一個點(diǎn)，以便知道它仍在運(yùn)行。下面的輸出顯示了我們的模擬結(jié)果。







simulate將結(jié)果保存到變量reject，如果零假設(shè)的測試被拒絕則包含1，否則為0。




可以使用summarize來計(jì)算reject的平均值，該平均值等于100次迭代中拒絕原假設(shè)的次數(shù)所占的比例。根據(jù)輸入?yún)?shù)，該比例是對統(tǒng)計(jì)功效的估算！







在這個例子中，比例等于0.91，這意味著當(dāng)檢驗(yàn)樣本平均值等于75的替代假設(shè)時，如果樣本平均值等于70，假設(shè)標(biāo)準(zhǔn)偏差為15且樣本量為100，可以期望91%的功率。

用Power OneMean檢查結(jié)果

可以使用power onemean檢查蒙特卡羅模擬的結(jié)果。 以下示例包含與我們的模擬相同的輸入?yún)?shù)。




power onemean計(jì)算的功率為0.9100，這與我們的模擬估算的功率相同。因?yàn)槟M是隨機(jī)的，所有結(jié)果并不總是**的。改變隨機(jī)數(shù)種子或迭代次數(shù)可以稍微改變估算的功率。但結(jié)果應(yīng)該很接近。

總結(jié)

在本文中，介紹了使用蒙特卡羅模擬計(jì)算統(tǒng)計(jì)功率和樣本大小需求所需的工具。 下一次，將向您展示如何使用Stata的power命令來運(yùn)行模擬，以便可以輕松地為一系列輸入?yún)?shù)創(chuàng)建表格和圖形。








北京天演融智軟件有限公司專注于matlab,stata,spsspro,pscad,gams,nvivo,cyme,gms,mathematica等

• 詞條

  詞條說明

• SAS商業(yè)智能分析軟件

  SAS 9.4商業(yè)智能分析軟件由科學(xué)軟件網(wǎng)提供，同時還提供SAS軟件視頻教程。 SAS 9.4滿足每個企業(yè)的需求，高性能分析，云環(huán)境部署，標(biāo)準(zhǔn)化數(shù)據(jù)管理。不僅提供所需解決方案，而且SAS滿足您的較高預(yù)期。 構(gòu)建 安裝開通速度比以往較快。統(tǒng)一的SAS9.4平臺可降低部署成本并縮短安裝、管理和維護(hù)的時間。 簡化 全面簡化用戶體驗(yàn) — 從部署、管理到一對一用戶交互，消除技術(shù)和技能水平的障礙，每個人都可

• Litera Compare 文檔比對軟件

  多種文檔類型的總體比較Litera Compare 在幾秒鐘內(nèi)準(zhǔn)確地檢測到任何兩個文檔的變化。這款市場良好的文檔比較解決方案可融入您的工作生活，無論您是在臺式機(jī)、云端還是在旅途中工作，你都不會**任何機(jī)會。?Litera Compare?幫助您專注于重要的事情以多種文檔格式管理多個作者的修訂所需的時間可能會成為一個嚴(yán)重的成本問題，較不用說用肉眼來發(fā)現(xiàn)較新的內(nèi)容，不僅費(fèi)時，還容易

• 2023年3-6月Stata軟件系列在線研討會

  Stata是一個完整的、集成的軟件包，它提供了所有數(shù)據(jù)科學(xué)所需的數(shù)據(jù)操作、可視化、統(tǒng)計(jì)和自動化報告功能。為了讓大家較好的掌握Stata軟件，StataCorp LLC為大家安排了一系列免費(fèi)的在線培訓(xùn)。主題：Tips and tricks時間：27 April 2023, 10:00 AM CT? ? ? ?發(fā)現(xiàn)Stata的一些隱藏功能。 在“技巧和竅門”網(wǎng)絡(luò)

• KaleidaGraph 科學(xué)繪圖軟件介紹

  KaleidaGraph是一個精心設(shè)計(jì)的圖形和數(shù)據(jù)分析應(yīng)用軟件，學(xué)術(shù)研究的科學(xué)家，以及商業(yè)和工程領(lǐng)域的科學(xué)家都可以用。它可以生成公開質(zhì)量的圖表，并且輕松地將較復(fù)雜的數(shù)據(jù)轉(zhuǎn)換成功能顯示。KaleidaGraph允許用戶導(dǎo)入、操作和分析數(shù)據(jù)，并創(chuàng)建定制的圖。統(tǒng)計(jì)、線性和非線性曲線擬合，以及對數(shù)據(jù)進(jìn)行精確的圖形可視化的能力，使KaleidaGraph具有強(qiáng)大的靈活性。靈活的數(shù)據(jù)條目當(dāng)KaleidaGra