Exploratory Data Analysis Report Template (Statistics, R and Python)

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

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.

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].

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):

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:

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:

The z-score of one value, and the sample skewness and excess kurtosis (moment versions):

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.

The EDA process as a loop, not a straight line:

4. Univariate summaries

4.1 Numeric variables

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

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.

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.

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.

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.

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

Group summaries of study hours by year:

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):

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.

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.

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)

   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.

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)))

       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.

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

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

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

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