From Analysis Script to Rendered Report
A GLMM analysis may run correctly in an analyst’s interactive R session and still fail when someone else tries to reproduce it. The cause may be an object created earlier but not included in the report, a file path that only works on one computer, or a figure that depends on an external resource. A Quarto report addresses this by bringing explanatory text, executable R code, and selected outputs into one document that can be rendered from a clean session.
The goal is not to print every line of code or every item in a model summary. A useful report explains the analysis, runs the necessary preparation and fitting steps in a clear order, and presents a focused set of results. A self-contained HTML file can embed supporting resources such as figures and styling, but it does not by itself guarantee that the analysis is reproducible or that private source data are safe to share.
The report should state the analysis population and relevant data-processing decisions, then show the model and the selected results. It should also distinguish descriptive summaries from model-based estimates. The earlier material on GLMM interpretation and reporting provides the statistical context; here the emphasis is on implementing that reporting workflow reliably.
Set Up the Document and Its Execution Order
A Quarto source file commonly has a .qmd extension. It begins with YAML metadata between lines containing three hyphens. For an HTML report, the metadata can specify the title and output format. The following settings request embedded resources:
---
title: "GLMM Analysis"
format:
html:
embed-resources: true
---
The document body can then be divided into sections such as “Data and analysis sample,” “Model,” “Results,” and “Limitations.” Executable R chunks are delimited by lines beginning with three backticks and an R chunk header, such as ```{r}. A chunk may have options on lines beginning with #|. For example, #| label: setup gives a chunk a meaningful label. Labels make render errors easier to locate and help keep the source organized.
Chunk order matters: a chunk can use an object only after an earlier chunk has created it. Put data import and preparation before model fitting, and model fitting before tables or plots that use the fitted model. Avoid relying on objects left over in an interactive workspace. Rendering should work when the report is executed in a fresh R session with only its declared inputs and installed packages available.
Worked Example: Create a Minimal Report Structure
Suppose a hypothetical analysis uses a comma-separated file called readmissions.csv. The file contains one row per patient, with readmission status, treatment assignment, age, and ward. The report should read the file, prepare variables, fit a model, and then present selected output.
Start the source with the YAML metadata shown above. Add an introductory paragraph describing the hypothetical population and the outcome, then add a chunk labelled setup. In that chunk, load the packages needed by the report and import the file using a project-relative path: dat = read.csv("data/readmissions.csv").
Next, create a preparation chunk that checks the required columns and defines the analysis variables. For instance, dat$readmitted = as.integer(dat$readmitted) and dat$ward = factor(dat$ward). If treatment is recorded as 0 and 1, define its labels explicitly before fitting, so that the reference category is unambiguous. A compact check such as sum(is.na(dat[c("readmitted", "treatment", "age", "ward")])) reports how many missing entries occur in those required columns; it does not, on its own, decide how missing observations should be handled.
The document’s execution order is therefore import, check and prepare, fit, and report. Each step is visible in the source and runs again during rendering. If the report fails because the input file is absent, the relative path makes the dependency identifiable; the correction is to include the permitted input file in the project or supply an appropriate path, not to rely on a file elsewhere in the analyst’s working directory.
State the Model Before Fitting It
Before a report runs a model, its analysis section should specify the response, link, predictors, and random-effects structure. For a hypothetical binary readmission outcome, let \(Y_{ij}\) indicate readmission for patient \(i\) in ward \(j\), let \(T_{ij}\) indicate treatment, and let \(A_{ij}\) be the patient’s raw age in years. A random-intercept logistic GLMM is
This specification assumes a binary outcome, a logit mean structure as written, and normally distributed ward random intercepts with mean zero and variance \(\sigma_b^2\). It also assumes conditional independence of patient outcomes given the ward effects and included predictors. The model’s interpretation depends on these assumptions and on an appropriate data structure: for example, each patient must be assigned to the correct ward, and the outcome must be coded consistently as 0 or 1. A report cannot make a misspecified model valid merely by rendering without errors.
With the R package lme4, the corresponding fitting call is glmer(readmitted ~ treatment + I(age - 55) + (1 | ward), family = binomial(link = "logit"), data = dat). The centering expression subtracts 55 once from raw age. In this parameterization, \(\beta_1\) is the conditional log odds ratio for treatment versus its reference category, holding age and the ward random intercept fixed.
Worked Example: Fit the Model and Report a Treatment Contrast
Suppose the hypothetical dataset contains 240 patients from 12 wards, with no missing values in the four required variables after the stated preparation. Assume the model has converged and that its illustrative output gives a treatment estimate of \(0.40\) with standard error \(0.15\). These are stipulated hypothetical model results, not a claim about a real dataset.
The treatment odds ratio conditional on the ward random intercept is \(\exp(0.40)=1.492\). Thus, at the same age and ward random-intercept value, the fitted odds of readmission for the treatment category are about 1.492 times those for the reference category. This is an odds comparison, not a probability increase of 49.2 percentage points.
For illustration, a Wald 95% interval on the log-odds scale is \(0.40 \pm 1.96(0.15)\). The margin is \(1.96(0.15)=0.294\), so the endpoints are \(0.106\) and \(0.694\). Exponentiating both endpoints gives the conditional odds-ratio interval \([\exp(0.106),\exp(0.694)] \approx [1.112, 2.002]\). The interval describes values compatible with the data and model under this approximate procedure; it does not assign a 95% probability to the fixed parameter.
A report can include a concise table with the treatment coefficient, standard error, and exponentiated estimate and interval. The table heading should identify the odds-ratio scale, the comparison, and the conditional interpretation. It should not label the exponentiated estimate as a risk ratio or imply that it is a population-average effect.
Select Tables and Figures Deliberately
A report is easier to review when it displays results that answer the stated question rather than a full software dump. For example, a coefficient table might include the treatment contrast and a small number of prespecified covariates, with estimates, uncertainty, and units or contrasts. The report can create such a table from the fitted object; the essential point is to label it accurately and check that the displayed values correspond to the model actually fitted.
A useful figure for this example could show observed readmission proportions by ward. It is descriptive, not a display of the model’s fitted conditional probabilities. The analysis should identify it as such, particularly because ward sizes may differ and the plot does not adjust for age or treatment.
Worked Example: Add a Descriptive Ward Figure
Suppose a prepared hypothetical dataset includes 12 wards. For each ward, calculate the number of patients and the observed proportion readmitted. If a ward has 20 patients and 3 readmissions, its observed proportion is \(3/20=0.15\), or 15%. If another ward has 40 patients and 10 readmissions, its proportion is \(10/40=0.25\), or 25%. The two proportions have different denominators, so the report should include ward sample sizes when presenting them.
In R, summarize with a grouped operation that calculates n = dplyr::n() and proportion = mean(readmitted), then plot ward against proportion using ggplot2 and geom_point(). Add the number of patients to the displayed data or figure labels so readers can see the denominators. The figure’s caption should say “Observed readmission proportion by ward,” not “Model-adjusted ward risk.”
This figure may reveal unusual proportions or wards with little information, prompting a check of the source data and model assumptions. It does not establish that ward differences are statistically significant, and it does not replace examination of model fit. Its role in the report is to provide a transparent descriptive view alongside the fitted model.
Why the Exponentiated Interval Is Reported This Way
The interval calculation above illustrates a general reporting principle: transform both endpoints when the reported parameter is a monotone transformation of the fitted coefficient. If a Wald interval for \(\beta_1\) is \([L,U]\), and the target is the odds ratio \(\exp(\beta_1)\), then the exponential function is strictly increasing. Therefore \(L \leq \beta_1 \leq U\) implies \(\exp(L) \leq \exp(\beta_1) \leq \exp(U)\). Here, \(L=0.106\) and \(U=0.694\), giving \([1.112,2.002]\) after rounding. Transforming only the point estimate, or adding and subtracting a standard error on the odds-ratio scale, would not produce this interval.
This justification concerns the transformation of the stated Wald interval; it does not establish that the Wald approximation is accurate for every GLMM. The report should identify the interval method used and avoid presenting approximate uncertainty as exact.
Render, Inspect, and Diagnose
Render the report from its project directory using Quarto’s render command, such as quarto render analysis.qmd. Then inspect the resulting HTML, not only the source or the R console. Check that the sample description agrees with the analyzed data, the table reflects the intended model, axes and captions are legible, and no output is cut off. For a report meant to travel as a single HTML file, open the rendered file after moving or copying it away from the project directory; if figures or styling disappear, the output is not functioning as intended for that use.
Common failures have different causes and require different remedies:
- Missing objects: A chunk uses an object created only in an interactive session. Reorder chunks or create the object in the report.
- Missing files: An input path works only from one working directory. Use a project-relative path and make the required data dependency explicit.
- Missing packages: A package is installed on the author’s computer but not in the rendering environment. Document package requirements and test in the intended environment.
- Model warnings or nonconvergence: Rendering can complete while the statistical fit remains questionable. Investigate the warning and model diagnostics; do not hide it merely to produce a clean-looking report.
- Unexpectedly large output: Printing the full dataset or many diagnostics can make the HTML unwieldy. Select outputs that answer the analysis question and keep the rest available in appropriate supplementary material.
A clean render is evidence that the document’s code ran in that environment, not proof that the model assumptions hold, that the analysis is correct, or that the result will reproduce identically with different software versions. Record the package and R environment as appropriate for the project. Also review privacy: embedding resources does not anonymize data, and including patient-level rows or sensitive details in code output can disclose them.
Alternatives and Limits
A Quarto HTML report is appropriate when readers need a navigable, reproducible analysis document that can include executable R chunks and selected outputs. An R script is a simpler alternative for code-only workflows, but it does not naturally combine explanations and rendered results. A PDF may suit print-oriented review, while a dashboard may suit interactive exploration; neither is automatically a substitute for a transparent analysis report. The choice depends on readers, access needs, and the intended use of the output.
Embedding resources improves portability of the rendered document, but the original data, source file, and software environment remain separate dependencies unless explicitly provided. Large resources can make the HTML file cumbersome. Further, a report is not a guarantee of independent verification: readers may need access to permitted data and package versions to reproduce the calculations. These limits should be stated rather than obscured by the phrase “self-contained.”
Check Your Understanding
Answer each question, then compare your reasoning with the solutions below.
- Why should a report be tested in a fresh R session rather than only in an interactive session?
- In the hypothetical model, what is the interpretation of \(\exp(\beta_1)\), and what does it not represent?
- A ward has 3 readmissions among 20 patients. What observed proportion should be reported, and what denominator should accompany it?
- Does embedded-resource HTML ensure that another analyst can rerun the model? Give one reason for your answer.
Solutions. (1) A fresh session exposes objects or settings that were created interactively but are absent from the report’s execution sequence. (2) It is the treatment-versus-reference conditional odds ratio at the same age and ward random-intercept value; it is not a probability difference or automatically a population-average effect. (3) The proportion is \(3/20=0.15\), or 15%, and it should be shown with \(n=20\). (4) No. Embedded resources concern portability of the rendered HTML; rerunning still requires the source, accessible input data, and a suitable software and package environment.