R语言中的流行病数据分析神器

大家在分析临床数据和流调数据的时候指定在为各种模型，表格构建煞费苦心。今天就给大家介绍一个解决这些问题的R包sjPlot。此包不仅可以实现三线表的绘制，同时可以进行模型结果的可视化展示、评估。首先我们看下包的安装：

install.packages("sjPlot")
install.packages("rstanarm")##贝叶斯回归模型

接下来通过实例来看下其主要的功能：

##主要包的载入

library(sjPlot)
library(sjmisc)
library(sjlabelled)
library(ggplot2)

1. 回归模型的三线表的绘制。

回归模型包括线性模型，广义线性模型，混合线性模型等的结果在此包中都可以进行展示，我们实例就以简单的线性模型为主。

##数据的载入
data('efc')
efc <-as_factor(efc,c161sex,c172code)##将efc中的c161sex,c172code两个变量转化为因子格式
##构建模型，创建表
m1 <-lm(barthtot ~ c160age + c12hour + c161sex + c172code, data = efc)
tab_model(m1)

上表可以看到每个变量显示的是具体的描述标签，但是并没有显示具体变量的名字。那么如何进行查看，赋值，具体每个变量的标签如下：

或者包中的函数get_label(efc)["c160age"]也可以获得变量的标签。

如果没有这个标签那或者关闭显示标签，么就会直接显示变量的名字。另外表生成的格式是html，这样方便进行网页展示。

##关闭变量标签
tab_model(m1,auto.label = F)

那么问题来了，单模型既然可以展示，为了进行两个模型的比较，也可以进行多个模型结果的展示来进行结果的对比。

m2 <- lm(neg_c_7 ~ c160age + c12hour + c161sex + e17age, data =efc)
tab_model(m1,m2)

为了适应广东使用者，此包还可以对需要显示的内容进行控制，其中show.est （显示预测值），show.ci （显示置信区间），show.se （显示标准误），show.std （显示系数），show.p （显示p值），show.stat（显示系数的统计学检验值），show.df （显示自由度）。通过进行设置可以进行对需要的值进行展示。

tab_model(m1,show.se = TRUE, show.std = TRUE, show.stat = TRUE,auto.label = F)

##置信区间的合并
tab_model(m1,show.se = TRUE, show.std = TRUE, show.stat = TRUE,auto.label = F,collapse.ci =T)

##P值用*代替
tab_model(m1,show.se = TRUE, show.std = TRUE, show.stat = TRUE,auto.label = F,collapse.ci =T, p.style = "stars")

##P值科学计数法保留两位有效数字
tab_model(m1,show.se = TRUE, show.std = TRUE, show.stat = TRUE,auto.label = F,collapse.ci =T, p.style = "scientific", digits.p = 2)

##仅展示特定的自变量
tab_model(m1,show.se = TRUE, show.std = TRUE, show.stat = TRUE,auto.label = F,collapse.ci =T, p.style = "scientific", digits.p = 2,rm.terms =c("c172code2", "c161sex2"))

##表中样式颜色设置
tab_model(m1,show.se = TRUE, show.std = TRUE, show.stat = TRUE,auto.label = F,collapse.ci =T, p.style = "scientific", digits.p = 2,rm.terms =c("c172code2", "c161sex2"), CSS = list(
    css.depvarhead = 'color: red;',
    css.centeralign = 'text-align: left;',
    css.firsttablecol = 'font-weight: bold;',
    css.summary = 'color: blue;'
  )
)

##添加CSS样式
tab_model(m1,show.se = TRUE, show.std = TRUE, show.stat = TRUE,auto.label = F,collapse.ci =T, p.style = "scientific", digits.p = 2,rm.terms =c("c172code2", "c161sex2"), CSS =css_theme("cells"))

2. 回归模型结果的可视化展示

#森林图展示
plot_model(m1)

##改变中线的颜色
plot_model(m1,vline.color = "red")

##进行排序
plot_model(m1,sort.est = TRUE)

##添加值
plot_model(m1,show.values = TRUE, value.offset = .3)

##去掉标签
plot_model(m1, auto.label= F)

###整体风格改变
plot_model(m1, colors= "Accent", show.values = TRUE,value.offset = .4,value.size = 4,dot.size= 3,line.size = 1.5,vline.color = "blue",width = 1.5)

#置信区间结果展示
plot_model(m1,type = "pred", terms = "c160age")##type还有eff,emm,int

#多个自变量的展示
plot_model(m1,type = "pred", terms = c("c12hour","c160age"))

##三变量
plot_model(m1,type = "pred", terms =c("c12hour","c160age","c161sex"))

##模型评估残差可视化。主要是展示预测与真实值之间的距离结果
plot_residuals(m1)

##去除不想看的变量
plot_residuals(m1,remove.estimates = c( "c161sex","c172code"))

##贝叶斯回归模型的展示
m <-stan_glm(barthtot ~ c160age + c12hour + c161sex + c172code, data =efc,chains=1)
plot_model(m)

###自定义外观
plot_model(m, bpe= "mean",bpe.style = "dot",prob.inner = .4,prob.outer = .8)

3. 基础信息汇总三线表

##构建数据表
data(efc)
varlabs <-get_label(efc)
# recveivefirst item of COPE-index scale
start <-which(colnames(efc) == "c82cop1")
# recveive lastitem of COPE-index scale
end <-which(colnames(efc) == "c90cop9")
 
# create dataframe with COPE-index scale
mydf <-data.frame(efc[, start:end])
colnames(mydf)<- varlabs[start:end]
tab_df(mydf)

###基础信息统计
tab_itemscale(mydf)

##相关性表创建
tab_corr(mydf,p.numeric = TRUE)

##展示PCA结果
tab_pca(mydf)

4. 基础信息的可视化展示

##频率分布图
plot_grpfrq(efc$e42dep,efc$c172code, geom.colors = "gs")

##设置样式
set_theme(
    geom.outline.color ="antiquewhite4",
    geom.outline.size = 1,
    geom.label.size = 2,
    geom.label.color = "grey50",
    title.color = "red",
    title.size = 1.5,
    axis.angle.x = 45,
    axis.textcolor = "blue",
    base = theme_bw()
)
plot_grpfrq(efc$e42dep,efc$c172code)

##自定义颜色
plot_grpfrq(efc$e42dep,efc$c172code,  title = NULL,
  geom.colors = c("cadetblue", "coral"),
  geom.size = 0.4
)

#坐标轴转换
plot_grpfrq(efc$e42dep,efc$c172code ,coord.flip = TRUE)

#散点图
plot_scatter(efc,e16sex, neg_c_7)

6. 李克特量表的可视化

首先介绍下这个量表的定义：是属评分加总式量表最常用的一种，属同一构念的这些项目是用加总方式来计分，单独或个别项目是无意义的。它是由美国社会心理学家李克特于1932年在原有的总加量表基础上改进而成的。该量表由一组陈述组成，每一陈述有"非常同意"、"同意"、"不一定"、"不同意"、"非常不同意"五种回答，分别记为5、4、3、2、1，每个被调查者的态度总分就是他对各道题的回答所得分数的加总，这一总分可说明他的态度强弱或他在这一量表上的不同状态。（转自百度百科）。

mydf <-find_var(efc, pattern = "cop", out = "df")
plot_likert(mydf)

#只列出总的结果
plot_likert(
  mydf,
  grid.range = c(1.2, 1.4),
  expand.grid = FALSE,
  values = "sum.outside",
  show.prc.sign = TRUE
)

##进行分组展示
plot_likert(mydf,groups = c(2, 1, 1, 1, 1, 2, 2, 2, 1))

##组名的展示
plot_likert(
  mydf,
  c(rep("A", 4), rep("B",5)),
  sort.groups = FALSE,
  grid.range = c(0.9, 1.1),
  geom.colors = "RdBu",
  rel_heights = c(6, 8),
  wrap.labels = 40,
  reverse.scale = TRUE
)

##其它同类数据的展示
# controllegend items
six_cat_example= data.frame(
  matrix(sample(1:6, 600, replace = TRUE), ncol= 6)
)
 
six_cat_example<-
  six_cat_example %>%
  dplyr::mutate_all( ~ ordered(., labels =c("+++", "++", "+", "-","--", "---")))
 
# Old default
plot_likert(
  six_cat_example,
  groups = c(1, 1, 1, 2, 2, 2),
  group.legend.options = list(nrow = 2, byrow =FALSE)
)

8.多项式拟合结果展示

sjp.poly(efc$c160age,efc$quol_5, 1)#其中1代表1次多项式

##自动进行多个方程的拟合
sjp.poly(efc$c160age,efc$quol_5, 1:4, show.scatter = FALSE)

##解决线性模型不显著的变量
fit <-lm(tot_sc_e ~ c12hour + e17age + e42dep, data = efc)
# inspectrelationship between predictors and response
plot_model(fit,type = "slope")

#"e17age" does not seem to be linear correlated to response
# try to findappropiate polynomial. Grey line (loess smoothed)
# indicatesbest fit. Looks like x^4 has the best fit,
# however, onlyx^3 has significant p-values.
sjp.poly(fit,"e17age", 2:4, show.scatter = FALSE)

欢迎互相学习交流！