Regression Analysis Explained: Coefficients, R² and a Real Example

import numpy as np
import pandas as pd
from datetime import date

# The ScholarsColab Research Methods Index (papers level): anonymised, updated daily
DATA = "https://scholarscolab.in/wp-json/scholarscolab/v1/research-methods-data?level=papers"
papers = pd.read_csv(DATA, comment="#")   # lines starting with # hold the citation

def p_text(p):
    return "< 0.001" if p < 0.001 else f"= {p:.3f}"

print(papers.shape)
papers.head()

THIS_YEAR = date.today().year
df = papers.dropna(subset=["year"]).copy()
df = df[df["year"] >= THIS_YEAR - 30]
df["age"] = THIS_YEAR - df["year"]                  # years since publication

df["log_citations"] = np.log(df["citations"] + 1)   # reduces the effect of a few huge values


df[["age", "citations", "log_citations"]].describe().round(2)

import statsmodels.formula.api as smf
import matplotlib.pyplot as plt

simple = smf.ols("log_citations ~ age", data=df).fit()
print(f"log_citations = {simple.params['Intercept']:.3f} + {simple.params['age']:.3f} × age"
      f"   (R² = {simple.rsquared:.3f})")

sample = df.sample(2000, random_state=1)            # plot a sample so the points stay readable

fig, ax = plt.subplots(figsize=(8, 4.5))
ax.scatter(sample["age"], sample["log_citations"], s=8, alpha=0.25)
x = np.linspace(0, 30, 50)
ax.plot(x, simple.params["Intercept"] + simple.params["age"] * x, color="green", linewidth=3)
ax.set_xlabel("Age of paper (years)")
ax.set_ylabel("log(1 + citations)")
plt.show()

model = smf.ols(
    'log_citations ~ age + review + title_words + title_colon'
    ' + C(venue_group, Treatment("Other journals")) + C(discipline)',
    data=df).fit()

def tidy(name):                                       # shorter names for the table

    name = name.replace('C(venue_group, Treatment("Other journals"))[T.', "venue: ")
    return name.replace("C(discipline)[T.", "discipline: ").rstrip("]")

table = pd.DataFrame({"coef": model.params, "p-value": model.pvalues})
table["% change in citations"] = (np.exp(table["coef"]) - 1) * 100
table.index = [tidy(i) for i in table.index]
print(f"R² = {model.rsquared:.3f}, n = {int(model.nobs):,}")
table.drop(index="Intercept").round(3)

ci = model.conf_int()
ci.index = [tidy(i) for i in ci.index]
show = [i for i in table.index if i.startswith("venue: ")] + ["review", "title_colon"]

pct = (np.exp(table.loc[show, "coef"]) - 1) * 100
low = (np.exp(ci.loc[show, 0]) - 1) * 100
high = (np.exp(ci.loc[show, 1]) - 1) * 100

fig, ax = plt.subplots(figsize=(8, 4))
ax.errorbar(pct, range(len(show)), xerr=[pct - low, high - pct], fmt="o", color="green")
ax.axvline(0, color="grey")
ax.set_yticks(range(len(show)), show)
ax.set_xlabel("% change in citations (95% confidence interval)")
ax.invert_yaxis()
plt.tight_layout()
plt.show()