多元统计分析成果-多元统计分析与ＳＡＳ实现

本章所用数据来自2016年全国流动人口动态监测调查，我们只使用该数据库中7个项目实施点的2335名流动人口数据来做模型示范，本节将数据集命名为exe11来演示倾向评分匹配。

以下示例说明了如何对干预组和对照组的观测值进行贪婪匹配，以便可以在随后的结果分析中将匹配的观察值用于估计干预效果。流动人口健康素养对其接受健康教育的影响，该数据集包含了流动人口的相关信息，接受健康教育种类，健康素养，性别，年龄，户口，婚姻状况，是否有工作，就业身份。

其变量赋值具体如下：

因变量：

health＿education：接受健康教育的种类（连续变量）

自变量：

HL：健康素养（1＝高健康素养，0＝低健康素养）

协变量：

gender：性别（1＝男性，0＝女性）

age：年龄（连续变量）

hukou：户口（1＝农业，0＝非农业）

marrige：婚姻状况（1＝在婚；0＝非在婚）

job：是否有工作（1＝有，0＝无）

employment＿status：就业身份（1＝雇员，0＝非雇员）

SAS程序：

SAS程序解释：

proc psmatch是应用倾向性评分的程序。此示例说明了使用proc psmatch程序将干预组中个体的观察结果与对照组中具有相似倾向得分的观察结果进行匹配。proc psmatch过程仅匹配那些倾向得分在region＝选项中指定的支持区域内的观测值。在这里，region＝treated选项要求仅将倾向得分位于已处理观测值定义的区域中的那些观测值用于匹配。默认情况下，该区域扩展为倾向得分对数的通用标准偏差的合并估计的0.25倍。匹配的观察值保存在输出数据集中，加上结果变量后，可用于提供干预效果的无偏估计。

psmodel语句指定逻辑回归模型，该模型为每个个体创建倾向分数，这是该个体高健康素养的概率。健康素养HL是二元干预指标变量，并且treated＝“1”将1标识为干预组。将性别、年龄、户口、婚姻状况、工作、就业身份变量纳入模型，因为它们被认为与分配有关。class语句指定分类变量。

match语句请求匹配并指定匹配条件。distance＝lps选项（默认设置）要求将倾向得分的对数用于计算观察对之间的差异。method＝greedy（k＝1）选项请求贪婪最近邻匹配，其中一个控制单元与治疗组中的每个单元匹配，这样，在该处理单元的所有可用对中，对内差异最小。caliper＝0.5选项指定匹配的卡尺要求。仅当两组中成对的单位的倾向得分的对数的差值小于或等于标准偏差的合并估计值的0.5倍时，才对单位进行匹配。

assess语句会生成一张表格并绘制图表，以汇总干预组和对照组之间指定变量的差异。可以使用这些结果来评估匹配如何在这些变量的分布中达到平衡。根据lps和var＝选项的要求，变量是倾向得分的对数以及协变量性别、年龄等。weight＝none选项禁止显示加权匹配观测值的差异。当proc psmatch将一个控制单元与每个处理单元匹配时，它将为所有匹配的处理单元和控制单元分配1的权重，因此加权匹配观测值和匹配观测值的结果相同。

output语句中的out（obs＝match）＝选项创建一个名为OutEx4的输出数据集，其中包含匹配的观察值。默认情况下，此数据集包括变量＿PS＿（提供倾向得分）和变量＿MATCHWGT＿（提供匹配的观察权重）。lps＝＿lps选项添加一个名为＿lps的变量，该变量提供倾向得分的对数，而matchid＝＿matchid选项添加一个名为＿matchid的变量，该变量标识匹配的观测值集。

需要说明的是该程序只能在SAS9.4 M5及以上版本运行。

SAS结果：

SAS部分结果输出如下：

图11-1　倾向值程序数据信息

图11-2　倾向值评分信息

图11-3　匹配信息

图11-4　干预组和对照组标准化均值差异

图11-5　标准化均值差异

图11-6　匹配集前5对数据信息

SAS结果解释：

图11-1显示了干预组和对照组的样本量分别为738和1256，支持域中倾向得分的上限和下限分别为0.5320和0.2257，支持域中干预组和对照组得样本量分别为738和1254。

图11-2显示了干预组和对照组的基本统计信息。这些统计数据是针对所有观测值，支持区域中的观测值以及匹配的观测值计算的。因为指定了region＝treated，所以这三组统计信息在干预组是相同的。

图11-3显示了匹配方法为1∶1贪婪最近邻匹配，卡尺为0.1160，匹配组数738，干预组和对照组中匹配观察值的数量都为738，所有匹配倾向得分对数的总绝对差为0.5861。

图11-4和图11-5显示了干预组和对照组之间变量的标准平均差，这些变量是针对所有观察值，支持域中的观察值以及匹配的观察值计算的。在匹配的观察中，所有变量的标准化均值差异相较于所有观察值均显著降低，且小于建议的上限0.25（由图11-5可知），匹配观测值中的处理后控制方差比为0.9741和1.0053，也在0.5到2的建议范围内（由图11-4可知）。

在对变量平衡性不满意的情况下，可以执行以下一项或多项操作来改善：选择另一组变量以适应倾向得分模型；修改匹配条件，或者选择另一种匹配方法，直到匹配结果可以很好地平衡变量，然后将匹配的观测值保存在输出数据集中，以用于后续结果分析。

output语句创建一个输出数据集OutEx4，其中包含匹配的观察值，并列出了前五个匹配集中的观测值。默认情况下，输出数据集包括变量＿PS＿（倾向值得分）和变量＿MATCHWGT＿（匹配的观察权重）。每个处理单元的权重为1。由于在match语句的method＝optimal选项中指定了k＝1，因此每个处理组都匹配一个对照组，因此每个匹配控制单元的权重也为1。matchid＝＿MatchID选项创建一个名为的变量＿MatchID（其标识匹配的观测值）。

SAS程序：

运行匹配后的结果分析的SAS程序如下：

proc ttest data＝exe11；

class HL；

var health＿education；

run；

SAS结果：

图11-7　t检验结果（匹配前）

SAS结果解释：

图11-7是对原始数据（匹配前）进行t检验的结果，p＜0.001显示高健康素养组和低健康素养组接受健康教育的种类具有显著性差异。尽管t检验显示了显著的影响，但这种影响可能与流动人口的年龄、性别、户籍、婚姻状况、工作、就业身份有关。

以下是控制了相关影响因素的SAS程序。

SAS程序：

proc glm data＝exe11；

class HL(ref＝ˈ0ˈ)gender hukou marrige jobs employment＿status；

model health＿education＝HL age gender hukou marrige jobs employment＿status/solution；

run；

SAS结果：(www.xing528.com)

图11-8　回归分析结果（匹配前）

SAS结果解释：

在控制了这些变量的情况下，结果仍然显示出流动人口健康素养对于其接受健康教育的促进作用，具体来说，高健康素养组接受健康教育的种类比低健康素养组多0.3870。

通过执行使用贪婪最近邻匹配的倾向得分匹配，以确保足够的协变量平衡。由于贪婪最近邻匹配的语句在上文已经给出，此处不再重复，提示可不使用assess语句。

以下语句使用t检验从匹配的观察值估计健康素养的效应。

SAS程序：

proc ttest data＝OutEx4a；

class HL；

var health＿education；

run；

SAS结果：

图11-9　t检验结果（匹配后）

SAS结果解释：

图11-9中，p＜0.001显示在匹配后高健康素养组和低健康素养组接受健康教育的种类具有显著性差异。

接下来进行对匹配观察值的回归分析，SAS程序如下。

SAS程序：

proc glm data＝OutEx4a；

class HL(ref＝ˈ0ˈ)gender hukou marrige jobs employment＿status；

model health＿education＝HL age gender hukou marrige jobs employment＿status/solution；

run；

SAS结果：

图11-10　回归分析结果（匹配后）

SAS结果解释：

图11-10中，控制年龄、性别、户籍、婚姻状况、工作、就业身份的影响，结果表明在匹配后流动人口健康素养对于其接受健康教育具有显著的促进作用，具体来说，高健康素养组接受健康教育的种类比低健康素养组多0.236。

对比匹配前和匹配后的回归结果可以看出，基于匹配观察的结果分析得出的结论与基于原始数据的结果分析得出的结论一致，即提高流动人口健康素养显著促进其接受更多类型的健康教育，虽然匹配前健康素养的效应（0.3870）略大于匹配后健康素养的效应（0.2362）。

接下来进行敏感性分析，以下语句计算health＿education每个匹配集中干预组与控制组之间的差异Diff。

SAS程序：

SAS结果：

图11-11　匹配集中健康教育差异

SAS结果解释：

图11-11列出了前十个匹配集中health＿education干预组和控制组的差异。

然后用以下语句执行符号秩检验。

SAS程序：

SAS结果：

图11-12　符号秩检验结果

SAS结果解释：

图11-12结果表明，干预组的接受健康教育的种类在0.0001水平上显著提升。

基于匹配观测值的敏感性分析部分中描述的方法，用以下语句计算符号秩统计量。

SAS程序：

SAS结果：

图11-13　符号秩统计量

SAS结果解释：

图11-13显示了符号秩检验的统计信息，使用此统计信息，以下语句将进行符号秩等级计算并显示p值，这些值对应于范围从1到2的Γ值。

SAS程序：

SAS结果：

图11-14　敏感性分析结果

SAS结果解释：

图11-14显示，在Γ＝1.65时，p值为0.0586，犯I类错误概率大于0.05。而一般来说，Γ＞2的情况下，可以认为研究的敏感性较低，可以说研究已免除隐藏性偏差的影响，因此该案例得出的结论对于未观察到的混杂因素所隐藏的偏见并不可靠，还需要进一步改善匹配的共变量。