In cuijj6/obstetrics:

knitr::opts_chunk$set(echo = TRUE)

R Markdown

source("test/main_pkgs.R")
load_all()
last_row <- . %>% .[nrow(.), ]

fix_height <- function(d) {
    height1 <- d$height1
    height2 <- d$height2

    I <- is.na(height1)
    height1[I] <- height2[I]

    I <- is.na(height2)
    height2[I] <- height1[I]

    d$height1 <- height1
    d$height2 <- height2
    d
}

BMI <- function(height, weight){
    weight/(height/100)^2
}

df <- read_xlsx("INPUT/GDM复发_final_20191013.xlsx") %>% data.table()
df[, GDM := factor(GDM_1 * 2 + GDM_2, 0:3, c("health", "new", "No Recurrence", "Recurrence"))]
types <- c("Recurrence", "No Recurrence")
df <- df[GDM %in% types, ] %>% fix_height()

df %<>% plyr::mutate(
    BMI_befo2 = BMI(height2, wght_b2), 
    BMI_deli2 = BMI(height2, wght_d2), 
    BMI_befo1 = BMI(height1, wght_b1), 
    BMI_deli1 = BMI(height1, wght_d1), 
    delvdate1 = ymd(delvdate1), 
    delvdate2 = ymd(delvdate2), 
    interval_deli = as.numeric(difftime(delvdate2, delvdate1, units = "days")/30),
    delta_BMI_befo = BMI_befo2 - BMI_befo1, 
    delta_BMI_deli = BMI_deli2 - BMI_deli1, 
    delta_wght2     = wght_d2 - wght_b2, 
    delta_wght_befo = wght_b2 - wght_b1, 
    delta_wght_deli = wght_d2 - wght_d1, 
    baby_age  = pweek2 * 7 + pday2,
    OGTT1_0_d = as.numeric(OGTT1_0_d)
    )
names(df)
# View(d_diff)

# c("GDM_1", "GDM_2", "group", "admin", "serinum", "name", "birth", 
# "age2", "delvdate2", "famhis2", "insulin2", "hypoins2", "gender2_1", "weight2_1", "brth2_1", "out2_1", "gender2_2", "weight2_2", "brth2_2", "out2_2", "ptimes2", "dtimes2", "pweek2", "pday2", "mode2", "wght_b2", "wght_d2", "height2", "diagn2", "hist2", "HBA2_d", "OGTT2_0_d", "OGTT2_1_d", "OGTT2_2_d", "CHOL2_d", "TG2_d", "HDL2_d", "LDL2_d", "uglu2_d", "pro2_d", "GLU2_d", "HBA2_a", 
# "OGTT2_0_a", "OGTT2_1_a", "OGTT2_2_a", "CHOL2_a", "TG2_a", "HDL2_a", "LDL2_a", "insu2_0_a", "uglu2_1_a", "uglu2_2_a", 
# # Index
# "age1", "delvdate1", "gender1_1", "weight1_1", "brth1_1", "out1_1", "gender1_2", "weight1_2", "brth1_2", "out1_2", "insulin1", "hypoins1", "ptimes1", "dtimes1", "pweek1", "pday1", "mode1", "wght_b1", "wght_d1", "height1", "diagn1", "pdiagn1", "hist1", "HBA1_d", "OGTT1_0_d", "OGTT1_1_d", "OGTT1_2_d", "OGTT1_3_d", "CHOL1_d", "TG1_d", "HDL1_d", "LDL1_d", "GDM")

df$his  <- tidy_his(df$famhis2)
df$mode <- tidy_mode(df$mode2)

# -d: pregnancy (only for OGTT)
# -a: after
# df[, .(height1,height2)]
varnames = c(
    "his", "mode",
    # "hypoins2",
    "weight2_1",
    "age2", "ptimes2", "dtimes2",
    # 2. 体重信息
    "wght_b2", "wght_d2",
    "interval_deli", 
    "delta_wght2",
    "delta_wght_befo",
    "delta_wght_deli",
    "baby_age",
    # BMI
    "BMI_befo2", "BMI_deli2",
    "delta_BMI_deli",
    "delta_BMI_befo",
    # 3. 测量指标
    "OGTT1_0_d", "OGTT1_1_d", "OGTT1_2_d", 
    "CHOL2_d", "TG2_d", "HDL2_d", "LDL2_d"
    ) %>% set_names(., .)
lst <- foreach(indicator = varnames) %do% {
    basinInfo(df, indicator)    
}

info <- melt_list(lst, "variales")
write_list2xlsx(list("basinInfo" = info), "Table.1 BasicInfo.xlsx")
info

info_sign <- info[pval < 0.05]

# INPUT data for logistics reg
data = df[, .SD, .SDcols = c("GDM", info_sign$variales) %>% setdiff(c("OGTT1_1_d"))]
data$GDM %<>% as.character() %>% factor()
m <- glm(GDM~., binomial(), data)

digit = 2; fmt3 = sprintf("%%.%df [%%.%df, %%.%df]", digit, digit, digit)
l_pr <- prLogisticDelta(GDM~., data = data, pattern="marginal")
pr <- l_pr$ci %>% set_colnames(c("y", "min", "max")) %>% data.table() %>% .[, sprintf(fmt3, y, min, max)]

variables <- l_pr$ci %>% rownames()

or   <- logistic.display(m)
d_or <- or$table %>% {.[seq(1, nrow(.)-1, 2), ]} %>% .[, c(2, 1, 3, 4)]%>% data.table() %>% {
    .[[3]] %<>% as.numeric()
    .[[4]] %<>% as.numeric()
    .
}

info_table2 <- cbind(variables, pr, d_or)
write_list2xlsx(list("Table1" = info_table2), "Table2. Logistics Regression.xlsx")

## 
s <- summary(m)
probabilities <- m %>% predict(data, type = "response")
pred <- ifelse(probabilities > 0.5, types[1], types[2])

tbl = data.table(ref = data$GDM, pred) %>% table()
precision = precision(tbl, "Recurrence")  # 预测阳性中 正确的比例
recall    = recall(tbl, "Recurrence")     # 样本中 正确预测的比例

FN  = tbl[2, 1]/sum(tbl[, 1])             # 在所有实际为阴性的样本中，被错误地判断为阳性之比率。
ACC = sum(diag(tbl))/sum(as.numeric(tbl))

library(prLogistic)
library(epiDisplay)

logistic.display(m)

cuijj6/obstetrics documentation built on Sept. 7, 2020, 12:04 a.m.