R语言 PSM之前之后使用Table1输出组间信息对比

2023-06-07  本文已影响0人  王叽叽的小心情

刚开始使用table1,后面使用的是tableone,先备忘录一下

library(tableone)  
library(Matching)  
library(gtools) # For the logit function  
library(lmtest)  
library(sandwich)
library(MatchIt)
library(dplyr)
library(table1) 
library(printr, quietly=TRUE) # to save table1 plot as csv
library(kableExtra)
library(magick)

pvalue <- function(x, ...) {
  # Construct vectors of data y, and groups (strata) g
  y <- unlist(x)
  g <- factor(rep(1:length(x), times=sapply(x, length)))
  if (is.numeric(y)) {
    # For numeric variables, perform a standard 2-sample t-test
    p <- t.test(y ~ g)$p.value
  } else {
    # For categorical variables, perform a chi-squared test of independence
    p <- chisq.test(table(y, g))$p.value
  }
  # Format the p-value, using an HTML entity for the less-than sign.
  # The initial empty string places the output on the line below the variable label.
  c("", sub("<", "&lt;", format.pval(p, digits=3, eps=0.001)))
}

# read data
df <- read.csv("E:\\panel_data.csv", header = 1)
df <- na.omit(df) 
# generate the treatment variable
df <- df %>% mutate(trent = ifelse(df$nrty > 2014, 1, 0))
fig_path = "E:\\DID20230606\\"

# before PSM
t_before <- table1(~ df[['ntl']] + poi_comm + poi_trans + lucc_forest + lucc_grass + 
                     poi_publi | treatment, data=df, overall=F, 
                   extra.col=list(`P-value`=pvalue))
# save the figure
t_before[1] %>% save_kable(file = paste(fig_path, "ntl", "_table_before_psm.png"))
# save the numbers
write.csv(as.data.frame(t_before), paste(fig_path, "ntl", "_table_before_psm.csv"), 
          row.names=TRUE)

# PSM
set.seed(50)
psm <- matchit(data = df,
               formula = treatment ~ poi_comm + poi_trans + lucc_forest + 
                 lucc_grass + poi_publi,
               method = "nearest",
               distance = "logit",
               replace = FALSE,
               caliper = 0.05)

# after matched
df_matched <- match.data(psm)

t_after <- table1( ~ ntl + poi_comm + poi_trans + lucc_forest + lucc_grass + 
                     poi_publi | treatment, data=df_matched, 
                   overall=F, extra.col=list(`P-value`=pvalue))
# save the figure
t_after[1] %>% save_kable(file = paste(fig_path, "ntl", "_table_after_psm.png"))
# save the numbers
write.csv(as.data.frame(t_after), paste(fig_path, "ntl", "_table_after_psm.csv"), 
          row.names=TRUE)

上一篇下一篇

猜你喜欢

热点阅读