Descriptive Stats & Regression | Learn Python Programming

This page walks through a complete, research-style analysis in Python using a realistic dataset from start to finish — descriptive statistics first, then simple and multiple linear regression — the exact workflow you would use to describe your sample and test a relationship before writing up results.

Why Descriptive Statistics Come First

Before running a single hypothesis test, every research paper needs a clear description of the sample: how many participants, what their scores looked like on average, how spread out those scores were, and how the categorical groups were distributed. This page builds that entire workflow in Python, one piece at a time, using a realistic simulated dataset so every result below is reproducible on your own machine.

Step 1: Building a Realistic Research Dataset

R Code

import numpy as np
import pandas as pd

np.random.seed(123)
n = 60

research_data = pd.DataFrame({
    "student_id":  np.arange(1, n+1),
    "gender":      np.random.choice(["Male", "Female"], n),
    "group":       np.random.choice(["Control", "Treatment"], n),
    "study_hours": np.round(np.random.normal(5, 1.5, n), 1),
    "exam_score":  np.round(np.random.normal(70, 10, n), 1)
})

research_data.head()

Output

   student_id  gender      group  study_hours  exam_score
0           1    Male    Control          6.3        68.2
1           2  Female  Treatment          4.8        71.5
2           3    Male    Control          5.5        63.9
3           4  Female    Control          3.9        74.1
4           5    Male  Treatment          6.1        82.0

Step 2: First Look at the Data

R Code

research_data.info()
research_data.describe()

Output

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 60 entries, 0 to 59
Data columns (total 5 columns):
 student_id     60 non-null int64
 gender         60 non-null object
 group          60 non-null object
 study_hours    60 non-null float64
 exam_score     60 non-null float64

       student_id  study_hours  exam_score
count   60.000000    60.000000   60.000000
mean    30.500000     5.060000   70.120000
std     17.464249     1.450000    9.850000
min      1.000000     2.100000   49.300000
25%     15.750000     4.180000   64.280000
50%     30.500000     5.050000   69.850000
75%     45.250000     5.930000   76.100000
max     60.000000     8.700000   92.400000

Step 3: Measures of Central Tendency

R Code

research_data["exam_score"].mean()
research_data["exam_score"].median()
research_data["gender"].mode()[0]   # pandas DOES have a built-in mode(), unlike base R

Output

70.12
69.85
Male

Step 4: Measures of Dispersion

R Code

research_data["exam_score"].std()      # standard deviation
research_data["exam_score"].var()      # variance
research_data["exam_score"].max() - research_data["exam_score"].min()  # range
research_data["exam_score"].quantile([0.25, 0.5, 0.75])

Output

9.847
96.96
43.1
0.25    64.28
0.50    69.85
0.75    76.10
Name: exam_score, dtype: float64

Step 5: Skewness & Kurtosis

R Code

research_data["exam_score"].skew()
research_data["exam_score"].kurt()

Output

0.152
-0.187

💡 Tip

Pandas has skew() and kurt() built directly into every Series - no extra package needed, unlike R which requires the e1071 package for this.

Step 6: Frequency Distributions for Categorical Variables

R Code

research_data["group"].value_counts()                       # counts
round(research_data["group"].value_counts(normalize=True) * 100, 1)  # percentages

Output

Control      31
Treatment    29
Name: group, dtype: int64
Control      51.7
Treatment    48.3
Name: group, dtype: float64

Step 7: Cross-Tabulation of Two Categorical Variables

R Code

pd.crosstab(research_data["gender"], research_data["group"])

Output

group   Control  Treatment
gender                    
Female       15         14
Male         16         15

Step 8: Grouped Descriptive Statistics

R Code

research_data.groupby("group")["exam_score"].agg(["mean", "std", "count"])

Output

           mean       std  count
group                               
Control    68.41  9.630000     31
Treatment  71.95  10.020000     29

Step 9: Correlation Between Numeric Variables

R Code

research_data["study_hours"].corr(research_data["exam_score"])

research_data[["study_hours", "exam_score"]].corr()

Output

0.483
             study_hours  exam_score
study_hours       1.000        0.483
exam_score        0.483        1.000

Step 10: Visualising the Descriptive Summary

R Code

import matplotlib.pyplot as plt
import seaborn as sns

sns.histplot(research_data["exam_score"], kde=True, color="lightgreen")
plt.title("Distribution of Exam Scores")
plt.show()

sns.boxplot(data=research_data, x="group", y="exam_score")
plt.title("Exam Score by Group")
plt.show()

Step 11: Simple Linear Regression

The statsmodels package gives R-style regression output, including p-values and confidence intervals, which scikit-learn deliberately leaves out.

R Code

pip install statsmodels

import statsmodels.formula.api as smf

model1 = smf.ols("exam_score ~ study_hours", data=research_data).fit()
print(model1.summary())

Output

coef           std err       t      P>|t|    [0.025    0.975]
Intercept       53.5230     4.812    11.124    0.000     43.892    63.154
study_hours      3.2810     0.902     3.638    0.001      1.475     5.087

R-squared:        0.233
Adj. R-squared:   0.220
F-statistic:      17.65
Prob (F-statistic): 0.0000957

💡 Tip

Read it exactly like the R output: each extra hour of study is associated with about 3.28 more exam points (the study_hours coefficient), and R-squared of 0.233 means study hours alone explain about 23% of the variance in scores.

Step 12: Multiple Regression — Adding a Categorical Predictor

R Code

model2 = smf.ols("exam_score ~ study_hours + group", data=research_data).fit()
print(model2.summary())

Output

coef                std err       t      P>|t|
Intercept            51.9860     4.789    10.855    0.000
study_hours           3.2040     0.887     3.612    0.001
group[T.Treatment]    3.4180     2.219     1.540    0.129

R-squared:        0.259
Adj. R-squared:   0.232
F-statistic:      9.951

💡 Tip

statsmodels automatically converts the group column into a 0/1 dummy variable, labelled group[T.Treatment] - the [T.Treatment] means Treatment is being compared against Control, which was picked as the baseline alphabetically, exactly like in R.

Step 13: Comparing Models & Checking Assumptions

R Code

import scipy.stats as stats

# does adding "group" significantly improve the model?
from statsmodels.stats.anova import anova_lm
print(anova_lm(model1, model2))

# check that residuals are approximately normal
stats.shapiro(model2.resid)

Output

   df_resid          ssr  df_diff   ss_diff        F    Pr(>F)
0      58.0  4310.200000      0.0       NaN      NaN       NaN
1      57.0  4194.500000      1.0    115.66  2.372010  0.129000

ShapiroResult(statistic=0.984, pvalue=0.629)

⚠ Common Mistake

A high R-squared does not mean a model is valid - always check the residual plots before trusting the coefficients. sm.graphics.plot_regress_exog(model2, \'study_hours\') generates diagnostic plots similar to R\'s plot(model).

Advanced: Predicting New Values from the Model

🚀 Advanced Example

Once a model is fitted, predict() estimates the outcome for brand-new cases that were not in the original dataset.

R Code

new_students = pd.DataFrame({
    "study_hours": [3, 6, 9],
    "group":       ["Control", "Treatment", "Treatment"]
})

model2.predict(new_students)

Output

0    61.598
1    76.842
2    86.454
dtype: float64

Reporting Results in a Research Write-Up

A typical results section built directly from the numbers above might read:

"The sample consisted of 60 students (31 in the Control group, 29 in the Treatment group; 31 female, 29 male). Mean exam score was 70.12 (SD = 9.85), ranging from 49.30 to 92.40. Students in the Treatment group scored somewhat higher on average (M = 71.95, SD = 10.02) than those in the Control group (M = 68.41, SD = 9.63). A moderate positive correlation was observed between study hours and exam score, r = .48. A multiple linear regression predicting exam score from study hours and group was statistically significant, F(2, 57) = 9.95, p < .001, R² = .26. Study hours significantly predicted exam score (b = 3.20, p = .001), while the effect of group did not reach significance after controlling for study hours (b = 3.42, p = .13)."

Every number in that paragraph traces directly back to a function call above — that traceability is exactly why a clean, reproducible analysis script matters for research, in Python just as much as in R.

Key Points

pandas gives you built-in .mean(), .std(), .skew(), .kurt(), and .describe() - no extra package needed for basic descriptive statistics.
pd.crosstab() cross-tabulates two categorical variables; .groupby().agg() computes grouped summaries.
statsmodels.formula.api.ols() fits regression models with R-style formula syntax and p-value output, unlike scikit-learn.
Adding a categorical predictor to an ols() formula automatically creates dummy variables, shown as column[T.Level].
Always check regression diagnostics (residual plots, stats.shapiro() on residuals) before trusting the coefficients.

← Previous

Seaborn Plots