Capstone: from a messy dataset to a reported result

Every lesson so far handed you a tidy dataset. Real data does not arrive that way. It comes with a typo where an age should be, a blank cell, a column read as text when you need a number, and a clinician who wants one answer, not a methods seminar. This capstone takes a single messy dataset from the raw rows to one sentence you could put in a paper. We use a dengue cohort, because dengue fills Malaysian wards every monsoon season and the analysis touches most of the course.

The clinical question is narrow on purpose. Among adults admitted with dengue, does a secondary infection, that is a repeat infection with a different serotype, raise the odds of progressing to severe dengue? And is any association we see just an effect of age? You met every piece of this on its own already. Here you connect them into one pipeline.

The shape of a real analysis

Before any code, hold the whole arc in your head. A complete analysis runs in close to the same order every time.

1. Import the data and get it into a data frame.
2. Inspect every column: types, ranges, and how many values are missing.
3. Clean the problems you found, and write down each decision.
4. Describe the sample with a baseline table before any model.
5. Choose the method from the study design and the variable types.
6. Run the model, crude first, then adjusted for the confounder.
7. Check the fit with a quick sanity test against the observed rates.
8. Interpret the estimate in plain clinical language, with its interval.
9. Report it as a methods-and-results sentence, and keep the script.

The order earns its keep. Cleaning before you describe, describing before you model, and checking before you interpret each catch a mistake the next step would otherwise bury.

The dataset

Picture an export from a ward register: 65 adults admitted with confirmed dengue over one season, three columns each. age in years, secondary recording whether serology showed a secondary infection, and severe flagging whether the patient met the criteria for severe dengue (1) or not (0). The export is typical of a hand-kept register. Someone mistyped an age, two ages are blank, and the infection column arrived as text rather than a number.

Import and inspect

The first move is never to model. It is to look. Build the data frame and ask what each column actually holds, because the answer decides everything downstream.

age <- c(43, 51, 67, 52, 68, 69, 44, 57, 62, 53, 58, 45, 52, 30, 51, 49, 41, 26, 39, 39, 50, 38, 60, 29, 24, 24, 30, 65, 48, 65, 57, 27, 42, 46, 58, 43, 33, 31, 20, 47, 57, 57, 50, 41, 69, 19, 31, 57, 51, 60, 19, 42, 55, 33, 18, 20, 22, 70, 28, 49, 57, 43, 999, NA, NA)
secondary <- c("yes", "no", "yes", "yes", "yes", "yes", "yes", "no", "no", "yes", "yes", "yes", "yes", "yes", "yes", "yes", "no", "no", "no", "yes", "no", "no", "yes", "no", "no", "yes", "yes", "yes", "no", "no", "yes", "no", "yes", "yes", "no", "yes", "no", "no", "no", "no", "yes", "no", "yes", "yes", "yes", "no", "no", "no", "yes", "no", "no", "yes", "yes", "yes", "no", "no", "no", "yes", "yes", "yes", "yes", "no", "yes", "no", "yes")
severe <- c(0, 1, 1, 1, 1, 0, 1, 0, 0, 1, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1)
dengue <- data.frame(age, secondary, severe, stringsAsFactors = FALSE)

# look before you model: what is actually in each column?

str(dengue)
summary(dengue$age)        # maximum is 999, not a human age
sum(is.na(dengue$age))     # 2 ages are missing
table(dengue$secondary)    # text "yes"/"no", not the 0/1 a model needs

import pandas as pd
import numpy as np
age = [43, 51, 67, 52, 68, 69, 44, 57, 62, 53, 58, 45, 52, 30, 51, 49, 41, 26, 39, 39, 50, 38, 60, 29, 24, 24, 30, 65, 48, 65, 57, 27, 42, 46, 58, 43, 33, 31, 20, 47, 57, 57, 50, 41, 69, 19, 31, 57, 51, 60, 19, 42, 55, 33, 18, 20, 22, 70, 28, 49, 57, 43, 999, np.nan, np.nan]
secondary = ["yes", "no", "yes", "yes", "yes", "yes", "yes", "no", "no", "yes", "yes", "yes", "yes", "yes", "yes", "yes", "no", "no", "no", "yes", "no", "no", "yes", "no", "no", "yes", "yes", "yes", "no", "no", "yes", "no", "yes", "yes", "no", "yes", "no", "no", "no", "no", "yes", "no", "yes", "yes", "yes", "no", "no", "no", "yes", "no", "no", "yes", "yes", "yes", "no", "no", "no", "yes", "yes", "yes", "yes", "no", "yes", "no", "yes"]
severe = [0, 1, 1, 1, 1, 0, 1, 0, 0, 1, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1]
dengue = pd.DataFrame({"age": age, "secondary": secondary, "severe": severe})

# look before you model: what is actually in each column?

print(dengue.dtypes)
print(dengue["age"].describe())       # max 999, impossible
print(dengue["age"].isna().sum())     # 2 missing
print(dengue["secondary"].value_counts())   # text, not 0/1

Three problems surface at once. The age summary reports a maximum of 999, which is no human age. Two ages are missing. And secondary is a character column of "yes" and "no", not the 0/1 a model reads. None of these is exotic. Almost every real dataset hides a few like them, and inspection is how you find them before they corrupt a result.

Make the cleaning decisions, and document them

Cleaning is a set of decisions, and each one needs a reason you could defend to a reviewer. Write them down as you make them. Three decisions cover this dataset.

Clean the data, then describe it

Apply the three decisions in code, then build a baseline table from the cleaned data. The baseline table, called Table 1 in most papers, summarises the sample by exposure group before any model runs.

age <- c(43, 51, 67, 52, 68, 69, 44, 57, 62, 53, 58, 45, 52, 30, 51, 49, 41, 26, 39, 39, 50, 38, 60, 29, 24, 24, 30, 65, 48, 65, 57, 27, 42, 46, 58, 43, 33, 31, 20, 47, 57, 57, 50, 41, 69, 19, 31, 57, 51, 60, 19, 42, 55, 33, 18, 20, 22, 70, 28, 49, 57, 43, 999, NA, NA)
secondary <- c("yes", "no", "yes", "yes", "yes", "yes", "yes", "no", "no", "yes", "yes", "yes", "yes", "yes", "yes", "yes", "no", "no", "no", "yes", "no", "no", "yes", "no", "no", "yes", "yes", "yes", "no", "no", "yes", "no", "yes", "yes", "no", "yes", "no", "no", "no", "no", "yes", "no", "yes", "yes", "yes", "no", "no", "no", "yes", "no", "no", "yes", "yes", "yes", "no", "no", "no", "yes", "yes", "yes", "yes", "no", "yes", "no", "yes")
severe <- c(0, 1, 1, 1, 1, 0, 1, 0, 0, 1, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1)
dengue <- data.frame(age, secondary, severe, stringsAsFactors = FALSE)
dengue$age[dengue$age > 120] <- NA                    # impossible value -> missing
dengue$sec <- ifelse(dengue$secondary == "yes", 1, 0)   # text column -> 0/1
clean <- dengue[!is.na(dengue$age), ]                  # keep complete ages

# clean first, then build Table 1 from the clean data

nrow(clean)                                  # 62 of 65 rows kept
table(clean$secondary, clean$severe)         # the 2x2 of exposure by outcome
tapply(clean$age, clean$secondary, mean)     # mean age 39.0 primary , 49.7 secondary
tapply(clean$age, clean$severe, mean)        # mean age 38.0 non-severe , 50.8 severe
tapply(clean$severe, clean$secondary, mean)  # severe rate 0.32 primary , 0.71 secondary

import pandas as pd
import numpy as np
age = [43, 51, 67, 52, 68, 69, 44, 57, 62, 53, 58, 45, 52, 30, 51, 49, 41, 26, 39, 39, 50, 38, 60, 29, 24, 24, 30, 65, 48, 65, 57, 27, 42, 46, 58, 43, 33, 31, 20, 47, 57, 57, 50, 41, 69, 19, 31, 57, 51, 60, 19, 42, 55, 33, 18, 20, 22, 70, 28, 49, 57, 43, 999, np.nan, np.nan]
secondary = ["yes", "no", "yes", "yes", "yes", "yes", "yes", "no", "no", "yes", "yes", "yes", "yes", "yes", "yes", "yes", "no", "no", "no", "yes", "no", "no", "yes", "no", "no", "yes", "yes", "yes", "no", "no", "yes", "no", "yes", "yes", "no", "yes", "no", "no", "no", "no", "yes", "no", "yes", "yes", "yes", "no", "no", "no", "yes", "no", "no", "yes", "yes", "yes", "no", "no", "no", "yes", "yes", "yes", "yes", "no", "yes", "no", "yes"]
severe = [0, 1, 1, 1, 1, 0, 1, 0, 0, 1, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1]
dengue = pd.DataFrame({"age": age, "secondary": secondary, "severe": severe})
dengue.loc[dengue["age"] > 120, "age"] = np.nan            # impossible value -> NA
dengue["sec"] = (dengue["secondary"] == "yes").astype(int)  # text column -> 0/1
clean = dengue.dropna(subset=["age"]).copy()              # keep complete ages

# clean first, then build Table 1 from the clean data

print(len(clean))                                       # 62 of 65 rows kept
print(pd.crosstab(clean["secondary"], clean["severe"]))   # the 2x2
print(clean.groupby("secondary")["age"].mean().round(1))  # 39.0 , 49.7
print(clean.groupby("severe")["age"].mean().round(1))     # 38.0 , 50.8
print(clean.groupby("secondary")["severe"].mean().round(2))  # 0.32 , 0.71

Sixty-two of the 65 patients remain after the three bad rows drop out. Read the baseline table before you fit anything.

Two patterns jump out. Patients with a secondary infection are older on average (49.7 against 39.0 years) and far more often severe (71% against 32%). Look at the outcome the other way and the same story holds: severe cases are older than non-severe cases (50.8 against 38.0 years). So age is tied to the exposure and to the outcome at the same time. That is the fingerprint of a confounder, and the baseline table flagged it before a single model ran.

Choose the method from the design

The method falls out of two facts you already hold. The outcome, severe dengue, is binary. The design is a cohort followed to that outcome, and you want to adjust for age. A binary outcome with one or more variables to adjust for points to logistic regression, the choice you made in the method-selection lesson (16.1). Logistic regression models the log odds of the outcome and returns an odds ratio for each term, holding the others fixed. Put age in the model and the secondary-infection effect comes back adjusted for age in a single fit. This is the regression form of the stratify-and-adjust idea from the confounding chapter, and it reads exactly like the model you built in the logistic-regression chapter.

Run the model: crude, then adjusted

Fit two models. The crude model uses the exposure alone, the way a raw 2x2 table would. The adjusted model adds age. Running both, and comparing them, is how you watch confounding move the estimate.

age    <- c(43, 51, 67, 52, 68, 69, 44, 57, 62, 53, 58, 45, 52, 30, 51, 49, 41, 26, 39, 39, 50, 38, 60, 29, 24, 24, 30, 65, 48, 65, 57, 27, 42, 46, 58, 43, 33, 31, 20, 47, 57, 57, 50, 41, 69, 19, 31, 57, 51, 60, 19, 42, 55, 33, 18, 20, 22, 70, 28, 49, 57, 43)
sec    <- c(1, 0, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 1, 1, 0, 0, 1, 0, 1, 1, 0, 1, 0, 0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 1, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 1, 1, 0)
severe <- c(0, 1, 1, 1, 1, 0, 1, 0, 0, 1, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 1, 1, 1, 1, 1, 1)

# crude uses exposure alone; adjusted adds age

crude <- glm(severe ~ sec, family = binomial)
adj   <- glm(severe ~ sec + age, family = binomial)
exp(coef(crude))["sec"]                 # crude OR 5.07
exp(confint.default(crude))["sec", ]    # 95% CI 1.72 to 14.97
exp(coef(adj))["sec"]                   # adjusted OR 3.41
exp(confint.default(adj))["sec", ]      # 95% CI 1.06 to 10.95
exp(coef(adj))["age"]                   # 1.06 per year of age
# quick check: model-predicted group risk should match the observed rate
tapply(predict(adj, type = "response"), sec, mean)  # 0.32 , 0.71
tapply(severe, sec, mean)                            # 0.32 , 0.71

import pandas as pd
import numpy as np
import statsmodels.formula.api as smf
age = [43, 51, 67, 52, 68, 69, 44, 57, 62, 53, 58, 45, 52, 30, 51, 49, 41, 26, 39, 39, 50, 38, 60, 29, 24, 24, 30, 65, 48, 65, 57, 27, 42, 46, 58, 43, 33, 31, 20, 47, 57, 57, 50, 41, 69, 19, 31, 57, 51, 60, 19, 42, 55, 33, 18, 20, 22, 70, 28, 49, 57, 43]
sec = [1, 0, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 1, 1, 0, 0, 1, 0, 1, 1, 0, 1, 0, 0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 1, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 1, 1, 0]
severe = [0, 1, 1, 1, 1, 0, 1, 0, 0, 1, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 1, 1, 1, 1, 1, 1]
df = pd.DataFrame({"age": age, "sec": sec, "severe": severe})

# crude uses exposure alone; adjusted adds age

crude = smf.logit("severe ~ sec", data=df).fit(disp=0)
adj   = smf.logit("severe ~ sec + age", data=df).fit(disp=0)
print(np.exp(crude.params["sec"]))                   # crude OR 5.07
print(np.exp(crude.conf_int().loc["sec"]).round(2))  # 1.72 , 14.97
print(np.exp(adj.params["sec"]))                     # adjusted OR 3.41
print(np.exp(adj.conf_int().loc["sec"]).round(2))    # 1.06 , 10.95
print(np.exp(adj.params["age"]))                     # 1.06 per year
df["fit"] = adj.predict(df)
print(df.groupby("sec")[["severe", "fit"]].mean().round(2))  # 0.32/0.32 , 0.71/0.71

The crude odds ratio is 5.07. On its own it says a secondary infection multiplies the odds of severe dengue about fivefold. The adjusted odds ratio, with age in the model, is 3.41. The estimate fell because age was doing part of the work the crude model had credited to infection.

A quick check before you trust it

Never report a model you have not sanity-checked. The fastest check compares the risk the model predicts in each group with the rate you actually saw. Average the fitted probabilities within each infection group and you get 0.32 for primary and 0.71 for secondary. The observed rates are 9/28 = 0.32 and 24/34 = 0.71. They match, which is what you expect when the grouping variable is in the model. A mismatch here would point straight to a coding error, such as the outcome reversed, and you would fix it before reading a single odds ratio.

Interpret the adjusted odds ratio

Put the adjusted result in clinical language. After accounting for age, adults with a secondary dengue infection had about 3.4 times the odds of severe disease compared with a primary infection (95% CI 1.06 to 10.95). The interval stays above 1, so the cohort is consistent with a real independent effect, though with 62 patients it is wide and its lower edge sits close to 1. The honest reading is a likely effect that this sample is too small to pin down tightly.

The gap between 5.07 and 3.41 is confounding by age, made into a number. Older patients were both more likely to carry a secondary infection and more likely to develop severe dengue. The age odds ratio is 1.06 per year, about 1.8 per decade. That shared link inflated the crude estimate, because the crude model had no way to tell the infection effect apart from the age effect riding along with it. Adjusting for age separated them and left the part that belongs to infection itself. Report both numbers. The size of the move from 5.07 to 3.41 is itself the evidence of how much age mattered.

Write the result the way a paper would

A reader does not want your code. They want a sentence or two that state what you did and what you found, with the numbers attached and the exclusions declared.

Make it reproducible

The analysis above is one script that runs from the raw rows to 3.41 without a manual step in between. That is the point of scripting it. Anyone with the script and the data can rerun it and land on the same answer, including you after you have forgotten the details. A few habits keep it that way.

• Script every step. Do every edit in code, never by hand in a spreadsheet. The 999 became NA on a line you can read, not in a cell you cannot.
• Set the seed before anything random. If a step bootstraps or simulates, call set.seed() in R or np.random.default_rng(seed) in Python first, so the run repeats exactly.
• Comment each decision next to the line that makes it. The reason a row dropped should live beside the code that drops it.
• Keep one dataset and one set of names. Carrying severe and sec through the whole script means a reviewer never has to relearn your coding.

Common mistakes

• Modelling before inspecting. Fit the model on the raw frame and the age of 999 quietly pulls the age coefficient toward nonsense. The bad value has to go before, not after, you run anything.
• Mean-imputing an impossible value. Replacing 999 with the mean age treats a known error as if it were data you measured. Recode it to missing instead, and let the analysis drop it.
• Reporting only the adjusted estimate. The crude and the adjusted belong together. Without the crude 5.07, a reader cannot see that age moved the result, which is half the finding.
• Reading exp(coef) as a risk ratio. It is an odds ratio. With a common outcome like severe dengue here, the odds ratio sits further from 1 than the risk ratio, so calling 3.41 a risk ratio overstates the effect.
• Skipping the fit check. If the predicted group risks had not matched 0.32 and 0.71, the odds ratio would have been meaningless, and only the check would have caught it.

Tips

• Run a one-line range or summary on every numeric column before you trust it. Impossible maximums and negative ages surface in seconds.
• Build Table 1 first. Group imbalance in the baseline table tells you which confounders to adjust for before you fit anything.
• Decide your confounders from clinical knowledge, here that age drives both dengue exposure history and severity, not by hunting for whichever variable shifts the estimate.
• Quote the crude and adjusted odds ratios side by side. The distance between them is the clearest evidence of confounding you can show a reader.
• Exponentiate coefficients and their confidence limits, never the other way round. The interval is built on the log-odds scale, then antilogged.

You can now run an analysis end to end

That is the whole course in one pipeline. You imported a messy file, inspected it, found and fixed an impossible value, a pair of missing ages, and a mistyped column. You described the sample, chose logistic regression from the design, ran the crude and adjusted models, checked the fit, read the odds ratio in clinical language, and wrote the sentence a paper would print. You did all of it in both R and Python on the same data, and you kept the script so the next person can reproduce it. The statistics you learned along the way, proportions and intervals, the chi-squared and t tests, confounding and stratification, and regression, are the tools. Knowing which one a real question needs, and carrying it cleanly from raw data to a defensible result, is the work. You can do that now. Well done.