My Coauthor Reran the Code and Got Different Results

Symptoms

• The same script appears to give a different coefficient, sample size, or sign on a coauthor’s machine.
• A report agrees on fitted predictions but disagrees on the displayed group coefficient.

What this usually means

“The same analysis” includes the input data, selected rows, formula, coding, defaults, random-number configuration, and environment. Compare these intermediate objects before comparing the final table. An opposite sign can represent an opposite reference contrast rather than different predictions.

Common causes

• A refreshed data file or join changes rows or duplicates participants.
• An added variable triggers missing-row exclusion.
• Factor references, units, contrasts, or package defaults differ.
• Random seeds, generator kinds, parallel execution, or software versions differ for stochastic analyses.

Run these checks

1. Match input checksums and record the data version without sharing private records unnecessarily.
2. Compare row counts, stable identifiers, filters, joins, and model-frame exclusions.
3. Compare formulas, model matrices, contrasts, reference levels, weights, and transformations.
4. Record sessionInfo(), package versions, and random-number settings. Start from a clean session.
5. Locate the earliest differing object, then rerun the agreed pipeline from the raw input.

What not to do

Do not average conflicting results or copy the preferred output into the manuscript. A seed alone does not freeze data or software, and identical displayed p-values are weak evidence that the workflows agree.

Treatment options

Create an explicit input-to-report script with saved provenance. Set reference levels and exclusions in code, make data validation failures visible, and lock relevant dependencies. Compare numerical outputs using justified tolerances rather than expecting every platform to produce byte-identical floating-point results.

Worked example

Eight synthetic rows contain one missing covariate. The script separates the sample effect from adjustment and reference coding:

Changing the reference gives identical fitted values: their maximum absolute difference is 0.00000000 at the displayed precision. The example also prints an input checksum, included IDs, the formula, contrast defaults, RNG settings, and session information. The checksum is for this generated CSV; do not treat it as a platform-independent data identity standard.

What to tell the reviewers

We reconciled input versions, included observations, model formulas, and reference coding. The discrepancy arose from a changed analysis sample and contrast direction. We regenerated the report from the agreed pipeline and recorded the data and software versions used.

See it in R and Python

Both languages construct the same deterministic observations and reproduce changes in sample size, adjustment, and contrast direction. Each prints language-appropriate provenance information.

Python dependencies: NumPy and SciPy. Install with python -m pip install numpy scipy.

# Synthetic workflow divergence; base R only. No external data.
d <- data.frame(id=1:8,arm=factor(rep(c('control','treatment'),each=4)),y=c(1,2,3,4,3,4,5,10),z=c(1,2,3,4,1,2,3,NA))
a <- lm(y~arm,d)
b <- lm(y~arm+z,d)
same <- d[complete.cases(d[c('y','arm','z')]),]
c <- lm(y~arm,same)
reverse <- transform(d,arm=relevel(arm,ref='treatment'))
r <- lm(y~arm,reverse)
cat(sprintf('Original model: n=%d, treatment-control=%.3f\n',nobs(a),coef(a)[2]))
cat(sprintf('Same formula on complete rows: n=%d, treatment-control=%.3f\n',nobs(c),coef(c)[2]))
cat(sprintf('Added z: n=%d, conditional treatment-control=%.3f\n',nobs(b),coef(b)[2]))
cat(sprintf('Reversed reference: control-treatment=%.3f; max fitted-value difference=%.8f\n',coef(r)[2],max(abs(fitted(a)-fitted(r)))))
# Stable audit manifest: record input checksum, identifiers, formula and session.
f <- tempfile(fileext='.csv'); write.csv(d,f,row.names=FALSE,quote=FALSE)
cat('Input MD5:',unname(tools::md5sum(f)),'\n'); unlink(f)
cat('Included IDs:',paste(d$id[as.integer(rownames(model.frame(a)))],collapse=','),'\n')
print(formula(a)); print(options('contrasts','na.action')); print(RNGkind()); sessionInfo()
stopifnot(nobs(a)==8,nobs(b)==7,max(abs(fitted(a)-fitted(r)))<1e-8)

# DASS Analysis Clinic Case 009; Python dependencies: numpy, scipy.
from pathlib import Path
import numpy as np
from scipy import stats, optimize
HERE = Path(__file__).resolve().parent

def ols(X, y):
    b = np.linalg.lstsq(X, y, rcond=None)[0]
    residual = y - X @ b
    df = len(y) - X.shape[1]
    cov = (residual @ residual / df) * np.linalg.inv(X.T @ X)
    se = np.sqrt(np.diag(cov))
    ci = np.column_stack((b-stats.t.ppf(.975, df)*se,b+stats.t.ppf(.975,df)*se))
    return b, cov, ci

def power(n, d):
    # Matches R power.t.test(strict=FALSE): rejection tail in effect direction.
    return stats.nct.sf(stats.t.ppf(.975,2*n-2),2*n-2,d*np.sqrt(n/2))

import platform, sys, hashlib
arm=np.repeat([0.,1.],4); y=np.array([1.,2.,3.,4.,3.,4.,5.,10.]); z=np.array([1.,2.,3.,4.,1.,2.,3.,np.nan])
X=np.column_stack((np.ones(8),arm)); full=ols(X,y); keep=~np.isnan(z)
same=ols(X[keep],y[keep]); adjusted=ols(np.column_stack((X[keep],z[keep])),y[keep]); reverse=ols(np.column_stack((np.ones(8),1-arm)),y)
print('All rows treatment-control:',full[0][1]); print('Complete rows treatment-control:',same[0][1]); print('Adjusted complete rows:',adjusted[0][1]); print('Reverse contrast:',reverse[0][1])
print('Included IDs:',np.arange(1,9)); print('Python:',sys.version,'platform:',platform.platform(),'NumPy:',np.__version__)
print('SHA256 of script:',hashlib.sha256(Path(__file__).read_bytes()).hexdigest())
assert np.allclose([full[0][1],same[0][1],adjusted[0][1],reverse[0][1]],[3,1.5,2,-3])
assert np.allclose(X@full[0],np.column_stack((np.ones(8),1-arm))@reverse[0])

Download Python script

Sources

• R Core Team. Fitting linear models: lm documentation.
• R Core Team. Collect information about the current R session.
• R Core Team. Random number generation.