组轨迹模型(GBTM)的SAS实现与编程技巧

2023-08-31 15:06 作者:统计_碎碎念 0人读过 | 我要投稿

今天，通过一个实例，来展示SAS实现组轨迹模型。

首先，SAS自带的Proc过程步没有实现组轨迹模型的程序，需要从网站上下载“Proc Traj”过程步进行安装，该程序其由Bobby L. Jones团队创作（这里再次跪拜大神，网址奉上：http://www.contrib.andrew.cmu.edu/~bjones/download.htm）。

一、SAS Proc Traj安装过程

查看自己的SAS版本（打开SAS就可以看到）

2.选择适合自己的版本，小编这里选择64位版本的trajv9-64.zip（网站网速可能比较慢，需要的小伙伴也可以后台私信小编）。

3.下载并解压文件夹，其中traj.dll为SAS过程步文件，trajplot.sas，trajplotnew.sas，dropoutplot.sas和trajtest.sas为traj.dll配套使用的宏程序。

4. 找到SAS过程步dll格式文件和宏程序按照路径。我的位置在“C:\Program Files\SAShome\SASFoundation\9.4\core\sasexe”将traj.dll文件复制进该文件夹；将trajplot.sas，trajplotnew.sas，dropoutplot.sas和trajtest.sas宏程序复制进“C:\Program Files\SAShome\SASFoundation\9.4\core\sasmacro”（SAS自带宏都储存在该文件夹下，自动调用）。

5.打开SAS尝试输入程序语句，发现Proc traj字样和Model字样都有了SAS过程步的颜色，即为按照成功。

二、组轨迹模型实战

某研究者对学生的心理状况进行了评估，共进行了7学年的随访。T1-T7为评估学年，O1- O7为心理评估得分，得分范围为0~10分，得分越高心理状况越好。

研究目的：在7个评估学年中，学生的心理发展轨迹。

SAS数据格式整理如下：

1.数据导入及转化操作

proc import datafile="&path_raw\实例数据.xlsx"

out=dataset dbms=xlsx replace;

getnames=yes;

run;

*宽形数据转化为长形数据以便绘图;

data plot1;

set dataset ;

array vv {*} O1-O7;

array tt {*} t1-t7;

do i = 1 to dim(vv);

idd=id;

time=tt(i);

value=vv(i);

output;

end;

keep time value idd;

run;

2.窥探数据分布，描述数据分布情况、缺失情况、发展情况。

proc means data =dataset n nmiss mean stddev

stderr min p25 median p75 max qrange clm;

var O1-O7 T1-T7;

run;

proc sgplot data = plot1;

vbox value / group=time ;

run;

ods html dpi=300;

proc sgplot data = plot1;

vline time / response=value group=idd ;

run;

3. SAS的语法结构

（1）组轨迹模型可适用的数据分布类型：PROC TRAJ可以拟合删失正态分布(CNORM)、β分布(BETA)、零膨胀Poisson分布(ZIP)和伯努利分布(LOGIT)。本研究因变量为心理评估得分，估拟合CNORM模型。

（2）PROC TRAJ基本语法结构：

PROC TRAJ /option调用组轨迹模型。

OUT= 输出组分配和成员概率数据集（可根据此数据集计算平均后验概率）。

OUTPLOT=轨迹图数据（根据此数据集可以进行个性化绘图）。

OUTSTAT=轨迹组的参数估计数据集。

OUTEST=参数和协方差矩阵估计数据集。

ITDETAIL 显示用于监控模型拟合进度的最小化迭代。

ci95m 估计95%置信区间。

ID：一般指研究对象唯一个案，这里可以为OUT=数据集中，想保留的变量。

VAR：因变量。

INDEP：测量因变量（VAR）时的自变量（如年龄、时间）。

MODEL：因变量分布（BETA、CNORM、ZIP、LOGIT），例如MODEL CNORM。

MIN：CNORM分布特有，指定删失值的最小值。

MAX：CNORM分布特有，指定删失值的最大值。

NGROUPS：拟合轨迹组数量。

ORDER：轨迹组的形态参数，每个轨迹组的多项式 (0=intercept, 1=linear, 2=quadratic, 3=cubic)。

（3）建立模型

考虑需要拟合多个组轨迹模型，为了模型拟合的方便性和便捷性，建议有基础的小伙伴采用宏程序处理。

小编构建了宏名称为%traj_M()的宏程序，包含data,id,avepp,EST,plot,stat,VAR,INDEP,NGROUPS,ORDER,title 11个宏参数。程序中%TRAJPLOTnew()为Proc traj 配套使用的绘图宏程序。

%macro traj_M(data,id,avepp,EST,plot,stat,VAR,INDEP,NGROUPS,ORDER,title);

PROC TRAJ DATA=&data. OUTPLOT=&plot. OUT=&avepp. OUTSTAT=&stat.

OUTEST=&EST. ITDETAIL ci95m;

ID &id.;

VAR &VAR.;

INDEP &INDEP.;

MODEL CNORM;

min 0;

MAX 10; NGROUPS &NGROUPS.; ORDER &ORDER.;

RUN;

proc sort data = &avepp.; by group;run;

ods html dpi=300;

%TRAJPLOTnew(&plot.,&stat.,"ZBMI",'Cnorm Model','BMI','Time(month)');

quit;

%mend;

（4）构建2~6组轨迹，轨迹形态分别为linear, quadratic, cubic的模型。也就是我们需要构建15个组轨迹模型。小编以拟合3组轨迹模型分别是cubic的示例：

%let data =dataset;/*调用的数据集名称*/

%let prefix=test_;/*输出数据集的前缀*/

%let VAR=O1-O7;/*因变量*/

%let INDEP=T1-T7;/*时点*/

%let id=ID;/*ID*/

%let postfix=3_333;/*输出数据集的后缀*/

%let n_grp=3;/*3组轨迹模型*/

/*1.拟合3-333轨迹模型*/

%traj_M(data=&data.,id=&id.,avepp=&prefix.avepp_&postfix.,EST=&prefix.est_&postfix.,plot=&prefix.plot_&postfix.,stat=&prefix.stat_&postfix.,VAR=&var.,INDEP=&INDEP.,NGROUPS=&n_grp.,ORDER=3 3 3,title="&postfix.")

%traj_M直接输出的结果和图形，我们需要利用这些结果和图形计算模型的评价指标，%Traj_Assess是小编写的计算模型拟合指标的宏程序，可以直接计算所有拟合指标：

数据集：&prefix.avepp_&postfix.结果