1 GLM模型简介
参数 family 规定回归模型的类型
链接函数 | Y的分布 | 回归类型 |
binomial(link = "logit") | 二元离散 二项分布 | 逻辑回归 |
gaussian(link = "identity") | 连续 正态分布 | 线性回归 |
Gamma(link = "inverse") | 非负 右偏误差分布 | Gamma回归 |
inverse.gaussian(link = "1/mu^2") | 连续的正值 逆高斯分布 | 逆高斯回归 |
poisson(link = "log") | 非负次数 离散型 泊松分布 | 泊松回归 |
quasi(link = "identity", variance = "constant") | 连续性 非负性 常数方差 | 平方根误差回归 |
quasibinomial(link = "logit") | 二元离散 非负性 | 逻辑回归 |
quasipoisson(link = "log") | 非负次数 离散型 | 广义线性模型 |
2 拟合glm模型
我们拟合常规的logistic回归
# 1 转化为因子
dat1$sex <- factor(dat1$sex)
# 2 拟合模型
mod <- glm(insomnia~sex,family = binomial,data=dat1)
# 3 总结模型
summary(mod)
可以知道,拟合的logistic回归模型,相对于Sex赋值为1的人群,Sex赋值为2的人群,OR=EXP(0.68023),P值<0.001
3 抓取回归的结果
接下来我们对其进行组合
为了方便接下来的工作,我们加载dplyr包
# 加载包
library("dplyr")
mod_dat <- mod %>% summary %>% coef %>% as.data.frame()
此时的结果长这个样子,已经有
Estimate:β回归系数
Std. Error:S.E (标准误)
z value: Z值
Pr(>|z|):P值
4 OR的计算
很明显,这个结果还不够,下面我们需要计算OR值和95%CI
mod_dat <- mod %>% summary %>% coef %>% as.data.frame() %>%
mutate(Beta=sprintf("%0.2f",Estimate)) %>%
mutate(S.E=sprintf("%0.2f",`Std. Error`)) %>%
mutate(Z=sprintf("%0.2f",`z value`)) %>%
mutate(P=ifelse(`Pr(>|z|)`<0.001,'<.001',sprintf("%0.3f",`Pr(>|z|)`))) %>%
mutate(OR=sprintf("%0.2f",exp(Estimate))) %>%
mutate(lower=sprintf("%0.2f",exp(Estimate - qnorm(0.975)*`Std. Error`))) %>%
mutate(upper=sprintf("%0.2f",exp(Estimate + qnorm(0.975)*`Std. Error`)))
mod_dat
到这里,其实结果就告一段落了,此时的结果为
5 OR (95%CI)的计算
Beta、S.E、Z、P、OR、lower和upper都有了,各位到这里就可以结束了
细心的同学可以发现,我在计算95%CI的时候用的是:
exp(Estimate - qnorm(0.975)*`Std. Error`)
这里为什么要这样去用呢,当计算的CI很贴近于1,你用1.96*S.E去计算CI,会得到Lower=1.00,更甚至会得到虽然Lower<1,Upper>1,但是P值是显著的,此时你会怀疑自己,或者怀疑软件了,本质就是因为1.96的精度不够。因此,我们采用qnorm(0.975)去计算
qnorm (0.975)= 1.959964,这个和1.96还是有差别的,用1.96去计算会得到更宽的区间,这也就解释了Lower<1,Upper>1,但是明明P值是显著的,这个结果了。
然后,我们想做的更卷一点
想想,如果,OR=1.000001,四舍五入是不是就变成了1.00,然后你又怀疑自己了,OR=1.00,但是P却显著,因此,我们继续做一个转化
mutate(OR=exp(Estimate)) %>%
mutate(OR=ifelse(OR<1.01 & OR>1 & `Pr(>|z|)`<0.05,1.01,#进行特殊情况的转化
ifelse(OR>0.99 & OR<1 & `Pr(>|z|)`<0.05,0.99,OR)))
怎么理解呢?
当OR<1.01并且OR>1,且P值显著时,那就将OR手动转为1.01,反之
当OR>0.99并且OR<1,且P值显著时,那就将OR手动转为0.99
虽然不太合适,但这样是不是可以解除误会了?
6 结果的组合
同样,我们将Lower和Upper也进行转化,我们上一套完整的代码
mod_dat <- mod %>% summary %>% coef %>% as.data.frame() %>%
mutate(Beta=sprintf("%0.2f",Estimate)) %>%
mutate(S.E=sprintf("%0.2f",`Std. Error`)) %>%
mutate(Z=sprintf("%0.2f",`z value`)) %>%
mutate(P=ifelse(`Pr(>|z|)`<0.001,'<.001',sprintf("%0.3f",`Pr(>|z|)`))) %>%
mutate(OR=exp(Estimate)) %>%
mutate(OR=ifelse(OR<1.01 & OR>1 & `Pr(>|z|)`<0.05,1.01,#进行特殊情况的转化
ifelse(OR>0.99 & OR<1 & `Pr(>|z|)`<0.05,0.99,OR))) %>%
mutate(lower=exp(Estimate - qnorm(0.975)*`Std. Error`)) %>%
mutate(lower=ifelse(lower>1 & lower<1.01 & `Pr(>|z|)`<0.05,1.01,#进行特殊情况的转化
ifelse(lower<1 & lower>0.99 & `Pr(>|z|)`<0.05,0.99,lower))) %>%
mutate(upper=exp(Estimate + qnorm(0.975)*`Std. Error`)) %>%
mutate(upper=ifelse(upper>1 & upper<1.01 & `Pr(>|z|)`<0.05,1.01,#进行特殊情况的转化
ifelse(upper<1 & upper>0.99 & `Pr(>|z|)`<0.05,0.99,upper))) %>%
mutate(ORCI=paste0(sprintf("%0.2f",OR)," (",
sprintf("%0.2f",lower)," - ",
sprintf("%0.2f",upper),")")) %>%
dplyr::select(Beta,S.E,Z,P,ORCI,OR,lower,upper) #仅保留的生成的新变量
这个结果就很清爽了
7 结果的优化
虽然清爽,但是我们还想更卷一些,看看哪里还可以卷,我们刚才的解释为:相对于Sex赋值为1的人群,Sex赋值为2的人群,OR=1.97, P值<0.001,很明显,这个结果,我们看不到参照的水平,那就继续优化
- 生成一行,变量名为sex,但是水平是1,并添加Reference字样的值
- 表格合并
7. 1参照水平
怎么获取参照水平?R语言默认以第一个水平为参照,因此只需要
levels(dat1[,"sex"])[1]
生成一行,并且将ORCI的值赋值为“1.00 (Reference)”
data.frame(Variables= c("sex",paste0(" ",levels(dat1[,"sex"])[1])),
Beta=c("",""),S.E=c("",""),Z=c("",""),P=c("",""),ORCI=c("","1.00 (Reference)"),
OR=c(NA,NA),lower=c(NA,NA),upper=c(NA,NA))
结果为
7.2 表格合并
用这个结果,再去和mod_dat(去除第一行的常数项)进行纵向合并就好了
rest_rowt <- rbind(data.frame(Variables= c("sex",paste0(" ",levels(dat1[,"sex"])[1])),
Beta=c("",""),S.E=c("",""),Z=c("",""),P=c("",""),ORCI=c("","1.00 (Reference)"),
OR=c(NA,NA),lower=c(NA,NA),upper=c(NA,NA)),
data.frame(Variables= paste0(" ",levels(dat1[,"sex"])[-1]),mod_dat[-1,])) %>% as_tibble()
rest_rowt
到这里,我们就卷完了
7.3 变量的批量处理
处理连续变量,我们就不需要额外的去生成一个新的参照水平的数据,直接就可以用了,最后我们把所有的变量进行合并输出就行了,见全部的代码
var=c("education","sex","huji","smoking","drinking")
cat_var=c("education","sex","huji","smoking","drinking")
df <- dat1
####因子转化####
for (fv in cat_var){
df[,fv] <- factor(df[,fv])
}
rest_row <- list()
for (i in 1:length(var)){
print(paste0("正在分析",var[i]),quote=F)
formula=paste0("insomnia~",var[i])
mod <- glm(formula,family = binomial,data=df)
mod_dat <- mod %>% summary %>% coef %>% as.data.frame() %>%
mutate(Beta=sprintf("%0.2f",Estimate)) %>%
mutate(S.E=sprintf("%0.2f",`Std. Error`)) %>%
mutate(Z=sprintf("%0.2f",`z value`)) %>%
mutate(P=ifelse(`Pr(>|z|)`<0.001,'<.001',sprintf("%0.3f",`Pr(>|z|)`))) %>%
mutate(OR=exp(Estimate)) %>%
mutate(OR=ifelse(OR<1.01 & OR>1 & `Pr(>|z|)`<0.05,1.01,#进行特殊情况的转化
ifelse(OR>0.99 & OR<1 & `Pr(>|z|)`<0.05,0.99,OR))) %>%
mutate(lower=exp(Estimate - qnorm(0.975)*`Std. Error`)) %>%
mutate(lower=ifelse(lower>1 & lower<1.01 & `Pr(>|z|)`<0.05,1.01,#进行特殊情况的转化
ifelse(lower<1 & lower>0.99 & `Pr(>|z|)`<0.05,0.99,lower))) %>%
mutate(upper=exp(Estimate + qnorm(0.975)*`Std. Error`)) %>%
mutate(upper=ifelse(upper>1 & upper<1.01 & `Pr(>|z|)`<0.05,1.01,#进行特殊情况的转化
ifelse(upper<1 & upper>0.99 & `Pr(>|z|)`<0.05,0.99,upper))) %>%
mutate(ORCI=paste0(sprintf("%0.2f",OR)," (",
sprintf("%0.2f",lower)," - ",
sprintf("%0.2f",upper),")")) %>%
dplyr::select(Beta,S.E,Z,P,ORCI,OR,lower,upper)
rownames (mod_dat)<-NULL
mod_dat <- mod_dat[-1,]
if (df[,var[i]] %>% is.factor ==T){
rest_row[i] <- rbind(data.frame(Variables= c(var[i],paste0(" ",levels(df[,var[i]])[1])),
Beta=c("",""),S.E=c("",""),Z=c("",""),P=c("",""),ORCI=c("","1.00 (Reference)"),
OR=c(NA,NA),lower=c(NA,NA),upper=c(NA,NA)),
data.frame(Variables= paste0(" ",levels(df[,var[i]])[-1]),mod_dat)) %>% list()
}
else if (df[,var[i]] %>% is.factor ==F){
rest_row[i] <- data.frame(Variables= var[i],mod_dat)%>% list()
}
}
final_single <- do.call(rbind,rest_row) %>% as_tibble()
结果见下图
7.4 结果的输出
结果导出word,结果
final_single <- do.call(rbind,rest_row) %>%
select(- OR,-lower,-upper)%>% as_tibble()
library("flextable")
c <- flextable(
final_single,
cwidth = 1,
cheight = 0.1,
defaults = list(),
theme_fun = theme_booktabs
)
save_as_docx(c,path="sing.DOCX")