# Stata 

```stata
###导入数据
use "D:\文件文档\word\搞得事\wmq meta\meta\再做\数据分析\first.dta" 
drop num year author country meanage age_rank samples cases men_cases men_samples women_cases women_samples p se _ES _seES _LCI _UCI _WT _rsample 

###数据类型转换
encode a, g(a1)

###亚组分析
gen p=cases/samples
gen se=sqrt(p*(1-p)/samples)
gen lnp=ln(p)
gen lnse=lnp(se)
metan p se, lcols(author) random xlabel(0.2,0.4,0.6,0.8) dp(5) 
#men
gen men_cases1 = real(men_cases)
gen men_samples1 = real(men_samples)
gen p_man=men_cases1/men_samples1
gen se_man=sqrt(p_man*(1-p_man)/men_samples1)
metan p_man se_man, lcols(author) random xlabel(0.2,0.4,0.6,0.8) dp(3)
#women
gen women_cases1 = real(women_cases)
gen women_samples1 = real(women_samples)
gen p_woman=women_cases1/women_samples1
gen se_woman=sqrt(p_woman*(1-p_woman)/women_samples1)
metan p_woman se_woman, lcols(author) random xlabel(0.2,0.4,0.6,0.8) dp(3)

###回归分析
metareg p tobacco_use alcohol_use areca_betel_nut_use mean_age, wsse(se) bsest(reml)
metareg p mean_age, wsse(se) bsest(reml)  graph
metareg p tobacco_use, wsse(se) bsest(reml)  graph
graph export "D:\文件文档\word\搞得事\wmq meta\meta\再做\数据分析\Graph.svg", as(svg) name("Graph")
metareg p alcohol_use, wsse(se) bsest(reml)  graph
graph export "D:\文件文档\word\搞得事\wmq meta\meta\再做\数据分析\alcohol_use.svg", as(svg) name("Graph")
metareg p areca_betel_nut_use, wsse(se) bsest(reml)  graph
graph export "D:\文件文档\word\搞得事\wmq meta\meta\再做\数据分析\areca_betel_nut_use.svg", as(svg) name("Graph")
```



# R

```R
library(tibble)

data_plot <- tibble(study,Cases,Samples,Prevalence,lower,upper)
data_plot$` ` <- paste(rep(" ", 40), collapse = " ")
data_plot$`95% CI` <- ifelse(is.na(data_plot$Prevalence), "",
                           sprintf("%.3f to %.3f",
                                   data_plot$Lower, data_plot$Upper))
#write.csv(data_plot, "data_plot_new.csv", row.names = FALSE)
data_plot <- read.csv("zhutu.csv",header = T)
#data_plot <-data_plot[,c(1,2,3,4,5,6,7)]

#summary和overall后移
pattern <- "[0-9]"
data_plot$Study <- ifelse(grepl(pattern, data_plot$Study),
                          paste("  ", data_plot$Study),
                          data_plot$Study)
data_plot$se <- (log(data_plot$Upper) - log(data_plot$Prevalence)+0.5)/1.96

# dt_tmp <- rbind(data_plot[-1, ], data_plot[1, ])
# dt_tmp[nrow(dt_tmp), 1] <- "Overall"

#upper.sum <- meta[["upper.common"]]
#lower.sum <- meta[["lower.common"]]
#names(data_plot)
#add_row(data_plot, study = "Overall", Pre = meta[["TE.common"]], lower = lower.sum, upper = upper.sum, .after = 68)
#install.packages("forestploter")

library(forestploter)

tm <- forest_theme(base_size = 12,
                   # Confidence interval point shape, line type/color/width
                   ci_pch = 15,
                   ci_col = "#762a83",
                   ci_lty = 1,
                   ci_lwd = 1,
                   ci_Theight = 0.02, # Set an T end at the end of CI 
                   # Reference line width/type/color
                   refline_lwd = 1,
                   refline_lty = "dashed",
                   refline_col = "grey20",
                   # Vertical line width/type/color
                   vertline_lwd = 1,
                   vertline_lty = "dashed",
                   vertline_col = "grey20",
                   # Change summary color for filling and borders
                   summary_fill = "#4575b4",
                   summary_col = "#4575b4",
                   # Footnote font size/face/color
                   footnote_cex = 1,
                   footnote_fontface = "bold",
                   footnote_col = "black")
library(grid)
# pdf("plotnew.pdf", width = 12, height = 18)
plotnew.pdf <- forest(data_plot[, c(1,2,3,4,5,9,10)],
               est = data_plot$Prevalence,
               lower = data_plot$Lower,
               upper = data_plot$Upper,
               ci_column = 6,
               ref_line = 0.03, #参考线
               sizes = data_plot$se,
               #arrow_lab = c(" ", "Prevalence"),
               #is_summary = c(rep(FALSE, nrow(dt_tmp)-1), TRUE),
               is_summary =!grepl(pattern, data_plot$Study),
               xlim = c(0, 0.5),
               ticks_at = c(0, 0.1,0.2,0.3,0.4,0.5),
               footnote = "
               
       Test of overall effect:  z = 26.204  p = 0.000    I² (%) = 99.3% p = 0.000",
       theme = tm)
plotnew.pdf <- add_border(plotnew.pdf,part = "header") 
plotnew.pdf <- edit_plot(plotnew.pdf,
               row = c(44,57,70,71),
               gp = gpar(fontface = "bold"))
library(ggplot2)
ggsave("plotnew.pdf",plotnew.pdf,width = 350, height = 490, units = "mm" )

# dev.new()
dev.off()
#data_plot <- data[ ,c(2,3,7,8)]
# rm(list=ls())
# fix(data_plot)
# edit(data_plot)
data_plot$Cases <- format(data_plot$Cases,justify="centre" )


###subgroup
dt <- read.csv("亚组率.csv",header = T)
dt$Subgroup <- ifelse(is.na(dt$Prevalence), 
                      dt$Subgroup,
                      paste0("   ", dt$Subgroup))
dt$` ` <- paste(rep(" ", 20), collapse = " ")
dt$`95% CI` <- ifelse(is.na(dt$Prevalence), "",
                           sprintf("%.3f to %.3f",
                                   dt$lower, dt$upper))
dt$se <- (log(dt$upper) - log(dt$Prevalence))/1.96
# dt$Prevalence <- ifelse(is.na(dt$Prevalence), "", dt$Prevalence)
# dt$lower <- ifelse(is.na(dt$lower), "", dt$lower)
# dt$upper <- ifelse(is.na(dt$upper), "", dt$upper)

p <- forest(dt[,c(1,2,5,6)],
            est = dt$Prevalence,
            lower = dt$lower, 
            upper = dt$upper,
            #sizes = dt$se,
            ci_column = 3,
            ref_line = 0.05,
            #arrow_lab = c("Placebo Better", "Treatment Better"),
            xlim = c(0, 0.3),
            ticks_at = c(0, 0.1,0.2),
            theme = tm)
            #footnote = "This is the demo data. Please feel free to change\nanything you want.")
p<- add_border(p,part = "header") 
p <- edit_plot(p,
                         row = c(1,4,11),
                         gp = gpar(fontface = "bold"))
library(ggplot2)
ggsave("p.pdf",p,width = 160, height = 120, units = "mm" )
```