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

功率和样本大小的计算是科学研究计划的重要组成部分。可以使用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命令来运行模拟，以便可以轻松地为一系列输入参数创建表格和图形。