Reproducible Research and Publication-Ready Outputs

2. Publication-Quality Tables — Beyond Copy-Paste

Research Question

"How do I create a regression table suitable for journal submission — with proper formatting, stars, and multiple model columns — directly from my Python output?"

Intuition

Manual tables are the enemy of reproducibility. The moment you copy a coefficient from a Python console into Excel, you have broken the link between your code and your paper. If you re-estimate the model (new data, different controls, corrected specification), you must manually update the table — and you will make errors.

The solution: generate tables programmatically. Python can produce LaTeX (for academic papers), HTML (for websites), and formatted text output directly from model results. The summary_col function in statsmodels creates multi-model comparison tables; for more control, extract coefficients and standard errors yourself and format them.

Python Implementation

import pandas as pd; import numpy as np
import statsmodels.api as sm; import statsmodels.formula.api as smf
from statsmodels.iolib.summary2 import summary_col
import warnings; warnings.filterwarnings('ignore')

df = sm.datasets.get_rdataset("mroz", "wooldridge").data.dropna(subset=['lwage'])

# ── Estimate four model variants ──────────────────────────────────
m1 = smf.ols("lwage ~ educ", data=df).fit()
m2 = smf.ols("lwage ~ educ + exper + expersq", data=df).fit()
m3 = smf.ols("lwage ~ educ + exper + expersq + age + kidslt6 + kidsge6", data=df).fit()
m4 = smf.ols("lwage ~ educ + exper + expersq + age + kidslt6 + kidsge6 + nwifeinc", data=df).fit()

# ── Method 1: summary_col (quick, good for exploration) ───────────
table = summary_col(
    [m1, m2, m3, m4],
    stars=True,
    model_names=['(1) Basic', '(2) +Experience', '(3) +Demographics', '(4) +Husband Inc'],
    info_dict={'N': lambda x: f"{int(x.nobs)}", 'R²': lambda x: f"{x.rsquared:.3f}"}
)
print(table)

# ── Method 2: Manual table (full control for journal submission) ──
def build_reg_table(models, var_list, model_names):
    """Build a formatted regression table string."""
    header = f"{'':<18}" + ''.join(f"{name:>14}" for name in model_names)
    lines = [header, '-' * len(header)]

    for var in var_list:
        row = f"{var:<18}"
        for m in models:
            if var in m.params.index:
                coef = m.params[var]; se = m.bse[var]; p = m.pvalues[var]
                stars = '***' if p < 0.01 else '**' if p < 0.05 else '*' if p < 0.10 else ''
                row += f" {coef:>8.4f}{stars}"
                row += f"\n{'':<18} ({se:>8.4f})"
            else:
                row += f" {'':>14}"
        lines.append(row)

    # Model stats
    n_obs = f"{'Observations':<18}" + ''.join(f" {int(m.nobs):>14}" for m in models)
    r_sq = f"{'R-squared':<18}" + ''.join(f" {m.rsquared:>14.3f}" for m in models)
    lines.extend([n_obs, r_sq, '', 'Robust standard errors in parentheses', '*** p<0.01, ** p<0.05, * p<0.10'])
    return '\n'.join(lines)

var_list = ['educ', 'exper', 'expersq', 'age', 'kidslt6', 'kidsge6', 'nwifeinc']
print("\n─── Publication-Ready Regression Table ───")
print(build_reg_table([m1, m2, m3, m4], var_list,
                       ['(1)', '(2)', '(3)', '(4)']))
            

─── Publication-Ready Regression Table ───
                          (1)           (2)           (3)           (4)
----------------------------------------------------------------------
educ                   0.1086***     0.1079***     0.0988***     0.0952***
                      (0.0144)      (0.0143)      (0.0146)      (0.0148)
exper                                0.0406***     0.0384***     0.0404***
                                    (0.0133)      (0.0134)      (0.0133)
expersq                             -0.0007*      -0.0006       -0.0007*
                                    (0.0004)      (0.0004)      (0.0004)
age                                               -0.0034       -0.0038
                                                  (0.0061)      (0.0062)
kidslt6                                           -0.0966       -0.0407
                                                  (0.0884)      (0.0952)
kidsge6                                           -0.0092       -0.0095
                                                  (0.0329)      (0.0328)
nwifeinc                                                        0.0026
                                                               (0.0033)
Observations              428           428           428           428
R-squared                0.117         0.145         0.157         0.167

Interpretation

The education coefficient is stable across specifications (0.095–0.109) — a strong signal of robustness. Progressive addition of controls increases R² from 0.117 to 0.167. The table is ready for a journal: coefficients with significance stars, standard errors in parentheses, observation counts, and fit statistics. No Excel, no manual entry, no copy-paste errors.

3. Publication-Quality Figures — Every Journal's Requirements

Research Question

"How do I create a figure that meets journal standards: proper resolution, readable fonts, consistent styling, and multi-panel layout?"

Intuition

Journals have specific figure requirements: 300+ DPI for raster images, vector formats preferred (PDF/SVG), readable font sizes (≥8pt at final size), consistent colour palettes that work in greyscale, and proper labelling. Most default matplotlib output fails these standards — default fonts are too small, default DPI is 100, and colours may be indistinguishable when printed.

Journal Figure Checklist

Python Implementation — Complete Publication Figure

import matplotlib.pyplot as plt
import seaborn as sns

plt.style.use('seaborn-v0_8-darkgrid')
plt.rcParams.update({
    'font.size': 10, 'axes.titlesize': 12, 'axes.labelsize': 11,
    'legend.fontsize': 9, 'figure.dpi': 150, 'savefig.dpi': 300,
    'savefig.bbox': 'tight', 'font.family': 'serif'
})

# ── Create a four-panel diagnostic figure ─────────────────────────
df = sm.datasets.get_rdataset("mroz", "wooldridge").data.dropna(subset=['lwage'])
m = smf.ols("lwage ~ educ + exper + expersq", data=df).fit()
resid = m.resid; fitted = m.fittedvalues

fig, axes = plt.subplots(2, 2, figsize=(12, 10))

# Panel (a): Q-Q Plot
from scipy import stats
sm.qqplot(resid, stats.norm, fit=True, line='45', ax=axes[0,0],
          markerfacecolor='steelblue', markeredgecolor='none', alpha=0.5)
axes[0,0].set_title('(a) Q-Q Plot of Residuals', fontweight='bold')

# Panel (b): Residuals vs Fitted
axes[0,1].scatter(fitted, resid, alpha=0.4, edgecolors='none', color='#2c7fb8')
axes[0,1].axhline(y=0, color='#d7191c', linestyle='--', linewidth=1)
axes[0,1].set_xlabel('Fitted Values'); axes[0,1].set_ylabel('Residuals')
axes[0,1].set_title('(b) Residuals vs Fitted', fontweight='bold')

# Panel (c): Coefficient Plot with CIs
coefs = m.params[['educ', 'exper', 'expersq']]; ses = m.bse[['educ', 'exper', 'expersq']]
y_pos = range(len(coefs))
axes[1,0].errorbar(coefs.values, y_pos, xerr=1.96*ses.values, fmt='o',
                    capsize=5, markersize=8, color='#2c7fb8', linewidth=2)
axes[1,0].axvline(x=0, color='gray', linestyle=':', alpha=0.5)
axes[1,0].set_yticks(y_pos); axes[1,0].set_yticklabels(coefs.index)
axes[1,0].set_xlabel('Coefficient Estimate'); axes[1,0].set_title('(c) Coefficient Plot (95% CI)', fontweight='bold')

# Panel (d): Density of Residuals vs Normal
from scipy.stats import gaussian_kde
x_range = np.linspace(resid.min(), resid.max(), 200)
axes[1,1].hist(resid, bins=30, density=True, alpha=0.5, color='steelblue', edgecolor='white')
kde = gaussian_kde(resid)
axes[1,1].plot(x_range, kde(x_range), 'b-', linewidth=2, label='Residuals')
axes[1,1].plot(x_range, stats.norm.pdf(x_range, resid.mean(), resid.std()),
              'r--', linewidth=2, label='Normal')
axes[1,1].set_xlabel('Residuals'); axes[1,1].set_title('(d) Residual Density vs Normal', fontweight='bold')
axes[1,1].legend(fontsize=8)

plt.suptitle('Figure 1: Diagnostic Plots for Mincer Wage Equation', fontweight='bold', fontsize=14, y=1.01)
plt.tight_layout()
plt.savefig('../../assets/images/m10-diagnostic-figure.png', dpi=300)
plt.savefig('../../assets/images/m10-diagnostic-figure.pdf')  # Vector format for journal
plt.show()

print("Figure saved at 300 DPI (PNG) and as vector PDF for journal submission.")
            

Figure saved at 300 DPI (PNG) and as vector PDF for journal submission.

Journal submission tip: Submit PDF for the manuscript (scalable, small file).
Provide high-res PNG/TIFF only if the journal specifically requires raster format.

4. Dynamic Documents — When Your Paper Updates Itself

Research Question

"If I update my dataset or change a specification, can my tables, figures, and manuscript update automatically?"

Intuition

A dynamic document embeds code directly in the manuscript. When you recompile the document, the code runs, producing updated tables and figures that are inserted automatically. No more "Did I remember to update Table 3 after I added that control variable?" The document and the analysis are one.

Three approaches, in order of complexity:

• Jupyter Notebook → LaTeX: Write your analysis in a notebook, export to LaTeX via nbconvert. Good for papers with a single main analysis.
• Quarto: The modern successor to R Markdown. Write in .qmd files mixing markdown and Python; render to HTML, PDF, or Word. Best choice for new projects.
• Papermill: Parameterise notebooks — run the same analysis with different parameters (regions, time periods, specifications) without duplicating code.

Python Implementation — Quarto Example

# Install Quarto (one-time)
# Download from https://quarto.org/docs/get-started/

# Create a new Quarto document
# quarto create project my-paper

# Render to PDF
quarto render my-paper.qmd --to pdf

# Render to HTML
quarto render my-paper.qmd --to html
            

# Example: my-paper.qmd contents
# ---
# title: "Returns to Education — Replication Package"
# author: "Your Name"
# format: pdf
# jupyter: python3
# ---
#
# ## Results
#
# ```{python}
# import statsmodels.formula.api as smf
# import statsmodels.api as sm
# df = sm.datasets.get_rdataset("mroz", "wooldridge").data.dropna(subset=['lwage'])
# model = smf.ols("lwage ~ educ + exper + expersq", data=df).fit()
# print(model.summary())
# ```
#
# The return to education is `{python} round(model.params['educ'], 3)` (p < 0.001).