分类变量的重复测量资料分析

1、分类变量的重复测量分类变量的重复测量资料分析资料分析流行病与卫生统计学教研室流行病与卫生统计学教研室沈毅沈毅2005.3.15 分类变量分类变量（categorical variable）又称为定性变量（qualitative variable），在工作中应用甚广。根据其不同的取值性质，又可分为3种类型： 第一种是名义刻度名义刻度（nominal scale）的分类变量，它是按事物属性分类的变量，如性别、职业等。在统计学上为了计算方便，将这些不同的属性进行数量化处理，如男性赋值为1，女性赋值为2。这种数值只是作为属性的代码，其间并无大小之分。 第二种为有序刻度有序刻度(ordinal scal

2、e)的分类变量，它是根据事物呈现出的程度或水平不同进行赋值。如临床化验结果用符号“-、+、+、+”，文化程度用“文盲、小学、中学、大学、研究生”来划分等级，在进行数量化处理时赋值1、2、3、。这里需要注意的是，1与2之差不一定等于2与3之差。 第三种是区间刻度区间刻度（interval scale)，如人口学统计中的年龄分组，“0-，10-，20-，”就是典型的例子。根据资料的性质，区间跨度有等距的，也有不等距的。 把分类变量作为反应变量进行重复观察的情形在工作中应用较广。在本书第九章第五节中介绍了二分类反应变量的重复测量资料分析方法。 本章将介绍分类反应变量重复测量资料的一般分析方法。主要介

3、绍加权最小二乘法分析方法。第一节一个总体的二分类反应重复测量资料的分析。第一节第一节 二分类反应重复测量二分类反应重复测量资料的分析资料的分析一、资料结构 设有n例受试者，反应变量y，分为 r水平（当r=2时称y为二分类反应变量)，有p个观测时间点，每一受试者可能出现的结果共有rp种。当r=2、p=3时的所有可能结果共有8种。当对y用（0，1）赋值方法时，n例受试者在这8种组合情况下出现的边际合计数为n000至n111。详细分配情况见下面示意图图11 .1。 边际合计数n000至n111是每种可能取值结果组合情况下的观察例数。总例数n为边际合计数之和。如果将上述示意图资料按独立观察资料整理各时

4、间点的频数分布时（表11.1),每一时间点的合计数都为n=n1+=n2+=n3+，总例数变成了3n而不是n。但实际是观察例数为n，观察次数为3n。一般统计学方法是以表11 .1的资料类型为出发点进行分析的。而本章介绍的方法是以图11 .1的边际合计数为出发点进行分析的。二、模型结构 受试者i在时间点j的反应变量为yij ，其取值概率为pj，（j=1,p）。单总体重复测量资料的任务为是要分析这p个概率p1，pp，是否相等。令反应频数的某种函数即反应函数（response function）Ff（m），建立的线性模型为： 式中X为dq维的设计矩阵，d=(r-1)p，为反应函数。B为q维未知参数。

5、由于每一个体在不同时间点的反应变量的取值是相关的，因此，在配合模型时应将这种相关性纳入模型中。对式（11.1)求参数估计值的加权最小二乘解法就是在考虑到这种相关性的基础上建立起来的。模型（11.1)的加权最小二乘解法求参数的估计值 的计算公式为 例11.1一项药物对风湿性关节炎的治疗效果观察。用一种药物对46名患风湿性关节炎的病人进行治疗，共三个疗程。每个疗程结束后评价其治疗效果。有效者记为yij=1，否则yij=0 i=1，,46; j=1,2,3。资料列于表11.2中。表中的T1,T2,T3分别代表三个疗程。本项研究的目的是探讨不同疗程的疗差别。三、配合线性模型的步骤 表11.2为资料的原

6、始记录形式，需要将其整理成边际频数表的格式后再配合模型。计算步骤介绍如下。 1.首先用下列SAS程序计算边际合计数 程序中的subj为受试者号，time1、time2、time3代表3个疗程。计算出的边际频数列于表11.3中。从而得到边际频数向量n及其函数边际频率向量m为： n= (6 16 2 4 2 4 6 6) M=(6/46 0.3478 0.0435 0.0870 0.0435 0.0870 0.1304 0.1304) 用nj及mj表示向量n及m中的元素。 2.给出转换矩阵A 本例的转换矩阵A的结构为 A的作用是将反应频率M转换成反应函数F。 3. 求出反应函数F 它是每一时间点反

7、应变量各个水平的边际概率。但它不包含该反应变量最后一个水平的反应概率。本例有3个时间点（time1 time2 time3），每个时间点的反应变量为r=2水平，故每一时间点只有一个反应函数。共有3个反应函数为： 如f1=0.6086957=0.1304348+0.3478261+0.0434783+0.0869565 以上分析可用SAS的CATMOD过程。程序如下： PROC CATMOD DATA=exmp11_1 order=DATA; RESPONSE MARGINALS; MODEL time1*time2*time3=_RESPONSE_/ONEWAY COV; REPEATED t

8、ime3/_RESPONSE_=time; RUN; 程序中的RESPONSE MARGINALS语句是指定反应函数为在MODEL中规定的反应变量的边际概率。输出结果见表11.5、 表11.6、 表11.7。 从输出的表11.5反应函数和协方差矩阵中看出，第一和第二疗程的有效率都为0.6087（60.87%），第三疗程的有效率为0.3478（34.78%）。从表11.6方差分析看，标记time的有效率分析，第三疗程的效果显著低于前两个疗程。对各疗程疗效一致性的假设检验可用下列公式计算。令对比矩阵C为： 从表11.7加权最小二乘估计值分析可见，有关TIME的两个参数估计值都为正值，反映药物在疗程

9、1、2期的效果优于第3期的效果。在程序中的REPEATED语句的下面加入CONTRAST语句：CONTRAST time1 VS time3_RESPONSE_2 1; 或CONTRAST time1 VS time3 ALL _PARMS 0 2 1;为了检验第一疗程与第三疗程疗效的相等性，所用的检验假设为： Ho：22 3或Ho : 2 2 3 0 ；H1：2 2 3 0 这是因为在参数化时用了关系式4 2 3（ 4代表第三疗程疗效），故有2 4等价于2 2 3 上述SAS程序是用每个病例的记录建立数据集配合模型用的，也可以用频数表资料（表11.3)配合模型。所用SAS程序为：得到输出结果

10、见表11.8。从表11.8对比分析看出，第一期的疗效与第三期疗效间有显著差异（P0.0161）。DATA exmpl 1_1; INPUT time1 time2 time3 count;CARDS; 1 1 1 6 1 1 0 16 1 0 1 2 1 0 0 4 0 1 1 2 0 1 0 4 0 0 1 6 0 0 0 6 ;PROC CATMOD DATA=exmpl 1_1 ORDER=DATA; RESPONSE MARGINALS; WEIGHT count;MODEL time1*time2*time3=_RESPONSE_/ONEWAY COV; REPEATED time3

11、/_RESPONSE_=time;RUN;第二节 多组资料的分析方法 当有多组资料时，统计分析所要回答的问题有三点： 资料是否存在组间差异。 反应是否随时间变化。 是否存在处理组与时间之间的交互效应。分析所用公式以及计算步骤与单总体的情形基本一致。 下面用一个具有二分类反应变量的2组重复测量资料为例加以介绍。 例11.2 对65岁及以上老年人参与某项社会活动的6年追踪调查。样本量为男性662人，女性1311人，分别在当年、第3年和第6年调查他们参与社会活动的情况。反应变量y的赋值规则为： y= 1 经常参加该项社会活动 0 不经常参加该项社会活动 资料列于表11.9中。 本例的分组变量为性别，

12、分组数a2。反应变量y的分类数：r2。重复测量时间点数p3。计算步骤如下。 1列出边际频数矩阵N及相应的边际频率矩阵M 每一组的边际频数占一列，即N(n1，n2 )。与之对应，每一组的边际频率也占频率矩阵中的一列，即M(ml，m2) 。 同时可以建立一个161维的边际频率向量L为： L(0.69 0.07 0.02 0.04 0.03 0.02 0.02 0.12 0.59 0.05 0.02 0.04 0.02 0.03 0.03 0.22） 2. 建立转换矩阵A 由于每一组都是3个反应变量，每一反应变量具有2个水平，所以所建立的转换矩阵A与例11 1的A完全相同。即3计算反应函数F 4分组

13、计算边际频率mg（g1，2）的协方差矩阵V（ml）和V（m2） 5两组反应频率的协方差矩阵V（M）本例的V（M）为1616维矩阵，对角元素为V（mg），非主对角元素为0。 6求F的协方差矩阵V（F） 本例的V（F）为66维矩阵：7定设计矩阵X 例11.2资料的设计矩阵（含交互效应）见表11.11故本例含有交互作用项的设计矩阵为： 8用式（11 .2）计算参数估计值（表11.12）9用对比矩阵C作两组比较的假设检验 C的结构为：利用式（11.6）得到231.38，自由度3，查附表3： 2界值表，有20.05（3）7.81，反映两组间的差异有统计学意义。以上计算过程所用SAS程序如下：DATA e

14、xmpl 1_2; input gender $ year0 year3 year6 count ;CARDS; f 1 1 1 904 f 1 1 0 88 f 1 0 1 25 f 1 0 0 51 f 0 1 1 33 f 0 1 0 22 f 0 0 1 30 f 0 0 0 158 m 1 1 1 391 m 1 1 0 36 m 1 0 1 12 m 1 0 0 143 m 0 1 1 15 m 0 1 0 21 m 0 0 1 18 m 0 0 0 143 ;PROC CATMOD ORDER=DATA; WEIGHT count; RESPONSE marginals; MOD

15、EL year0*year3*year6=gender| _RESPONSE_/PRED=FREQ COV; REPEATED year; RUN;输出结果见表11.13表11.17。 从输出表11 .17方差分析结果看，性别（gender）和年份（year）都有显著的统计学意义（P0.0001）；但性别*年份（gender*year）间交互作用无显著的统计学意义（P0.6476）。 因此，在模型中剔除交互作用项后从新配合下列模型： MODEL year0*year3*year6genderRESPONSE/FREDFREQ COV； 得到输出结果见表11 .18，表11.19。 输出结果表11.19表明，性别间的差异和年份间的差异都具有统计学意义。由于无性别关年份交互作用，所以不同性别受试者对参与某项社会活动的时间趋势是类似的。 为进一步验证这一结论，可对不同时间的效应进行两两比较。提出的检验假设为： H0: year3-year0=year6-year3, H1:year3-year0year6-year3,由于在参数化中year0 , year3的效应是用3和4表示的,year6的效应则为34。故上述假设用参数表示为： 在表11.19加权最小二乘估计值分析中已给出了检