---
title: "From SAS descriptive jobs to R"
author: "John Ehrlinger"
date: today
vignette: >
  %\VignetteIndexEntry{From SAS descriptive jobs to R}
  %\VignetteEngine{quarto::html}
  %\VignetteEncoding{UTF-8}
format:
  html:
    toc: true
---

```{r}
#| label: setup
#| include: false
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
```

## Who this is for

This guide is for a SAS biostatistician who knows the `tp.dc.*` templates and
wants to run the same descriptive work in an hvtiR study. "R equivalent" means
the same numbers when the method is the same, and an explicit difference when
the R job makes a different choice. It does not mean translating every SAS
statement line by line.

## Install once

Install the hvtiR package set, then check it before starting a study job:

```r
pak::pak("ehrlinger/hvtiR")
hvtiR::status()
```

## Open the study first

Open the study's `.Rproj` file in RStudio before creating or rendering a job.
The project establishes the study root; do not change it manually. Attach the
study packages, then confirm that the project has already been adopted and
that its registered data still match the manifest:

```r
library(hvtiRutilities)
library(hvtiRtemplates)

study_root()
study_status()
verify_manifest()
```

The descriptive jobs can be the first jobs in a study. At this stage there may
be no modeled endpoint or analysis type. In the examples below, `cohort` names
the subject being described and `eda` names the stage of work.

## The map

| SAS job | R job | template | engine |
|:--|:--|:--|:--|
| `descriptive/dc.tables.ods.sas` (`%desc_tab`) | `dc-tables` | `10_descriptive/dc-tables.qmd` | `hv_tbl_summary()`; `hv_correlation_table()` and `hv_correlation_matrix()` |
| `descriptive/dc.gfup.sas` | `dc-gfup` | `10_descriptive/dc-gfup.qmd` | `proc_means()` |
| `graphs/dp.trends.R` | `dc-trends` | `10_descriptive/dc-trends.qmd` | `hv_trends()` |
| `tp.dp.DescriptiveSummary.qmd`, `tp.dp.EDA_barplots_scatterplots*.R` | `dc-eda` | `10_descriptive/dc-eda.qmd` | `hv_eda()`, through `hv_eda_pages()` |

## Five steps for every job

1. Scaffold the job. For example,
   `job <- open_job("dc", "cohort", "eda", qualifier = "tables")` writes
   `descriptive/dc.tables.cohort.eda.qmd` in an adopted legacy study and opens
   it in the editor. `study_dir("descriptive")` resolves the actual directory,
   including `10_descriptive/` in a numbered-layout study.
2. Work the `EDIT:` markers, rendering `render_job(job)` as you go. Unfinished
   work renders as a draft, as far as the markers you have worked.
3. Stop when every marker is gone and the draft banner with it. Render the
   accepted result with `render_job(job, final = TRUE)`, which stops instead
   of drafting if a marker remains.
4. Compare the result with the SAS `.lst`, section by section.
5. Keep the authored QMD and its HTML report flat in the job's working
   directory: `descriptive/` for the computation and EDA jobs, or
   `graphs/` for trends. The tables job writes its Word file under
   `documents/cohort-eda/`; the plot jobs write PNG files under
   `graphs/cohort-eda/`. `study_dir()` resolves each directory in either study
   layout.

The filename records the choices in the call. `cohort-eda` is the shared job
set, `dc` is the descriptive-computation prefix, and `tables` selects the table
template within that prefix. These are organizational labels. When no outcome
or model has been chosen, do not invent one merely to name an EDA job.

## The draft banner and `EDIT:` markers

Every line a study author has to change carries an `EDIT:` marker. Markers
left in a job do not stop its render; while any remain, the report opens with a
banner:

```
DRAFT -- this job is unfinished

Unresolved markers remain. The numbers below are not a result.

5 unresolved EDIT: marker(s) remain in this job:
  - <!-- EDIT: name the job this replaces (descriptive/dc.tables.*) and say whether it
  - # EDIT: the variables to report, grouped. Each name is a section heading and
  ...
```

The render log carries the same list as a warning. The banner says the
numbers are not a result because some of them may still come from the
template's placeholder choices rather than your study.

The banner does not mean every section will run. A section whose marker is not
yet worked still reads the template's placeholders, and most of those stop with
their own error rather than draw a table from nothing. A new `dc-tables` job,
for instance, stops with "Unknown table variable(s)" until `GROUPS` names your
columns, and before that if the `eda` analysis set has not been written. Work
the markers from the top of the file down, and each render gets further.

Each listed line is the start of a marker. Read the comment around it, make the
change it asks for, delete the marker line, and render again. To see what is
left from the RStudio Console:

```r
grep("EDIT:", readLines(job), value = TRUE)
```

The guard scans the whole file, so a marker in the prose, inside an HTML
comment, counts as much as one in a code chunk. When the last marker goes, so
does the banner.

For a final render, such as the one that goes to an investigator, make an
unfinished job stop rather than render as a draft:

```r
render_job(job, final = TRUE)
```

`render_job()` runs Quarto from the job's own directory, so the same call works
for adopted and numbered study layouts. A final render fails at the
`setup` chunk with the unresolved list.

## Which data a job reads

A SAS job names its dataset in a `set` statement, so the answer is in the code.
An R job takes it from `_study.yml`, which can offer the same study more than
one way. The descriptive jobs choose with two settings in their
`edit-study-choices` chunk, and every rendered report prints what they chose:

```
Data read: dataset `builtr` (builtr_20261007.parquet), 2514 rows, 135 columns.
```

The file named is the dated version registering `builtr.sas7bdat` made, which
is what R reads.

`_study.yml` can declare three kinds of data:

```yaml
built: built.sas7bdat            # the study dataset
additional_datasets:
  builtr:                        # another file, read as it stands
    built: builtr.sas7bdat
analysis_sets:
  eda:                           # columns kept and rows excluded, written by R
    ...
```

| to read | `DATASET` | `ANALYSIS_SET` | exclusions shown |
|:--|:--|:--|:--|
| the study dataset | `"built"` | `NULL` | none applied |
| an additional dataset, such as a column subset written for R | its name, `"builtr"` | `NULL` | none applied |
| an analysis set | `"built"` | its name, `"eda"` | yes, in order |

Three things catch people:

- **An additional dataset is not an analysis set.** Putting `"builtr"` in
  `ANALYSIS_SET` stops with "No analysis set `builtr`". The two are separate
  lists in `_study.yml`.
- **An analysis set must be written before a job can read it.** Run
  `hvtiRdatabuild::write_analysis_set("eda")` once from the study root. If the
  study dataset or the set's declaration changes afterwards, the job stops and
  names that same call, so it never describes a cohort that has since moved.
- **An analysis set is always cut from the study dataset.** Pairing one with a
  `DATASET` other than `"built"`, or its older name `"study"`, stops rather
  than quietly using the study dataset instead.

Before any of this, `verify_manifest()` checks every file's checksum against
`manifest.yaml`. A stop on that first line means a registered file changed, not
that a setting is wrong. A dataset rebuilt since it was registered does not
stop the job: it reads the registered version, `verify_manifest()` lists the
rebuilt file as `PENDING`, and the report says a newer file is waiting, with
the `hvtiRutilities::update_manifest()` call that registers it.

## Descriptive tables and correlations: `dc-tables`

The SAS job commonly calls `%desc_tab` once for categorical variables and once
for continuous variables. The R job supplies both sets to one table. SAS
comment banners become the names in `GROUPS`; `by=` and `byvalue=` become `BY`;
and `perctile=15 50 85` becomes `continuous_stat = "both"` with percentiles 15
and 85. A `by=` group gets a p-value, `COMPARE = "pvalue"`, as `%desc_tab`
printed one. Three things differ on purpose: a categorical variable with sparse
cells is tested with Fisher's exact test, where `%desc_tab` always used
chi-square; two groups of a continuous variable are compared with the Wilcoxon
rank-sum test, exact in small groups, where `%desc_tab` reported the
Kruskal-Wallis chi-square approximation; and the low and high percentiles use
SAS's default definition rather than the `PCTLDEF=4` `%desc_tab` set, so they
can differ in small groups.
`migrate_job()` writes these choices and lists the differences in its report.
The report shows the table in the report style, with each variable's column
name; the Word file holds it in the manuscript style, under `MANUSCRIPT_CAPTION`
and shaded in `MANUSCRIPT_FILL`. A migrated job that shows both statistics gets
`MANUSCRIPT_LAYOUT <- "describe"`, the manuscript layout that draws them.
The `dc.tables.ods_<topic>.sas` correlation variant becomes `CORR` in the same
job.

```sas
%desc_tab(vartype=category,   varlist=female diabetes nyha_pr);
%desc_tab(vartype=continuous, varlist=age bmi, perctile=15 50 85);
```

```{r dc-tables, eval=requireNamespace("hvtiRtables", quietly = TRUE) && requireNamespace("hvtiPlotR", quietly = TRUE)}
d <- hvtiPlotR::sample_correlation_data(n = 80)
d$female <- rep(0:1, length.out = nrow(d))
groups <- list(
  Demography = c("female"),
  Measurements = c("a1c", "glucose", "creatinine")
)
hvtiRtables::hv_tbl_summary(
  d,
  groups = groups,
  binary = "female",
  continuous = c("a1c", "glucose", "creatinine"),
  continuous_stat = "both",
  percentiles = c(15, 85),
  compare = "none"
)
```

```{r dc-correlations, eval=requireNamespace("hvtiRtables", quietly = TRUE) && requireNamespace("hvtiPlotR", quietly = TRUE)}
hvtiRtables::hv_correlation_table(
  d,
  vars = c("glucose", "creatinine"),
  with = "a1c"
)
```

## Goodness of follow-up: `dc-gfup`

The SAS `proc means` and `proc univariate` calls become `proc_means()`. The R
job reports follow-up distributions overall, by vital status, and among
survivors. It deliberately omits the patient-level listing: a rendered report
should answer the completeness question without exposing an identifier.

```sas
proc means data=built n mean std min max;
  class dead;
  var iv_dead iv_fup;
run;
```

```{r dc-gfup}
if (requireNamespace("hvtiRutilities", quietly = TRUE)) {
  d <- hvtiRutilities::generate_survival_data(n = 80)
  d$iv_fup <- pmin(d$iv_dead, 0.85 * d$iv_dead)
  hvtiRutilities::proc_means(
    d,
    vars = c("iv_dead", "iv_fup"),
    class = "dead",
    stats = c("n", "nmiss", "mean", "std", "min", "p15", "median", "p85", "max")
  )
}
```

## Trends over operation year: `dc-trends`

The shipped `dc-trends` template replaces hand-built annual summaries and
per-series `smooth.spline(df = 6)` calls with `hv_trends()`. Its smoother is
LOESS, not the SAS-era smoothing spline, so compare the old and new figures by
shape rather than pointwise equality.

```r
fit <- smooth.spline(year, annual_mean, df = 6)
```

```{r dc-trends}
if (requireNamespace("hvtiPlotR", quietly = TRUE)) {
  d <- hvtiPlotR::sample_eda_data(n = 120)
  trend_data <- data.frame(year = d$year, value = d$ef, group = "EF")
  plot(hvtiPlotR::hv_trends(trend_data))
}
```

## EDA report: `dc-eda`

The older `Function_DataPlotting()` loops become one `hv_eda()` call per entry
in `VARS`. Character and factor columns produce bars. A numeric column produces
bars only when it has at most six distinct non-missing values and every value
is a non-negative whole number no greater than six; otherwise it produces a
scatter panel. The template lays those panels out in printable pages.

A legacy EDA report drew no follow-up panels, so `migrate_job()` writes its
`dc-eda` job with `SECTIONS <- c("continuous", "percent", "count")`: the
continuous, percent and count pages, each with its table. Add `"followup"` to
draw the follow-up panels as well.

```r
for (v in variables) Function_DataPlotting(data, v, "year")
```

```{r dc-eda, eval=requireNamespace("hvtiPlotR", quietly = TRUE)}
d <- hvtiPlotR::sample_eda_data(n = 120)
variables <- c(male = "Male", ef = "Ejection fraction")
panels <- lapply(names(variables), function(v) {
  plot(hvtiPlotR::hv_eda(d, x_col = "year", y_col = v, y_label = variables[[v]]))
})
panels[[1]]
panels[[2]]
```

## Where the numbers differ

Percentiles use SAS `QNTLDEF=5`, equivalent to R's `type = 2`;
`proc_means()` and the table engine already make that choice. Pin `type = 2`
in any hand-written `quantile()` call. `N` is the non-missing count, not always
the number of rows. Continuous-variable p-values are non-parametric.
Correlation intervals are 68% Fisher intervals by default. Finally,
`%desc_tab` wrote separate categorical and continuous RTF files; the R job
writes one DOCX table.

## What is not ported

The templates do not reproduce RTF output; they produce DOCX and self-contained
HTML instead. They also do not reproduce identified patient listings. Use the
distribution summaries to assess completeness, then investigate any necessary
patient-level records inside the governed study environment rather than in a
portable report.