## Created by EmpowerStats @ Tue, 16 Sep 25 19:46:50 +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_184_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: 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,CUNBAE,CUNBAE.Z); 
sink(); 
vname<-c(NA,'SEX','SEX.0','SEX.1','AGE','NAFLD','NAFLD.0','NAFLD.1','BMI','MONTH','BMI2','CUNBAE','CUNBAE.Z')[-1]; 
vlabel<-c(NA,'SEX','  0','  1','AGE','NAFLD','  0','  1','BMI','MONTH','BMI2','CUNBAE','CUNBAE Z-score')[-1]; 
varused4this <- c('SEX','AGE','NAFLD','BMI','MONTH','BMI2','CUNBAE','CUNBAE.Z'); 
pkgs<-c('survival','gdata','Hmisc','Matrix','survey'); 
for (g in pkgs) {  
if (!(g %in% rownames(installed.packages()))) install.packages(g,repos='https://cloud.r-project.org'); 
}
library(survival); 
library(gdata); 
library(Hmisc); 
library(Matrix); 
library(survey); 
WD <- EmpowerStatsR; rm(EmpowerStatsR); gc(); 
title<-'Cox 模型'; 
wd.subset=''; 
svy<- 0; 
weights.var<- NA; 
colvname<-'NAFLD'; colvlv<- 2; 
xvname<-c('CUNBAE.Z'); 
sxf<-c(0); 
xlv<-c(0); 
par1<-1; 
timevar<-'MONTH'; timevlv<- 0; 
vname.start<- NA; 
vname.stop<- NA; 
par3<- NA; 
bvar<- NA; 
dec<-4; 
##R package## survival gdata Hmisc Matrix survey ##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=" "))
}
runcoxfml<-function(fml,fm12,ofname,xxstr) {
 nx<-length(xvname)
 if (svy==1) {
   print(summary(svy.DSN));
   tmp.cox<-svycoxph(formula(fml),svy.DSN, method=mth,weights=WD1$weights,data=WD1); 
 } else {
   if (!is.na(weights.var)) {tmp.cox<-coxph(formula(fml),method=mth,weights=WD1$weights,robust=TRUE,data=WD1); 
   } else {tmp.cox<-coxph(formula(fml),method=mth,data=WD1); }
 }
 print(summary(tmp.cox))
 if (is.na(weights.var) & svy==0) {
  tmp.zph <- cox.zph(tmp.cox); print(tmp.zph)
  martingale<-resid(tmp.cox,type="martingale")
  deviance<-resid(tmp.cox,type="deviance")
  if (mth != "exact") {
    res.score<-resid(tmp.cox,type="score")
    if (!is.matrix(res.score)) {res.score<-matrix(res.score,ncol=1); colnames(res.score)<-"score";
    } else {colnames(res.score)<-paste("score",colnames(res.score),sep=".")
    }
    resid<-cbind(martingale,deviance,res.score)
    if (nrow(resid)==nrow(WD1)) resid<-cbind(WD1,resid)
    write.table(resid,file=paste(ofname,"_resid.xls",sep=""),col.names=TRUE,row.names=FALSE,sep="\t")
    schoenfeld<-resid(tmp.cox,type="schoenfeld")
    event.time<-as.numeric(rownames(schoenfeld))
    schoenfeld<-cbind(event.time,schoenfeld)
    write.table(schoenfeld,file=paste(ofname,"_schoenfeld.xls",sep=""),col.names=TRUE,row.names=FALSE,sep="\t")
    scho.time <-cbind(round(tmp.zph$x,4), as.numeric(rownames(tmp.zph$y)), round(tmp.zph$y,4))
    colnames(scho.time)<-c("time-transform", timevar, colnames(tmp.zph$y))
    write.table(scho.time,file=paste(ofname,"resd_time.xls",sep="_"),col.names=TRUE,row.names=FALSE,sep="\t") 
    for (i in (1:nx)) {
      png(paste(ofname,xvname[i],"martingale.png",sep="_"),width=720,height=560)
      plot(WD1[,xvname[i]],martingale, xlab=xb.utf8[i], ylab="Martingale residual"); lines(lowess(WD1[,xvname[i]],martingale),col="red")
      dev.off()
    }
  }
  if (!strata) {
    cname<-colnames(tmp.zph$y); np<-length(cname); nr<-round(sqrt(np),0); nc<-ceiling(np/nr)
	png(paste(ofname,"beta_time.png",sep="_"),width=960,height=720)
    par(mfrow=c(nr,nc))
    for (i in 1:np) try(plot(tmp.zph[i]))
	dev.off()  
  }
 }
 if (fml2>"") {
   sfit<-survfit(formula(fml2),weights=WD1$weights,data=WD1)
   tmp.miny<-floor(min(sfit$lower,na.rm=TRUE)*10)/10
   png(paste(ofname,"_survfit.png",sep=""),width=720,height=560)
   if (!is.na(vname.stop)) {
     plot(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)
   } else {
     plot(sfit,ylim=c(tmp.miny,1),ylab="Survival",xlab=timevb.utf8,main=colvb.utf8)
   }
   dev.off()
 }
 ss<-summary(tmp.cox); print(fml); print(ss);
 coe<-ss$coefficients; if (!is.matrix(coe)) coe<-matrix(coe,nrow=1)
 xxname<-rownames(coe)
 if (is.na(weights.var) & svy==0) { 
   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]))
   if (!is.matrix(coe0)) coe0<-matrix(coe0,nrow=1)
   coe<-cbind(numfmt(coe0,4),pvformat(coe[,5],4))
   if (!is.matrix(coe)) coe<-matrix(coe,nrow=1)
   zph<-numfmt(tmp.zph$table,4)
   oo1<-cbind(xxstr[[1]],c(ss$n,ss$nevent))
   oo2<-rbind(xxstr[[2]],cbind(xxname,coe))
   oo3<-paste(paste(names(ss$rsq),"=",numfmt(ss$rsq,4)),collapse=", ")
   oo4<-paste(xxstr[[4]],ss$loglik[2], ", AIC=", round(extractAIC(tmp.cox)[2],1))
   oo5<-numfmt(rbind(ss$logtest,ss$waldtest,ss$sctest),5)
   oo5<-cbind(c(" ","Likelihood ratio test","Wald test","Score(log rank)test"),rbind(colnames(oo5),oo5))
   oo6<-cbind(c("",rownames(zph)),rbind(colnames(zph),zph))
   if (length(unlist(strsplit(survStr,",")))<3) {
     hc<-rcorrcens(formula(paste(survStr,"~I(-1*predict(tmp.cox))")))
     oo7<-rbind(c(" ",colnames(hc)),cbind("Model",round(hc,4)))
   } else {oo7<-NULL;}	
 } else {
    #R363 is different from R343 with output of this
    coe0<-cbind(coe[,1],coe[,4],coe[,5],coe[,2],1/coe[,2], exp(coe[,1]-1.96*coe[,4]), exp(coe[,1]+1.96*coe[,4]))
    if (!is.matrix(coe0)) coe0<-matrix(coe0,nrow=1)
    coe<-cbind(numfmt(coe0,4),pvformat(coe[,6],4))
    if (!is.matrix(coe)) coe<-matrix(coe,nrow=1)
    oo1<-cbind(xxstr[[1]],c(ss$n,ss$nevent))
    oo2<-rbind(xxstr[[2]],cbind(xxname,coe))
 }
 ww<-"";
 if (svy==1) ww<-c(ww, "</br></br>Survey design (check .log or .lst for details)");
 ww<-c(ww,"</br></br>model: ", fml);
 ww<-c(ww,"</br></br><table border=3>",mat2htmltable(oo1),"</table>")
 ww<-c(ww,"</br>",xxstr[[3]],"<table border=3>",mat2htmltable(oo2),"</table>")
 if (is.na(weights.var) & svy==0) {
   ww<-c(ww,paste("</br>",oo3))
   ww<-c(ww,paste("</br>",oo4))
   ww<-c(ww,"</br><table border=3>",mat2htmltable(oo5),"</table>")
   ww<-c(ww,"</br>",xxstr[[5]],"<table border=3>",mat2htmltable(oo6),"</table></br>")
   if (!is.null(oo7)) {
     ww<-c(ww,"</br>Somers' Rank Correlation for Censored Data</br><table border=3>",mat2htmltable(oo7),"</table></br>")
     ww<-c(ww,"C: Harrell\'s C index; Dxy = 2(C - 0.5), Dxy is a U-statistic.")
     ww<-c(ww,"</br>SD: Standard error of Dxy, so, SE. for C = SD/2 </br>")
   }
 }  
 if (!strata & (nx>1) & svy==0 & is.na(weights.var)) {
   tmp.stp<-step(tmp.cox)
   ww<-c(ww,"</br>Step wise model selection<table border=3>",mat2htmltable(as.matrix(tmp.stp$anova)), "</table></br>")
   cos<-summary(tmp.stp)$coefficients
   if (length(cos)>0) {
     if (!is.matrix(cos)) cos<-matrix(cos,nrow=1)
     xxname<-rownames(cos)
     cos<-cbind(cos[,1],cos[,3],cos[,4],cos[,2],1/cos[,2], exp(cos[,1]-1.96*cos[,3]), exp(cos[,1]+1.96*cos[,3]), cos[,5])
     cos<-numfmt(cos,4)
     oos<-rbind(xxstr[[2]],cbind(xxname,cos))
     ww<-c(ww,"</br>Final model:<table border=3>",mat2htmltable(oos),"</table></br>")
   }
 }
 return(ww)
}
vlabelN<-(substr(vlabel,1,1)==" ");
vlabelZ<-vlabel[vlabelN];vlabelV<-vlabel[!vlabelN]
vnameV<-vname[!vlabelN];vnameZ<-vname[vlabelN]
sxf <- toupper(sxf)
if (!is.na(weights.var)) {weights<-WD[,weights.var];} else {weights<-1;}
WD<-cbind(WD,weights);
allvname<-c(xvname,colvname,bvar,vname.start,vname.stop,timevar,"weights"); 
allvname<-allvname[!is.na(allvname)]
WD<-WD[,allvname]; WD<-WD[(apply(is.na(WD),1,sum)==0),]
if (!is.na(par3)) {
 if (is.numeric(par3)) {tp<-par3;} else {tp<-as.numeric(unlist(strsplit(par3," ")));}
 if (sum(is.na(tp))==0) bvar="TIME.PHASE.4LMCOX";
}
xv1<-xvname; xv1[xlv>2]<-paste("factor(",xvname[xlv>2],")",sep="")
survStr<-""; fml<-""; fml2<-""
strata<-FALSE
if (!is.na(bvar) & bvar=="TIME.PHASE.4LMCOX") {
  if (sum(sxf=="S")==0) {
    WD<-getLandmarkDF(colvname, 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]));
    }
  } else if (is.na(vname.start)) {
    WD<-survSplit(Surv(WD[,timevar], WD[,colvname])~., data=WD, cut=tp, episode="timeSeg", id="tmp.ID")
	xv2<-xv1
	xv2[sxf=="S"]<-paste(xv2[sxf=="S"],"strata(timeSeg)",sep=":")
	fml<-paste("Surv(tstart, tstop, event) ~", paste(xv2,collapse="+"))
	if (!is.na(vname.stop)) fml<-paste(fml,"+strata(",vname.stop,")",sep=""); 
	survStr<-"Surv(tstart, tstop, event)"
	bvar<-NA
	strata<-TRUE
  }
}
if (!is.na(vname.stop)) {
  strata<-TRUE
  bvb<-vlabelV[match(vname.stop,vnameV)]; if (is.na(bvb)) bvb<-vname.stop; 
  bv.lv<-levels(factor(WD[,vname.stop])); 
  bv.lb<-vlabelZ[match(paste(vname.stop,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)
  sxf<-sxf[xvname!=vname.stop];xlv<-xlv[xvname!=vname.stop];xvname<-xvname[xvname!=vname.stop]
}
if (!is.na(bvar)) {
  stvb<-vlabelV[match(bvar,vnameV)]; if (is.na(stvb)) stvb<-bvar; 
  stv.lv<-levels(factor(WD[,bvar]));
  stv.lb<-vlabelZ[match(paste(bvar,stv.lv,sep="."),vnameZ)]
  stv.lb[is.na(stv.lb)]<-stv.lv[is.na(stv.lb)]
  nstg<-length(stv.lv)
  sxf<-sxf[xvname!=bvar];xlv<-xlv[xvname!=bvar];xvname<-xvname[xvname!=bvar]
}
mth<-"efron"; if (!is.na(par1)) mth<-c("efron","breslow","exact")[par1]
colvb=vlabelV[match(colvname,vnameV)]; if (is.na(colvb)) colvb<-colvname;
colvb.utf8<-colvb; Encoding(colvb.utf8)<-"UTF-8"
xb<-vlabelV[match(xvname,vnameV)]; xb[is.na(xb)]<-xvname[is.na(xb)]
xb.utf8<-xb; Encoding(xb.utf8)<-"UTF-8"
timevb<-timevar;
if (!is.na(timevar)) {
 timevb<-vlabelV[match(timevar,vnameV)]; if (is.na(timevb)) timevb<-timevar;
}
if (!is.na(vname.start)) {
 tvb1<-vlabelV[match(vname.start,vnameV)]; if (is.na(tvb1)) tvb1<-vname.start;
 timevb<-paste(tvb1,"-",timevb);
}
timevb.utf8<-timevb; Encoding(timevb.utf8)<-"UTF-8"
if (survStr=="") {
  if (!is.na(vname.start)) {
    fml<-paste("Surv(",vname.start,",",timevar,",",colvname,")~",paste(xv1,collapse="+"))
    fml2<-paste("Surv(",vname.start,",",timevar,",",colvname,")");
  } else {
    if (!is.na(timevar)) fml<-paste("Surv(",timevar,",",colvname,")~",paste(xv1,collapse="+"))
    fml2<-paste("Surv(",timevar,",",colvname,")");  
  }
  survStr<-fml2;
  if (!is.na(vname.stop) & fml>"") {
    fml<-paste(fml,"+strata(",vname.stop,")",sep=""); 
	fml2<-paste(fml2,"~factor(",vname.stop,")")
  } else {fml2<-paste(fml2,"~1");}
  xinstrstr<-""; fmla1<-""
  if (!is.na(bvar)) {
    fmla0<-paste(fml,"+factor(",bvar,")");
    if (sum(sxf=="S")>0) {
      xintrstr<-paste(paste(xv1[sxf=="S"],paste("factor(",bvar,")",sep=""),sep="*"),collapse="+")
      fmla1<-paste(fmla0,"+",xintrstr)
    }
  }
}
w<-c("<!DOCTYPE html><html lang='zh'><head><meta charset='utf-8'></head><body>")
w<-c(w,"<h2>Cox Model</h2>")
sink(paste(ofname,".lst",sep=""))
if (svy==1) summary(svy.DSN)
 

xx1<-c("观察数（人数）","事件发生数")
xx2<-c(" ","回归系数(coef)","回归系数标准误","z","危险比(HR=exp(coef))","1/危险比","95%区间下限","95%区间上限","P　值")
xx3<-"比例灾害模型"
xx4<-"Log Likelihood （对数似然值）="
xx5<-"比例灾害模型诊断"
xxname<-list(xx1,xx2,xx3,xx4,xx5)

w<-c(w,paste("</br>结局变量: ",colvb))
w<-c(w,paste("</br>时间变量: ",timevb,"</br>"))

if (fml>"") {
 if (!is.na(bvar)) {
  for (s in (1:nstg)) {
    WD1<-WD[WD[,bvar]==stv.lv[s],]
    attach(WD1)
    w<-c(w,paste("</br>",stvb,": ",stv.lb[s],"</br>"))
	w<-c(w,runcoxfml(fml,fml2,paste(ofname,bvar,stv.lv[s],sep="_"),xxname))
	detach(WD1)
  }
  w<-c(w,"</br>Total:</br>")
  WD1<-WD;attach(WD1)  
  w<-c(w,runcoxfml(fmla0,fml2,ofname,xxname))
  if (fmla1>"") {
    w<-c(w,"</br>Total: 交互作用模型</br>")
    w<-c(w,runcoxfml(fmla1,"",ofname,xxname))
  }	
 } else {
   WD1<-WD;attach(WD1)  
   w<-c(w,runcoxfml(fml,fml2,ofname,xxname))
 }
}
sink()
w<-c(w,wd.subset)
if (mth == "exact") w<-c(w,"Tier method: exact, HR is OR")
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)