Some Friendly Advice for Data Analysis

Who am I to give advice?

• Statistical Consultant at UVA Library StatLab since 2013
• Have consulted with hundreds of students on statistical methods, data wrangling, creating visualizations and debugging code
• Have seen where students get stuck and what they struggle with

Agenda

• Share some practical tips on data analysis
• Occasional stats refresher
• Perhaps help someone save time and avoid pain and stress

Make a statistical plan

• Some students design a study with no plan of analysis
• Consequently, the structure of the data they collect is too complicated for their statistical comfort
• If you’re designing a study that requires a statistical analysis, identify the statistical analysis before you collect the data

Example of a statistical plan

Sample Size

• If you know your analysis plan, you may want to think about your sample size.
• Example: I suspect US quarters do not land tails 50% of the time. I’m going to flip a quarter a bunch of times and run a one-sample test for proportions.
• How many times do I need to flip a quarter?
• The answer depends on three things:
  • Hypothesized meaningful effect size. What do I think the real probability of landing tails is?
  • Desired power. What probability do I want to have of correctly rejecting a null hypothesis of fair coin?
  • Acceptable significance level. What probability am I willing to accept for incorrectly rejecting the null hypothesis of fair coin?
• There is no single correct answer to sample size planning.

Sample size calculation

This type of calculation is easy to do in R. We need to flip the coin about 1048 times assuming an effect size of 0.55 versus 0.50, power of 0.90, and significance level of 0.05.

library(pwr)
pwr.p.test(h = ES.h(0.5, 0.55), sig.level = 0.05, power = 0.90)

     proportion power calculation for binomial distribution (arcsine transformation) 

              h = 0.1001674
              n = 1047.232
      sig.level = 0.05
          power = 0.9
    alternative = two.sided

Think carefully about sample size. Don’t assume <insert number> is enough.

Precision

• Another goal for an analysis may be to make a precise estimate.
• Example: I suspect US quarters do not land tails 50% of the time. I’m going to flip a quarter a bunch of times and calculate a 95% confidence interval on the observed proportion of tails.
• How many times do I need to flip a quarter?
• The answer depends on two things:
  • Hypothesized probability. What do I think the real probability of landing tails is?
  • Desired width of confidence interval. How precise do I want to make my measurement.
• Again, no single correct answer.

Sample size calculation

This type of calculation is also easy to do in R. We need to flip the coin about 2374 times assuming proportion of tails will be 0.55 and our CI width is 0.04 (ie, \(\pm\) 0.02).

library(presize)
prec_prop(p = 0.55, conf.width = 0.04, method = "agresti-coull")

     sample size for a proportion with Agresti-Coull confidence interval. 

     p      padj        n conf.width conf.level       lwr       upr
1 0.55 0.5499192 2373.139       0.04       0.95 0.5299192 0.5699192

NOTE: padj is the adjusted proportion, from which the ci is calculated.

Think carefully about sample size. Don’t assume <insert number> is enough.

Cleaning Data

• Precious time is often spent cleaning data instead of analyzing it because of how data was entered in Excel or Google Sheets
• This usually ties back to not knowing the statistical analysis for the study design
• This time can be reduced by adhering to a few simple rules
• Broman, K. W., & Woo, K. H. (2018). Data Organization in Spreadsheets. The American Statistician, 72(1), 2–10.
  https://doi.org/10.1080/00031305.2017.1375989

How to Excel

• Organize the data as a single rectangle (with observations as rows, variables as columns, single header row, one row per observation)
• Do not include calculations in the raw data files
• Do not use font color or highlighting as data
• Write dates as YYYY-MM-DD
• Put just one thing in a cell

Not Good Excel

This will require cleaning before analyzing in R, Python, SPSS, etc.

Good Excel

This is ready for analysis in R, Python, SPSS, etc.

Not Good Excel

This will require cleaning before analyzing in R, Python, SPSS, etc.

Good Excel

This is ready for analysis in R, Python, SPSS, etc.

Get to Know Your Data

• Resist the urge to go directly to modeling/machine learning
• Tabulate categorical variables and cross-tabulate multiple categories
• Check range and distribution of numeric variables
• Determine where you have missingness

Logistic regression gone bad

Why are there singularities? Why are there NAs?

m <- glm(y ~ z + x + g + x:g, data = d, family = binomial)
summary(m)

Call:
glm(formula = y ~ z + x + g + x:g, family = binomial, data = d)

Coefficients: (2 not defined because of singularities)
            Estimate Std. Error z value Pr(>|z|)  
(Intercept)   0.1942     0.3609   0.538   0.5906  
z                 NA         NA      NA       NA  
xb            0.6168     0.5575   1.106   0.2686  
xc           -0.5306     0.4647  -1.142   0.2535  
gm            1.0868     0.6211   1.750   0.0802 .
xb:gm        -0.7991     0.9127  -0.876   0.3813  
xc:gm             NA         NA      NA       NA  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 199.84  on 147  degrees of freedom
Residual deviance: 186.56  on 143  degrees of freedom
  (2 observations deleted due to missingness)
AIC: 196.56

Number of Fisher Scoring iterations: 4

Simple investigation

Calling summary() on the data frame reveals no variability in predictor z along with two missing values.

summary(d)

 y      x      g             z    
 0:62   a:54   f:106   Min.   :1  
 1:88   b:47   m: 44   1st Qu.:1  
        c:49           Median :1  
                       Mean   :1  
                       3rd Qu.:1  
                       Max.   :1  
                       NA's   :2  

A cross-tabulation of x and g reveals a 0 count when x = "c" and g = "m".

xtabs(~ x + g, data = d)

   g
x    f  m
  a 31 23
  b 26 21
  c 49  0

No need to include z in model or to fit an interaction between x and g.

Absolute vs Relative Change

• Think carefully whether change is absolute (additive) or relative (multiplicative)
• In 2015 the World Health Organization (WHO) announced that eating 50g of processed meat a day was associated with an 18% increase in the risk of bowel cancer.
  • This is relative change. If the risk of bowel cancer without consuming 50g of processed meat per day is 0.001, then a relative increase of 18% raises the risk to 0.001 x 1.18 = 0.00118.
• In 2014, 6.94% of U.S. households had at least one motorcycle. In 2018, that figure rose to 8.02%. [source]
  • Can report as absolute change: the percent of households with at least one motorcycle increased by 1.08% (8.02% - 6.94%).
  • Can report as relative change: motorcycle ownership increased by 8.02/6.94 = 1.1556, or 15.56%.

Probability vs Odds

• Probability (p) ranges from 0 to 1.
• Odds = p/(1 - p) and range from 0 to infinity.
• If the probability of rain tomorrow is 0.2, the odds are 0.2/(1 - 0.2) = 1/4.
• If the probability of rain tomorrow is 0.8, the odds are 0.8/(1 - 0.8) = 4.
• An Odds Ratio is a ratio of odds and is often expressed as a relative increase.
• If the odds of cancer for males in 1/1000, and the odds of cancer for females is 1/1200, then the odds ratio is (1/1000)/(1/1200) = 1.2. The odds of cancer for males is 20% higher.

Odds Ratios are Relative

• A paper published in 2024 announced that coffee drinkers had a 20% lower risk of developing Parkinson disease
• The odds ratio was 0.80 with a 95% CI of 0.67 - 0.95
• If baseline probability of developing Parkinson disease is 0.001…
  • Odds = 0.001/(1 - 0.001) \(\approx\) 0.001
  • Odds for coffee drinkers: 0.001 * 0.8 \(\approx\) 0.0008
• Relative change can often sound more dramatic than it really is

Let’s get Normal

It’s ok to not be Normal

• Your data do not need to be normally distributed
• In linear statistical modeling/regression, the assumption of normality is on the residuals, not the variables in the model
• “For the purpose of estimating the regression line (as compared to predicting individual data points), the assumption of normality is typically barely important at all.” (Gelman et al., 2020, p. 155)
• There is no need to test for normality before doing t-tests (Rasch et al, 2011)
• Normality tests are sensitive to sample size; use with caution (or not at all)

Normality “tests”

shapiro.test(x)

    Shapiro-Wilk normality test

data:  x
W = 0.90239, p-value = 0.1703

shapiro.test(x)

    Shapiro-Wilk normality test

data:  x
W = 0.99903, p-value = 0.02258

Means and Medians

Some distributions

Mean and Median

Correlation Review

Correlation near 0

Correlations near 0 do not always mean there is no relationship. See the Datasaurus Dozen. Plot your data if calculating correlations.

The Quest for Significance

• Data analysis usually involves statistical inference: use sample statistics to infer something about population parameters
• We often calculate statistics for different groups and compare the statistics
• The big question: is the difference statistically significant?
• The question is usually answered with a p-value

p-value Review

• Let’s say we sample 100 females and 100 male UVA first-year students and ask them if they floss their teeth every day
• 35/100 females say “yes”; 28/100 males say “yes”
• Difference: 0.35 - 0.28 = 0.07
• Assume there is no difference in proportion who floss (i.e., if we could survey all first year students the proportion who floss would be exactly equal between groups). We usually call this the Null hypothesis.
• What’s the probability of seeing a difference of 0.07 or larger (in absolute value) in our sample if our assumption is true? That’s the p-value
• The p-value for this hypothesis test: 0.36
• Based on the p-value, we might fail to reject the Null hypothesis and conclude there is no “statistically significant” difference between the groups

How to perform test in R

prop.test(x = c(35,28), n = c(100,100))

    2-sample test for equality of proportions with continuity correction

data:  c(35, 28) out of c(100, 100)
X-squared = 0.8342, df = 1, p-value = 0.3611
alternative hypothesis: two.sided
95 percent confidence interval:
 -0.06838886  0.20838886
sample estimates:
prop 1 prop 2 
  0.35   0.28 

Statistical Significance

• Significance is usually achieved by crossing an arbitrary threshold, usually 0.05
• However, a p-value less than 0.05 is not proof of an effect or phenomenon
• Many statisticians are now calling to abolish the 0.05 threshold and the phrase “statistically significant”
• In 2019 The American Statistician published a special issue with 43 articles discussing p-values
• Wasserstein, R. L., Schirm, A. L., & Lazar, N. A. (2019). Moving to a World Beyond “p<0.05.” The American Statistician, 73(sup1), 1–19.
  https://doi.org/10.1080/00031305.2019.1583913

Handle with Caution

• p-values aren’t going anywhere just yet; know how to interpret them
• They do not measure the probability the null hypothesis is true or that the data were produced by random chance alone
• They do not justify proof of something because they’re less than 0.05
• They do not prove the null hypothesis because they’re greater than 0.05
• They do not measure the size of an effect or the importance of a result. A smaller p-value does not mean a “bigger” effect
• A p-value is the probability of the data, or more extreme data, assuming a null hypothesis is true.

p-Hacking

• The quest for p < 0.05 can tempt one to test unplanned hypotheses about their data
  • There was no difference in flossing between the males and females, but what about students from VA versus students out-of-state?
  • There was no difference in VA students vs out-of-state students, but what about students who use Macs vs students who use Windows?
• This is called p-hacking, where one performs many tests and only reports those with p < 0.05 as if those tests were planned from the beginning
• This can lead to false positive results that are wrong and don’t replicate in future experiments

p-Hacking example

Simulate data with absolutely no relationship, and then hack at the data until significance is achieved. The final result is a false positive.

# simulate data
set.seed(666)
n <- 200
y <- rnorm(n)
x1 <- rnorm(n)
x2 <- rnorm(n)
g1 <- sample(c("a","b"), size = n, replace = TRUE)
g2 <- sample(c("i","j","k"), size = n, replace = TRUE)
# first test
t.test(y ~ g1)$p.value

[1] 0.2258509

# second test
t.test(y ~ g1, subset = g2 == "j")$p.value

[1] 0.122574

# third test
t.test(y ~ g1, subset = g2 == "j" & x2 < 1)$p.value

[1] 0.06973587

# fourth test
t.test(y ~ g1, subset = g2 == "j" & x2 < 1 & x1 < 1)$p.value

[1] 0.03239254

Confidence Intervals

• Slightly better than p-values but no less slippery
• Measure uncertainty in estimated statistics
• Difference in proportions between females and males who floss: 0.07
• 95% confidence interval: [-0.07,0.21]
• Since the interval overlaps 0 we’re not sure which group flosses more, though it would seem to be females
• In general, confidence intervals are more informative than p-values. They tell us about direction and magnitude.

Confidence, not Probability

• Unfortunately, confidence intervals are not probability intervals
• We can’t say that there’s a 95% probability that the difference in proportions is between -0.07 and 0.21
• We call it “confidence” because we are confident in the method
• The method: if we take thousands of samples and calculate thousands of 95% confidence intervals, about 95% of the intervals should contain the true value

Confidence Interval Demo in R

Below we simulate asking 100 females and 100 males if they floss, assuming the true difference in proportions is 0.10. We then calculate the confidence interval on the difference in proportions and check if it contains 0.10. About 95% of the confidence intervals contain the true value of 0.10.

ci <- replicate(n = 1000, expr = {
  f_samp <- sample(c("yes", "no"), size = 100, replace = TRUE, prob = c(0.3, 0.7))
  m_samp <- sample(c("yes", "no"), size = 100, replace = TRUE, prob = c(0.2, 0.8))
  f_floss <- sum(f_samp == "yes")
  m_floss <- sum(m_samp == "yes")
  test <- prop.test(x = c(f_floss, m_floss), n = c(100,100))
  test$conf.int[1] < 0.1 & test$conf.int[2] > 0.1})

# proportion of CI that contain difference of 0.1
mean(ci)

[1] 0.964

Probability, not Confidence

• If we want probability intervals, we need to do Bayesian statistics
• Bayesian statistics produce credibility intervals that are interpreted as probability intervals
• A Bayesian analysis of our fictional flossing experiment returns a 95% credibility interval of [-0.06, 0.20]
• We can say, “there is a 95% probability that the difference in flossing proportions is between -0.06 and 0.20”
• Which method should we use? Which method is better? Which method gives the right answer?

Statistics is Subjective

Students are frequently left with the impression that statistics is a monolith, firmly rooted in mathematics and thus somehow pure — while in truth the subject is highly contentious and as closely aligned to philosophy as it is to math.

Darren Dahly, from his article Statistical Reform

What to do?

Wasserstein, et al. (2019) make the following recommendations:

• Accept uncertainty. Be thoughtful, open, and modest. (ATOM)
  • “Treat statistical results as being much more incomplete and uncertain than is currently the norm”
  • Evaluate research based on the importance of the questions being asked and the methods used to answer them, rather than the outcomes obtained
  • Describe all analyses performed and present all findings obtained, without regard to statistical significance
  • Understand and clearly express the limitations of your work

What about AI?

• It is amazing and worth exploring, but not always right.
• Anecdote: I asked ChatGPT the following on 2024-10-15:
  • I would like to run an experiment where I ask males and females whether or not they floss every day. How many subjects should I sample if I suspect the difference in proportion who answer yes is 0.1, and I want to have 0.90 power to reject the null hypothesis of no difference at a significance level of 0.05?
• It interpreted 0.1 as Cohen’s h, not a absolute difference in proportions, and gave a wrong answer. (I gave a bad prompt.) I tried to correct it:
  • It looks like you treated the effect size as Cohen’s h. Can you re-run the calculation assuming one group is 0.4 and the other group is 0.5, and thus the effect size of 0.1 is the difference in proportions?
• It told me: “It seems like I can’t do more advanced data analysis right now. Please try again later.” Link to exchange

How to do this in R

In case you’re interested…

power.prop.test(p1 = 0.4, p2 = 0.5, sig.level = 0.05, power = 0.9)

     Two-sample comparison of proportions power calculation 

              n = 518.0372
             p1 = 0.4
             p2 = 0.5
      sig.level = 0.05
          power = 0.9
    alternative = two.sided

NOTE: n is number in *each* group

AI Tips

• Ask for explanations of statistical concepts (“Explain the central limit theorem to me like I’m 5”)
• Ask it to interpret statistical output (copy and paste statistical results and ask for a one paragraph summary)
• Use it to help write scripts in Python or R (“Can you write a Python function to extract an email address from a sentence?”)
• ChatGPT gave me these tips: Link to exchange
• Remain skeptical and treat it like fire

Thanks

I hope you enjoyed.

Link to slides: https://is.gd/btpGXD

For help and advice with your data analysis: statlab@virginia.edu

References

Appendix: Bayesian analysis in R

Code to run Bayesian analysis of fictional flossing data

library(rstanarm)
# put data into data frame
d <- data.frame(yes = c(28, 35), 
                no = c(100-28, 100-35), 
                sex = factor(c("male", "female"), 
                             levels = c("male", "female")))
# fit logistic regression model
m <- stan_glm(cbind(yes, no) ~ sex, data = d, 
              family = binomial)
# extract the posterior sample
sims <- as.matrix(m)
# sample of male probabilities
male <- plogis(sims[,1])
# sample of female probabilities
female <- plogis(rowSums(sims))
# 95% credibility interval
quantile(female - male, probs = c(0.025, 0.975))