HARKing: Exploration Is Honest Until You Rewrite When the Hypothesis Was Born
On this page
อ่านฉบับภาษาไทย (Thai version)
Abstract
A subgroup result found while exploring data is a hypothesis, not a confirmation. HARKing, hypothesizing after the results are known, presents such a hypothesis as if it had been stated before the data were seen. The rewrite hides how many analyses were examined, so readers cannot judge the role of chance. With 20 independent looks at the 5% level and no true effect anywhere, the chance of at least one P value below 0.05 is 0.64. In a simulated hip-fracture registry, a declared list of twenty subgroup splits of the effect of early surgery on delirium was analysed with inverse probability weighting; five reached P below 0.05. Only frailty changes the odds ratio in the generating model; the other four carry frailty and different baseline risks, which shift the risk ratio. This article concludes that exploration is legitimate when it is labelled, that every subgroup examined should be reported, and that a post hoc finding needs new data before it is called confirmed.
A striking subgroup in a hip-fracture registry
A national hip-fracture registry holds 20,000 older adults. For each patient it records whether surgery took place within the first day after admission, called early surgery here, and whether postoperative delirium followed. The registry is fictional and its data are simulated.
The 950 patients in a small randomised trial inside the registry are set aside. An analyst studies the other 19,050, whose surgery timing was chosen by their clinicians. Frail patients and those with dementia or on anticoagulants tended to wait longer, so the analysis weights the comparison to balance such differences. Overall, the weighted risk ratio of delirium for early surgery is 0.92 (95% CI 0.86 to 0.97), which would be a modest benefit if the weighting removed all confounding.
Then an exploratory look by blood albumin, a protein often read as a marker of nutrition and general health, finds something striking. In patients whose albumin was above 37.8 g/L, the top third of recorded values, the weighted risk ratio is 0.68 (95% CI 0.57 to 0.83). In the bottom third it is 0.96.
Six months later, the draft introduction reads "We hypothesised that well-nourished patients would benefit most from early surgery". Nobody had written that hypothesis down before the data were opened.
Exploration is how hypotheses are born
Looking at data without a fixed hypothesis is a normal and valuable part of research. An exploratory analysis searches the data for patterns worth studying. A confirmatory analysis tests a hypothesis that was stated, with its analysis plan, before the data were seen.
A subgroup analysis estimates an effect separately within groups of patients defined by a characteristic measured before treatment, such as age band or frailty. When the effect, measured on a stated scale such as the risk ratio or the odds ratio, truly differs between such groups, the characteristic is an effect modifier, and the situation is called effect modification. Finding a candidate modifier by exploring a registry is a sound way to start. The trouble begins only when the paper misstates when the idea arose.
HARKing: rewriting when the hypothesis was born
HARKing stands for hypothesizing after the results are known. Kerr defined it as presenting a post hoc hypothesis, one based on or informed by the results, in a research report as if it were an a priori hypothesis, one stated in advance [1]. Rubin described three forms: building a hypothesis from the results, retrieving one from the literature afterwards, and dropping an advance hypothesis that the results did not support [2]. All three hide when and why each hypothesis entered the paper.
The draft in the opening scene shows the first form. The data led to the albumin result, and the result led to the hypothesis. The manuscript tells the story in reverse: hypothesis, then study, then result. A reader who trusts that order will treat a finding from a search as if it had survived a test.
Why it misleads: multiplicity
When early surgery has no effect, a test at the 5% level still gives P below 0.05 with probability 0.05. Multiplicity is the problem of running many such tests: the chance that at least one is falsely positive grows with their number. The probability of at least one false positive across a whole set of tests is the family-wise error rate. HARKing hides the size of that set, because the paper reports one hypothesis and one test.
Hand example: twenty looks at a world with no effect
Hand example. An analyst looks at 20 subgroups in a world where early surgery has no effect on anyone. Treat the 20 looks as independent, and count a look as a "finding" when its P value is below 0.05.
-
One look with no effect
\[ 1 - 0.05 = 0.95 \]
A single look gives P of 0.05 or more with probability 0.95.
-
Twenty looks with no finding
\[ 0.95^{20} = 0.358 \]
Independent probabilities multiply, so all 20 looks stay at 0.05 or above with probability 0.358.
-
At least one finding
\[ 1 - 0.95^{20} = 1 - 0.358 = 0.64 \]
The chance of at least one false finding is the complement, 0.64.
-
Expected number of findings
\[ 20 \times 0.05 = 1 \]
On average the analyst collects one false finding per 20 looks.
Result: With 20 independent looks and no true effect anywhere, the family-wise error rate is 0.64, or 64%, not 5%.
Real subgroup looks overlap, so they are not independent and the exact figure differs. It never falls below the 0.05 of a single look.
Two errors that often travel with HARKing
One changes the analysis, the other its wording.
p-hacking: choosing the analysis by its P value
p-hacking is trying analyses until one crosses the significance threshold and reporting only that one. Simmons and colleagues showed by simulation how a few undisclosed, ordinary choices, such as which outcome to report or when to stop recruiting, can push the false-positive rate far above the nominal 5% [3]. p-hacking bends the analysis toward a result, while HARKing bends the hypothesis to fit the result.
Causal spin: wording an association as an effect
Causal spin words an association as an effect. A weighted registry analysis finds less delirium after early surgery, and the discussion says that early surgery prevents delirium. That is a claim about what would happen under an intervention.
Such claims hold only if the assumptions behind the estimate hold, such as no unmeasured confounding. Hernán argues that the remedy is not to avoid causal words but to state the causal question openly, with the assumptions under which the estimate answers it [4].
The garden of forking paths
Multiplicity can arise even when only one analysis is run, a situation Gelman and Loken called the garden of forking paths [5]. Each choice of subgroup, cut point, outcome definition or adjustment set is a fork. When those choices follow the data, the one analysis that is run was in effect picked by the data. Its P value then overstates the evidence, even though nothing else was tried.
The registry offers many such forks: age cut at 80, 85 or 90 years or grouped in bands, and albumin split into thirds or cut at 35 or 30 g/L. Each gives a different P value.
The widget below lets you walk such a garden in a simulated world where early surgery has no effect. The panel counts the paths you have taken, written $k$, and how many reached P below 0.05. It shows that count beside $0.05 \times k$, the number expected by chance, and beside $1 - 0.95^{k}$, the chance of at least one false finding if the $k$ paths were independent.
Twenty declared splits in the simulated registry
This article analyses the same 19,050 patients again with a declared list of 20 splits, written down before this reanalysis ran but after the albumin look. The list shows the whole family that the draft left out. It includes the albumin thirds, and for that split the reanalysis is a disclosed re-look, not a confirmation, because a confirmation needs data the hypothesis did not come from. Unlike the hand example, early surgery has a real effect in the registry, and one characteristic, frailty, changes its odds ratio.
ASA grade 3 or higher, used in several splits, marks severe systemic disease on the anaesthetic grading of physical status. Each level gets its own weighted risk ratio of delirium for early versus later surgery. This is an average treatment effect (ATE): the risk if all its patients had early surgery divided by the risk if all had later surgery. Each split gets one P value for a difference in risk ratio between its levels, from a test that compares their log risk ratios.
The weighting is inverse probability of treatment weighting (IPTW): each patient counts as the inverse of their estimated probability of the surgery timing they had. The weighted groups are then balanced on the characteristics in the propensity model, the model that estimates those probabilities. A separate guide to IPTW explains how the weights are built and checked.
All 20 declared splits
| Split | Risk ratio (95% CI) by level (albumin in g/L) | P for a difference in risk ratio between levels |
|---|---|---|
| Sex | men: 1.00 (0.89 to 1.12); women: 0.88 (0.83 to 0.94) | 0.0540 |
| Age, cut at 80 years | under 80: 0.88 (0.80 to 0.97); 80 or over: 0.92 (0.87 to 0.98) | 0.4035 |
| Age, cut at 85 years | under 85: 0.89 (0.83 to 0.95); 85 or over: 0.93 (0.87 to 0.99) | 0.4511 |
| Age, cut at 90 years | under 90: 0.90 (0.85 to 0.95); 90 or over: 0.94 (0.88 to 1.02) | 0.2930 |
| Age, three bands | under 75: 0.88 (0.77 to 1.01); 75 to 84: 0.90 (0.84 to 0.97); 85 or over: 0.93 (0.87 to 0.99) | 0.7982 |
| Age, by decade | under 70: 0.96 (0.78 to 1.18); 70 to 79: 0.86 (0.78 to 0.95); 80 to 89: 0.91 (0.85 to 0.98); 90 or over: 0.94 (0.88 to 1.02) | 0.5039 |
| ASA grade | 1 or 2: 0.90 (0.77 to 1.07); 3 or higher: 0.91 (0.87 to 0.97) | 0.8911 |
| Anticoagulant | no: 0.93 (0.88 to 0.98); yes: 0.95 (0.76 to 1.19) | 0.8310 |
| Dementia | no: 0.86 (0.80 to 0.93); yes: 0.96 (0.92 to 1.01) | 0.0160 |
| Frailty | not frail: 0.65 (0.58 to 0.73); frail: 1.02 (0.97 to 1.07) | below 0.001 |
| Albumin, thirds | 34.0 or less: 0.96 (0.90 to 1.04); above 34.0 to 37.8: 0.90 (0.71 to 1.13); above 37.8: 0.68 (0.57 to 0.83) | 0.0043 |
| Albumin, cut at 35 g/L | under 35: 0.97 (0.89 to 1.06); 35 or more: 0.74 (0.65 to 0.84) | 0.0006 |
| Albumin, cut at 30 g/L | under 30: 0.94 (0.84 to 1.06); 30 or more: 0.89 (0.80 to 0.98) | 0.4366 |
| Albumin recorded | missing: 0.94 (0.88 to 1.01); recorded: 0.90 (0.82 to 0.98) | 0.3644 |
| Sex by age 80 | men under 80: 0.95 (0.81 to 1.11); men 80 or over: 1.03 (0.91 to 1.17); women under 80: 0.86 (0.77 to 0.96); women 80 or over: 0.88 (0.83 to 0.94) | 0.1131 |
| Sex by ASA grade | men ASA 1 or 2: 1.07 (0.76 to 1.50); men ASA 3 or higher: 0.99 (0.90 to 1.08); women ASA 1 or 2: 0.82 (0.73 to 0.91); women ASA 3 or higher: 0.89 (0.83 to 0.95) | 0.0533 |
| Sex by dementia | men without dementia: 0.88 (0.78 to 1.00); men with dementia: 1.08 (1.01 to 1.15); women without dementia: 0.85 (0.78 to 0.93); women with dementia: 0.91 (0.86 to 0.96) | below 0.001 |
| Dementia by ASA grade | neither: 0.83 (0.71 to 0.97); ASA 3 or higher only: 0.87 (0.80 to 0.95); dementia only: 0.94 (0.81 to 1.08); both: 0.97 (0.93 to 1.02) | 0.0581 |
| Anticoagulant by ASA grade | neither: 0.87 (0.79 to 0.95); ASA 3 or higher only: 0.94 (0.89 to 1.00); anticoagulant only: 1.65 (0.89 to 3.07); both: 0.86 (0.75 to 0.99) | 0.0985 |
| Age 80 by dementia | under 80 without dementia: 0.82 (0.71 to 0.94); under 80 with dementia: 0.95 (0.85 to 1.05); 80 or over without dementia: 0.88 (0.81 to 0.96); 80 or over with dementia: 0.96 (0.92 to 1.01) | 0.0654 |
What the generating model says
Because the data are simulated, the truth is known. In the model that generated delirium, early surgery multiplies the odds of delirium by 0.50 in patients who are not frail and by 1.00 in frail patients. No other characteristic enters the surgery term of that model, so frailty is the only modifier of the odds ratio for early surgery.
Across the whole simulated population behind the registry, trial patients included, early surgery lowers the true risk of delirium from 0.157 to 0.092 in patients who are not frail, a risk ratio of 0.59. In frail patients the true risk is 0.561 with or without early surgery, a risk ratio of 1.00. The frailty split, estimated in the patients outside the trial, recovers this pattern with P below 0.001. Its risk ratio is 0.65 (95% CI 0.58 to 0.73) in patients who are not frail and 1.02 (0.97 to 1.07) in frail patients, intervals that include both true values.
Why four other splits also reached P below 0.05
Four more splits reached P below 0.05: dementia, albumin in thirds, albumin cut at 35 g/L, and sex by dementia. They are not chance findings, yet in the generating model none of dementia, albumin or sex changes the odds ratio for early surgery. Frailty makes dementia more likely and lowers albumin, so these splits sort frail patients unevenly between their levels. Levels with more frail patients, such as dementia or low albumin, show risk ratios nearer 1, because more of their patients gain nothing from early surgery.
The risk-ratio scale adds to these gaps. The generating model fixes an odds ratio, and a fixed odds ratio gives a risk ratio close to itself when delirium is uncommon and closer to 1 when delirium is common. Dementia raises the risk of delirium without early surgery and higher albumin lowers it, so their levels truly differ in risk ratio even among patients who are not frail, which is effect modification on the risk-ratio scale. Whether an effect differs between groups depends on the scale it is measured on, a question taken further in a separate article on interaction and scale.
Read as a story, the albumin split invites the nutrition hypothesis that the draft wrote into its introduction. In this simulated world that story is wrong. Albumin does not change the odds ratio for early surgery in any patient; its split stands in for frailty and for differences in baseline risk. The subgroup was real as a pattern and false as an explanation.
Near misses and striking single levels
Two splits only just missed the threshold: sex, with P 0.0540, and sex by ASA grade, with P 0.0533. Sex plays no part in how delirium was generated and is unrelated to frailty and to every other cause of delirium in this world, so men and women share the same true risk ratio. The gap between men (1.00) and women (0.88) is chance alone. A slightly different fork could move such a split across 0.05, and a HARKed paper would then explain it.
Men with dementia have a risk ratio of 1.08 (95% CI 1.01 to 1.15), an interval wholly above 1. Yet early surgery cannot raise anyone's risk of delirium in this world, because its true odds ratio is 0.50 or 1.00 for every patient. With dozens of level intervals in the table, a few are expected to miss the truth by chance alone, and this is one of them.
Stata: one loop over the 20 declared splits
local splits sex age_80 age_85 age_90 age_band3 age_decade asa3 anticoag dementia frail
local splits `splits' albumin_tertile albumin_35 albumin_30 albumin_recorded
local splits `splits' sex_age80 sex_asa3 sex_dementia dementia_asa3 anticoag_asa3 age80_dementia
local lv_sex male female
local lv_age_80 lt80 ge80
local lv_age_85 lt85 ge85
local lv_age_90 lt90 ge90
local lv_age_band3 lt75 75to84 ge85
local lv_age_decade lt70 70to79 80to89 ge90
local lv_asa3 no yes
local lv_anticoag no yes
local lv_dementia no yes
local lv_frail no yes
local lv_albumin_tertile t1 t2 t3
local lv_albumin_35 lt35 ge35
local lv_albumin_30 lt30 ge30
local lv_albumin_recorded no yes
local lv_sex_age80 male_lt80 male_ge80 female_lt80 female_ge80
local lv_sex_asa3 male_asa12 male_asa3 female_asa12 female_asa3
local lv_sex_dementia male_nodem male_dem female_nodem female_dem
local lv_dementia_asa3 neither asa3_only dementia_only both
local lv_anticoag_asa3 neither asa3_only anticoag_only both
local lv_age80_dementia lt80_nodem lt80_dem ge80_nodem ge80_dem
local drop_sex female
local drop_asa3 asa3
local drop_anticoag anticoag
local drop_dementia dementia fd
local drop_frail frail fd
local drop_sex_age80 female
local drop_sex_asa3 female asa3
local drop_sex_dementia female dementia fd
local drop_dementia_asa3 dementia fd asa3
local drop_anticoag_asa3 anticoag asa3
local drop_age80_dementia dementia fd
* ---- IPTW risk ratio of delirium within each level, and the P value for a difference between levels ----
* within a split the levels are independent samples: with b the log risk ratio of a level and w = 1/SE^2,
* the Wald chi-squared statistic sum w*b^2 - (sum w*b)^2/sum w has k - 1 df (for two levels it is the
* square of the z statistic (b1 - b2)/sqrt(SE1^2 + SE2^2))
local ps_all age_c age_c2 female frail dementia fd asa3 anticoag
scalar sc_sig = 0
local i = 0
display as text _newline %-18s "split" %-15s "level" %7s "n" %8s "RR" " 95% CI"
foreach s of local splits {
local i = `i' + 1
local terms : list ps_all - drop_`s'
local k : word count `lv_`s''
scalar sc_sw = 0
scalar sc_swb = 0
scalar sc_swb2 = 0
forvalues j = 1/`k' {
local lv : word `j' of `lv_`s''
quietly teffects ipw (delirium) (surg24 `terms', logit) if g_`s' == `j', pomeans
pomrr
matrix RES = nullmat(RES) \ (`i', `j', e(N), exp(scalar(sc_b)), scalar(sc_lb), scalar(sc_ub))
scalar sc_w = 1/scalar(sc_se)^2
scalar sc_sw = scalar(sc_sw) + scalar(sc_w)
scalar sc_swb = scalar(sc_swb) + scalar(sc_w)*scalar(sc_b)
scalar sc_swb2 = scalar(sc_swb2) + scalar(sc_w)*scalar(sc_b)^2
display as text %-18s cond(`j' == 1, "`s'", "") %-15s "`lv'" %7.0f e(N) %8.2f exp(scalar(sc_b)) ///
" " %4.2f scalar(sc_lb) " to " %4.2f scalar(sc_ub)
}
scalar sc_p = chi2tail(`k' - 1, scalar(sc_swb2) - scalar(sc_swb)^2/scalar(sc_sw))
matrix PV = nullmat(PV) \ scalar(sc_p)
scalar sc_sig = scalar(sc_sig) + (scalar(sc_p) < 0.05)
display as text %-18s "" "P for a difference between levels = " %6.4f scalar(sc_p) ///
cond(scalar(sc_p) < 0.05, " (P < 0.05)", "")
}
display as text _newline "splits with P < 0.05: " scalar(sc_sig) " of " `i'. display as text _newline %-18s "split" %-15s "level" %7s "n" %8s "RR" " 95% CI"
split level n RR 95% CI
. foreach s of local splits {
2. local i = `i' + 1
3. local terms : list ps_all - drop_`s'
4. local k : word count `lv_`s''
5. scalar sc_sw = 0
6. scalar sc_swb = 0
7. scalar sc_swb2 = 0
8. forvalues j = 1/`k' {
9. local lv : word `j' of `lv_`s''
10. quietly teffects ipw (delirium) (surg24 `terms', logit) if g_`s' == `j', pomeans
11. pomrr
12. matrix RES = nullmat(RES) \ (`i', `j', e(N), exp(scalar(sc_b)), scalar(sc_lb), scalar(sc_ub))
13. scalar sc_w = 1/scalar(sc_se)^2
14. scalar sc_sw = scalar(sc_sw) + scalar(sc_w)
15. scalar sc_swb = scalar(sc_swb) + scalar(sc_w)*scalar(sc_b)
16. scalar sc_swb2 = scalar(sc_swb2) + scalar(sc_w)*scalar(sc_b)^2
17. display as text %-18s cond(`j' == 1, "`s'", "") %-15s "`lv'" %7.0f e(N) %8.2f exp(scalar(sc_b)) ///
> " " %4.2f scalar(sc_lb) " to " %4.2f scalar(sc_ub)
18. }
19. scalar sc_p = chi2tail(`k' - 1, scalar(sc_swb2) - scalar(sc_swb)^2/scalar(sc_sw))
20. matrix PV = nullmat(PV) \ scalar(sc_p)
21. scalar sc_sig = scalar(sc_sig) + (scalar(sc_p) < 0.05)
22. display as text %-18s "" "P for a difference between levels = " %6.4f scalar(sc_p) ///
> cond(scalar(sc_p) < 0.05, " (P < 0.05)", "")
23. }
sex male 5665 1.00 0.89 to 1.12
female 13385 0.88 0.83 to 0.94
P for a difference between levels = 0.0540
age_80 lt80 8824 0.88 0.80 to 0.97
ge80 10226 0.92 0.87 to 0.98
P for a difference between levels = 0.4035
age_85 lt85 12075 0.89 0.83 to 0.95
ge85 6975 0.93 0.87 to 0.99
P for a difference between levels = 0.4511
age_90 lt90 14780 0.90 0.85 to 0.95
ge90 4270 0.94 0.88 to 1.02
P for a difference between levels = 0.2930
age_band3 lt75 5787 0.88 0.77 to 1.01
75to84 6288 0.90 0.84 to 0.97
ge85 6975 0.93 0.87 to 0.99
P for a difference between levels = 0.7982
age_decade lt70 3325 0.96 0.78 to 1.18
70to79 5499 0.86 0.78 to 0.95
80to89 5956 0.91 0.85 to 0.98
ge90 4270 0.94 0.88 to 1.02
P for a difference between levels = 0.5039
asa3 no 7674 0.90 0.77 to 1.07
yes 11376 0.91 0.87 to 0.97
P for a difference between levels = 0.8911
anticoag no 15798 0.93 0.88 to 0.98
yes 3252 0.95 0.76 to 1.19
P for a difference between levels = 0.8310
dementia no 14024 0.86 0.80 to 0.93
yes 5026 0.96 0.92 to 1.01
P for a difference between levels = 0.0160 (P < 0.05)
frail no 10899 0.65 0.58 to 0.73
yes 8151 1.02 0.97 to 1.07
P for a difference between levels = 0.0000 (P < 0.05)
albumin_tertile t1 4471 0.96 0.90 to 1.04
t2 4409 0.90 0.71 to 1.13
t3 4270 0.68 0.57 to 0.83
P for a difference between levels = 0.0043 (P < 0.05)
albumin_35 lt35 5494 0.97 0.89 to 1.06
ge35 7656 0.74 0.65 to 0.84
P for a difference between levels = 0.0006 (P < 0.05)
albumin_30 lt30 1242 0.94 0.84 to 1.06
ge30 11908 0.89 0.80 to 0.98
P for a difference between levels = 0.4366
albumin_recorded no 5900 0.94 0.88 to 1.01
yes 13150 0.90 0.82 to 0.98
P for a difference between levels = 0.3644
sex_age80 male_lt80 2647 0.95 0.81 to 1.11
male_ge80 3018 1.03 0.91 to 1.17
female_lt80 6177 0.86 0.77 to 0.96
female_ge80 7208 0.88 0.83 to 0.94
P for a difference between levels = 0.1131
sex_asa3 male_asa12 2223 1.07 0.76 to 1.50
male_asa3 3442 0.99 0.90 to 1.08
female_asa12 5451 0.82 0.73 to 0.91
female_asa3 7934 0.89 0.83 to 0.95
P for a difference between levels = 0.0533
sex_dementia male_nodem 4114 0.88 0.78 to 1.00
male_dem 1551 1.08 1.01 to 1.15
female_nodem 9910 0.85 0.78 to 0.93
female_dem 3475 0.91 0.86 to 0.96
P for a difference between levels = 0.0000 (P < 0.05)
dementia_asa3 neither 6024 0.83 0.71 to 0.97
asa3_only 8000 0.87 0.80 to 0.95
dementia_only 1650 0.94 0.81 to 1.08
both 3376 0.97 0.93 to 1.02
P for a difference between levels = 0.0581
anticoag_asa3 neither 6435 0.87 0.79 to 0.95
asa3_only 9363 0.94 0.89 to 1.00
anticoag_only 1239 1.65 0.89 to 3.07
both 2013 0.86 0.75 to 0.99
P for a difference between levels = 0.0985
age80_dementia lt80_nodem 7411 0.82 0.71 to 0.94
lt80_dem 1413 0.95 0.85 to 1.05
ge80_nodem 6613 0.88 0.81 to 0.96
ge80_dem 3613 0.96 0.92 to 1.01
P for a difference between levels = 0.0654
. display as text _newline "splits with P < 0.05: " scalar(sc_sig) " of " `i'
splits with P < 0.05: 5 of 20R: the same loop with WeightIt
splits <- list(
sex = list(lv = c("male", "female"), drop = "female"),
age_80 = list(lv = c("lt80", "ge80"), drop = character(0)),
age_85 = list(lv = c("lt85", "ge85"), drop = character(0)),
age_90 = list(lv = c("lt90", "ge90"), drop = character(0)),
age_band3 = list(lv = c("lt75", "75to84", "ge85"), drop = character(0)),
age_decade = list(lv = c("lt70", "70to79", "80to89", "ge90"), drop = character(0)),
asa3 = list(lv = c("no", "yes"), drop = "asa3"),
anticoag = list(lv = c("no", "yes"), drop = "anticoag"),
dementia = list(lv = c("no", "yes"), drop = c("dementia", "fd")),
frail = list(lv = c("no", "yes"), drop = c("frail", "fd")),
albumin_tertile = list(lv = c("t1", "t2", "t3"), drop = character(0)),
albumin_35 = list(lv = c("lt35", "ge35"), drop = character(0)),
albumin_30 = list(lv = c("lt30", "ge30"), drop = character(0)),
albumin_recorded = list(lv = c("no", "yes"), drop = character(0)),
sex_age80 = list(lv = c("male_lt80", "male_ge80", "female_lt80", "female_ge80"), drop = "female"),
sex_asa3 = list(lv = c("male_asa12", "male_asa3", "female_asa12", "female_asa3"), drop = c("female", "asa3")),
sex_dementia = list(lv = c("male_nodem", "male_dem", "female_nodem", "female_dem"),
drop = c("female", "dementia", "fd")),
dementia_asa3 = list(lv = c("neither", "asa3_only", "dementia_only", "both"), drop = c("dementia", "fd", "asa3")),
anticoag_asa3 = list(lv = c("neither", "asa3_only", "anticoag_only", "both"), drop = c("anticoag", "asa3")),
age80_dementia = list(lv = c("lt80_nodem", "lt80_dem", "ge80_nodem", "ge80_dem"), drop = c("dementia", "fd"))
)
# ---- IPTW risk ratio of delirium within each level, and the P value for a difference between levels ----
# within a split the levels are independent samples: with b the log risk ratio of a level and w = 1/SE^2,
# the Wald chi-squared statistic sum w*b^2 - (sum w*b)^2/sum w has k - 1 df (for two levels it is the
# square of the z statistic (b1 - b2)/sqrt(SE1^2 + SE2^2))
ps_all <- c("age_c", "age_c2", "female", "frail", "dementia", "fd", "asa3", "anticoag")
cat(sprintf("\n%-18s%-15s%7s%8s 95%% CI\n", "split", "level", "n", "RR"))
res <- lapply(names(splits), function(s) {
sp <- splits[[s]]
f_ps <- reformulate(setdiff(ps_all, sp$drop), response = "surg24")
rows <- lapply(seq_along(sp$lv), function(j) {
dj <- d0[!is.na(d0[[paste0("g_", s)]]) & d0[[paste0("g_", s)]] == j, ]
W <- weightit(f_ps, data = dj, method = "glm", estimand = "ATE")
fit <- glm_weightit(delirium ~ surg24, data = dj, weightit = W, family = quasipoisson(link = "log"))
b <- coef(fit)[["surg24"]]; se <- sqrt(vcov(fit)["surg24", "surg24"])
cat(sprintf("%-18s%-15s%7d%8.2f %4.2f to %4.2f\n", if (j == 1) s else "", sp$lv[j], nrow(dj),
exp(b), exp(b - z * se), exp(b + z * se)))
data.frame(split = s, level = sp$lv[j], n = nrow(dj), rr = exp(b), lb = exp(b - z * se),
ub = exp(b + z * se), b = b, se = se)
})
r <- do.call(rbind, rows)
w <- 1 / r$se^2
p <- pchisq(sum(w * r$b^2) - sum(w * r$b)^2 / sum(w), df = nrow(r) - 1, lower.tail = FALSE)
cat(sprintf("%-18sP for a difference between levels = %6.4f%s\n", "", p,
if (p < 0.05) " (P < 0.05)" else ""))
list(levels = r, p = p)
})
names(res) <- names(splits)
sig <- sum(sapply(res, function(x) x$p < 0.05))
cat(sprintf("\nsplits with P < 0.05: %d of %d\n", sig, length(res))) > cat(sprintf("\n%-18s%-15s%7s%8s 95%% CI\n", "split",
+ "level", "n", "RR"))
split level n RR 95% CI
> res <- lapply(names(splits), function(s) {
+ sp <- splits[[s]]
+ f_ps <- reformulate(setdiff(ps_all, sp$drop), response = "surg24")
+ rows <- lapply(seq_along(sp$lv), function(j) {
+ dj <- d0[!is.na(d0[[paste0("g_", s)]]) & d0[[paste0("g_",
+ s)]] == j, ]
+ W <- weightit(f_ps, data = dj, method = "glm", estimand = "ATE")
+ fit <- glm_weightit(delirium ~ surg24, data = dj, weightit = W,
+ family = quasipoisson(link = "log"))
+ b <- coef(fit)[["surg24"]]
+ se <- sqrt(vcov(fit)["surg24", "surg24"])
+ cat(sprintf("%-18s%-15s%7d%8.2f %4.2f to %4.2f\n",
+ if (j == 1)
+ s
+ else "", sp$lv[j], nrow(dj), exp(b), exp(b - z *
+ se), exp(b + z * se)))
+ data.frame(split = s, level = sp$lv[j], n = nrow(dj),
+ rr = exp(b), lb = exp(b - z * se), ub = exp(b + z *
+ se), b = b, se = se)
+ })
+ r <- do.call(rbind, rows)
+ w <- 1/r$se^2
+ p <- pchisq(sum(w * r$b^2) - sum(w * r$b)^2/sum(w), df = nrow(r) -
+ 1, lower.tail = FALSE)
+ cat(sprintf("%-18sP for a difference between levels = %6.4f%s\n",
+ "", p, if (p < 0.05)
+ " (P < 0.05)"
+ else ""))
+ list(levels = r, p = p)
+ })
sex male 5665 1.00 0.89 to 1.12
female 13385 0.88 0.83 to 0.94
P for a difference between levels = 0.0540
age_80 lt80 8824 0.88 0.80 to 0.97
ge80 10226 0.92 0.87 to 0.98
P for a difference between levels = 0.4035
age_85 lt85 12075 0.89 0.83 to 0.95
ge85 6975 0.93 0.87 to 0.99
P for a difference between levels = 0.4511
age_90 lt90 14780 0.90 0.85 to 0.95
ge90 4270 0.94 0.88 to 1.02
P for a difference between levels = 0.2930
age_band3 lt75 5787 0.88 0.77 to 1.01
75to84 6288 0.90 0.84 to 0.97
ge85 6975 0.93 0.87 to 0.99
P for a difference between levels = 0.7982
age_decade lt70 3325 0.96 0.78 to 1.18
70to79 5499 0.86 0.78 to 0.95
80to89 5956 0.91 0.85 to 0.98
ge90 4270 0.94 0.88 to 1.02
P for a difference between levels = 0.5039
asa3 no 7674 0.90 0.77 to 1.07
yes 11376 0.91 0.87 to 0.97
P for a difference between levels = 0.8911
anticoag no 15798 0.93 0.88 to 0.98
yes 3252 0.95 0.76 to 1.19
P for a difference between levels = 0.8310
dementia no 14024 0.86 0.80 to 0.93
yes 5026 0.96 0.92 to 1.01
P for a difference between levels = 0.0160 (P < 0.05)
frail no 10899 0.65 0.58 to 0.73
yes 8151 1.02 0.97 to 1.07
P for a difference between levels = 0.0000 (P < 0.05)
albumin_tertile t1 4471 0.96 0.90 to 1.04
t2 4409 0.90 0.71 to 1.13
t3 4270 0.68 0.57 to 0.83
P for a difference between levels = 0.0043 (P < 0.05)
albumin_35 lt35 5494 0.97 0.89 to 1.06
ge35 7656 0.74 0.65 to 0.84
P for a difference between levels = 0.0006 (P < 0.05)
albumin_30 lt30 1242 0.94 0.84 to 1.06
ge30 11908 0.89 0.80 to 0.98
P for a difference between levels = 0.4366
albumin_recorded no 5900 0.94 0.88 to 1.01
yes 13150 0.90 0.82 to 0.98
P for a difference between levels = 0.3644
sex_age80 male_lt80 2647 0.95 0.81 to 1.11
male_ge80 3018 1.03 0.91 to 1.17
female_lt80 6177 0.86 0.77 to 0.96
female_ge80 7208 0.88 0.83 to 0.94
P for a difference between levels = 0.1131
sex_asa3 male_asa12 2223 1.07 0.76 to 1.50
male_asa3 3442 0.99 0.90 to 1.08
female_asa12 5451 0.82 0.73 to 0.91
female_asa3 7934 0.89 0.83 to 0.95
P for a difference between levels = 0.0533
sex_dementia male_nodem 4114 0.88 0.78 to 1.00
male_dem 1551 1.08 1.01 to 1.15
female_nodem 9910 0.85 0.78 to 0.93
female_dem 3475 0.91 0.86 to 0.96
P for a difference between levels = 0.0000 (P < 0.05)
dementia_asa3 neither 6024 0.83 0.71 to 0.97
asa3_only 8000 0.87 0.80 to 0.95
dementia_only 1650 0.94 0.81 to 1.08
both 3376 0.97 0.93 to 1.02
P for a difference between levels = 0.0581
anticoag_asa3 neither 6435 0.87 0.79 to 0.95
asa3_only 9363 0.94 0.89 to 1.00
Warning in weightit(f_ps, data = dj, method = "glm", estimand = "ATE") :
Some extreme weights were generated. Examine them with `summary()` and maybe
trim them with `trim()`.
anticoag_only 1239 1.65 0.89 to 3.07
both 2013 0.86 0.75 to 0.99
P for a difference between levels = 0.0985
age80_dementia lt80_nodem 7411 0.82 0.71 to 0.94
lt80_dem 1413 0.95 0.85 to 1.05
ge80_nodem 6613 0.88 0.81 to 0.96
ge80_dem 3613 0.96 0.92 to 1.01
P for a difference between levels = 0.0654
> names(res) <- names(splits)
> sig <- sum(sapply(res, function(x) x$p < 0.05))
> cat(sprintf("\nsplits with P < 0.05: %d of %d\n",
+ sig, length(res)))
splits with P < 0.05: 5 of 20 What was done versus what may be claimed
The cure for HARKing is not silence about exploratory results but wording each result at the level its design supports.
Honest wording by what was done
| What was done | What it supports | Honest wording |
|---|---|---|
| Looked through the data with no prior hypothesis and noticed a pattern | A description of the data | "In an exploratory look, the weighted risk ratio was lower in patients with higher albumin." |
| Ran an analysis chosen after seeing the data | A new hypothesis, to be tested in other data | "A post hoc, exploratory analysis suggested that ..." |
| Tested a hypothesis written down, with its analysis plan, before the data were seen | A confirmatory test of that hypothesis | "We tested the prespecified hypothesis that ..." |
| Analysed every subgroup on a declared list | A family of tests, read together | "Of 20 subgroup analyses on a declared list, 5 had P below 0.05; all 20 are reported." |
| Estimated an effect from observational data under stated assumptions | An effect estimate, if the assumptions hold | "Assuming no unmeasured confounding, early surgery was estimated to ..." |
| Randomised the intervention | The most direct evidence about an effect | "Early surgery reduced delirium in this trial." |
Safeguards that keep exploration honest
Preregistration is a time-stamped record of hypotheses and analysis plans, made before the data are seen [6]. It does not forbid exploration. It makes the boundary between planned and unplanned analysis visible, and changes to the plan are reported rather than hidden.
Trial registries and published protocols let readers compare what was planned with what was reported. One such comparison found outcomes reported selectively and primary outcomes changed [7]. An observational analysis can keep the same discipline with an analysis plan dated before the outcome data are opened.
Credibility criteria for subgroup claims give a checklist [8]. Was the subgroup one of a few specified in advance, with its expected direction? Does a test of interaction, a test of whether the effect differs between levels, make chance an unlikely explanation? Is the effect independent of other subgroup effects, and consistent across studies and related outcomes?
The dementia and albumin splits in the registry fail the independence question, because each carries frailty. Reasoning from the causal question first also helps [4]: name the effect of interest and the characteristics that might modify it, with a reason, before any subgroup is examined. A post hoc finding that survives all of this is still a hypothesis until new data confirm it.
Common misreadings and their fixes
-
"All post hoc analysis is worthless."
The fault in the opening scene lies in the rewritten introduction, not in the exploration that found the albumin pattern.
Fix: Post hoc analysis is how new hypotheses are found. The error is presenting it as if it had been planned; label it, report everything examined, and test it in new data.
-
"The split with the smallest P value shows where the effect differs."
Five splits reached P below 0.05, yet only frailty changes the odds ratio for early surgery. The dementia and albumin splits carry frailty and different baseline risks, which shift the risk ratio.
Fix: Ask what else differs between the levels, check whether the scale of the effect measure could create the gap, and prefer a modifier named in advance with a reason.
-
"P just above 0.05 means no difference, and P just below means a real one."
The sex split, with P 0.0540, reflects chance alone, while the dementia split, with P 0.0160, reflects frailty and a higher baseline risk. The ASA split, with P 0.8911, shows no gap even though higher ASA grades carry a higher baseline risk of delirium, so the 0.05 threshold neither marks the splits where frailty or baseline risk shifts the risk ratio nor clears the ones where only chance is at work. In this world sex has the same true risk ratio in men and women, while dementia and albumin truly differ in risk ratio because each carries frailty.
Fix: Report the estimate and its interval for every split, and let the design and the mechanism decide what deserves a follow-up study.
-
"A multiplicity correction fixes HARKing."
A Bonferroni correction, which compares each P value with 0.05 divided by the number of tests, controls the family-wise error rate only over the tests that are reported.
Fix: Disclose every analysis examined, because a correction works only on a visible family.
What to do in your own analysis
- Before opening the outcome data, write down the hypotheses, each subgroup with its expected direction and the test of interaction, and date or register the plan.
- Label every analysis in the paper as prespecified or post hoc, and report every subgroup examined, including those that showed nothing.
- Judge a subgroup by a test of interaction and by its estimates with their intervals, not by separate P values within levels.
- Consider modelling a continuous candidate modifier, such as age or albumin, as continuous rather than choosing a cut point after seeing the data.
Glossary
- HARKing (การตั้งสมมติฐานหลังรู้ผล)
- Hypothesizing after the results are known: presenting a hypothesis formed after seeing the results as if it had been stated in advance.
- exploratory analysis
- An analysis that searches the data for patterns worth studying; it generates hypotheses.
- confirmatory analysis
- An analysis that tests a hypothesis stated, with its analysis plan, before the data were seen.
- subgroup analysis
- An estimate of an effect within groups of patients defined by a characteristic measured before treatment.
- effect modification
- A true difference in the effect of a treatment between groups of patients, on a stated scale.
- multiplicity
- The growth in the chance of at least one false positive as more tests are run.
- family-wise error
- The probability of at least one false positive across a whole set of tests.
- p-hacking (การปั่นค่า p)
- Trying analyses until one crosses a significance threshold and reporting only that one.
- garden of forking paths (ทางแยกของการวิเคราะห์)
- Data-dependent analysis choices that inflate false positives even when only one analysis is run.
- preregistration (การลงทะเบียนแผนการวิเคราะห์ล่วงหน้า)
- A time-stamped record of hypotheses and analysis plans made before the data are seen.
- average treatment effect (ATE) (ผลเฉลี่ยของการรักษาในประชากรทั้งหมด (ATE))
- The average effect of treatment had everyone in a group been treated versus had no one been treated, on a stated scale such as the risk ratio.
- inverse probability of treatment weighting (IPTW) (การถ่วงน้ำหนักด้วยส่วนกลับของความน่าจะเป็นในการได้รับการรักษา)
- Weighting each patient by the inverse probability of the treatment they received, so the weighted groups are balanced on the characteristics in the model for that probability.
References
- Kerr NL. HARKing: hypothesizing after the results are known. Pers Soc Psychol Rev. 1998;2(3):196-217. doi:10.1207/s15327957pspr0203_4 https://doi.org/10.1207/s15327957pspr0203_4
- Rubin M. When does HARKing hurt? Identifying when different types of undisclosed post hoc hypothesizing harm scientific progress. Rev Gen Psychol. 2017;21(4):308-320. doi:10.1037/gpr0000128 https://doi.org/10.1037/gpr0000128
- Simmons JP, Nelson LD, Simonsohn U. False-positive psychology: undisclosed flexibility in data collection and analysis allows presenting anything as significant. Psychol Sci. 2011;22(11):1359-1366. doi:10.1177/0956797611417632 https://doi.org/10.1177/0956797611417632
- Hernán MA. The C-word: scientific euphemisms do not improve causal inference from observational data. Am J Public Health. 2018;108(5):616-619. doi:10.2105/AJPH.2018.304337 https://doi.org/10.2105/AJPH.2018.304337
- Gelman A, Loken E. The statistical crisis in science. Am Sci. 2014;102(6):460-465. doi:10.1511/2014.111.460 https://doi.org/10.1511/2014.111.460
- Nosek BA, Ebersole CR, DeHaven AC, Mellor DT. The preregistration revolution. Proc Natl Acad Sci U S A. 2018;115(11):2600-2606. doi:10.1073/pnas.1708274114 https://doi.org/10.1073/pnas.1708274114
- Chan A-W, Hrobjartsson A, Haahr MT, Gotzsche PC, Altman DG. Empirical evidence for selective reporting of outcomes in randomized trials: comparison of protocols to published articles. JAMA. 2004;291(20):2457-2465. doi:10.1001/jama.291.20.2457 https://doi.org/10.1001/jama.291.20.2457
- Sun X, Briel M, Walter SD, Guyatt GH. Is a subgroup effect believable? Updating criteria to evaluate the credibility of subgroup analyses. BMJ. 2010;340:c117. doi:10.1136/bmj.c117 https://doi.org/10.1136/bmj.c117
Key takeaways
- Exploring data to find new hypotheses is legitimate; HARKing is the undisclosed rewrite of when a hypothesis was born.
- With 20 independent looks at the 5% level and no true effect anywhere, the chance of at least one P below 0.05 is 0.64.
- In the simulated registry, 5 of 20 declared splits reached P below 0.05, yet only frailty changes the odds ratio for early surgery; the dementia and albumin splits carry frailty and different baseline risks, which shift the risk ratio.
- Data-dependent choices of subgroup, cut point or adjustment set inflate false findings even when only one analysis is run.
- Preregister the plan, label post hoc analyses, report every subgroup examined, and confirm a new finding in new data.
Related in the wiki: [[iptw-guide]]