## Created by EmpowerStats @ Wed, 17 Sep 25 17:26:28 +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_213_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('MASS','gdata','survival'); 
for (g in pkgs) {  
if (!(g %in% rownames(installed.packages()))) install.packages(g,repos='https://cloud.r-project.org'); 
}
library(MASS); 
library(gdata); 
library(survival); 
WD <- EmpowerStatsR; rm(EmpowerStatsR); gc(); 
title<-'阈值效应分析-女'; 
WD<-subset(WD, ((!is.na(SEX) & (SEX == 0)))); 
wd.subset<-paste("Use subset of data: (!is.na(SEX) & (SEX == 0))"); 
weights.var<- NA; 
yvname<-c('NAFLD'); 
ydist<-c('breslow'); 
ylink<-c('0'); 
ylv<-c(2); 
par1<-'自动寻找最佳拐点'; 
avname<-c('ALT.2','AST.2','TP.1','ALB.1','GLB.1','BUN','CR','UA','TC1','TG1','HDL1','LDL1','GLU1','SBP1','DBP1'); 
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,0); 
bvar<- NA; 
xvname<-c('CUNBAE'); 
xlv<-c(0); 
chk<- 0; 
cox<- 1; 
timevar<-'MONTH'; timevlv<- 0; 
vname.start<- NA; 
subjvname<- NA; 
gee.TYPE<-NA; 
prn<-1; 
dec<-1; 
##R package## MASS gdata 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)
}
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=" "))
}
setcox<-function(fmlx, i, wdtmp) {
  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]))
  }
  return(mdl)
}
mdl2oo<-function(mdl, xxname) {
 if (substr(mdl[1],1,5)!="Error") {
  decp<-dec+2; if (decp>4) decp<-4
  gs<-summary(mdl); print(mdl$formula); print(gs)
  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,decp);}
  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));}
 return(o1)
}
removeNA<-function(i,j,wdf) {
 vvv<-c(yvname[i],xvname[j],avname,subjvname,bvar,vname.start,timevar); 
 vvv<-vvv[!is.na(vvv)]; vvv<-vvv[vvv>" "]
 tmp<-is.na(wdf[,vvv]); 
 return(wdf[apply(tmp,1,sum)==0,])
}
get.tpval<-function(i,j,g,wdtmp,tppmin=NA,tppmax=NA) {
  if (sum(is.na(wdtmp))>0) {
    if (is.na(g)) {wdtmp<-removeNA(i,j,WD);
    } else if (g<nblv) {wdtmp0<-WD[WD[,bvar]==blv[g],]; wdtmp<-removeNA(i,j,wdtmp0);
    } else {wdtmp<-removeNA(i,j,WD); }
  }
  if (is.na(g)) {fmladj1<-fmladj;
  } else if (g<nblv) {fmladj1<-fmladj;
  } else {fmladj1<-paste(fmladj,"+factor(",bvar,")",sep="");}
  xTMP <- wdtmp[,xvname[j]]
  tmp.ss<-seq(0.05,0.95,0.05)
  tp<-quantile(xTMP,probs=tmp.ss,na.rm=TRUE) 
  tmp.llk<-rep(NA,length(tmp.ss))
  fml<-paste(xvname[j],"+tmp.X",fmladj1)
  if (!is.na(tppmin) & !is.na(tppmax)) {
    tp2.min = tppmin; tp2.max = tppmax;
  } else {
    for (k in (1:length(tmp.ss))) {
      tmp.X<-(xTMP > tp[k])*(xTMP-tp[k]); wdtmp1<-cbind(wdtmp,tmp.X)
      tmp.mdl<-setcox(fml, i, wdtmp1);
      tmp.llk[k]<-tmp.mdl$loglik[2];
      rm(wdtmp1, tmp.X)
    }
    tp1<-tmp.ss[which.max(tmp.llk)]
    tp2.min = tp1 - 0.04
    tp2.max = tp1 + 0.04
    if (tp2.min<0.05) {tp2.min=0.05}
    if (tp2.max>0.95) {tp2.max=0.95}
  }	
  tp.pctlrange<-quantile(xTMP,probs=c(tp2.min,tp2.max),na.rm=TRUE)
  tp.range<-unique(xTMP[xTMP>tp.pctlrange[1] & xTMP<tp.pctlrange[2]])
  while (length(tp.range)>5) {
    tmp.pct3<-quantile(tp.range,probs=c(0,0.25,0.5,0.75,1),type=3)
    tmp.llk3<-rep(NA,3)
    for (k in (2:4)) {
      tmp.X<-(xTMP>tmp.pct3[k])*(xTMP-tmp.pct3[k]); wdtmp1<-cbind(wdtmp,tmp.X)
      tmp.mdl<-setcox(fml, i, wdtmp1);
      tmp.llk3[k-1]<-tmp.mdl$loglik[2];
      rm(wdtmp1, tmp.X)
     }
     tmp.min3<-which.max(tmp.llk3)
     tp.range<-tp.range[tp.range>=tmp.pct3[tmp.min3] & tp.range<=tmp.pct3[tmp.min3+2]]
  }
  if (length(tp.range)>0) {
    if (length(tp.range)==1) {tp.val=tp.range[1];} else {
      tmp.llk<-rep(NA,length(tp.range))
      for (k in (1:length(tp.range))) {
         tmp.X<-(xTMP>tp.range[k])*(xTMP-tp.range[k]); wdtmp1<-cbind(wdtmp,tmp.X)
         tmp.mdl<-setcox(fml, i, wdtmp1);
         tmp.llk[k]<-tmp.mdl$loglik[2];
         rm(wdtmp1, tmp.X)
      }
      tp.val<-tp.range[which.max(tmp.llk)]
    }
  } else { tp.val<-tp.pctlrange[1];}
  return(round(tp.val,dec));
}
get2lines<-function(i,j,g,tp.value) {
  if (is.na(g)) {fmladj1<-fmladj;wdtmp<-removeNA(i,j,WD);
  } else if (g<nblv) {fmladj1<-fmladj;wdtmp0<-WD[WD[,bvar]==blv[g],];wdtmp<-removeNA(i,j,wdtmp0);
  } else {fmladj1<-paste(fmladj,"+factor(",bvar,")",sep="");wdtmp<-removeNA(i,j,WD);}
  xTMP<-wdtmp[,xvname[j]]
  tmp.X1<-(xTMP<=tp.value)*(xTMP-tp.value)
  tmp.X2<-(xTMP> tp.value)*(xTMP-tp.value)
  wdtmp1<-cbind(wdtmp,xTMP,tmp.X1,tmp.X2)
  fml0<-paste("xTMP+tmp.X2",fmladj1)
  fml1<-paste("tmp.X1+tmp.X2",fmladj1)
  fml2<-paste("xTMP",fmladj1)
  fmlp<-paste("xTMP+tmp.X2")
  tmpn<-nrow(wdtmp); 
  tmp.mdl0<-setcox(fml0,i,wdtmp1); tmp.mdl1<-setcox(fml1,i,wdtmp1); tmp.mdl2<-setcox(fml2,i,wdtmp1)
  m2<-mdl2oo(tmp.mdl2,"xTMP")
  m1<-mdl2oo(tmp.mdl1,c("tmp.X1","tmp.X2"))
  m0<-mdl2oo(tmp.mdl0,"tmp.X2")
  plrt<-pvformat(anova(tmp.mdl0,tmp.mdl2)$"P(>|Chi|)"[2],3)
  plrti<-""; plrtx<-""
  if (!is.na(g) & g==nblv) {
    fmli<-paste("xTMP+xTMP*factor(",bvar,")",fmladj)
    fmlx<-paste("tmp.X1+tmp.X2+tmp.X1*factor(",bvar,")+tmp.X2*factor(",bvar,")",fmladj)
    tmp.mdli<-setcox(fmli,i,wdtmp1);
    tmp.mdlx<-setcox(fmlx,i,wdtmp1);
    plrti<-paste("P-interaction:", pvformat(anova(tmp.mdli,tmp.mdl2)$"P(>|Chi|)"[2],3))
    plrtx<-paste("P-interaction:", pvformat(anova(tmp.mdlx,tmp.mdl1)$"P(>|Chi|)"[2],3))
  }
  oo<-list(c(plrti,m2,plrtx,tp.value,m1,m0,plrt),tmpn)
  return(oo)
}
get3lines<-function(i,j,g,tp.value) {
  if (is.na(g)) {fmladj1<-fmladj;wdtmp<-removeNA(i,j,WD);
  } else if (g<nblv) {fmladj1<-fmladj;wdtmp0<-WD[WD[,bvar]==blv[g],];wdtmp<-removeNA(i,j,wdtmp0);
  } else {fmladj1<-paste(fmladj,"+factor(",bvar,")",sep="");wdtmp<-removeNA(i,j,WD);}
  xTMP<-wdtmp[,xvname[j]]; tp1<-tp.value[1]; tp2<-tp.value[2]
  tmp.X1<- (xTMP< tp1)*(xTMP-tp1)
  tmp.X2<-((xTMP>=tp1) & (xTMP<=tp2))*(xTMP-tp1)
  tmp.X3<- (xTMP> tp2)*(xTMP-tp2)
  tmp.B1<- (xTMP>tp2)
  wdtmp1<-cbind(wdtmp,xTMP,tmp.X1,tmp.X2,tmp.X3,tmp.B1)
  fml0<-paste("xTMP+tmp.X1+tmp.X3+tmp.B1",fmladj1)
  fml1<-paste("tmp.X1+tmp.X2+tmp.X3+tmp.B1",fmladj1)
  fml2<-paste("xTMP",fmladj1)
  tmpn<-nrow(wdtmp)
  tmp.mdl0<-setcox(fml0,i,wdtmp1); tmp.mdl1<-setcox(fml1,i,wdtmp1); tmp.mdl2<-setcox(fml2,i,wdtmp1);
  m2<-mdl2oo(tmp.mdl2,"xTMP")
  m1<-mdl2oo(tmp.mdl1,c("tmp.X1","tmp.X2","tmp.X3"))
  m0<-mdl2oo(tmp.mdl0,c("tmp.X1","tmp.X3"))
  plrt<-pvformat(anova(tmp.mdl0,tmp.mdl2)$"P(>|Chi|)"[2],3)
  plrti<-""; plrtx<-""
  if (!is.na(g) & g==nblv) {
    fmli<-paste("xTMP+xTMP*factor(",bvar,")",fmladj)
    fmlx<-paste("tmp.X1+tmp.X2+tmp.X3+tmp.X1*factor(",bvar,")+tmp.X2*factor(",bvar,")+tmp.X3*factor(",bvar,")+tmp.B1",fmladj)
    tmp.mdli<-setcox(fmli,i,wdtmp1);
    tmp.mdlx<-setcox(fmlx,i,wdtmp1);
    plrti<-paste("P-interaction:", pvformat(anova(tmp.mdli,tmp.mdl2)$"P(>|Chi|)"[2],3))
    plrtx<-paste("P-interaction:", pvformat(anova(tmp.mdlx,tmp.mdl1)$"P(>|Chi|)"[2],3))
  }
  oo<-list(c(plrti,m2,plrtx,paste(tp.value,collapse=", "),m1,m0,plrt),tmpn)
  return(oo)
}
getci4tp<-function(i,j,g,tp0=NA) {
  set.seed(123456)
  if (is.na(g)) {wdt<-removeNA(i,j,WD);
  } else if (g<nblv) {wdt<-WD[WD[,bvar]==blv[g],];wdt<-removeNA(i,j,wdt);
  } else {wdt<-removeNA(i,j,WD);}
  nnwd<-nrow(wdt); tp.vv<-rep(NA,1000)
  if (!is.na(tp0)) {
    tpp0 = sum(wdt[,xvname[j]] < tp0)/length(wdt[,xvname[j]]) 
    tppmin = max(tpp0 - 0.09, 0.05) 
	tppmax = min(tpp0 + 0.09, 0.95)
  } else {
    tppmin = NA; tppmax = NA
  }
  for (s in (1:1000)) {
    WDi<-wdt[sample(1:nnwd,nnwd,replace=T),]
    tp.vv[s]<-get.tpval(i,j,NA,WDi, tppmin, tppmax); rm(WDi)
  }
  tpci<-paste(quantile(tp.vv,probs=c(0.025,0.975)),collapse=", ")
  return(tpci);
}
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]
allvname<-c(yvname,xvname,bvar,avname,subjvname,vname.start,timevar); 
allvname<-allvname[!is.na(allvname)]
WD<-WD[,allvname];
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>Spline smoothing only applies for continuous variables")
  if (!is.na(subjvname)) saf<-rep(0,length(saf))
}
if (sum(xlv>0)>0) w<-c(w,"Categorical exposure variables were ignored")
xvname<-xvname[xlv==0];
if (!is.na(subjvname)) WD<-WD[order(WD[,subjvname]),];
fmladj<-""; avb=""; smoothav<-0; 
if (length(avname)>0) {
  avb<-vlabelV[match(avname,vnameV)];
  avname_ <- avname 
  smoothavi<-((saf=="s" | saf=="S") & alv==0)
  smoothav<-sum(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="")
  fmladj<-paste("+",paste(avname_,collapse="+"))
}
if (is.na(bvar)) {
  blvb<-"N"; blvb_<-"N"; nblv<-1; blbl<-"";
} else {
  blbl<-vlabelV[match(bvar,vnameV)]; if (is.na(blbl)) blbl<-bvar;
  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(blbl,blvb_,sep="="),"Total");
  blvb_<-c(blvb_,"Total")
  WD<-WD[!is.na(WD[,bvar]),]
}
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)]
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]),]
  }
}
colprn<-prn
clogit<-0; 
if (!is.na(subjvname) & is.na(timevar)) clogit<-1
if (is.na(par1)) par1<-"";
if (is.numeric(par1)) {tp.vv<-par1;
} else {
  tmp<-as.numeric(strsplit(par1," ")[[1]]);  tp.vv<-c(tmp[!is.na(tmp)],NA)
}
prn<-ifelse(!is.na(bvar), "S", ifelse(nx>1 & ny==1, "X", "Y"));

if (length(tp.vv)>2) tp.vv<-tp.vv[1:2]
ntp<-length(tp.vv); 
getci<-FALSE
 

if (ntp==1) {
  cc0<-c("模型 I","&nbsp&nbsp一条直线效应");
  cc0<-c(cc0,"模型 II","&nbsp&nbsp折点(K)","&nbsp&nbsp &lt K 段效应 1","&nbsp&nbsp &gt K 段效应 2","&nbsp&nbsp 2与1的效应差")
  if (is.na(tp.vv[1]) & chk==1 & clogit==0) getci<-TRUE;
} else {
  cc0<-c("模型 I","&nbsp&nbsp一条直线效应");
  cc0<-c(cc0,"模型 II","&nbsp&nbsp折点(K1,K2)","&nbsp&nbsp &lt K1 段效应 1","&nbsp&nbsp K1-K2 段效应 2","&nbsp&nbsp &gt K2 段效应 3")
  cc0<-c(cc0,"&nbsp&nbsp 1与2的效应差","&nbsp&nbsp 3与2的效应差")
}
cc0<-c(cc0,"对数似然比检验");
if (getci) cc0<-c(cc0,"折点的95可信区间");
if (clogit==1 & chk==1) w<-c(w,"</br>Bootstrap for 95% CI of turning point was ignored for conditional logistic regression")
sink(paste(ofname,".lst",sep=""))
nn<-c("Outcome","Exposure",blvb);
if (prn=="Y") {
  for (j in 1:nx) {
    tt<-cc0;
    for (i in 1:ny) {
      if (is.na(tp.vv[1])) {tp.v<-get.tpval(i,j,NA,NA);} else {tp.v<-tp.vv;}
      if (ntp==1) tmpij<-get2lines(i,j,NA,tp.v);  
      if (ntp==2) tmpij<-get3lines(i,j,NA,tp.v);  
      if (getci) {tt<-cbind(tt,c(tmpij[[1]],getci4tp(i,j,NA,tp.v)));} else {tt<-cbind(tt,tmpij[[1]]);}
      nn<-rbind(nn,c(yb[i],xb[j],tmpij[[2]]))
    }
    tt<-rbind(c("Outcome: ",yb),tt)
    w<-c(w,paste("</br>For exposure:",xb[j]))
    w<-c(w,"</br><table border=3>", mat2htmltable(tt), "</table>")
  }
}
if (prn=="X") {
  for (i in 1:ny) {
    tt<-cc0;
    for (j in 1:nx) {
      if (is.na(tp.vv[1])) {tp.v<-get.tpval(i,j,NA,NA);} else {tp.v<-tp.vv;}
      if (ntp==1) tmpij<-get2lines(i,j,NA,tp.v);  
      if (ntp==2) tmpij<-get3lines(i,j,NA,tp.v);  
      if (getci) {tt<-cbind(tt,c(tmpij[[1]],getci4tp(i,j,NA,tp.v)));} else {tt<-cbind(tt,tmpij[[1]]);}
      nn<-rbind(nn,c(yb[i],xb[j],tmpij[[2]]))
    }
    tt<-rbind(c("Exposure: ",xb),tt)
    w<-c(w,paste("</br>For outcome:",yb[i]))
    w<-c(w,"</br><table border=3>", mat2htmltable(tt), "</table>")
  }
}
if (prn=="S") {
  for (i in 1:ny) {
    tt<-cc0;
    for (j in 1:nx) {
      nnij<-c(yb[i],xb[j])
      for (g in 1:nblv) {
        if (is.na(tp.vv[1])) {tp.v<-get.tpval(i,j,g,NA);} else {tp.v<-tp.vv;}
        if (ntp==1) tmpij<-get2lines(i,j,g,tp.v);  
        if (ntp==2) tmpij<-get3lines(i,j,g,tp.v);  
        if (getci) {tt<-cbind(tt,c(tmpij[[1]],getci4tp(i,j,g,tp.v)));} else {tt<-cbind(tt,tmpij[[1]]);}
        nnij<-c(nnij,tmpij[[2]])
      }
      nn<-rbind(nn,nnij)
    }
    tt<-rbind(c(blbl,blvb_),tt)
    w<-c(w,paste("</br>For outcome:",yb[i]))
    w<-c(w,paste("</br>For Exposure:",xb[j]))
    w<-c(w,"</br><table border=3>", mat2htmltable(tt), "</table>")
  }
}
sink()
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=" "))
if (length(avname)==0) avb1<-"None";
w<-c(w,paste(c("</br>调整变量:",paste(avb1,collapse="; ")),collapse=" "))
if (smoothav>0) w<-c(w,". Restricted cubic spline smoothing 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)