## Created by EmpowerStats @ Mon, 15 Sep 25 16:29:18 +0800## 
#******************** Regarding ALL Following R Code   ******************** 
#*** COPYRIGHT (c) 2010, 2021 X&Y Solutions, ALL RIGHT RESERVED *********** 
#******************* www.EmpowerStats.com ********************************* 
#************************************************************************** 
Sys.setlocale(category = 'LC_ALL', locale = 'English_United States.1252'); 
.libPaths(file.path(R.home(),'library')); 
library(doBy); 
options(timeout=600); 
library(plotrix); 
library(stringi); 
library(stringr); 
library(survival); 
library(rms); 
library(nnet); 
library(car); 
library(mgcv); 
pdfwd<-6; pdfht<-6; 
load('D:/EmpowerRCH/Analysis/zfg2/longitudinal 2010.Rdata'); 
if (length(which(ls()=='EmpowerStatsR'))==0) EmpowerStatsR<-get(ls()[1]); 
names(EmpowerStatsR)<-toupper(names(EmpowerStatsR)); 
originalVNAME<-names(EmpowerStatsR); 
ofname<-'zfg2_176_tbl'; 
attach(EmpowerStatsR); 
sink(paste(ofname,'_datastep.lst',sep='')); 
print('Creating new variable: BMI2'); 
BMI2<-BMI*BMI; 
summary(BMI2); 
EmpowerStatsR<-cbind(EmpowerStatsR,BMI2); 
print('Creating new variable: CUNBAE'); 
CUNBAE<-rep(NA, times=nrow(EmpowerStatsR)); 
tmp <- is.na(CUNBAE) & (SEX==0); 
tmp <- (!is.na(tmp) & tmp); 
CUNBAE[tmp]<-((0.503*AGE)+10.689+(3.172*BMI)-(0.026*BMI2)+(0.181*BMI)-(0.02*BMI*AGE)-(0.005*BMI2)+(0.00021*BMI2*AGE)-44.988)[tmp]; 
tmp <- is.na(CUNBAE) & (SEX==1); 
tmp <- (!is.na(tmp) & tmp); 
CUNBAE[tmp]<-((0.503*AGE)+(3.172*BMI)-(0.026*BMI2)-(0.02*BMI*AGE)+(0.00021*BMI2*AGE)-44.988)[tmp]; 
summary(CUNBAE); 
EmpowerStatsR<-cbind(EmpowerStatsR,CUNBAE); 
rm(BMI2,CUNBAE); 
sink(); 
vname<-c(NA,'SEX','SEX.0','SEX.1','AGE','NAFLD','NAFLD.0','NAFLD.1','BMI','BMI2','CUNBAE')[-1]; 
vlabel<-c(NA,'SEX','  0','  1','AGE','NAFLD','  0','  1','BMI','BMI2','CUNBAE')[-1]; 
varused4this <- c('SEX','AGE','NAFLD','BMI','BMI2','CUNBAE'); 
pkgs<-c('pROC','survivalROC','plyr','rms','rmda','Hmisc'); 
for (g in pkgs) {  
if (!(g %in% rownames(installed.packages()))) install.packages(g,repos='https://cloud.r-project.org'); 
}
library(pROC); 
library(survivalROC); 
library(plyr); 
library(rms); 
library(rmda); 
library(Hmisc); 
WD <- EmpowerStatsR; rm(EmpowerStatsR); gc(); 
title<-'诊断试验与ROC分析0'; 
wd.subset=''; 
colvname<-'NAFLD'; colvlv<- 2; 
xvname<-c('CUNBAE','BMI'); 
xlv<-c(0,0); 
cox<- 0; 
timevar<- NA; 
vname.start<- NA; 
par3<- NA; 
par1<-2; 
avname<-c(); 
chk<- 0; 
bvar<- NA; 
par2<- NA; 
dec<-4; 
##R package## pROC survivalROC plyr rms rmda Hmisc ##R package##;
pvformat<-function(p,dec) {
  pp <- sprintf(paste("%.",dec,"f",sep=""),as.numeric(p))
  if (is.matrix(p)) {pp<-matrix(pp, nrow=nrow(p)); colnames(pp)<-colnames(p);rownames(pp)<-rownames(p);}
  lw <- paste("<",substr("0.00000000000",1,dec+1),"1",sep="");
  pp[as.numeric(p)<(1/10^dec)]<-lw
  return(pp)
}
numfmt<-function(p,dec) {
  if (is.list(p)) p<-as.matrix(p)
  pp <- sprintf(paste("%.",dec,"f",sep=""),as.numeric(p))
  if (is.matrix(p)) {pp<-matrix(pp, nrow=nrow(p));colnames(pp)<-colnames(p);rownames(pp)<-rownames(p);}
  pp[as.numeric(p)>10000000]<- "inf."
  pp[is.na(p) | gsub(" ","",p)==""]<- ""
  pp[p=="-Inf"]<-"-Inf"
  pp[p=="Inf"]<-"Inf"
  return(pp)
}
rocplot<-function(x,xb) { 
 xlm<-c(0,1);
 x1 <- 1-x$specificities; y1<-x$sensitivities;
 if (x$percent) {x1<-x1/100; y1<-y1/100;}
 plot(x1, y1, xlim=c(0, 1), ylim=c(0, 1), xlab="1 - Specificity", ylab="Sensitivity",
  asp=0.9, mar=c(4, 4, 2, 2)+.1, mgp=c(2.5, 1, 0), col=par("col"),
  lty=par("lty"), lwd=2, type="l", main=paste(xb, "\nAUC =", round(x$auc,3)));
 abline(a=0,b=1,col="darkgrey",lty=1,lwd=1);
}
rocplot2<-function(x,xb,mainb) {
 n=length(x); xlm<-c(0,1); ccol<-(1:n)
 x1 <- 1-x[[1]]$specificities; y1<-x[[1]]$sensitivities; auc<-x[[1]]$auc
 if (x[[1]]$percent) {x1<-x1/100; y1<-y1/100;}
 if (n==1) {mb<-paste(xb[1],"\nAUC =", round(auc,3));
 } else {mb<-mainb;}
 plot(x1, y1, xlim=xlm, ylim=xlm,
 xlab="1 - Specificity", ylab="Sensitivity",
 asp=0.9, mar=c(4, 4, 2, 2)+.1, mgp=c(2.5, 1, 0),
 col=1, lty=1, lwd=2, type="l", main=mb);
 abline(a=0,b=1,col="darkgrey",lty=1,lwd=1);
 for (k in (2:n)) {
   auc<-c(auc,x[[k]]$auc)
   par(new=T)
   xk <- 1-x[[k]]$specificities; yk<-x[[k]]$sensitivities
   if (x[[k]]$percent) {xk<-xk/100; yk<-yk/100;}
   plot(xk, yk, xlim=xlm, ylim=xlm, xlab="", ylab="",
   asp=0.9, mar=c(4, 4, 2, 2)+.1, mgp=c(2.5, 1, 0),
   col=k, lty=1, lwd=2, type="l", main="");
 }
 if (n>1) {
   mlen<-min(max(nchar(xb)),15); 
   mspace<-paste(rep(" ",mlen),collapse="")
   xb1<-substr(paste(xb,mspace),1,mlen)
   tauc<-paste(xb1,pvformat(auc,3),sep=": ")
   ord<-order(auc,decreasing=TRUE)
   tauc<-tauc[ord]; clr<-(1:n)[ord]
   legend(0.6,0.4,tauc,title="AUC",col=clr,lty=1,lwd=1,bty="n")
 }
}
diagss<-function(d,x,xb,xlb,ylb,pngf) {
  cmp<- !is.na(d) & !is.na(x); d<-d[cmp]; x<-x[cmp]
  d0<-d; d<-2-d;  ndL <- length(levels(factor(d))); nxL <- length(levels(factor(x)));  
  if (ndL!=2 | nxL<2) return (c(xb,rep(" ",6)))
  lr <- rep(NA,nxL)
  if (nxL==2) {x<-2-x; xlb<-c(xlb[2],xlb[1])}
  t0<-table(x,d);  ctot<-apply(t0,2,sum);  rtot<-apply(t0,1,sum)
  for (i in (1:nrow(t0))) lr[i] <- (t0[i,1]/ctot[1]) / (t0[i,2]/ctot[2])
  tt <- rbind(c(xb, " "," "," "),cbind(xlb,t0,format(round(lr,dec),nsmall=dec)))
  if (nxL==2) {
    sen <- t0[1,1]/ctot[1]; spe<- t0[2,2]/ctot[2]
    ppv <- t0[1,1]/rtot[1]; npv<- t0[2,2]/rtot[2]
    acc<-(t0[1,1]+t0[2,2])/(ctot[1]+ctot[2])
    plr<-sen/(1-spe); nlr<-(1-sen)/spe;  nnd<-1/(sen+spe-1); dor<-plr/nlr
    tt<-cbind(tt,rbind(format(round(c(sen,spe,ppv,npv,acc,plr,nlr,dor,nnd),dec),nsmall=dec),matrix(" ",nrow=nxL,ncol=9)))
  } else {
    tmp.glm<-glm(d0~factor(x),family=binomial(link="logit"))
    roc0<-roc(d0,predict(tmp.glm))
    a.tss<-roc2xls(roc0,d0)
    colnames(a.tss)<-c("threshold","specificity","sensitivity","accuracy","positive-LLR","negative-LLR","diagnose-OR","N-for-diagnose","postive-pv","negative-pv","a","b","c","d")
    bb<-format(round(a.tss[which.max(a.tss[,2]+a.tss[,3]),c(3,2,9,10,4,5,6,7,8)],dec),nsmall=dec)
    xlsfname<-paste(pngf,"_roc.xls",sep="")
    write.table(a.tss,file=xlsfname,row.names=FALSE,col.names=TRUE,sep="\t",append=FALSE,quote=FALSE)
    png(paste(ofname,"_raw.png",sep="")); rocplot(roc0,xb); dev.off()
    pdf(paste(ofname,"_raw.pdf",sep=""), width=pdfwd, height=pdfht, family="Helvetica"); rocplot(roc0,xb); dev.off()
    tt<-cbind(tt,rbind(bb,matrix(" ",nrow=(nrow(tt)-1),ncol=9)))
  }
  return(tt)
}
roc2xls<-function(roc0,d) {
  a.tss<-t(coords(roc0,"all",ret=c("threshold", "specificity", "sensitivity"), transpose=TRUE))
  bdac<-table(d); t.bd<-bdac[1]; t.ac<-bdac[2]
  t.a<-a.tss[,3]*t.ac; t.c<-t.ac-t.a; t.d<-a.tss[,2]*t.bd; t.b<-t.bd-t.d;
  acc<-(t.a+t.d)/(t.ac+t.bd)
  plr<-(t.a/(t.a+t.c))/(t.b/(t.b+t.d)) 
  nlr<-(t.c/(t.a+t.c))/(t.d/(t.b+t.d))
  dor<-plr/nlr
  nnd<-1/(t.a/(t.a+t.c)-t.b/(t.b+t.d))
  ppv<-t.a/(t.a+t.b)
  npv<-t.d/(t.c+t.d)
  a.tss<-cbind(a.tss,acc,plr,nlr,dor,nnd,ppv,npv,t.a,t.b,t.c,t.d)
  return(a.tss)
}
glm2formula<-function(mdl) {
  coe<-summary(mdl)$coefficients; 
  tmp<-rownames(coe);
  tmp<-gsub(")","=",tmp)
  tmp[substr(tmp,1,7)=="factor("]<-paste(substr(tmp[substr(tmp,1,7)=="factor("],7,99),")",sep="")
  tmp<-paste("*",tmp,sep="")
  tmp[1]<-""; 
  tme<-gsub(" ","",numfmt(coe[,1],5)); 
  tme[as.numeric(tme)>0]<-paste("+",tme[as.numeric(tme)>0],sep="")
  tme<-paste(tme,tmp,sep="",collapse=" ")
  if (substr(tme,1,1)=="+") tme<-substr(tme,2,999);
  return(paste("logit(",colvname,") = ",tme,sep=""))
}
diagrocx1<-function(d,x,xb,pngfname) {
  cmp<- !is.na(d) & !is.na(x); d<-d[cmp]; x<-x[cmp]
  ndL <- length(levels(factor(d)));
  tt<-c(xb,format(c(sum(d),length(d)-sum(d)),nsmall=0))
  if ((min(x)==max(x)) | ndL!=2) return (c(tt,rep(" ",3)))
  roc0<-roc(d,x)
  a.tss<-roc2xls(roc0,d)
  colnames(a.tss)<-c("threshold","specificity","sensitivity","accuracy","positive-LLR","negative-LLR","diagnose-OR","N-for-diagnose","postive-pv","negative-pv","a","b","c","d")
  b.tss<-a.tss[which.max(a.tss[,2]+a.tss[,3]),]
  xlsfname<-paste(pngfname,"_roc.xls",sep="")
  write.table(a.tss,file=xlsfname,row.names=FALSE,col.names=TRUE,sep="\t",append=FALSE,quote=FALSE)
  png(paste(pngfname,"_raw.png",sep="")); rocplot(roc0,xb);  dev.off()
  pdf(paste(pngfname,"_raw.pdf",sep=""), width=pdfwd, height=pdfht, family="Helvetica"); rocplot(roc0,xb); dev.off()
  roc1<-try(roc(d,x,smooth=bt.smooth,ci=TRUE,boot.n=bt.times,ci.alpha=0.95))
  if (length(roc1)>2) {
   if (bt.smooth) {
     png(paste(pngfname,".png",sep="")); rocplot(roc1,xb);  dev.off()   
     pdf(paste(pngfname,".pdf",sep=""), width=pdfwd, height=pdfht, family="Helvetica"); rocplot(roc1,xb); dev.off()
     sens.ci<-try(ci.se(roc1,specificities=seq(0,1,0.05)))
     if (length(sens.ci)>1) {
       png(paste(pngfname,"_ci.png",sep=""))
       plot(roc1,grid=TRUE,legacy.axes=TRUE,main=paste(xb, "\nAUC = ",round(roc1$auc,3)))
       plot(sens.ci,type="shape",col="lightblue")
       plot(sens.ci,type="bars")
       dev.off()
       pdf(paste(pngfname,"_ci.pdf",sep=""), width=pdfwd, height=pdfht, family="Helvetica"); 
       plot(roc1,grid=TRUE,legacy.axes=TRUE,main=paste(xb, "\nAUC = ",round(roc1$auc,3)))
       plot(sens.ci,type="shape",col="lightblue")
       plot(sens.ci,type="bars")
       dev.off()
     }
   }
   aucci<-c(numfmt(c(roc1$auc,roc1$ci[c(1,3)],b.tss[1:10]),dec),b.tss[11:14])
  } else {aucci<-rep(NA,17);}
  return(list(c(tt,aucci),roc0))
}
diagrocx2<-function(d,x1,x1b,x2,x2b,pngfname,k) {
  cmp<- !is.na(d) & !is.na(x1) & !is.na(x2); d<-d[cmp];x1<-x1[cmp];x2<-x2[cmp]
  ndL <- length(levels(factor(d)));
  tt1<-c(x1b,format(c(sum(d),length(d)-sum(d)),nsmall=0))
  tt2<-c(x2b,format(c(sum(d),length(d)-sum(d)),nsmall=0))
  if (ndL!=2) return (rbind(c(tt1,rep(" ",3)),c(tt2,rep(" ",3))))
  tt<-c("Test/model","D+","D-","AUC","AUC.low","AUC.upp");
  cnm<-c("threshold", "specificity", "sensitivity","accuracy","positive-LLR","negative-LLR")
  cnm<-c(cnm,"diagnose-OR","N-for-diagnose","postive-pv","negative-pv","a","b","c","d")
  tt<-c(tt,cnm)
  if (min(x1)!=max(x1)) {
    roc0<-roc(d,x1)
    a.tss1<-roc2xls(roc0,d)
    colnames(a.tss1)<-cnm
    b.tss<-a.tss1[which.max(a.tss1[,2]+a.tss1[,3]),]
    png(paste(pngfname,"_raw1.png",sep="")); rocplot(roc0,x1b); dev.off()  
    pdf(paste(pngfname,"_raw1.pdf",sep=""), width=pdfwd, height=pdfht, family="Helvetica"); rocplot(roc0,x1b); dev.off() 
    roc1<-try(roc(d,x1,smooth=bt.smooth,ci=TRUE,boot.n=bt.times,ci.alpha=0.95))
    if (length(roc1)>2) {
      if (bt.smooth) {
        png(paste(pngfname,"_1.png",sep="")); rocplot(roc1,x1b); dev.off() 
        pdf(paste(pngfname,"_1.pdf",sep=""), width=pdfwd, height=pdfht, family="Helvetica"); rocplot(roc1,x1b); dev.off() 
        sen1.ci<-try(ci.se(roc1,specificities=seq(0,1,0.05)))
        if (length(sen1.ci)>1) {
          png(paste(pngfname,"_ci1.png",sep=""))
          plot(roc1,xlim=c(1.0,0.0),grid=TRUE,legacy.axes=TRUE,main=paste(x1b, "\n AUC = ",round(roc1$auc,3)))
          plot(sen1.ci,type="shape",col="lightblue")
          plot(sen1.ci,type="bars")
          dev.off()
          pdf(paste(pngfname,"_ci1.pdf",sep=""), width=pdfwd, height=pdfht, family="Helvetica");
          plot(roc1,xlim=c(1.0,0.0),grid=TRUE,legacy.axes=TRUE,main=paste(x1b, "\n AUC = ",round(roc1$auc,3)))
          plot(sen1.ci,type="shape",col="lightblue")
          plot(sen1.ci,type="bars")
          dev.off()
        }
      }
      auc1<-c(tt1,numfmt(c(roc1$auc,roc1$ci[c(1,3)],b.tss[1:10]),dec),b.tss[11:14])
    } else {auc1<-rep(NA,length(tt));}
    tt1=cbind(tt,auc1)
  }
  if (min(x2)!=max(x2)) {
    roc02<-roc(d,x2)
    a.tss2<-roc2xls(roc02,d)
    colnames(a.tss2)<-cnm
    b.tss<-a.tss2[which.max(a.tss2[,2]+a.tss2[,3]),]
    png(paste(pngfname,"_raw2.png",sep="")); rocplot(roc02,x2b);  dev.off() 
    pdf(paste(pngfname,"_raw2.pdf",sep=""), width=pdfwd, height=pdfht, family="Helvetica"); rocplot(roc02,x2b);  dev.off() 
    roc2<-try(roc(d,x2,smooth=bt.smooth,ci=TRUE,boot.n=bt.times,ci.alpha=0.95))
    if (length(roc2)>2) {
      if (bt.smooth) {
        png(paste(pngfname,"_2.png",sep="")); rocplot(roc2,x2b); dev.off()      
        pdf(paste(pngfname,"_2.pdf",sep=""), width=pdfwd, height=pdfht, family="Helvetica");rocplot(roc2,x2b); dev.off()      
        sen2.ci<-try(ci.se(roc2,specificities=seq(0,1,0.05)))
        if (length(sen2.ci)>1) {
          png(paste(pngfname,"_ci2.png",sep=""))
          plot.roc(roc2,xlim=c(1.0,0.0),grid=TRUE,legacy.axes=TRUE,main=paste(x2b, "\n AUC = ",round(roc2$auc,3)))
          plot(sen2.ci,type="shape",col="lightblue")
          plot(sen2.ci,type="bars")
          dev.off()
          pdf(paste(pngfname,"_ci2.pdf",sep=""), width=pdfwd, height=pdfht, family="Helvetica");
          plot.roc(roc2,xlim=c(1.0,0.0),grid=TRUE,legacy.axes=TRUE,main=paste(x2b, "\n AUC = ",round(roc2$auc,3)))
          plot(sen2.ci,type="shape",col="lightblue")
          plot(sen2.ci,type="bars")
          dev.off()
        }
      }
      auc2<-c(tt2,numfmt(c(roc2$auc,roc2$ci[c(1,3)],b.tss[1:10]),dec),b.tss[11:14])
    } else {auc2<-rep(NA,length(tt));}
    tt1=cbind(tt1,auc2)
  }
  if ((min(x1)!=max(x1)) & (min(x2)!=max(x2))) {
    if (nrow(a.tss1)==nrow(a.tss2)) {
      a.tss<-cbind(a.tss1,a.tss2)
      colnames(a.tss)<-c(paste(colnames(a.tss1),"1",sep="."), paste(colnames(a.tss2),"2",sep="."))
      write.table(a.tss,file=paste(pngfname,"_roc.xls",sep=""),row.names=FALSE, col.names=TRUE, sep="\t",append=FALSE,quote=FALSE)
    } else {
      write.table(a.tss1,file=paste(pngfname,"_roc1.xls",sep=""),row.names=FALSE, col.names=TRUE, sep="\t",append=FALSE,quote=FALSE)
      write.table(a.tss2,file=paste(pngfname,"_roc2.xls",sep=""),row.names=FALSE, col.names=TRUE, sep="\t",append=FALSE,quote=FALSE)
    }
    if (nbg>1) {tmp.mb<-paste(bvb,bv.lb[k],sep=": ");} else {tmp.mb<-paste("Compare 2 models");}
    png(paste(pngfname,"_2raw.png",sep="")); rocplot2(list(roc0,roc02),c(x1b,x2b),tmp.mb); dev.off()
    pdf(paste(pngfname,"_2raw.pdf",sep=""), width=pdfwd, height=pdfht, family="Helvetica");
    rocplot2(list(roc0,roc02),c(x1b,x2b),tmp.mb); dev.off()
    if (length(roc1)>2 & length(roc2)>2) {
      p<-roc.test(roc1,roc2,reuse.auc=FALSE,boot.n=bt.times)$p.value
      ccp<-c("P(compare)",rep(" ",2),format(round(p,dec),nsmall=dec))
      ccp<-c(ccp,rep(" ",length(tt)-length(ccp)))
      tt1=cbind(tt1,ccp) 
      if (bt.times>0) {
        png(paste(pngfname,"_2smooth.png",sep="")); rocplot2(list(roc1,roc2),c(x1b,x2b),tmp.mb); dev.off()
        pdf(paste(pngfname,"_2smooth.pdf",sep=""), width=pdfwd, height=pdfht, family="Helvetica");
          rocplot2(list(roc1,roc2),c(x1b,x2b),tmp.mb); dev.off()
      }
    } 
  }
  if (!is.matrix(tt1)) {
    return(rep(" ",length(tt)));
  } else {
    if (ncol(tt1)>1) {return(tt1[,-1])} else {return(rep(" ",length(tt)))}
  }
}
plotnomogram <- function(tmp.nom) {
   tmp.nom.x<-tmp.nom[length(tmp.nom)][[1]]$x
   x.mid<-(max(tmp.nom.x)+min(tmp.nom.x))/2
   x.mid.dis<-abs(tmp.nom.x-x.mid)
   x.setblank<-(1:length(tmp.nom.x))[-c(1,which(x.mid.dis==min(x.mid.dis)),length(tmp.nom.x))]
   tmp.nom[length(tmp.nom)][[1]]$fat[x.setblank]<-""
   lbl.nchar<-max(nchar(c(names(tmp.nom),"linear.predictor")))
   lblfac<-0.2
   if (lbl.nchar>20) lblfac<-0.3
   if (lbl.nchar>30) lblfac<-0.4
   plot(tmp.nom,xfrac=lblfac); 
}
mdl2oddsratio<-function(mdlobj, rr="OR") {
  tmps <- summary(mdlobj)$coefficients
  tmpc <- cbind(tmps[,1],tmps[,1]-1.96*tmps[,2],tmps[,1]+1.96*tmps[,2])
  tmpor<- cbind(exp(tmpc),tmps[,4])
  colnames(tmpor) <- c(rr,"Low 95%CI","High 95%CI", "P value")
  rownames(tmpor) <- rownames(tmps)
  return(tmpor)
}
mat2htmltable<-function(mat, rcname=FALSE) {
  if (rcname) mat<-cbind(c("",rownames(mat)),rbind(colnames(mat),mat))
  t1<- apply(mat,1,function(z) paste(z,collapse="</td><td>"))
  t2<- paste("<tr><td>",t1,"</td></tr>")
  return(paste(t2,collapse=" "))
}
vlabelN<-(substr(vlabel,1,1)==" ");
vlabelZ<-vlabel[vlabelN];vlabelV<-vlabel[!vlabelN]
vnameV<-vname[!vlabelN];vnameZ<-vname[vlabelN]
bt.smooth<-chk; bt.times<-500;
if (is.na(par2) || par2=="") {
  study.type <- "cohort";
} else {
  if (par2>=1) par2 <- par2/100
  if (par2>0) {study.type<-"case-control"; population.r <- par2;} else {study.type="cohort";}
}  
if (length(avname)==0 & par1==3) par1==2
if (length(avname)>0) av<-WD[,avname]
if (!is.na(bvar)) {
 bv<-WD[,bvar]
 bvb<-vlabelV[match(bvar,vnameV)]; if (is.na(bvb)) bvb<-bvar; 
 bv.lv<-levels(factor(WD[,bvar])); 
 bv.lb<-vlabelZ[match(paste(bvar,bv.lv,sep="."),vnameZ)]
 bv.lb[is.na(bv.lb)]<-bv.lv[is.na(bv.lb)]
 nbg<-length(bv.lv); 
} else {bv<-rep(1,nrow(WD)); bv.lv<-rep(1,nrow(WD)); nbg<-1; bvb<-"ALL";}
nx<-length(xvname); xb<-vlabelV[match(xvname,vnameV)]; xb[is.na(xb)]<-xvname[is.na(xb)];
xv<-WD[,xvname]
if (nx==1) xv<-matrix(xv,ncol=1)
colv.lv<-levels(factor(WD[,colvname]))
yb<-vlabelV[match(colvname,vnameV)]; if (is.na(yb)) yb<-colvname; 
ylb<-c("0","1")
if (length(colv.lv)==2 & colv.lv[1]==0 & colv.lv[2]==1) {
  yv<-WD[,colvname]; 
} else { 
  if (length(colv.lv)>=2) {
    yv<-(WD[,colvname]==colv.lv[length(colv.lv)])*1; yv[is.na(WD[,colvname])]<-NA;
    ylb<-c("0",paste(colvname,"=",colv.lv[length(colv.lv)],sep=""))
  } else {par1==0;}
}
ofname1<-ofname;
tmp.ss<-rep(NA,13); if (nbg>1) tmp.ss<-rep(NA,14)
tmp.roc<-rep(NA,20)
dec<-4; rocc<-list(NA)
sink(paste(ofname1,".txt",sep=""))
tmp.fml<-NULL
w.00<-""
tmp1<-c("Test","1","0","ROC area(AUC)","95%CI low","95%CI upp","Best threshold","Specificity","Sensitivity")
tmp1<-c(tmp1,"Accuracy","Positive-LR","Negative-LR","Diagnose-OR","N-for-diagnose","Postive-pv","Negative-pv","a","b","c","d")
for (k in (1:nbg)) {
  ncxplot<-0
  if (nbg>1) {
    print(paste(bvb,"=",bv.lb[k]))
    w.00<-c(w.00,"</br>", paste(bvb,"=",bv.lb[k]), "</br></br>")
  }  
  if (par1==1) {
    dcas <- list()
    ndca <- 0
    dcaname <- NULL
    for (i in (1:nx)) {  
	  if (is.na(bvar)) {
	    x.tmp<-xv[,i]; y.tmp<-yv;
	  } else {
	    x.tmp<-xv[bv==bv.lv[k],i]; y.tmp<-yv[bv==bv.lv[k]]
	  }
      cmp<-(!is.na(x.tmp) & !is.na(y.tmp))
      x.tmp<-x.tmp[cmp]; y.tmp<-y.tmp[cmp]; WDTMP <- as.data.frame(cbind(y.tmp,x.tmp))
      singleChk.x<-(max(x.tmp)!=min(x.tmp))
      if (length(levels(factor(y.tmp)))==2  & singleChk.x) {
        pngf<-paste(ofname1,xvname[i],sep="_")
        if (nbg>1) pngf<-paste(pngf,bv.lv[k],sep="_")
        if (xlv[i]==0) {
          tmp.xb<-xb[i]; if (nbg>1) tmp.xb<-paste(tmp.xb," (",bvb,"=",bv.lb[k],")",sep="")
          tmp<-diagrocx1(y.tmp,x.tmp,tmp.xb,pngf)
          tmp.roc<-rbind(tmp.roc,tmp[[1]])
          ncxplot<-ncxplot+1
          rocc[[ncxplot]]<-tmp[[2]]
          if (ncxplot==1) {xxb<-tmp.xb;} else {xxb<-c(xxb,tmp.xb);}
          ndca <- ndca + 1
          dcas[[ndca]] <- decision_curve(y.tmp~x.tmp, data=WDTMP, study.design="cohort", policy="opt-in", bootstraps = 50)
          write.table(dcas[[ndca]]$derived.data,file=paste(pngf,xvname[i],"dca.xls",sep="_"),col.names=TRUE,row.names=FALSE,sep="\t")
          if (is.null(dcaname)) {dcaname <- xvname[i];} else {dcaname<-c(dcaname, xvname[i]);}
        } else {
          xi.lv<-levels(factor(x.tmp))
          tmp.xlb<-vlabelZ[match(paste(xvname[i],xi.lv,sep="."),vnameZ)]
          tmp.xlb[is.na(tmp.xlb)]<-xi.lv[is.na(tmp.xlb)]
          tmp<-diagss(y.tmp,x.tmp,xb[i],tmp.xlb,ylb,pngf)
          if (nbg>1) tmp<-cbind(bv.lb[k],tmp)
          tmp.ss<-rbind(tmp.ss,tmp)
        }
      } 
      rm(WDTMP)
    }
    if (ncxplot>1) {
      pngf<-ofname1; if (nbg>1) pngf<-paste(pngf,bv.lv[k],sep="_")
      if (nbg>1) {tmp.mb<-paste(bvb,bv.lb[k],sep=": ");} else {tmp.mb<-paste("ROC curve for", yb);}
      png(paste(pngf,"_rocs.png",sep=""),width=640,height=640); rocplot2(rocc,xxb,tmp.mb); dev.off()
      pdf(paste(pngf,"_rocs.pdf",sep=""),width=pdfwd, height=pdfht, family="Helvetica");rocplot2(rocc,xxb,tmp.mb); dev.off()
    }
    if (is.matrix(tmp.ss)) {
      tmp0<-c("Test","1","0","Likelihood ratio(LR)","Sensitivity","Specificity","Positive pv","Negative pv","Accuracy","Positive-LR","Negative-LR","Diagnosis OR","Number for diagnose")
      if (nbg>1) tmp0<-c(bvb,tmp0)
      tmp.ss<-rbind(tmp0,tmp.ss[-1,])
      w.00<-c(w.00,"</br>Sensitivity and Specificity for diagnositic test</br><table border=3>",mat2htmltable(tmp.ss),"</table>")
    }
    if (is.matrix(tmp.roc)) {
      tmp.roc<-rbind(tmp1,tmp.roc[-1,])
      tmp.roc1<-tmp.roc[,(1:9)]
      tmp.roc2<-tmp.roc[,-(2:6)]
      w.00<-c(w.00,"</br>ROC analysis for continuous predictor</br><table border=3>",mat2htmltable(tmp.roc1),"</table>")
      w.00<-c(w.00,"</br>Best threshold analysis</br><table border=3>",mat2htmltable(tmp.roc2),"</table>")
    }
    if (ndca>0) {
      png(paste(pngf,"_dca.png",sep=""),width=960,height=840); 
      if (study.type=="cohort") {
        plot_decision_curve(dcas, curve.names = dcaname, col = (2:(ndca+1)), confidence.intervals = FALSE)
      } else {
        plot_decision_curve(dcas, curve.names = dcaname, col = (2:(ndca+1)), confidence.intervals = FALSE, 
        study.design = "case-control", population.prevalence = population.r)
      }  
      dev.off()
      pdf(paste(pngf,"_dca.pdf",sep=""),width=pdfwd*1.5, height=pdfht*1.5, family="Helvetica");
      if (study.type=="cohort") {
        plot_decision_curve(dcas, curve.names = dcaname, col = (2:(ndca+1)), confidence.intervals = FALSE)
      } else {
        plot_decision_curve(dcas, curve.names = dcaname, col = (2:(ndca+1)), confidence.intervals = FALSE, 
        study.design = "case-control", population.prevalence = population.r)
      }  
      dev.off()         
    }
  }
  if (par1==2) {
    WD0<-cbind(yv,xv); colnames(WD0)<-c(colvname,xvname);
    if (nbg>1) WD0<-WD0[bv==bv.lv[k],]
    cmp<-(apply(is.na(WD0),1,sum)==0)
    WD0<-as.data.frame(WD0[cmp,])
    if (length(levels(factor(WD0[,1])))==2) {      
      xv1.1<-xvname
      xv1.1[xlv>2]<-paste("factor(",xv1.1[xlv>2],")",sep="")
      singleChk<- apply(cbind(WD0[,xvname],1),2,function(x) return(max(x,na.rm=TRUE)-min(x,na.rm=TRUE)!=0))
      fml.1<-paste(colvname,"~",paste(xv1.1[singleChk],collapse="+"))
      tmp.glm<-try(glm(formula(fml.1),family=binomial(link="logit"),data=WD0))
      pngf<-paste(ofname1,"model",sep="_")
      if (nbg>1) pngf<-paste(pngf,bv.lv[k],sep="_")
      if (substr(tmp.glm[[1]][1],1,5)!="Error") {
        print(summary(tmp.glm)); print(mdl2oddsratio(tmp.glm))
        tmp.lrm<-try(lrm(formula(fml.1),data=WD0))
        if (substr(tmp.lrm[[1]][1],1,5)!="Error") {
          tmp.dst<-datadist(WD0)
          options(datadist='tmp.dst')
          tmp.nom<-nomogram(tmp.lrm, fun=function(x)1/(1+exp(-x)),funlabel=yb)
          print("Nomogram"); print(tmp.nom)
          png(paste(pngf,"_nom.png",sep=""),width=960,height=840); plotnomogram(tmp.nom); dev.off()
          pdf(paste(pngf,"_nom.pdf",sep=""),width=pdfwd*1.5, height=pdfht*1.5, family="Helvetica");plotnomogram(tmp.nom); dev.off()
        }
        x.tmp<-predict(tmp.glm)
        y.tmp<-WD0[,1]
        tmp.xb<-"Model"; 
        tmpi<-glm2formula(tmp.glm)
        if (nbg>1) {
          tmp.xb<-paste(tmp.xb," (",bvb,"=",bv.lb[k],")",sep="")
          tmpi<-paste("</br>",bvb,"=",bv.lb[k],":</br>",tmpi)
        }
        tmp.fml<-c(tmp.fml,tmpi)
        tmp<-diagrocx1(y.tmp,x.tmp,tmp.xb,pngf)
        tmp.roc<-cbind(tmp.roc,tmp[[1]])
        
        dca0 <- decision_curve(formula(fml.1), data=WD0, study.design="cohort", policy="opt-in", bootstraps = 50)
        png(paste(pngf,"_dca.png",sep=""),width=960,height=840); 
        if (study.type=="cohort") {
          plot_decision_curve(dca0, curve.names="Model",  confidence.intervals = FALSE, col="red")
        } else {
          plot_decision_curve(dca0, curve.names="Model",  confidence.intervals = FALSE, col="red",
          study.design = "case-control", population.prevalence = population.r)
        }  
        dev.off()
        pdf(paste(pngf,"_dca.pdf",sep=""),width=pdfwd*1.5, height=pdfht*1.5, family="Helvetica");
        if (study.type=="cohort") {
          plot_decision_curve(dca0, curve.names="Model",  confidence.intervals = FALSE, col="red")
        } else {
          plot_decision_curve(dca0, curve.names="Model",  confidence.intervals = FALSE, col="red", 
          study.design = "case-control", population.prevalence = population.r)
        }  
        dev.off()
        write.table(dca0$derived.data,file=paste(pngf,"_dca.xls",sep=""),col.names=TRUE,row.names=FALSE,sep="\t")
      } else {tmp.roc<-"Model error!";}
      rm(WD0)
    }
    if (is.matrix(tmp.roc)) {
      tmp.roc<-cbind(tmp1,tmp.roc[,-1])
    } else {w.00<-c(w.00,"</br>",tmp.roc);}
  }
  if (par1==3) {
    WD0<-cbind(yv,xv,av); colnames(WD0)<-c(colvname,xvname,avname)
    if (nbg>1) WD0<-WD0[bv==bv.lv[k],]
    cmp<-(apply(is.na(WD0),1,sum)==0)
    WD0<-as.data.frame(WD0[cmp,])
    if (length(levels(factor(WD0[,1])))==2) {
      xv1<-xvname
      xv1[xlv!=0]<-paste("factor(",xv1[xlv!=0],")",sep="")
      singleChk.xv1 <- apply(cbind(WD0[,xvname],1),2,function(x) return(max(x)-min(x)!=0))
      xv1<- xv1[singleChk.xv1];
      fml<-paste(colvname,"~",paste(xv1,collapse="+"))
      tmp.glm<-try(glm(formula(fml),family=binomial(link="logit"),data=WD0))
      if (substr(tmp.glm[[1]][1],1,5)!="Error") {
	    or.tb1<- numfmt(mdl2oddsratio(tmp.glm),4)
        print(summary(tmp.glm)); print(or.tb1)
        x.tmp<-predict(tmp.glm)
        y.tmp<-WD0[,1]
        av1<-avname
        av1[alv!=0]<-paste("factor(",av1[alv!=0],")",sep="")
        singleChk.av1 <- apply(cbind(WD0[,avname],1),2,function(x) return(max(x)-min(x)!=0))
        av1<- av1[singleChk.av1];
        fml2<-paste(colvname,"~",paste(av1,collapse="+"))
        tmp.glm2<-try(glm(formula(fml2),family=binomial(link="logit"),data=WD0))
        if (substr(tmp.glm2[[1]][1],1,5)!="Error") {
          or.tb2<- numfmt(mdl2oddsratio(tmp.glm2),4)
          print(summary(tmp.glm2)); print(or.tb2)
          x.tmp2<-predict(tmp.glm2)
          pngf<-paste(ofname1,"model2",sep="_")
          if (nbg>1) pngf<-paste(pngf,bv.lv[k],sep="_")
          tmp.x1b<-"Model1"; if (nbg>1) tmp.x1b<-paste(tmp.x1b,"(",bvb,"=",bv.lb[k],")",sep="")
          tmp.x2b<-"Model2"; if (nbg>1) tmp.x2b<-paste(tmp.x2b,"(",bvb,"=",bv.lb[k],")",sep="")
          tmp.fml<-glm2formula(tmp.glm)
          tmp.fml2<-glm2formula(tmp.glm2)
          tmp.roc<-cbind(tmp.roc,diagrocx2(y.tmp,x.tmp,tmp.x1b,x.tmp2,tmp.x2b,pngf,k))

          impr<-improveProb(x.tmp>=as.numeric(tmp.roc[7,2]), x.tmp2>=as.numeric(tmp.roc[7,3]), y.tmp) 
          print(table(x.tmp>=as.numeric(tmp.roc[7,2]), x.tmp2>=as.numeric(tmp.roc[7,3]), y.tmp));
		  print(impr); 
		  w.nrp<-c("Increase for events (1)", "Increase for non-events (2)", "Decrease for events (3)", "Decrease for non-events (4)")
          w.nrp<-cbind(w.nrp, numfmt(unlist(impr[c("pup.ev","pup.ne","pdown.ev","pdown.ne")]), dec))
          nri<-matrix(unlist(impr[c("nri", "nri.ev", "nri.ne", "idi", "se.nri", "se.nri.ev", "se.nri.ne", "se.idi")]),ncol=2)
          nri<-cbind(nri, nri[,1]/nri[,2])
          nri<-cbind(nri, 2-2*pnorm(abs(nri[,3])))
          nri<-cbind(nri, nri[,1]-1.96*nri[,2], nri[,1]+1.96*nri[,2])
          w.nri<-numfmt(nri,dec)
          w.nri[,4]<-pvformat(nri[,4],4)
          w.nri<-rbind(c("Estimate","SE","Z","P value","95% Lower","95% Upper"), w.nri)
          w.nri<-cbind(c("","NRI (1-3+4-2)","NRI for events (1-3)","NRI for non-events (4-2)","IDI"),w.nri)
                   
          dca1 <- decision_curve(formula(fml), data=WD0, study.design="cohort", policy="opt-in", bootstraps = 50)
          dca2 <- decision_curve(formula(fml2), data=WD0, study.design="cohort", policy="opt-in", bootstraps = 50)
          png(paste(pngf,"_dca.png",sep=""),width=960,height=840); 
          if (study.type=="cohort") {
            plot_decision_curve(list(dca1, dca2), curve.names = c("Model 1", "Model 2"), col = c("blue", "red"), confidence.intervals = FALSE)
          } else {
            plot_decision_curve(list(dca1, dca2), curve.names = c("Model 1", "Model 2"), col = c("blue", "red"), confidence.intervals = FALSE,
            study.design = "case-control", population.prevalence = population.r)
          }  
          dev.off()
          pdf(paste(pngf,"_dca.pdf",sep=""),width=pdfwd*1.5, height=pdfht*1.5, family="Helvetica");
          if (study.type=="cohort") {
            plot_decision_curve(list(dca1, dca2), curve.names = c("Model 1", "Model 2"), col = c("blue", "red"), confidence.intervals = FALSE)
          } else {
            plot_decision_curve(list(dca1, dca2), curve.names = c("Model 1", "Model 2"), col = c("blue", "red"), confidence.intervals = FALSE,
            study.design = "case-control", population.prevalence = population.r)
          }  
          dev.off()         
          write.table(dca1$derived.data,file=paste(pngf,"_dca1.xls",sep=""),col.names=TRUE,row.names=FALSE,sep="\t")
          write.table(dca2$derived.data,file=paste(pngf,"_dca2.xls",sep=""),col.names=TRUE,row.names=FALSE,sep="\t")
        } else {tmp.roc<-"Model 2 error!";}
      } else {tmp.roc<-"Model 1 error!";}
      rm(WD0)
    }
    if (is.matrix(tmp.roc)) {
      tmp.roc<-cbind(tmp1,tmp.roc[,-1])
      w.00<-c(w.00,"</br>Predictive model:</br>","Model 1：",tmp.fml,"</br>","Model 2: ",tmp.fml2)
      w.00<-c(w.00,"</br></br>Compare 2 predictive models (ROC Curve)</br><table border=3>",mat2htmltable(tmp.roc),"</table>")
      w.00<-c(w.00,"</br></br>Based on best threshold in the above table:</br>Net Reclassification Index, Integrated Discrimination Index")
	  w.00<-c(w.00,"</br><table border=3>",mat2htmltable(w.nrp),"</table>")	  
	  w.00<-c(w.00,"</br><table border=3>",mat2htmltable(w.nri),"</table>")
      w.00<-c(w.00,"</br>Model 1</br><table border=3>",mat2htmltable(or.tb1,TRUE),"</table>")
      w.00<-c(w.00,"</br>Model 2</br><table border=3>",mat2htmltable(or.tb2,TRUE),"</table>")	  
    } else {w.00<-c(w.00,"</br>",tmp.roc);}
  }
}
sink()
w<-c("<!DOCTYPE html><html lang='zh'><head><meta charset='utf-8'></head><body>")
w<-c(w,"<h2>",title,"</h2>")
 

w<-c(w,"</br>结局变量:", yb, "</br>")
if (par1==2) {
  w <-c(w,"</br>Predictive model:</br>",tmp.fml)
  tmp1<-c("检测项目","1","0","ROC面积(AUC)","95% 区间下限","95% 区间上限","最佳阈值","特异度","敏感度")
  tmp1<-c(tmp1,"准确度","阳性似然比","阴性似然比","诊断比值比","诊断需要检测数","阳性预测值","阴性预测值","a","b","c","d");
  tmp.roc<-cbind(tmp1,tmp.roc[,-1])
  w<-c(w,"</br></br>Model ROC curve and best threshold analysis</br><table border=3>",mat2htmltable(tmp.roc),"</table>")
} else {
  w<-c(w,w.00)
  w<-c(w,"</br>Sensitivity: 敏感度; Specificity: 特异度; Positive pv: 阳性预测值; Negative pv: 阴性预测值; Accuracy: 准确度; Positive-LR: 阳性似然比; Negative-LR: 阴性似然比; Diagnosis OR: 诊断比值比; Number for diagnose: 诊断需要检测数 </br>")
}
if (bt.smooth) w<-c(w,paste("</br>AUC 可信区间与显著性检验采用非参数的重复采样方法(Bootstrap resampling times=",bt.times,")"))
w<-c(w,"</br></br>最佳阈值取敏感度+特异度最大的分界值。各分界点对应的敏感度特异保存在ROC输出文件（.xls）里")
w<-c(w,wd.subset)
w<-c(w,paste("</br></br>Created by EmpowerStats (www.empowerstats.com) and R on",Sys.Date()))
w<-c(w,"</body></html>")
fileConn<-file(paste(ofname1,".htm",sep="")); writeLines(w, fileConn)