Basic Data Analysis with R

Part 1. Know Your Data: What is YRBSS?

Before you analyse anything, understand what you have

Load packages, read the data

library(tidyverse)   # data wrangling + ggplot2
library(gtsummary)   # publication-ready tables
library(janitor)     # clean_names(), tabyl()

yrbss <- read_csv("yrbss.csv")

Why this order matters:

• tidyverse first - gives you read_csv(), mutate(), filter(), plus pipes
• gtsummary next - for the descriptive table you’ll build later
• janitor last - small helpers; loaded last so its tabyl() doesn’t get masked

Inspect the structure first

glimpse(yrbss)

Rows: 20,000
Columns: 11
$ record      <dbl> 931897, 333862, 36253, 1095530, 1303997, 261619, 926649, 1…
$ age         <fct> 15 years old, 17 years old, 18 years old or older, 15 year…
$ sex         <fct> Female, Female, Male, Male, Male, Male, Male, Male, Male, …
$ grade       <fct> 10th, 12th, 11th, 10th, 9th, 9th, 11th, 12th, 12th, 10th, …
$ race4       <fct> White, White, Hispanic/Latino, Black or African American, …
$ race7       <chr> "White", "White", "Hispanic/Latino", "Black or African Ame…
$ bmi         <dbl> 17.1790, 20.2487, NA, 27.9935, 24.4922, NA, 20.5435, 19.25…
$ stweight    <dbl> 54.43, 57.15, NA, 85.73, 66.68, NA, 70.31, 58.97, 123.38, …
$ obese       <fct> Not obese, Not obese, NA, Not obese, Not obese, NA, Not ob…
$ bmi_cat     <fct> Underweight, Normal, NA, Overweight, Normal, NA, Normal, N…
$ bmi_missing <lgl> FALSE, FALSE, TRUE, FALSE, FALSE, TRUE, FALSE, FALSE, FALS…

Why we always do this first:

glimpse() is the analyst’s first move on any new dataset. In one line you learn:

• How many rows and columns is this what you expected?
• Variable names and types <dbl>, <chr>, <fct>?
• The first few values do they look sensible?

Catching a wrongly-typed variable here (e.g. age stored as text) saves hours of debugging later.

Audit missing data - never skip this

colSums(is.na(yrbss))

     record         age         sex       grade       race4       race7 
          0         205         231         406         710        1187 
        bmi    stweight       obese     bmi_cat bmi_missing 
       6458        6557        6458        6458           0 

What this tells you:

• bmi and stweight each have ~6,400 missing (~32%) - these are students who didn’t report height/weight
• age, sex, grade, race4 are <2% missing - safe for descriptive tables

Variable types determine your analysis

yrbss |>
  summarise(across(everything(), class)) 

# A tibble: 1 × 11
  record  age   sex   grade race4 race7 bmi   stweight obese bmi_cat bmi_missing
  <chr>   <chr> <chr> <chr> <chr> <chr> <chr> <chr>    <chr> <chr>   <chr>      
1 numeric fact… fact… fact… fact… char… nume… numeric  fact… factor  logical    

The four R types you need to recognise:

Why type matters: gtsummary, ggplot2, and every statistical test treat numeric and categorical variables differently. If grade is stored as <chr> instead of <fct>, your tables sort alphabetically ("10th" before "9th") - wrong.

Descriptive Statistics: The four dimensions of any variable

Summarise every variable before you model anything

Describe all variables before any modelling. A one-line summary per variable prevents major errors downstream.

What gtsummary does for you

• Auto-detects variable types: continuous vs categorical
• Continuous → median (IQR) by default · Categorical → n (%)
• add_p() automatically picks the correct statistical test
• add_overall() adds an all-participants column
• Output-ready for HTML, Word, and PDF via as_gt() or as_flex_table()
• Bold or italicise labels in one pipe step

Stratified tbl_summary() code

yrbss |>
  select(sex, grade, race4, bmi) |>
  tbl_summary(by = sex) |>     # split columns by sex
  add_p() |>                   # auto-pick the right test
  add_overall() |>             # add an "Overall" column
  bold_labels() |>
  italicize_levels()

The by = argument is the magic. It splits the table into one column per group and - combined with add_p() - automatically chooses:

• Wilcoxon for skewed continuous (like BMI)
• Chi-square for categorical
• Fisher’s exact if any expected cell < 5

Customising labels and statistics

yrbss |>  select(sex, bmi, stweight) |>
  tbl_summary( by = sex,
    statistic = list(
      all_continuous()  ~ "{median} ({p25}, {p75})",
      all_categorical() ~ "{n} ({p}%)"),
    digits = list(all_continuous() ~ 1),
    label  = list(
      bmi      ~ "Body Mass Index (kg/m²)",
      stweight ~ "Body weight (kg)"),
    missing = "ifany", missing_text = "Missing, n") |>
  add_p(test = list(all_continuous() ~ "wilcox.test")) |>
  add_overall() |>  bold_labels() |>
  modify_caption("**Table 1. Participant characteristics**")

How to read a gtsummary table

Exploratory Data Analysis

• Why we plot every variable first?

  • Visualise before you conclude

“Exploratory data analysis can never be the whole story, but nothing else can serve as the foundation stone.”
- John Tukey, 1977

BMI distribution (histogram)

med_bmi <- median(yrbss$bmi, na.rm = TRUE)
ggplot(yrbss, aes(x = bmi)) +
  geom_histogram(binwidth = 2, fill = "#0D9488",
                 colour = "white") +
  geom_vline(xintercept = med_bmi,
             colour = "#EF4444", linewidth = 1,
             linetype = "dashed") +
  annotate("text", x = med_bmi + 1, y = 1700,
           label = paste("Median =",
                         round(med_bmi, 1)),
           colour = "#EF4444", hjust = 0) +
  labs(title = "Distribution of BMI",
       x = "BMI (kg/m²)", y = "Count")

What to look for in any histogram:

• Centre: where is the bulk of the data?
• Spread: narrow peak or wide spread?
• Shape: symmetric (bell), skewed (long tail), or bimodal (two peaks)?
• Outliers: values far from the bulk that need investigation

For BMI: clear right-skew (long tail to the right) → use median, not mean.

BMI density curves by sex

ggplot(yrbss |> filter(!is.na(sex)),
       aes(x = bmi, fill = sex)) +
  geom_density(alpha = 0.45, colour = NA) +
  scale_fill_manual( values = c(Female = "#1B3A6B",
               Male   = "#F59E0B")) +
  labs(title = "BMI density by sex",
       x = "BMI (kg/m²)", y = "Density", fill = NULL) +
  theme(legend.position = "top")

Reading two density curves side-by-side:

A density plot is a smoothed histogram - the area under each curve sums to 1. So you can compare two groups of different sizes on the same axes.

• Both curves are right-skewed - same shape, no surprises there
• Male curve sits to the right - males have systematically higher BMI
• Similar widths → similar IQRs (similar spread)

This is what a single number like “mean BMI” would hide entirely.

Box plots: comparing groups

ggplot(yrbss |> filter(!is.na(grade)),
       aes(x = grade, y = bmi, fill = grade)) +
  geom_boxplot() +
  scale_fill_manual(values = c("#1B3A6B", "#0D9488","#F59E0B","#EF4444" )) +
  labs(title = "BMI by Grade", x = NULL, 
       y = "BMI (kg/m²)") +
  theme(legend.position = "none")

Anatomy of a box plot:

• Box = IQR (25th-75th percentile)
• Line in box = median
• Whiskers = 1.5 × IQR
• Points beyond whiskers = outlier candidates

Stacked bars: composition within groups

yrbss |> filter(!is.na(grade), !is.na(sex)) |>
  count(grade, sex) |>
  group_by(grade) |>
  mutate(pct = n / sum(n) * 100) |>
  ggplot(aes(x = grade, y = pct, fill = sex)) +
  geom_col(position = "fill") +
  geom_text(aes(label = sprintf("%.0f%%", pct)),
    position = position_fill(vjust = 0.5),
    colour = "white", fontface = "bold", size = 3.5) +
  scale_fill_manual(values = c("#1B3A6B","#0D9488")) +
  scale_y_continuous(labels = scales::percent) +
  labs(title = "Sex distribution by grade",
       subtitle = "Each bar sums to 100%",
       x = "Grade", y = "Proportion", fill = "Sex")

Stacked vs dodged: pick the right one:

• Stacked (position = "fill") - best for proportions within a group; each bar sums to 100%.
• Dodged (position = "dodge") - best for absolute counts between groups; bars side-by-side.

Use stacked when relative composition matters (here: sex balance per grade). Use dodged when raw counts matter.

What to Report

Mean vs median

Table or graph: when to use which

Quick-reference reporting card

Linear regression: BMI predicted by weight

model1 <- lm(bmi ~ stweight, data = yrbss)
summary(model1)

Call:
lm(formula = bmi ~ stweight, data = yrbss)

Residuals:
     Min       1Q   Median       3Q      Max 
-11.4075  -1.6816  -0.2079   1.4228  24.4669 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 6.560047   0.092984   70.55   <2e-16 ***
stweight    0.251396   0.001336  188.21   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.611 on 13439 degrees of freedom
  (6559 observations deleted due to missingness)
Multiple R-squared:  0.725, Adjusted R-squared:  0.7249 
F-statistic: 3.542e+04 on 1 and 13439 DF,  p-value: < 2.2e-16

Tidy regression output with broom

broom::tidy(model1, conf.int = TRUE)

# A tibble: 2 × 7
  term        estimate std.error statistic p.value conf.low conf.high
  <chr>          <dbl>     <dbl>     <dbl>   <dbl>    <dbl>     <dbl>
1 (Intercept)    6.56    0.0930       70.6       0    6.38      6.74 
2 stweight       0.251   0.00134     188.        0    0.249     0.254

broom::glance(model1)

# A tibble: 1 × 12
  r.squared adj.r.squared sigma statistic p.value    df  logLik    AIC    BIC
      <dbl>         <dbl> <dbl>     <dbl>   <dbl> <dbl>   <dbl>  <dbl>  <dbl>
1     0.725         0.725  2.61    35421.       0     1 -31973. 63951. 63974.
# ℹ 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>

broom turns model output into a tidy data frame easy to filter, format, and pass to gtsummary or ggplot2.

• tidy() - one row per coefficient (estimate, SE, p-value, CI)
• glance() - one row per model (R², AIC, n)
• augment() - adds fitted values + residuals to your data

Multiple regression: adding sex and age

model2 <- lm(bmi ~ stweight + sex + race4, data = yrbss)
summary(model2)

Call:
lm(formula = bmi ~ stweight + sex + race4, data = yrbss)

Residuals:
     Min       1Q   Median       3Q      Max 
-13.7772  -1.3761  -0.1224   1.1726  25.1075 

Coefficients:
                                Estimate Std. Error t value Pr(>|t|)    
(Intercept)                     6.115324   0.086190  70.951   <2e-16 ***
stweight                        0.279391   0.001225 228.160   <2e-16 ***
sexMale                        -2.809706   0.041274 -68.075   <2e-16 ***
race4Black or African American  0.061472   0.058664   1.048    0.295    
race4Hispanic/Latino            0.557191   0.056076   9.936   <2e-16 ***
race4White                     -0.491424   0.052357  -9.386   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.214 on 13022 degrees of freedom
  (6972 observations deleted due to missingness)
Multiple R-squared:  0.8025,    Adjusted R-squared:  0.8024 
F-statistic: 1.058e+04 on 5 and 13022 DF,  p-value: < 2.2e-16

Publication table for the regression

gtsummary::tbl_regression(model2,
  label = list(stweight ~ "Body weight (kg)",
    sex      ~ "Sex", race4      ~ "Ethinicity")) |>
  bold_labels() |>
  add_glance_source_note(
    include = c(r.squared, adj.r.squared, AIC, nobs))

Publication table for the regression

Logistic regression: when the outcome is yes / no

So far our outcome (bmi) has been continuous. But many real questions are about a binary outcome:

• Is the student overweight? (yes / no)
• Did the patient survive? (yes / no)
• Did the household adopt the intervention? (yes / no)

Why we can’t just use linear regression here:

• A linear model can predict probabilities below 0 or above 1 - meaningless
• The relationship between predictors and probability is S-shaped, not straight
• The variance of a 0/1 outcome isn’t constant - violates the equal-variance assumption.

The fix: transform the probability so the linear part of the model has no upper or lower bound. That transformation is called the logit.

The logit link - modelling probabilities

We don’t model the probability \(p\) directly. We model its log-odds:

\[\text{logit}(p) = \log\!\left(\frac{p}{1-p}\right) = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \dots\]

Why this works

• p / (1 - p) converts probability into odds (ranges 0 to ∞)
• Taking log() makes it range from −∞ to +∞
• Now the right-hand side can take any value - no more impossible predictions
• Inverting the logit gives back a probability between 0 and 1

In R: glm(y ~ x, family = binomial) fits this model. The default link is logit. You almost never need to change it.

yrbss <- yrbss |>
  mutate(overweight = factor(ifelse(bmi >= 25, "Yes", "No")))

# Then model:
glm(overweight ~ age + sex + grade + race4,
    data = yrbss,
    family = binomial)

Call:  glm(formula = overweight ~ age + sex + grade + race4, family = binomial, 
    data = yrbss)

Coefficients:
                   (Intercept)                 age13 years old  
                      -0.97364                        -1.08840  
               age14 years old                 age15 years old  
                      -0.58921                        -0.33808  
               age16 years old                 age17 years old  
                      -0.13499                         0.06457  
      age18 years old or older                         sexMale  
                       0.29058                         0.21353  
                     grade10th                       grade11th  
                      -0.09245                        -0.11765  
                     grade12th  race4Black or African American  
                      -0.17767                         0.37905  
          race4Hispanic/Latino                      race4White  
                       0.32089                        -0.07208  

Degrees of Freedom: 12991 Total (i.e. Null);  12978 Residual
  (7008 observations deleted due to missingness)
Null Deviance:      15540 
Residual Deviance: 15280    AIC: 15310

Odds and odds ratios - the interpretation key

Odds = probability of event / probability of no event.

Odds ratio (OR) = odds in group A / odds in group B.

Step 2 - Fit the model

model_log <- glm(overweight ~ age + sex + grade + race4,
                 data   = yrbss, family = binomial)
summary(model_log)

Call:
glm(formula = overweight ~ age + sex + grade + race4, family = binomial, 
    data = yrbss)

Coefficients:
                               Estimate Std. Error z value Pr(>|z|)    
(Intercept)                    -0.97364    0.50687  -1.921   0.0547 .  
age13 years old                -1.08840    0.66961  -1.625   0.1041    
age14 years old                -0.58921    0.50904  -1.157   0.2471    
age15 years old                -0.33808    0.50531  -0.669   0.5035    
age16 years old                -0.13499    0.50453  -0.268   0.7890    
age17 years old                 0.06457    0.50466   0.128   0.8982    
age18 years old or older        0.29058    0.50613   0.574   0.5659    
sexMale                         0.21353    0.03941   5.419 6.01e-08 ***
grade10th                      -0.09245    0.06773  -1.365   0.1723    
grade11th                      -0.11765    0.08165  -1.441   0.1496    
grade12th                      -0.17767    0.09343  -1.902   0.0572 .  
race4Black or African American  0.37905    0.05855   6.474 9.52e-11 ***
race4Hispanic/Latino            0.32089    0.05651   5.679 1.36e-08 ***
race4White                     -0.07208    0.05488  -1.313   0.1890    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 15542  on 12991  degrees of freedom
Residual deviance: 15278  on 12978  degrees of freedom
  (7008 observations deleted due to missingness)
AIC: 15306

Number of Fisher Scoring iterations: 4

Reading the raw summary() output

The output has four parts you should always check:

Step 3 - Get odds ratios with broom

broom::tidy() does the exponentiation, the CIs, and tidies everything into a data frame in one call:

broom::tidy(model_log, exponentiate = TRUE,   # log-odds → ORs
            conf.int     = TRUE,   # add 95% CI
            conf.level   = 0.95)

# A tibble: 14 × 7
   term                  estimate std.error statistic p.value conf.low conf.high
   <chr>                    <dbl>     <dbl>     <dbl>   <dbl>    <dbl>     <dbl>
 1 (Intercept)              0.378     0.507    -1.92    0.055    0.13      0.986
 2 age13 years old          0.337     0.67     -1.62    0.104    0.088     1.27 
 3 age14 years old          0.555     0.509    -1.16    0.247    0.211     1.62 
 4 age15 years old          0.713     0.505    -0.669   0.503    0.274     2.07 
 5 age16 years old          0.874     0.505    -0.268   0.789    0.336     2.53 
 6 age17 years old          1.07      0.505     0.128   0.898    0.41      3.09 
 7 age18 years old or o…    1.34      0.506     0.574   0.566    0.513     3.88 
 8 sexMale                  1.24      0.039     5.42    0        1.15      1.34 
 9 grade10th                0.912     0.068    -1.36    0.172    0.798     1.04 
10 grade11th                0.889     0.082    -1.44    0.15     0.758     1.04 
11 grade12th                0.837     0.093    -1.90    0.057    0.697     1.01 
12 race4Black or Africa…    1.46      0.059     6.47    0        1.30      1.64 
13 race4Hispanic/Latino     1.38      0.057     5.68    0        1.23      1.54 
14 race4White               0.93      0.055    -1.31    0.189    0.836     1.04 

Now the estimate column is the OR, and conf.low / conf.high give the 95 % CI - ready to paste into a manuscript.

Publication table with gtsummary