A Data Analysis with pandas¶
From a Messy File to a Defensible Finding¶
STAT 540 — Computing in Statistics Yen-Yi Ho · Department of Statistics, University of South Carolina
Companion notebook to the lecture slides A Data Analysis with pandas.
Lecture 1 taught the pandas vocabulary — Series, DataFrame, loc, groupby, merge.
This notebook is about putting that vocabulary to work on a real question, in the order a
real analysis actually happens.
The data is the tips dataset: every check a waiter recorded over about ten weeks
(Bryant & Smith, 1995). You met a simulated version of it in Lecture 1. This is the
real thing — and it arrives in the state real data arrives in.
The brief¶
A restaurant owner hands you one file and three questions:
- What is a typical tip here? Does the "15 to 20 percent" norm hold?
- Does the tip rate depend on the size of the bill?
- Do smokers, men, women, weekends — any of the groups servers argue about — tip differently?
Deliverable: a short written finding for each, with the evidence behind it and an honest statement of what the data cannot settle.
The eight stages¶
| Stage | The question it answers | |
|---|---|---|
| 1 | Load and inspect | What did I actually get? |
| 2 | Clean | What is wrong with it, and what will I do about it? |
| 3 | Derive | What do I need that is not in the file? |
| 4 | One variable | How does the outcome behave on its own? |
| 5 | Two variables | What moves with what? |
| 6 | Compare groups | Is the difference real, or is it something else? |
| 7 | Visualize | What is the one picture that shows it? |
| 8 | Report | What can I claim, and what can I not? |
Stages 1–3 take most of the time. That is normal, and it is not wasted time.
Setup¶
The cell below writes the data file next to this notebook if it is not already there, so the notebook runs anywhere. It is not part of the lesson — collapse it and move on.
import base64, zlib, pathlib
_BLOB = """
eNqNWsmOHMcRvfsr+jDHciEjcj8KtgYWoOGBog48ETQ5BAfi9AicoWH66x1LVldu1eJFIMRgZcb2
4r3I/vD5/sMf7x4+Li9PL++/vPv3w5cvy8vDn8vz/X+X58enP+6/Lh/ff6f/9Xi/PD/8j/5z//U/
91/fnZ9e7v/2D2OMW27QriEtdrWw3L3/cr+8elp++3Ze/vlwPtM/x4Xtkl1uwKwmLrDSn2/vH4vl
m8/fvi5fvp0/fF6ALdEBfTGukBZcjdks394/L7+9f1k+Nh8lU8irc3Q4md5thm/+9fvr5Vf5phii
9cvpBvzq7ak2lS8210S+pl9jpMPxcs0zOfTTm8YSAJebuKInf/Zbsuf7J8UfCIHORrd65LNdqmN0
2kytmKZEptasLpzU+buffv1ZPHre4+nFI/rLG7CrBzW8xOiFw1m5DnTgDVDMOfBhD/xyetNbWstB
ghXiiWy9WT61SapsTY6SJEyLayJ/+/WhDZP3kiIJ6GbY1ofTw8kfxDUh2cEeoyqTEiLEzJe0dDKH
KKKGiCyfu1SCy8sNHZzIl2gvH3zZy61kXC9IB1s5uMuhGgVyF8Lq/PL3nCXfMzOQpCClT+OnZ2pM
fnn16ufX5WtcOtavGBbf1c4QlSznkpFGb1rdxnDikOoQuML4YPnY76+a6sISYuAQU6V1jdX467kQ
kRyGUgrT+lanMUuB2bRsBXOeFIynxrJupY4lx7f4kbsfK3eNB4mfd3RD3xVLHcDMEJFWKsK9mb/P
qi9wciN3ia0izabV2aWfg5ha0BDeTlGMevkGCcUY7cDvF9xjoo0MXAiOP+Z6DCMcaTCMvsWRTvJJ
3GqmcllDE0BOpmDvzUbHnroC5FykFRkQY2ju1zhrUGDOphKXS++2V8tiRRAPaw6XThtDx+iKhJpU
Le6woq22W5IUh1zh1qyqUhR/s9nS8bYkrflaohsStEauqrTno68+pIAI/ud0aiqBbNnrCmBMTIps
pkK2DTfUTluTa5C+SI5or/M06eoU6Eo3zqyEgba071hSUWueMgJc0E1J1TXvonQvDSVf3+s0fPDH
0QXpk4nmEp1MU6LMhaFzo5daJgxF7dzzU6m86noITuYRRm3da+DCUZFccAnYw6T9cJebmCU0VPC2
Cc34RS0+OnTHZ6n81790CO2kj5xtvzgAVkLB58Anu4K7w6kFDCI3G+7+Pg+JU7AXf/0W6K7i0Uiv
GWEy7eitiFRJL4MVTTauZDQ95dIbBrGMMiqDadGqniC4fHh6/JNvaVFCQ4yqxei9XsvpIGNBGp2J
YQGPoTWtVTvGXRp3E0jV+mIWRbnL5nSZwVdqDLRokyQQax5X+yR9BfxhmrGOCQUBxd1muXukyODk
CoSu+TSWRQNeUeE1sU8Q9hSdhnta5Zsm9LA0AGL2Uh1+Y1zSCiOHc1LhlHpsvzfwQp+lYS3nh4Jf
RaivSx7wxC4il4fdIHvgWyZoz5D1Pns61q6VEVBHQOgqiDLTjBVqKKF6IGyYYOgKH+XpTXf03GMZ
l8cjGeJQpryTs92xXPFKapBHX+9MbYjUDiedA6cOpob6sUoibU3Zz81wVrssLc4e8yeJn95uc6ot
Mmowb7jIWFeZ6R0LZ0ZJtUPVH1Xf7MOvuI2i0lCK4mB2qN8mCm3xnG0Hndc1FYpGPhmsev3pStMS
x+BxCXGRcXMRP5XrBXaDkkmehWTd1sWpBUErQIBdsY0MwcosZLVA9aNE+9MlOrDBHwJFnRSaCZWe
+d5xsZJGlKnkssL5dsXTKJFArmhidcXzhGYFK/WD0Pf2OCSCypWoE+d2nw8b0ynuQMiSHc8sZkeB
aXZs0vEkkBq5KMdkB2VtMVRQOgVIBTMan+0MG8rMWenXxBwlTgl3ESGXEgfeRdweIgVIbIj6skYL
7rizwTvhlt7/lXL1inuhNGtXjE17GTE17vr8QAtytIHq6O6O2ggmKPPmkSQ0+OiWxEW1eASdoZEZ
dQKzk+g44WY4JRcljKw5HQ+EU2F8U+wBTxXhEnepb/j+PoaVgoDukbLQCn9cORCF1GShwfGQmQHv
m5xnmmIbqicHV6gMRkuW2tRxV4+LIZ1FgeHB+l4NDzjCcs4mhi87pq7mR0TkROBIC+wVO2g1RBUI
EJR/lwt2yA2eZYuE2HKjdOBQK8ms49xAd8OBVrMEo54ytm2AYfYGZXCqYFKB7YmiAwHtwBkJqULO
vRSK1JVlHan2zKUV3KHMNqg3tP1yqdfFCFKtlEBizbkzbbuZ2LLckvO31eEUmUD3ar7ov8dDjgdi
adjxvygeUC5K9CSWRc4VxsPbHMgazm0hNXAy9oa+J/O8QoeRNqKqCowc8qP9bKFQHCFgpdcJ0b71
mWtZoYyByePdgfowvqwPfGAQs+5wB+KcyBSe0RYPdRREJxdMvmz0SkcPVDCLvFR89/t8HqYpL+FI
UFNUqC5itUEd6jEYLTMOYxAKw4P8NDBgdgQjn9pI9DrWJTTiMYFykTyHfFVAIl7Q5HhZ53VDyfju
8Yo+QB73VmqBN8JHC3MUwWF4vwvMN6brU5DllRS/rfaiw4yKcXPCshN3c8ZtrMKSGd4HxlmvoIi8
hsN8hcXLVEm80rNHA9frGFW9JjzDx5qrndrdGX+R4mI5fFXRTDZiTgPtuE2u8QxeTBFNdKneBQ/z
EY1SB7bKeMBGNNzMmbCQ/XZkDKv0KLwlsBoi4jQFh8If8jQ5E8TRpxuT6lX5bJ6Ciu7LbvGywW34
rAxKDbnJQ2H0MocRgusiNf13oNCBW9AzeLsrBDTzxtwR4NA9/b5TfdOtDQ33CwVd1iLuIOjuojb4
dY2JvBm2amOgftgUfFkMMfwQ77rt9iL6wIRG36LAibDNB8qtqCyljbnmyufJBg6dLgli6vcdYtmC
hm4UExzyb+1x3TMRNdurjfLdtiN49oYfAk7t1wbYA33CUQVjO0HUsQqnCi/EpiN7bEEVY/1mtFKX
hUZh6++mfseKDNriLDiCOWZHLFVZzDM0R7iiNqLuY4w+Zd5d3inaJatuTkxsMncez0UmCFRehim6
aOnJ5VAGURFsPrUstFnZallb1OVT19RtwRjtgFyneMRRWbMi73aBd6PX6JiqVMcE2F5n/foOITne
VtqT12Wj6xXLOBrCkQhDfiNhSeI7mdirv+RkcoRUaau3kw1r2J/eKeyXZ5wBwgK/aHsWxs51W/Sa
AKakgjdVgHyeGILViQW+Wf1MN/Oihpyuk3gP0WnEgp36unVZKU/oftl/oC6Lo+x0DYxgXJrd69iP
NHr90TNdWd4VURLlBQnjlPeWOFr5qpNXzLArg/7dFrIRzQtmsgX5WEEckw7SlCa2T1KDeE/C4q2o
HM/KoFW9jRhTub15PpM52i1eFzCyP4tu+uykAWKhz5fsZeVpKFvZCPqy3cjT1WGJD+rTPNM36tor
I81wq/IPJww7vj2A974krR8oO4HuNbhVlkYFUZAdrPWzdw61NPKINv19Q/2ixDOXf6oi28DZ/qWg
rPYqWWivHtANI1MvrJnVp6m396MasuV3Kp4dfzzSs2Blw7EdfAUVkZkO6O9Z8CI9B7NC5lPud1jD
m3pWGmq1VQ8ZVkgi6yJXg7+2j7sw7xp1Jr9JQJ4Ekat6B7wJaDvVG1mWduH44RWg7HFN1QHnSe6k
7YEXOo3dgLOcOf5lANsJTZv+Boqpl808SONlkM4WG+yGA+4jf7DN8dokYSOoatfFWGUQQ5xn0QDd
zyqGoSy/D8PZ28PwiyqvBAi2fZP2R//OjELdE18PwvF+zUi56NCbvBWqJWgE9ZcfQrysm0dw25fn
uE1wHB9fm/IPMpvl9wsA8yXE/wF3CT/3
"""
_p = pathlib.Path("tips_messy.csv")
if not _p.exists():
_p.write_bytes(zlib.decompress(base64.b64decode("".join(_BLOB.split()))))
print(_p, "-", _p.stat().st_size, "bytes")
tips_messy.csv - 10108 bytes
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
pd.set_option("display.max_rows", 12)
pd.set_option("display.precision", 3)
# A two-colour palette, checked for colour-vision deficiency separation.
BLUE, ORANGE, GREY = "#2a78d6", "#eb6834", "#52514e"
plt.rcParams.update({
"figure.figsize": (7, 3.2), "font.size": 10, "figure.dpi": 110,
"axes.spines.top": False, "axes.spines.right": False,
"axes.grid": True, "grid.alpha": 0.25, "axes.axisbelow": True,
})
print("pandas", pd.__version__)
pandas 2.3.3
Two ground rules¶
1. Never modify the raw file. Read it, copy it, transform the copy. The raw file is your evidence; if you overwrite it you cannot reproduce or check anything.
2. Every cleaning decision is a research decision. Dropping a row, recoding a label, filling a missing value — each one changes your answer. Write each down as you make it, because you will have to defend it.
raw = pd.read_csv("tips_messy.csv")
raw.shape
(247, 9)
247 rows and 9 columns. The owner said "about 240 checks", so the row count is plausible — but we were promised 8 columns, not 9.
Always check shape first. If it surprises you, stop and find out why before doing
anything else.
raw.head()
| check_id | total_bill | tip | sex | smoker | day | time | size | server_note | |
|---|---|---|---|---|---|---|---|---|---|
| 0 | C0004 | $23.68 | 3.31 | Male | No | Sun | Dinner | 2.0 | NaN |
| 1 | C0083 | $10.07 | 1.83 | Female | No | Thur | lunch | 1.0 | NaN |
| 2 | C0241 | $27.18 | 2.00 | Female | Yes | Sat | dinner | 2.0 | NaN |
| 3 | C0081 | $19.44 | 3.00 | Male | Yes | THUR | Lunch | 2.0 | NaN |
| 4 | C0235 | $15.53 | 3.00 | Male | Y | Sat | Dinner | 2.0 | NaN |
Three problems are visible in the first five rows alone: a $ inside total_bill,
a server_note column nobody mentioned, and inconsistent capitalisation in day and
time.
1.2 info() — the most informative single call¶
raw.info()
<class 'pandas.core.frame.DataFrame'> RangeIndex: 247 entries, 0 to 246 Data columns (total 9 columns): # Column Non-Null Count Dtype --- ------ -------------- ----- 0 check_id 247 non-null object 1 total_bill 247 non-null object 2 tip 244 non-null float64 3 sex 247 non-null object 4 smoker 247 non-null object 5 day 247 non-null object 6 time 247 non-null object 7 size 243 non-null float64 8 server_note 5 non-null object dtypes: float64(2), object(7) memory usage: 17.5+ KB
1.3 What dtypes is telling you¶
raw.dtypes
check_id object total_bill object tip float64 sex object smoker object day object time object size float64 server_note object dtype: object
total_bill is text, not a number. pandas read it as object/str because of the
dollar signs. raw["total_bill"].mean() would fail outright — or, worse, in some
versions concatenate the strings.
A numeric column arriving as text is the single most common data-loading defect, and
dtypes catches it in one line.
# proof, if you want to see it fail
try:
raw["total_bill"].mean()
except TypeError as e:
print("TypeError:", str(e)[:70])
TypeError: Could not convert string '$23.68$10.07$27.18$19.44 $15.53 $15.77$7.25
1.4 describe() on what little is numeric¶
raw.describe()
| tip | size | |
|---|---|---|
| count | 244.000 | 243.000 |
| mean | 1.771 | 2.576 |
| std | 11.350 | 0.956 |
| min | -99.000 | 1.000 |
| 25% | 2.000 | 2.000 |
| 50% | 2.960 | 2.000 |
| 75% | 3.603 | 3.000 |
| max | 10.000 | 6.000 |
Only two columns are numeric — and tip has a minimum of −99, which is impossible
for a tip.
That is a sentinel value: a code someone chose to mean "missing". −99, −999, 9999,
"N/A", ".", −1 — every institution has its own, and none of them are missing to
pandas until you say so. describe() finds them, because an impossible minimum or
maximum sticks out immediately.
1.5 Inspect every categorical column¶
for c in ["sex", "smoker", "day", "time"]:
print(f"{c}: {sorted(raw[c].dropna().unique())}\n")
sex: ['F', 'Female', 'M', 'MALE', 'Male', 'female', 'male'] smoker: ['N', 'No', 'Y', 'Yes', 'no', 'yes'] day: [' Fri', ' Sat', ' Sun', ' Thur', 'FRI', 'Fri', 'Fri ', 'SAT', 'SUN', 'Sat', 'Sat ', 'Sun', 'Sun ', 'THUR', 'Thur', 'Thur ', 'fri', 'sat', 'sun', 'thur'] time: ['DINNER', 'Dinner', 'Lunch', 'dinner', 'lunch']
Seven spellings for two sexes. Twenty for four days. A groupby("day") right now would
produce twenty groups and a meaningless answer.
1.6 Missing values and duplicates¶
print(raw.isna().sum())
print("\nduplicated rows:", raw.duplicated().sum())
check_id 0 total_bill 0 tip 3 sex 0 smoker 0 day 0 time 0 size 4 server_note 242 dtype: int64 duplicated rows: 3
raw.loc[~raw["server_note"].isna(), "server_note"]
62 comp 93 comp 98 comp 171 comp 244 comp Name: server_note, dtype: object
1.7 The five questions to ask any new dataset¶
| Question | One line that answers it | |
|---|---|---|
| 1 | How big is it? | df.shape |
| 2 | What type is each column? | df.dtypes / df.info() |
| 3 | What is missing? | df.isna().sum() |
| 4 | Are there repeated rows? | df.duplicated().sum() |
| 5 | What values does each category take? | df[c].value_counts() |
Run all five on every dataset, every time, before you analyze anything. They take thirty seconds and they will save you hours.
1.8 The defect inventory¶
| # | Defect | Fix |
|---|---|---|
| 1 | server_note is empty and unused |
drop the column |
| 2 | total_bill is text with $ |
strip, convert to float |
| 3 | tip uses −99 for missing |
replace with NaN |
| 4 | sex, smoker have many spellings |
harmonise |
| 5 | day, time have case/space noise |
strip, title-case |
| 6 | 3 duplicated rows | drop_duplicates |
| 7 | 6 missing tips, 4 missing sizes | decide, then document |
Write the inventory before you write any cleaning code. It becomes your checklist, and later it becomes a paragraph of your report.
d = raw.copy()
d.shape
(247, 9)
From here on raw is untouched evidence and d is the working copy. If a cleaning step
goes wrong you re-run from raw — you never re-download.
2.2 Defect 1: a column nobody needs¶
raw["server_note"].value_counts(dropna=False)
server_note NaN 242 comp 5 Name: count, dtype: int64
d = d.drop(columns="server_note")
d.columns.tolist()
['check_id', 'total_bill', 'tip', 'sex', 'smoker', 'day', 'time', 'size']
242 blanks and 5 comp flags — not enough to analyze. We drop it, and we say so in the
report. Silently discarding a column is how findings become irreproducible.
2.3 Defect 2: money stored as text¶
raw["total_bill"].head(5).tolist()
['$23.68', '$10.07', '$27.18', '$19.44', ' $15.53 ']
d["total_bill"].isna().sum() # 0
np.int64(0)
d["total_bill"] = (d["total_bill"]
.str.strip() # kill the padding
.str.replace(r"[$,]", "", regex=True) # kill $ and thousands ,
.astype(float))
d["total_bill"].head(5).tolist()
[23.68, 10.07, 27.18, 19.44, 15.53]
Chain the .str methods, then astype. Removing , matters the moment a bill reaches
$1,000 — and it costs nothing to handle now.
Always check what a type conversion did¶
print("new NaNs:", d["total_bill"].isna().sum())
print("range:", d["total_bill"].min(), "-", d["total_bill"].max())
new NaNs: 0 range: 3.07 - 50.81
astype(float) raises if a value will not parse, so a clean run means every value
converted. But pd.to_numeric(..., errors="coerce") — which you will often want, because
it survives one bad row — turns failures into NaN silently. That is why you count
them afterwards either way.
A range of $3.07 to $50.81 for restaurant checks is plausible. Good.
2.4 Defect 3: the sentinel value¶
print("rows coded -99:", (raw["tip"] == -99).sum())
d["tip"] = d["tip"].replace(-99, np.nan)
print("tip min now:", d["tip"].min())
print("tip NaNs now:", d["tip"].isna().sum())
rows coded -99: 3 tip min now: 1.0 tip NaNs now: 6
Six missing tips in total: three written as -99 and three left blank. Two codings for
the same thing in one column — entirely typical.
2.5 Defect 4: one category, many spellings¶
sorted(raw["smoker"].unique())
['N', 'No', 'Y', 'Yes', 'no', 'yes']
You could write a dictionary with every spelling in it. That works until next month's
export contains "YES " and quietly produces a NaN.
Better: normalise first, then map. Strip the whitespace, upper-case it, take the first letter — now there are only two possible values left to map, and the mapping is robust to spellings you have not seen.
d["sex"] = d["sex"].str.strip().str.upper().str[0].map({"M": "Male", "F": "Female"})
d["smoker"] = d["smoker"].str.strip().str.upper().str[0].map({"N": "No", "Y": "Yes"})
d["day"] = d["day"].str.strip().str.title()
d["time"] = d["time"].str.strip().str.title()
for c in ["sex", "smoker", "day", "time"]:
print(f"{c}: {sorted(d[c].unique())}")
sex: ['Female', 'Male'] smoker: ['No', 'Yes'] day: ['Fri', 'Sat', 'Sun', 'Thur'] time: ['Dinner', 'Lunch']
The one dangerous property of map¶
map returns NaN for anything not in the dictionary — silently. A spelling you did
not anticipate does not raise an error; it becomes a missing value and disappears from
your groupby. So after every map on a category, count the NaNs.
d[["sex", "smoker"]].isna().sum()
sex 0 smoker 0 dtype: int64
Zero. Nothing escaped the mapping.
2.6 Defect 5: duplicate rows¶
d[d.duplicated(keep=False)].sort_values("check_id")
| check_id | total_bill | tip | sex | smoker | day | time | size | |
|---|---|---|---|---|---|---|---|---|
| 20 | C0012 | 35.26 | 5.0 | Female | No | Sun | Dinner | 4.0 |
| 48 | C0012 | 35.26 | 5.0 | Female | No | Sun | Dinner | 4.0 |
| 28 | C0091 | 28.97 | 3.0 | Male | Yes | Fri | Dinner | 2.0 |
| 53 | C0091 | 28.97 | 3.0 | Male | Yes | Fri | Dinner | 2.0 |
| 172 | C0204 | 16.40 | 2.5 | Female | Yes | Thur | Lunch | 2.0 |
| 173 | C0204 | 16.40 | 2.5 | Female | Yes | Thur | Lunch | 2.0 |
before = len(d)
d = d.drop_duplicates()
print(f"{before} -> {len(d)} rows")
247 -> 244 rows
Three exact duplicates removed.
Be careful here. Two different checks could legitimately have identical bills,
tips, and party sizes. What proves these are duplicates is the repeated check_id.
Without an identifier, drop_duplicates() on the full row would destroy real data.
2.7 Defect 6: missing values — count first¶
print(d.isna().sum())
print("\nas a percentage:")
print((100 * d[["tip", "size"]].isna().mean()).round(1))
check_id 0 total_bill 0 tip 6 sex 0 smoker 0 day 0 time 0 size 4 dtype: int64 as a percentage: tip 2.5 size 1.6 dtype: float64
About 2.5% of tips and 1.7% of sizes. Small — but "small" is not a decision.
Then decide¶
| Option | When it is right | Here? |
|---|---|---|
| Drop the rows | few, and not systematic | yes |
| Impute | many, and you must keep the sample size | no |
| Keep as its own category | the missingness is itself informative | no |
tip is our outcome variable. Imputing an outcome means inventing the very thing we
are trying to measure, so that is out. Under 3% missing, so we drop — and we report that
we dropped 10 rows.
Check who you are about to drop¶
This is the step people skip. If every missing tip were on a Saturday, dropping them would bias the Saturday estimate.
gone = d[d[["tip", "size"]].isna().any(axis="columns")]
print(gone["day"].value_counts(), "\n")
print(gone[["check_id", "total_bill", "tip", "day", "size"]])
day
Sat 5
Sun 3
Thur 1
Fri 1
Name: count, dtype: int64
check_id total_bill tip day size
18 C0067 16.45 NaN Sat 2.0
25 C0129 11.38 NaN Thur 2.0
34 C0061 20.29 NaN Sat 2.0
65 C0233 11.61 3.39 Sat NaN
66 C0216 12.90 NaN Sat 2.0
118 C0107 20.49 4.06 Sat NaN
147 C0054 9.94 NaN Sun 2.0
178 C0164 13.81 2.00 Sun NaN
179 C0101 11.35 NaN Fri 2.0
209 C0190 23.10 4.00 Sun NaN
Spread across all four days, with ordinary-looking bills. Dropping is defensible.
d = d.dropna(subset=["tip", "size"]).copy()
d["size"] = d["size"].astype(int)
d.shape
(234, 8)
234 checks from the original 247: −3 duplicates, −10 incomplete.
2.8 Put the whole thing in a function¶
Everything above, in one reusable place.
def load_tips(path="tips_messy.csv"):
"""Read the register export and return analysis-ready data."""
d = pd.read_csv(path).drop(columns="server_note")
d["total_bill"] = (d["total_bill"].str.strip()
.str.replace(r"[$,]", "", regex=True).astype(float))
d["tip"] = d["tip"].replace(-99, np.nan)
d["sex"] = d["sex"].str.strip().str.upper().str[0].map({"M": "Male", "F": "Female"})
d["smoker"] = d["smoker"].str.strip().str.upper().str[0].map({"N": "No", "Y": "Yes"})
d["day"] = d["day"].str.strip().str.title()
d["time"] = d["time"].str.strip().str.title()
d = d.drop_duplicates().dropna(subset=["tip", "size"]).copy()
d["size"] = d["size"].astype(int)
return d.reset_index(drop=True)
Why a function? It runs the same way every time, it is one line to re-run when next month's export arrives, and it is the honest answer to "what did you do to the data?"
2.9 Validate — state your expectations as code¶
tips = load_tips()
assert tips["total_bill"].between(1, 500).all(), "implausible bill"
assert tips["tip"].ge(0).all(), "negative tip"
assert tips["tip"].le(tips["total_bill"]).all(), "tip exceeds bill"
assert set(tips["sex"]) == {"Male", "Female"}, "unexpected sex"
assert set(tips["smoker"]) == {"No", "Yes"}, "unexpected smoker"
assert set(tips["day"]) == {"Thur", "Fri", "Sat", "Sun"}, "unexpected day"
assert tips["size"].between(1, 20).all(), "implausible party size"
assert not tips.isna().any().any(), "missing values remain"
assert tips["check_id"].is_unique, "duplicate check_id"
print(f"all checks passed on {len(tips)} rows")
all checks passed on 234 rows
These assertions are the cheapest insurance in data analysis. They cost one line each and they fail loudly next month when the register starts exporting something new.
Notice that they encode subject-matter knowledge, not just types: a tip cannot exceed the bill, a party cannot have 400 people. Those are the errors that type checks miss.
print(tips.dtypes, "\n")
tips.head()
check_id object total_bill float64 tip float64 sex object smoker object day object time object size int64 dtype: object
| check_id | total_bill | tip | sex | smoker | day | time | size | |
|---|---|---|---|---|---|---|---|---|
| 0 | C0004 | 23.68 | 3.31 | Male | No | Sun | Dinner | 2 |
| 1 | C0083 | 10.07 | 1.83 | Female | No | Thur | Lunch | 1 |
| 2 | C0241 | 27.18 | 2.00 | Female | Yes | Sat | Dinner | 2 |
| 3 | C0081 | 19.44 | 3.00 | Male | Yes | Thur | Lunch | 2 |
| 4 | C0235 | 15.53 | 3.00 | Male | Yes | Sat | Dinner | 2 |
2.10 Which checks did we lose?¶
We went from 247 rows to 234. The owner will ask which ones — and answering that is a comparison between two columns living in two different DataFrames.
lost = pd.Index(raw["check_id"]).difference(tips["check_id"])
print(len(lost), "checks dropped:")
print(lost.tolist())
10 checks dropped: ['C0054', 'C0061', 'C0067', 'C0101', 'C0107', 'C0129', 'C0164', 'C0190', 'C0216', 'C0233']
That is R's setdiff(raw$check_id, tips$check_id).
Comparing two columns from two frames¶
Wrap a column in pd.Index(...) and you get the set operations. The Index type carries
them; a plain Series does not.
| R | pandas | Returns |
|---|---|---|
a %in% b |
a.isin(b) |
Boolean mask, aligned to a |
which(a %in% b) |
np.flatnonzero(a.isin(b)) |
positions in a |
a[a %in% b] |
a[a.isin(b)] |
the matching values |
match(a, b) |
pd.Index(b).get_indexer(a) |
position in b, -1 if absent |
match(b, a) |
pd.Index(a).get_indexer(b) |
position in a, -1 if absent |
intersect(a, b) |
pd.Index(a).intersection(b) |
in both |
setdiff(a, b) |
pd.Index(a).difference(b) |
in a, not b |
setdiff(b, a) |
pd.Index(b).difference(a) |
in b, not a |
union(a, b) |
pd.Index(a).union(b) |
either |
a = pd.Series(["a", "b", "c"])
b = pd.Series(["b", "f", "g", "a", "i"])
print("a %in% b ", a.isin(b).tolist())
print("which(a %in% b) ", np.flatnonzero(a.isin(b)))
print("a[a %in% b] ", a[a.isin(b)].tolist())
print("match(a, b) ", pd.Index(b).get_indexer(a))
print("match(b, a) ", pd.Index(a).get_indexer(b))
print("intersect(a, b) ", pd.Index(a).intersection(b).tolist())
print("setdiff(a, b) ", pd.Index(a).difference(b).tolist())
print("setdiff(b, a) ", pd.Index(b).difference(a).tolist())
print("union(a, b) ", pd.Index(a).union(b).tolist())
a %in% b [True, True, False] which(a %in% b) [0 1] a[a %in% b] ['a', 'b'] match(a, b) [ 3 0 -1] match(b, a) [ 1 -1 -1 0 -1] intersect(a, b) ['a', 'b'] setdiff(a, b) ['c'] setdiff(b, a) ['f', 'g', 'i'] union(a, b) ['a', 'b', 'c', 'f', 'g', 'i']
Read match(a, b) → [3, 0, -1] element by element: "a" sits at position 3 of b,
"b" at position 0, and "c" is absent, so -1. R would say [4, 1, NA] — the same
answer, 1-based, with NA instead of -1.
In Python
-1is a legal position meaning "last element". Filter on>= 0before usingget_indexeroutput to index anything, or you will silently pick up the wrong row.
pos = pd.Index(b).get_indexer(a)
ok = pos >= 0
print("all positions :", pos)
print("safe to use :", pos[ok])
print("the values :", b.iloc[pos[ok]].tolist()) # .iloc -- these are POSITIONS
all positions : [ 3 0 -1] safe to use : [3 0] the values : ['a', 'b']
Comparison between two Series aligns on the index, so it requires identical labels — and columns from two different frames rarely have them.
isin, get_indexer, and the Index set operations all sidestep this: they compare
values as a set, ignoring both the index and the length. That is precisely why they
are the right tools for cross-DataFrame work.
print("rows dropped :", len(raw) - len(tips))
print("ids dropped :", len(lost))
rows dropped : 13 ids dropped : 10
Thirteen rows disappeared but only ten check numbers. The other three were the
duplicate rows — and a duplicate shares its check_id with a row that survived, so that
id never left the set.
Set operations work on values, not on rows. When the two counts disagree, the disagreement is itself informative: here it confirms exactly the decomposition we expect, 13 = 10 incomplete + 3 duplicated.
And which rows were they?¶
difference finds which ids; isin pulls those rows back. The two operations pair up.
raw.loc[raw["check_id"].isin(lost), ["check_id", "total_bill", "tip", "size", "day"]]
| check_id | total_bill | tip | size | day | |
|---|---|---|---|---|---|
| 18 | C0067 | $16.45 | -99.00 | 2.0 | Sat |
| 25 | C0129 | $11.38 | NaN | 2.0 | Thur |
| 34 | C0061 | $20.29 | NaN | 2.0 | Sat |
| 65 | C0233 | $11.61 | 3.39 | NaN | Sat |
| 66 | C0216 | $12.90 | -99.00 | 2.0 | Sat |
| 118 | C0107 | $20.49 | 4.06 | NaN | Sat |
| 147 | C0054 | $9.94 | -99.00 | 2.0 | Sun |
| 178 | C0164 | $13.81 | 2.00 | NaN | Sun |
| 179 | C0101 | $11.35 | NaN | 2.0 | Fri |
| 209 | C0190 | $23.10 | 4.00 | NaN | Sun |
tips["tip_pct"] = tips["tip"] / tips["total_bill"]
tips[["total_bill", "tip", "tip_pct"]].head()
| total_bill | tip | tip_pct | |
|---|---|---|---|
| 0 | 23.68 | 3.31 | 0.140 |
| 1 | 10.07 | 1.83 | 0.182 |
| 2 | 27.18 | 2.00 | 0.074 |
| 3 | 19.44 | 3.00 | 0.154 |
| 4 | 15.53 | 3.00 | 0.193 |
This is the most consequential decision in the whole analysis, and it is made before
any computation. Analyze tip and you will conclude that large parties tip well; analyze
tip_pct and you will conclude the opposite. Both are true. They answer different
questions.
3.2 Several columns at once with assign¶
tips = tips.assign(
bill_pp = lambda d: d["total_bill"] / d["size"],
is_weekend = lambda d: d["day"].isin(["Sat", "Sun"]),
bill_q = lambda d: pd.qcut(d["total_bill"], 4, labels=["Q1", "Q2", "Q3", "Q4"]),
)
tips[["total_bill", "size", "bill_pp", "is_weekend", "bill_q"]].head()
| total_bill | size | bill_pp | is_weekend | bill_q | |
|---|---|---|---|---|---|
| 0 | 23.68 | 2 | 11.840 | True | Q3 |
| 1 | 10.07 | 1 | 10.070 | False | Q1 |
| 2 | 27.18 | 2 | 13.590 | True | Q4 |
| 3 | 19.44 | 2 | 9.720 | False | Q3 |
| 4 | 15.53 | 2 | 7.765 | True | Q2 |
assign takes lambda d: ... so each new column can refer to the frame being built, and
it returns a new frame rather than mutating in place — which is what lets you chain it.
3.3 Why bin the bill?¶
tips.groupby("bill_q", observed=True)["total_bill"].agg(["count", "min", "max"]).round(2)
| count | min | max | |
|---|---|---|---|
| bill_q | |||
| Q1 | 61 | 3.07 | 13.42 |
| Q2 | 56 | 13.51 | 17.82 |
| Q3 | 58 | 17.89 | 24.52 |
| Q4 | 59 | 24.55 | 50.81 |
qcut makes four groups of roughly equal count. We will use these quartiles over and
over to ask the question that separates a good analysis from a bad one:
Is this group difference still there once I compare like with like?
print(tips.shape)
print(tips.columns.tolist())
(234, 12) ['check_id', 'total_bill', 'tip', 'sex', 'smoker', 'day', 'time', 'size', 'tip_pct', 'bill_pp', 'is_weekend', 'bill_q']
tips["tip_pct"].describe().round(4)
count 234.000 mean 0.160 std 0.061 min 0.036 25% 0.127 50% 0.154 75% 0.191 max 0.710 Name: tip_pct, dtype: float64
Median 15.4%, mean 16.0%. Mean above median means a right skew — a few very large tip rates pulling the average up. The maximum of 71% confirms it.
4.2 A picture beats the six numbers¶
fig, ax = plt.subplots()
ax.hist(tips["tip_pct"], bins=30, color=BLUE, edgecolor="white", linewidth=0.5)
ax.axvline(tips["tip_pct"].median(), color=ORANGE, lw=2)
ax.text(tips["tip_pct"].median() + .012, 22,
f"median {tips['tip_pct'].median():.1%}", color=ORANGE)
ax.set_xlabel("tip as a share of the bill")
ax.set_ylabel("number of checks")
ax.set_title("One tall cluster, and a long thin tail")
plt.show()
fig, axes = plt.subplots(1, 2, sharex=True, sharey=True, figsize=(9, 3.2))
for ax, (lab, col) in zip(axes, [("Lunch", BLUE), ("Dinner", ORANGE)]):
sub = tips.loc[tips["time"] == lab, "tip_pct"]
ax.hist(sub, bins=np.arange(0, 0.75, 0.025), color=col,
edgecolor="white", linewidth=0.5)
ax.axvline(sub.median(), color=GREY, lw=2, linestyle="--")
ax.set_title(f"{lab} (n = {len(sub)}, median {sub.median():.1%})", fontsize=10)
ax.set_xlabel("tip as a share of the bill")
axes[0].set_ylabel("number of checks")
fig.tight_layout()
plt.show()
4.3 Does the 15–20% norm hold?¶
print(f"within 15-20%: {tips['tip_pct'].between(0.15, 0.20).mean():.1%}")
print(f"below 15% : {(tips['tip_pct'] < 0.15).mean():.1%}")
print(f"above 20% : {(tips['tip_pct'] > 0.20).mean():.1%}")
within 15-20%: 38.9% below 15% : 45.3% above 20% : 15.8%
print("highest tip rates:")
print(tips.nlargest(4, "tip_pct")[["check_id", "total_bill", "tip", "tip_pct", "size"]])
print("\nlowest tip rates:")
print(tips.nsmallest(4, "tip_pct")[["check_id", "total_bill", "tip", "tip_pct", "size"]])
highest tip rates:
check_id total_bill tip tip_pct size
211 C0173 7.25 5.15 0.710 2
57 C0179 9.60 4.00 0.417 2
230 C0068 3.07 1.00 0.326 1
199 C0184 23.17 6.50 0.281 4
lowest tip rates:
check_id total_bill tip tip_pct size
141 C0238 32.83 1.17 0.036 2
106 C0103 44.30 2.50 0.056 3
23 C0058 26.41 1.50 0.057 2
210 C0001 16.99 1.01 0.059 2
Look at the total_bill column in each. The highest rates are all on small bills;
the lowest rates are all on large ones. A $5.15 tip on a $7.25 bill is 71%, but
in dollars it is unremarkable.
Two tails, one pattern — and it is a strong hint about Stage 5.
4.5 Do we drop the outliers?¶
No. And it is worth being explicit about why.
- They are not errors. A 71% tip is unusual, not impossible, and our validation
already confirmed
tip <= total_bill. - They are not noise — they are the clearest cases of the very pattern we are about to find.
- Dropping points because they are inconvenient is how analyses become dishonest.
What we do instead: report the median alongside the mean, since the median barely moves, and keep every row. Exercise 5 asks you to check that the conclusions survive their removal — which is the principled version of the same worry.
fig, ax = plt.subplots()
ax.scatter(tips["total_bill"], tips["tip"], s=22, color=BLUE,
alpha=0.65, edgecolor="white", linewidth=0.4)
ax.set_xlabel("total bill ($)"); ax.set_ylabel("tip ($)")
ax.set_title("Upward — and fanning out to the right")
plt.show()
fig, ax = plt.subplots()
ax.scatter(tips["total_bill"], tips["tip"], s=22, color=BLUE,
alpha=0.65, edgecolor="white", linewidth=0.4)
# --- the regression line
b1, b0 = np.polyfit(tips["total_bill"], tips["tip"], 1)
r = tips["total_bill"].corr(tips["tip"])
xs = np.linspace(tips["total_bill"].min(), tips["total_bill"].max(), 100)
ax.plot(xs, b0 + b1*xs, color=ORANGE, lw=2.2,
label=f"tip = {b0:.2f} + {b1:.3f} × bill ($r^2$ = {r**2:.2f})")
ax.legend(frameon=False, loc="upper left")
ax.set_xlabel("total bill ($)"); ax.set_ylabel("tip ($)")
ax.set_title("Tip rises with the bill — but from a positive base")
ax.axhline(b0, color=GREY, lw=1, linestyle=":")
ax.annotate(f"intercept ${b0:.2f}", xy=(tips["total_bill"].min(), b0),
xytext=(26, 0.9), color=GREY, fontsize=9,
arrowprops=dict(arrowstyle="->", color=GREY, lw=1))
plt.show()
5.2 Putting a number on it¶
tips[["total_bill", "tip", "size", "tip_pct"]].corr().round(3)
| total_bill | tip | size | tip_pct | |
|---|---|---|---|---|
| total_bill | 1.000 | 0.674 | 0.594 | -0.340 |
| tip | 0.674 | 1.000 | 0.488 | 0.341 |
| size | 0.594 | 0.488 | 1.000 | -0.141 |
| tip_pct | -0.340 | 0.341 | -0.141 | 1.000 |
tip and total_bill correlate at 0.68 — strong, and exactly what you would expect.
Now look at the last column. tip_pct correlates with total_bill at −0.34. Bigger
bills get bigger tips and smaller tip rates, at the same time. That is not a
contradiction, and the next cell explains why.
5.3 A line through the cloud¶
b1, b0 = np.polyfit(tips["total_bill"], tips["tip"], 1)
print(f"tip = {b0:.2f} + {b1:.3f} x bill")
tip = 0.92 + 0.105 x bill
Read it out loud: a tip starts at about 92 cents and adds about 10.5 cents per dollar of bill.
That positive intercept is the whole story. If tipping were a pure percentage, the line would pass through the origin.
5.4 The mechanism¶
If $\text{tip} = 0.92 + 0.105 \times \text{bill}$, then dividing through by the bill:
$$\text{tip\%} = \frac{0.92}{\text{bill}} + 0.105$$
The fixed 92 cents is a shrinking share of a growing bill.
bills = np.array([10, 20, 30, 40, 50])
pd.DataFrame({
"bill": bills,
"fixed part": (b0 / bills).round(3),
"rate part": b1.round(3),
"predicted tip %": (b0 / bills + b1).round(3),
})
| bill | fixed part | rate part | predicted tip % | |
|---|---|---|---|---|
| 0 | 10 | 0.092 | 0.105 | 0.197 |
| 1 | 20 | 0.046 | 0.105 | 0.151 |
| 2 | 30 | 0.031 | 0.105 | 0.135 |
| 3 | 40 | 0.023 | 0.105 | 0.128 |
| 4 | 50 | 0.018 | 0.105 | 0.123 |
The falling tip rate is not a change in generosity — it is arithmetic. That is a genuinely useful thing to be able to tell the owner, because it means the servers arguing about "cheap big tables" are arguing about a property of division.
5.5 Seeing it directly¶
cent = tips.groupby("bill_q", observed=True)["total_bill"].mean()
qm = tips.groupby("bill_q", observed=True)["tip_pct"].mean()
fig, ax = plt.subplots()
ax.scatter(tips["total_bill"], tips["tip_pct"], s=18, color=BLUE,
alpha=0.5, edgecolor="white", linewidth=0.4)
ax.plot(cent, qm, color=ORANGE, lw=2, marker="o", ms=8, label="quartile mean")
ax.legend(frameon=False)
ax.set_xlabel("total bill ($)"); ax.set_ylabel("tip rate")
ax.set_title("Tip rate falls as the bill grows")
plt.show()
tips.groupby("bill_q", observed=True)["tip_pct"].agg(["count", "mean", "median"]).round(4)
| count | mean | median | |
|---|---|---|---|
| bill_q | |||
| Q1 | 61 | 0.184 | 0.167 |
| Q2 | 56 | 0.166 | 0.155 |
| Q3 | 58 | 0.157 | 0.155 |
| Q4 | 59 | 0.133 | 0.139 |
tips.groupby("bill_q", observed=True)["tip"].agg(["mean", "std"]).round(2)
| mean | std | |
|---|---|---|
| bill_q | ||
| Q1 | 1.93 | 0.70 |
| Q2 | 2.63 | 0.72 |
| Q3 | 3.22 | 0.95 |
| Q4 | 4.29 | 1.70 |
The standard deviation of the dollar tip more than doubles from Q1 to Q4. Large bills are not just bigger — they are less predictable.
A statistician calls this heteroscedasticity. It is what makes the scatter fan out to the right, and it matters if you later fit a regression, because ordinary least squares assumes it is absent.
Stating Finding 2 for the owner¶
A server who turns a $20 table is likely to see about 15% of it. A server who turns a $40 table sees about 13% — but that is $5.20 instead of $3.00. The rate falls; the money rises.
This is why the choice of outcome in Stage 3 mattered so much.
tips.groupby("smoker")["tip_pct"].agg(["count", "mean", "median"]).round(4)
| count | mean | median | |
|---|---|---|---|
| smoker | |||
| No | 146 | 0.159 | 0.156 |
| Yes | 88 | 0.163 | 0.154 |
Smokers 16.3%, non-smokers 15.9%. A gap of 0.4 percentage points.
If you stopped here you would report "no real difference between smokers and non-smokers." That conclusion would be wrong — not because the means are wrong, but because the mean is the wrong summary.
6.2 Look at the spread¶
tips.groupby("smoker")["tip_pct"].agg(["mean", "std", "min", "max"]).round(4)
| mean | std | min | max | |
|---|---|---|---|---|
| smoker | ||||
| No | 0.159 | 0.039 | 0.057 | 0.266 |
| Yes | 0.163 | 0.087 | 0.036 | 0.710 |
The standard deviation among smoking-section parties is more than twice that of non-smokers. Same centre, completely different behaviour.
fig, ax = plt.subplots(figsize=(5.6, 3.2))
groups = [tips.loc[tips["smoker"] == v, "tip_pct"] for v in ["No", "Yes"]]
bp = ax.boxplot(groups, patch_artist=True, widths=0.5,
medianprops=dict(color="white", lw=2),
flierprops=dict(markersize=4, markeredgecolor=GREY))
for patch, c in zip(bp["boxes"], [BLUE, ORANGE]):
patch.set_facecolor(c); patch.set_edgecolor(c)
ax.set_xticks([1, 2]); ax.set_xticklabels(["No", "Yes"]) # not boxplot(labels=)
ax.set_xlabel("smoking section"); ax.set_ylabel("tip rate")
ax.set_title("Same centre, very different spread")
plt.show()
6.3 Quantifying the tails¶
tips.groupby("smoker")["tip_pct"].apply(
lambda s: pd.Series({"under 10%": (s < 0.10).mean(),
"over 25%": (s > 0.25).mean()})).unstack().round(3)
| under 10% | over 25% | |
|---|---|---|
| smoker | ||
| No | 0.062 | 0.014 |
| Yes | 0.193 | 0.080 |
Finding 3¶
Smoking-section parties are three times as likely to tip under 10% (19.3% vs 6.2%) and nearly six times as likely to tip over 25%. The averages match; the behaviour does not.
The lesson: compare distributions, not just centres. Two groups with identical means can be completely different, and a table of means will never tell you.
6.4 Sex: a raw difference¶
tips.groupby("sex")[["total_bill", "tip_pct"]].mean().round(4)
| total_bill | tip_pct | |
|---|---|---|
| sex | ||
| Female | 18.300 | 0.167 |
| Male | 20.911 | 0.157 |
Women tip 16.7%, men 15.7% — a full percentage point. A tempting headline.
But look at the first column before you write it. Women in this sample pay smaller bills ($18.30 vs $20.91). And we know from Finding 2 that smaller bills carry higher tip rates for reasons that have nothing to do with who is paying.
Two competing explanations:
sex ─────────────────────────────> tip_pct (a real sex effect)
sex ──> total_bill ──> tip_pct (or just Finding 2 again)
The fix: compare men and women at the same bill size. That is what stratification means.
6.5 Stratifying¶
tips.pivot_table(index="bill_q", columns="sex", values="tip_pct",
aggfunc="mean", observed=True).round(4)
| sex | Female | Male |
|---|---|---|
| bill_q | ||
| Q1 | 0.187 | 0.182 |
| Q2 | 0.175 | 0.159 |
| Q3 | 0.152 | 0.158 |
| Q4 | 0.137 | 0.131 |
piv = tips.pivot_table(index="bill_q", columns="sex", values="tip_pct",
aggfunc="mean", observed=True)
cnt = tips.pivot_table(index="bill_q", columns="sex", values="tip_pct",
aggfunc="size", observed=True)
err = tips.pivot_table(index="bill_q", columns="sex", values="tip_pct",
aggfunc="sem", observed=True)
xs = np.arange(len(piv)); w = 0.38
fig, ax = plt.subplots()
for i, (lab, col) in enumerate([("Female", BLUE), ("Male", ORANGE)]):
off = (i - 0.5) * w
ax.bar(xs + off, piv[lab], w, color=col, label=lab,
yerr=err[lab], capsize=3, error_kw=dict(ecolor=GREY, lw=1))
for x, v, n in zip(xs + off, piv[lab], cnt[lab]):
ax.text(x, 0.008, f"n={n}", ha="center", va="bottom",
fontsize=7.5, color="white")
ax.set_xticks(xs); ax.set_xticklabels(piv.index)
ax.set_ylim(0, 0.235)
ax.yaxis.set_major_formatter(lambda v, p: f"{v:.0%}")
ax.legend(frameon=False, ncol=2, loc="upper right")
ax.set_xlabel("bill quartile"); ax.set_ylabel("mean tip rate")
ax.set_title("Within every bill quartile, the sex gap disappears")
ax.set_xticklabels(
[f"{q}\n(${tips.loc[tips['bill_q']==q,'total_bill'].min():.0f}–"
f"${tips.loc[tips['bill_q']==q,'total_bill'].max():.0f})" for q in piv.index])
plt.show()
Finding 4¶
Within every bill quartile the gap collapses to a fraction of a point, and it changes direction twice. The apparent sex difference was a bill-size difference wearing a disguise.
The general lesson¶
A difference between two groups is evidence about those groups only if everything else is comparable. Before you believe any group difference:
| Ask | pandas |
|---|---|
| Are the groups the same size? | df.groupby(g).size() |
| Do they differ on the main driver? | df.groupby(g)[driver].mean() |
| Does the gap survive stratification? | df.pivot_table(index=bin, columns=g) |
| Is it spread rather than centre? | df.groupby(g)[y].std() |
6.6 Day and time: check the cross-tab FIRST¶
pd.crosstab(tips["day"], tips["time"], margins=True)
| time | Dinner | Lunch | All |
|---|---|---|---|
| day | |||
| Fri | 11 | 7 | 18 |
| Sat | 82 | 0 | 82 |
| Sun | 73 | 0 | 73 |
| Thur | 1 | 60 | 61 |
| All | 167 | 67 | 234 |
Stop. Saturday and Sunday are dinner-only. Thursday is lunch except for a single check. There is almost no overlap between the two variables.
print(tips.groupby("time")["tip_pct"].mean().round(4))
print()
print(tips.groupby("day")["tip_pct"].mean().round(4))
time Dinner 0.159 Lunch 0.164 Name: tip_pct, dtype: float64 day Fri 0.167 Sat 0.152 Sun 0.167 Thur 0.161 Name: tip_pct, dtype: float64
In this dataset, "lunch tips better than dinner" and "Thursday tips better than Saturday" are the same sentence. The design does not let you separate the day from the meal.
No amount of cleverness fixes this — it is a property of the data collection, not of your method. The honest move is to say so in the report.
This is exactly why crosstab belongs before groupby, not after.
6.7 What we can still say about days¶
tips.groupby("day").agg(n=("tip", "size"),
bill=("total_bill", "mean"),
pct=("tip_pct", "mean")).round(3)
| n | bill | pct | |
|---|---|---|---|
| day | |||
| Fri | 18 | 17.474 | 0.167 |
| Sat | 82 | 20.691 | 0.152 |
| Sun | 73 | 21.648 | 0.167 |
| Thur | 61 | 17.786 | 0.161 |
Saturday is lowest at 15.2% — and Saturday also has the second-highest average bill, so Finding 2 may be all we are looking at. With only 18 checks on Friday, the day-to-day differences are small relative to the noise.
6.8 One table for the owner¶
tips.pivot_table(index="time", columns="smoker", values="tip_pct",
aggfunc=["mean", "count"], margins=True).round(3)
| mean | count | |||||
|---|---|---|---|---|---|---|
| smoker | No | Yes | All | No | Yes | All |
| time | ||||||
| Dinner | 0.158 | 0.160 | 0.159 | 102 | 65 | 167 |
| Lunch | 0.161 | 0.170 | 0.164 | 44 | 23 | 67 |
| All | 0.158 | 0.163 | 0.160 | 146 | 88 | 234 |
6.9 transform: removing the bill effect in one line¶
Stratifying with pivot_table gives a grid you read by eye. transform does the same
adjustment arithmetically, so you can keep analysing the adjusted values.
tips["pct_vs_peers"] = (tips["tip_pct"]
- tips.groupby("bill_q", observed=True)["tip_pct"].transform("mean"))
print("mean deviation from same-quartile peers:")
print(tips.groupby("sex")["pct_vs_peers"].mean().round(4))
print()
print(tips.groupby("smoker")["pct_vs_peers"].mean().round(4))
mean deviation from same-quartile peers: sex Female 0.004 Male -0.002 Name: pct_vs_peers, dtype: float64 smoker No -0.003 Yes 0.004 Name: pct_vs_peers, dtype: float64
transform broadcasts each quartile's mean back onto every row, so subtracting it strips
the bill-size effect out of every check at once.
The residual sex gap is under a quarter of a percentage point — the same conclusion as the stratified table, in one line instead of a grid.
Stage 7 — One Figure That Tells the Story¶
fig, axes = plt.subplots(2, 2, figsize=(10, 6))
ax = axes[0, 0]
ax.hist(tips["tip_pct"], bins=25, color=BLUE, edgecolor="white", linewidth=.4)
ax.axvline(tips["tip_pct"].median(), color=ORANGE, lw=1.8)
ax.set_title("1. Median tip 15.4%, long right tail")
ax.set_xlabel("tip rate")
ax = axes[0, 1]
ax.scatter(tips["total_bill"], tips["tip_pct"], s=12, color=BLUE, alpha=.45)
ax.plot(cent, qm, color=ORANGE, lw=1.8, marker="o", ms=5)
ax.set_ylim(0, .45)
ax.set_title("2. Rate falls as the bill grows")
ax.set_xlabel("bill ($)")
ax = axes[1, 0]
g = [tips.loc[tips["smoker"] == v, "tip_pct"] for v in ["No", "Yes"]]
bp = ax.boxplot(g, patch_artist=True, widths=.45, showfliers=False,
medianprops=dict(color="white", lw=1.6))
for p, c in zip(bp["boxes"], [BLUE, ORANGE]):
p.set_facecolor(c); p.set_edgecolor(c)
ax.set_xticks([1, 2]); ax.set_xticklabels(["non-smoking", "smoking"])
ax.set_title("3. Same centre, double the spread")
ax = axes[1, 1]
piv = tips.pivot_table(index="bill_q", columns="sex", values="tip_pct",
aggfunc="mean", observed=True)
xs = np.arange(len(piv))
ax.bar(xs - .19, piv["Female"], .36, color=BLUE, label="Female")
ax.bar(xs + .19, piv["Male"], .36, color=ORANGE, label="Male")
ax.set_xticks(xs); ax.set_xticklabels(piv.index)
ax.set_ylim(0, .285); ax.legend(frameon=False, ncol=2, loc="upper center")
ax.set_title("4. No sex gap within bill quartiles")
fig.tight_layout()
fig.savefig("tips_summary.png", dpi=200, bbox_inches="tight")
plt.show()
Figure rules that matter:
- One message per panel, and the message goes in the panel title — not in a caption the reader has to hunt for.
- A legend whenever there are two series, and never colour alone to carry meaning.
- Panel 3 hides the outlying points (
showfliers=False) so the boxes stay readable — which is fine because we showed them in full earlier. Hiding data silently is not. - These two colours were checked for colour-vision-deficiency separation. Roughly 1 in 12 men cannot distinguish a red/green pair.
Stage 8 — Report¶
The findings¶
1. The norm is a myth here. The median check tips 15.4%, but only 39% of parties fall in the 15–20% band and 45% tip less than 15%.
2. The tip rate falls with the bill, from 18.5% on the smallest quarter of checks to 13.3% on the largest. The cause is arithmetic, not generosity: tips behave like $0.92 + 10.5% of the bill, and the fixed part shrinks as a share. In dollars, large tables are still worth far more.
3. Smoking-section parties are not worse tippers — they are more variable ones. The means differ by 0.4 points, but they are three times as likely to tip under 10% and six times as likely to tip over 25%.
4. There is no sex difference. The raw one-point gap disappears entirely once men and women are compared at the same bill size.
What we cannot say¶
An analysis is judged as much by its stated limits as by its findings.
- Day and meal cannot be separated. Saturday and Sunday are dinner-only; Thursday is lunch-only. Any "weekend effect" here is also a "dinner effect".
- One restaurant, one server, about ten weeks. Nothing here generalises to restaurants in general.
- 234 of 247 checks. We dropped 10 incomplete records and 3 duplicates. The dropped rows were spread across days, but we cannot prove they were missing at random.
- Nothing here is causal. We observed what happened; we changed nothing.
- Small cells. Friday has 18 checks. Treat its average as a hint, not a number.
The whole analysis, start to finish¶
tips_final = (load_tips("tips_messy.csv")
.assign(tip_pct = lambda d: d["tip"] / d["total_bill"],
bill_pp = lambda d: d["total_bill"] / d["size"],
bill_q = lambda d: pd.qcut(d["total_bill"], 4,
labels=["Q1", "Q2", "Q3", "Q4"])))
headline = tips_final.groupby("bill_q", observed=True)["tip_pct"].mean()
by_smoker = tips_final.groupby("smoker")["tip_pct"].agg(["mean", "std"])
by_sex = tips_final.pivot_table(index="bill_q", columns="sex",
values="tip_pct", observed=True)
print(headline.round(4), "\n")
print(by_smoker.round(4), "\n")
print(by_sex.round(4))
bill_q
Q1 0.184
Q2 0.166
Q3 0.157
Q4 0.133
Name: tip_pct, dtype: float64
mean std
smoker
No 0.159 0.039
Yes 0.163 0.087
sex Female Male
bill_q
Q1 0.187 0.182
Q2 0.175 0.159
Q3 0.152 0.158
Q4 0.137 0.131
Twelve lines, from a broken CSV to every number in the report. That is the goal: an analysis someone else can re-run and get your answer.
A checklist you can reuse¶
| Stage | Do not skip |
|---|---|
| Load | shape, dtypes, isna().sum(), duplicated().sum(), value_counts() on every category |
| Clean | work on a copy; write the defect inventory first |
count NaNs after every map and every type conversion |
|
assert your expectations; check who you are dropping |
|
| Derive | choose the outcome deliberately; bin the main driver |
| Explore | median beside mean; look before you summarise |
| Compare | crosstab before groupby; stratify every raw gap; spread, not just centre |
| Report | state the limits as plainly as the findings |
References¶
- Bryant, P.G. & Smith, M.A. (1995) Practical Data Analysis: Case Studies in Business
Statistics. Irwin — the original source of the
tipsdata. - McKinney, W. Python for Data Analysis, 3rd ed. — Data cleaning and preparation · Data aggregation and group operations · Plotting and visualization
- Wickham, H. (2014) "Tidy Data", Journal of Statistical Software 59(10). doi:10.18637/jss.v059.i10
- Irizarry, R.A. Introduction to Data Science — on confounding and stratification.
- pandas user guide