原文链接：http://tecdat.cn/?p=10809

原文出处：拓端数据部落公众号

简介

本文档用于比较六个不同统计软件程序（SAS，Stata，HLM，R，SPSS和Mplus）的两级多层（也称分层或层次）线性模型的过程和输出。 

下面介绍的六个模型都是两级分层模型的变体，也称为多级模型，这是混合模型的特殊情况。此比较仅对完全嵌套的数据有效（不适用于交叉或其他设计的数据，可以使用混合模型进行分析）。尽管HLM软件的网站声明可以用于交叉设计，但这尚未得到确认。下面的SAS，Stata，R，SPSS和Mplus中使用的过程是其多层次或混合模型过程的一部分，并且可以扩展为非嵌套数据。

但是出于比较的目的，我们将仅研究完全嵌套的数据集。除了HLM（完全由GUI运行）以外，所有程序的下面都包含用于每个模型的代码/语法。我们提供了HLM和SPSS的屏幕截图。此外，每个模型均以分层格式和混合格式指定。尽管模型的这两个表达式是等效的，但一些研究领域更倾向于可视化层次结构，因为它更容易看到层次之间的分离，而另一些研究领域则更喜欢混合格式，在其中容易区分固定效果和随机效果。

模型注意事项将预测变量添加到本文档讨论的六个模型中时，我们选择以均值居中为中心，这意味着我们从每个受试者的得分中减去了该变量的总体均值。

正如Enders＆Tofighi（2007）所详细讨论的那样，以总体平均值为中心，而不是以组平均值（每个组的平均值均以该组中受试者的得分为准）为中心，并不适合所有模型。 。使用哪种居中方法的选择应由所询问的具体研究问题决定。另一个考虑因素是这些程序使用的估计方法来产生参数估计，即最大似然（ML）或受限最大似然（REML）。每种都有自己的优点和缺点。ML更适合不平衡的数据，但是会产生偏差的结果。REML是无偏的，但是在将两个嵌套模型与似然比检验进行比较时，不能使用REML。

两种方法将产生相同的固定效应估计，但它们对随机效应的估计却有所不同（Albright＆Marinova，2010）。正如我们将在下面讨论的模型中看到的那样，这两种方法产生的结果非常相似，并且不会极大地影响随机因素的p值。但是，重要的是要意识到，方法的选择会影响随机因素的估计，标准误差和p值，并且可能会影响宣布随机因素是否重要的​​决策。SAS，HLM，R和SPSS默认使用REML，而Stata和Mplus使用ML。在本文档中的Stata示例中，我们告诉Stata使用REML以便将输出与其他四个程序进行比较。

类内相关系数

我们还报告了每种模型的类内相关系数（ICC）ρ。ICC是结果变量中方差的比例，由分层模型的分组结构解释。它是根据组级别误差方差与总误差方差之比来计算的：

其中，是2级残差的方差，是1级残差的方差。换句话说，与总的无法解释的方差（方差之内和之间）相比，ICC报告了模型中任何可归因于分组变量的预测变量无法解释的变化量。

示例数据集

流行的数据集由来自不同班级的学生组成，并且由于每个学生都属于一个唯一的班级，因此它是一个嵌套设计。因变量是“流行”，它是一个自评的流行度，范围为0-10。预测指标包括学生级别的性别（二分法）和Extrav（连续的自我评价的外向得分），以及班级的Texp（多年的老师经验， 是连续的）。

仅截距模型（无条件模型）

无条件混合模型规范类似于单因素方差分析，其总体均值和类效应。但是，我们将其视为随机效应（均值为零的正态分布变量），而不是像方差分析中那样的固定因子效应。因此，我们将估计值解释为每个类别的平均数在总体平均人气得分附近的方差。

估算值是每个班级的“大众”平均值的平均值，而不是研究中所有学生的平均值。如果数据完全平衡（即每个班级的学生人数相同），则无条件模型的结果将与方差分析程序的结果相同。

SAS结果

 需要“ covtest”选项来报告方差分量估计的标准误差。另外， 需要指定非结构化协方差矩阵类型，这是HLM和R默认情况下使用的类型，我们在这里使用它进行比较。SAS的输出等于Hox的书表2.1中的结果。我们可以得出结论，各类别之间的平均人气得分为5.078，并且各类别之间的差异（1.221）比不同类别之间的差异（0.702）多。当我们为该模型计算ICC时，将对此进行进一步讨论。

Stata结果

 Stata的xtmixed命令需要因变量，后跟“ ||”  指定固定变量和随机变量之间的分隔。我们必须包括方差选项以查看输出中方差分量的估计值，以及reml选项以使用受限的最大似然估计。还要注意，Stata不会输出随机分量估计的p值，但是可以通过置信区间中是否包含零来确定有效值。这些结果与SAS的结果完全匹配

HLM结果

HLM报告方差组件的标准偏差，而不是标准误差。同样，对于随机效应，他 仅报告截距的卡方统计量和p值。这些结果与其他程序的结果相同。

R结果

R报告方差成分（例如HLM）的标准偏差，而lme4软件包报告固定效应的t统计量。 

 SPSS结果

 屏幕截图：

 

 需要在“随机”窗口中指定非结构化协方差类型。这些结果与其他程序和本文得出的结果相同。请注意，像SAS和Mplus一样，SPSS报告方差分量的标准误差，而HLM和R报告标准差。我们无法得出结论，哪个更适合报告，但是差异不会影响这些参数的p值。

 因为这是一个无条件模型，所以我们不需要指定任何WITHIN或BETWEEN变量。下面列出了在MODEL语句中列出变量的标准。在以下各节中，我们将看到前三个示例：

1.％WITHIN％– 1级固定因子（非随机斜率）2.具有潜在斜率变量的％WITHIN％– 1级随机因子3.％BETWEEN％– 2级固定因子4.在任一个陈述–在学生水平上测得的变量，但具有1级和2级方差估计 。

上表显示了Mplus输出底部的“模型结果”部分的结果。Mplus确实会报告每个估计的p值，并且所有估计都与其他程序的p值匹配，但随机截距的方差估计相差约0.007。这种差异是由于Mplus使用ML估计这一事实造成的。尽管存在这种差异，但我们看不到任何变量的重要性发生变化。

 汇总

总体而言，这六个程序对于仅截取模型产生了非常相似的结果（唯一的差异发生在随机效应的Mplus估计中）。唯一的区别是他们如何报告随机方差估计的精度。此模型的ICC等于：

这告诉我们，“流行”课程总变化的大约三分之一可以由每个学生所在的班级解释。

具有一个固定的Level-1因子的随机截距（非随机斜率）模型

该模型增加了一个学生级别的固定因子Extrav，即自我报告的外向得分。混合模型看起来像是基于带有协变量Extrav的类的ANCOVA，但请记住，我们仍然认为这是随机效应，而不是固定效应。因此， 估算值与ANCOVA程序所得出的估算值不同。 

 在此数据的实际应用中，Extrav应该具有固定的效果而不是随机的效果是没有意义的，因为学生外向性的水平应随班级而变化。但是，出于比较这四个程序的目的，我们仍然希望调查一个具有一个学生级别固定因子的案例。

SAS结果

现在，我们对Extrav的固定效果进行了估算。学生报告的外向得分每增加一个单位，他们的受欢迎度得分就会增加0.486。这些结果等于使用REML的其他程序的结果。

Stata结果 

当我们向Stata中的模型添加预测变量时，我们添加了cov（un）选项，指定了非结构化协方差矩阵。我们将 Extraversion变量放在“ ||”之前 表示它是一个固定因子（具有非随机斜率）。这些结果与其他程序的结果相同。

 HLM结果

 

R结果

 SPSS结果

Mplus结果

 现在，我们在VARIABLE语句的WITHIN选项中包括居中的Extrav变量。对于内部MODEL规范，我们必须使用“ ON”选项，以告知Mplus Extrav是固定的1级因子。  可以看到由于使用ML估计而不是REML，许多估计和估计的标准误差（以及t统计量）存在细微差异。由于方差的估计值与其他程序不同，因此Mplus报告的ICC与下面报告的有所不同。

汇总

对于此模型，前五个程序的结果完全相同，而Mplus的估算值相差很小。此模型的ICC大于无条件模型的ICC（正如预期的那样，因为我们通过添加固定因子来控制某些学生水平的变化）：

 

 使用一个学生水平的固定因子，“流行”总变化的几乎一半可以由该学生的班级和学生水平的固定因子“外向”解释

一个一级因子的随机截距和斜率模型

 

该模型包含Extrav的随机斜率，这意味着我们允许回归方程的斜率随类而变化。该模型比以前的模型更适合于所使用的变量，因为可以直观地假设外向因类而异。

SAS结果

随机Extrav斜率的估计值很重要（p值为0.003），因此我们可以说学生的外向性得分随班级的不同而变化。这些结果与其他程序的结果完全匹配，除了固定效果的t统计量存在一些细微差异。

 Stata结果

 HLM结果

SPSS结果

 Mplus结果

 这次我们在WITHIN语句中包括一个潜在的斜率变量，以将Extrav指定为随机因子，该变量告诉Mplus不要在数据集中寻找“ randoms1”，因为没有观察到它。 可以将此变量的输出解释为Extrav的随机斜率分量。我们必须这样做，因为Mplus是为结构方程模型设计的，其多级模型功能是其潜在潜伏分析程序的改编。

汇总

总体而言，前五个程序对该模型产生相同的结果，而Mplus再次由于ML估计而相差很小。此模型的ICC为：

 通过将Extrav的影响从固定变为随机，由于我们正在考虑在学生级别上更多随机变化，因此ICC会略有增加。

两个1级因子的随机斜率模型

 

 对于此模型，我们包括第二个学生级别的变量Sex，该变量也具有随机斜率。这意味着我们既要考虑学生的性别，又要考虑他们的外向得分，并且允许这两个因素的斜率随班级而变化。

SAS结果

 

 在此输出中，我们可以看到性别确实对学生自我报告的知名度有重大影响（p值<0.001）。对Sex的固定估计意味着，在Extrav不变的情况下，女学生（Sex = 1）的普遍得分比男学生（基线组，Sex = 0）高1.244。

SAS不喜欢在该模型中，Sex的估计方差非常接近零，因此没有输出标准误差或p值。因为非常接近于零，所以我们可以得出结论，性别不会因类别而显着变化。 

Stata结果

 Stata在运行该模型时引用了一个错误：标准误差计算失败，这意味着未计算随机效应的标准误差。我们发现通过删除cov（un）选项，不会出现此错误。但是，该输出中的所有估计均与其他程序不同，因此我们选择使用非结构化协方差矩阵规范报告输出。我们不确定这是否是在Stata中运行此类模型的常见问题，但重要的是要意识到它会发生。

HLM结果

这些估计大致等于其他计划的结果，但随机性别影响的估计除外。由于这种影响非常接近于零，因此程序不会报告完全相同的值，但是所有结果都表明该值远非重要。

 R结果

SPSS结果

 Mplus结果

 这次，我们在WITHIN语句中包括了两个潜在的斜率变量，以将Extrav和Sex指定为随机因子。我们可以将“ randoms1”的输出解释为Extrav的估计，将“ randoms2”的输出解释为Sex的估计。

Mplus针对此模型的输出所得出的估计值与先前模型中的其他程序相距甚远。我们看到，由于模型必须估计更多随机参数，因此估计程序（ML与REML）之间的差异变得更加明显。但是，Mplus同意其他程序的观点，即“性别”的随机方差部分以外的所有估计值都非常重要。

 汇总

对于方差非常接近零的随机效应，六个程序以不同的方式处理估计值。SAS和Stata无法报告随机效应的标准误差或p值，而其他变量的估计值和标准误差均具有相当大的差异。Mplus结果也显示出比以前的模型更大的差异。此模型的ICC为：

同样，当我们在模型中添加另一个学生级别的效果（包括随机斜率）时，ICC略有增加。 

一个2级因子和两个随机1级因子（无交互）

 这是我们看到的第一个具有2级（班级）变量的模型：教师的多年经验（Texp），也是以均值为中心的。正如Enders和Tofighi（2007）指出的那样，级别2变量的唯一居中选项是均值居中。  无法对均值中心Texp进行分组，因为它已经在班级水平上进行了度量，这意味着“分组均值”将等于原始值。

在分层格式中， 可以看到它具有固定的斜率系数，并且对于每个类j都是唯一的。该模型在教师的经验和学生水平的变量之间没有任何相互作用。如果我们有理由相信Texp不会缓和Sex和Extrav对Popular的影响，那么我们将使用此模型，这意味着我们的学生水平变量的斜率是相同的，无论学生是否有新教师或新教师。一位拥有多年经验的人。

SAS结果

 

 现在，我们在固定效果表中看到了Texp，估计值为0.089，p值很大。这意味着，在使学生的性别和性取向得分保持不变的情况下，每增加一年的教师经验，该学生的热门得分就会增加0.089。

同样，我们看到SAS无法处理随机性别效应的很小变化。因此，没有报告标准误差，z统计量或p值。

Stata结果

 与以前的模型一样，我们收到一个错误，告诉我们Stata无法计算方差分量的标准误差。但是，这些估计值与其他程序的估计值大致相同。

HLM结果

 

 

 这些估计值与其他程序的结果略有不同  。

R结果

SPSS结果 

 

 Mplus结果

对于Level-2因子，我们在BETWEEN语句中包括Ctexp。我们再次看到这些估计数与其他五个计划的输出有微小出入

 汇总

使用REML的五个程序的输出实质上是相等的，仅相差几千个单位。与以前的模型一样，最大的差异出现在随机性别效应的方差估计中，因为它非常接近零。

请注意，此模型的ICC比以前的模型有所降低（= 0.542）：

 请记住，ICC是衡量 所在的班级可以解释多少无法解释的变化的方法。通过添加班级级别的预测变量，我们可以解释不同班级中较大比例的变化。因此，与没有任何2级预测变量的模型相比，该模型的随机截距存在较少的变异，因此ICC也较低。

具有相互作用的一个2级因子和两个随机1级因子

 

 这是我们在班级变量Texp与学生级变量Sex和Extrav之间进行跨级交互的唯一模型。例如，如果我们想找出具有更多经验的教师是否比新教师对学生的外向性或性别与他们自我报告的知名度之间的关系有不同的影响，则可以使用此模型。换句话说，教师的经历是否适度了性格外向或性别对受欢迎程度的影响？

 可以看到，在分层格式中，Texp在三个方程式的每个方程式内都有一个斜率系数。这与混合模型中的交互项有关，即通过外向的教师体验和按性别的教师体验。

SAS结果

 在固定效果表中，有两个交互作用项，其中一个（）远不重要，p值> 0.5。 

在随机方差分量表中，我们看到外向随机斜率的估计值和性别随机斜率的估计值与零没有显着差异。这意味着没有证据表明这两个因素实际上在该模型中因类别而异。 

 Stata结果

Stata无法自动识别变量之间的交互项，因此我们必须为两个跨级别的交互手动创建变量（请参见上面的代码中的gen语句）。当我们使用带有非结构化协方差矩阵选项的xtmixed命令运行时，Stata给出了一个错误，指出Hessian不是负半定性，一致性错误，并且没有产生任何输出。 

HLM结果

这些估计值大致等于其他程序的结果。

R结果 

SPSS结果 

 对于SPSS 19而言，此模型实在太多了。对于具有非结构化协方差矩阵的更复杂的模型，其他程序可能会运行更有效的算法，因此优于SPSS。

 Mplus结果

 现在，我们在BETWEEN模型部分中包含两个ON语句，以指示与教师经验的跨层次交互。同样，我们发现与其他输出之间的细微差异，但Mplus同意Texp和Sex之间的固定交互作用不显着，而Extrav和Sex的随机组成也并不重要。

 汇总

加上两个跨层交互项，Stata和SPSS无法使用非结构化协方差选项运行模型。这并不是说不应该将它们用于这种类型的分析，但是在向具有非结构化协方差矩阵的模型中添加更复杂的参数时，应谨慎使用。

与以前的模型一样，SAS，HLM和R的结果相对接近相等，而Mplus的估计略有不同。另外，ICC与模型5几乎完全相同，这意味着交互作用项不会改变按类别说明的差异比例：

总结

进行比较的目的是调查来自六个不同统计软件程序的嵌套两级层次模型的过程和结果可能存在的差异。总体而言，我们发现SAS，Stata（带有reml选项），HLM，R和SPSS产生的实际估计值之间没有太大差异。Mplus使用另一种估算方法ML，这导致其估算值与其他估算值有所不同。另外，重要的是要注意以下几点：

1.对于方差估计非常接近零的随机效应，SAS无法产生标准误差或p值。其他三个程序在估计这些参数方面的差异与其他效果相比更大。

2. Stata和SPSS无法处理最复杂的模型，该模型包含两个跨级别的交互项。建议使用其他程序来分析复杂模型并指定非结构化协方差矩阵。

此外，我们研究了每种模型中类内相关系数的值。通过添加1级预测因子，ICC有所增加。但是，当我们添加2级预测变量时，ICC会大大降低，甚至比无条件模型更低。这是由于在类级别添加了预测变量时，无法解释的Level-2变异（随机截距项）减少了。

尽管本文档可以用作为嵌套数据集运行各种两级分层模型的指南，但我们强烈建议读者仅在适合回答您的特定研究问题时使用这些模型。在确定固定因素和随机因素之间，以及对于中心平均值为1的总体平均值或组平均值时，必须谨慎使用。 

参考文献

Enders, Craig K. and Tofighi, Davood (2007). “Centering Predictor Variables in Cross-Sectional Multilevel Models: A New Look at an Old Issue.” Psychological Methods, vol. 12, pg. 121-138.

Hox, Joop J. (2010). Multilevel Analysis (2nd ed.). New York: Routledge

---

最受欢迎的见解

1.用SPSS估计HLM层次线性模型模型

2.R语言线性判别分析（LDA），二次判别分析（QDA）和正则判别分析（RDA）

3.基于R语言的lmer混合线性回归模型

4.R语言Gibbs抽样的贝叶斯简单线性回归仿真分析

5.在r语言中使用GAM（广义相加模型）进行电力负荷时间序列分析

6.使用SAS，Stata，HLM，R，SPSS和Mplus的分层线性模型HLM

7.R语言中的岭回归、套索回归、主成分回归：线性模型选择和正则化

8.R语言用线性回归模型预测空气质量臭氧数据

9.R语言分层线性模型案例

一、练习

• 列出价格大于6000的国产汽车的价格

  sysuse auto, clear
  list price if price > 6000 & foreign == 0
  

• 给出1978年维修记录少于3次或产地为国外的汽车价格和重量的描述性统计信息

  sum price weight if (rep78 < 3 | foreign == 1)
  

• 本数据中有多少辆国产汽车价格大于6000？

  count if price > 6000 & foreign == 0
  

• 列出 price wei len mpg turn foreign 变量的均值，标准差，中位数，最大值，最小值（对于这种要求，将一系列统计量以一张表格的形式展现出来更合适）

  tabstat price wei len mpg turn foreign, stat(mean sd median max min)
  

• 按国产和非国产为标准分类对price wei 进行描述性统计

  by foreign : sum price wei
  

• 将数据按照价格进行升序和降序排列

  gsort price
  gsort -price
  

• 删除1978年维修记录缺失的数据

  drop if rep78 == .
  

• 将字符型变量make转换成数值型变量cenmake

  encode make, gen(cenmake)
  

• 将price 转化成字符型变量

  tostring price, replace
  

• 将price 转化成数值型变量

  destring price, replace
  

• 生成一个虚拟变量，其中价格大于6000为1，否则表示为0

  gen vari = 1 if price > 6000
  replace vari = 0 if price <= 6000
  //这样比较麻烦，可以写成如下语句
  gen dum_price = price > 6000
  

• 给出price wei len mpg 相关系数矩阵

  pwcorr price wei len mpg
  

• 画出price wei len mpg 的相关系数矩阵散点图（graph是最需要掌握的画图命令）

  graph matrix price wei len mpg
  

• 画出price频率直方图

  histogram price, frequency
  

• 画出price和weight的散点图

  scatter price weight
  

• 给price 和weight散点图加上price和wei的拟合线

  twoway (scatter price weight) (lfit price weight)
  

---

二、线性回归模型

• 做一个简单的OLS回归，其中被解释变量为price,解释变量为weight foreign

1. 初步探索

cd D:\Stata_Projects  //进入存储数据的工作目录
use "food.dta", clear  //使用这个数据集

* 得到数据之后，一般先浏览，检查数据，并看变量的描述性统计信息
* 画出food_exp与income的散点图，并保存
twoway (scatter food_exp income)
* 保存图片, 如果有了就替代
graph save food1, replace 
* 一步到位
twoway (scatter food_exp income), saving(food1, replace) 

* 更丰富的画图功能, 使用help graph就可以查看到
twoway (scatter food_exp income) ///  /* 基本语句，画什么的散点图 */
	   (lfit food_exp income), ///    /*拟合线*/
       ylabel(0(100)600)         ///  /* Y轴的刻度从0到600，间隔为100 */           
	   xlabel(0(5)35)            ///  /* X轴刻度0-35，间隔5 */
	   title(食物消费数据)   /// /* 图片名字 */
	   xtitle(每周收入)   ///
       ytitle(每周家庭食物花费)
* 还有更多的画图内容可以学习，图片的保存等也可以直接点击图片操作

2. 简单OLS回归(最小二乘法)

* food_ex = coef * income + ... 对income进行回归
regress food_exp income

* 计算y的拟合值，这里的y是food_exp'; xb是计算线性拟合值
predict yhat, xb  //calculate linear prediction
* 计算残差项的拟合值
predict ehat, residuals 
browse

* 做预测，预测其他的收入情况在该模型中的食物支出
* 例如，预测家庭每周收入2000$的家庭，每周食物支出
set obs 41 //本数据有40个观测值，我们多加一个观测值，但是没有定义具体观测值内容，它目前是缺失
replace income = 20 in 41 //将收入第41个观测值替换成20百
predict yhat1, xb
list yhat1 in 41

* 保存dta文件
* 直接点save菜单，或者输入下面命令
save "newfood.dta", replace

---

三、二次式模型

$y={\beta }_1+{\beta }_2\ast x^2+e$

* 二次式模型
use "br.dta",clear  
sum //sum, detail

desc

gen sqft2 =sqft^2 //生成变量
regress price sqft2  //二次模型
predict priceq, xb
br price sqft sqft2 priceq

* 画出拟合曲线
twoway (scatter price sqft) 	   ///     	/* 基本的点 */
      (lfit price sqft,            ///		/* 拟合线 */
       sort lwidth(medthick))	   /*设定线的厚度 */
graph save br_quad1, replace       //保存图片 

twoway (scatter price sqft) 	  ///     	/* 基本的点 */
       (line price sqft,          ///		/* 折线 */
       sort lwidth(medthick))	  /*设定线的厚度 */
graph save br_quad2, replace      //保存图片	
	
twoway (scatter price sqft) 	  ///     	/* 基本的点 */
       (line priceq sqft,         ///		/* 折线 */
       sort lwidth(medthick))     /*设定线的厚度 */
graph save br_quad3, replace      //保存图片	
graph combine "br_quad1" "br_quad2" "br_quad3"
*	比较三个图片区别	

---

四、对数线性模型

$lny={\beta }_1+{\beta }_{2}x+e$

histogram price

graph save price, replace 
	
gen lnprice=ln(price)        //price 存在一些值特别大,尾巴很长,为了把数据平滑一点,取对数
histogram lnprice,percent    //更接近正态分布

reg lnprice sqft    //对数线性模型回归
predict lpricef, xb //得到对数形式的拟合值

* 画出拟合曲线
generate pricef = exp(lpricef)
twoway (scatter price sqft)  
       (line pricef sqft, sort lwidth(medthick))
graph save br_loglin, replace

---

五、指示变量回归

首先进一步学习创建虚拟变量/factor variables/指示变量：

• 方法1. 设定是否满足某个条件来生成，条件为真，取值 $1$，否则取 $0$：

  gen bh = (bedrooms>4) // 卧室数量大于四为大房子，取值为1 
  

  

• 方法2. $tabulate$ 创建指示变量。让每个可能的 $baths$ 取值都设定一个独立的虚拟变量：

  tabulate region, gen(reg)  // Create indicator variables for region, called reg1, reg2, ...
  tabulate baths, gen(bath)  // 为每个baths取值计数，gen(bath)就会产生5个指示变量bath1-5
  

  

---

六、练习

use "utown.dta", clear
describe  //拿到数据的规定动作，了解数据内容和描述性统计信息
sum  

* 按照房子是否位于大学附近进行分类画图	
histogram price if utown==0, width(12) start(130) percent  ///
          xtitle(House prices ($1000) in Golden Oaks)      ///
	      xlabel(130(24)350) legend(off)
graph save utown_0, replace

histogram price if utown==1, width(12) start(130) percent  ///
          xtitle(House prices ($1000) in University Town)  ///
	      xlabel(130(24)350) legend(off)
graph save utown_1, replace
graph combine "utown_0" "utown_1", col(1) iscale(1) 
*把两张图放在一张图上，col(#)分成几列显示，iscale()文本和标记的大小
graph save combined ,replace

histogram每一个语句的意思需要好好理解，例如：

• histogram price if utown==0：意思是做房屋不位于大学附近的房子的价格的直方图；
• width(12)是柱子宽度为 $12$ 个刻度；
• start(130)价格从 $130$ 开始画（price最小值为 $134$，这样x轴不会从 $0$ 开始）
• percent表示纵轴为占比；
• xlabel表示 $x$ 轴刻度从 $130$ 到 $350$，单位刻度(即间隔刻度)为 $25$；
• legend(off)表示不显示图例

label define utownlabel 0 "Golden Oaks" 1 "University Town"  //定义utown值的label
label value utown utownlabel //将utown的label也按照值的label定义

• label variable varname 给变量定标签
• label define labelname 给变量的值定义标签
• label value varlist(一列变量名) labelname 将取值的标签定义给变量

* by(utown, cols(1))  按照utown分组画图	
histogram price, by(utown, cols(2))  			 
          start(130) percent                	 
          xtitle(House prices ($1000))     		 
          xlabel(130(24)350) legend(off)
graph save combined2, replace

• 检验独立社区房价的描述性统计信息
  sum price if utown == 0
sum price if utown == 1
  
• 复习，用 bysort 分组统计
  bysort utown:sum price
  
• 另一种写法
  by utown,sort:sum price
  

$by$ and $bysort$ are really the same command; $bysort$ is just $by$ with the $sort$ option.