# Exploratory Data Analysis Report: [Your data set or research question] > A report template for exploratory data analysis (EDA): the first, open-minded look at a data set before any formal > test. It walks you through purpose, data overview, univariate summaries (centre, spread and shape), frequency tables, > bivariate summaries, three charts, an outlier screen and a findings section that separates what you noticed from what > you will later test. It is for students learning descriptive statistics, master's and PhD students meeting their > thesis data for the first time, and analysts who need a quick, honest profile of a new table. It speeds up learning > by fixing the order of the questions, showing one variable fully worked, one half worked and one left for you, and by > ending with self-test questions whose answers sit in a separate key. The worked examples use fictional data; replace > them with your own. > Source basis: Tukey (1977) on exploratory data analysis, the NIST/SEMATECH e-Handbook of Statistical Methods > (chapter 1, EDA), the GAISE College Report (ASA 2016), and the R, pandas and SciPy documentation, all written up in > original words and examples. ## Document control | Field | Value | |---|---| | Course or project | [STAT 101 assignment / thesis chapter 3 / client data review] | | Learner or analyst | [Your name] | | Data set | [Lakemoor University study habits survey, 20 students (fictional) / your data set, file name and source] | | Data version or extract date | [v1, extracted YYYY-MM-DD] | | Software and versions | R 4.5.1 (base and stats); Python 3.14 with numpy 2.5.4, pandas 3.1.0, scipy 1.18.1 | | Date | YYYY-MM-DD | | Version | [0.1] | ## 1. Purpose and questions **Purpose.** [Example: Before testing whether study time predicts exam performance in the Lakemoor University STAT 101 cohort, describe the data, check its quality, and find features (skew, outliers, group differences) that the later analysis plan must handle.] EDA is the stage where you look at the data with as few assumptions as possible. John Tukey, who named the field, compared it to detective work: you collect clues and form hypotheses. Confirmatory analysis is the courtroom: a hypothesis written down in advance is tested with a method chosen in advance. Both are needed, but they must not be mixed. If you find a pattern while exploring and then test it on the same data, the p-value is too optimistic, because you chose the test after seeing the pattern. | # | Question | Why it matters | Answered in section | |---|---|---|---| | *Example:* Q1 | What is a typical weekly study time, and how variable is it? | Sets expectations for the main analysis | 3 | | *Example:* Q2 | Are there unusual values or data entry errors? | Outliers can dominate means and correlations | 7 | | *Example:* Q3 | Does exam score rise with study hours, and is the pattern linear? | Decides between correlation, regression or a transformation | 6 | | *Example:* Q4 | Do the year groups differ in study time? | Possible confounder for Q3 | 6 | | [Q5] | [Your question] | [Why] | [ ] | ## 2. Data overview **Source and unit.** [Example: an anonymous survey of 20 STAT 101 students at Lakemoor University (fictional), one row per student, collected in week 8; exam score from the midterm.] **Sample size:** n = 20. **Missing values:** none. **Known quality issues:** [none / list]. | Variable | Description | Measurement level | Values or unit | Missing | |---|---|---|---|---| | *Example:* id | Student identifier | Nominal (key) | 1 to 20 | 0 | | *Example:* year | Year of study | Ordinal | 1, 2, 3 | 0 | | *Example:* mode | Usual way of getting to campus | Nominal | Bike, Bus, Car, Walk | 0 | | *Example:* hours | Self-reported study time per week | Ratio | hours, to 0.5 h | 0 | | *Example:* score | Midterm exam score | Interval (treated as) | 0 to 100 points | 0 | | [variable] | [description] | [level] | [values] | [n] | The measurement level decides the summaries you may use. For nominal variables you can count and give proportions and the mode. For ordinal variables you can also give the median and percentiles. For interval and ratio variables you can use means and standard deviations, and for ratio variables only, statements like "twice as much" make sense. The example data (fictional): | id | year | mode | hours | score | | id | year | mode | hours | score | |---|---|---|---|---|---|---|---|---|---|---| | 1 | 1 | Bus | 3.5 | 50 | | 11 | 2 | Car | 8.5 | 67 | | 2 | 1 | Walk | 6.0 | 58 | | 12 | 2 | Bike | 9.0 | 65 | | 3 | 1 | Bus | 7.5 | 66 | | 13 | 2 | Walk | 11.0 | 75 | | 4 | 1 | Car | 8.0 | 61 | | 14 | 2 | Bus | 4.5 | 55 | | 5 | 1 | Walk | 5.0 | 57 | | 15 | 3 | Car | 13.0 | 74 | | 6 | 1 | Bus | 9.0 | 70 | | 16 | 3 | Walk | 8.0 | 69 | | 7 | 1 | Bike | 10.0 | 68 | | 17 | 3 | Bus | 7.0 | 62 | | 8 | 2 | Walk | 12.5 | 78 | | 18 | 3 | Bike | 6.0 | 60 | | 9 | 2 | Bus | 6.5 | 64 | | 19 | 3 | Walk | 10.0 | 73 | | 10 | 2 | Bus | 7.0 | 59 | | 20 | 3 | Bus | 28.0 | 81 | ## 3. How EDA works: the summaries and what they mean Every numeric variable has three features worth describing: its **centre** (a typical value), its **spread** (how far values sit from the centre) and its **shape** (symmetric or skewed, one peak or several, heavy or light tails). On top of that you look for **unusual values**. There are two families of summaries. The **moment-based** family (mean, standard deviation, skewness, kurtosis) uses every value and is efficient for symmetric data, but a single extreme value can move it a lot. The **order-based** family (median, quartiles, interquartile range, five-number summary) uses the sorted positions and is **resistant**: moving the largest value from 28 to 280 changes none of them. A good EDA reports both and asks why they disagree when they do. The mean and the sample standard deviation: $$ \bar{x} = \frac{1}{n}\sum_{i=1}^{n} x_i, \qquad s = \sqrt{\frac{1}{n-1}\sum_{i=1}^{n}\left(x_i - \bar{x}\right)^2} $$ The divisor $n-1$ makes $s^2$ an unbiased estimate of the population variance. R's `sd()` and pandas' `.std()` use $n-1$; NumPy's `np.std()` uses $n$ unless you pass `ddof=1`. On the example data this is the difference between 5.122 and 4.992 hours. The interquartile range and the 1.5 IQR fences used in boxplots: $$ \text{IQR} = Q_3 - Q_1, \qquad \text{lower fence} = Q_1 - 1.5\,\text{IQR}, \qquad \text{upper fence} = Q_3 + 1.5\,\text{IQR} $$ The z-score of one value, and the sample skewness and excess kurtosis (moment versions): $$ z_i = \frac{x_i - \bar{x}}{s}, \qquad g_1 = \frac{\frac{1}{n}\sum (x_i - \bar{x})^3}{\left[\frac{1}{n}\sum (x_i - \bar{x})^2\right]^{3/2}}, \qquad g_2 = \frac{\frac{1}{n}\sum (x_i - \bar{x})^4}{\left[\frac{1}{n}\sum (x_i - \bar{x})^2\right]^{2}} - 3 $$ Many programs report the adjusted versions $G_1 = g_1\sqrt{n(n-1)}/(n-2)$ and $G_2 = \frac{n-1}{(n-2)(n-3)}\left[(n+1)g_2 + 6\right]$. pandas `.skew()` and `.kurt()` give $G_1$ and $G_2$; SciPy's `stats.skew()` and `stats.kurtosis()` give $g_1$ and $g_2$ by default (`bias=True`). In small samples the two differ noticeably, so always say which one you report. **Quantiles are not unique.** With a finite sample, "the 25th percentile" can be defined in several ways. R offers nine types; the default is type 7 (linear interpolation between order statistics at position $1 + p(n-1)$). The pandas default `interpolation="linear"` and NumPy's default `method="linear"` are the same rule as R type 7, and so is Excel's `QUARTILE.INC`. SPSS and Minitab use type 6 (position $p(n+1)$, the same as Excel `QUARTILE.EXC`), and R's `fivenum()` and many textbooks use Tukey's hinges (medians of the lower and upper halves). On the example data $Q_1$ is 6.375 (type 7), 6.125 (type 6) and 6.25 (Tukey hinge). None is wrong; name the rule you use. | Symbol | Meaning | |---|---| | $n$ | Sample size (number of observations) | | $x_i$ | Value of observation $i$ | | $\bar{x}$ | Sample mean | | $s$ | Sample standard deviation, divisor $n-1$ | | $Q_1, Q_3$ | First and third quartiles (25th and 75th percentiles) | | IQR | Interquartile range, $Q_3 - Q_1$ | | $z_i$ | Standardised score of observation $i$ | | $g_1, G_1$ | Sample skewness: moment version and adjusted version | | $g_2, G_2$ | Sample excess kurtosis (normal = 0): moment version and adjusted version | | $r$ | Pearson correlation coefficient | The EDA process as a loop, not a straight line: ```mermaid flowchart TD A["Write purpose and questions"] --> B["Data overview: unit, n, variables, levels, missing"] B --> C["Univariate: centre, spread, shape"] C --> D["Charts: histogram and boxplot"] D --> E{"Unusual values?"} E -- "Yes" --> F["Check source: error or real?"] F --> G["Log decision in cleaning log"] G --> C E -- "No" --> H["Bivariate: scatter, crosstab, group summaries"] H --> I["Write findings as questions"] I --> J["Confirmatory analysis on a plan or new data"] ``` ## 4. Univariate summaries ### 4.1 Numeric variables | Variable | n | Mean | SD | Min | Q1 | Median | Q3 | Max | IQR | Skewness $G_1$ | Excess kurtosis $G_2$ | |---|---|---|---|---|---|---|---|---|---|---|---| | *Example:* hours | 20 | 9.00 | 5.12 | 3.5 | 6.375 | 8.0 | 10.0 | 28.0 | 3.625 | 2.86 | 10.37 | | *Example:* score | 20 | 65.60 | 8.09 | 50 | 59.75 | 65.5 | 70.75 | 81 | 11.0 | [ ] | [ ] | | [variable] | [ ] | [ ] | [ ] | [ ] | [ ] | [ ] | [ ] | [ ] | [ ] | [ ] | [ ] | Quartiles: R type 7 (the R, pandas and NumPy default). **Reading the study hours row.** The mean (9.0 hours) is above the median (8.0 hours), the maximum (28) is far above $Q_3$ (10), and the skewness is strongly positive (2.86). All three clues say the same thing: the distribution has a long right tail, driven mostly by one student. For a skewed variable, the median and IQR describe a typical student better than the mean and SD. The excess kurtosis of 10.4 says the tail is far heavier than a normal distribution's, again because of one value: kurtosis is very sensitive to single extreme values, so read it as a flag rather than a precise measurement. ### 4.2 Categorical variables: frequency tables | Mode of transport | Count | Percent | Cumulative percent | |---|---|---|---| | *Example:* Bus | 8 | 40% | 40% | | *Example:* Walk | 6 | 30% | 70% | | *Example:* Bike | 3 | 15% | 85% | | *Example:* Car | 3 | 15% | 100% | | **Total** | **20** | **100%** | | Sort nominal categories by frequency (as above); keep ordinal categories in their natural order (Year 1, 2, 3). With small samples, give counts as well as percentages: "15%" sounds like a lot until the reader learns it is 3 students. Cumulative percent only makes sense for ordinal variables or for the sorted list of categories, as here. | Year of study | Count | Percent | |---|---|---| | *Example:* 1 | 7 | 35% | | *Example:* 2 | 7 | 35% | | *Example:* 3 | 6 | 30% | ## 5. Charts **Histogram of study hours** (bin width 4 hours). Most students study 4 to 12 hours a week; the bar at 28 to 32 hours is a single student, separated from the rest by an empty gap. A gap like that is the visual signature of an outlier or of a second group. ```vega-lite { "$schema": "https://vega.github.io/schema/vega-lite/v5.json", "title": "Weekly study hours, Lakemoor STAT 101 (n = 20, fictional)", "width": 380, "height": 220, "data": {"values": [ {"hours": 3.5}, {"hours": 6.0}, {"hours": 7.5}, {"hours": 8.0}, {"hours": 5.0}, {"hours": 9.0}, {"hours": 10.0}, {"hours": 12.5}, {"hours": 6.5}, {"hours": 7.0}, {"hours": 8.5}, {"hours": 9.0}, {"hours": 11.0}, {"hours": 4.5}, {"hours": 13.0}, {"hours": 8.0}, {"hours": 7.0}, {"hours": 6.0}, {"hours": 10.0}, {"hours": 28.0} ]}, "mark": "bar", "encoding": { "x": {"field": "hours", "bin": {"step": 4, "extent": [0, 32]}, "title": "Study hours per week"}, "y": {"aggregate": "count", "title": "Number of students"} } } ``` The bin counts behind the chart (from `np.histogram` with edges 0, 4, ..., 32) are 1, 8, 8, 2, 0, 0, 0, 1: one student below 4 hours, eight in 4 to under 8, eight in 8 to under 12, two in 12 to under 16 and one in 28 to 32 (the 28-hour student). Each bin includes its lower edge and excludes its upper edge, so a value of exactly 8 counts in the 8 to 12 bin. **Boxplot of study hours by year.** The box spans $Q_1$ to $Q_3$, the line inside is the median, the whiskers reach the most extreme values inside the 1.5 IQR fences, and points beyond the fences are drawn separately. ```vega-lite { "$schema": "https://vega.github.io/schema/vega-lite/v5.json", "title": "Study hours by year of study (fictional)", "width": 300, "height": 220, "data": {"values": [ {"year": "Year 1", "hours": 3.5}, {"year": "Year 1", "hours": 6.0}, {"year": "Year 1", "hours": 7.5}, {"year": "Year 1", "hours": 8.0}, {"year": "Year 1", "hours": 5.0}, {"year": "Year 1", "hours": 9.0}, {"year": "Year 1", "hours": 10.0}, {"year": "Year 2", "hours": 12.5}, {"year": "Year 2", "hours": 6.5}, {"year": "Year 2", "hours": 7.0}, {"year": "Year 2", "hours": 8.5}, {"year": "Year 2", "hours": 9.0}, {"year": "Year 2", "hours": 11.0}, {"year": "Year 2", "hours": 4.5}, {"year": "Year 3", "hours": 13.0}, {"year": "Year 3", "hours": 8.0}, {"year": "Year 3", "hours": 7.0}, {"year": "Year 3", "hours": 6.0}, {"year": "Year 3", "hours": 10.0}, {"year": "Year 3", "hours": 28.0} ]}, "mark": {"type": "boxplot", "extent": 1.5}, "encoding": { "x": {"field": "year", "type": "nominal", "title": "Year of study"}, "y": {"field": "hours", "type": "quantitative", "title": "Study hours per week"} } } ``` The medians rise from 7.5 (Year 1) to 8.5 (Year 2) and 9.0 hours (Year 3), and the boxes overlap heavily. Year 3 contains the 28-hour value, drawn as a separate point: within Year 3 the quartiles are 7.25 and 12.25, so the upper fence is 19.75 and 28 lies beyond it. **Scatter plot of score against hours, with a least-squares line.** The chart shows how one point can bend a trend line, which is the main bivariate lesson in this data set. ```vega-lite { "$schema": "https://vega.github.io/schema/vega-lite/v5.json", "title": "Midterm score against study hours (fictional)", "width": 380, "height": 240, "data": {"values": [ {"hours": 3.5, "score": 50}, {"hours": 6.0, "score": 58}, {"hours": 7.5, "score": 66}, {"hours": 8.0, "score": 61}, {"hours": 5.0, "score": 57}, {"hours": 9.0, "score": 70}, {"hours": 10.0, "score": 68}, {"hours": 12.5, "score": 78}, {"hours": 6.5, "score": 64}, {"hours": 7.0, "score": 59}, {"hours": 8.5, "score": 67}, {"hours": 9.0, "score": 65}, {"hours": 11.0, "score": 75}, {"hours": 4.5, "score": 55}, {"hours": 13.0, "score": 74}, {"hours": 8.0, "score": 69}, {"hours": 7.0, "score": 62}, {"hours": 6.0, "score": 60}, {"hours": 10.0, "score": 73}, {"hours": 28.0, "score": 81} ]}, "layer": [ { "mark": {"type": "point", "filled": true, "size": 60}, "encoding": { "x": {"field": "hours", "type": "quantitative", "title": "Study hours per week"}, "y": {"field": "score", "type": "quantitative", "title": "Midterm score", "scale": {"zero": false}} } }, { "transform": [{"regression": "score", "on": "hours"}], "mark": {"type": "line", "color": "#c0392b"}, "encoding": { "x": {"field": "hours", "type": "quantitative"}, "y": {"field": "score", "type": "quantitative"} } } ] } ``` The fitted line (computed in R and Python below) is score = 54.23 + 1.26 × hours. Without student 20 it becomes score = 43.05 + 2.72 × hours: the slope more than doubles. Among the 19 other students, the points lie close to a steep straight line; student 20 sits far to the right and well below where that line would predict (about 119 points, which is impossible on a 100-point exam). That is a sign that the relationship flattens at high study time, or that the 28 hours is a misreport, or both. ## 6. Bivariate summaries | Pair | Summary | Value | Note | |---|---|---|---| | *Example:* hours and score | Pearson r | 0.800 | Pulled down by student 20; 0.939 without that student | | *Example:* hours and score | Spearman rho | 0.945 | Rank-based, so student 20 counts as "the most hours", not "28 hours" | | *Example:* year and hours | Group medians (Year 1, 2, 3) | 7.5, 8.5, 9.0 | Year 3 mean 12.0 is inflated by one value | | *Example:* mode and year | Crosstab | see below | Bus is the largest group in every year | | [pair] | [summary] | [value] | [note] | Group summaries of study hours by year: | Year | n | Mean | Median | SD | |---|---|---|---|---| | *Example:* 1 | 7 | 7.00 | 7.5 | 2.29 | | *Example:* 2 | 7 | 8.43 | 8.5 | 2.73 | | *Example:* 3 | 6 | 12.00 | 9.0 | 8.22 | The Year 3 SD (8.22) is three times that of the other years, but the Year 3 median is only half an hour above Year 2. If you compared year groups with a test of means, one student would drive the result. Note that as a question for the confirmatory stage, not as a finding. Crosstab of transport mode by year (counts): | Mode | Year 1 | Year 2 | Year 3 | Total | |---|---|---|---|---| | *Example:* Bike | 1 | 1 | 1 | 3 | | *Example:* Bus | 3 | 3 | 2 | 8 | | *Example:* Car | 1 | 1 | 1 | 3 | | *Example:* Walk | 2 | 2 | 2 | 6 | | **Total** | **7** | **7** | **6** | **20** | With cells of 1 to 3 students, no percentage comparison across years is meaningful. Small cells are a finding in themselves: they tell the analysis plan that transport mode cannot be used as a grouping factor without merging categories or collecting more data. ## 7. Outlier screening Two common screening rules: - **Tukey fences:** flag values below $Q_1 - 1.5\,\text{IQR}$ or above $Q_3 + 1.5\,\text{IQR}$ (use 3 IQR for "far out" values). The rule is resistant: the fences are built from quartiles, which the outlier itself barely moves. - **z-score rule:** flag $|z| > 3$ (some use 2.5). This rule is not resistant: an extreme value inflates the SD, which shrinks its own z-score. In samples of 10 or fewer observations, $|z|$ cannot even reach 3: the largest possible value is $(n-1)/\sqrt{n}$, which is 2.85 for $n = 10$. A resistant variant uses the median and the median absolute deviation (MAD) instead of the mean and SD. | Variable | Rule | Lower limit | Upper limit | Flagged ids | Checked against source | Decision | |---|---|---|---|---|---|---| | *Example:* hours | 1.5 IQR (type 7) | 0.94 | 15.44 | 20 (28 h) | Survey form says 28; student confirmed "about 4 hours a day" | Keep; run main analysis with and without | | *Example:* hours | \|z\| > 3 | | | 20 (z = 3.71) | as above | as above | | *Example:* score | 1.5 IQR (type 7) | [ ] | [ ] | [ ] | [ ] | [ ] | | [variable] | [rule] | [ ] | [ ] | [ ] | [ ] | [ ] | **Decision rules for flagged values.** Correct a value only when the source shows it is a recording error (a typo, a unit mix-up such as minutes typed as hours). Delete it only if it is impossible (a score of 140 on a 100-point exam) and cannot be corrected. Keep a value that is real but extreme, and show how much it matters by running the key summaries with and without it (a sensitivity analysis). Never drop a value because it spoils a result. Write every decision, with its reason and date, in `academic/stats-data-cleaning-log`. ## 8. Code and output (R and Python) ### R (base R and stats only) The R block uses only base R and the `stats` package, so it runs on any fresh installation. Skewness and kurtosis are computed by hand because base R has no function for them. ```r df <- data.frame( id = 1:20, year = c(1,1,1,1,1,1,1,2,2,2,2,2,2,2,3,3,3,3,3,3), mode = c("Bus","Walk","Bus","Car","Walk","Bus","Bike","Walk","Bus","Bus", "Car","Bike","Walk","Bus","Car","Walk","Bus","Bike","Walk","Bus"), hours = c(3.5,6,7.5,8,5,9,10,12.5,6.5,7,8.5,9,11,4.5,13,8,7,6,10,28), score = c(50,58,66,61,57,70,68,78,64,59,67,65,75,55,74,69,62,60,73,81) ) h <- df$hours summary(h) c(mean = mean(h), sd = sd(h), IQR = IQR(h)) # sd() uses n - 1 quantile(h, c(.25, .75)) # type 7 (default) quantile(h, c(.25, .75), type = 6) # SPSS / Minitab style fivenum(h) # Tukey hinges n <- length(h); m <- mean(h) g1 <- sum((h - m)^3) / n / (sum((h - m)^2) / n)^1.5 G1 <- g1 * sqrt(n * (n - 1)) / (n - 2) c(skew_g1 = g1, skew_G1 = G1) q <- quantile(h, c(.25, .75)); f <- q + c(-1.5, 1.5) * IQR(h) f # 1.5 IQR fences df$id[h < f[1] | h > f[2]] # flagged ids z <- (h - mean(h)) / sd(h); round(max(z), 2) table(df$mode) round(prop.table(table(df$mode)), 2) aggregate(hours ~ year, data = df, FUN = function(x) c(mean = mean(x), median = median(x), sd = sd(x))) cor(df$hours, df$score) cor(df$hours, df$score, method = "spearman") coef(lm(score ~ hours, data = df)) boxplot.stats(h)$out # boxplot's own outlier list (uses hinges) g2 <- sum((h - m)^4) / n / (sum((h - m)^2) / n)^2 - 3 G2 <- ((n + 1) * g2 + 6) * (n - 1) / ((n - 2) * (n - 3)) c(kurt_g2 = g2, kurt_G2 = G2) ``` ```text Min. 1st Qu. Median Mean 3rd Qu. Max. 3.500 6.375 8.000 9.000 10.000 28.000 mean sd IQR 9.000000 5.122191 3.625000 25% 75% 6.375 10.000 25% 75% 6.125 10.000 [1] 3.50 6.25 8.00 10.00 28.00 skew_g1 skew_G1 2.637158 2.855982 25% 75% 0.9375 15.4375 [1] 20 [1] 3.71 Bike Bus Car Walk 3 8 3 6 Bike Bus Car Walk 0.15 0.40 0.15 0.30 year hours.mean hours.median hours.sd 1 1 7.000000 7.500000 2.291288 2 2 8.428571 8.500000 2.729905 3 3 12.000000 9.000000 8.221922 [1] 0.8004005 [1] 0.9446344 (Intercept) hours 54.225878 1.263791 [1] 28 kurt_g2 kurt_G2 7.669229 10.372623 ``` ### Python (pandas, numpy, scipy) The Python block builds the same data frame and prints the same summaries, plus the correlations and slopes with and without student 20 used in sections 5 and 6. ```python import numpy as np, pandas as pd from scipy import stats df = pd.DataFrame({ "id": range(1, 21), "year": [1,1,1,1,1,1,1,2,2,2,2,2,2,2,3,3,3,3,3,3], "mode": ["Bus","Walk","Bus","Car","Walk","Bus","Bike","Walk","Bus","Bus", "Car","Bike","Walk","Bus","Car","Walk","Bus","Bike","Walk","Bus"], "hours": [3.5,6,7.5,8,5,9,10,12.5,6.5,7,8.5,9,11,4.5,13,8,7,6,10,28], "score": [50,58,66,61,57,70,68,78,64,59,67,65,75,55,74,69,62,60,73,81], }) h = df["hours"] print(df[["hours", "score"]].describe().round(2)) print("skewness (adjusted):", round(h.skew(), 3), " excess kurtosis:", round(h.kurt(), 3)) print("scipy skew (bias=True):", round(stats.skew(h), 3), " kurtosis:", round(stats.kurtosis(h), 3)) print("SD n-1:", round(h.std(), 3), " SD n:", round(h.std(ddof=0), 3)) print("Q1, Q3 linear (R type 7):", h.quantile([.25, .75]).tolist()) print("Q1, Q3 weibull (R type 6):", np.quantile(h, [.25, .75], method="weibull").tolist()) q1, q3 = h.quantile([.25, .75]); iqr = q3 - q1 print("IQR:", iqr, " fences:", q1 - 1.5*iqr, q3 + 1.5*iqr) z = (h - h.mean()) / h.std() print("|z| > 3:", df.loc[z.abs() > 3, "id"].tolist(), " max z:", round(z.max(), 2)) print(df["mode"].value_counts()) print(pd.crosstab(df["mode"], df["year"], margins=True)) print(df.groupby("year")["hours"].agg(["count", "mean", "median", "std"]).round(2)) print("Pearson r:", round(h.corr(df["score"]), 3), " Spearman rho:", round(h.corr(df["score"], method="spearman"), 3)) d2 = df[df["id"] != 20] print("Pearson r without id 20:", round(d2["hours"].corr(d2["score"]), 3)) print("slope, intercept:", np.round(np.polyfit(h, df["score"], 1), 3)) print("slope, intercept without id 20:", np.round(np.polyfit(d2["hours"], d2["score"], 1), 3)) print(np.histogram(h, bins=range(0, 36, 4))) ``` ```text hours score count 20.00 20.00 mean 9.00 65.60 std 5.12 8.09 min 3.50 50.00 25% 6.38 59.75 50% 8.00 65.50 75% 10.00 70.75 max 28.00 81.00 skewness (adjusted): 2.856 excess kurtosis: 10.373 scipy skew (bias=True): 2.637 kurtosis: 7.669 SD n-1: 5.122 SD n: 4.992 Q1, Q3 linear (R type 7): [6.375, 10.0] Q1, Q3 weibull (R type 6): [6.125, 10.0] IQR: 3.625 fences: 0.9375 15.4375 |z| > 3: [20] max z: 3.71 mode Bus 8 Walk 6 Car 3 Bike 3 Name: count, dtype: int64 year 1 2 3 All mode Bike 1 1 1 3 Bus 3 3 2 8 Car 1 1 1 3 Walk 2 2 2 6 All 7 7 6 20 count mean median std year 1 7 7.00 7.5 2.29 2 7 8.43 8.5 2.73 3 6 12.00 9.0 8.22 Pearson r: 0.8 Spearman rho: 0.945 Pearson r without id 20: 0.939 slope, intercept: [ 1.264 54.226] slope, intercept without id 20: [ 2.717 43.051] (array([1, 8, 8, 2, 0, 0, 0, 1]), array([ 0, 4, 8, 12, 16, 20, 24, 28, 32])) ``` **Where R and Python agree and differ, and why.** | Quantity | R | Python | Reason | |---|---|---|---| | SD | `sd()` = 5.122 | `Series.std()` = 5.122; `np.std()` = 4.992 | pandas and R divide by $n-1$; NumPy divides by $n$ unless `ddof=1` | | Quartiles | `quantile()` type 7: 6.375, 10 | `Series.quantile()` linear: 6.375, 10 | Same rule (R type 7 = pandas and NumPy "linear"); type 6 gives 6.125 | | Boxplot hinges | `fivenum()` 6.25, 10 | not built in | Tukey hinges are a third rule; base R `boxplot()` uses them | | Skewness | no base function; $g_1$ = 2.637, $G_1$ = 2.856 by hand | pandas `.skew()` = 2.856 ($G_1$); `stats.skew()` = 2.637 ($g_1$) | Adjusted versus moment formula | | Excess kurtosis | by hand: $g_2$ = 7.669, $G_2$ = 10.373 | pandas `.kurt()` = 10.373; `stats.kurtosis()` = 7.669 | Adjusted versus moment formula; both subtract 3 | ## 9. Worked examples ### Worked example A (fully solved): five-number summary and fences for study hours **Step 1. Sort the 20 values.** 3.5, 4.5, 5, 6, 6, 6.5, 7, 7, 7.5, 8, 8, 8.5, 9, 9, 10, 10, 11, 12.5, 13, 28. **Step 2. Median.** With $n = 20$ (even), the median is the mean of the 10th and 11th values: $(8 + 8)/2 = 8.0$. **Step 3. Quartiles, type 7.** The position of $Q_1$ is $1 + 0.25 \times 19 = 5.75$: three quarters of the way from the 5th value (6) to the 6th (6.5), so $Q_1 = 6 + 0.75 \times 0.5 = 6.375$. The position of $Q_3$ is $1 + 0.75 \times 19 = 15.25$: between the 15th and 16th values, both 10, so $Q_3 = 10$. **Step 4. Five-number summary.** Minimum 3.5, $Q_1$ 6.375, median 8.0, $Q_3$ 10.0, maximum 28. **Step 5. Fences.** IQR $= 10 - 6.375 = 3.625$; $1.5 \times 3.625 = 5.4375$. Lower fence $6.375 - 5.4375 = 0.9375$, upper fence $10 + 5.4375 = 15.4375$. Only 28 lies outside, so student 20 is flagged. **Step 6. z-score check.** $z_{20} = (28 - 9.0)/5.122 = 3.71$, so the z-rule flags the same student. **Step 7. Interpretation.** Typical weekly study time is about 8 hours (median), with the middle half of students between about 6.4 and 10 hours. The distribution is right-skewed (mean 9.0 above median 8.0, $G_1$ = 2.86) because of one student who reports 28 hours. Report the median and IQR as the headline and the mean and SD in a table, and check the 28-hour value against its source. **Step 8. Report sentence.** "Weekly study time was right-skewed (*Mdn* = 8.0 hours, IQR = 3.6, range 3.5 to 28; *M* = 9.0, *SD* = 5.1, *n* = 20). One student reported 28 hours, beyond the upper 1.5 IQR fence (15.4 hours); the value was confirmed with the student and retained." ### Worked example B (half solved): study hours in Year 3 only Year 3 values: 13, 8, 7, 6, 10, 28 ($n = 6$). 1. Sorted values: 6, 7, 8, 10, 13, 28. 2. Median: mean of the 3rd and 4th values $= (8 + 10)/2 =$ [ ]. 3. $Q_1$ (type 7): position $1 + 0.25 \times 5 = 2.25$, so $Q_1 = 7 + 0.25 \times (8 - 7) =$ [ ]. 4. $Q_3$ (type 7): position $1 + 0.75 \times 5 =$ [ ], so $Q_3 = 10 + 0.75 \times (13 - 10) =$ [ ]. 5. IQR $=$ [ ]; upper fence $=$ [ ]. Is 28 flagged within Year 3? [ ] 6. Mean $= 72/6 = 12.0$. Why is the mean 3 hours above the median? [ ] 7. The Year 3 SD is 8.22. Would you compare years using means and SDs or medians and IQRs? Why? [ ] ### Exercise C (unsolved): midterm score Using the 20 scores in section 2: 1. Give the five-number summary of score (type 7 quartiles). 2. Compute the IQR and the 1.5 IQR fences. Are any scores flagged? 3. Compare the mean and the median. What do they suggest about the shape? 4. Check your answers with `quantile(df$score)` in R or `df["score"].describe()` in Python. 5. Write a one-sentence description of the score distribution in the style of Example A, step 8. ## 10. Findings and questions for confirmatory analysis | # | Observation (exploratory) | Evidence | Question for confirmatory analysis | Planned method | |---|---|---|---|---| | *Example:* F1 | Study hours are right-skewed with one extreme value | Mean 9.0 > median 8.0; $G_1$ = 2.86; id 20 at 28 h | Does the main result hold with and without id 20? | Sensitivity analysis; `academic/stats-data-transformations` (log hours) | | *Example:* F2 | Score rises with hours, possibly flattening at high hours | r = 0.80, rho = 0.945; slope 2.72 without id 20 vs 1.26 with | Is the relation linear, or does it level off? | Regression with log(hours) or a spline, pre-specified | | *Example:* F3 | Year 3 studies a little more in the median | Medians 7.5, 8.5, 9.0; heavy overlap | Is year a confounder for hours and score? | Include year as a covariate; `academic/stats-ancova` | | *Example:* F4 | Transport groups are tiny | Cells of 1 to 3 | None with these data; collect more or merge categories | No test | | [F5] | [observation] | [evidence] | [question] | [method] | **Separating exploration from confirmation.** Three practical ways to keep the two apart: 1. **Write down the confirmatory questions before you collect the confirmatory data** (a pre-analysis plan or preregistration, see `academic/pre-analysis-plan` and `science/preregistration`). 2. **Split the data:** explore on a random part (for example 30%) and confirm on the rest, which you do not look at until the plan is fixed. 3. **Label clearly:** if neither is possible, report the patterns as exploratory, without p-values, or with a sentence saying that the hypotheses were generated by the same data. **Summary paragraph.** [Example: The 20 students studied a median of 8 hours a week (IQR 3.6). One student reported 28 hours; the value was confirmed and kept. Midterm scores (*M* = 65.6, *SD* = 8.1) were roughly symmetric with no flagged values. Scores rose with study hours (Spearman rho = 0.95), and the pattern appeared to flatten at high study time. These observations are exploratory; the analysis plan will test the hours-score relation with year as a covariate and with a sensitivity analysis for the extreme value.] ## 11. Reporting checklist - [ ] Purpose and questions written before the analysis. - [ ] Unit of observation, sample size n and data source stated. - [ ] Every variable listed with its measurement level and missing count. - [ ] Centre, spread and shape reported for every key numeric variable, with both mean/SD and median/IQR when skewed. - [ ] Quantile rule and skewness formula named (for example "type 7 quartiles; adjusted skewness $G_1$"). - [ ] Frequency tables with counts and percentages for categorical variables. - [ ] At least one chart per key variable, with a one-sentence caption that says what it shows. - [ ] Outlier rule stated in advance; each flagged value checked against the source and the decision logged. - [ ] Bivariate summaries match the pair of measurement levels; n given for every group. - [ ] Findings written as questions; no confirmatory claims or p-values for patterns found here. - [ ] Code and software versions included so that every number can be reproduced. ## 12. Common mistakes and misconceptions | Mistake or misconception | Why it is wrong | Do this instead | |---|---|---| | "The mean is the typical value." | In skewed data the mean is pulled toward the long tail (here 9.0 vs median 8.0) | Report the median and IQR for skewed variables, and say why | | Reporting only the mean and SD | Two very different distributions can share a mean and SD | Always look at a histogram or boxplot as well | | Deleting every point beyond the 1.5 IQR fences | The rule flags about 0.7% of values even from a perfect normal distribution; in skewed data it flags many real values | Treat fences as a screen; check the source; keep real values and run a sensitivity analysis | | Assuming R, Python, Excel and SPSS give the same quartiles | There are at least nine quantile rules; defaults differ across programs | Name the rule; use the same rule throughout | | Mixing `np.std()` (divisor n) with `sd()` (divisor n - 1) | Small but real difference (4.99 vs 5.12 here) | Use `ddof=1` in NumPy for a sample SD | | Testing a pattern on the same data that suggested it | The p-value ignores the search that found the pattern, so false positives rise | Write it as a question; test it on new or held-out data or with a pre-specified plan | | Reading a correlation of 0.80 as "strong linear relation" without a scatter plot | One point can raise or lower r a lot (0.80 with, 0.94 without student 20) | Always plot; report Spearman as well when outliers or curvature appear | | Using percentages alone for small groups | "15%" may be 3 people | Give counts and n | | Choosing histogram bins to make the picture look normal | Bin width can hide gaps or invent bumps | Try two or three bin widths; report the one that shows the data honestly | | Treating a high kurtosis value as a precise measure of "peakedness" | Kurtosis mainly measures tail weight and is dominated by the most extreme values | Use it as a flag for heavy tails or outliers, together with a plot | ## 13. Self-test questions 1. A variable has mean 42 and median 31. What does that suggest about its shape, and which summary of centre would you report as the headline? 2. Why is the 1.5 IQR rule more resistant to outliers than the $|z| > 3$ rule? 3. A colleague's pandas output gives $Q_1$ = 6.375 for study hours and your SPSS output gives 6.125. Who is wrong? 4. NumPy gives an SD of 4.99 and R gives 5.12 for the same data. Explain the difference and how to make them agree. 5. During EDA you notice that Year 3 students score higher than Year 1 students and you run a t-test on the same data, getting p = .02. What is the problem with reporting this as a confirmed finding? 6. Which summary and chart would you use for (a) transport mode, (b) year of study, (c) exam score by transport mode? ## 14. My notes and next steps | Item | Entry | |---|---| | Confidence with descriptive statistics (1-5) | [ ] | | Confidence with choosing charts (1-5) | [ ] | | Confidence with outlier decisions (1-5) | [ ] | | Open questions | [What I am still unsure about] | | Mistakes to log | [Link to `academic/stats-error-log` entry] | | Next review dates | YYYY-MM-DD (2 days), YYYY-MM-DD (1 week), YYYY-MM-DD (1 month) | | Next template | [`academic/stats-data-transformations` / `academic/stats-correlation` / `academic/stats-analysis-plan`] | ## 15. Sources and further reading Structure aligned with the exploratory data analysis approach of Tukey (1977), the EDA chapter of the NIST/SEMATECH e-Handbook of Statistical Methods and the GAISE College Report recommendations; all text, data and examples here are original. - Tukey, J. W. (1977). *Exploratory Data Analysis*. Addison-Wesley. - NIST/SEMATECH e-Handbook of Statistical Methods, chapter 1 Exploratory Data Analysis, https://www.itl.nist.gov/div898/handbook/ (accessed 2026-10-25). Public domain. - Hyndman, R. J., and Fan, Y. (1996). Sample quantiles in statistical packages. *The American Statistician*, 50(4), 361–365. - Joanes, D. N., and Gill, C. A. (1998). Comparing measures of sample skewness and kurtosis. *Journal of the Royal Statistical Society, Series D (The Statistician)*, 47(1), 183–189. - Wickham, H., and Stryjewski, L. (2012). *40 years of boxplots*. Technical report, had.co.nz. - GAISE College Report ASA Revision Committee (2016). *Guidelines for Assessment and Instruction in Statistics Education: College Report 2016*. American Statistical Association. - Further reading (by name only): de Smith, M. J. (2024). *Statistical Analysis Handbook* (2024 edition, StatsRef); Diez, D. M., Çetinkaya-Rundel, M., and Barr, C. D., *OpenIntro Statistics* (CC BY-SA). - Software documentation: R `quantile()`, `fivenum()`, `boxplot.stats()`; pandas `DataFrame.describe()`, `Series.skew()`, `Series.kurt()`; SciPy `stats.skew()`, `stats.kurtosis()`. ## Answer key **Worked example B.** (2) Median = 9.0. (3) $Q_1$ = 7.25. (4) Position 4.75, $Q_3$ = 12.25. (5) IQR = 5.0; upper fence $12.25 + 7.5 = 19.75$; yes, 28 is flagged within Year 3 as well. (6) The single value of 28 pulls the mean upward; the median, based on the middle two values, does not move. Without 28, the Year 3 mean would be 8.8. (7) Medians and IQRs, because one value dominates the Year 3 mean and SD; with groups of 6 or 7, one student can change the story. If you report means, show them with and without the extreme value. **Exercise C.** (1) Minimum 50, $Q_1$ = 59.75, median 65.5, $Q_3$ = 70.75, maximum 81. (2) IQR = 11.0; fences $59.75 - 16.5 = 43.25$ and $70.75 + 16.5 = 87.25$; no score is flagged. (3) Mean 65.6 and median 65.5 are almost equal, which suggests a roughly symmetric distribution (the histogram confirms one peak, no gaps). (4) R gives 59.75, 65.50, 70.75; pandas gives the same (25% 59.75, 50% 65.50, 75% 70.75). (5) "Midterm scores were roughly symmetric (*M* = 65.6, *SD* = 8.1; *Mdn* = 65.5, IQR = 11.0, range 50 to 81, *n* = 20), with no values beyond the 1.5 IQR fences." **Self-test.** 1. Mean well above the median suggests right skew (a long tail of high values). Report the median as the headline, with the IQR, and give the mean and SD in a table. 2. The fences are built from quartiles, which barely move when one value becomes extreme. The z-rule uses the mean and SD, which the extreme value itself inflates, so it partly hides itself (masking); with several outliers this gets worse. 3. Neither. pandas uses linear interpolation (R type 7) and SPSS uses type 6. Both are valid quantile definitions; say which one you used and keep it consistent. 4. NumPy's `np.std()` divides by $n$ (population formula); R's `sd()` divides by $n-1$ (sample formula). Use `np.std(x, ddof=1)` or pandas `.std()` to match R. 5. The hypothesis was suggested by the same data used to test it. The p-value assumes the test was chosen in advance, so it understates the chance of a false positive. Report the difference as exploratory, or test it on new or held-out data with a plan written beforehand. 6. (a) Nominal: frequency table with counts and percentages, bar chart. (b) Ordinal: frequency table in natural order, bar chart; median and percentiles if useful. (c) Numeric by groups: group summaries (n, median, IQR, mean, SD) and side-by-side boxplots (with tiny groups, a strip plot of the raw points is better). ## Revision history | Version | Date | Author | Change summary | |---|---|---|---| | 0.1 | YYYY-MM-DD | [Your name] | First draft from the template | | 0.2 | YYYY-MM-DD | [Your name] | *Example:* Replaced Lakemoor example with thesis survey data; outlier decisions logged |