df16_2 <- foreign::read.spss("datasets/例16-02.sav",
to.data.frame = T,
use.value.labels = F,
reencode = "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" ...31 p-for-trend/ p-for-interaction/ per-1-sd R语言实现
本篇主要介绍P for trend、p for interaction、per 1 sd的R语言实现，关于每一项的具体含义，可参考文中给出的链接，或者自己搜索学习。
这几个概念在统计学课本是找不到的，但是在临床研究的SCI论文中经常见到，所以有必要学习它，了解它，实现它！
31.1 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群文件中，需要的加群下载即可。
数据一共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.08400.09279918就是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))
## 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组的，因为第一组是参考，所有组都是和第一组进行比较。
31.2 p for interaction
p for interaction是交互作用的P值，关于其含义可以参考松哥统计的这篇文章：p for interaction是什么
目前计算P for interaction两种方法： 1. 对于数值与等级或二分类，可以直接模型中增加相乘项【如x1×X2】，然后看交互项有无意义。 2. 而对于多项分类【如血型】，产生哑变量后，相乘则会产生多个交互项，此时不能整体判断交互作用是否有意义。我们可以先构建一个无交互作用项的模型，再构建一个有交互作用项的模型。然后采用似然比检验（likelihood ratio test）进行比较有个模型差异，则可以判定交互项整体是否有意义。
31.2.1 方法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)
##
## 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，交互作用项是没有统计学意义的。
31.2.2 方法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，也就是交互作用项没有统计学意义。
31.3 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.0673 7.15 0.00000000285 1.42 1.85
## 2 weight.scaled 1.13 0.0680 1.75 0.0864 0.986 1.29结果给出了P值，OR值以及95%的可信区间。