p for trend/ p for interaction/ per 1 sd 的R语言实现
本篇主要介绍per 1 sd的R语言实现,关于每一项的具体含义,可参考文中给出的链接,或者自己搜索学习。这几个概念在统计学课本是找不到的,但是在临床研究的SCI论文中经常见到,所以有必要学习它,了解它,实现它!
本篇主要介绍P for trend
、p for interaction
、per 1 sd
的R语言实现,关于每一项的具体含义,可参考文中给出的链接,或者自己搜索学习。
这几个概念在统计学课本是找不到的,但是在临床研究的SCI论文中经常见到,所以有必要学习它,了解它,实现它!
P for trend
P for trend
是线性趋势检验的P值,用于反映自变量和因变量是否存在线性趋势关系。线性趋势检验,之前介绍过Cochran Armitage
检验,不过是针对分类变量的。
今天要介绍的P for trend
主要是针对连续型变量的。
关于p for trend
具体含义和数值型变量分箱的方法,大家可以参考医咖会的文章:p for trend是个啥
把连续性变量转换为分类变量(在R里转变为因子),设置哑变量,进行回归分析,即可得到OR值及95%的可信区间;把转换好的分类变量当做数值型,进行回归分析,即可得到P for trend
使用之前逻辑回归的例子演示,来自孙振球版医学统计学第4版,电子版和配套数据均放在QQ群文件中,需要的加群下载即可。(加我微信问群号的一律拉黑!)
df16_2 <- foreign::read.spss("../000统计学/例16-02.sav",
to.data.frame = T,
use.value.labels = F,
reencode = "utf-8")
## re-encoding from utf-8
str(df16_2)
## 'data.frame': 54 obs. of 11 variables:
## $ .... : num 1 2 3 4 5 6 7 8 9 10 ...
## $ x1 : num 3 2 2 2 3 3 2 3 2 1 ...
## $ x2 : num 1 0 1 0 0 0 0 0 0 0 ...
## $ x3 : num 0 1 0 0 0 1 1 1 0 0 ...
## $ x4 : num 1 1 1 1 1 1 0 1 0 1 ...
## $ x5 : num 0 0 0 0 0 0 0 1 0 0 ...
## $ x6 : num 0 0 0 0 1 0 0 0 0 0 ...
## $ x7 : num 1 1 1 1 1 2 1 1 1 1 ...
## $ x8 : num 1 0 0 0 1 1 0 0 1 0 ...
## $ y : num 0 0 0 0 0 0 0 0 0 0 ...
## $ PGR_1: num 1 0 0 0 1 1 0 0 0 0 ...
## - attr(*, "variable.labels")= Named chr [1:11] "" "" "" "" ...
## ..- attr(*, "names")= chr [1:11] "...." "x1" "x2" "x3" ...
数据一共11列,第1列是编号,第2-9列是自变量,第10列是因变量。
具体说明:
-
x1:年龄,小于45岁是1,45-55是2,55-65是3,65以上是4;
-
x2:高血压病史,1代表有,0代表无;
-
x3:高血压家族史,1代表有,0代表无;
-
x4:吸烟,1代表吸烟,0代表不吸烟;
-
x5:高血脂病史,1代表有,0代表无;
-
x6:动物脂肪摄入,0表示低,1表示高
-
x7:BMI,小于24是1,24-26是2,大于26是3;
-
x8:A型性格,1代表是,0代表否;
-
y:是否是冠心病,1代表是,0代表否
这里的x1~y
虽然是数值型,但并不是真的代表数字大小,只是为了方便标识,
年龄x1
应该是数值型的,但是为了方便解释逻辑回归的意义,我们对它进行了分箱处理,也就是把它转换为了分类变量。数值型变量进行分箱,是回归分析中计算p for trend的第一步
此时x1
是数值型,我们直接进行逻辑回归,得到的P值就是 p for trend
f <- glm(y ~ x1 + x2,
data = df16_2,
family = binomial())
broom::tidy(f)
## # A tibble: 3 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) -2.22 1.03 -2.15 0.0313
## 2 x1 0.712 0.423 1.68 0.0928
## 3 x2 1.08 0.625 1.73 0.0840
0.0928
就是x1
的p for trend
,而且还是校正了x2
这个变量之后的p for trend
,是不是很简单?
此时如果我们把x1
变成因子型,那在进行回归分析时会自动进行哑变量编码,就可以得到几个组的OR值和95%的可信区间,关于R语言中分类变量进行回归分析时常用的一些编码方法,强烈你看一下这篇推文:R语言分类变量进行回归分析的编码方案。
# 变为因子型
df16_2$x1.f <- factor(df16_2$x1)
# 把因子放入自变量
f <- glm(y ~ x1.f + x2,
data = df16_2,
family = binomial())
broom::tidy(f,conf.int=T,exponentiate=T)
## # A tibble: 5 × 7
## term estimate std.error statistic p.value conf.low conf.high
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 0.200 1.10 -1.47 0.142 0.0104 1.24
## 2 x1.f2 2.32 1.19 0.704 0.481 0.289 49.3
## 3 x1.f3 4.48 1.26 1.19 0.233 0.485 102.
## 4 x1.f4 9.42 1.63 1.38 0.169 0.508 438.
## 5 x2 2.94 0.639 1.69 0.0918 0.854 10.7
这样就得到了x1.f
中4/3/2分别和1进行比较的OR值和95%的可信区间。当然你写函数提取也行:
# OR值
exp(coef(f))
## (Intercept) x1.f2 x1.f3 x1.f4 x2
## 0.200000 2.319343 4.476753 9.415697 2.936212
# OR值的95%的可信区间
exp(confint(f))
## Waiting for profiling to be done...
## 2.5 % 97.5 %
## (Intercept) 0.01043892 1.240092
## x1.f2 0.28932733 49.284803
## x1.f3 0.48497992 102.335996
## x1.f4 0.50766137 437.541812
## x2 0.85353009 10.723068
这样就得到了每个组的OR值和95%的可信区间,可以看到没有第1组的,因为第一组是参考,所有组都是和第一组进行比较。
p for interaction
p for interaction
是交互作用的P值,关于其含义可以参考松哥统计的这篇文章:p for interaction是什么
目前计算P for interaction两种方法:
-
对于数值与等级或二分类,可以直接模型中增加相乘项【如x1×X2】,然后看交互项有无意义。
-
而对于多项分类【如血型】,产生哑变量后,相乘则会产生多个交互项,此时不能整体判断交互作用是否有意义。我们可以先构建一个无交互作用项的模型,再构建一个有交互作用项的模型。然后采用似然比检验(likelihood ratio test)进行比较有个模型差异,则可以判定交互项整体是否有意义。
方法1
假如探索年龄(x1)和BMI(x7)之间对因变量y有没有交互作用,我们首先新建一列相乘列,然后进行回归分析。
# 新建1列
df16_2$x17 <- df16_2$x1 * df16_2$x7
str(df16_2)
## 'data.frame': 54 obs. of 13 variables:
## $ .... : num 1 2 3 4 5 6 7 8 9 10 ...
## $ x1 : num 3 2 2 2 3 3 2 3 2 1 ...
## $ x2 : num 1 0 1 0 0 0 0 0 0 0 ...
## $ x3 : num 0 1 0 0 0 1 1 1 0 0 ...
## $ x4 : num 1 1 1 1 1 1 0 1 0 1 ...
## $ x5 : num 0 0 0 0 0 0 0 1 0 0 ...
## $ x6 : num 0 0 0 0 1 0 0 0 0 0 ...
## $ x7 : num 1 1 1 1 1 2 1 1 1 1 ...
## $ x8 : num 1 0 0 0 1 1 0 0 1 0 ...
## $ y : num 0 0 0 0 0 0 0 0 0 0 ...
## $ PGR_1: num 1 0 0 0 1 1 0 0 0 0 ...
## $ x1.f : Factor w/ 4 levels "1","2","3","4": 3 2 2 2 3 3 2 3 2 1 ...
## $ x17 : num 3 2 2 2 3 6 2 3 2 1 ...
## - attr(*, "variable.labels")= Named chr [1:11] "" "" "" "" ...
## ..- attr(*, "names")= chr [1:11] "...." "x1" "x2" "x3" ...
# 进行逻辑回归
f <- glm(y ~ x1 + x7 + x17,
family = binomial(),
data = df16_2
)
summary(f)
##
## Call:
## glm(formula = y ~ x1 + x7 + x17, family = binomial(), data = df16_2)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.7967 -0.9092 -0.6507 1.1651 1.7708
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.3953 2.6477 -0.149 0.881
## x1 -0.4867 1.1142 -0.437 0.662
## x7 -1.1509 1.7095 -0.673 0.501
## x17 0.9249 0.7489 1.235 0.217
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 74.786 on 53 degrees of freedom
## Residual deviance: 62.508 on 50 degrees of freedom
## AIC: 70.508
##
## Number of Fisher Scoring iterations: 5
结果中显示x17
的P值(p for interaction)是:0.217,交互作用项是没有统计学意义的。
方法2
# 先构建一个没有交互项的逻辑回归模型
f1 <- glm(y ~ x1 + x7,
family = binomial(),
data = df16_2)
# 再构建一个有交互作用的逻辑回归模型
f2 <- glm(y ~ x1 + x7 + x17,
family = binomial(),
data = df16_2)
# 似然比检验
lmtest::lrtest(f1,f2)
## Likelihood ratio test
##
## Model 1: y ~ x1 + x7
## Model 2: y ~ x1 + x7 + x17
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 3 -32.216
## 2 4 -31.254 1 1.9238 0.1654
结果显示P(p for interaction)=0.1654,也就是交互作用项没有统计学意义。
per 1 sd
关于什么是per 1 sd
,可以参考松哥统计的这篇文章:per 1 sd
Per 1 sd的实现,其实就是把原始数据进行标准化,另存为一个新的变量X,新变量X因为是被标准化后的数据,因此其均数和标准差为0和1。然后让x进入模型进行分析。请问大家此时x每增加1个单位,效应量增加的风险为HR。因为标准差为1,此时x增加1个单位,就是Per 1 sd。1=Per 1 sd。就是自变量每增加1个标准差。
为了方便演示,我们新建一列数据weight
,然后进行标准化,再进行逻辑回归。
# 新建一列weight
df16_2$weight <- rnorm(54, 70,11)
# 进行标准化
df16_2$weight.scaled <- scale(df16_2$weight)
# 进行逻辑回归
f <- glm(y ~ weight.scaled, data = df16_2)
broom::tidy(f,conf.int=T,exponentiate=T)
## # A tibble: 2 × 7
## term estimate std.error statistic p.value conf.low conf.high
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 1.62 0.0690 6.98 0.00000000529 1.41 1.85
## 2 weight.scaled 1.05 0.0696 0.702 0.486 0.916 1.20
结果给出了P值,OR值以及95%的可信区间。
![Logo](https://devpress.csdnimg.cn/79de2bf0b7994defa4242ef90d5513fa.jpg)
开放原子开发者工作坊旨在鼓励更多人参与开源活动,与志同道合的开发者们相互交流开发经验、分享开发心得、获取前沿技术趋势。工作坊有多种形式的开发者活动,如meetup、训练营等,主打技术交流,干货满满,真诚地邀请各位开发者共同参与!
更多推荐
所有评论(0)