I am using 2025 NFL team statistics from nflverse. Each row shows one team’s stats for one game. The variables I will use are team, opponent, passing yards, and passing touchdowns. Team and opponent are categorical variables. Passing yards and passing touchdowns are numerical variables. Passing yards is how many yards a team gains through passing, and passing touchdowns is how many touchdowns it scores through passes. I want to see if teams with more passing yards also have more passing touchdowns.
# Load tidyverse for reading data and making graphslibrary(tidyverse)
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr 1.2.1 ✔ readr 2.2.0
✔ forcats 1.0.1 ✔ stringr 1.6.0
✔ ggplot2 4.0.3 ✔ tibble 3.3.1
✔ lubridate 1.9.5 ✔ tidyr 1.3.2
✔ purrr 1.2.2
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
# Read the 2025 NFL team statistics CSVnfl_data <- readr::read_csv("https://github.com/nflverse/nflverse-data/releases/download/stats_team/stats_team_week_2025.csv")
Rows: 570 Columns: 138
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (7): team, season_type, game_id, opponent_team, fg_made_list, fg_misse...
dbl (131): season, week, completions, attempts, passing_yards, passing_tds, ...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
# Show the first six rows of the datasethead(nfl_data)
# Keep only regular-season gamesnfl_regular <- nfl_data %>%filter(season_type =="REG")
# Keep the four variables needed for this projectnfl_clean <- nfl_regular %>%select(team, opponent_team, passing_yards, passing_tds)# Preview the selected columnshead(nfl_clean)
# A tibble: 6 × 4
team opponent_team passing_yards passing_tds
<chr> <chr> <dbl> <dbl>
1 ARI NO 163 2
2 ATL TB 298 1
3 BAL BUF 209 2
4 BUF BAL 394 2
5 CAR JAX 154 1
6 CHI MIN 210 1
# Count missing values in the four selected columnssum(is.na(nfl_clean))
[1] 0
I kept only regular-season games and selected the four variables needed for my project. I checked these columns for missing values and found none, so I did not need to remove any rows for missing data.
Linear Regression
# Model passing touchdowns using passing yardsfit1 <-lm(passing_tds ~ passing_yards, data = nfl_clean)# Show the regression resultssummary(fit1)
Call:
lm(formula = passing_tds ~ passing_yards, data = nfl_clean)
Residuals:
Min 1Q Median 3Q Max
-2.4797 -0.7243 -0.0973 0.6305 3.8443
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.2617666 0.1450744 -1.804 0.0717 .
passing_yards 0.0077884 0.0006147 12.671 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.021 on 542 degrees of freedom
Multiple R-squared: 0.2285, Adjusted R-squared: 0.2271
F-statistic: 160.5 on 1 and 542 DF, p-value: < 2.2e-16
The positive slope means that teams with more passing yards tend to have more passing touchdowns. For every additional 100 passing yards, the model predicts about 0.78 more passing touchdowns on average.
The p-value for passing yards is less than 2e-16, which is below 0.05. This suggests a statistically significant relationship between passing yards and passing touchdowns. The adjusted R-squared is 0.2271, or about 22.7%. This means the model explains about 22.7% of the variation in passing touchdowns, so passing yards alone does not explain most of the differences between games.
Diagnostic Plots
# Show diagnostic plots to check the regression modelplot(fit1)
The Residuals vs Fitted plot has diagonal bands because touchdowns are whole numbers. Its red line stays near zero but bends slightly. In the Q-Q plot, most points follow the line in the middle but move away at the ends. The Scale-Location plot suggests that the spread of the errors is not completely constant. Most points have low leverage, while a few games have higher leverage. These plots show limitations of the linear model, so its predictions should be treated cautiously.
Visualization
# Keep the Ravens, Bills, and Chiefs for the graphnfl_teams <- nfl_clean %>%filter(team %in%c("BAL", "BUF", "KC"))
# Compare passing yards and touchdowns for three teamsggplot(nfl_teams, aes(x = passing_yards, y = passing_tds, color = team)) +geom_point() +labs(title ="Passing Yards and Touchdowns in 2025",x ="Passing Yards per Game",y ="Passing Touchdowns per Game",color ="Team",caption ="Source: nflverse, 2025 regular season") +scale_color_manual(values =c("BAL"="purple", "BUF"="blue", "KC"="red")) +theme_bw()
Conclusion
I cleaned the data by keeping only regular-season games and selecting team, opponent, passing yards, and passing touchdowns. There were no missing values in these four columns. I used all teams for the regression, but narrowed the graph to the Ravens (BAL), Bills (BUF), and Chiefs (KC) so it was easier to read.
Each dot shows one team’s stats for one game. Purple represents the Ravens, blue represents the Bills, and red represents the Chiefs. The graph shows that similar passing yard totals can come with different numbers of touchdowns. One Chiefs game had over 350 passing yards but no passing touchdowns. The regression also showed a positive relationship, but passing yards alone explained only about 23% of the variation in passing touchdowns.
I would have liked to include rushing touchdowns to get a fuller picture of each team’s offense. For this project, I focused on passing yards and passing touchdowns.
