R/epitool_core.R

In R4VN: Health Data Analysis and Publication-Ready Reporting

# =============================================================================
# R4VN EpiTool Studio - Core epidemiology engines
# File: R/epitool_core.R
# =============================================================================

.r4vn_epi_num <- function(x, digits = 2L) {
  ifelse(is.finite(x), formatC(x, format = "f", digits = digits), NA_character_)
}

.r4vn_epi_pct <- function(x, digits = 1L) {
  ifelse(is.finite(x), paste0(formatC(100 * x, format = "f", digits = digits), "%"), NA_character_)
}

.r4vn_epi_ci_text <- function(est, low, high, digits = 2L) {
  paste0(.r4vn_epi_num(est, digits), " (95% CI ",
         .r4vn_epi_num(low, digits), "\u2013", .r4vn_epi_num(high, digits), ")")
}

.r4vn_epi_safe_date <- function(x) {
  if (inherits(x, "Date")) return(x)
  if (inherits(x, c("POSIXct", "POSIXlt"))) return(as.Date(x))
  if (is.numeric(x)) {
    # Excel dates are common in field files.
    z <- suppressWarnings(as.Date(x, origin = "1899-12-30"))
    if (mean(!is.na(z)) > .8) return(z)
  }
  ch <- trimws(as.character(x))
  out <- suppressWarnings(as.Date(ch))
  bad <- is.na(out) & nzchar(ch)
  if (any(bad)) {
    formats <- c("%d/%m/%Y", "%m/%d/%Y", "%d-%m-%Y", "%Y/%m/%d",
                 "%d.%m.%Y", "%d %b %Y", "%d %B %Y")
    for (fmt in formats) {
      ii <- which(bad)
      if (!length(ii)) break
      z <- suppressWarnings(as.Date(ch[ii], format = fmt))
      ok <- !is.na(z)
      out[ii[ok]] <- z[ok]
      bad[ii[ok]] <- FALSE
    }
  }
  out
}

.r4vn_epi_safe_datetime <- function(x, tz = "UTC") {
  if (inherits(x, "POSIXct")) return(x)
  if (inherits(x, "POSIXlt")) return(as.POSIXct(x))
  if (inherits(x, "Date")) return(as.POSIXct(x, tz = tz))
  ch <- trimws(as.character(x))
  suppressWarnings(as.POSIXct(ch, tz = tz))
}

.r4vn_epi_var_label <- function(x, fallback = "") {
  z <- attr(x, "label", exact = TRUE)
  if (is.null(z) || !length(z) || is.na(z[1]) || !nzchar(as.character(z[1]))) fallback else as.character(z[1])
}

.r4vn_epi_positive <- function(x, event = NULL) {
  if (is.logical(x)) return(!is.na(x) & x)
  if (!is.null(event) && length(event) && !is.na(event[1]) && nzchar(as.character(event[1]))) {
    return(!is.na(x) & as.character(x) == as.character(event[1]))
  }
  if (is.numeric(x) || is.integer(x)) return(!is.na(x) & x == 1)
  ch <- trimws(tolower(as.character(x)))
  yes <- c("1","yes","y","true","positive","pos","exposed","case","disease",
           "ill","confirmed","co","c\u00f3","duong tinh","d\u01b0\u01a1ng t\u00ednh","benh","b\u1ec7nh")
  hit <- ch %in% yes
  if (any(hit, na.rm = TRUE)) return(!is.na(ch) & hit)
  lv <- unique(ch[!is.na(ch) & nzchar(ch)])
  if (length(lv) == 2L) return(!is.na(ch) & ch == lv[2L])
  rep(FALSE, length(x))
}

.r4vn_epi_binom_ci <- function(x, n, conf.level = .95) {
  if (!is.finite(x) || !is.finite(n) || n <= 0 || x < 0 || x > n) return(c(NA_real_, NA_real_))
  as.numeric(stats::binom.test(x, n, conf.level = conf.level)$conf.int[1:2])
}

.r4vn_epi_2x2 <- function(a, b, c, d, conf.level = .95, correction = .5) {
  vals <- c(a = a, b = b, c = c, d = d)
  if (any(!is.finite(vals)) || any(vals < 0)) {
    stop("All four cell counts must be non-negative finite numbers.", call. = FALSE)
  }
  if (sum(vals) <= 0) stop("The 2 x 2 table is empty.", call. = FALSE)

  risk_e <- if ((a+b) > 0) a/(a+b) else NA_real_
  risk_u <- if ((c+d) > 0) c/(c+d) else NA_real_
  rd <- risk_e-risk_u
  z <- stats::qnorm(1-(1-conf.level)/2)

  corrected <- any(vals == 0)
  cc <- if (corrected) correction else 0
  aa <- a+cc; bb <- b+cc; cc2 <- c+cc; dd <- d+cc

  rr <- (aa/(aa+bb))/(cc2/(cc2+dd))
  se_rr <- sqrt(1/aa - 1/(aa+bb) + 1/cc2 - 1/(cc2+dd))
  rr_ci <- exp(log(rr) + c(-1,1)*z*se_rr)

  or <- aa*dd/(bb*cc2)
  se_or <- sqrt(1/aa + 1/bb + 1/cc2 + 1/dd)
  or_ci <- exp(log(or) + c(-1,1)*z*se_or)

  se_rd <- sqrt(risk_e*(1-risk_e)/max(a+b,1) + risk_u*(1-risk_u)/max(c+d,1))
  rd_ci <- pmax(-1, pmin(1, rd+c(-1,1)*z*se_rd))

  mat <- matrix(c(a,b,c,d), nrow = 2, byrow = TRUE)
  fisher <- tryCatch(stats::fisher.test(mat), error=function(e) NULL)
  chi <- tryCatch(suppressWarnings(stats::chisq.test(mat, correct=FALSE)), error=function(e) NULL)

  list(
    cells=vals, risk_exposed=risk_e, risk_unexposed=risk_u,
    risk_difference=rd, risk_difference_ci=rd_ci,
    rr=rr, rr_ci=rr_ci, or=or, or_ci=or_ci,
    attributable_fraction_exposed=if (is.finite(rr) && rr != 0) (rr-1)/rr else NA_real_,
    vaccine_effectiveness=if (is.finite(rr)) 1-rr else NA_real_,
    nnt=if (is.finite(rd) && rd != 0) 1/abs(rd) else Inf,
    fisher_p=if (is.null(fisher)) NA_real_ else fisher$p.value,
    chisq_p=if (is.null(chi)) NA_real_ else chi$p.value,
    corrected=corrected, correction=correction, conf.level=conf.level
  )
}

.r4vn_epi_2x2_table <- function(x) {
  data.frame(
    Measure = c("Risk among exposed","Risk among unexposed","Risk difference",
                "Risk ratio","Odds ratio","Attributable fraction among exposed",
                "Protective / vaccine effectiveness","NNT / NNH",
                "Fisher exact p-value","Pearson chi-square p-value"),
    Estimate = c(
      .r4vn_epi_pct(x$risk_exposed),
      .r4vn_epi_pct(x$risk_unexposed),
      paste0(.r4vn_epi_pct(x$risk_difference)," (95% CI ",
             .r4vn_epi_pct(x$risk_difference_ci[1])," to ",
             .r4vn_epi_pct(x$risk_difference_ci[2]),")"),
      .r4vn_epi_ci_text(x$rr,x$rr_ci[1],x$rr_ci[2]),
      .r4vn_epi_ci_text(x$or,x$or_ci[1],x$or_ci[2]),
      .r4vn_epi_pct(x$attributable_fraction_exposed),
      .r4vn_epi_pct(x$vaccine_effectiveness),
      if (is.finite(x$nnt)) .r4vn_epi_num(x$nnt,1) else "Inf",
      if (is.finite(x$fisher_p)) format.pval(x$fisher_p,digits=3,eps=.001) else NA,
      if (is.finite(x$chisq_p)) format.pval(x$chisq_p,digits=3,eps=.001) else NA
    ),
    check.names=FALSE, stringsAsFactors=FALSE
  )
}

.r4vn_epi_measure <- function(type, numerator, denominator, multiplier=100, conf.level=.95) {
  type <- match.arg(type,c("prevalence","risk","attack_rate","case_fatality","rate"))
  if (!is.finite(numerator) || !is.finite(denominator) || numerator < 0 || denominator <= 0) {
    stop("Numerator must be non-negative and denominator must be positive.", call.=FALSE)
  }
  if(type != "rate" && numerator > denominator) {
    stop("For a proportion-based measure, numerator cannot exceed denominator.",call.=FALSE)
  }
  est <- numerator/denominator
  if (type=="rate") {
    z <- stats::qnorm(1-(1-conf.level)/2)
    se <- sqrt(numerator)/denominator
    ci <- pmax(0, est+c(-1,1)*z*se)
  } else ci <- .r4vn_epi_binom_ci(numerator, denominator, conf.level)
  labels <- c(prevalence="Prevalence", risk="Cumulative incidence / risk",
              attack_rate="Attack rate", case_fatality="Case fatality ratio",
              rate="Incidence / mortality rate")
  list(type=type,label=unname(labels[type]),estimate=est*multiplier,
       ci=ci*multiplier,raw=est,multiplier=multiplier,numerator=numerator,
       denominator=denominator,conf.level=conf.level)
}

.r4vn_epi_diag <- function(tp, fp, fn, tn, conf.level=.95) {
  vals <- c(tp=tp,fp=fp,fn=fn,tn=tn)
  if (any(!is.finite(vals)) || any(vals < 0)) stop("Diagnostic cell counts must be non-negative.",call.=FALSE)
  metric <- function(num,den) {
    est <- if (den>0) num/den else NA_real_
    ci <- .r4vn_epi_binom_ci(num,den,conf.level)
    c(est=est,low=ci[1],high=ci[2])
  }
  se <- metric(tp,tp+fn); sp <- metric(tn,tn+fp)
  ppv <- metric(tp,tp+fp); npv <- metric(tn,tn+fn)
  acc <- metric(tp+tn,sum(vals)); prev <- metric(tp+fn,sum(vals))
  lr_pos <- if (is.finite(se[1]) && is.finite(sp[1]) && sp[1] < 1) se[1]/(1-sp[1]) else Inf
  lr_neg <- if (is.finite(se[1]) && is.finite(sp[1]) && sp[1] > 0) (1-se[1])/sp[1] else Inf
  list(cells=vals,sensitivity=se,specificity=sp,ppv=ppv,npv=npv,
       accuracy=acc,prevalence=prev,lr_pos=lr_pos,lr_neg=lr_neg,
       diagnostic_or=if (is.finite(lr_pos)&&is.finite(lr_neg)&&lr_neg!=0) lr_pos/lr_neg else Inf,
       youden=se[1]+sp[1]-1,conf.level=conf.level)
}

.r4vn_epi_diag_table <- function(x) {
  pctci <- function(z) paste0(.r4vn_epi_pct(z[1])," (95% CI ",.r4vn_epi_pct(z[2]),"\u2013",.r4vn_epi_pct(z[3]),")")
  data.frame(
    Measure=c("Sensitivity","Specificity","Positive predictive value","Negative predictive value",
              "Accuracy","Observed prevalence","LR+","LR-","Diagnostic OR","Youden index"),
    Estimate=c(pctci(x$sensitivity),pctci(x$specificity),pctci(x$ppv),pctci(x$npv),
               pctci(x$accuracy),pctci(x$prevalence),.r4vn_epi_num(x$lr_pos),
               .r4vn_epi_num(x$lr_neg),.r4vn_epi_num(x$diagnostic_or),.r4vn_epi_num(x$youden)),
    stringsAsFactors=FALSE
  )
}

.r4vn_epi_profile <- function(data) {
  if (is.null(data) || !is.data.frame(data)) {
    return(list(n=0L,p=0L,dates=character(),binary=character(),categorical=character(),
                numeric=character(),case_candidates=character(),geo_candidates=character(),
                id_candidates=character(),lat_candidates=character(),lon_candidates=character(),
                missing=NA_real_))
  }
  nms <- names(data); low <- tolower(nms)
  nlev <- vapply(data,function(x) length(unique(x[!is.na(x)])),integer(1))
  pick <- function(pattern) nms[grepl(pattern,low,perl=TRUE)]
  list(
    n=nrow(data), p=ncol(data),
    dates=nms[vapply(data,function(x) inherits(x,c("Date","POSIXct","POSIXlt")),logical(1)) |
              grepl("date|ngay|ng\u00e0y|onset|time",low)],
    binary=nms[nlev==2L],
    categorical=nms[vapply(data,function(x) is.factor(x)||is.character(x)||is.logical(x),logical(1)) & nlev<=30],
    numeric=nms[vapply(data,is.numeric,logical(1))],
    case_candidates=pick("case|status|benh|b\u1ec7nh|ill|outcome|disease|confirmed|probable|suspected"),
    geo_candidates=pick("province|district|commune|ward|site|location|tinh|t\u1ec9nh|huyen|huy\u1ec7n|xa$|x\u00e3$"),
    id_candidates=pick("^id$|case.?id|patient.?id|record.?id|stt|code|ma.?ca|m\u00e3.?ca"),
    lat_candidates=pick("^lat$|latitude|vi.?do|v\u0129.?do"),
    lon_candidates=pick("^lon$|^lng$|longitude|kinh.?do"),
    missing=sum(is.na(data))/max(1,nrow(data)*ncol(data))
  )
}

.r4vn_epi_health <- function(data,id=NULL,onset=NULL,report=NULL,age=NULL) {
  if (is.null(data)||!is.data.frame(data)) stop("A data frame is required.",call.=FALSE)
  n <- nrow(data)
  issues <- data.frame(Issue=character(),Severity=character(),N=integer(),Percent=numeric(),stringsAsFactors=FALSE)
  add <- function(issue,severity,count) {
    if (!is.finite(count)||count<=0) return(NULL)
    data.frame(Issue=issue,Severity=severity,N=as.integer(count),
               Percent=if(n>0)100*count/n else NA_real_,stringsAsFactors=FALSE)
  }
  push <- function(z) if(!is.null(z)) issues <<- rbind(issues,z)
  push(add("Records with at least one missing value","Info",sum(!stats::complete.cases(data))))
  if (!is.null(id)&&nzchar(id)&&id%in%names(data)) {
    x <- data[[id]]
    push(add("Duplicate IDs","Critical",sum(!is.na(x)&duplicated(x))))
    push(add("Missing IDs","Critical",sum(is.na(x)|trimws(as.character(x))=="")))
  }
  onset_date <- NULL
  if (!is.null(onset)&&nzchar(onset)&&onset%in%names(data)) {
    onset_date <- .r4vn_epi_safe_date(data[[onset]])
    push(add("Missing onset dates","Warning",sum(is.na(onset_date))))
    push(add("Onset dates in the future","Critical",sum(!is.na(onset_date)&onset_date>Sys.Date())))
  }
  if (!is.null(report)&&nzchar(report)&&report%in%names(data)) {
    report_date <- .r4vn_epi_safe_date(data[[report]])
    push(add("Missing report dates","Warning",sum(is.na(report_date))))
    push(add("Report dates in the future","Warning",sum(!is.na(report_date)&report_date>Sys.Date())))
    if (!is.null(onset_date)) push(add("Report date before onset date","Critical",
      sum(!is.na(onset_date)&!is.na(report_date)&report_date<onset_date)))
  }
  if (!is.null(age)&&nzchar(age)&&age%in%names(data)) {
    z <- suppressWarnings(as.numeric(data[[age]]))
    push(add("Impossible age (<0 or >120)","Critical",sum(!is.na(z)&(z<0|z>120))))
  }
  if (!nrow(issues)) issues <- data.frame(Issue="No rule-based issues detected",Severity="OK",N=0L,Percent=0)
  penalty <- sum(ifelse(issues$Severity=="Critical",pmin(issues$Percent,20),
                 ifelse(issues$Severity=="Warning",pmin(issues$Percent,10)/2,0)),na.rm=TRUE)
  list(table=issues,readiness=max(0,100-penalty),
       critical=sum(issues$Severity=="Critical"),warning=sum(issues$Severity=="Warning"))
}

.r4vn_epi_bucket <- function(x, interval=c("day","week","month","hour")) {
  interval <- match.arg(interval)
  if (interval=="hour") {
    dt <- .r4vn_epi_safe_datetime(x)
    return(as.POSIXct(cut(dt,breaks="hour"),tz="UTC"))
  }
  d <- .r4vn_epi_safe_date(x)
  if (interval=="day") return(d)
  if (interval=="month") return(as.Date(format(d,"%Y-%m-01")))
  as.Date(cut(d,breaks="week",start.on.monday=TRUE))
}

.r4vn_epi_curve <- function(data,onset,group=NULL,interval="day") {
  if (!onset%in%names(data)) stop("Onset variable not found.",call.=FALSE)
  b <- .r4vn_epi_bucket(data[[onset]],interval)
  keep <- !is.na(b)
  if (!any(keep)) stop("No usable onset dates were found.",call.=FALSE)
  if (is.null(group)||!nzchar(group)||!group%in%names(data)) {
    tab <- as.data.frame(table(b[keep]),stringsAsFactors=FALSE)
    names(tab) <- c("Time","Cases"); tab$Time <- as.character(tab$Time)
    tab$Cases <- as.integer(tab$Cases); tab$Group <- "Cases"
  } else {
    g <- as.character(data[[group]]); g[is.na(g)|!nzchar(g)] <- "Missing"
    tab <- as.data.frame(table(b[keep],g[keep]),stringsAsFactors=FALSE)
    names(tab) <- c("Time","Group","Cases"); tab <- tab[tab$Cases>0,,drop=FALSE]
    tab$Time <- as.character(tab$Time); tab$Cases <- as.integer(tab$Cases)
  }
  tab
}

.r4vn_epi_exposure_screen <- function(data,outcome,exposures,event=NULL,exposed_event=NULL,conf.level=.95) {
  if (!outcome%in%names(data)) stop("Outcome variable not found.",call.=FALSE)
  exposures <- intersect(exposures,names(data))
  if (!length(exposures)) stop("Select at least one exposure variable.",call.=FALSE)
  yraw <- data[[outcome]]; y <- .r4vn_epi_positive(yraw,event)
  rows <- lapply(exposures,function(v) {
    xraw <- data[[v]]; x <- .r4vn_epi_positive(xraw,exposed_event)
    ok <- !is.na(yraw)&!is.na(xraw)
    a <- sum(ok&x&y); b <- sum(ok&x&!y); c <- sum(ok&!x&y); d <- sum(ok&!x&!y)
    fit <- tryCatch(.r4vn_epi_2x2(a,b,c,d,conf.level),error=function(e)NULL)
    if (is.null(fit)) return(NULL)
    data.frame(Exposure=.r4vn_epi_var_label(xraw,v),Variable=v,
      Exposed_cases=a,Exposed_total=a+b,AR_exposed=fit$risk_exposed,
      Unexposed_cases=c,Unexposed_total=c+d,AR_unexposed=fit$risk_unexposed,
      RR=fit$rr,RR_low=fit$rr_ci[1],RR_high=fit$rr_ci[2],
      OR=fit$or,OR_low=fit$or_ci[1],OR_high=fit$or_ci[2],
      p=fit$fisher_p,Corrected=fit$corrected,stringsAsFactors=FALSE)
  })
  rows <- Filter(Negate(is.null),rows)
  if (!length(rows)) stop("No analyzable exposure tables were produced.",call.=FALSE)
  out <- do.call(rbind,rows); out <- out[order(out$RR,decreasing=TRUE,na.last=TRUE),,drop=FALSE]
  rownames(out) <- NULL; out
}

.r4vn_epi_stratified <- function(data,outcome,exposure,strata,event=NULL,exposed_event=NULL,conf.level=.95) {
  needed <- c(outcome,exposure,strata)
  if (any(!needed%in%names(data))) stop("Outcome, exposure, or strata variable not found.",call.=FALSE)
  yraw <- data[[outcome]]; xraw <- data[[exposure]]; sraw <- data[[strata]]
  y <- .r4vn_epi_positive(yraw,event); x <- .r4vn_epi_positive(xraw,exposed_event)
  ok <- !is.na(yraw)&!is.na(xraw)&!is.na(sraw)
  s <- droplevels(as.factor(sraw[ok])); y <- y[ok]; x <- x[ok]
  lev <- levels(s)
  arr <- array(0,dim=c(2,2,length(lev)),dimnames=list(Exposure=c("No","Yes"),Outcome=c("No","Yes"),Stratum=lev))
  rows <- lapply(seq_along(lev),function(i) {
    z <- s==lev[i]; a<-sum(z&x&y); b<-sum(z&x&!y); c<-sum(z&!x&y); d<-sum(z&!x&!y)
    arr[,,i] <<- matrix(c(d,b,c,a),nrow=2,byrow=FALSE)
    fit <- tryCatch(.r4vn_epi_2x2(a,b,c,d,conf.level),error=function(e)NULL)
    if(is.null(fit)) return(data.frame(Stratum=lev[i],a=a,b=b,c=c,d=d,RR=NA,RR_low=NA,RR_high=NA,OR=NA,OR_low=NA,OR_high=NA))
    data.frame(Stratum=lev[i],a=a,b=b,c=c,d=d,RR=fit$rr,RR_low=fit$rr_ci[1],RR_high=fit$rr_ci[2],OR=fit$or,OR_low=fit$or_ci[1],OR_high=fit$or_ci[2])
  })
  mh <- tryCatch(stats::mantelhaen.test(arr,conf.level=conf.level),error=function(e)NULL)
  crude <- .r4vn_epi_2x2(sum(x&y),sum(x&!y),sum(!x&y),sum(!x&!y),conf.level)
  list(strata=do.call(rbind,rows),crude=crude,
       mh_or=if(is.null(mh))NA_real_ else unname(mh$estimate),
       mh_ci=if(is.null(mh))c(NA,NA) else as.numeric(mh$conf.int[1:2]),
       mh_p=if(is.null(mh))NA_real_ else mh$p.value,conf.level=conf.level)
}

.r4vn_epi_surveillance <- function(data,date,count=NULL,interval=c("week","month"),baseline_years=3L,k=2) {
  interval <- match.arg(interval)
  if(!date%in%names(data)) stop("Date variable not found.",call.=FALSE)
  d <- .r4vn_epi_safe_date(data[[date]]); ok <- !is.na(d)
  if(!any(ok)) stop("No usable dates were found.",call.=FALSE)
  value <- if(is.null(count)||!nzchar(count)||!count%in%names(data)) rep(1,nrow(data)) else suppressWarnings(as.numeric(data[[count]]))
  value[!is.finite(value)] <- 0
  year <- if(interval=="month") as.integer(format(d,"%Y")) else as.integer(format(d,"%G"))
  period <- if(interval=="month") as.integer(format(d,"%m")) else as.integer(format(d,"%V"))
  agg <- stats::aggregate(value[ok],list(year=year[ok],period=period[ok]),sum)
  names(agg)[3] <- "cases"; agg <- agg[order(agg$year,agg$period),]
  agg$baseline_mean <- agg$baseline_sd <- agg$threshold <- NA_real_
  agg$signal <- FALSE
  for(i in seq_len(nrow(agg))) {
    hist <- agg$cases[agg$period==agg$period[i] & agg$year<agg$year[i] & agg$year>=agg$year[i]-baseline_years]
    if(length(hist)>=2) {
      mu <- mean(hist); sdv <- stats::sd(hist)
      agg$baseline_mean[i] <- mu; agg$baseline_sd[i] <- sdv
      agg$threshold[i] <- mu+k*sdv
      agg$signal[i] <- is.finite(agg$threshold[i]) && agg$cases[i] > agg$threshold[i]
    }
  }
  agg
}

.r4vn_epi_ss_oneprop <- function(p=.5,d=.05,confidence=.95,deff=1,finite_population=NULL) {
  z <- stats::qnorm(1-(1-confidence)/2)
  n <- z^2*p*(1-p)/d^2*deff
  if(!is.null(finite_population)&&is.finite(finite_population)&&finite_population>0) {
    n <- n/(1+(n-1)/finite_population)
  }
  ceiling(n)
}

.r4vn_epi_ss_twoprop <- function(p1,p2,power=.8,alpha=.05,ratio=1) {
  if(any(!is.finite(c(p1,p2,power,alpha,ratio)))||p1<=0||p1>=1||p2<=0||p2>=1||ratio<=0) {
    stop("Invalid sample-size inputs.",call.=FALSE)
  }
  z1 <- stats::qnorm(1-alpha/2); z2 <- stats::qnorm(power)
  pbar <- (p1+ratio*p2)/(1+ratio)
  n1 <- ((z1*sqrt((1+1/ratio)*pbar*(1-pbar)) +
          z2*sqrt(p1*(1-p1)+p2*(1-p2)/ratio))^2)/(p1-p2)^2
  c(group1=ceiling(n1),group2=ceiling(n1*ratio))
}

.r4vn_epi_html_table <- function(x,digits=3) {
  if(is.null(x)||!NROW(x)) return("<p>No results.</p>")
  y <- as.data.frame(x,stringsAsFactors=FALSE)
  y[] <- lapply(y,function(z) if(is.numeric(z)) formatC(z,format="fg",digits=digits) else as.character(z))
  esc <- function(s) {
    s <- gsub("&","&amp;",s,fixed=TRUE); s <- gsub("<","&lt;",s,fixed=TRUE)
    s <- gsub(">","&gt;",s,fixed=TRUE); s
  }
  h <- paste0("<table class='epi-table'><thead><tr>",paste0("<th>",esc(names(y)),"</th>",collapse=""),"</tr></thead><tbody>")
  for(i in seq_len(nrow(y))) h <- paste0(h,"<tr>",paste0("<td>",esc(as.character(y[i,])),"</td>",collapse=""),"</tr>")
  paste0(h,"</tbody></table>")
}