Do Stroke Units Shorten Hospital Stay?

A worked random-effects meta-analysis of 9 randomised trials — and a lesson in reading heterogeneity, not hiding it

Evidence-synthesis worked example · prepared by Abdelwahab Ghazy, MD Portfolio: add your link · ORCID: add your ID · dr.abdelwahab.ghazy@gmail.com

In one breath. Pooling nine randomised trials gives a headline “stroke units save ~15 days,” but the confidence interval crosses zero and the trials disagree almost completely (I² ≈ 97%). That disagreement is not noise — the benefit scales with how sick the patients are. Where strokes are severe and stays are long, stroke units help enormously; where stays are already short, they barely move the needle. The professional deliverable here is not a single pooled number — it is the explanation for why one pooled number would mislead.

1 How to read this document

New to R? Read this once and everything below will click. Already fluent? Skip straight to the next section.

Each grey box is a code chunk — the exact instructions R runs. Under most chunks sits a purple Plain-English walkthrough that translates the code, piece by piece. A little recurring grammar unlocks all of it:

• A function is a command that does a job; you call it by its name followed by round brackets, e.g. forest(...). Whatever you put inside the brackets are its arguments — the settings you hand the command.
• <- means “save this as.” m <- metacont(...) runs metacont and keeps the result in a box named m; afterwards we just write m to reuse it.
• $ reaches inside a result to grab one piece: m$TE.random = “the pooled estimate stored inside m.”
• |> is the pipe. It takes whatever is on its left and feeds it into the function on its right, so a chain reads left-to-right instead of inside-out.
• Some sentences contain a snippet written as an r followed by a short expression, all wrapped in back-ticks. That is inline code: R replaces it with the live number when the report is built, so the words and the numbers can never disagree.
• Packages (meta, metafor, …) are toolboxes we load once at the top; each hands us ready-made functions so we never reinvent them.

2 The clinical question

We compared organised stroke-unit care — a dedicated ward with a specialist multidisciplinary team — against conventional care on a general ward, measuring length of hospital stay, in days.

PICO — P stroke inpatients · I organised stroke-unit care · C conventional general-ward care · O length of stay (days).

Lower is better: a negative Mean Difference means the stroke unit sends patients home sooner.

3 The data

Nine randomised trials assembled by Normand (1999) [1]. Each row is one trial, with the mean, SD, and sample size for both arms.

data(dat.normand1999, package = "metafor")
d <- dat.normand1999   # source = trial label; 1 = stroke unit, 2 = conventional

kable(
  d[, c("source", "n1i", "m1i", "sd1i", "n2i", "m2i", "sd2i")],
  col.names = c("Trial", "N", "Mean", "SD", "N", "Mean", "SD"),
  caption = "Length of hospital stay (days) by trial and arm."
) |>
  add_header_above(c(" " = 1, "Stroke unit" = 3, "Conventional" = 3)) |>
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = FALSE)

Length of hospital stay (days) by trial and arm.

	Stroke unit			Conventional
Trial	N	Mean	SD	N	Mean	SD
Edinburgh	155	55	47	156	75	64
Orpington-Mild	31	27	7	32	29	4
Orpington-Moderate	75	64	17	71	119	29
Orpington-Severe	18	66	20	18	137	48
Montreal-Home	8	14	8	13	18	11
Montreal-Transfer	57	19	7	52	18	4
Newcastle	34	52	45	33	41	34
Umea	110	21	16	183	31	27
Uppsala	60	30	27	52	23	20

Plain-English walkthrough.

• data(dat.normand1999, package = "metafor") loads a ready-made example dataset that ships inside the metafor toolbox.
• d <- dat.normand1999 copies it into a short name, d, that is quicker to type from here on.
• d[, c("source", "n1i", ...)] picks out just the columns we want to display, and kable(...) turns them into a tidy table.
• col.names = c(...) gives the columns human-friendly headers, and caption = ... puts a title above the table.
• The |> pipe then passes that table on to add_header_above() (which adds the “Stroke unit / Conventional” spanning row) and to kable_styling() (the striped, hover-highlight look).

How to read the columns. The suffix 1 is the intervention arm (stroke unit); 2 is the control arm (conventional care). Within each arm: n = sample size, m = mean length of stay, sd = standard deviation. (The trailing i in n1i/m1i/… is simply metafor’s naming convention and carries no analytic meaning.)

Two things are already visible and will matter later: conventional-care stay ranges from 18 to 137 days — a nearly ten-fold spread — and three of the trials come from a single centre (Orpington), deliberately split by stroke severity.

4 The analysis

We pool the trials as a Mean Difference under a random-effects model, estimating the between-study variance by REML and using the Hartung–Knapp adjustment for the confidence interval [3] — the modern default that guards against over-precise intervals when trials are few and heterogeneous.

m <- metacont(
  n1i, m1i, sd1i,               # stroke-unit arm
  n2i, m2i, sd2i,               # conventional arm
  studlab = source, data = d,
  sm               = "MD",      # effect measure: Mean Difference
  random           = TRUE,      # random-effects model
  common           = FALSE,
  method.tau       = "REML",
  method.random.ci = "HK",      # Hartung-Knapp confidence interval
  prediction       = TRUE,
  label.e = "Stroke unit", label.c = "Conventional"
)
summary(m)

##                          MD               95%-CI %W(random)
## Edinburgh          -20.0000 [-32.4744;  -7.5256]       11.0
## Orpington-Mild      -2.0000 [ -4.8271;   0.8271]       11.7
## Orpington-Moderate -55.0000 [-62.7656; -47.2344]       11.4
## Orpington-Severe   -71.0000 [-95.0223; -46.9777]        9.6
## Montreal-Home       -4.0000 [-12.1539;   4.1539]       11.4
## Montreal-Transfer    1.0000 [ -1.1176;   3.1176]       11.7
## Newcastle           11.0000 [ -8.0620;  30.0620]       10.3
## Umea               -10.0000 [-14.9237;  -5.0763]       11.6
## Uppsala              7.0000 [ -1.7306;  15.7306]       11.4
## 
## Number of studies: k = 9
## Number of observations: o = 1158 (o.e = 548, o.c = 610)
## 
##                            MD              95%-CI     t p-value
## Random effects model -15.1060 [-36.3187;  6.1067] -1.64  0.1392
## Prediction interval           [-78.8739; 48.6618]              
## 
## Quantifying heterogeneity (with 95%-CIs):
##  tau^2 = 684.6462 [292.9181; 2889.4513]; tau = 26.1657 [17.1149; 53.7536]
##  I^2 = 96.7% [95.2%; 97.7%]; H = 5.46 [4.54; 6.58]
## 
## Test of heterogeneity:
##       Q d.f.  p-value
##  238.92    8 < 0.0001
## 
## Details of meta-analysis methods:
## - Inverse variance method
## - Restricted maximum-likelihood estimator for tau^2
## - Q-Profile method for confidence interval of tau^2 and tau
## - Calculation of I^2 based on Q
## - Hartung-Knapp adjustment for random effects model (df = 8)
## - Prediction interval based on t-distribution (df = 8)

Plain-English walkthrough — metacont(). The name is short for meta-analysis of a continuous outcome. It takes each trial’s six summary numbers and pools them into one result. Reading its arguments in order:

• n1i, m1i, sd1i then n2i, m2i, sd2i — the six columns it needs: sample size, mean, and standard deviation for the stroke-unit arm (1), then the same three for the conventional arm (2).
• studlab = source — which column holds the trial name to print on each row.
• data = d — the table where all those columns live.
• sm = "MD" — the summary measure. “MD” = Mean Difference; because both arms are measured in the same unit (days), we can simply subtract them.
• random = TRUE — produce the random-effects result, which allows the true effect to genuinely differ from trial to trial.
• common = FALSE — hide the fixed-effect (“common effect”) result, which would pretend every trial is estimating one identical number.
• method.tau = "REML" — the recipe for estimating how much the trials differ (the between-study variance, τ²). REML is the modern default.
• method.random.ci = "HK" — use the Hartung–Knapp confidence interval, which stays honest (appropriately wide) when trials are few and disagree.
• prediction = TRUE — also compute the prediction interval: the range a brand-new trial might plausibly fall in.
• label.e / label.c — readable names for the experimental and control arms, used later on the plot.

Then summary(m) prints everything metacont just calculated.

5 The forest plot

forest(
  m,
  sortvar             = TE,
  leftcols            = c("studlab", "n.e", "mean.e", "sd.e", "n.c", "mean.c", "sd.c"),
  leftlabs            = c("Trial", "Total", "Mean", "SD", "Total", "Mean", "SD"),
  rightcols           = c("effect", "ci", "w.random"),
  rightlabs           = c("MD", "95% CI", "Weight"),
  label.e             = "Stroke unit", label.c = "Conventional",
  col.diamond         = "darkred",  col.diamond.lines = "black",
  col.square          = "steelblue3", col.square.lines = "steelblue3",
  common = FALSE, random = TRUE, overall = TRUE, overall.hetstat = TRUE,
  prediction          = TRUE,                # 95% prediction interval (red bar)
  spacing             = 1.2,
  colgap              = unit(3, "mm"),
  colgap.forest.left  = unit(0, "mm"),
  colgap.forest.right = unit(1, "mm"),
  xlab                = "Favours stroke unit          Favours conventional",
  squaresize = 0.8,
  digits = 1, digits.mean = 1, digits.sd = 1,
  digits.pval = 2, digits.pval.Q = 2, digits.I2 = 1, digits.tau2 = 1,
  fontsize = 11, fs.heading = 12, fs.study = 11, fs.random = 11,
  fs.hetstat = 10, fs.axis = 10, fs.xlab = 12, fs.smlab = 12,
  col.hetstat = "gray50"
)

Forest plot of the nine trials. Squares are trial-level mean differences (size <U+221D> weight); the red diamond is the random-effects pooled estimate and the red bar is the 95% prediction interval.

Plain-English walkthrough — forest(). forest(m, ...) draws the classic forest plot from the pooled object m. There are many arguments, but they only do a handful of jobs:

• What to show, and where. leftcols / leftlabs list the columns (and their headers) down the left side; rightcols / rightlabs do the same on the right (the effect, its CI, and each trial’s weight). sortvar = TE orders the trials by effect size.
• The overall result. random = TRUE, common = FALSE, overall, and overall.hetstat show the random-effects diamond and the heterogeneity statistics; prediction = TRUE adds the red prediction-interval bar.
• Colours & sizes. col.diamond, col.square, squaresize, spacing, and the colgap* arguments are pure cosmetics — the diamond colour, the box colour, how big the boxes are, and the gaps between columns.
• Rounding. Every digits* argument sets how many decimals a given number shows (e.g. digits = 1 → one decimal place for the effects).
• Font sizes. fontsize is the baseline; each fs.* fine-tunes one element (heading, study rows, heterogeneity line, axis, labels).
• xlab = "Favours stroke unit … Favours conventional" labels which side of the vertical line is the good-news side.

In short: the first few arguments decide the plot’s content, and the long tail of col.*, fs.*, and digits.* simply polishes its appearance.

6 What the numbers are telling us

res <- data.frame(
  Metric = c("Pooled MD (random effects)", "95% confidence interval",
             "95% prediction interval", "I² (inconsistency)",
             "τ² (between-trial variance)", "Cochran's Q (p-value)"),
  Value = c(
    sprintf("%.1f days", m$TE.random),
    sprintf("%.1f to %.1f days", m$lower.random, m$upper.random),
    sprintf("%.1f to %.1f days", m$lower.predict, m$upper.predict),
    sprintf("%.0f%%", m$I2 * 100),
    sprintf("%.1f", m$tau2),
    format.pval(m$pval.Q, digits = 2, eps = 0.001)
  )
)
kable(res, col.names = c("", "Result"),
      caption = "Random-effects summary.") |>
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)

Random-effects summary.

	Result
Pooled MD (random effects)	-15.1 days
95% confidence interval	-36.3 to 6.1 days
95% prediction interval	-78.9 to 48.7 days
I<U+00B2> (inconsistency)	97%
<U+03C4><U+00B2> (between-trial variance)	684.6
Cochran’s Q (p-value)	<0.001

Plain-English walkthrough. This chunk builds the summary table by hand so we control exactly what appears.

• data.frame(Metric = ..., Value = ...) makes a two-column table: the labels on the left, the formatted results on the right.
• sprintf("%.1f days", m$TE.random) is text formatting: %.1f means “a number with one decimal,” so the pooled estimate prints as, e.g., -15.1 days. Each value is pulled straight from the model with $ (m$TE.random, m$I2, m$tau2, …), so the table can never drift away from the analysis.
• format.pval(..., eps = 0.001) prints very small p-values neatly as < 0.001 instead of a long string of zeros.
• kable() |> kable_styling() renders and styles the table, just as before.

The pooled estimate of -15.1 days looks impressive, but two signals stop us from reporting it:

1. The 95% confidence interval (-36.3 to 6.1) crosses zero — the pooled effect is not statistically distinguishable from “no difference.”
2. I² is 97% [5] — the trials essentially do not agree. The 95% prediction interval [7], -78.9 to 48.7 days, makes this concrete: a new trial could plausibly land anywhere from a large benefit to a large harm.

When a prediction interval is that wide, a single pooled number is not a finding — it is an average of things that don’t belong together. So we don’t report it. We ask why the trials disagree.

7 Hunting for the cause

The trial labels hint at the answer. Three trials come from one centre (Orpington) split by stroke severity; the rest span different countries and health systems, and length of stay depends heavily on both case-mix and how each system manages discharge. Because conventional-care stay ranges from 18 to 137 days, a natural question is whether baseline stay explains the disagreement. We test it with a meta-regression.

mr <- metareg(m, ~ m2i)   # does control-arm mean explain the between-trial MD?
mr

## 
## Mixed-Effects Model (k = 9; tau^2 estimator: REML)
## 
## tau^2 (estimated amount of residual heterogeneity):     18.4346 (SE = 17.2577)
## tau (square root of estimated tau^2 value):             4.2936
## I^2 (residual heterogeneity / unaccounted variability): 72.19%
## H^2 (unaccounted variability / sampling variability):   3.60
## R^2 (amount of heterogeneity accounted for):            97.31%
## 
## Test for Residual Heterogeneity:
## QE(df = 7) = 19.1945, p-val = 0.0076
## 
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 7) = 68.4411, p-val < .0001
## 
## Model Results:
## 
##          estimate      se     tval  df    pval    ci.lb    ci.ub      
## intrcpt   12.5103  3.6203   3.4556   7  0.0106   3.9497  21.0709    * 
## m2i       -0.5505  0.0665  -8.2729   7  <.0001  -0.7079  -0.3932  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Plain-English walkthrough — metareg(). A meta-regression asks whether some trait of the trials explains why their results disagree.

• metareg(m, ~ m2i) — take the pooled object m and test one predictor. The ~ m2i part is an R formula; read it as “explain the effect using m2i” (the control-arm mean — i.e. baseline length of stay).
• Printing mr reports the slope (how much the effect shifts per extra baseline day), its p-value, and R² — the share of the between-trial disagreement that this one predictor accounts for.

Baseline stay explains almost all of the between-trial variation (R² ≈ 97%, slope p < 0.001): trials where patients would otherwise stay a long time (severe strokes) show the biggest reductions, while trials with already-short stays show almost none. The heterogeneity was never noise — it was baseline-dependence.

Read this before believing it. Regressing the Mean Difference on the control-group mean is exploratory, not confirmatory. The control mean is part of the outcome itself (MD = intervention − control), so the two are mathematically linked and part of this association is built in. Treat it as hypothesis-generating, not proof.

8 A cleaner test: the Orpington severity gradient

The caveat above has a clean answer. Instead of baseline stay — which is tangled with the outcome — we use an independent moderator: pre-randomisation stroke severity. The three Orpington trials are ideal, because they hold the hospital, team, and health system constant and vary only by severity.

orp <- subset(d, grepl("Orpington", source))
orp$severity <- c("Mild", "Moderate", "Severe")
orp$MD <- orp$m1i - orp$m2i    # stroke unit - conventional

kable(orp[, c("severity", "m1i", "m2i", "MD")],
      col.names = c("Stroke severity", "Stroke-unit mean",
                    "Conventional mean", "Mean difference (days)"),
      caption = "Orpington trials — same centre, same team, severity the only difference.") |>
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)

Orpington trials <u+2014> same centre, same team, severity the only difference.</u+2014>

	Stroke severity	Stroke-unit mean	Conventional mean	Mean difference (days)
2	Mild	27	29	-2
3	Moderate	64	119	-55
4	Severe	66	137	-71

Plain-English walkthrough.

• subset(d, grepl("Orpington", source)) — keep only the rows whose source name contains the text “Orpington.” grepl() is a text search that returns TRUE/FALSE for every row, and subset() keeps the TRUE ones.
• orp$severity <- c("Mild", "Moderate", "Severe") — add a new column by hand, labelling those three rows.
• orp$MD <- orp$m1i - orp$m2i — compute each trial’s mean difference directly (stroke-unit mean minus conventional mean) — the very subtraction metacont does internally.
• kable(...) |> kable_styling(...) — display it as a tidy table.

Holding the system fixed, the mean difference goes from -2 days (mild) → -55 days (moderate) → -71 days (severe). A clean, monotonic dose–response on an independent clinical characteristic — exactly the confirmatory signal the baseline-stay regression could only hint at.

This is the defensible version of the story: severity is clinically meaningful, available in the source trials, and not mathematically bound to the outcome. In a real review, we would pre-specify it as the moderator and extract it during data collection.

9 Sensitivity: is any single trial driving the result?

Before drawing conclusions, we check that no one trial is doing all the work, by re-pooling with each trial removed in turn (leave-one-out).

forest(metainf(m, pooled = "random"),
       col.square = "steelblue3", col.diamond = "darkred",
       fontsize = 11)

Leave-one-out analysis: the pooled estimate re-computed with each trial omitted.

Plain-English walkthrough — metainf(). metainf(m, pooled = "random") re-runs the random-effects pool nine times — each run leaving one trial out (“inf” is short for influence). Wrapping it in forest(...) draws one row per omission, so you can see at a glance whether dropping any single trial would change the conclusion.

The pooled estimate stays in a narrow band (roughly −9 to −18 days) and never crosses to favouring conventional care, so no single trial is fabricating the signal. The two omissions that pull the estimate toward zero are the Orpington-Moderate and Orpington-Severe trials — which is not fragility, but the same message again: the high-severity trials carry the benefit.

10 A note on small-study effects

With only 9 trials, a funnel plot and Egger’s test are underpowered and easily misread — guidance is not to test for small-study effects below ~10 studies, and severe heterogeneity makes an asymmetry test uninterpretable here anyway [5]. We therefore state the limitation explicitly rather than present a plot that cannot support a conclusion. Naming this is the honest, rigorous move.

11 What we learned

• The naive headline — “stroke units save 15 days” — is not supported: the CI crosses zero and I² is 97%.
• The disagreement is not noise. The benefit depends on baseline severity: largest where conventional stays are longest (severe strokes), minimal where they are already short — confirmed cleanly within the Orpington centre.
• The result is robust to leave-one-out, and small-study-effect testing is correctly withheld at this sample size.
• Transferable method habits: (1) pick the effect measure from the outcome (MD for a shared scale); (2) default to random-effects + Hartung–Knapp;
  3. read the prediction interval, not just the CI; (4) explain heterogeneity before quoting any pooled number; and (5) prefer a moderator that is independent of the outcome.

12 Real-world validation

This teaching result echoes the definitive evidence: the Cochrane Stroke Unit Trialists’ Collaboration review found that organised stroke-unit care improves outcomes, with effects that vary by patient severity and care setting [4] — reassurance that the analytic story here points the same way as the full body of trials.

13 Reproducibility & data availability

• Data: dat.normand1999 from the metafor package (Normand 1999) [1;6] — no private data; the analysis is fully reproducible from this document.
• Code: every figure and number is generated in-line from the fitted model (no hard-coded results). Use the Code ▸ Download Rmd button (top-right) to get the source.
• Environment: exact package versions are recorded in the appendix below.

14 References

1. Normand SL. Meta-analysis: formulating, evaluating, combining, and reporting. Stat Med. 1999;18(3):321–359.
2. Balduzzi S, Rücker G, Schwarzer G. How to perform a meta-analysis with R: a practical tutorial. Evid Based Ment Health. 2019;22(4):153–160. (the meta package)
3. IntHout J, Ioannidis JPA, Borm GF. The Hartung–Knapp–Sidik–Jonkman method for random-effects meta-analysis is straightforward and considerably outperforms the standard DerSimonian–Laird method. BMC Med Res Methodol. 2014;14:25.
4. Stroke Unit Trialists’ Collaboration. Organised inpatient (stroke unit) care for stroke. Cochrane Database Syst Rev. 2013;(9):CD000197.
5. Higgins JPT, Thompson SG. Quantifying heterogeneity in a meta-analysis. Stat Med. 2002;21(11):1539–1558.
6. Viechtbauer W. Conducting meta-analyses in R with the metafor package. J Stat Softw. 2010;36(3):1–48.
7. IntHout J, Ioannidis JPA, Rovers MM, Goeman JJ. Plea for routinely presenting prediction intervals in meta-analysis. BMJ Open. 2016;6:e010247.

Appendix: session information

Plain-English walkthrough. sessionInfo() prints the exact R version and the version of every package used. Including it lets anyone rerun this analysis with the same tools and get the same numbers — a hallmark of trustworthy work.

sessionInfo()

## R version 4.5.2 (2025-10-31)
## Platform: x86_64-apple-darwin20
## Running under: macOS Sequoia 15.7.7
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.5-x86_64/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.5-x86_64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
## 
## locale:
## [1] C
## 
## time zone: Africa/Cairo
## tzcode source: internal
## 
## attached base packages:
## [1] grid      stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
## [1] kableExtra_1.4.0 knitr_1.51       meta_8.3-0       metadat_1.4-0   
## [5] metabook_0.2-0  
## 
## loaded via a namespace (and not attached):
##  [1] sass_0.4.10         generics_0.1.4      xml2_1.5.2         
##  [4] stringi_1.8.7       lattice_0.22-7      hms_1.1.4          
##  [7] lme4_1.1-37         digest_0.6.39       magrittr_2.0.5     
## [10] evaluate_1.0.5      RColorBrewer_1.1-3  CompQuadForm_1.4.4 
## [13] fastmap_1.2.0       jsonlite_2.0.0      Matrix_1.7-4       
## [16] purrr_1.2.2         viridisLite_0.4.2   scales_1.4.0       
## [19] textshaping_1.0.4   numDeriv_2016.8-1.1 jquerylib_0.1.4    
## [22] reformulas_0.4.1    Rdpack_2.6.4        cli_3.6.6          
## [25] rlang_1.2.0         rbibutils_2.3       splines_4.5.2      
## [28] cachem_1.1.0        yaml_2.3.12         otel_0.2.0         
## [31] tools_4.5.2         tzdb_0.5.0          nloptr_2.2.1       
## [34] minqa_1.2.8         metafor_4.8-0       dplyr_1.2.1        
## [37] ggplot2_4.0.3       mathjaxr_1.8-0      boot_1.3-32        
## [40] vctrs_0.7.3         R6_2.6.1            lifecycle_1.0.5    
## [43] stringr_1.6.0       MASS_7.3-65         pkgconfig_2.0.3    
## [46] pillar_1.11.1       bslib_0.10.0        gtable_0.3.6       
## [49] glue_1.8.1          Rcpp_1.1.1-1.1      systemfonts_1.3.1  
## [52] xfun_0.57           tibble_3.3.1        tidyselect_1.2.1   
## [55] rstudioapi_0.18.0   farver_2.1.2        htmltools_0.5.9    
## [58] nlme_3.1-168        svglite_2.2.2       rmarkdown_2.31     
## [61] readr_2.2.0         compiler_4.5.2      S7_0.2.0