## Created by EmpowerStats @ Tue, 16 Sep 25 21:50:29 +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_202_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); 
print('Creating new variable: ALT.CX'); 
ALT.CX<- 0+(ALT.2>=40); 
ALT.CX[is.na(ALT.2)]<- NA; 
summary(ALT.CX); 
EmpowerStatsR<-cbind(EmpowerStatsR,ALT.CX); 
rm(BMI2,ALT.2,AST.2,TP.1,ALB.1,GLB.1,CUNBAE,CUNBAE.Z,ALT.CX); 
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','ALT.CX','ALT.CX.0','ALT.CX.1')[-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','ALT categorical','  <40','  >=40')[-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','ALT.CX'); 
pkgs<-c('mgcv','geepack','survival'); 
for (g in pkgs) {  
if (!(g %in% rownames(installed.packages()))) install.packages(g,repos='https://cloud.r-project.org'); 
}
library(mgcv); 
library(geepack); 
library(survival); 
WD <- EmpowerStatsR; rm(EmpowerStatsR); gc(); 
title<-'交互作用检验ALT'; 
wd.subset=''; 
weights.var<- NA; 
yvname<-c('NAFLD'); 
ydist<-c('breslow'); 
ylink<-c('0'); 
ylv<-c(2); 
colvname<-'ALT.CX'; colvlv<- 2; 
avname<-c('AST.2','TP.1','ALB.1','GLB.1','BUN','CR','UA','TC1','TG1','HDL1','LDL1','GLU1','SBP1','DBP1','SEX'); 
saf<-c(0,0,0,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,0,0,2); 
svname<-c(); 
xvname<-c('CUNBAE.Z'); 
xlv<-c(0); 
bvar<- NA; 
chk<- 1; 
cox<- 1; 
timevar<-'MONTH'; timevlv<- 0; 
vname.start<- NA; 
subjvname<- NA; 
gee.TYPE<-NA; 
prn<-1; 
dec<-1; 
##R package## mgcv geepack survival ##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=" "))
}
chkMatrixSinglar <- function(mat) {
  oo<-rep(0,ncol(mat))
  tmp<-mat[apply(is.na(mat),1,sum)==0,]
  oo[apply(tmp,2,function(z) length(levels(as.factor(z)))==1)]<-1
  rr<-cor(tmp)
  for (i in (2:ncol(rr))) if (sum(rr[1:(i-1),i]>0.999)>0) oo[i]<-1
  if (sum(oo)>0) print(rr)
  return(oo)
}
mdl2oo<-function(mdl, mdlo, xxname, mdl2) {
  gs<-summary(mdl); print(mdl$formula); print(gs)
  gs2<-summary(mdlo);  print(gs2)
  gsparm <- gs$coefficients;tmpn <- gs$n; 
  pintr<- 1-pchisq(2*(gs$loglik[2]-gs2$loglik[2]),gs$logtest[2]-gs2$logtest[2]) 
  pintr<-pvformat(pintr,4)
  if (!is.na(mdl2[[1]][1])) {
    gs3<-summary(mdl2); print(gs3)
    pintr2<- 1-pchisq(2*(gs3$loglik[2]-gs2$loglik[2]),gs3$logtest[2]-gs2$logtest[2]) 
    pintr<-paste(pintr, " (", pvformat(pintr2,4), " #)",sep="")
  }
  gsp<-gsparm[match(xxname,rownames(gsparm)),]
  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,4);}
  o1<-paste(numfmt(exp(beta),dec)," (",numfmt(exp(ci1),dec),", ",numfmt(exp(ci2),dec),")",sep="")
  o1<-paste(o1,pvv)
  return(list(o1,pintr))
}
mdl2bb<-function(mdl,xxname) {
  gs<-summary(mdl); print(mdl$formula); print(gs)
  gsparm <- gs$coefficients;
  bb<-round(gsparm[match(xxname,rownames(gsparm)),1],4)
  return(bb)
}
vv2fml<-function(vvname,vvlv) {
  vvname_ <- vvname
  vvname_[vvlv>0]<-paste("factor(",vvname[vvlv>0],")",sep="")
  fml<-paste(vvname_,collapse="+")
  return(fml)
}
addvia<-function(vvname,vvlv,iavname) {    
  vvname_ <- vvname
  vvname_[vvlv>0]<-paste("factor(",vvname[vvlv>0],")",sep="")
  fml0<-paste(c(vvname_,paste("factor(",iavname,")",sep="")),collapse="+")
  vvname_2<-vvname
  vvname_2<-paste(vvname_2,"*factor(",iavname,")",sep="")
  vvname_2[vvlv>0]<-paste("factor(",vvname[vvlv>0],")*factor(",iavname,")",sep="")
  for (z in (1:length(vvname_))) {fml0<-c(fml0,paste(c(vvname_[-z],vvname_2[z]),collapse="+"));}
  vvia<-rep(FALSE,length(vvname))
  for (i in (1:ny)) {
    for (j in (1:nx)) {
        xxname<- xvianame_[[j]]; if (xlv[j]>0) xxname<-xxname[-1]
        xxname1<-paste(xxname,collapse="+")
        fml<-paste(fmlcox0,yvname[i],")~",xxname1,"+",fml0,fmlst); bb<-NA
        for (k in (1:length(fml))) {
          cox0<-try(coxph(formula(fml[k]),weights=wdtmp$weights,data=wdtmp,method=ydist[i]))
          bb<-rbind(bb,mdl2bb(cox0, xxname));
        }
        bb<-bb[-1,]; tmprname_=paste("Model",1:length(fml)); rownames(bb)<-tmprname_; 
        tmpmdls_<-cbind(tmprname_,fml);  print(bb); print(tmpmdls_);
        bba<-(apply((bb< matrix(rep(bb[1,]-abs(bb[1,])*0.1,nrow(bb)),nrow=nrow(bb),byrow=TRUE)),1,sum)>0)
        bbu<-(apply((bb> matrix(rep(bb[1,]+abs(bb[1,])*0.1,nrow(bb)),nrow=nrow(bb),byrow=TRUE)),1,sum)>0)
        vvia <- (vvia | (bba | bbu)[-1])
    }  
  }
  fmlx<-paste(c(paste("factor(",iavname,")",sep=""),vvname_[!vvia],vvname_2[vvia]),collapse="+")
  return(list(fmlx,vvia,vvname[vvia]))
}
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"];
}}
colprn<-prn; if (colprn==0) colprn<-1;
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]
sink(paste(ofname,".lst",sep=""))
Total_<-1; WD<-cbind(WD,Total_)
bvname<-bvar;
if (is.na(bvar)) {bvar<-"Total_"; bv<-WD[,bvar];bvb<-bvar;}
blv<-levels(factor(WD[,bvar])); nlb<-length(blv)
if (!is.na(bvname)) {
  bvb<-vlabelV[match(bvar,vnameV)][1]; if (is.na(bvb)) bvb<-bvar
  blb<-vlabelZ[match(paste(bvname,".",blv,sep=""),vnameZ)]
  blb[is.na(blb)]<-blv[is.na(blb)]
  stvlb<-c(paste(bvb,": ",blb,sep=""),"Total")
}
nx<-length(xvname); xb<-vlabelV[match(xvname,vnameV)]; xb[is.na(xb)]<-xvname[is.na(xb)]
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)]
clv<-levels(factor(WD[,colvname])); nlv<-length(clv)
colvb<-vlabelV[match(colvname,vnameV)]; if (is.na(colvb)) colvb<-colvname;
clb<-vlabelZ[match(paste(colvname,".",clv,sep=""),vnameZ)]
clb[is.na(clb)]<-clv[is.na(clb)]
xvname_ <- xvname
xvname_[xlv>0]<-paste("factor(",xvname[xlv>0],")",sep="")
xiaxname_ <- paste(xvname_,"+",xvname_,"*",colvname,sep="")
xvianame_<-list(NA); xvialbl_<-list(NA); 
for (j in (1:nx)) {
  tmpnewxv_<-NA; tmpnewxm_<-NA; tmpnewxb1<-""; tmpnewxb2<-"";
  if (xlv[j]==0) {
    for (k in (1:nlv)) {
      tmpnewxv_<-cbind(tmpnewxv_,WD[,xvname[j]]*(WD[,colvname]==clv[k]))      
      tmpnewxm_<-c(tmpnewxm_,paste(xvname[j],".",k,sep=""))
      tmpnewxb1<-c(tmpnewxb1,xb[j]);tmpnewxb2<-c(tmpnewxb2,clb[k])
    }
  } else {
    xjlv <- levels(factor(WD[,xvname[j]])); nxjlv<-length(xjlv)
    xjlb <- vlabelZ[match(paste(xvname[j],".",xjlv,sep=""),vnameZ)] 
    xjlb[is.na(xjlb)]<-xjlv[is.na(xjlb)]
    tmpxcc_<-rep(NA,nrow(WD)); tmpcode_<-0
    for (k in (1:nlv)) {
      for (z in (1:nxjlv)) {
        tmprzk_<-(WD[,xvname[j]]==xjlv[z] & WD[,colvname]==clv[k])
        tmpxcc_[tmprzk_]<-tmpcode_; tmpcode_<-tmpcode_+1
        tmpnewxv_<-cbind(tmpnewxv_,tmprzk_*1)
        tmpnewxm_<-c(tmpnewxm_,paste(xvname[j],".",z,".",colvname,".",k,sep=""))
        tmpnewxb1<-c(tmpnewxb1,xjlb[z]);tmpnewxb2<-c(tmpnewxb2,clb[k])
      }
    }
    tmpnewxv_<-cbind(tmpnewxv_,tmpxcc_);
    tmpnewxm_<-c(tmpnewxm_,paste(xvname[j],".",colvname,sep=""))
  }
  tmpnewxv_<-tmpnewxv_[,-1]; tmpnewxm_<-tmpnewxm_[-1]; 
  tmpnewxb1<-tmpnewxb1[-1]; tmpnewxb2<-tmpnewxb2[-1]
  colnames(tmpnewxv_)<-tmpnewxm_
  WD<-cbind(WD,tmpnewxv_)
  if (xlv[j]==0) {xvianame_[[j+1]]<-tmpnewxm_;} else {xvianame_[[j+1]]<-tmpnewxm_[-length(tmpnewxm_)];}
  xvialbl_[[j+1]]<-list(tmpnewxb1,tmpnewxb2)
}
xvianame_<-xvianame_[-1]; xvialbl_<-xvialbl_[-1]
cond.logit<-0; fmlst<-"";
if (!is.na(vname.start)) {
  fmlcox0 <- paste("Surv(",vname.start,",",timevar,",",sep="")
  cmp0<-(!is.na(vname.start) & !is.na(timevar))
} else {
  if (is.na(timevar)) {
    timevar<-"tmp.time"; tmp.time<-rep(1,nrow(yv)); WD<-cbind(WD,tmp.time);
    ydist<-rep('exact',ny); cond.logit<-1;
    fmlcox0 <- "Surv(tmp.time,";
    fmlst<-paste("+strata(",subjvname,")",sep="")
  } else {
    fmlcox0 <- paste("Surv(",timevar,",",sep="");
    cmp0<-(!is.na(timevar))
  }
}  
allvname<-c(yvname,xvname,colvname,bvar,avname,svname,subjvname,vname.start,timevar,"weights"); 
allvname<-allvname[!is.na(allvname) & allvname>""]
cmpd<-WD[,allvname];
wdtmp<-WD[(apply(is.na(cmpd),1,sum)==0),]
if (!is.na(subjvname)) wdtmp<-wdtmp[order(wdtmp[,subjvname]),]
fmlxx<-""; fmlzz<-""; fmltt<-"Crude"; fmlvv<-list(""); 
if (length(avname)>0) {
  avb<-vlabelV[match(avname,vnameV)]; avb[is.na(avb)]<-avname[is.na(avb)]
  avbstr<-paste(avb,collapse=", ")
  fmlxx<-c(fmlxx,vv2fml(avname,alv))
  fmlzz<-c(fmlzz,avbstr)
  fmltt<-c(fmltt,"Model I")
  fmlvv[[length(fmlvv)+1]]<-avname
  if (chk==1) {
    tmp<-addvia(avname,alv,colvname)
    if (sum(tmp[[2]])>0) {
      fmlxx<-c(fmlxx,tmp[[1]]); 
      str1<- paste(vlabelV[match(tmp[[3]],vnameV)],collapse=", ")
      fmlzz<-c(fmlzz,paste(avbstr, " and the interaction terms for following variables:", str1));
      fmltt<-c(fmltt,"Model I*")
      fmlvv[[length(fmlvv)+1]]<-avname
    }
  }
}
if (length(svname)>0) {
  svbstr<-paste(sb,collapse=", ")
  fmlxx<-c(fmlxx,vv2fml(svname,slv))
  fmlzz<-c(fmlzz,svbstr)
  fmltt<-c(fmltt,"Model II")
  fmlvv[[length(fmlvv)+1]]<-svname
  if (chk==1) {
    tmp<-addvia(svname,slv,colvname)
    if (sum(tmp[[2]])>0) {
      fmlxx<-c(fmlxx,tmp[[1]]); 
      str1<- paste(vlabelV[match(tmp[[3]],vnameV)],collapse=", ")
      fmlzz<-c(fmlzz,paste(svbstr, " and the interaction terms for following variables:", str1));
      fmltt<-c(fmltt,"Model II*")
      fmlvv[[length(fmlvv)+1]]<-svname
    }
  }
}
nmdl<-length(fmlxx)
if (nlb>1) {stmp_<-c(1:nlb,0);} else {stmp_<-0;}
oo1<-list(NA,NA,NA,NA); oo2<-list(NA,NA,NA,NA,NA,NA,NA);
for (s in stmp_) {
  for (i in (1:ny)) {
    for (j in (1:nx)) {
      ifelse(xlv[j]>0, tmpvxz<-paste(xvname[j],".",colvname,sep=""), tmpvxz<-colvname) 
      wdtmp0<-WD[(apply(is.na(WD[,c(yvname[i],xvname[j],colvname)]),1,sum)==0),]
      numi<-tapply(!is.na(wdtmp0[,yvname[i]]),as.factor(wdtmp0[,tmpvxz]),sum)
      oo1<-rbind(oo1,list(s,i,j,numi))
      for (k in (1:nmdl)) {
        tmpvn<-c(yvname[i],xvname[j],colvname);
        if (fmlvv[[k]][1]>"") tmpvn<-c(tmpvn,fmlvv[[k]])
        wdtmp<-WD[(apply(is.na(WD[,tmpvn]),1,sum)==0),]
        if (s>0) wdtmp<-wdtmp[wdtmp[,bvar]==blv[s],]
        nn<-tapply(!is.na(wdtmp[,yvname[i]]),as.factor(wdtmp[,tmpvxz]),sum)
        xxname<- xvianame_[[j]];
        if (xlv[j]>0) {
          xxname<-xxname[-1]
          fmli<-paste(fmlcox0,yvname[i],")~",paste(xxname,collapse="+"), fmlst,sep="")
        } else{  
          fmli<-paste(fmlcox0,yvname[i],")~",paste(xxname,collapse="+"),"+factor(",colvname,")", fmlst,sep="")
        }
        fmlo<-paste(fmlcox0,yvname[i],")~",xvname_[j],"+factor(",colvname,")",fmlst,sep="")
        fml2<-paste(fmlcox0,yvname[i],")~",xiaxname_[j],"+factor(",colvname,")",fmlst,sep="")
        if (fmlxx[k]>"  ") {
          fmli<-paste(fmli,"+",fmlxx[k],sep=""); 
          fmlo<-paste(fmlo,"+",fmlxx[k],sep="")
          fml2<-paste(fml2,"+",fmlxx[k],sep="")
        }
        if (s==0 & !is.na(bvname)) {
          fmli<-paste(fmli,"+factor(",bvname,")",sep="")
          fmlo<-paste(fmlo,"+factor(",bvname,")",sep="")
          fml2<-paste(fml2,"+factor(",bvname,")",sep="")
        }
        if (min(apply(wdtmp[,c(xxname,colvname)],2,function(m) length(unique(m))))>1) {
          coxi<-try(coxph(formula(fmli),weights=wdtmp$weights,data=wdtmp,method=ydist[i]))
          coxo<-try(coxph(formula(fmlo),weights=wdtmp$weights,data=wdtmp,method=ydist[i]))
          cox2<-NA; if (nlv>2) cox2<-try(coxph(formula(fml2),weights=wdtmp$weights,data=wdtmp,method=ydist[i]))
          mdli<-mdl2oo(coxi,coxo,xxname,cox2)
        } else {
          mdli<-list(rep("",times=length(xxname)),"");
        }
        oo2<-rbind(oo2,list(s,i,j,k,mdli[[1]],mdli[[2]],nn))
      }
    }
  }
}

ncontx<-sum(xlv==0); ncatgx<-nx-ncontx
w<-""
if (ncontx>0) {
  tt1<-c("Model",paste(colvname,": ",clb,sep=""),"P interaction")
  nn1<-c("Y","Strata","X","Model",clb)
  for (i in (1:ny)) {
    if (ny>1) tt1<-rbind(tt1,c(paste("For ",yb[i],sep=""),rep(" ",nlv+1)))
    for (s in stmp_) {
       if (nlb>1) {
         ifelse(s>0, sslb<-stvlb[s], sslb<-stvlb[nlb+1])
         tt1<-rbind(tt1,c(sslb,rep(" ",nlv+1)))
       } else {sslb<-"Total";}
       for (j in (1:nx)) {
         if (xlv[j]==0) {
           if (nx>1) tt1<-rbind(tt1,c(xb[j],rep(" ",nlv+1)))
           for (k in (1:nmdl)) {
             oo2row <-(oo2[,1]==s & oo2[,2]==i & oo2[,3]==j & oo2[,4]==k)
             oo2row[is.na(oo2row)]<-FALSE
             mdloo<-oo2[oo2row,][[5]]; pia<-oo2[oo2row,][[6]];mdlnn<-oo2[oo2row,][[7]]
             tt1<-rbind(tt1,c(fmltt[k],mdloo,pia))
             nn1<-rbind(nn1,c(yb[i],sslb,xb[j],fmltt[k],mdlnn))
           }
         }
       }
    }
  }
}

if (ncatgx>0) {
  tt2<-c("Exposure",colvb,"N")
  if (nmdl==1) mdlk<-1
  if (nmdl==2) mdlk<-2
  if (nmdl==3) {if (fmlxx[3]=="Model II") {mdlk<-c(2,3);} else {mdlk<-3;}}
  if (nmdl==4) {
    if (fmlxx[4]=="Model II") mdlk<-c(3,4)
    if (fmlxx[4]=="Model II*") mdlk<-c(2,4)
  }
  if (nmdl==5) mdlk<-c(3,5)
  mdlz<-(1:nmdl)[-mdlk]; nmdlz<-length(mdlz)
  tt2<-c(tt2,fmltt[mdlk])
  if (nmdlz>0) {tt3<-c("Exposure",colvb,fmltt[mdlz]);} else tt3<-NA
  tbncol<-length(tt2)
  nn2<-c("Y","Strata","X","Model","Nobs","Tot N")
  for (i in (1:ny)) {
    if (ny>1) {
      tt2<-rbind(tt2,c(paste("For ",yb[i],sep=""),rep(" ",tbncol-1)))
      if (nmdlz>0) tt3<-rbind(tt3,c(paste("For ",yb[i],sep=""),rep(" ",nmdlz+1)))
    }
    for (s in stmp_) {
       if (nlb>1) {
         ifelse(s>0, sslb<-stvlb[s], sslb<-stvlb[nlb+1])
         tt2<-rbind(tt2,c(sslb,rep(" ",tbncol-1)))
         if (nmdlz>0) tt3<-rbind(tt3,c(sslb,rep(" ",nmdlz+1)))
       } else {sslb<-"Total";}
       for (j in (1:nx)) {
         if (xlv[j]>0) {
           oo1row <-(oo1[,1]==s & oo1[,2]==i & oo1[,3]==j)
           oo1row[is.na(oo1row)]<-FALSE
           mdlcc<-cbind(c(xb[j],xvialbl_[j][[1]][[1]]),c(" ",xvialbl_[j][[1]][[2]]))
           if (nmdlz>0) mdlcx<-mdlcc
           mdlcc<-cbind(mdlcc,c(" ",oo1[oo1row,][[4]]))
           pia<-c("P interaction"," "," ")
           if (nmdlz>0) pix<-c("P interaction"," ")
           for (k in mdlk) {
             oo2row <-(oo2[,1]==s & oo2[,2]==i & oo2[,3]==j & oo2[,4]==k)
             oo2row[is.na(oo2row)]<-FALSE
             mdloo<-c("","Ref.",oo2[oo2row,][[5]]); 
             pia<-c(pia,oo2[oo2row,][[6]]);
             mdlnn<-oo2[oo2row,][[7]]; sumnn<-sum(mdlnn)
             mdlcc<-cbind(mdlcc,mdloo)             
             nn2<-rbind(nn2,c(yb[i],sslb,xb[j],fmltt[k],paste(mdlnn,collapse=";"),sumnn))
           }
           tt2<-rbind(tt2,mdlcc,pia)
           if (nmdlz>0) {
             for (k in mdlz) {
               oo2row <-(oo2[,1]==s & oo2[,2]==i & oo2[,3]==j & oo2[,4]==k)
               oo2row[is.na(oo2row)]<-FALSE
               mdloo<-c("","Ref.",oo2[oo2row,][[5]]); 
               pix<-c(pix,oo2[oo2row,][[6]]);
               mdlnn<-oo2[oo2row,][[7]]; sumnn<-sum(mdlnn)
               mdlcx<-cbind(mdlcx,mdloo)             
               nn2<-rbind(nn2,c(yb[i],sslb,xb[j],fmltt[k],paste(mdlnn,collapse=";"),sumnn))
             }
             tt3<-rbind(tt3,mdlcx,pix)
           }
         }
       }
    }
  }
}
sink()
w<-c("<!DOCTYPE html><html lang='zh'><head><meta charset='utf-8'></head><body>")
w<-c(w,paste("<h2>", title, "</h2>"))
if (colprn<=2) {prnopt<-"HR (95%CI) Pvalue ";} else {prnopt<-"HR (95%CI) ";}
if (cond.logit==1) {w<-c(w, "Conditional logistic regression</br>");
} else {w<-c(w, "Proportional Hazards Model (Cox model)</br>");}

 

if (ncontx>0) {
  w<-c(w,"</br><table border=3>", mat2htmltable(tt1), "</table>")
  w<-c(w,"</br>Results in table:", prnopt)
  if (nlv>2) w<-c(w,"</br># interaction trend test (treat ", colvb, " as continuous)")
  w<-c(w,paste(c("</br>结果变量:",yb),collapse=" "))
  w<-c(w,paste(c("</br>危险因素:",xb[xlv==0]),collapse=" "))
  w<-c(w,paste(c("</br>效应修饰因子:",colvb),collapse=" "))
  for (k in (1:nmdl)) {
    if (fmlxx[k]>"") w<-c(w,paste("</br>",fmltt[k]," 调整了: ",fmlzz[k],sep=""))
  }
  w<-c(w,"</br></br>各模型所用的样本量</br><table border=3>", mat2htmltable(nn1), "</table>")
}
if (ncatgx>0) {
  w<-c(w,"</br><table border=3>", mat2htmltable(tt2), "</table>")
  w<-c(w,"</br>Results in table:", prnopt)
  if (nlv>2) w<-c(w,"</br># interaction trend test (treat ", colvb, " as continuous)")
  w<-c(w,paste(c("</br>结果变量:",yb),collapse=" "))
  w<-c(w,paste(c("</br>危险因素:",xb[xlv>0]),collapse=" "))
  w<-c(w,paste(c("</br>效应修饰因子:",colvb),collapse=" "))
  for (k in mdlk) {
    if (fmlxx[k]>"") w<-c(w,paste("</br>",fmltt[k]," 调整了: ",fmlzz[k],sep=""))
  }
  if (!is.na(tt3[1])) {
    w<-c(w,"</br><table border=3>", mat2htmltable(tt3), "</table>")
    w<-c(w,"</br>Results in table:", prnopt)
    w<-c(w,paste(c("</br>结果变量:",yb),collapse=" "))
    w<-c(w,paste(c("</br>危险因素:",xb[xlv>0]),collapse=" "))
    w<-c(w,paste(c("</br>效应修饰因子:",colvb),collapse=" "))
    for (k in mdlz) {
      if (fmlxx[k]>"") w<-c(w,paste("</br>",fmltt[k]," 调整了: ",fmlzz[k],sep=""))
    }
  }
  w<-c(w,"</br></br>各模型所用的样本量</br><table border=3>", mat2htmltable(nn2), "</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)