## Created by EmpowerStats @ Wed, 17 Sep 25 19:08: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_215_tbl'; 
g2idx<-function(x) {
 g<-levels(factor(x));  y<-rep(NA,length(x))
 for (i in 1:length(g)) y[x==g[i]]<-i;
 y[is.na(y)]<-length(g)+1;
 return(y);
}
zscorevname<-function(xvname, gvname, mu, std, ref, stratum=c()) {
 tmpd<-EmpowerStatsR[,xvname]; tmpg_<-rep(1,length(tmpd));
 ifelse(gvname>"", tmpd<-cbind(tmpd,EmpowerStatsR[,gvname]), tmpd<-cbind(tmpd,tmpg_));
 ng<-length(levels(factor(tmpd[,2])))
 if (ref>"") {
   tmpd0<-EmpowerStatsR[eval(parse(text=ref)),xvname]; tmpg0_<-rep(1,length(tmpd0));
   ifelse(gvname>"", tmpd0<-cbind(tmpd0,EmpowerStatsR[eval(parse(text=ref)),gvname]), tmpd0<-cbind(tmpd0,tmpg0_));
 } else { tmpd0<-tmpd; }
 ng0<-length(levels(factor(tmpd0[,2])))
 if (ng!=ng0) tmpd0<-tmpd
 tmpd[,2]<-g2idx(tmpd[,2]); tmpd0[,2]<-g2idx(tmpd0[,2])
 nmu<-length(mu)
 if (nmu==0 | nmu!=length(std) | !is.numeric(mu) | !is.numeric(std) | nmu!=ng) mu<-c()
 if (length(stratum)>0) {
   if (length(mu) != length(stratum)) mu <- c()
 }
 if (length(mu)==0) {
   mu<-c(tapply(tmpd0[,1],tmpd0[,2], mean, na.rm=TRUE))
   std<-c(tapply(tmpd0[,1],tmpd0[,2], sd, na.rm=TRUE))
 } else {
   if (length(stratum)>0) {
      ord <- match(stratum, as.numeric(levels(factor(tmpd0[,2]))))
      mu <- mu[ord]
      std <- std[ord]
   }
 }
 print("Mean:"); print(mu); print("SD:"); print(std);
 y<-(tmpd[,1]-mu[tmpd[,2]])/std[tmpd[,2]];
 return(y);
} 
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); 
print('Creating new variable: CUNBAE.Z'); 
CUNBAE.Z<- zscorevname("CUNBAE","",c(),c(),"",c()); 
summary(CUNBAE.Z); 
EmpowerStatsR<-cbind(EmpowerStatsR,CUNBAE.Z); 
rm(BMI2,ALT.2,AST.2,TP.1,ALB.1,GLB.1,CUNBAE,CUNBAE.Z); 
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','CUNBAE.Z')[-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','CUNBAE Z-score')[-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','CUNBAE.Z'); 
pkgs<-c('survival','gdata'); 
for (g in pkgs) {  
if (!(g %in% rownames(installed.packages()))) install.packages(g,repos='https://cloud.r-project.org'); 
}
library(survival); 
library(gdata); 
WD <- EmpowerStatsR; rm(EmpowerStatsR); gc(); 
title<-'多个回归方程-CUNBAE45'; 
wd.subset=''; 
weights.var<- NA; 
yvname<-c('NAFLD'); 
ydist<-c('breslow'); 
ylink<-c('0'); 
ylv<-c(2); 
avname<-c('ALT.2','AST.2','TP.1','ALB.1','GLB.1','BUN','CR','UA','TC1','TG1','HDL1','LDL1'); 
saf<-c(0,0,0,0,0,0,0,0,0,0,0,0); 
alv<-c(0,0,0,0,0,0,0,0,0,0,0,0); 
svname<-c('ALT.2','AST.2','TP.1','ALB.1','GLB.1','BUN','CR','UA','TC1','TG1','HDL1','LDL1','GLU1','SBP1','DBP1'); 
sdf<-c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0); 
slv<-c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0); 
prn<-1; 
xvname<-c('CUNBAE.Z'); 
sxf<-c(3); 
xlv<-c(0); 
par1<-1; 
chk<- 1; 
cox<- 1; 
timevar<-'MONTH'; timevlv<- 0; 
vname.start<- NA; 
subjvname<- NA; 
gee.TYPE<-NA; 
bvar<- NA; 
colvname<- NA; 
par3<-1; 
dec<-1; 
##R package## survival gdata ##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)
}
varstats<-function(var,vlvl,dec) {
  if (length(vlvl)==1 & vlvl[1]==0) {
    return(paste(numfmt(mean(var,na.rm=TRUE),dec),numfmt(sd(var,na.rm=TRUE),dec),sep="+"))
  } else {
    a<-table(var)
    b<-matrix(paste(a, " (", numfmt(a/sum(a)*100,dec), "%)",sep=""),ncol=1)
    return(c(" ",b[match(vlvl,names(a))]))
  }
}
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=" "))
}
mdl2oo<-function(i, j, m, k, b, wdtmp0) {
  if (par1==4) {fmlx<-paste(xvname_[j],fmlm[m],fmlm4[[m]][j]);} else {fmlx<-paste(xvname_[j],fmlm[m]);}
  if (b==nblv & !is.na(bvar)) fmlx<-paste(fmlx,"+factor(",bvar,")",sep="")
  if (k==nclv & !is.na(colvname)) fmlx<-paste(fmlx,"+factor(",colvname,")",sep="")
  wdtmp<-removeNA(i,j,m,wdtmp0)
  xxname<-xxname_[[j]]
  if (clogit==1) {
    fml<-paste(yvname[i],"~",fmlx,"+strata(",subjvname,")",sep="")
    mdl<-try(clogit(formula(fml),weights=wdtmp$weights,data=wdtmp,na.action=na.omit,method="exact"))
  } else {
    fml<-paste(fmlt,yvname[i],")~",fmlx);
    mdl<-try(coxph(formula(fml),weights=wdtmp$weights,data=wdtmp,na.action=na.omit,method=ydist[i]))
  }
  if (substr(mdl[1],1,5)!="Error") {
    gs<-summary(mdl); print(mdl$formula); print(gs)
    tmpn<-gs$n; gsparm <- gs$coefficients
    gsp<-gsparm[match(xxname,rownames(gsparm)),]
    if (length(xxname)==1) {beta<-gsp[1]; se<-gsp[3]; pv<-gsp[5]; 
    } else {beta<-gsp[,1]; se<-gsp[,3]; pv<-gsp[,5]; }
    ci1<- beta-1.96*se; ci2<- beta+1.96*se
    pvx<-substr(rep("****",length(pv)),1,(pv<=0.05)+(pv<=0.01)+(pv<=0.001)) 
    if (colprn==3) {pvv<-pvx;} else {pvv<-pvformat(pv,dec+2);}
    o1<-paste(numfmt(exp(beta),dec)," (",numfmt(exp(ci1),dec),", ",numfmt(exp(ci2),dec),")",sep="")
    o1<-paste(o1,pvv); o1[is.na(beta)]<-NA;
  } else {o1<-rep(NA,length(xxname)); tmpn<-nrow(wdtmp);}
  if (length(xxname)>1) {
    o1[is.na(o1) & substr(xxname,1,7)=="factor("]<-"1.0"
    o1[is.na(o1)]<-"";
  }
  return(list(o1,tmpn))
}
recodevar <- function (var,oldcode,newcode) {
  tmp.v <- var
  nc.tmp <- length(oldcode)
  for (i in (1:nc.tmp)) {tmp.v[(var==oldcode[i])]=newcode[i]}
  if (is.factor(tmp.v)) {tmp.v1<-as.numeric(as.character(tmp.v))} else {tmp.v1<-as.numeric(tmp.v)}
  rm(tmp.v);  return(tmp.v1)
} 
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; }
  return(outvar)
}
removeNA<-function(i,j,m,wdf) {
 vvv<-c(yvname[i],adjvv[[m]],subjvname,colvname,bvar,vname.start,timevar); 
 if (j<=nx) {vvv<-c(vvv,xvname[j]);} else {vvv<-c(vvv,xvname);}
 vvv<-vvv[!is.na(vvv)]; vvv<-vvv[vvv>" "]
 tmp<-is.na(wdf[,vvv]); 
 return(wdf[apply(tmp,1,sum)==0,])
}
if (!is.na(bvar)) {
 if (bvar=="TIME.PHASE.4LMCOX") {
  WD<-getLandmarkDF(yvname, WD, timevar, tp);
  tmptimeb <- vlabel[vname==timevar];
  vname<-c(vname,"TIME.PHASE.4LMCOX", "TIME.PHASE.4LMCOX.1");
  vlabel<-c(vlabel, paste(tmptimeb, "segment"), paste("  <=", tp[1]));
  if (length(tp)>1) {
    vname<-c(vname, "TIME.PHASE.4LMCOX.2", "TIME.PHASE.4LMCOX.3");
    vlabel<-c(vlabel,  paste("  >", tp[1], " & <=", tp[2],sep=""),paste("  >", tp[2]));
  } else {
    vname<-c(vname, "TIME.PHASE.4LMCOX.2");
    vlabel<-c(vlabel,  paste("  >", tp[1]));
  }
  bv <- WD[,"TIME.PHASE.4LMCOX"];
}}
if (!is.na(colvname)) {
 if (colvname=="TIME.PHASE.4LMCOX") {
  WD<-getLandmarkDF(yvname, WD, timevar, tp);
  tmptimeb <- vlabel[vname==timevar];
  vname<-c(vname,"TIME.PHASE.4LMCOX", "TIME.PHASE.4LMCOX.1");
  vlabel<-c(vlabel, paste(tmptimeb, "segment"), paste("  <=", tp[1]));
  if (length(tp)>1) {
    vname<-c(vname, "TIME.PHASE.4LMCOX.2", "TIME.PHASE.4LMCOX.3");
    vlabel<-c(vlabel,  paste("  >", tp[1], " & <=", tp[2],sep=""),paste("  >", tp[2]));
  } else {
    vname<-c(vname, "TIME.PHASE.4LMCOX.2");
    vlabel<-c(vlabel,  paste("  >", tp[1]));
  }
  colv <- WD[,"TIME.PHASE.4LMCOX"];
}}
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];
w<-c("<!DOCTYPE html><html lang='zh'><head><meta charset='utf-8'></head><body>")
w<-c(w,paste("<h2>", title, "</h2>"))
if (length(avname)>0) { 
 if (sum((saf=="s" | saf=="S") & alv>0)>0) w<-c(w,"</br>Restrict cubic spline smoothing only applies for continuous variables")
}
if (length(svname)>0) {
 if (sum((sdf=="s" | sdf=="S") & slv>0)>0) w<-c(w,"</br>Restrict cubic spline smoothing only applies for continuous variables")
}
allvname<-c(yvname,xvname,colvname,bvar,avname,svname,subjvname,vname.start,timevar,"weights"); 
allvname<-allvname[!is.na(allvname)]
WD<-WD[,allvname];
if (!is.na(subjvname)) WD<-WD[order(WD[,subjvname]),]
sxf<-as.numeric(sxf);sxf[is.na(sxf)]<-0;
if (sum(sxf>1 & xlv>0)>0) w<-c(w,"Categorizing only applies to continuous variables");
if (sum(sxf>1 & xlv==0)>0) {
     t.xname<-NA;t.xlv<-NA; nx<-length(xvname)
     for (i in 1:nx) {
       if (sxf[i]>1 & xlv[i]==0) {
          tmp.Xi<- rankvar(WD[,xvname[i]],sxf[i])
          tmp.newcode <- tapply(WD[,xvname[i]],tmp.Xi,function(z) median(z,na.rm=TRUE))
          tmp.low <- tapply(WD[,xvname[i]],tmp.Xi,function(z) min(z,na.rm=TRUE))
          tmp.upp <- tapply(WD[,xvname[i]],tmp.Xi,function(z) max(z,na.rm=TRUE))
          tmp.Xi2<- recodevar(tmp.Xi,(1:sxf[i])-1,tmp.newcode)
          tmp.Xi<-cbind(tmp.Xi,tmp.Xi2)
          tmp.NM<-paste(xvname[i],c("grp","grp.cont"),sep=".")
          colnames(tmp.Xi)<-tmp.NM
          WD<-cbind(WD,tmp.Xi)
          t.xname<-c(t.xname,tmp.NM)
          t.xlv<-c(t.xlv,sxf[i],0)
          vnameV<-c(vnameV,tmp.NM)
          vlabelV<-c(vlabelV,paste(vlabelV[vnameV==xvname[i]],c("group","group trend")))
          vnameZ<-c(vnameZ,paste(tmp.NM[1],(1:sxf[i])-1,sep="."))
          vlabelZ<-c(vlabelZ,paste(tmp.low,"-",tmp.upp))
       } else {
          t.xname<-c(t.xname,xvname[i]); t.xlv<-c(t.xlv,xlv[i])
       }
     }
     xvname<-t.xname[-1]; xlv<-t.xlv[-1];
}
if (!is.na(subjvname)) {
 if (length(avname)>0) saf<-rep(0,length(saf)); 
 if (length(svname)>0) sdf<-rep(0,length(sdf)); 
 WD<-WD[order(WD[,subjvname]),];
}
ny=length(yvname); nx=length(xvname); 
xb<-vlabelV[match(xvname,vnameV)]; xb[is.na(xb)]<-xvname[is.na(xb)]
yb<-vlabelV[match(yvname,vnameV)]; yb[is.na(yb)]<-yvname[is.na(yb)]
xvname_ <- xvname
xvname_[xlv>0]<-paste("factor(",xvname[xlv>0],")",sep="")
xxname_<-list(NA); xxlbl_<-list(NA); xxlvl_<-list(NA)
for (j in (1:nx)) {
  if (xlv[j]==0) {
    xxname_[[j+1]]<-xvname[j];xxlbl_[[j+1]]<-xb[j];xxlvl_[[j+1]]<-0
  } else {
    xxlvl_[[j+1]]<-levels(factor(WD[,xvname[j]]))
    tmp<-paste(xvname[j],".",xxlvl_[[j+1]],sep="")
    xxlbl_[[j+1]]<-c(xb[j],vlabelZ[match(tmp,vnameZ)])
    xxlbl_[[j+1]]<-paste(c("",rep("&nbsp&nbsp",length(xxlbl_[[j+1]])-1)),xxlbl_[[j+1]])
    xxname_[[j+1]]<-c(xvname[j],paste("factor(",xvname[j],")",xxlvl_[[j+1]],sep=""))
  }
}
xxname_<-xxname_[-1]; xxlbl_<-xxlbl_[-1]; xxlvl_<-xxlvl_[-1];
if (nx==1) par1<-1; 
if (is.na(par1) | par1<1) par1<-1;
if (par1>1) {
  tmp1<-xxname_[[1]]; tmp2<-xxlbl_[[1]]
  for (j in 2:nx) {tmp1<-c(tmp1,xxname_[[j]]); tmp2<-c(tmp2,xxlbl_[[j]]);}
  xxname_[[nx+1]]<-tmp1; xxlbl_[[nx+1]]<-tmp2;
  xvname_<-c(xvname_,paste(xvname_,collapse="+"))
}
contx<-(sum(xlv>0)==0)
if (par1==3 & !is.na(bvar)) {w<-c(w,"</br>Column stratified variable was ignored"); bvar<-NA; bvname<-NA;}

fmlm<-" "; fmlb<-"Non-adjusted"; tmp<-""; adjvv<-list(NA); adjvb<-"None"; fmlp<-"glm";
fmlm4<-list(rep("",nx))
na=0; avb=""; smoothav<-0; nadjm<-0
if (length(avname)>0) { 
  na<-length(avname)
  avb<-vlabelV[match(avname,vnameV)];
  avname_ <- avname 
  smoothavi<-((saf=="s" | saf=="S") & alv==0)
  smoothav<-sum(smoothavi)
  smoothavname<-avname[smoothavi]
  avname_[smoothavi]<-paste("rcs(",avname[smoothavi],")",sep="")
  avb1<-avb
  avb1[smoothavi]<-paste(avb[smoothavi],"(Smooth)",sep="")
  avname_[alv>0]<-paste("factor(",avname[alv>0],")",sep="")
  fmlm<-c(fmlm,paste("+",paste(avname_,collapse="+")))
  fmlb<-c(fmlb,"Adjust")
  nadjm<-nadjm+1; tmp<-c(tmp,"I"); adjvv[[nadjm+1]]<-avname; 
  adjvb<-c(adjvb, paste(avb1, collapse="; "))
  fmlp<-c(fmlp,ifelse(!is.na(subjvname), "gee", ifelse(smoothav>0, "gam", "glm")))
  fmlm4<-c(fmlm4, list(rep("",nx)))  
}
ns=0; svb=""; smoothsv<-0
if (length(svname)>0) { 
  svb<-vlabelV[match(svname,vnameV)];
  svname_ <- svname 
  smoothsvi<-((sdf=="s" | sdf=="S") & slv==0)
  smoothsv<-sum(smoothsvi)
  smoothsvname<-svname[smoothsvi]
  svname_[smoothsvi]<-paste("rcs(",svname[smoothsvi],")",sep="")
  svb1<-svb
  svb1[smoothsvi]<-paste(svb[smoothsvi],"(Smooth)",sep="")
  svname_[slv>0]<-paste("factor(",svname[slv>0],")",sep="")
  if (par1==4) {
    fmlm<-c(fmlm, paste(fmlm[length(fmlm)],"+")); 
	fmlm4<-c(fmlm4, list(svname_));
  } else {
    fmlm<-c(fmlm,paste("+",paste(svname_,collapse="+"))); 
	fmlm4<-c(fmlm4, list(rep("",nx)));
  } 
  fmlb<-c(fmlb,"Adjust")
  nadjm<-nadjm+1; tmp<-c(tmp,"II"); adjvv[[nadjm+1]]<-svname
  adjvb2 <- paste(svb1, collapse="; ")
  if (par1 == 4) {
     if (avb>"") {adjvb2<-paste(avb1, "AND one in (", adjvb2,") respectively for each exposure X");
     } else {     adjvb2<-paste("ONE of (", adjvb2,") respectively for each exposure X");}
  }
  adjvb<-c(adjvb, adjvb2)
  fmlp<-c(fmlp,ifelse(!is.na(subjvname), "gee", ifelse(smoothsv>0, "gam", "glm")))
}
if (chk==0 & length(fmlm)>1) {
  fmlm<-fmlm[-1]; fmlb<-fmlb[-1]; tmp<-tmp[-1]; adjvv<-adjvv[-1]; adjvb<-adjvb[-1]; fmlp<-fmlp[-1]; fmlm4<-fmlm4[-1];
}
if (nadjm>1) fmlb<-paste(fmlb,tmp)
nmdl<-length(fmlm)
if (is.na(bvar) & !is.na(colvname) & nmdl==1 & par1!=3) {if ((ny==1) | (nx==1 & contx)) {bvar<-colvname; colvname<-NA;}}
if (is.na(colvname)) {
  nclv<-1; clvb<-"Total"; clvb_<-"Total"
} else {
  clv<-levels(factor(WD[,colvname])); nclv<-length(clv)+1
  clvb_<-vlabelZ[match(paste(colvname,".",clv,sep=""),vnameZ)]; clvb_[is.na(clvb_)]<-clv[is.na(clvb_)];
  clvb<-c(paste(vlabelV[vnameV==colvname],clvb_,sep="="),"Total");
  clvb_<-c(clvb_,"Total")
  WD<-WD[!is.na(WD[,colvname]),]
} 
if (is.na(bvar)) {
  blvb<-"Total"; blvb_<-"Total";nblv<-1;
} else {
  blv<-levels(factor(WD[,bvar])); nblv<-length(blv)+1
  blvb_<-vlabelZ[match(paste(bvar,".",blv,sep=""),vnameZ)]; blvb_[is.na(blvb_)]<-blv[is.na(blvb_)];
  blvb<-c(paste(vlabelV[vnameV==bvar],blvb_,sep="="),"Total");
  blvb_<-c(blvb_,"Total")
  WD<-WD[!is.na(WD[,bvar]),]
}
if (!is.na(vname.start)) {
  fmlt<-paste("Surv(,",vname.start,",",timevar,",",sep="");
  WD<-WD[!is.na(WD[,vname.start]) & !is.na(WD[,timevar]),]
} else {
  if (!is.na(timevar)) {
    fmlt<-paste("Surv(",timevar,",",sep="");
    WD<-WD[!is.na(WD[,timevar]),]
  }
}
aa<-c(1,2,3,4)
for (i in 1:4) {
      for (j in c(1:4)[-i]) {
          for (k in c(1:4)[-c(i,j)]) aa<-rbind(aa,c(i,j,k,c(1:4)[-c(i,j,k)]))
      }
}
if (nx==1) par1<-1; 
if (is.na(par3) | par3<1) par3<-1;
rord<-aa[par3,]
colprn<-prn; if (colprn>3) {colprn<-1; w<-c(w,"</br>Output format was reset to option 1");}
if (!is.na(bvar)) {prn<-"S";
} else {
  if (par3>1) {
    rordc<-ifelse(rord[4]==1,3,4)
    prn<-c("G","Y","M","X")[rord[rordc]]
    if (prn=="X") {
      if (!contx & nx>1) prn<-ifelse(nmdl>1, "M", ifelse(ny>nx & contx, "X", "Y"))
      if (!contx & nx==1) prn<-"CX"
      if (par1==2) {tmp<-ifelse(rordc==4,ifelse(rord[3]==1,2,3),2); prn<-c("G","Y","M","X")[rord[tmp]];}
    }
  } else {
    prn<-ifelse(nmdl>1, "M", ifelse(ny>nx & contx, "X", "Y"))
    if (par1==2) prn<-ifelse(nmdl>1, "M", "Y")
  }
  if (par1==3) prn<-"UM"
}
clogit<-0; 
if (!is.na(subjvname) & is.na(timevar)) clogit<-1
sink(paste(ofname,".lst",sep=""))
if (par1==2) {xbgn<-nx+1; xend<-nx+1;} else {xbgn<-1; xend<-nx;}
if (prn=="Y") {
  tt<-c(0,0,0,0,"Exposure",yb); nn<-c(0,0,0,0,yb); 
  for (k in (1:nclv)) {
    wdtmp0<-WD;
    if (!is.na(colvname)) {
      if (k<nclv) wdtmp0<-WD[WD[,colvname]==clv[k],];
      print(paste("Stratified by",colvname, ":", clvb[k]))
    } 
    for (m in 1:nmdl) {
      for (j in (xbgn:xend)) {
        colj<-cbind(k,0,m,j,xxlbl_[[j]])
        nnj <-c(k,0,m,j)
        for (i in (1:ny)) {
          tmpooi<-mdl2oo(i,j,m,k,0,wdtmp0)
          colj<-cbind(colj,tmpooi[[1]]); nnj<-c(nnj,tmpooi[[2]])
        }
        tt<-rbind(tt,colj); nn<-rbind(nn,nnj)
      }
    }
  }
} 
if (prn=="S") {
  tt<-c(0,0,0,0,"Exposure",blvb); nn<-c(0,0,0,0,blvb);
  for (k in (1:nclv)) {
    wdtmp0<-WD;
    if (!is.na(colvname)) {
      if (k<nclv) wdtmp0<-WD[WD[,colvname]==clv[k],];
      print(paste("Stratified by",colvname, ":", clvb[k]))
    } 
    for (i in (1:ny)) {
      for (m in 1:nmdl) {
        for (j in (xbgn:xend)) {
          colj<-cbind(k,i,m,j,xxlbl_[[j]]);
          nnj <- c(k,i,m,j)
          for (b in (1:nblv)) {
            print(paste("Stratified by",bvar, ":", blvb[b]))
            if (b<nblv) {wdtmp1<-wdtmp0[wdtmp0[,bvar]==blv[b],];} else {wdtmp1<-wdtmp0;}
            tmpooi<-mdl2oo(i,j,m,k,b,wdtmp1)
            colj<-cbind(colj,tmpooi[[1]]); nnj<-c(nnj,tmpooi[[2]])
          }
          tt<-rbind(tt,colj); nn<-rbind(nn,nnj)
        }
      }
    }
  }
} 
if (prn=="M") {
  tt<-c(0,0,0,0,"Exposure",fmlb);  nn<-c(0,0,0,0,fmlb)
  for (k in (1:nclv)) {
    wdtmp0<-WD; 
    if (!is.na(colvname)) {
      if (k<nclv) wdtmp0<-WD[WD[,colvname]==clv[k],];
      print(paste("Stratified by",colvname, ":", clvb[k]))
    } 
    for (i in 1:ny) {
      for (j in xbgn:xend) {
        colj<-cbind(k,i,0,j,xxlbl_[[j]]); nnj<-c(k,i,0,j)
        for (m in 1:nmdl) {
          tmpooi<-mdl2oo(i,j,m,k,0,wdtmp0)
          colj<-cbind(colj,tmpooi[[1]]); nnj<-c(nnj,tmpooi[[2]])
        }
        tt<-rbind(tt,colj); nn<-rbind(nn,nnj)
      }
    }
  }
}
if (prn=="X") {
  tt<-c(0,0,0,0,"Outcome",xb); nn<-c(0,0,0,0,xb); 
  for (k in (1:nclv)) {
    wdtmp0<-WD;
    if (!is.na(colvname)) {
      if (k<nclv) wdtmp0<-WD[WD[,colvname]==clv[k],];
      print(paste("Stratified by",colvname, ":", clvb[k]))
    } 
    for (m in 1:nmdl) {
      for (i in (1:ny)) {
        colj<-cbind(k,i,m,0,yb[i])
        nnj <-c(k,i,m,0)
        for (j in (1:nx)) {
          tmpooi<-mdl2oo(i,j,m,k,0,wdtmp0)
          colj<-cbind(colj,tmpooi[[1]]); nnj<-c(nnj,tmpooi[[2]])
        }
        tt<-rbind(tt,colj); nn<-rbind(nn,nnj)
      }
    }
  }
} 
if (prn=="CX") {
  tt<-c(0,0,0,0,"Outcome",xxlbl_[[1]][-1]); nn<-c(0,0,0,0,xb[1]); 
  nxl<-length(xxlbl_[[1]])-1
  for (k in (1:nclv)) {
    wdtmp0<-WD;
    if (!is.na(colvname)) {
      if (k<nclv) wdtmp0<-WD[WD[,colvname]==clv[k],];
      print(paste("Stratified by",colvname, ":", clvb[k]))
    } 
    for (m in 1:nmdl) {
      for (i in (1:ny)) {
          colj<-c(k,i,m,0,yb[i]); nnj <-c(k,i,m,0)
          tmpooi<-mdl2oo(i,j,m,k,0,wdtmp0)
          colj<-c(colj,tmpooi[[1]][-1]); nnj<-c(nnj,tmpooi[[2]])
          tt<-rbind(tt,colj); nn<-rbind(nn,nnj)
      }
    }
  }
} 
if (prn=="UM") {
  tt<-c(0,0,0,0,"Exposure","Univariable","Multivariable");  nn<-c(0,0,0,0,"Univariable","Multivariable")
  for (k in (1:nclv)) {
    wdtmp0<-WD; 
    if (!is.na(colvname)) {
      if (k<nclv) wdtmp0<-WD[WD[,colvname]==clv[k],];
      print(paste("Stratified by",colvname, ":", clvb[k]))
    } 
    for (i in 1:ny) {
      for (m in 1:nmdl) {
        colm<-rep(NA,6); nnm<-rep(NA,5)
        for (j in 1:(nx+1)) {
          colj<-cbind(k,i,m,j,xxlbl_[[j]]); nnj<-c(k,i,m,j)
          tmpooi<-mdl2oo(i,j,m,k,0,wdtmp0)
          colj<-cbind(colj,tmpooi[[1]]); nnj<-c(nnj,tmpooi[[2]])
          if (j<=nx) {colm<-rbind(colm,colj); nnm<-rbind(nnm,nnj);
          } else {colm<-cbind(colm[-1,],tmpooi[[1]]); nnm<-cbind(nnm[-1,],tmpooi[[2]]);}
        }
        tt<-rbind(tt,colm); nn<-rbind(nn,nnm)
      }
    }
  }
}
sink()
if (!contx & prn!="CX") rord<-c(rord[rord!=4],4)
if (prn=="X" | prn=="CX") rord<-rord[rord!=4]
if (prn=="Y") rord<-rord[rord!=2]
if (prn=="M") rord<-rord[rord!=3]
if (nx==1 & contx) rord<-rord[rord!=4]
if (nmdl==1) rord<-rord[rord!=3]
if (ny==1) rord<-rord[rord!=2]
if (is.na(colvname)) rord<-rord[rord!=1]
if (length(rord)==0) rord<-1
nrr<-length(rord)
for (i in nrr:1) {nn<-nn[order(as.numeric(nn[,rord[i]])),];tt<-tt[order(as.numeric(tt[,rord[i]])),];}
if (nn[2,4]>0) {nn[,4]<-c("Exposure",xb[as.numeric(nn[-1,4])]); } else {nn<-nn[,-4];}
if (nn[2,3]>0) {nn[,3]<-c("Model",fmlb[as.numeric(nn[-1,3])]); } else {nn<-nn[,-3];}
if (nn[2,2]>0) {nn[,2]<-c("Outcome",yb[as.numeric(nn[-1,2])]); } else {nn<-nn[,-2];}
if (!is.na(colvname)) {nn[,1]<-c(vlabelV[vnameV==colvname],clvb_[as.numeric(nn[-1,1])]);} else {nn<-nn[,-1];}
tb<-matrix(as.numeric(tt[,c(1:4)]),ncol=4); 
if (!is.na(colvname)) {tt[,1]<-c(vlabelV[vnameV==colvname],clvb[tb[-1,1]]);}
if (ny>1) {tt[,2]<-c("Outcome",yb[tb[-1,2]]);}
if (nmdl>1) {tt[,3]<-c("Model",fmlb[tb[-1,3]]);}
nrr1<-nrr-1; oo<-tt[1,]; nc<-ncol(tt)-5; nr<-nrow(tt)
for (i in 2:nr) {
  if (nrr>1) {
    for (j in 1:nrr1) {
      if (tb[i,rord[j]]!=tb[i-1,rord[j]]) oo<-rbind(oo,c(rep(tt[i,rord[j]],5),rep(" ",nc)))
    }
  }
  oo<-rbind(oo,tt[i,])
}
if (rord[nrr]!=4 & tt[1,rord[nrr]]!="0") {oo<-cbind(oo[,rord[nrr]],oo[,-(1:5)]);} else {oo<-oo[,-(1:4)]}
w<-c(w,"</br><table border=3>", mat2htmltable(oo), "</table>")
 

w<-c(w,"</br>表中数据：HR (95% CI) Pvalue *P<0.05 **P<0.01 ***P<0.001")
w<-c(w,paste(c("</br>结果变量:",paste(yb,collapse="; ")),collapse=" "))
w<-c(w,paste(c("</br>暴露变量:",paste(xb,collapse="; ")),collapse=" "))
for (m in 1:nmdl) w<-c(w,paste("</br>",fmlb[m],"model adjust for:", adjvb[m]))
if (smoothav>0 | smoothsv>0) w<-c(w,". Restricted cubic spline were applied")
if (clogit==0) {
  tmpv<-c(vname.start,timevar); tmpv<-tmpv[!is.na(tmpv)];
  w<-c(w,paste(c("</br>Cox model 时间变量:", tmpv)))
}
if (clogit==1) w<-c(w,paste(c("</br>条件 Logistic 回归, 配对组编号变量:", subjvname)))
w<-c(w,"</br></br>各模型所用的样本量</br><table border=3>", mat2htmltable(nn), "</table>")
w<-c(w,wd.subset)
w<-c(w,paste("</br>此表用易侕统计软件 (www.empowerstats.com) 和R软件生成，生成日期：",Sys.Date()))
w<-c(w,"</body></html>")
fileConn<-file(paste(ofname,".htm",sep="")); writeLines(w, fileConn)