使用蒙特卡罗模拟计算功率，*1部分: 基础知识

时间：2019-12-20

功率和样本大小的计算是科学研究计划的重要组成部分。可以使用Stata的power命令来计算许多常用统计测试的功率和样本大小需求。但对于更复杂的模型是没有简单公式的，如多层/纵向模型和结构方程模型(SEMs)。蒙特卡罗模拟是计算复杂模型的功率和样本大小要求的一种方法，Stata提供了执行此操作所需要的所有工具。甚至可以将模拟集成到Stata的power命令中，可以轻松地为一系列参数值创建自定义表格和图形。

比如，下面的自定义程序power simmixed模拟纵向模型的功率，假设参与者的数量（级别2）为100到500，每次增加100，每个参与者有5到6个观察值(级别1)。power simmixed还创建了一个表(未显示)和下图来显示模拟的结果。

我的同事和我写了一系列文章展示如何做到这一点。在今天的文章中，将介绍使用模拟计算功率和样本大小要求所需的基本工具。在*二篇文章中，将展示如何将模拟集成到Stata的power命令中。然后，我们将展示线性回归，逻辑回归，多级/纵向模型和结构方程模型的具体示例。

基本思路

统计功率是当零假设为假时拒绝零假设的概率。功率的计算是基于一组假设，例如样本大小，alpha水平和特定的替代假设。例如，我们可能希望计算t检验的功率，假设零假设的样本均值为70，替代假设为75，样本大小为100，alpha水平为0.05。

使用蒙特卡罗模拟计算功率的基本步骤是

1. 生成假设替代假设为真(例如，均值=75)的数据集。

2. 使用数据集测试零假设（例如，测试均值= 70）。

3. 保存测试结果（例如，“拒绝”或“未拒绝”）。

4. 多次重复步骤1-3（通常为1,000或更多）。

零假设被拒绝的比例是我们对统计能力的估算。在上面的例子中，我们可能在1000次迭代中观察到834次“拒绝”，这使我们的估算功率为0.834或83.4％。

要执行这些步骤，需要熟悉Stata的一些编程工具。以下是本文中介绍的主题列表。如果您熟悉其中一些主题，可以单击下面的链接跳到不熟悉的主题。

主题列表

标量和局部宏

创建伪随机数据集

存储模型输出

如何使用程序创建简单程序

如何使用程序创建有用的程序

如何使用模拟多次运行程序

用power onemean来检查结果

标量和局部宏

标量和本地宏是模拟的重要工具，因为它们允许将数字临时存储在内存中。例如，可以通过输入将数字1存储到名为i的标量中

稍后通过输入来引用此标量

如果对标量和变量使用相同的名称，Stata将会变得混乱。可以通过为标量使用一的名称来避免这种混淆，或者可以使用scalar()函数来引用标量。

还可以使用本地宏存储数字。例如，可以通过输入将数字1存储到名为i的本地宏。

然后，可以通过在其前面引用左单引号（通常位于“1”键左侧的键盘上）并使用右单引号（通常位于Enter键左侧）来引用本地宏。）。

本地宏通常用于定义模拟的输入参数，许多Stata命令的结果存储为标量。

创建伪随机数据集

您还需要生成随机数来进行模拟。在这里，将展示一些常用的随机数函数，并参考Stata Functions Reference Manual以获取完整列表。首先清除Stata的内存，将随机数种子设置为15，然后使用set obs告诉Stata我们想要创建一个包含200个观察值的数据集。

可以使用runiform()函数在区间(0,1)内生成均匀分布的随机数。

可以使用类似的runiformint(a,b)函数在区间[a,b]上生成均匀分布的随机整数变量。例如，可以输入以下命令，在区间[18,65]内生成一个人的年龄的随机值:

可以使用rbinomial（n，p）函数生成二项式(n,p)随机变量，其中n是试验次数，p是成功概率。例如，可以通过输入以下命令为女性生成随机指示器：

还可以使用rnormal(m,s)函数从正常密度生成随机值，其中均值等于m，标准差等于s。例如，可以使用以下命令生成重量和高度的变量。在这里，*了72千克的平均值和15的标准偏差，以及170厘米的平均值和10的标准偏差。

在这个例子中，我们独立地生成了weight和height。但是像weight和height这样的变量很可能是相关的，可以使用drawnorm来生成相关的变量。

在下面的示例中，平均值存储在矩阵m中，标准偏差存储在矩阵s中，矩阵C中存储变量之间的相关性。然后，可以将这些矩阵作为参数包含在drawingorm选项中以创建相关变量height和weight.。

上述均值和标准差的估算值与我们的输入参数相似，下面的估算相关系数为0.5049，接近我们在上面的相关矩阵中*的0.5的值。

存储模型输出

许多Stata命令将其结果存储为标量，宏和矩阵。运行命令后，可以通过输入return list查看存储结果的列表。还可以在估算命令（如regress）后输入ereturn list。

在下面的例子中，我使用ttest命令来测试平均权重等于70的零假设。

输入return list显示Stata存储在内存中的标量列表。

上述均值和标准差的估算值与我们的输入参数相似，下面的估算相关系数为0.5049，接近我们在上面的相关矩阵中*的0.5的值。

存储输出模型

许多Stata命令将其结果存储为标量，宏和矩阵。运行命令后，可以通过输入return list查看存储结果的列表。还可以在估算命令（如regress）后键入ereturn list。

在下面的示例中，我使用ttest命令来测试平均权重等于70的零假设。

输入return list会显示Stata存储在内存中的标量列表。

可以将这些标量中的任意一个存储到另一个标量中。例如，可以通过输入将存储在r(p)中的双面p值存储到名为pvalue的标量中

如何使用程序创建一个简单的程序

您将运行生成随机数据的命令，并在运行模拟时多次测试零假设。还需要为模拟定义输入参数并返回假设检验的结果。有效地完成这些任务的一种方法是使用program定义自己的Stata程序。下面的代码块定义了一个名为myprogram的程序，它接受输入参数n()，显示n的值，并返回n的值。

让我们逐行考虑这个代码块。

**行是capture program drop myprogram。定义程序后，必须先使用program drop将程序从内存中删除，然后才能修改和重新定义程序。因为很可能在完成之前会多次修改此程序，所以在重新定义程序之前，它将节省运行program drop的时间。**次运行此代码块时不会定义程序，所以program drop将返回一个错误。在program drop之前输入capture将捕获此错误并允许代码块继续运行。

*二行是程序program myprogram, rclass。开始定义程序myprogram。定义程序后，可以输入myprogram，Stata将在program和end之间运行所有命令。选项rclass告诉Stata我们想要使用return返回程序中的值。

*三行，15.1版本告诉Stata希望使用Stata 15.1中编写的功能来运行程序。可以在Stata Programming Reference Manual中了解有关Stata版本控制的更多信息。

*四行是syntax, n(integer)。定义了程序的语法。该程序要求用户输入逗号，后跟n()的整数。n()的值存储在名为n的本地宏中，可以将此本地宏称为“n”。

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

*六行是return scalar N = `n’。此行指示程序以名为N的标量形式返回n的值。
最后一行是end。这告诉Stata你已经完成了程序myprogram的定义。

让我们运行myprogram，看看它做了什么。

为n()输入值50，Stata显示结果n = 50。
可以输入return list，并看到myprogram以标量r（N）的形式返回n的值。

非常有用！

如何使用程序创建有用的程序

让我们定义一个名为simttest的程序，它根据输入参数生成一个随机数据集，测试零假设，并返回我们的假设检验结果。下面的代码块仅使用注释来定义程序。

现在，让我们在下面的代码块中填写详细信息。

syntax的定义包括一个必需参数，n()和四个可选参数，alpha(),m0(),ma(),和sd()。可选参数包含在方括号中。所有输入参数都包含对输入值的限制。例如，n()必须是整数，可选参数必须是实数。可选参数还包括一个默认值，如果用户在输入程序名称时未*输入参数，则使用该值。例如，除非用户*不同的值，否则将为m0分配值0。

注意，在此示例中，n()是样本大小，alpha()是alpha级别，m0()是假设零假设的平均值，ma()是假设备选假设的平均值，sd()是标准偏差。

然后drawnorm使用输入参数从均值为`ma’、标准差为`sd’的正态分布中生成一个观测值为`n’的样本。接下来，ttest测试样本均值等于`m0′的原假设。

可以通过输入return scalar reject =（r（p）<`alpha'）来返回假设检验的结果。回想一下ttest将双边p值存储在标量r(p)中。当r(p)小于`alpha’*的alpha级别时，我们的程序使用标量拒绝返回值1，否则返回0。

现在，可以输入simttest和输入参数来运行模拟。

通过输入return list，可以看到模拟的结果。

如何使用模拟多次运行程序

simttest程序将执行模拟的一次迭代所需的所有操作。接下来，需要一种方法来多次运行simttest并收集结果。下面的代码块显示了如何使用simulate来完成这两项任务。

参数reject=r(reject)告诉simulate将r(reject)中返回的结果保存到名为reject的变量中。选项reps(100)指示模拟运行程序100次。选项seed(12345)设置随机数种子，以便我们的结果可以重现。

冒号之后是simttest以及我们模拟的输入参数。某些模拟需要很长时间才能运行，并且simulate会在结果窗口中显示一个点，以便知道它仍在运行。下面的输出显示了我们的模拟结果。

simulate将结果保存到变量reject，如果零假设的测试被拒绝则包含1，否则为0。

可以使用summarize来计算reject的平均值，该平均值等于100次迭代中拒绝原假设的次数所占的比例。根据输入参数，该比例是对统计功效的估算！

在这个例子中，比例等于0.91，这意味着当检验样本平均值等于75的替代假设时，如果样本平均值等于70，假设标准偏差为15且样本量为100，可以期望91%的功率。

用Power OneMean检查结果

可以使用power onemean检查蒙特卡罗模拟的结果。以下示例包含与我们的模拟相同的输入参数。

power onemean计算的功率为0.9100，这与我们的模拟估算的功率相同。因为模拟是随机的，所有结果并不总是**的。改变随机数种子或迭代次数可以稍微改变估算的功率。但结果应该很接近。

总结

在本文中，介绍了使用蒙特卡罗模拟计算统计功率和样本大小需求所需的工具。下一次，将向您展示如何使用Stata的power命令来运行模拟，以便可以轻松地为一系列输入参数创建表格和图形。

18510103847

词条
词条说明
2023年1-2月Stata软件系列在线培训
Stata是一个完整的、集成的软件包，它提供了所有数据科学所需的数据操作、可视化、统计和自动化报告功能。为了让大家更好的掌握Stata软件，StataCorp LLC为大家安排了一系列免费的在线培训。Ready. Set. Go Stata培训时间：5 January 2023, 1:00 PM CT
TOUGH 系列建模软件介绍
TOUGH3TOUGH3 V1.0使得在PC机上使用单一代码对裂隙多孔介质中非等温多相流动进行预测建模的模拟成为可能。凭借其改进的功能、可用性和输出选项，TOUGH3 是一种更强大、更高效的工具，用于解决与地下流动建模相关的各种科学和实际应用中的困难和计算要求高的问题。注意：可执行文件版本可用于Windows系统。Mac 或 Linux系统，请使用源代码版本并编译代码。 TO
东城Matlab电话
Matlab，作为一款由美国MathWorks公司出品的商业数学软件，其在科学计算、数据分析、可视化以及交互式程序设计等领域广泛应用。科学软件网作为一家引进国研软件并提供软件服务的网站，为国内高校、科研院所和研发型企业事业单位提供了Matlab等科研软件及相关服务。### Matlab的基本概述Matlab全称Matrix Laboratory，是建立在矩阵计算基础上的**编程语言和环境。它集成了
Latent GOLD 潜在类和有限混合软件
Latent GOLD是一个功能强大的潜在类和有限混合程序，具有友好的点击式界面（GUI）。有两个附加选项可用于扩展基本版本的程序。 Advanced / Syntax add-on通过使用包括直观的LG-equations的Syntax命令语言，为高级用户提供更多控制。此附加组件还包含更高级的GUI建模功能，如Latent（Hidden）Markov和Multilevel模型。Choice ad