## Created by EmpowerStats @ Mon, 15 Sep 25 16:11:41 +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_164_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: ALT.2'); 
ALT.2<-ALT; 
ALT.2[(is.na(ALT))]<-16; 
summary(ALT.2); 
EmpowerStatsR<-cbind(EmpowerStatsR,ALT.2); 
print('Creating new variable: AST.2'); 
AST.2<-AST; 
AST.2[(is.na(AST))]<-21; 
summary(AST.2); 
EmpowerStatsR<-cbind(EmpowerStatsR,AST.2); 
print('Creating new variable: TP.1'); 
TP.1<-TP; 
TP.1[(is.na(TP))]<-73.9; 
summary(TP.1); 
EmpowerStatsR<-cbind(EmpowerStatsR,TP.1); 
print('Creating new variable: ALB.1'); 
ALB.1<-ALB; 
ALB.1[(is.na(ALB))]<-44.5; 
summary(ALB.1); 
EmpowerStatsR<-cbind(EmpowerStatsR,ALB.1); 
print('Creating new variable: GLB.1'); 
GLB.1<-GLB; 
GLB.1[(is.na(GLB))]<-29.3; 
summary(GLB.1); 
EmpowerStatsR<-cbind(EmpowerStatsR,GLB.1); 
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,ALT.2,AST.2,TP.1,ALB.1,GLB.1,CUNBAE); 
sink(); 
vname<-c(NA,'SEX','SEX.0','SEX.1','AGE','NAFLD','NAFLD.0','NAFLD.1','ALT','AST','TP','ALB','GLB','BUN','CR','UA','GLU1','TC1','TG1','HDL1','LDL1','BMI','SBP1','DBP1','MONTH','BMI2','ALT.2','AST.2','TP.1','ALB.1','GLB.1','CUNBAE')[-1]; 
vlabel<-c(NA,'SEX','  0','  1','AGE','NAFLD','  0','  1','ALT','AST','TP','ALB','GLB','BUN','CR','UA','GLU1','TC1','TG1','HDL1','LDL1','BMI','SBP1','DBP1','MONTH','BMI2','ALT','AST','TP','ALB','GLB','CUNBAE')[-1]; 
varused4this <- c('SEX','AGE','NAFLD','ALT','AST','TP','ALB','GLB','BUN','CR','UA','GLU1','TC1','TG1','HDL1','LDL1','BMI','SBP1','DBP1','MONTH','BMI2','ALT.2','AST.2','TP.1','ALB.1','GLB.1','CUNBAE'); 
pkgs<-c('mgcv','gdata','survival','rms'); 
for (g in pkgs) {  
if (!(g %in% rownames(installed.packages()))) install.packages(g,repos='https://cloud.r-project.org'); 
}
library(mgcv); 
library(gdata); 
library(survival); 
library(rms); 
WD <- EmpowerStatsR; rm(EmpowerStatsR); gc(); 
title<-'平滑曲线拟合'; 
wd.subset=''; 
weights.var<- NA; 
yvname<-c('NAFLD'); 
ydist<-c('breslow'); 
ylink<-c('0'); 
ylv<-c(2); 
par1<-1; 
xvname<-c('ALT.2','AST.2','TP.1','ALB.1','GLB.1','BUN','CR','UA','TC1','TG1','HDL1','LDL1','GLU1','SBP1','DBP1'); 
sxf<-c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0); 
xlv<-c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0); 
svname<-c('CUNBAE','BMI'); 
sdf<-c(0,0); 
slv<-c(0,0); 
par3<-2; 
timevar<-'MONTH'; timevlv<- 0; 
vname.start<- NA; 
vname.stop<- NA; 
subjvname<- NA; 
colvname<- NA; 
chk<- 1; 
dec<-4; 
##R package## mgcv gdata survival rms ##R package##;
getNumber<-function(str, n) {
  str<-substr(str,2,nchar(str)-1)
  for (i in (1:nchar(str))) {if (substr(str,i,i)==",") {p=i; break}; }
  ifelse(n==1,return(substr(str,1,p-1)),return(substr(str,p+1,nchar(str))))
}
legLocate<-function(x,y) {
  x[is.infinite(y)]<-NA
  y[is.infinite(y)]<-NA
  xmin<-min(x,na.rm=TRUE); xmax<-max(x,na.rm=TRUE)
  ymin<-min(y,na.rm=TRUE); ymax<-max(y,na.rm=TRUE)
  yoff<-(ymax-ymin); tmp<-table(cut(x,3),cut(y,4))
  tmp.r=which.min(tmp[,4]);tmp.c=4
  if (tmp[2,1]==0) {tmp.r=2;tmp.c=1}
  if (tmp[1,1]==0) {tmp.r=1;tmp.c=1}
  if (tmp[3,1]==0) {tmp.r=3;tmp.c=1}
  if (tmp[2,4]==0) {tmp.r=2;tmp.c=4}
  if (tmp[1,4]==0) {tmp.r=1;tmp.c=4}
  if (tmp[3,4]==0) {tmp.r=3;tmp.c=4}
  pos.y<-colnames(tmp)[tmp.c];  pos.x<-rownames(tmp)[tmp.r];  pct<-0.15
  if (tmp.c==4) {
     if (min(tmp[,4])>0) {pct<-0.3}
     ymax<-ymax+yoff*pct; legy<-ymax;ymin<-ymin-yoff*0.1
  } 
  if (tmp.c==1) {
     if (min(tmp[,1])>0) {pct<-0.3}
     legy<-as.numeric(getNumber(pos.y,2));ymin<-ymin-yoff*pct;ymax=ymax+yoff*0.1
  } 
  legx<-as.numeric(getNumber(pos.x,1))
  return(cbind(xmin,xmax,ymin,ymax,legx,legy))
}
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=" "))
}
cox2htmltable<-function(mdl) {
  ss<-summary(mdl)
  oo1<-cbind(c("N","Number of events"),c(ss$n,ss$nevent))
  coe<-ss$coefficients
  coe<-cbind(coe[,c(1,3,4,5)],coe[,2],1/coe[,2],exp(coe[,1]-1.96*coe[,3]),exp(coe[,1]+1.96*coe[,3]))
  coe<-format(coe,digits=4,nsmall=4)
  cname<-c(" ","coef","se(coef)","z","P.value","exp(coef)","1/exp(coef)","95%CI low","95%CI upp")
  oo2<-rbind(cname,cbind(rownames(coe),coe))
  oo3<-paste(paste(names(ss$rsq),"=",format(ss$rsq,digits=4,nsmall=4)),collapse=", ")
  oo4<-paste("Log Likelihood =",ss$loglik[2])
  oo5<-format(rbind(ss$logtest,ss$waldtest,ss$sctest),digits=4,nsmall=4)
  oo5<-cbind(c(" ","Likelihood ratio test","Wald test","Score(log rank)test"),rbind(colnames(oo5),oo5))
  oo<-c("</br><table border=3>",mat2htmltable(oo1),"</table>")
  oo<-c(oo,"</br>Proportional Harzard Model<table border=3>",mat2htmltable(oo2),"</table>")
  oo<-c(oo,paste("</br>",oo3))
  oo<-c(oo,paste("</br>",oo4))
  oo<-c(oo,"</br><table border=3>",mat2htmltable(oo5),"</table>")
  return(oo)
}
cox2pred<-function(mdl,yi,xi) {
  pred<-predict(mdl,type="terms",se.fit=TRUE)
  mfit<-NA; sfit<-NA;
  if (xi==0) {kb=1; ke=ns;} else {kb=xi; ke=xi;}
  for (k in (kb:ke)) {mfit<-cbind(mfit,pred$fit[,sxx[k]]); sfit<-cbind(sfit,pred$se.fit[,sxx[k]]);}
  tname<-svname[kb:ke];
  if (nx>0) {mfit<-cbind(mfit,pred$fit[,xvv]); sfit<-cbind(sfit,pred$se.fit[,xvv]); tname<-c(tname,xvname);}
  mfit<-matrix(mfit[,-1],ncol=length(tname)); colnames(mfit)<-paste(tname,".fit",sep="");
  sfit<-matrix(sfit[,-1],ncol=length(tname)); colnames(sfit)<-paste(tname,".se",sep="");
  oo<-cbind(wd,mfit,sfit)
  return(oo)
}
cox2surv<-function(mdl,yi,xi,gi=NA) {
    if (!is.na(gi)) {
      svi <-sort(WD[WD[,colvname]==colv.lv[gi],svname[xi]])
      tmpf<-paste(ofname,yvname[yi],svname[xi],colvname,colv.lv[gi],"Surv",sep="_")
    } else {
      svi <-sort(WD[,svname[xi]])
      tmpf<-paste(ofname,yvname[yi],svname[xi],"Surv",sep="_")
    }
    newd <-matrix(rep(refData, length(svi)),ncol=length(refData),byrow=TRUE)
    newd <-cbind(0, predtime, newd)
    colnames(newd) <-c(yvname[yi], timevar, refName)
    newd[,svname[xi]] <- svi
    fv<-predict(mdl,newdata=data.frame(newd),type="survival",se=TRUE) 
    fs<-cbind(fv$fit, fv$fit-2*fv$se.fit, fv$fit+2*fv$se.fit)
    colnames(fs) <- c("Survival", "SurvCI.LOW", "SurvCI.UPP")
    ooi<-cbind(newd,fs) 
    if (!is.na(gi)) {ooi<-cbind(ooi, colv.lv[gi]); colnames(ooi)[ncol(ooi)]<-colvname;}
    ooi<-ooi[,-1]
	ooi <-ooi[!duplicated(ooi[,svname[xi]]),]
	
    png(paste(tmpf,"png",sep="."),width=720,height=560)
    plot(ooi[,svname[xi]], ooi[,"Survival"], type="p", pch=19, ylim=c(0,1), xlab=sb[xi], ylab=paste("Survival at",timevar,":",predtime),main=yb[yi])
    points(ooi[,svname[xi]], ooi[,"SurvCI.LOW"], pch=1)
    points(ooi[,svname[xi]], ooi[,"SurvCI.UPP"], pch=1)
    dev.off()
    pdf(paste(tmpf,"pdf",sep="."),width=pdfwd, height=pdfht, family="GB1");
    plot(ooi[,svname[xi]], ooi[,"Survival"], type="p", pch=19, ylim=c(0,1), xlab=sb[xi], ylab=paste("Survival at",timevar,":",predtime),main=yb[yi])
    points(ooi[,svname[xi]], ooi[,"SurvCI.LOW"], pch=1)
    points(ooi[,svname[xi]], ooi[,"SurvCI.UPP"], pch=1)
    dev.off()
    return(ooi)
}
coxpred2pngs<-function(mat,yi,xi) {
  if (xi==0) {kb=1; ke=ns;} else {kb=xi; ke=xi;}
  px<-c(20,1:9); gg<-"";
  tmp.ylab<-paste("Log RR for",yb.utf8[yi])
  for (k in (kb:ke)) {
     cname1<-paste(svname[k],".fit",sep=""); cname2<-paste(svname[k],".se",sep="")
     mat<-mat[order(mat[,svname[k]]),]
     y.tmp<-mat[,cname1]; se.tmp<-mat[,cname2]; 
     y.low<-y.tmp-1.96*se.tmp; y.upp<-y.tmp+1.96*se.tmp
     x.tmp<-mat[,svname[k]]; 
     pngf<-paste(ofname,yvname[yi],svname[k],"smooth.png",sep="_")
     pdff<-paste(ofname,yvname[yi],svname[k],"smooth.pdf",sep="_")
     pngf0<-paste(ofname,yvname[yi],svname[k],"smooth1.png",sep="_")
     pdff0<-paste(ofname,yvname[yi],svname[k],"smooth1.pdf",sep="_")
     if (is.na(colvname)) {
       if (chk==1) {tmp.col<-rainbow(2);} else {tmp.col<-rep("black",2);}
       xy<-legLocate(c(x.tmp,x.tmp),c(y.low,y.upp))
       png(pngf,width=720,height=560)
       if (slv[k]==0) {
        plot(y.tmp~x.tmp,ylim=c(xy[3],xy[4]),xlim=c(xy[1],xy[2]),col=tmp.col[1],type="p", pch=20, ylab="", xlab="")
        par(new=TRUE); 
        plot(y.low~x.tmp,ylim=c(xy[3],xy[4]),xlim=c(xy[1],xy[2]),col=tmp.col[2],type="p", pch=".", ylab="", xlab="")
        par(new=TRUE); 
        plot(y.upp~x.tmp,ylim=c(xy[3],xy[4]),xlim=c(xy[1],xy[2]),col=tmp.col[2],type="p", pch=".", ylab=tmp.ylab, xlab=sb.utf8[k])
        dev.off()
        png(pngf0,width=720,height=560)
        plot(y.tmp~x.tmp,ylim=c(xy[3],xy[4]),xlim=c(xy[1],xy[2]),col=tmp.col[1],type="l", lty=1, lwd=2, ylab="", xlab="")
        par(new=TRUE); 
        plot(y.low~x.tmp,ylim=c(xy[3],xy[4]),xlim=c(xy[1],xy[2]),col=tmp.col[2],type="l", lty=3, lwd=1, ylab="", xlab="")
        par(new=TRUE); 
        plot(y.upp~x.tmp,ylim=c(xy[3],xy[4]),xlim=c(xy[1],xy[2]),col=tmp.col[2],type="l", lty=3, lwd=1, ylab=tmp.ylab, xlab=sb.utf8[k])
       } else {
        plotCI(x.tmp,y=y.tmp,li=y.low,ui=y.upp,pch=20,lwd=1,col=tmp.col[1],xlim=c(xy[1],xy[2]),ylim=c(xy[3],xy[4]), ylab=tmp.ylab, xlab=sb.utf8[k])
        lines(x.tmp,y.tmp,lty=2)
       }
       dev.off()

       pdf(pdff,width=pdfwd, height=pdfht, family="GB1");
       if (slv[k]==0) {
        plot(y.tmp~x.tmp,ylim=c(xy[3],xy[4]),xlim=c(xy[1],xy[2]),col=tmp.col[1],type="p", pch=20, ylab="", xlab="")
        par(new=TRUE); 
        plot(y.low~x.tmp,ylim=c(xy[3],xy[4]),xlim=c(xy[1],xy[2]),col=tmp.col[2],type="p", pch=".", ylab="", xlab="")
        par(new=TRUE); 
        plot(y.upp~x.tmp,ylim=c(xy[3],xy[4]),xlim=c(xy[1],xy[2]),col=tmp.col[2],type="p", pch=".", ylab=tmp.ylab, xlab=sb.utf8[k])
        dev.off()
        pdf(pdff0,width=pdfwd, height=pdfht, family="GB1");
        plot(y.tmp~x.tmp,ylim=c(xy[3],xy[4]),xlim=c(xy[1],xy[2]),col=tmp.col[1],type="l", lty=1, lwd=2, ylab="", xlab="")
        par(new=TRUE); 
        plot(y.low~x.tmp,ylim=c(xy[3],xy[4]),xlim=c(xy[1],xy[2]),col=tmp.col[2],type="l", lty=3, lwd=1, ylab="", xlab="")
        par(new=TRUE); 
        plot(y.upp~x.tmp,ylim=c(xy[3],xy[4]),xlim=c(xy[1],xy[2]),col=tmp.col[2],type="l", lty=3, lwd=1, ylab=tmp.ylab, xlab=sb.utf8[k])
       } else {
        plotCI(x.tmp,y=y.tmp,li=y.low,ui=y.upp,pch=20,lwd=1,col=tmp.col[1],xlim=c(xy[1],xy[2]),ylim=c(xy[3],xy[4]), ylab=tmp.ylab, xlab=sb.utf8[k])
        lines(x.tmp,y.tmp,lty=2)
       }
       dev.off()

     } else {
       if (chk==1) {tmp.col<-rainbow(ncg);tmp.col1<-c("red","blue")} else {tmp.col<-rep("black",ncg);tmp.col1<-tmp.col}
       for (b in (1:ncg)) {
         y00<-y.tmp[mat[,colvname]==colv.lv[b]]; x00<-x.tmp[mat[,colvname]==colv.lv[b]]
         y01<-y.low[mat[,colvname]==colv.lv[b]]; y02<-y.upp[mat[,colvname]==colv.lv[b]];
         xy1<-legLocate(c(x00,x00),c(y01,y02))
         pngf1<-paste(ofname,yvname[yi],svname[k],colvname,colv.lv[b],"smooth.png",sep="_")
         pdff1<-paste(ofname,yvname[yi],svname[k],colvname,colv.lv[b],"smooth.pdf",sep="_")

         png(pngf1,width=720,height=560)
         if (slv[k]==0) {
          plot(y00~x00,ylim=c(xy1[3],xy1[4]),xlim=c(xy1[1],xy1[2]),col=tmp.col1[1], type="p", pch=20, ylab="", xlab="")
          par(new=TRUE); 
          plot(y01~x00,ylim=c(xy1[3],xy1[4]),xlim=c(xy1[1],xy1[2]),col=tmp.col1[2], type="p", pch=1, ylab="", xlab="")
          par(new=TRUE); 
          plot(y02~x00,ylim=c(xy1[3],xy1[4]),xlim=c(xy1[1],xy1[2]),col=tmp.col1[2], type="p", pch=1, ylab=tmp.ylab, xlab=sb.utf8[k])
         } else {
          plotCI(x00,y=y00,li=y01,ui=y02,pch=20,lwd=1,col=tmp.col[1],xlim=c(xy1[1],xy1[2]),ylim=c(xy1[3],xy1[4]),ylab=tmp.ylab, xlab=sb.utf8[k])
          lines(x00,y00,lty=2)
         }
         dev.off()

         pdf(pdff1,width=pdfwd, height=pdfht, family="GB1");
         if (slv[k]==0) {
          plot(y00~x00,ylim=c(xy1[3],xy1[4]),xlim=c(xy1[1],xy1[2]),col=tmp.col1[1], type="p", pch=20, ylab="", xlab="")
          par(new=TRUE); 
          plot(y01~x00,ylim=c(xy1[3],xy1[4]),xlim=c(xy1[1],xy1[2]),col=tmp.col1[2], type="p", pch=1, ylab="", xlab="")
          par(new=TRUE); 
          plot(y02~x00,ylim=c(xy1[3],xy1[4]),xlim=c(xy1[1],xy1[2]),col=tmp.col1[2], type="p", pch=1, ylab=tmp.ylab, xlab=sb.utf8[k])
         } else {
          plotCI(x00,y=y00,li=y01,ui=y02,pch=20,lwd=1,col=tmp.col[1],xlim=c(xy1[1],xy1[2]),ylim=c(xy1[3],xy1[4]),ylab=tmp.ylab, xlab=sb.utf8[k])
          lines(x00,y00,lty=2)
         }
         dev.off()

       }
       xy<-legLocate(x.tmp,y.tmp)
       if (slv[k]==0) {
        png(pngf0,width=720,height=560)
        for (b in (1:ncg)) {
         y0<-y.tmp[mat[,colvname]==colv.lv[b]]; x0<-x.tmp[mat[,colvname]==colv.lv[b]]
         if (b>1) par(new=TRUE)
         plot(y0~x0,ylim=c(xy[3],xy[4]),xlim=c(xy[1],xy[2]),col=tmp.col[b],type="p", pch=px[b], ylab=tmp.ylab, xlab=sb.utf8[k])
        }
        legend(xy[5],xy[6],colv.lb.utf8,title=colvb.utf8,pch=px[1:ncg],bty="n",col=tmp.col)
        dev.off()

        pdf(pdff0,width=pdfwd, height=pdfht, family="GB1");
        for (b in (1:ncg)) {
         y0<-y.tmp[mat[,colvname]==colv.lv[b]]; x0<-x.tmp[mat[,colvname]==colv.lv[b]]
         if (b>1) par(new=TRUE)
         plot(y0~x0,ylim=c(xy[3],xy[4]),xlim=c(xy[1],xy[2]),col=tmp.col[b],type="p", pch=px[b], ylab=tmp.ylab, xlab=sb.utf8[k])
        }
        legend(xy[5],xy[6],colv.lb.utf8,title=colvb.utf8,pch=px[1:ncg],bty="n",col=tmp.col)
        dev.off()

       }

       png(pngf,width=720,height=560)
       for (b in (1:ncg)) {
         y0<-y.tmp[mat[,colvname]==colv.lv[b]]; x0<-x.tmp[mat[,colvname]==colv.lv[b]]
         tmp.ord<-order(x0); x00<-x0[tmp.ord]; y00<-y0[tmp.ord];
         if (b>1) par(new=TRUE)
         plot(y00~x00,ylim=c(xy[3],xy[4]),xlim=c(xy[1],xy[2]),col=tmp.col[b],type="l", lty=b, lwd=2, ylab=tmp.ylab, xlab=sb.utf8[k])
         rm(tmp.ord,x00,y00)
       }
       legend(xy[5],xy[6],colv.lb.utf8,title=colvb.utf8,lty=(1:ncg),bty="n",col=tmp.col)
       dev.off()

       pdf(pdff,width=pdfwd, height=pdfht, family="GB1");
       for (b in (1:ncg)) {
         y0<-y.tmp[mat[,colvname]==colv.lv[b]]; x0<-x.tmp[mat[,colvname]==colv.lv[b]]
         tmp.ord<-order(x0); x00<-x0[tmp.ord]; y00<-y0[tmp.ord];
         if (b>1) par(new=TRUE)
         plot(y00~x00,ylim=c(xy[3],xy[4]),xlim=c(xy[1],xy[2]),col=tmp.col[b],type="l", lty=b, lwd=2, ylab=tmp.ylab, xlab=sb.utf8[k])
         rm(tmp.ord,x00,y00)
       }
       legend(xy[5],xy[6],colv.lb.utf8,title=colvb.utf8,lty=(1:ncg),bty="n",col=tmp.col)
       dev.off()

     }
     gg<-c(gg,paste("<td>",yb[yi]," vs. ",sb[k],"</br><a href=\"",pngf,"\" target=_BLANK><img src=\"",pngf,"?time=\" width=320,height=320></a></td>",sep=""))
  }
  return(gg)
}
coxsurv2png<-function(ss,yi,xi) {
  if (chk==1) {tmp.col<-rainbow(2);} else {tmp.col<-rep("black",2);}
  y.tmp<- as.numeric(ss[,"Survival"])
  x.tmp<- as.numeric(ss[,svname[xi]])
  xy<-legLocate(x.tmp,y.tmp)
  tmpf<-paste(ofname,yvname[yi],svname[xi],"Survival",sep="_")
  png(paste(tmpf,"png",sep="."),width=720,height=560)
  for (b in 1:ncg) {
    y00 <- y.tmp[ss[,colvname] == colv.lv[b]]
    x00 <- x.tmp[ss[,colvname] == colv.lv[b]]
    if (b>1) par(new=TRUE)
    plot(y00~x00,lty=b,lwd=2,type="l",xlim=c(xy[1],xy[2]),ylim=c(0,1),ylab=paste("Survival at",timevar,":",predtime),xlab=sb.utf8[xi])
    rm(x00,y00)
  }
  legend(xy[5],xy[6],colv.lb.utf8,title=colvb.utf8,lty=(1:ncg),bty="n",col=tmp.col)
  dev.off()
  pdf(paste(tmpf,"pdf",sep="."),width=pdfwd, height=pdfht, family="GB1");
  for (b in 1:ncg) {
    y00 <- y.tmp[ss[,colvname] == colv.lv[b]]
    x00 <- x.tmp[ss[,colvname] == colv.lv[b]]
    if (b>1) par(new=TRUE)
    plot(y00~x00,lty=b,lwd=2,type="l",xlim=c(xy[1],xy[2]),ylim=c(0,1),ylab=paste("Survival at",timevar,":",predtime),xlab=sb.utf8[xi])
    rm(x00,y00)
  }
  legend(xy[5],xy[6],colv.lb.utf8,title=colvb.utf8,lty=(1:ncg),bty="n",col=tmp.col)
  dev.off()
  return(tmpf);
}
if (!is.na(weights.var)) {weights<-WD[,weights.var];} else {weights<-1;}
WD<-cbind(WD,weights);
vlabelN<-(substr(vlabel,1,1)==" ");
vlabelZ<-vlabel[vlabelN];vlabelV<-vlabel[!vlabelN]
vnameV<-vname[!vlabelN];vnameZ<-vname[vlabelN]
ny<-length(yvname); yb<-vlabelV[match(yvname,vnameV)]; yb[is.na(yb)]<-yvname[is.na(yb)]
ns<-length(svname); sb<-vlabelV[match(svname,vnameV)]; sb[is.na(sb)]<-svname[is.na(sb)]
yb.utf8<-yb; Encoding(yb.utf8)<-"UTF-8"
sb.utf8<-sb; Encoding(sb.utf8)<-"UTF-8"
ssf<-rep("",ns); ssf[sdf>0]<-paste(", ",sdf[sdf>0],sep="")
sxx<-paste("rcs(",svname,ssf,")",sep="")
sxx[slv>0]<-paste("factor(", svname[slv>0], ")",sep="")
sxterms<-paste("rcs(",svname,")",sep="")
sxterms[slv>0]<-paste("factor(", svname[slv>0], ")",sep="")
ncg<-1
if (!is.na(colvname)) {
  colv.lv<-levels(factor(WD[,colvname]));  
  ncg<-length(colv.lv);
  colv.lb<-vlabelZ[match(paste(colvname,colv.lv,sep="."),vnameZ)]
  colv.lb[is.na(colv.lb)]<-colv.lv[is.na(colv.lb)]
  colvb<-vlabelV[match(colvname,vnameV)]; if (is.na(colvb)) colvb<-colvname;
  colvb.utf8<-colvb; Encoding(colvb.utf8)<-"UTF-8"; colv.lb.utf8<-colv.lb; Encoding(colv.lb.utf8)<-"UTF-8"
} 
nx<-0
if (length(xvname)>0) {
  if (!is.na(colvname)) {xvname<-xvname[xvname!=colvname];}
  nx<-length(xvname); 
} else {
  xlv<-c();
}
if (nx>0) {
  xb<-vlabelV[match(xvname,vnameV)]; xb[is.na(xb)]<-xvname[is.na(xb)];
  xvv<-xvname; xvv[xlv>2]<-paste("factor(",xvname[xlv>2],")",sep="")
  xv1<-paste(xvv,collapse="+")
}
if (is.na(par1)) par1<-1
if (ny!=ns & par1==2) par1<-1
if (par1==3) {
  sxStr<-paste(sxx,collapse="+"); if (nx>0) sxStr<-paste(sxStr,"+",xv1,sep="")
} else {sxStr=sxx; if (nx>0) sxStr<-paste(sxStr,"+",xv1,sep=""); }
w<-c("<!DOCTYPE html><html lang='zh'><head><meta charset='utf-8'></head><body>")
w<-c(w,paste("<h2>", title, "</h2>"))
wtab<-"</br></br>Cox model</br>"
wpng<-"</br><table>";
cond.logit<-0;
refData = NULL; 
if (!is.na(vname.start)) {
  tmp.xx0<-c(vname.start,timevar); tmptvname<-paste(vname.start,",",timevar);
} else {
  if (is.na(timevar) & !is.na(subjvname)) {
    subj<-WD[,subjvname]; 
    timevar<-"tmp.time"; tmp.time<-rep(1,nrow(WD)); WD<-cbind(WD,tmp.time,subj);
    sxStr<-paste(sxStr,"+strata(subj)");
    ydist<-rep('exact',ny); cond.logit<-1;
  }
  tmp.xx0<-timevar; tmptvname<-timevar;
  if (cond.logit == 0) {
    refData <-c(refData, apply(cbind(1,WD[,c(svname,xvname[xlv==0])]),2,function(z) mean(z, na.rm=TRUE))[-1])
    if (sum(xlv>0)>0) refData <-c(refData, apply(cbind(1,WD[,xvname[xlv>0]]),2,function(z) min(z, na.rm=TRUE))[-1])
    refName <- c(svname, xvname[xlv==0], xvname[xlv>0])
    predtime <-as.numeric(vname.stop); if (is.na(predtime)) predtime <- median(WD[,timevar],na.rm=TRUE)
  }  
}
for (i in (1:ny)) {
  if (par1!=3) {
    wtmp<-"";
    if (par1==2) {jstart<-i; jstop<-i;} else {jstart<-1; jstop<-ns;}
    for (j in (jstart:jstop)) {
      tmp.xx<-c(tmp.xx0,yvname[i],svname[j])
      if (nx>0) tmp.xx<-c(tmp.xx,xvname)
      if (ncg>1) tmp.xx<-c(tmp.xx,colvname)
      if (cond.logit==1) tmp.xx<-c(tmp.xx,subjvname)
      wdd<-WD[,tmp.xx];
      wdd<-wdd[apply(is.na(wdd),1,sum)==0,]
      fml<-paste("Surv(",tmptvname,",",yvname[i],")~",sxStr[j],sep="")
      wtab<-c(wtab,paste("</br></br>Outcome:",yb[i]))
      wtab<-c(wtab,paste("</br>Exposure:",sb[j]))
      pp<-NULL; ss<-NULL;
      if (ncg>1) {
         for (b in (1:ncg)) {
            wd<-wdd[wdd[,colvname]==colv.lv[b],]
            tmp.cox<-coxph(formula(fml),weights=wd$weights,data=wd,method=ydist[i])
            pp<-rbind(pp,cox2pred(tmp.cox,i,j));
            wtab<-c(wtab,paste("</br></br>Sub-group:",colvb,"=",colv.lb[b]))
            wtab<-c(wtab,cox2htmltable(tmp.cox))
            ss<-rbind(ss,cox2surv(tmp.cox,i,j,b));
         }
		 tmpf<-paste(coxsurv2png(ss,i,j),"png",sep=".");
		 wtmp<-c(wtmp,paste("</br>",yb[i]," vs. ",sb[j],"</br><a href=\"",tmpf,"\" target=_BLANK><img src=\"",tmpf,"?time=\" width=320,height=320></a>",sep=""))
      } else {
        wd<-wdd
        tmp.cox<-coxph(formula(fml),weights=wd$weights,data=wd,method=ydist[i])
        pp<-cox2pred(tmp.cox,i,j)
        wtab<-c(wtab,cox2htmltable(tmp.cox))
        ss<-rbind(ss,cox2surv(tmp.cox,i,j));
      }
      xf<-paste(ofname,yvname[i],svname[j],"rcs.xls",sep="_");
      write.table(pp,file=xf,row.names=FALSE,col.names=TRUE,sep="\t",append=FALSE,quote=FALSE)
      wtmp<-c(wtmp,coxpred2pngs(pp,i,j))
      xs<-paste(ofname,yvname[i],svname[j],"surv.xls",sep="_");
      write.table(ss,file=xs,row.names=FALSE,col.names=TRUE,sep="\t",append=FALSE,quote=FALSE)	  
    }
    wpng<-c(wpng,"<tr>",wtmp,"</tr>")
  } else {
      tmp.xx<-c(tmp.xx0,yvname[i],svname);
      if (nx>0) tmp.xx<-c(tmp.xx,xvname)
      if (ncg>1) tmp.xx<-c(tmp.xx,colvname)
      wdd<-WD[,tmp.xx];
      wdd<-wdd[apply(is.na(wdd),1,sum)==0,]
      fml<-paste("Surv(",tmptvname,",",yvname[i],")~",sxStr[1],sep="")
      wtab<-c(wtab,paste("</br></br>Outcome:",yb[i]))
      pp<-NULL; ss<-NULL;
      if (ncg>1) {
         for (b in (1:ncg)) {
            wd<-wdd[wdd[,colvname]==colv.lv[b],]
            tmp.cox<-coxph(formula(fml),weights=wd$weights,data=wd,method=ydist[i])
            pp<-rbind(pp,cox2pred(tmp.cox,i,0));
            wtab<-c(wtab,paste("</br></br>Sub-group:",colvb,"=",colv.lb[b]))
            wtab<-c(wtab,cox2htmltable(tmp.cox))
         }
      } else {
         wd<-wdd
         tmp.cox<-coxph(formula(fml),weights=wd$weights,data=wd,method=ydist[i])
         pp<-cox2pred(tmp.cox,i,0)
         wtab<-c(wtab,cox2htmltable(tmp.cox))
      }
      xf<-paste(ofname,yvname[i],"rcs.xls",sep="_");
      write.table(pp,file=xf,row.names=FALSE,col.names=TRUE,sep="\t",append=FALSE,quote=FALSE)
      wtmp<-coxpred2pngs(pp,i,0)
      wpng<-c(wpng,"<tr>",wtmp,"</tr>")
	  xs<-paste(ofname,yvname[i],"surv.xls",sep="_");
  }
}
wpng<-c(wpng,"</table>")
w<-c(w,wpng,wtab)
w<-c(w,wd.subset)
w<-c(w,paste("</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)