## Created by EmpowerStats @ Mon, 15 Sep 25 16:14:21 +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_165_tbl'; 
rankvar <- function(var, num) {
  qprobs <- 1/num
  if (num>2) {for (i in (2:(num-1))) {qprobs <- c(qprobs, 1/num * i) }} 
  outvar <- rep(0, times=length(var))
  outvar[is.na(var)] <- NA
  cutpoints <- quantile(var,probs=qprobs, na.rm=TRUE)
  for (k in (1:length(cutpoints))) {outvar[var>=cutpoints[k]] <- k}
  tmp<-c(min(var,na.rm=TRUE),cutpoints,max(var,na.rm=TRUE))
  names(tmp)<-c("Min",names(cutpoints),"Max")
  print(tmp)
  return(outvar)
} 
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); 
print('Creating new variable: CUNBAE.Q4'); 
CUNBAE.Q4<- rankvar(CUNBAE,4); 
summary(CUNBAE.Q4); 
EmpowerStatsR<-cbind(EmpowerStatsR,CUNBAE.Q4); 
rm(BMI2,CUNBAE,CUNBAE.Q4); 
sink(); 
vname<-c(NA,'SEX','SEX.0','SEX.1','AGE','NAFLD','NAFLD.0','NAFLD.1','BMI','MONTH','BMI2','CUNBAE','CUNBAE.Q4','CUNBAE.Q4.0','CUNBAE.Q4.1','CUNBAE.Q4.2','CUNBAE.Q4.3')[-1]; 
vlabel<-c(NA,'SEX','  0','  1','AGE','NAFLD','  0','  1','BMI','MONTH','BMI2','CUNBAE','CUNBAE quartile','  Q1','  Q2','  Q3','  Q4')[-1]; 
varused4this <- c('SEX','AGE','NAFLD','BMI','MONTH','BMI2','CUNBAE','CUNBAE.Q4'); 
pkgs<-c('survival','gdata','gridExtra','ggplot2','tidyr','broom','survMisc','survminer','survRM2'); 
for (g in pkgs) {  
if (!(g %in% rownames(installed.packages()))) install.packages(g,repos='https://cloud.r-project.org'); 
}
library(survival); 
library(gdata); 
library(gridExtra); 
library(ggplot2); 
library(tidyr); 
library(broom); 
library(survMisc); 
library(survminer); 
library(survRM2); 
WD <- EmpowerStatsR; rm(EmpowerStatsR); gc(); 
title<-'KM生存曲线-CUNBAE'; 
wd.subset=''; 
svy<- 0; 
weights.var<- NA; 
colvname<-'NAFLD'; colvlv<- 2; 
bvar<-'CUNBAE.Q4'; blv<- 4; 
timevar<-'MONTH'; timevlv<- 0; 
vname.start<- NA; 
par2<- NA; 
par1<- NA; 
par3<-'Auto'; 
avname<-c(); 
dec<-4; 
##R package## survival gdata gridExtra ggplot2 tidyr broom survMisc survminer survRM2 ##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)
}

getLandmarkDF <- function(yvname, WD, timevar, tp) {
  ny<-length(yvname); tp<-c(0, tp, max(WD[,timevar],na.rm=TRUE)); ntp<-length(tp)-1;
  tmpWD<-WD[!is.na(WD[,timevar]),]
  for (i in (1:ntp)) {
    wd.i<-tmpWD[tmpWD[,timevar]>tp[i],]
    for (j in (1:ny)) {wd.i[wd.i[,yvname[j]]==1 & wd.i[,timevar]>tp[i+1],yvname[j]]<-0;}
    wd.i[wd.i[,timevar]>tp[i+1],timevar]<-tp[i+1]
    wd.i<-cbind(wd.i, TIME.PHASE.4LMCOX=i)
    if (i==1) {wd.p<-wd.i;} else {wd.p<-rbind(wd.p,wd.i);}
  }
  return(wd.p)
}
mat2htmltable<-function(mat) {
  t1<- apply(mat,1,function(z) paste(z,collapse="</td><td>"))
  t2<- paste("<tr><td>",t1,"</td></tr>")
  return(paste(t2,collapse=" "))
}

allname<-c(colvname, timevar, vname.start, avname, bvar)
allname<-allname[!is.na(allname)]
WD<-WD[apply(is.na(WD[,allname]),1,sum)==0, allname]
WD[,colvname]<-(WD[,colvname]==1)*1

if (!is.na(par1) & is.na(vname.start)) {
  if (is.numeric(par1)) {tp<-par1;
  } else {
    tmp<-as.numeric(strsplit(par1," ")[[1]]);  tp<-c(tmp[!is.na(tmp)],NA)
  }
  if (length(tp)>2) tp<-tp[1:2]
  tp0<-tp
  WD<-getLandmarkDF(colvname, WD, timevar, tp);
  tmptimeb <- vlabel[vname==timevar];
  tptimeb <- paste(tmptimeb, "<=", tp[1]);
  if (length(tp)>1) {
    tptimeb<-c(tptimeb,  paste(tmptimeb, ">", tp[1], "& <=", tp[2]),paste(tmptimeb, ">", tp[2]));
  } else {
    tptimeb<-c(tptimeb,  paste(tmptimeb, ">", tp[1]));
  }
  ntp<-length(tp)+1; pngfsuff<-paste("_seg",(1:ntp),sep=""); tp<-c(0,tp,max(WD[,timevar],na.rm=TRUE))
} else {
  ntp<-1; WD<-cbind(WD,TIME.PHASE.4LMCOX=1);pngfsuff<-""; tp<-c(0,max(WD[,timevar],na.rm=TRUE))
}

vlabelN<-(substr(vlabel,1,1)==" ");
vlabelZ<-vlabel[vlabelN];vlabelV<-vlabel[!vlabelN]
vnameV<-vname[!vlabelN];vnameZ<-vname[vlabelN]
colvb<-vlabelV[match(colvname,vnameV)]; if (is.na(colvb)) colvb<-colvname;  
colvb.utf8<-colvb; Encoding(colvb.utf8)<-"UTF-8"
fmlx<-"1"; fml1<-"1"
if (!is.na(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)]
  bvb.utf8<-bvb; Encoding(bvb.utf8)<-"UTF-8"; bv.lb.utf8<-bv.lb; Encoding(bv.lb.utf8)<-"UTF-8"
  nbg<-length(bv.lv)
  fmlx<-paste("factor(",bvar,")", sep="");
  fml1<-fmlx;
  wdnew <- matrix(as.numeric(bv.lv),ncol=1); colnames(wdnew)<-bvar; 
  if (length(avname)>0) {    
    av1<-avname; av1[alv>2]<-paste("factor(",avname[alv>2],")",sep="")
    fmlx<-paste(fmlx, paste(av1, collapse="+"), sep="+")
    for (i in 1:length(avname)) {
      if (alv[i] == 0) {wdnew <- cbind(wdnew, mean(WD[, avname[i]]));
	  } else {wdnew <- cbind(wdnew, min(WD[, avname[i]]));}
    }	
    colnames(wdnew)<-c(bvar,avname)
  }
  wdnew<-as.data.frame(wdnew)	
}  
fml0<-"";
if (!is.na(vname.start)) {
  fml0<-paste("Surv(,",vname.start,",",timevar,",",colvname,")~")
  timevb<-paste(vlabelV[match(vname.start,vnameV)],"--",vlabelV[match(timevar,vnameV)]);
  if (is.na(timevb)) timevb<-paste(vname.start,"--",timevar);
} else {
  if (!is.na(timevar)) fml0<-paste("Surv(", timevar, ",", colvname, ")~"); 
  timevb<-vlabelV[match(timevar,vnameV)];
  if (is.na(timevb)) timevb<-timevar;
}
timevb.utf8<-timevb; Encoding(timevb.utf8)<-"UTF-8"
fml<-paste(fml0, fmlx);
fml1<-paste(fml0, fml1);

tau <- NULL; armv<-"";
if (!is.na(par3) & !is.na(bvar)) {
  if (blv == 2) {
    if (max(WD[,bvar])!=1 | min(WD[,bvar]) !=0) {
	   tmp.GROUP<- (WD[,bvar]==max(WD[,bvar]))*1
	   tmp.GROUP[is.na(WD[,bvar])]<-NA
	   WD<-cbind(WD, tmp.GROUP)
	   armv<-"tmp.GROUP" 
	} else {armv<-bvar;}
    if (is.numeric(par3)) tau <- par3 	
  }
}

w<-c("<!DOCTYPE html><html lang='zh'><head><meta charset='utf-8'></head><body>")
w<-c(w,"<h2>Kaplan Meier Survival Curve</h2>")
w<-c(w,paste("</br>Outcome:", colvb))
w<-c(w,paste("</br>Time:", timevb))

surv.mat <- NULL
tval<-NULL; if (!is.na(par2)) tval<-par2;
if (!is.na(colvname) & (fml>"")) {
 sink(paste(ofname,".lst",sep="")) 
 for (k in (1:ntp)) {
   if (ntp>1) {
     w<-c(w,paste("</br></br>Landmark analysis for time segment:", tptimeb[k]), "</br>")
     print(paste("Landmark analysis for time segment:", tptimeb[k]))
   }	 
   pngf<-paste(ofname,pngfsuff[k],"_",c(1:4,0,"b"),".png",sep="")
   pdff<-paste(ofname,pngfsuff[k],"_",c(1:4,0,"b"),".pdf",sep="")
   wdtmp<-WD[WD[,"TIME.PHASE.4LMCOX"]==k,]
   tlim<-c(tp[k], tp[k+1])
   oo.med<-NULL
   if (length(avname)>0 | !is.na(bvar)) {
     tmp.cox <-coxph(formula(fml), data=wdtmp); print(tmp.cox)
     tmp.sfit<-survfit(tmp.cox, newdata=wdnew)
     tmp.med0<-surv_median(survfit(formula(fml1), data=wdtmp))
     tmp.med0<-cbind(tmp.med0[,1], numfmt(tmp.med0[,(2:4)],3))
     tmp.med<-surv_median(tmp.sfit)
     tmp.med<-cbind(tmp.med0[,1], numfmt(tmp.med[,(2:4)],3))
     oo.med<- rbind(colnames(tmp.med), c("Unadjusted","","",""), tmp.med0)
     oo.med<- rbind(oo.med, c("Adjusted","","",""), tmp.med)
	 if (!is.na(bvar)) tmp.sfitx<-survfit(formula(fml1), data=wdtmp)
   } else {
     tmp.sfit<-survfit(formula(fml), data=wdtmp)
   }	   
   print(tmp.sfit)  
   if (is.na(bvar)) {
     ss<-summary(tmp.sfit)
     cname<-c(timevb,"N.Risk","N.Event","N.Censor","Survival","95%CI Low","95%CI Upp")
     oo<-cbind(ss$time,ss$n.risk,ss$n.event,ss$n.censor, numfmt(cbind(ss$surv,ss$lower,ss$upper),4))
     if (is.null(surv.mat)) surv.mat<-cname
     surv.mat<-rbind(surv.mat, oo)
     oo<-rbind(cname,oo)
     oo2<-rbind(names(ss$table),round(ss$table,3))
     tmp.miny<-floor(min(tmp.sfit$lower,na.rm=TRUE)*10)/10

     png(pngf[1],width=720,height=560)
     if (!is.na(par2)) {
       plot(tmp.sfit,ylim=c(tmp.miny,1),xlim=tlim, ylab="Survival",xlab=timevb.utf8,main=colvb.utf8, mark.time=FALSE, xaxt="n")
	   axis(1, at=seq(0,tlim[2], by=par2))
	 } else {plot(tmp.sfit,ylim=c(tmp.miny,1),xlim=tlim, ylab="Survival",xlab=timevb.utf8,main=colvb.utf8, mark.time=FALSE)
     }	 
     dev.off()
     pdf(pdff[1],width=pdfwd, height=pdfht, family="GB1");
     if (!is.na(par2)) {
       plot(tmp.sfit,ylim=c(tmp.miny,1),xlim=tlim, ylab="Survival",xlab=timevb.utf8,main=colvb.utf8, mark.time=FALSE, xaxt="n")
	   axis(1, at=seq(0,tlim[2], by=par2))
	 } else {plot(tmp.sfit,ylim=c(tmp.miny,1),xlim=tlim, ylab="Survival",xlab=timevb.utf8,main=colvb.utf8, mark.time=FALSE)
     }		 
     dev.off()

     png(pngf[2],width=720,height=560)
     if (!is.na(par2)) {
       plot(tmp.sfit,mark.time=FALSE,fun="cumhaz",xlim=tlim,ylab="Cumulative Hazard",xlab=timevb.utf8,main=colvb.utf8, xaxt="n")
	   axis(1, at=seq(0,tlim[2], by=par2))
	 } else {plot(tmp.sfit,mark.time=FALSE,fun="cumhaz",xlim=tlim,ylab="Cumulative Hazard",xlab=timevb.utf8,main=colvb.utf8)
     }	 
     dev.off()
     pdf(pdff[2],width=pdfwd, height=pdfht, family="GB1");
     if (!is.na(par2)) {
       plot(tmp.sfit,mark.time=FALSE,fun="cumhaz",xlim=tlim,ylab="Cumulative Hazard",xlab=timevb.utf8,main=colvb.utf8, xaxt="n")
	   axis(1, at=seq(0,tlim[2], by=par2))
	 } else {plot(tmp.sfit,mark.time=FALSE,fun="cumhaz",xlim=tlim,ylab="Cumulative Hazard",xlab=timevb.utf8,main=colvb.utf8)
     } 
     dev.off()

     png(pngf[3],width=720,height=560)
     if (!is.na(par2)) {
       plot(tmp.sfit,mark.time=FALSE,fun="event",xlim=tlim,ylab="Cumulative events (1-y)",xlab=timevb.utf8,main=colvb.utf8, xaxt="n")
       axis(1, at=seq(0,tlim[2], by=par2))
	 } else {plot(tmp.sfit,mark.time=FALSE,fun="event",xlim=tlim,ylab="Cumulative events (1-y)",xlab=timevb.utf8,main=colvb.utf8)
	 }  
	 dev.off()
     pdf(pdff[3],width=pdfwd, height=pdfht, family="GB1");
     if (!is.na(par2)) {
       plot(tmp.sfit,mark.time=FALSE,fun="event",xlim=tlim,ylab="Cumulative events (1-y)",xlab=timevb.utf8,main=colvb.utf8, xaxt="n")
       axis(1, at=seq(0,tlim[2], by=par2))
	 } else {plot(tmp.sfit,mark.time=FALSE,fun="event",xlim=tlim,ylab="Cumulative events (1-y)",xlab=timevb.utf8,main=colvb.utf8)
	 } 
     dev.off()
     ow<-""
	 
   } else {
     sv<-summary(tmp.sfit)
     surv<-NULL; surv.cname<-NULL
     for (b in 1:nbg) {
       surv.cname<-c(surv.cname, paste(trimws(bv.lb[b]), c("Surv", "95%CI low", "95%CI upp"), sep=": "))
       surv<- cbind(surv, sv$surv[,b], sv$lower[,b], sv$upper[,b])
     }
     surv<-numfmt(surv,4)
     cname<-c(timevb,"N.Risk","N.Event","N.Censor", surv.cname)
     oo<-cbind(sv$time,sv$n.risk,sv$n.event,sv$n.censor, surv)
     if (nbg==2) {
       sdiff    <-sv$surv[,1]-sv$surv[,2]
       sdiff.se <-sqrt(apply(sv$std.err*sv$std.err, 1, sum))
       oo<-cbind(oo,round(sdiff,6),round(sdiff.se,6))
       cname<-c(cname,"Surv.diff","Surv.diff.se")
     }	 
     if (is.null(surv.mat)) surv.mat<-cname
     surv.mat<-rbind(surv.mat, oo)
     oo<-rbind(cname, oo)
     	
     tmp.miny<-floor(min(tmp.sfit$lower,na.rm=TRUE)*10)/10

     png(pngf[1],width=720,height=560)
     if (!is.na(par2)) {
       plot(tmp.sfit,lty=1,col=(2:(nbg+1)), ylim=c(tmp.miny,1.05),xlim=tlim, ylab="Survival", xlab=timevb.utf8, main=colvb.utf8, mark.time=FALSE, xaxt="n")
       axis(1, at=seq(0,tlim[2], by=par2))
     } else {plot(tmp.sfit,lty=1,col=(2:(nbg+1)), ylim=c(tmp.miny,1.05),xlim=tlim, ylab="Survival", xlab=timevb.utf8, main=colvb.utf8, mark.time=FALSE)	 
	 }
     legend("topright",bv.lb.utf8,lty=1,col=(2:(nbg+1)),title=bvb.utf8)
     dev.off()
     pdf(pdff[1],width=pdfwd, height=pdfht, family="GB1");
     if (!is.na(par2)) {
       plot(tmp.sfit,lty=1,col=(2:(nbg+1)), ylim=c(tmp.miny,1.05),xlim=tlim, ylab="Survival", xlab=timevb.utf8, main=colvb.utf8, mark.time=FALSE, xaxt="n")
       axis(1, at=seq(0,tlim[2], by=par2))
     } else {plot(tmp.sfit,lty=1,col=(2:(nbg+1)), ylim=c(tmp.miny,1.05),xlim=tlim, ylab="Survival", xlab=timevb.utf8, main=colvb.utf8, mark.time=FALSE)	 
	 }
     legend("topright",bv.lb.utf8,lty=1,col=(2:(nbg+1)),title=bvb.utf8)
     dev.off()
     
     png(pngf[2],width=720,height=560)
     if (!is.na(par2)) {
	   plot(tmp.sfit,lty=1,col=(2:(nbg+1)),mark.time=FALSE,fun="cumhaz",xlim=tlim, ylab="Cumulative Hazard",xlab=timevb.utf8,main=colvb.utf8,xaxt="n")
       axis(1, at=seq(0,tlim[2], by=par2))
     } else {plot(tmp.sfit,lty=1,col=(2:(nbg+1)),mark.time=FALSE,fun="cumhaz",xlim=tlim, ylab="Cumulative Hazard",xlab=timevb.utf8,main=colvb.utf8)
     }
     legend("topleft",bv.lb.utf8,lty=1,col=(2:(nbg+1)),title=bvb.utf8)
     dev.off()
     pdf(pdff[2],width=pdfwd, height=pdfht, family="GB1");
     if (!is.na(par2)) {
	   plot(tmp.sfit,lty=1,col=(2:(nbg+1)),mark.time=FALSE,fun="cumhaz",xlim=tlim, ylab="Cumulative Hazard",xlab=timevb.utf8,main=colvb.utf8,xaxt="n")
       axis(1, at=seq(0,tlim[2], by=par2))
     } else {plot(tmp.sfit,lty=1,col=(2:(nbg+1)),mark.time=FALSE,fun="cumhaz",xlim=tlim, ylab="Cumulative Hazard",xlab=timevb.utf8,main=colvb.utf8)
     }
     legend("topleft",bv.lb.utf8,lty=1,col=(2:(nbg+1)),title=bvb.utf8)
     dev.off()

     png(pngf[3],width=720,height=560)
     if (!is.na(par2)) {
       plot(tmp.sfit,lty=1,col=(2:(nbg+1)),mark.time=FALSE,fun="event",xlim=tlim, ylab="Cumulative events (1-y)",xlab=timevb.utf8,main=colvb.utf8,xaxt="n")
       axis(1, at=seq(0,tlim[2], by=par2))
     } else {plot(tmp.sfit,lty=1,col=(2:(nbg+1)),mark.time=FALSE,fun="event",xlim=tlim, ylab="Cumulative events (1-y)",xlab=timevb.utf8,main=colvb.utf8)
     }
     legend("topleft",bv.lb.utf8,lty=1,col=(2:(nbg+1)),title=bvb.utf8)
     dev.off()
     pdf(pdff[3],width=pdfwd, height=pdfht, family="GB1");
     if (!is.na(par2)) {
       plot(tmp.sfit,lty=1,col=(2:(nbg+1)),mark.time=FALSE,fun="event",xlim=tlim, ylab="Cumulative events (1-y)",xlab=timevb.utf8,main=colvb.utf8,xaxt="n")
       axis(1, at=seq(0,tlim[2], by=par2))
     } else {plot(tmp.sfit,lty=1,col=(2:(nbg+1)),mark.time=FALSE,fun="event",xlim=tlim, ylab="Cumulative events (1-y)",xlab=timevb.utf8,main=colvb.utf8)
     }
     legend("topleft",bv.lb.utf8,lty=1,col=(2:(nbg+1)),title=bvb.utf8)
     dev.off()
	
     if (length(avname)==0) {	 
       ggsurv0<-ggsurvplot(tmp.sfitx,data=wdtmp,pval=TRUE,conf.int=FALSE,xlim=tlim,xlab=timevb.utf8,legend.title=bvb.utf8,legend.labs=bv.lb.utf8,break.time.by=tval,
       risk.table=TRUE,risk.table.y.text.col=T,risk.table.height=0.25,risk.table.y.text=FALSE,ncensor.plot=TRUE,ncensor.plot.height = 0.25)
       png(pngf[5],width=720,height=720); print(ggsurv0);  dev.off()
       pdf(pdff[5],width=pdfwd, height=pdfht*1.2, family="GB1"); print(ggsurv0); dev.off()
    
       ggsurv<-ggsurvplot(tmp.sfitx,data=wdtmp,pval=TRUE,conf.int=FALSE,xlim=tlim,xlab=timevb.utf8,legend.title=bvb.utf8,legend.labs=bv.lb.utf8,break.time.by=tval,
       risk.table=TRUE,risk.table.height=0.25,linetype="strata",color="black")
       png(pngf[6],width=720,height=720); print(ggsurv);  dev.off()
       pdf(pdff[6],width=pdfwd, height=pdfht*1.2, family="GB1"); print(ggsurv); dev.off()
    
       svd<-survdiff(formula(fml), data=wdtmp)
       ow<-c("Log rank test: implements the G-rho family of Harrington and Fleming (1982), with weights on each death of S(t)^rho, where S is the Kaplan-Meier estimate of survival. With rho = 0 this is the log-rank or Mantel-Haenszel test.")
       tmp.svdmat<-rbind(c(" ","N","Observed","Expected"),cbind(names(svd$n),svd$n,svd$obs,numfmt(svd$exp,4)))
       ow<-c(ow,"</br></br><table border=3>",mat2htmltable(tmp.svdmat),"</table></br>") 
       ow<-c(ow,"Chisq=",numfmt(svd$chisq,4), " on ", length(svd$n)-1, " degree of freedom, p=", pvformat(1-pchisq(svd$chisq,length(svd$n)-1),4))
     } else {
       ss<-summary(tmp.cox);  coe<-ss$coefficients;
       coe0<-cbind(coe[,1],coe[,3],coe[,4],coe[,2],1/coe[,2], exp(coe[,1]-1.96*coe[,3]), exp(coe[,1]+1.96*coe[,3]))
       coe<-cbind(numfmt(coe0,4),pvformat(coe[,5],4))
       ow<-c(" ","coef.","se(coef.)","z","HR=exp(coef)","1/HR","95%CI lower","95%CI upper","P　value")
       ow<-rbind(ow,cbind(rownames(coe),coe)) 
	   ow<-paste("</br><table border=3>",mat2htmltable(ow),"</table></br>") 
     }
     if (armv > "") {
       covxx<-NULL; if (length(avname)>0) covxx <- wdtmp[,avname]
	   tauk <- tau;
	   if (is.numeric(tauk)) tauk <- min(max(wdtmp[,timevar],na.rm=TRUE), tauk)
       rs0<-rmst2(wdtmp[,timevar], wdtmp[,colvname], wdtmp[,armv], tau=tauk, covariates=NULL)
       print(rs0)
       t1<-c("Group","<strong>RMST</strong>","se","95%CI lower","95%CI upper","<strong>RMTL</strong>","se","95%CI lower","95%CI upper")
       t1<-rbind(t1,c("arm=1", numfmt(rs0$RMST.arm1$result[1,],3), numfmt(rs0$RMST.arm1$result[2,],3)))
       t1<-rbind(t1,c("arm=0", numfmt(rs0$RMST.arm0$result[1,],3), numfmt(rs0$RMST.arm0$result[2,],3)))
       t2<-rs0$unadjusted.result
       t2[,(1:3)]<-numfmt(t2[,(1:3)],3)
       t2[,4]<-pvformat(t2[,4],4)
       t3<-NULL
       if (!is.null(covxx)) {
         rs1<-rmst2(wdtmp[,timevar], wdtmp[,colvname], wdtmp[,armv], tau=tauk, covariates=covxx)
         print(rs1)
         ta<-rs1$adjusted.result
         ta<-paste(numfmt(ta[,1],3), " (", numfmt(ta[,2],3), ",", numfmt(ta[,3],3), ") ", pvformat(ta[,4],4), sep="")
         t2<-cbind(rownames(t2), paste(t2[,1], " (", t2[,2], ",", t2[,3], ") ", t2[,4], sep=""), ta)
         t2<-rbind(c("","Unadjusted","Adjusted"), t2)
         tmp1<-rs1$RMST.difference.adjusted
         tmp2<-rs1$RMST.ratio.adjusted
         tmp3<-rs1$RMTL.ratio.adjusted
         tmp1<-cbind(rownames(tmp1),numfmt(tmp1[,(1:3)],3),pvformat(tmp1[,4],4),"",numfmt(tmp1[,(5:6)],3))
         tmp2<-cbind(rownames(tmp2),numfmt(tmp2[,(1:3)],3),pvformat(tmp2[,4],4),numfmt(tmp2[,(5:7)],3))
         tmp3<-cbind(rownames(tmp3),numfmt(tmp3[,(1:3)],3),pvformat(tmp3[,4],4),numfmt(tmp3[,(5:7)],3))
         t3<-c("","coef","se","z","P value","exp(coef)","95%CI lower","95%CI upper")
         t3<-rbind(t3,c("Model: difference of RMST", rep("",7)),tmp1)
         t3<-rbind(t3,c("Model: ratio of RMST", rep("",7)),tmp2)
         t3<-rbind(t3,c("Model: ratio of RMTL", rep("",7)),tmp3)
       } else {	 
         t2<-cbind(c("",rownames(t2)), rbind(colnames(t2),t2))
       }	 
       rw<-c("</br><strong>Restricted mean survival time (RMST) analysis</strong></br>")
       rw<-c(rw,"Restricted to time point (tau=)", round(rs0$tau,3), "</br>")
       rw<-c(rw, paste("arm = 1 refer to", bvar, "=", max(WD[,bvar])), "</br>")
       rw<-c(rw,"</br>Unadjusted RMST (restricted mean survival time) and RMTL (restricted mean time lost)")
       rw<-c(rw,"<table border=3>",mat2htmltable(t1),"</table>")
       rw<-c(rw,"</br></br>Between group comparison<table border=3>",mat2htmltable(t2),"</table>")
       if (!is.null(t3)) rw<-c(rw,"</br></br>Model summary<table border=3>",mat2htmltable(t3),"</table>")
     }
  }
  w<-c(w,ow)     
  if (armv >"") w<-c(w,rw)
  w<-c(w,"</br>Survival table<table border=3>",mat2htmltable(oo),"</table></br>") 
  if (is.matrix(oo.med)) {w<-c(w,"</br>Median survival time<table border=3>",mat2htmltable(oo.med),"</table></br>");} 
  if (!is.null(surv.mat)) write.table(surv.mat, file=paste(ofname,".xls",sep=""), col.names=FALSE, row.names=FALSE, sep="\t");
 }
 if (!is.na(bvar) & ntp>1 & length(avname)==0) {
   WD<-cbind(WD, tgrp= WD[,"TIME.PHASE.4LMCOX"]*10+ WD[,bvar])
   fml2<-paste(fml0,"tgrp")
   tmp.sfit<-survfit(formula(fml2), data=WD)
   print(tmp.sfit)
   png(paste(ofname,".png",sep=""),width=720,height=560)
   plot(tmp.sfit,lty=(1:nbg),col=(2:(nbg+1)), ylim=c(tmp.miny,1.05), ylab="Survival", xlab=timevb.utf8, main=colvb.utf8)
   legend("topright",bv.lb.utf8 ,lty=(1:nbg),col=(2:(nbg+1)),title=bvb.utf8)
   abline(v=tp0, lty=1)
   dev.off()
   pdf(paste(ofname,".pdf",sep=""),width=pdfwd, height=pdfht, family="GB1");
   plot(tmp.sfit,lty=(1:nbg),col=(2:(nbg+1)), ylim=c(tmp.miny,1.05), ylab="Survival", xlab=timevb.utf8, main=colvb.utf8)
   legend("topright",bv.lb.utf8 ,lty=(1:nbg),col=(2:(nbg+1)),title=bvb.utf8)
   abline(v=tp0, lty=1)
   dev.off()

   fml3<-paste(fml0,bvar)
   tmp.sfit<-survfit(formula(fml3), data=WD)
   print(tmp.sfit)
   ggsurv<-ggsurvplot(tmp.sfit,data=WD,conf.int=TRUE,xlab=timevb.utf8,risk.table=TRUE,ncensor.plot=TRUE,color="black",break.time.by=tval)
   png(paste(ofname,"_0.png",sep=""),width=720,height=720); print(ggsurv);  dev.off()
   pdf(paste(ofname,"_0.pdf",sep=""),width=pdfwd, height=pdfht*1.5, family="GB1"); print(ggsurv); dev.off()  
 }
 sink()
}
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(ofname,".htm",sep="")); writeLines(w, fileConn)
dev.off()