library(openxlsx)
library(tidyverse)
library(pROC)
library(nricens)
library(PredictABEL)
library(geepack)
library(gee)
adam<- read.xlsx("C:/Users/99405/Desktop/增长轨迹/xlsx数据/12.17sjydata6.xlsx")

pre <- read.xlsx("C:/Users/99405/Desktop/整理/预测分析.xlsx")
adam <- adam %>% mutate(rebp_am_resultclass2=if_else(BP_AM_result==3,1,0))
adam <- adam %>% mutate(adophysical=if_else(cl3bmi==0,0,1))
adam1 <- adam %>% select(where(~sum(is.na(.)) <= nrow(adam)/2))
adam <- sjydata %>% select(overweight,underweight,BP_AM_result,bp_ch_result,matches("weizm\\d+"),matches("\\w+x\\d_\\d"),contains("slope"))
bpwei <- read.xlsx("C:/Users/99405/Desktop/rebp预测.xlsx")
bplen <- read.xlsx("C:/Users/99405/Desktop/lenpred.xlsx")
bpweipred <- bpwei %>% left_join(adam[c("childid","rebp_am_resultclass2")],by="childid") 
bplenpred <- bplen %>% left_join(adam[c("childid","rebp_am_resultclass2")],by="childid") 


#所用代码####
fit1 <- glm(overweight~weizm24,family=binomial,data=adam,x=TRUE)
adam1$pfit1 <- predict(fit1,adam1,type="response")
roc_result1 <- roc(adam1[["overweight"]], adam1[["pfit1"]], ci = TRUE, auc = TRUE)
aucciresult1 <- ci.auc(roc_result1)
aucresult1 <- tibble(auc=aucciresult1[2],
                     Lower_CI = aucciresult1[1],
                     Upper_CI = aucciresult1[3],
                     model=paste(fit1$formula,collapse = " "))


#批量化建模
#*所用模型公式####
models <- list(
#1.LEN WEI BMI 预测能力的差异
  #wei
  formula = as.formula(overweight ~ weizm6),
  formula = as.formula(overweight ~ weizm6 + Slope_WEI_240),
  formula = as.formula(overweight ~ weizm6 + Slope_WEI_240 + weirrerm0_6),
  #len
  formula = as.formula(overweight~lenzm6),
  formula = as.formula(overweight ~ lenzm6 + Slope_LEN_240),
  formula = as.formula(overweight ~ lenzm6 + Slope_LEN_240 + lenrrerm0_6),
  #bmi
  formula = as.formula(overweight~bmizm6),
  formula = as.formula(overweight ~ bmizm6 + Slope_BMI_240),
  formula = as.formula(overweight ~ bmizm6 + Slope_BMI_240 + bmirrerm0_6),
#2.不同时点6，12，18，24月预测能力的差异
  #6月
  formula = as.formula(overweight ~ weizm6),
  formula = as.formula(overweight ~ weizm6 + Slope_WEI_240),
  formula = as.formula(overweight ~ weizm6 + Slope_WEI_240 + weirrerm0_6),
  #12月
  formula = as.formula(overweight ~ weizm12),
  formula = as.formula(overweight ~ weizm12 + Slope_WEI_480),
  formula = as.formula(overweight ~ weizm12 + Slope_WEI_480 + weirrerm6_12),
  #18月
  formula = as.formula(overweight ~ weizm18),
  formula = as.formula(overweight ~ weizm18 + Slope_WEI_720),
  formula = as.formula(overweight ~ weizm18 + Slope_WEI_720 + weirrerm12_18),
  #24月
  formula = as.formula(overweight ~ weizm24),
  formula = as.formula(overweight ~ weizm24 + Slope_WEI_960),
  formula = as.formula(overweight ~ weizm24 + Slope_WEI_960 + weirrerm18_24),
#3.不同指标预测能力的差异
  #12月 
  formula = as.formula(overweight~weizm12),
  formula = as.formula(overweight~Slope_WEI_480),
  formula = as.formula(overweight~weirrerm6_12),
  #24月 
  formula = as.formula(overweight~weizm24),
  formula = as.formula(overweight~Slope_WEI_960),
  formula = as.formula(overweight~weirrerm18_24)
)



#*循环并储存结果####
# 初始化一个空的数据框来存储结果
aucresult <- tibble()

# 使用循环构建模型、进行预测、计算 ROC 曲线和 AUC，然后将结果添加到数据框
for (i in seq_along(models)) {
  # 构建模型
  fit <- glm(models[[i]], family = binomial(), data = adam1, x = TRUE)
  
  # 进行预测
  pred_column <- paste("pfit", i, sep = "")
  adam1[[pred_column]] <- predict(fit, adam1, type = "response")
  
  # 计算 ROC 曲线和 AUC
  roc_result <- roc(adam1[["overweight"]], adam1[[pred_column]], ci = TRUE, auc = TRUE)
  auc_ci_result <- ci.auc(roc_result)
  
  # 将结果添加到数据框
  auc_result <- tibble(
    auc = auc_ci_result[2],
    Lower_CI = auc_ci_result[1],
    Upper_CI = auc_ci_result[3],
    model = paste(models[[i]],collapse = " "))
  
  aucresult <- rbind(aucresult, auc_result)
}


#输出为xlsx
write.xlsx(aucresult,"C:/Users/99405/Desktop/aucresult.xlsx")



adam1$rebp_am_resultclass2


#血压
# 出生体重和出生身长切割为分类变量
adam1$weizm0
cols_to_cut <- c("weizm0","lenzm0","bmizm0")
adam1 <- adam1 %>%
  mutate(across(any_of(cols_to_cut), ~cut(., breaks = c(min(.,na.rm=T)-1,quantile(., probs = c(0.25, 0.50, 0.75,1),na.rm = TRUE)),labels=F), .names = "{.col}class4"))
adam1 <- adam1 %>%
  mutate(across(any_of(cols_to_cut), ~cut(., breaks = c(min(.,na.rm=T)-1,quantile(., probs = c(1/3,2/3,1),na.rm = TRUE)),labels=F), .names = "{.col}class3"))
adam1 <- adam1 %>%
  mutate(across(any_of(cols_to_cut), ~cut(., breaks = c(min(.,na.rm=T)-1,quantile(., probs = c(0.1,0.3,0.5,0.7,0.9,1),na.rm = TRUE)),labels=F), .names = "{.col}class6"))



model1 <- gee(formula = rebp_am_resultclass2 ~ weizm0+ ag117 + M_edu + F_edu + M_job + F_job + mother_age + ab106 + parity_c2 + pcaWI_C3 + treat_chris + ac203 ,id=ar06,subset=weizm0class4,family = binomial,corstr = "independence",data=adam1)
model2 <- glm(formula = rebp_am_resultclass2 ~ weizm0 ,family = binomial,data=adam1)

class(adam1$weizm0class4)

table(adam1$weizm0class3)
table(adam1$lenzm0class4)
table(adam1$weix9_1class3)

groupadam <- adam1 %>% group_nest(weizm0class4,na.rm=F) 
groupadam <- groupadam[1:4,]
result <-groupadam %>% mutate(model=map(data,~glm(rebp_am_resultclass2 ~ weizm0 ,data=.x,family=binomial)))
result <- result %>% mutate(result=map(model,tidy)) %>% select(weizm0class4,result) %>% unnest(result)

groupadam <- adam1 %>% group_nest(weizm0class3,na.rm=F) 
groupadam <- groupadam[1:3,]
result2 <-groupadam %>% mutate(model=map(data,~glm(rebp_am_resultclass2 ~ weizm0 ,data=.x,family=binomial)))
result2 <- result2 %>% mutate(result=map(model,tidy)) %>% select(weizm0class3,result) %>% unnest(result)

groupadam <- adam1 %>% group_nest(weizm0class6,na.rm=F) 
groupadam <- groupadam[1:5,]
result3 <-groupadam %>% mutate(model=map(data,~glm(rebp_am_resultclass2 ~ weizm0 ,data=.x,family=binomial)))
result3 <- result3 %>% mutate(result=map(model,tidy)) %>% select(weizm0class6,result) %>% unnest(result)



rebp_am_resultclass2 ~ weizm0+ 
groupadam <- adam1 %>% group_nest(weizm0class3,na.rm=F) 
groupadam <- groupadam[1:3,]
result2 <-groupadam %>% mutate(model=map(data,~gee(rebp_am_resultclass2 ~ weizm0+ag117 + M_edu + F_edu + M_job + F_job + mother_age + ab106 + parity_c2 + pcaWI_C3 + treat_chris + ac203 ,id=ar06,corstr = "independence",data=.x,family=binomial)))
result2 <- result2 %>% mutate(result=map(model,tidy)) %>% select(weizm0class3,result) %>% unnest(result)

# wei

groupadam <- adam1 %>% group_nest(weizm0class3,na.rm=F) 
groupadam <- groupadam[1:3,]
result2 <-groupadam %>% mutate(model=map(data,~glm(rebp_am_resultclass2 ~ weizm0 ,data=.x,family=binomial)))
result2a <- result2 %>% mutate(result=map(model,tidy)) %>% select(weizm0class3,result) %>% unnest(result)
result2 <-groupadam %>% mutate(model=map(data,~glm(rebp_am_resultclass2 ~ weizm0+ag117 + M_edu + F_edu + M_job + F_job + mother_age + ab106 + parity_c2 + pcaWI_C3 + treat_chris + ac203 ,data=.x,family=binomial)))
result2b <- result2 %>% mutate(result=map(model,tidy)) %>% select(weizm0class3,result) %>% unnest(result)

result2a <- result2a %>% filter(term=="weizm0")
result2b <- result2b %>% filter(term=="weizm0")
weiresult2 <- result2a %>% left_join(result2b,by="weizm0class3")


# len
groupadam <- adam1 %>% group_nest(lenzm0class3,na.rm=F) 
groupadam <- groupadam[1:3,]
result2 <-groupadam %>% mutate(model=map(data,~glm(rebp_am_resultclass2 ~ lenzm0 ,data=.x,family=binomial)))
result2a <- result2 %>% mutate(result=map(model,tidy)) %>% select(lenzm0class3,result) %>% unnest(result)
result2 <-groupadam %>% mutate(model=map(data,~glm(rebp_am_resultclass2 ~ lenzm0+ag117 + M_edu + F_edu + M_job + F_job + mother_age + ab106 + parity_c2 + pcaWI_C3 + treat_chris + ac203  ,data=.x,family=binomial)))
result2b <- result2 %>% mutate(result=map(model,tidy)) %>% select(lenzm0class3,result) %>% unnest(result)

result2a <- result2a %>% filter(term=="lenzm0")
result2b <- result2b %>% filter(term=="lenzm0")
lenresult2 <- result2a %>% left_join(result2b,by="lenzm0class3")


# bmi
groupadam <- adam1 %>% group_nest(bmizm0class3,na.rm=F) 
groupadam <- groupadam[1:3,]
result2 <-groupadam %>% mutate(model=map(data,~glm(rebp_am_resultclass2 ~ bmizm0 ,data=.x,family=binomial)))
result2a <- result2 %>% mutate(result=map(model,tidy)) %>% select(bmizm0class3,result) %>% unnest(result)
result2 <-groupadam %>% mutate(model=map(data,~glm(rebp_am_resultclass2 ~ bmizm0+ag117 + M_edu + F_edu + M_job + F_job + mother_age + ab106 + parity_c2 + pcaWI_C3 + treat_chris + ac203  ,data=.x,family=binomial)))
result2b <- result2 %>% mutate(result=map(model,tidy)) %>% select(bmizm0class3,result) %>% unnest(result)

result2a <- result2a %>% filter(term=="bmizm0")
result2b <- result2b %>% filter(term=="bmizm0")
bmiresult2 <- result2a %>% left_join(result2b,by="bmizm0class3")



write.xlsx(weiresult2,"C:/Users/99405/Desktop/weizm0class3.xlsx")
write.xlsx(lenresult2,"C:/Users/99405/Desktop/lenzm0class3.xlsx")
write.xlsx(bmiresult2,"C:/Users/99405/Desktop/bmizm0class3.xlsx")


quantile(adam1$weizm0, probs = c(1/3,2/3,1),na.rm=T)
quantile(adam1$bmizm0, probs = c(1/3,2/3,1),na.rm=T)      

bppredvars <- colnames(bpweipred)[-c(1:4,39)]
roclist <- list()

for (nm in bppredvars) {
  roc_result <- roc(bpweipred$rebp_am_resultclass2, bpweipred[[nm]], ci = TRUE, auc = TRUE)
  roclist[nm] <- list(roc_result)}  

aucresult <- tibble()
for(i in bppredvars){
  result <- auc(roclist[[i]])
  resultci <- ci.auc(roclist[[i]])
  idx <- match(i, bppredvars) 
  aucresult[idx,"varname"] <- i
  aucresult[idx,"auv"] <- as.numeric(result)
  aucresult[idx,"auc_lci"] <-resultci[1]
  aucresult[idx,"auc_uci"] <-resultci[3]}


write.xlsx(aucresult,"C:/Users/99405/Desktop/aucsinglewei.xlsx")

bppredvars <- colnames(bplenpred)[-c(1,36)]
roclist <- list()

for (nm in bppredvars) {
  roc_result <- roc(bplenpred$rebp_am_resultclass2, bplenpred[[nm]], ci = TRUE, auc = TRUE)
  roclist[nm] <- list(roc_result)}  

aucresult <- tibble()
for(i in bppredvars){
  result <- auc(roclist[[i]])
  resultci <- ci.auc(roclist[[i]])
  idx <- match(i, bppredvars) 
  aucresult[idx,"varname"] <- i
  aucresult[idx,"auv"] <- as.numeric(result)
  aucresult[idx,"auc_lci"] <-resultci[1]
  aucresult[idx,"auc_uci"] <-resultci[3]}


write.xlsx(aucresult,"C:/Users/99405/Desktop/aucsinglelen.xlsx")

adam$bp_am_resultclass2

bc <- adam %>% select(childid,bp_am_resultclass2)
pre <- pre %>% left_join(bc,by="childid")

bppredvars <- colnames(pre)[-1]
roclist <- list()

for (nm in bppredvars) {
  roc_result <- roc(pre$bp_am_resultclass2, pre[[nm]], ci = TRUE, auc = TRUE)
  roclist[nm] <- list(roc_result)}  

aucresult <- tibble()
for(i in bppredvars){
  result <- auc(roclist[[i]])
  resultci <- ci.auc(roclist[[i]])
  idx <- match(i, bppredvars) 
  aucresult[idx,"varname"] <- i
  aucresult[idx,"auv"] <- as.numeric(result)
  aucresult[idx,"auc_lci"] <-resultci[1]
  aucresult[idx,"auc_uci"] <-resultci[3]}


write.xlsx(aucresult,"C:/Users/99405/Desktop/preauc.xlsx")
d