如果把 Z 个应变量分开, 逐个进行回归分析, 就是建立 Z 个回归方程 :Y i = X, 这与逐个一样, 没有同时考虑 X 对 Z 个 Y 的作用现在做多元回归是把这 Z 个回归方程合并成一个方程如果数据本身不是对同一个指标多次测量, 而是多个应变量, 只要这些应变量反映的性质相同, 变化

Size: px

Start display at page:

Download "如果把 Z 个应变量分开, 逐个进行回归分析, 就是建立 Z 个回归方程 :Y i = X, 这与逐个一样, 没有同时考虑 X 对 Z 个 Y 的作用现在做多元回归是把这 Z 个回归方程合并成一个方程如果数据本身不是对同一个指标多次测量, 而是多个应变量, 只要这些应变量反映的性质相同, 变化"

泺付谢
7 years ago
Views:

1 多元分析与重复测量数据分析 ( 一 ) 多元 T 2 检验广义估计方程 (GEE) 多元回归 (1) 易尔统计软件三种方法具体操作 (2) 三种分析方法统计检验效率比较本文通过模拟比较不同统计方法, 进行多元分析的统计检验效率, 得出 GEE 多元回归分析效率最高, 并介绍了如何使用易尔统计软件很方便地进行 GEE 多元回归分析临床研究数据常有以下两种 : 1 对 N 个研究对象观测了 Z 项指标 ( 结局变量 ), 如比较某降血脂药用药前后总胆固醇低密度脂蛋白甘油三酯的变化, 结局变量有 3 个同时对这 Z 个结局变量进行分析被称为多元分析 (multivariate analysis) 2 对 N 个研究对象某一指标 ( 结局变量 ) 观测了 Z 次这 Z 次可能是 Z 个观测员或是在 Z 种情况下所做的观测, 如服某降血压药后小时分别测量的血压值对这 Z 次重复测量结果的分析就是重复测量数据的分析 Z 次测量也可以称为 Z 个结局变量上两种情况可以总称为多元 ( 多应变量 ) 分析一分析方法 : 1) 单个逐个比较每个应变量一般的统计方法是将 Z 个结局变量分开来分析, 即做 Z 个单变量的分析如分别比较用某降血脂药前后总胆固醇高密度脂蛋白甘油三酯的变化, 做 3 个 T 检验或回归分析 ; 分别比较用某降血压药后小时血压变化, 做 4 个 T 检验或回归分析这种将多个因变量拆分成单元进行分析的缺点是 :(1) 增加了 I 类错误, 如 α 定在 0.05, 一次检验有 5% 可能拒绝真的无效假设, 三次检验就有约 15% 的可能拒绝一个真的无效假设 (2) 将多因变量拆分进行单元分析的统计检验效率低多元分析同时对 Z 个变量进行检验, 不存在多次检验需要校正 p 值的问题, 且同时考虑到对 Z 个结局变量的变化, 检验效率高 2) Hotelling T2 检验同时比较 Z 个应变量单个变量两样本均数的比较用, 多个变量一起比较可用 Hotelling T 2 检验其原理与基本相似这个方法只适用于比较两组满足多元正态分布的 Z 个连续性变量 3) 广义估计方程 (GEE) 多元回归

2 如果把 Z 个应变量分开, 逐个进行回归分析, 就是建立 Z 个回归方程 :Y i = X, 这与逐个一样, 没有同时考虑 X 对 Z 个 Y 的作用现在做多元回归是把这 Z 个回归方程合并成一个方程如果数据本身不是对同一个指标多次测量, 而是多个应变量, 只要这些应变量反映的性质相同, 变化方向一致, 如血脂中的总胆固醇低密度脂蛋白甘油三酯三项指标, 肺功能的一秒肺活量与最大肺活量, 呼吸道症状的咳嗽咳痰气喘气短等, 也可以看成是重复测量数据也就是把 {Y1,Y2,, Yi 看成是对 Y 的重复测量首先把数据进行转换, 原来一个研究对象一条记录, 有 Z 个应变量 (Y), 转换成一个研究对象 Z 条记录, 每条一个应变量 (Y), 在每条记录后增加一个指示变量 (I) 表示该应变量值的来源然后建立回归模型 :Y= X + I 其中 I 是 Z 个指示变量这个模型假定 X 对不同的应变量 (Y) 的作用是一致的, 不同应变量的截距可以不同 ( 通过 I 表现 ), 但变化大小一样因为一般线性回归与 Logistic 回归都要求观察记录相互独立, 即研究个体间无相关关系但此时同一个人有多条记录, 这些记录间就有相关关系, 不符合独立性要求因此需要用广义估计方程 (Generalized Estimate Equation) 简称 GEE 这个方法适合于连续性变量, 也适用于两分类变量 ( 用 logistic 回归 ), 及其它分布类型 ( 用相应的联系函数 ) 如果原来的 Z 个应变量单位不同, 在单个回归方程中回归系数差异很大, 可通过调整相应的应变量的单位, 或进行标准化等方法, 尽量使得所研究的危险因素对每个应变量的回归系数接近二三种方法在易尔统计软件上的具体操作 : 1 ) 单个 : (1) 在基本统计菜单下调用对均数的比较与检验, 再选两样本比较的, 如下图所示 :

(2) 在两样本比较的窗口, 通过拖拉 ( 与双击 ) 方法, 选择分析变量, 同样方法选择分组变量, 如下图所示, 点击运行即可软件自动对每个变量进行, 输出结果摘录如下 : 两样本均数比较的按 GROUP 分层的 INFO 的样本量均数标准差 N Mean Sd 0 28 11.89286 2.819762 1 12 9.75000 2.

3 (2) 在两样本比较的窗口, 通过拖拉 ( 与双击 ) 方法, 选择分析变量, 同样方法选择分组变量, 如下图所示, 点击运行即可软件自动对每个变量进行, 输出结果摘录如下 : 两样本均数比较的按 GROUP 分层的 INFO 的样本量均数标准差 N Mean Sd Two Sample t-test data: INFO by GROUP t = , df = 38, p-value = alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: sample estimates:

mean in group 0 mean in group 1 11.89286 9.

4 mean in group 0 mean in group )Hotelling T2 检验同时比较 Z 个应变量检验 (1) 在基本统计菜单下调用多元分析, 再选两样本多元比较的 Hotelling T 平方 (2) 在两样本多元比较的 T 平方检验窗口, 通过拖拉 ( 与双击 ) 方法, 选择分析变量, 同样方法选择分组变量, 如下图所示, 点击运行即可

5 输出结果如下 : Hotelling T 平方检验 : INFO SIMIL ARITH PICT 分组变量 : GROUP 分组值 : 0, 1 N1 N2 T 平方 F 自由度 1 自由度 2 P 值 )GEE 多元回归分析 : (1) 在基本统计菜单下调用多元分析, 再选广义估计方程多应变量回归 (2) 在广义估计方程多元回归分析窗口, 通过拖拉 ( 与双击 ) 方法, 选择应变量, 同样方法选择危险因素 ( 同上述的分组变量 ), 如下图所示, 点击运行即可

易尔统计软件自动转换数据, 自动对每个应变量的单个回归系数, 与总的回归系数进行估计, 并评价单个回归系数之间是否有显著差异, 用一个总的回归系数表达危险因素的作用是否合适给出结果如下例所示 : 使用广义估计方程作多应变量回归分析 : 单变量 ( 结局 ) 回归分析 Call: geeglm(formula = tmp.

6 易尔统计软件自动转换数据, 自动对每个应变量的单个回归系数, 与总的回归系数进行估计, 并评价单个回归系数之间是否有显著差异, 用一个总的回归系数表达危险因素的作用是否合适给出结果如下例所示 : 使用广义估计方程作多应变量回归分析 : 单变量 ( 结局 ) 回归分析 Call: geeglm(formula = tmp.y ~ +GROUP.1.INFO + GROUP.1.SIMIL + GROUP.1.ARITH + GROUP.1.PICT + factor(tmp.yidx), family = gaussian(link = "identity"), data = WD, id = tmp.id, corstr = "independence")

7 Coefficients: Estimate Std.err Wald Pr(> W ) (Intercept) < 2e-16 *** GROUP.1.INFO * GROUP.1.SIMIL ** GROUP.1.ARITH * GROUP.1.PICT * factor(tmp.yidx) e-08 *** factor(tmp.yidx) factor(tmp.yidx) e-12 *** --- Signif. codes: 0 '***' '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Estimated Scale Parameters: Estimate Std.err (Intercept) Correlation: Structure = independencenumber of clusters: 40 Maximum cluster size: 4 分别为 4 个应变量单独的回归系数使用 GEE 作多变量 ( 结局 ) 回归分析 : Call: geeglm(formula = tmp.y ~ GROUP.1 + factor(tmp.yidx), family = gaussian(link = "identity"), data = WD, id = tmp.id, corstr = "independence") Coefficients: Estimate Std.err Wald Pr(> W ) (Intercept) < 2e-16 *** GROUP *** factor(tmp.yidx) e-11 *** factor(tmp.yidx) factor(tmp.yidx) e-15 *** --- Signif. codes: 0 '***' '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Estimated Scale Parameters: Estimate Std.err (Intercept) Correlation: Structure = independencenumber of clusters: 40 Maximum cluster size: 4 综合 4 个应变量的回归系数与 p 值检验危险因素 :GROUP 对不同结局变量的回归系数之间是否有显著性差异

8 Analysis of 'Wald statistic' Table Model 1 tmp.y ~ +GROUP.1.INFO + GROUP.1.SIMIL + GROUP.1.ARITH + GROUP.1.PICT + factor(tmp.yidx) Model 2 tmp.y ~ GROUP.1 + factor(tmp.yidx) Df X2 P(> Chi ) 适合用多变量回归系数来表示 GROUP 危险因素的总的效应通过比较两个模型似然值判断 4 个单个的回归系数是否有显著差异, 用一个总的回归系数是否合适使用易尔统计软件做 GEE 多元回归分析的优点 : 不用进行数据转换, 软件自动编程进行转换不用逐个运行单个回归模型, 检验回归系数差异是否显著, 软件自动编程进行比较判断 R 软件中的 geeglm 程序包要求数据是按分析对象排序的, 否则即使在调用过程中定义了 id= 与 corstr=, 这个定义没起作用, 结果并没有调整同一个 id 的多条记录的相关性, 得出的结论有误很多人在使用 R 的 geeglm 过程中很容易忽略了这个问题, 得出的结果是不正确的易尔软件自动编程, 不会忽略这个问题三三种方法统计检验效率比较如果你的数据适合三种方法, 哪种更有效率呢? 统计检验效率是指能检测出差异的能力临床研究中常常因为经费人力时间等原因, 样本量不够大, 希望能用小的样本检出较小的差异, 这就是怎么能提高统计检验效率的问题下面通过模拟, 对三种方法进行比较, 得出使用广义估计方程 (GEE) 对多应变量数据转换后进行回归分析, 统计检验效率最高而且随着变量数增加, 效率增加 Hotelling T 2 检验效率与应变量数关系不大 Hotelling T 2 检验与 GEE 多元回归两种方法随应变量之间的相关性增加时, 统计检验效率都会下降, 但 Hotelling T 2 检验降低更明显具体模拟方法 : (1) 假定两组数据, 各组 500 例 ; (2) 分别假定两个变量 X1 X2, 三个变量 X1 X2 X3, 四个变量 X1 X2 X3 X4; (3) 假定第一组 X1 X2 X3 X4 的均数分别是 , 第二组 X1 X2 X3 X4 的均数分别是

9 (4) 假定 X1 X2 X3 X4 的方差均是 10, 分别假定协方差是 , 即分别假定 X1 X2 X3 X4 相互间的相关系数分别是 ; 用随机抽样方法, 对上述每种情况分别模拟 1000 次, 对每次模拟出来的数据进行单个变量的 Hotelling T 2 检验 GEE 回归分析, 如得出的 p 值小于 0.05 表示检验成功检出差异统计 1000 个模拟数据中, 有多少个 p 值小于 0.05 结果如下 : 两个应变量模拟结果变量间相关系数 Hotelling GEE 多元回单变量 1 单变量 2 T 2 归三个应变量模拟结果变量间相关系数 Hotelling T 2 GEE 多元回归单变量 1 单变量 2 单变量四个应变量模拟结果变量间相关系数 Hotelling T 2 GEE 多元回归单变量 1 单变量 2 单变量 3 单变量附 R 模拟程序 : ## Power Simulation for comparing HotellingT2 test, MVGee and single t test library(geepack) library(mass) library(icsnp)

10 ## 2 variables simupower2<-function(mean, mdiff, sigma, r) { for (i in (1:1000)) { d1<-mvrnorm(500,mean, sigma) d2<-mvrnorm(500,mean+mdiff, sigma) p.t2<-hotellingst2(d1,d2)$p.value dd1<rbind(cbind(id=(1:500),grp=0,v=1,y=d1[,1]),cbind(id=(1:500),grp=0,v=2,y=d1[,2])) dd2<rbind(cbind(id=(501:1000),grp=1,v=1,y=d2[,1]),cbind(id=(501:1000),grp=1,v=2,y=d2[,2])) dd<-data.frame(rbind(dd1,dd2)) dd<-dd[with(dd, order(id)), ] gee<geeglm(y~grp+factor(v),id=id,corstr="independence",data=dd,family=gaussian(link="ident ity")) p.gee<-summary(gee)$coefficients[2,4] ppi<-c(2,r,p.t2,p.gee,t.test(d1[,1],d2[,1])$p.value,t.test(d1[,2],d2[,2])$p.value) if (i==1) {pp<-ppi else {pp<-rbind(pp,ppi) return(pp) ## 3 variables simupower3<-function(mean, mdiff, sigma, r) { for (i in (1:1000)) { d1<-mvrnorm(500,mean, sigma) d2<-mvrnorm(500,mean+mdiff, sigma) p.t2<-hotellingst2(d1,d2)$p.value dd1<rbind(cbind(id=(1:500),grp=0,v=1,y=d1[,1]),cbind(id=(1:500),grp=0,v=2,y=d1[,2]),cbind( id=(1:500),grp=0,v=3,y=d1[,3])) dd2<rbind(cbind(id=(501:1000),grp=1,v=1,y=d2[,1]),cbind(id=(501:1000),grp=1,v=2,y=d2[,2]), cbind(id=(501:1000),grp=1,v=3,y=d2[,3])) dd<-data.frame(rbind(dd1,dd2)) dd<-dd[with(dd, order(id)), ] gee<geeglm(y~grp+factor(v),id=id,corstr="independence",data=dd,family=gaussian(link="ident ity")) p.gee<-summary(gee)$coefficients[2,4] ppi<-c(3,r,p.t2,p.gee,t.test(d1[,1],d2[,1])$p.value, t.test(d1[,2],d2[,2])$p.value, t.test(d1[,3],d2[,3])$p.value) if (i==1) {pp<-ppi else {pp<-rbind(pp,ppi) return(pp) ## 4 variables simupower4<-function(mean, mdiff, sigma, r) { for (i in (1:1000)) { d1<-mvrnorm(500,mean, sigma) d2<-mvrnorm(500,mean+mdiff, sigma) p.t2<-hotellingst2(d1,d2)$p.value dd1<rbind(cbind(id=(1:500),grp=0,v=1,y=d1[,1]),cbind(id=(1:500),grp=0,v=2,y=d1[,2]),cbind( id=(1:500),grp=0,v=3,y=d1[,3]),cbind(id=(1:500),grp=0,v=4,y=d1[,4])) dd2<rbind(cbind(id=(501:1000),grp=1,v=1,y=d2[,1]),cbind(id=(501:1000),grp=1,v=2,y=d2[,2]), cbind(id=(501:1000),grp=1,v=3,y=d2[,3]),cbind(id=(501:1000),grp=1,v=4,y=d2[,4]))

11 dd<-data.frame(rbind(dd1,dd2)) dd<-dd[with(dd, order(id)), ] gee<geeglm(y~grp+factor(v),id=id,corstr="independence",data=dd,family=gaussian(link="ident ity")) p.gee<-summary(gee)$coefficients[2,4] ppi<c(4,r,p.t2,p.gee,t.test(d1[,1],d2[,1])$p.value,t.test(d1[,2],d2[,2])$p.value,t.test(d1 [,3],d2[,3])$p.value,t.test(d1[,4],d2[,4])$p.value) if (i==1) {pp<-ppi else {pp<-rbind(pp,ppi) return(pp) pw2.2 <- simupower2(c(2,3),0.2,matrix(c(10,2,2,10),2,2),0.2) pw2.4 <- simupower2(c(2,3),0.2,matrix(c(10,4,4,10),2,2),0.4) pw2.6 <- simupower2(c(2,3),0.2,matrix(c(10,6,6,10),2,2),0.6) pw2.8 <- simupower2(c(2,3),0.2,matrix(c(10,8,8,10),2,2),0.8) pw3.2 <- simupower3(c(2,3,4),0.2,matrix(c(10,2,2,2,10,2,2,2,10),3,3),0.2) pw3.4 <- simupower3(c(2,3,4),0.2,matrix(c(10,4,4,4,10,4,4,4,10),3,3),0.4) pw3.6 <- simupower3(c(2,3,4),0.2,matrix(c(10,6,6,6,10,6,6,6,10),3,3),0.6) pw3.8 <- simupower3(c(2,3,4),0.2,matrix(c(10,8,8,8,10,8,8,8,10),3,3),0.8) pw4.2 <- simupower4(c(2,3,4,5),0.2,matrix(c(10,2,2,2,2,10,2,2,2,2,10,2,2,2,2,10),4,4),0.2) pw4.4 <- simupower4(c(2,3,4,5),0.2,matrix(c(10,4,4,4,4,10,4,4,4,4,10,4,4,4,4,10),4,4),0.4) pw4.6 <- simupower4(c(2,3,4,5),0.2,matrix(c(10,6,6,6,6,10,6,6,6,6,10,6,6,6,6,10),4,4),0.6) pw4.8 <- simupower4(c(2,3,4,5),0.2,matrix(c(10,8,8,8,8,10,8,8,8,8,10,8,8,8,8,10),4,4),0.8) pw2<-rbind(pw2.2,pw2.4,pw2.6,pw2.8) pw3<-rbind(pw3.2,pw3.4,pw3.6,pw3.8) pw4<-rbind(pw4.2,pw4.4,pw4.6,pw4.8) pp2 <- ifelse(pw2[,-c(1,2)]<0.05,1,0) ps2 <- apply(pp2,2,x<-function(x) tapply(x,factor(pw2[,2]),sum)) colnames(ps2) <- c("t2","gee","v1","v2") pp3 <- ifelse(pw3[,-c(1,2)]<0.05,1,0) ps3 <- apply(pp3,2,x<-function(x) tapply(x,factor(pw3[,2]),sum)) colnames(ps3) <- c("t2","gee","v1","v2","v3") pp4 <- ifelse(pw4[,-c(1,2)]<0.05,1,0) ps4 <- apply(pp4,2,x<-function(x) tapply(x,factor(pw4[,2]),sum)) colnames(ps4) <- c("t2","gee","v1","v2","v3","v4") ps2 ps3 ps4 作者联系地址 :changzhong_chen@dfci.harvard.edu

1

015 年 ( 第四届 ) 全国大学生统计建模大赛参赛论文论文名称 : 打车软件影响下的西安市出租车运营市场研究和统计分析参赛学校 : 西安理工大学参赛队员 : 刘二嫚余菲王赵汉指导教师 : 王金霞肖燕婷提交日期 :015 年 6 月 9 日 1