如何用R語言畫森林圖展示Logistic回歸分析的結(jié)果

    

就采用表格的形式展示Logistic回歸分析的結(jié)果，上述表格把有統(tǒng)計學(xué)意義的結(jié)果進行了加粗，使得讀者看起來不那么費勁。那么，有沒有更加直觀的方法展示回歸結(jié)果呢？當(dāng)然有，那就是森林圖。近年來，越來越多文獻用森林圖來展示回歸的結(jié)果。接下來我們一起來學(xué)習(xí)一下如何用R作森林圖。

森林圖展示的數(shù)據(jù)通常是Logistic回歸分析的系數(shù)和95%置信區(qū)間以及顯著性檢驗的P值，那么如何獲得這些結(jié)果呢？

logistic回歸分析的代碼

data(Affairs,package = "AER")
df<-Affairs
df$ynaffairs<-ifelse(df$affairs>0,1,0)
df$ynaffairs<-factor(df$ynaffairs,
                     levels = c(0,1),
                     labels = c("No","Yes"))
fit.full<-glm(ynaffairs~gender+age+yearsmarried+
                children+religiousness+education+occupation+rating,
              data=df,family = binomial())

fit.result<-summary(fit.full)
df1<-fit.result$coefficients
df2<-confint(fit.full)
df3<-cbind(df1,df2)
df4<-data.frame(df3[-1,c(1,4,5,6)])
df4$Var<-rownames(df4)
colnames(df4)<-c("OR","Pvalue","OR_1","OR_2","Var")
df5<-df4[,c(5,1,2,3,4)]
df5$OR_mean<-df5$OR
df5$OR<-paste0(round(df5$OR,2),
               "(",
               round(df5$OR_1,2),
               "~",
               round(df5$OR_2,2),
               ")")
df5$Pvalue<-round(df5$Pvalue,3)
write.csv(df5,file = "forestplot_example.csv",
          quote = F,row.names = F)

導(dǎo)出數(shù)據(jù)以后需要自己手動添加一行， 最終作圖的數(shù)據(jù)如下

這里準(zhǔn)備數(shù)據(jù)的過程稍微有些繁瑣了，不知道大家有沒有簡便的方法呢？歡迎留言討論呀！

首先是安裝

install.packages("forestplot")

讀入數(shù)據(jù)并作圖

library(forestplot)
fp<-read.csv("forestplot_example.csv",header=T)

forestplot(labeltext=as.matrix(fp[,1:3]),
           mean=fp$OR_mean,
           lower=fp$OR_1,
           upper=fp$OR_2,
           zero=0,
           boxsize=0.2,
           graph.pos=2)

 

接下來是簡單的美化

forestplot(labeltext=as.matrix(fp[,1:3]),
           mean=fp$OR_mean,
           lower=fp$OR_1,
           upper=fp$OR_2,
           zero=0,
           boxsize=0.2,
           lineheight = unit(7,'mm'),
           colgap=unit(2,'mm'),
           lwd.zero=1.5,
           lwd.ci=2, 
           col=fpColors(box='#458B00',
                        summary='#8B008B',
                        lines = 'black',
                        zero = '#7AC5CD'),
           xlab="OR",
           lwd.xaxis =1,
           txt_gp = fpTxtGp(ticks = gpar(cex = 0.85),
                            xlab  = gpar(cex = 0.8),
                            cex = 0.9),
           lty.ci = "solid",
           title = "Forestplot", 
           line.margin = 0.08,
           graph.pos=2)