Soccer Code Logic.Rmd

---
title: "Soccer Regression Project"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

# We are going to start out with some EDA

I started by narrowing down the dataset to quantitative, per-game variables to prevent bias in the number of games a team played during the season.
```{r}
library(tidyverse)

womens_ncaa_soccer_20182019 <- read.csv("womens_ncaa_soccer_20182019.csv")

soccer_qunant <- womens_ncaa_soccer_20182019 %>% 
    mutate(goal_diff = goals - ga) %>% 
  select(team,
         team_games,
         corners_gp,
         assists_gp,
         fouls_gp,
         pk_pct,
         points_gp,
         save_pct,
         saves_gp,
         gpg,
         sog_pct,
         win_pct,
         sog_gp,
         goal_diff,
         season)
```

Constructing a correlation matrix of relevant per-game variables to determine the strength of relationships between the variables in the dataset.
```{r}

# install.packages("gt")
library(gt)
# install.packages("xtable")
library(xtable)

corr_table <- soccer_qunant %>%
  select(goal_diff,
         fouls_gp,
         pk_pct,
         save_pct,
         saves_gp,
         corners_gp,
         sog_pct,
         sog_gp) %>% 
  cor()
corr_table <- round(corr_table, 2)
corr_table[upper.tri(corr_table)] <- ""
print(xtable(corr_table), type = "html")
```


In this matrix I can explore which variables are strongly related and which might not be as relevant to my regression.

Looking at pairs plots for some strongly related variables.
```{r}
womens_ncaa_soccer_20182019 %>% 
  select(win_pct, 
         points_gp,
         shots_gp,
         saves_gp) %>% 
  pairs()

```
All of these variables look to be linearly related to one another, with varying levels of strength.

Exploring scatterplots of certain variables with win percent.
```{r}
womens_ncaa_soccer_20182019 <- 
  womens_ncaa_soccer_20182019 %>% 
  mutate(goal_diff = goals - ga)

# install.packages("ggpubr")
library(ggpubr)

womens_ncaa_soccer_20182019 %>% 
  ggplot(aes(x  = saves_gp,
             y = goal_diff)) +
  geom_point() +
    theme_bw() +
  labs(title = "Goal Differential and Saves Per Game",
       subtitle = "Women's NCAA Soccer 2018 & 2019",
       x = "Saves Per Game",
       y = "Season Goal Differential") +
  stat_cor(method = "pearson",
           label.x = 6,
           label.y = 60)

womens_ncaa_soccer_20182019 %>% 
  ggplot(aes(x  = fouls_gp,
             y = goal_diff)) +
  geom_point() + 
    theme_bw() +
  labs(title = "Goal Differential and Fouls Per Game",
       subtitle = "Women's NCAA Soccer 2018 & 2019",
       x = "Fouls Per Game",
       y = "Season Goal Differential") +
  stat_cor(method = "pearson",
           label.x = 10,
           label.y = 70)

womens_ncaa_soccer_20182019 %>% 
  ggplot(aes(x  = corners_gp,
             y = goal_diff)) +
  geom_point() + 
    theme_bw() +
  labs(title = "Goal Differential and Corners Per Game",
       subtitle = "Women's NCAA Soccer 2019",
       x = "Corners Per Game",
       y = "Season Goal Differential",
       caption = "corner data only available only for 2019 season") +
  stat_cor(method = "pearson",
           label.x = 4.5,
           label.y = -60)

womens_ncaa_soccer_20182019 %>% 
  ggplot(aes(x  = pk_pct,
             y = goal_diff)) +
  geom_jitter(height = 0) + 
    theme_bw() +
  labs(title = "Goal Differential and Penalty Kick Percentage",
       subtitle = "Women's NCAA Soccer 2018 & 2019",
       x = "Penalty Kick Percentage",
       y = "Season Goal Differential") +
  stat_cor(method = "pearson",
           label.x = 0,
           label.y = 70)

womens_ncaa_soccer_20182019 %>% 
  ggplot(aes(x  = save_pct,
             y = goal_diff)) +
  geom_point() + 
    theme_bw() +
  labs(title = "Goal Differential and Save Percentage",
       subtitle = "Women's NCAA Soccer 2018 & 2019",
       x = "Save Percentage",
       y = "Season Goal Differential") +
  stat_cor(method = "pearson",
           label.x = 0.75,
           label.y = -60)

```
It looks like points per game, shots per game, and assists per game have a pretty strong linear relationship with win percent. Saves per game has a moderate linear relationship with win percent. Fouls per game does not appear to have any sort of relevant relationship with win percent. Assists per game and points per game are clearly very strongly correlated. Corners per game and assists per game seem to have a weak relationship.

DZ: A more consice way to do this is to use `facet_wrap`:

```{r}
womens_ncaa_soccer_20182019 %>% 
  select(goal_diff, saves_gp, fouls_gp, corners_gp, pk_pct, save_pct) %>% 
  gather("predictor", "value", -goal_diff) %>% 
  ggplot(aes(x  = value,
             y = goal_diff)) +
  geom_point() +
  theme_bw() +
  facet_wrap(~predictor, scales = "free") +
  labs(title = "Goal Differential and Saves Per Game",
       subtitle = "Women's NCAA Soccer 2018 & 2019",
       x = "Saves Per Game",
       y = "Season Goal Differential") +
  stat_cor(method = "pearson")
```

Exploring the effect of season
```{r}
womens_ncaa_soccer_20182019 %>% 
  ggplot(aes(x = pk_pct,
             y = win_pct,
             color = as.factor(season))) +
  geom_point()
```
There does not appear to be a difference between the seasons in the relationships between penalty kick percentage and win percent.


Looking at the distributions of some variables that were having interesting effects.
```{r}
hist(womens_ncaa_soccer_20182019$fouls)

hist(womens_ncaa_soccer_20182019$psatt)
```
Fouls appears to be pretty normally distributed, while penalty shot attempts is right skewed with many teams having zero penalty kick attempts.

# Starting to fit a potential regression

I started by including the per-game predictors that had a correlation of more than 0.7 or less than -0.7 with win percent.
```{r}
init_fit <- lm(win_pct ~ assists_gp + gaa + 
                 points_gp + gpg + shots_gp + 
                 sog_gp,
               data = womens_ncaa_soccer_20182019)
summary(init_fit)
```
With this initial model, we can see that coefficient estimates are very small and many p-values are very high.


I plotted the residuals against some of the predictor variables to try and determine if there were clear interactions among these variables that need to be accounted for in the model.
```{r}
plot(init_fit$residuals~
       womens_ncaa_soccer_20182019$assists_gp, 
     color = womens_ncaa_soccer_20182019$season)
plot(init_fit$residuals~
       womens_ncaa_soccer_20182019$points_gp)
plot(init_fit$residuals~
       womens_ncaa_soccer_20182019$sog_gp)
plot(init_fit$residuals~
       womens_ncaa_soccer_20182019$shots_gp)
```
None of these plots show evidence of an interaction being necessary. There do seem to be some outliers, but the strong majority of the data is randomly scattered.

I began to narrow my model by removing the predictor with the highest p-value. After Monday's lecture, I realize that the model building method I used is very formulaic and doesn't create the best possible fit for the data. 
```{r}
second_fit <- lm(win_pct ~ assists_gp + gaa + 
                   points_gp + shots_gp + sog_gp,
               data = womens_ncaa_soccer_20182019)
summary(second_fit)
```

I again removed the variable with the highest p-value from the model.
```{r}
third_fit <- lm(win_pct ~ assists_gp + gaa + 
                  points_gp + sog_gp,
                 data = womens_ncaa_soccer_20182019)
summary(third_fit)
```

After talking with Caleb, I took the variable about goals against out of the model because it is a very obvious predictor of win percent.
```{r}
fourth_fit <- lm(win_pct ~ assists_gp + points_gp +
                   sog_gp,
                 data = womens_ncaa_soccer_20182019)
summary(fourth_fit)
```
Removing goals against caused a large drop in the Mutiple R-Squared for the model and caused the assists per game variable to lose its significance.

Looking at some residuals by variable from this fourth fit.
```{r}
plot(fourth_fit$residuals~
       womens_ncaa_soccer_20182019$assists_gp)
plot(fourth_fit$residuals~
       womens_ncaa_soccer_20182019$points_gp)
plot(fourth_fit$residuals~
       womens_ncaa_soccer_20182019$sog_gp)
```
We are seeing again that there is no clear evidence that an interaction is necessary for this model.

Even though the residual plots did not indicate presence of an interaction, I wanted to explore the reactions that may be present
```{r}
experimental_fit <- lm(win_pct ~ assists_gp + 
                         points_gp + sog_gp + 
                         assists_gp * points_gp,
                 data = womens_ncaa_soccer_20182019)
summary(experimental_fit)
plot(experimental_fit)
```
The interaction term between assists per game and points per game has a very small p-value. The residual plots for this model look very good.

Exploring a full interaction model to see if other interactions may be significant.
```{r}
interaction_fit <- lm(win_pct ~
                      assists_gp*points_gp*sog_gp, 
                      data = womens_ncaa_soccer_20182019)
summary(interaction_fit)
```
Some of the interactions between variables did have significance, but they did not contribute to the model as  a whole very much and they could be accounting for too much random noise in the data.


# Testing models

I'm setting up training and testing sets for my model. We are using 2018 data as the training set and 2019 for the testing set. This is logical because the 2018 data could have been useful before the 2019 season to predict.
```{r}
womens_ncaa_soccer_2018 <- womens_ncaa_soccer_20182019 %>% 
  filter(season == 2018)

womens_ncaa_soccer_2019 <- womens_ncaa_soccer_20182019 %>% 
  filter(season == 2019)
```

I'm fitting two candidate models on the 2018 season data
```{r}
candidate_model_1 <- lm(win_pct ~ assists_gp + points_gp +
                          sog_gp,
                data = womens_ncaa_soccer_2018)
summary(candidate_model_1)
plot(candidate_model_1)

candidate_model_2 <- lm(win_pct ~ assists_gp + points_gp +
                          sog_gp + points_gp * assists_gp,
                        data = womens_ncaa_soccer_2018)
summary(candidate_model_2)
plot(candidate_model_2)
```

Calculating the MSE for the candidate models.
```{r}
model_1_preds <- predict(candidate_model_1, 
                         newdata = womens_ncaa_soccer_2019)
model_1_mse <- mean((model_1_preds -
                       womens_ncaa_soccer_2019$win_pct)^2)
model_1_mse

model_2_preds <- predict(candidate_model_2,
                         newdata = womens_ncaa_soccer_2019)
model_2_mse <- mean((model_2_preds - 
                       womens_ncaa_soccer_2019$win_pct)^2)
model_2_mse

```
The MSE for candidate model 2 (with the interaction included) has a slightly smaller MSE of the two models. Because win percent has values between 0 and 1, these small values do not necessarily indicate an incredibly strong predictive value. The MSE needs to be compared among models.

# Rethinking things with Caleb

Over the weekend Caleb and I worked independently of one another to work with the data and begin fitting regression models. Caleb mutated the dataset to make a new variable, goal differential.
```{r}
womens_ncaa_soccer_20182019 <- 
  womens_ncaa_soccer_20182019 %>% 
  mutate(goal_diff = goals - ga)
```
This variable has units that are more easily interpretable, and it is very highly correlated with win percent. It encompasses multiple columns of data as well.

#Doing some EDA with goal differential

Exploring the distribution of goal differential
```{r}
library(rvest)

soccer_rankings <- read_html("https://www.ncaa.com/rankings/soccer-women/d1/ncaa-womens-soccer-rpi") %>% 
  html_table()

school_conf <- soccer_rankings[[1]] %>% 
  select("School", "Conference")

womens_ncaa_soccer_20182019 <- 
  womens_ncaa_soccer_20182019 %>% 
  left_join(school_conf, by = c("team" = "School"))

womens_ncaa_soccer_20182019 %>% 
  ggplot(aes(x = goal_diff)) +
  geom_histogram(binwidth = 10, 
                 color = "black",
                 fill = "#CC002B") +
  theme_bw() +
  labs(title = "Season Goal Differential",
       subtitle = "Women's NCAA Soccer 2018 & 2019",
       x = "Goal Differential",
       y = "count")
```
We can see that goal differential is pretty normally distributed

Exploring goal ratio as a possible predictor?
```{r}
womens_ncaa_soccer_20182019 <- 
  womens_ncaa_soccer_20182019 %>% 
  mutate(goal_ratio = goals / ga)

womens_ncaa_soccer_20182019 %>% 
  ggplot(aes(x = team_games,
             y = goal_ratio)) +
  geom_point()
```
Very confusing interpretation and still has a relationship with number of games played.

I'm looking at scatterplots and correlation values for goal differential and several per-game predictor variables.
```{r}
womens_ncaa_soccer_20182019 %>% 
  select(goal_diff, assists_gp, save_pct, corners_gp,
         fouls_gp, pk_pct, sog_gp) %>% 
  pairs()

womens_ncaa_soccer_20182019 %>% 
  select(goal_diff, assists_gp, save_pct, corners_gp,
         fouls_gp, pk_pct, sog_gp) %>% 
  cor()

womens_ncaa_soccer_20182019 %>% 
  ggplot(aes(x = assists_gp,
             y = goal_diff)) + 
  geom_point()

womens_ncaa_soccer_20182019 %>% 
  ggplot(aes(x = save_pct,
             y = goal_diff)) + 
  geom_point()

womens_ncaa_soccer_20182019 %>% 
  ggplot(aes(x = corners_gp,
             y = goal_diff)) + 
  geom_point()

womens_ncaa_soccer_20182019 %>% 
  ggplot(aes(x = fouls_gp,
             y = goal_diff)) + 
  geom_point()

womens_ncaa_soccer_20182019 %>% 
  ggplot(aes(x = pk_pct,
             y = goal_diff)) + 
  geom_point()

womens_ncaa_soccer_20182019 %>% 
  ggplot(aes(x = sog_gp,
             y = goal_diff)) + 
  geom_point()
```





#Models with goal differential as the response

This is the model that Caleb came to over the weekend. I think using corners per game as a predictor could be a problem because it is missing from 2018 data so we wouldn't be able to calculate the MSe between 2018 and 2019. However, from a practical standpoint, we can expect that future season data would include corner data so it would be useful in the future.
```{r}
caleb_fit <- lm(goal_diff ~ corners_gp + save_pct + 
                  pk_pct + assists_gp + sog_gp,
                data = womens_ncaa_soccer_20182019)
summary(caleb_fit)
plot(caleb_fit)
```
In the residual plots, we can see that this model is a pretty great fit but there are a few outliers that are strongly influencing the data.

We tried removing the two datapoints that were consistent outliers in the model from the dataset to see what happened.
```{r}
nooutlier <- womens_ncaa_soccer_20182019 %>% 
  filter(team != "Alcorn") %>% 
  filter(team != "Chicago St.")
nooutlierfit <- lm(goal_diff ~ corners_gp + save_pct +
                     pk_pct + assists_gp + sog_gp,
                   data = nooutlier)
summary(nooutlierfit)
plot(nooutlierfit)
```
Removing these outliers only changes the model very slightly, but the residual plots are more ideal. Both of these teams had a win percent of zero, and there are two other teams also with a win percent of zero. Arbitrarily cutting out data is obviously not ideal, but could there be justification for removing all the teams that had a win percent of zero?


#Testing data for goal differential model without corners

```{r}

candidate_model_3 <- lm(goal_diff ~ assists_gp + points_gp +
                          sog_gp,
                data = womens_ncaa_soccer_2018)
summary(candidate_model_3)

candidate_model_4 <- lm(goal_diff ~ assists_gp + points_gp +
                          sog_gp + points_gp * assists_gp,
                        data = womens_ncaa_soccer_2018)
summary(candidate_model_4)

                         

candidate_model_5 <- lm(goal_diff ~ save_pct + 
                  pk_pct  + sog_gp + saves_gp + fouls_gp,
                data = womens_ncaa_soccer_2018)
summary(candidate_model_5)

plot(candidate_model_5)
```
Model 5 has a stronger multiple R-squared value.

Calculating MSE for potential models with goal differential as response using 2019 season data as testing data.
```{r}
model_3_preds <- predict(candidate_model_3, 
                         newdata = womens_ncaa_soccer_2019)
model_3_mse <- mean((model_3_preds -
                       womens_ncaa_soccer_2019$goal_diff)^2)
model_3_mse

model_4_preds <- predict(candidate_model_4,
                         newdata = womens_ncaa_soccer_2019)
model_4_mse <- mean((model_4_preds - 
                       womens_ncaa_soccer_2019$goal_diff)^2)
model_4_mse

model_5_preds <- predict(candidate_model_5,
                         newdata = womens_ncaa_soccer_2019)
model_5_mse <- mean((model_5_preds - 
                       womens_ncaa_soccer_2019$goal_diff)^2)
model_5_mse

```
We can see that model 5 provides the best fit by quite a large margin.

I want to plot the two years data in side by side scatterplots with lines using the models with and without corners.
```{r}
model_1819 <- lm(goal_diff ~ save_pct + 
                  pk_pct + assists_gp + sog_gp,
                data = womens_ncaa_soccer_20182019)
model_19corners <- lm(goal_diff ~ corners_gp + save_pct + 
                  pk_pct + assists_gp + sog_gp,
                data = womens_ncaa_soccer_20182019)
preds_1819_18 <- predict(model_1819, womens_ncaa_soccer_2018)
preds_1819_19 <- predict(model_1819,
                         womens_ncaa_soccer_2019)
preds_19corners <- predict(model_19corners, womens_ncaa_soccer_2019)

womens_ncaa_soccer_20182019 %>% 
  filter(season == 2018) %>% 
  ggplot() +
  geom_point(aes(x = assists_gp, y = goal_diff)) +
  geom_smooth(aes(x = assists_gp, y = preds_1819_18)) 

womens_ncaa_soccer_20182019 %>% 
  filter(season == 2019) %>% 
  ggplot() +
  geom_point(aes(x = assists_gp, y = goal_diff)) +
  geom_smooth(aes(x = assists_gp, y = preds_1819_19)) +
  geom_smooth(aes(x = assists_gp, y = preds_19corners), color = "red")



```

#cross validation looking at corners for 2019 data

Starting by generating two subsets of the 2019 data
```{r}
n_teams <- nrow(womens_ncaa_soccer_2019)

train_nums <- sample(n_teams, 
                     n_teams / 2, 
                     replace = FALSE)

test_nums <- (1:n_teams)[-train_nums]

train_2019 <- womens_ncaa_soccer_2019[train_nums,]

test_2019 <- womens_ncaa_soccer_2019[test_nums,]
```

Setting up candidate models
```{r}
candidate_model_2019_1 <- lm(goal_diff ~ assists_gp +
                               points_gp +
                               sog_gp,
                             data = train_2019)

candidate_model_2019_2 <- lm(goal_diff ~ assists_gp +
                               points_gp +
                               sog_gp + 
                               points_gp * assists_gp,
                             data = train_2019)

candidate_model_2019_3 <- lm(goal_diff ~ save_pct + 
                               pk_pct + assists_gp +
                               sog_gp,
                             data = train_2019)

candidate_model_2019_4 <- lm(goal_diff ~ save_pct + 
                               pk_pct + assists_gp +
                               sog_gp + corners_gp,
                             data = train_2019)
summary(candidate_model_2019_4)

candidate_model_2019_32 <- lm(goal_diff ~ save_pct + 
                               pk_pct + sog_pct  + saves_gp + fouls_gp +
                               sog_gp,
                             data = train_2019)

candidate_model_2019_42 <- lm(goal_diff ~ save_pct + 
                               pk_pct + sog_pct  + saves_gp + fouls_gp +
                               sog_gp + corners_gp,
                             data = train_2019)

model_5 <- lm(goal_diff ~ save_pct + pk_pct + assists_gp +
                sog_gp + corners_gp, 
              data = womens_ncaa_soccer_2019)
summary(model_5)

womens_ncaa_soccer_20182019 <- womens_ncaa_soccer_20182019 %>% 
  mutate(save_pct100 = save_pct * 100) %>% 
  mutate(pk_pct100 = pk_pct * 100)


womens_ncaa_soccer_2019 <- womens_ncaa_soccer_2019 %>% 
  mutate(save_pct100 = save_pct *100) %>% 
  mutate(pk_pct100 = pk_pct * 100)

final_model2019 <- lm(goal_diff~ pk_pct100 + save_pct100 + saves_gp + fouls_gp +
                    sog_gp + corners_gp, data = womens_ncaa_soccer_2019)
summary(final_model2019)

final_model <- lm(goal_diff~ pk_pct100 + save_pct100 + saves_gp + fouls_gp +
                    sog_gp, data = womens_ncaa_soccer_20182019)
summary(final_model)
# install.packages("GGally")
library(GGally)

ggcoef(final_model2019,
       exclude_intercept = TRUE,
       vline = TRUE,
       vline_color = "red") +
  theme_bw() +
  labs(title = "Coefficient Estimates",
       subtitle = "Women's NCAA Soccer 2019",
       caption = "corner data only available for 2019 season
       Adjusted R-Squared = 0.8754
       rMSE from five-fold Cross-Validation = 6.87",
       x = "Estimate",
       y = "Covariate")

ggcoef(final_model,
       exclude_intercept = TRUE,
       vline = TRUE,
       vline_color = "red") +
  theme_bw() +
  labs(title = "Coefficient Estimates",
       subtitle = "Women's NCAA Soccer 2018 & 2019",
       caption = "Adjusted R-Squared = 0.8824,
       rMSE with 2018 as Training and 2019 as Testing = 7.07",
       x = "Estimate",
       y = "Covariate")

```

Calculating MSE
```{r}
model_2019_1_preds <- predict(candidate_model_2019_1,
                              newdata = test_2019)
model_2019_1_mse <- mean((model_2019_1_preds - 
                            test_2019$goal_diff)^2)

model_2019_2_preds <- predict(candidate_model_2019_2,
                              newdata = test_2019)
model_2019_2_mse <- mean((model_2019_2_preds - 
                            test_2019$goal_diff)^2)

model_2019_3_preds <- predict(candidate_model_2019_3,
                              newdata = test_2019)
model_2019_3_mse <- mean((model_2019_3_preds - 
                            test_2019$goal_diff)^2)

model_2019_4_preds <- predict(candidate_model_2019_4,
                              newdata = test_2019)
model_2019_4_mse <- mean((model_2019_4_preds - 
                            test_2019$goal_diff)^2)

model_2019_32_preds <- predict(candidate_model_2019_32,
                              newdata = test_2019)
model_2019_32_mse <- mean((model_2019_32_preds - 
                            test_2019$goal_diff)^2)

model_2019_42_preds <- predict(candidate_model_2019_42,
                              newdata = test_2019)
model_2019_42_mse <- mean((model_2019_42_preds - 
                            test_2019$goal_diff)^2)

model_2019_1_mse
model_2019_2_mse
model_2019_3_mse
model_2019_4_mse

model_2019_32_mse
model_2019_42_mse

summary(candidate_model_2019_42)

plot(candidate_model_2019_42)
womens_ncaa_soccer_2019 <- womens_ncaa_soccer_2019 %>% 
  mutate(save_pct100 = save_pct *100)
candidate_model_422 <- lm(goal_diff ~  
                               pk_pct + sog_pct + save_pct100 + saves_gp + fouls_gp +
                               sog_gp + corners_gp,
                          data = womens_ncaa_soccer_2019)
summary(candidate_model_422)

sqrt(model_2019_32_mse)
sqrt(model_2019_42_mse)

confint(candidate_model_422)
```

```{r}
womens_ncaa_soccer_2018 <- 
  womens_ncaa_soccer_2018 %>% 
  mutate(save_pct100 = save_pct * 100)

final_model <- lm(goal_diff~ pk_pct + save_pct100 + saves_gp + fouls_gp +
                    sog_gp, data = womens_ncaa_soccer_2018)
summary(final_model)
final_model_preds <- predict(final_model, newdata = womens_ncaa_soccer_2019)
final_model_mse <- mean((final_model_preds - womens_ncaa_soccer_2019$goal_diff)^2)

final_model_mse
sqrt(final_model_mse)
```



```{r}
library(ggdendro)

soccer_correlation <- 
  womens_ncaa_soccer_20182019 %>% 
  na.omit() %>% 
  select(!"X") %>% 
  select(!"X.1") %>% 
  select(-season) %>% 
  select(!"save_pct100") %>% 
  select(!"goal_ratio") %>% 
  select_if(is.numeric) %>%
  cor()

soccer_distribution <- 1-abs(soccer_correlation)

soccer_distribution %>% 
  as.dist() %>% 
  hclust(method = "complete") %>% 
  as.dendrogram() %>% 
  ggdendro::ggdendrogram(rotate = T) +
  labs(title = "Clustering with Complete Linkage",
       subtitle = "Women's NCAA Soccer 2018 & 2019")
```

```{r}
n_distinct(womens_ncaa_soccer_20182019$team)
summary(womens_ncaa_soccer_20182019$team_games)
hist(womens_ncaa_soccer_20182019$team_games)

dim(womens_ncaa_soccer_20182019)
```