# Partially Constrained Group Variable Selection to Adjust for Complementary
Unit Performance in American College Football

A. Skripnikova CONTACT A. Skripnikov. Email<EMAIL_ADDRESS>aNew College
of Florida, Sarasota, FL 34243, USA

###### Abstract

Given the importance of accurate team rankings in American college football
(CFB) - due to heavy title and playoff implications - strides have been made
to improve evaluation metrics across statistical categories, going from basic
averages (e.g. points scored per game) to metrics that adjust for a team’s
strength of schedule, but one aspect that hasn’t been emphasized is the
complementary nature of American football. Despite the same team’s offensive
and defensive units typically consisting of separate player sets, some aspects
of your team’s defensive (offensive) performance may affect the complementary
side: turnovers forced by your defense could lead to easier scoring chances
for your offense, while your offense’s ability to control the clock may help
your defense. For 2009-2019 CFB seasons111Data and source code files are made
available at
https://github.com/UsDAnDreS/SUBMISSIONS_OffenseST_DefenseST_GLM_LASSO_adjusted_rankings,
we incorporate natural splines with group penalty approaches to identify the
most consistently influential features of complementary football in a data-
driven way, conducting partially constrained optimization in order to
additionally guarantee the full adjustment for strength of schedule and
homefield factor. We touch on the issues arising due to reverse-causal nature
of certain within-game dynamics, discussing several potential remedies.
Lastly, game outcome prediction performances are compared across several
ranking adjustment approaches for method validation purposes.

###### keywords:

group penalty, LASSO, natural splines, regularized estimation, reverse
causality, sports statistics

##  1 Introduction

American football, while requiring immense physical ability, is an example of
an extremely strategic sport. Requiring full mental engagement on every play
from all 11 players your team puts on the field, it underlines the importance
of tactical preparation and accurate evaluation of the opposing team in
advance of the match-up. Besides that, in American college football
specifically, the ability to objectively evaluate team performance is pivotal
for the purposes of rankings, with the latter fully dictating the college
teams that get into the championship playoffs and high-profile bowl games,
including all the financial benefits to come with it. Historically, classical
averages have been used to evaluate team’s performance in a certain
statistical category [2]. For example, American football team’s offensive
output from yardage standpoint is traditionally measured in total yards the
team gains per game; its tendency to turn the ball over to the opposing team
is typically captured in their average turnover count per contest; team’s pure
scoring ability is ubiquitously described by an average amount of points
scored per game. While useful, such basic averages are extremely flawed
whenever trying to objectively evaluate and compare teams within the setup of
Football Bowl Subdivision (FBS) of American college football.

Football Bowl Subdivision (FBS), as of year 2019, consists of $131$
participating college football programs, split into $11$ conferences. Each
team is typically scheduled to play only $12$ opponents during the regular
season, with most of the games taking place against the teams from the same
conference, hence leaving the ranking picture incomplete in terms of relative
strengths for teams across the entire FBS. Certain conferences could end up
being especially strong, or weak, during a particular season, leading to
uneven quality of opposition which can’t be reflected in the calculating each
team’s basic per-game averages alone. That gave rise to what’s now known as
”strength of schedule” adjustment, originating in [15] via an offense-defense
model where points scored by a team are assumed to be a function that team’s
offensive strength and opposing team’s defensive strength, therefore adjusting
for quality of the opponent. Lastly, given the importance of home-field
advantage [10, 25], in addition to adjusting for strength of scheduled
opponents [15] also accounted for whether a team had a home- or road-heavy
schedule.

In this work, besides the strength of schedule and homefield adjustments, we
wanted to also incorporate the complementary nature of American football.
Unlike most other sports, in American football offense and defense are played
by separate, in most cases non-overlapping, units of players, that can’t be on
the field at the same time. It might lead to an assumption that performance of
a team’s offense could be treated independently of that same team’s defense.
Nonetheless, historically it’s been shown that the two sides are likely to
complement each other, e.g. turnovers forced by your defense could lead to
easier scoring chances for your offense, while your offense’s ability to
control the clock may help your defense. That concept has been termed
”complementary football”, and focus of this work was on determining the most
critical features of collegiate American football in affecting the
complementary side’s performance. Specifically, we utilized a variety of
regularized estimation techniques that impose a penalty for the purposes of
variable selection, implementing partially constrained optimization with the
goal of guaranteeing strength of schedule and homefield adjustment. Somewhat
related ideas were studied in [18, 4] but in application to modelling the
dependence between goals scored and allowed by a team in a soccer match. To
the best of our knowledge, there hasn’t been any involved research done on
complementary impacts an American football offense can have on its own
defense, and vice versa, after having adjusted for strength of schedule and
homefield factor.

The remainder of the paper is structured as follows. Section 2 lays out the
details of data collection and cleaning, formulates all regularized estimation
approaches for the purpose of finding the most critical complementary football
features in affecting certain statistical categories. The results of feature
selection, ensuing ranking adjustments, along with their impact on predicting
game outcomes, are presented in Section 3. Moreover, the issue of endogeneity
is brought up, resulting from reverse causal within-game dynamics of a
football game, and several partial solutions discussed, main one involving the
use of solely efficiency-based statistics. Section 4 contains concluding
remarks and discussions of future work.

##  2 Materials and methods

###  2.1 Data collection and cleaning

Game-by-game data for 2009-2019 FBS seasons was obtained via web scraping from
two primary sources: www.sports-reference.com/cfb/ and www.cfbstats.com/. The
former website provides detailed scoring and defensive play, and has webpage
format conducive to seamless scraping of analysis-ready data tables into R
Statistical Software [22] \- primary tool used for statistical analysis in
this work - via $rvest$ package [27]. The latter website helped augment the
data with more statistical measurements, although it required a more profound
web scraping effort while using $RSelenium$ package [14]. Certain defensive
play statistics (e.g. tackles, sacks, forced fumbles) were lacking for ”non-
major” teams, meaning colleges not in the FBS. Despite rare occurrence of such
observations (only about $7\%$), we wanted to incorporate all relevant games
into our analysis and, instead of removing them, we performed data imputation
[11] via linear regression. We regressed the missing statistics - e.g. sacks -
on the other statistical categories that were available for non-major teams
(e.g. yards, touchdowns, punts, etc), utilizing observations for all FBS teams
in fitting this regression, to later make predictions of that missing
statistic for non-major teams.

Having obtained game-by-game data on teams’ performances in a multitude of
statistical categories, we converted it into a format where each observation
describes the cumulative performance throughout the game whenever one team’s
offensive side (offense or offensive special teams) - and, consequently,
opponent’s defensive side (defense or defensive special teams) - was on the
field. E.g. for this team’s offense, it would include yards gained, touchdowns
scored, passes completed, field goals scored, etc, while for opponent’s
defense it would show total tackles, tackles for loss, sacks, quarterback
hurries, etc. To clarify, ”offensive special teams” imply the field goal
kicking and kick/punt returning player sets, while ”defensive special teams”
represent players participating in field goal blocking and kick/punt coverage.
That data format results into two observations per each game, one representing
plays produced when first team’s offensive side (and second team’s defensive
side) was on the field, second - when second team’s offensive side (and first
team’s defensive side) was on the field. This implies potential within-game
dependence introduced by such pairs of observations, but after having
calculated Intraclass Correlation Coefficients [3] across the models
considered in this paper, virtually none of those correlations ended up being
over 0.10 in magnitude, pointing to independence as a reasonable assumption
for our data.

###  2.2 Variable selection methodology outline

Below we introduce the modeling notation and main variable selection
approaches to help adjust team’s offensive (or defensive) quality for the
strength of schedule, homefield factor, and complementary nature of American
football.

####  2.2.1 Notation

Assume we have a total of $n$ teams in the league, and we are interested in a
particular statistical category $y$, be it yards (per game), turnovers,
points, etc. Let $y_{ij}$ denote the performance of team $i$ against team $j$
in category $y$, $n_{i}$ \- total number of opponents team $i$ has had through
the season, $i,j=1,\dots,n$. For example, if $y=$ {total yards gained per
game}, then $y_{ij}$ would correspond to the amount of yards team $i$ gained
when playing against team $j$, and, by symmetry, the amount of yards team $j$
allowed in the game against team $i$. In case teams $i$ and $j$ had to face
each other $L_{ij}$ times during the season, we add another index $l$. Then,
$y_{ijl}$ corresponds to performance of team $i$ in category $x$ during its
$l^{th}$ meeting with team $j$, $i,j=1,\dots,n,\ l=1,\dots,L_{ij}$.

Next, let $h_{ijl}$ denote a homefield indicator for $l^{th}$ game between
teams $i$ and $j$, taking on value $1$ if team $i$ is at home, $0$ if the game
site is neutral, $-1$ if team $i$ is at home. Such numerical encoding was
intuitive (typically easier to play at home than at a neutral site, and at a
neutral site than on the road) and also got confirmed by running a dummy-
variable encoding scheme, having shown increases (in points, touchdowns
scored, yards gained) for home and decreases for road games compared to the
neutral site baseline.

To incorporate adjustment for strength of the schedule, we introduce concepts
of offensive (defensive) worth of the ”league-average opponent”, and offensive
(defensive) margin for a team. In regards to a particular statistical category
$y$, one can define the league-average opponent via two parameters - offensive
and defensive worth. E.g. for points per game, offensive (defensive) worth of
the average opponent is the average points per game scored (allowed) by all
teams across all games that could have been played against one another
throughout the course of the season. Due to symmetry (team $i$ scoring
$y_{ij}$ points against team $j$ is equivalent to team $j$ allowing $y_{ij}$
points to team $i$), both offensive and defensive worth represent the same
value for the league-average opponent, which we denote as $\mu$. Now, for each
team $i$ we can posit parameters capturing two aspects of its performance
within a statistical category - offensive margin $\alpha_{i}$ and defensive
margin $\beta_{i}$. Offensive (defensive) margin describes by how much a team
would outperform the aforementioned defensive (offensive) worth $\mu$ of the
average opponent. The main assumption when adjusting for strength of schedule
is that performance of team $i$ against team $j$ in category $y$ is
attributable to both the offensive margin $\alpha_{i}$ of team $i$ and
defensive margin $\beta_{j}$ of team $j$.

Lastly, presuming that we consider $C$ complementary football statistics,
let’s use $x_{c,jil},\ c=1,\dots,C,$ to denote the value of $c^{th}$ statistic
that’s complementary to the $y_{ijl}$, meaning that $x_{c,jil}$ is obtained
when the defense (complementary unit for the offense) of $i^{th}$ team and
offense (complementary unit for the defense) of $j^{th}$ team were on the
field during their $l^{th}$ game of the season between these two teams.

####  2.2.2 Natural cubic splines

To model potentially non-linear effects of complementary football features, as
a well-known method we utilized natural cubic splines [16], where one uses a
mixture of piece-wise cubic and linear polynomials, smoothly connected at a
set of $K$ knots placed across the range of the explanatory variable. It
results into each complementary statistic $x_{c}$ being represented by a set
of basis functions $N_{1}(x_{c}),N_{2}(x_{c}),\dots,N_{K-1}(x_{c})$. For more
detail, see [16], keeping in mind that the intercept basis function
$N(x_{c})=1$ for each individual complementary statistic $x_{c}$ in our case
is omitted from the basis due to being folded into the overall model’s
intercept. We chose to use $K=5$ knots placed at $0.00$-, $0.25$-, $0.50$-,
$0.75$\- and $1.00$-quantiles, providing just enough flexibility to capture
any clear non-linearity, while decreasing chances of overfitting and low
interpretability that come with overly flexible fits. That results into each
complementary statistic $x_{c}$ being represented by four basis functions,
with its partial effect on response calculated via a linear combination of
these functions.

####  2.2.3 Partially constrained group penalization

Natural splines approach leads to each complementary football feature being
represented by a group of several parameters, where inclusion of a feature
into the model implies including the entire said parameter group, while
exclusion indicates setting that whole group of parameters to zero. Such
setting warrants use of group penalization approaches [28, 17, 26]. As we also
want to guarantee adjustment for homefield effect and strength of schedule,
certain parameters will be left unpenalized, resulting into the following
general format of partially constrained group penalization criteria:

$\min_{\begin{subarray}{c}\mu,\delta,\\{\alpha_{i}\\}_{1}^{n},\\{\beta_{j}\\}_{1}^{n}\\\
\\{\boldsymbol{\gamma}_{c}\\}_{1}^{K-1},\ c=1,\dots,C\end{subarray}}\
\sum_{i,j=1}^{n}\sum_{l=1}^{L_{ij}}(y_{ijl}-(\mu+\delta
h_{ijl}+\alpha_{i}+\beta_{j}+\sum_{k=1}^{K-1}\gamma_{1,k}N_{k}(x_{1,jil})+\dots+\sum_{k=1}^{K-1}\gamma_{C,k}N_{k}(x_{C,jil})))^{2}+\\\
\lambda\sum_{c=1}^{C}P_{c}(||\boldsymbol{\gamma}_{c}||_{2}),$ (1)

where $\boldsymbol{\gamma}_{c}=(\gamma_{c,1},\dots,\gamma_{c,K-1})$, capturing
all the natural cubic spline coefficients for $c^{th}$ complementary
statistic;
$||\boldsymbol{\gamma}_{c}||_{2}=\sqrt{\sum_{k=1}^{K-1}\gamma_{c,k}^{2}}$ \-
the Eucledian ($L_{2}$) norm; $P_{c}(\cdot)$ indicates a specific group
penalty function to be utilized for the coefficients of $c^{th}$ complementary
statistic; $\lambda$ \- tuning parameter responsible for strength of penalty
(larger $\lambda$ implies larger penalty for high values of
$||\boldsymbol{\gamma}_{c}||_{2}$). Here, the penalty is imposed solely on
coefficients $\boldsymbol{\gamma}_{1},\dots,\boldsymbol{\gamma}_{C}$ for
complementary football statistics, with the goal of inferring the most
critical ones for statistical adjustment purposes. This partially constrained
optimization criteria guarantees inclusion of homefield effect ($\delta$) and
strength of schedule adjustment effects ($\alpha_{i},\beta_{j},\
i,j=1,\dots,n$).

We considered an array of various group penalties, to subsequently identify
the most consistently selected features regardless of a particular penalty
choice. First, we utilized a classic group least absolute shrinkage and
selection operator (group LASSO, [28]),
$P_{c}(||\boldsymbol{\gamma}_{c}||_{2})=||\boldsymbol{\gamma}_{c}||_{2}$,
which is a natural extension of regular LASSO [24] and selects variables in a
grouped manner. Due to applying the same strength of penalty ($\lambda$) to
all coefficients, regular LASSO was shown to struggle with efficiency [12] and
model selection [20]. Given the shared nature of the penalty function, similar
issues persist for group LASSO. In case of regular LASSO, adaptive LASSO [30]
method was proposed to address this issue via introducing coefficient-specific
weights, subsequently yielding consistent variable selection and estimators
satisfying oracle property. The exact same intuition was carried out in [26]
to implement adaptive group LASSO. In our case, it leads to penalty function
$P_{c}(||\boldsymbol{\gamma}_{c}||_{2})=w_{c}||\boldsymbol{\gamma}_{c}||_{2},\
w_{c}=1/||\hat{\boldsymbol{\gamma}}^{LS}_{c}||^{\gamma}_{2}$, where
$\hat{\boldsymbol{\gamma}}^{LS}_{c}$ denotes the least squares estimate for
the effect of $c^{th}$ complementary statistic resulting from optimizing
criteria (1), but with the penalty term
($\lambda\sum_{c=1}^{C}P_{c}(||\boldsymbol{\gamma}_{c}||_{2})$) excluded.
While adaptive LASSO poses a two-stage approach, with estimation of weights
$\\{w_{c}\\}$ preceding the actual penalization, non-convex penalties present
a single-stage alternative that achieves analogous bias reduction for large
regression coefficients. We used two of the most popular non-convex penalties
- smoothly clipped absolute deviations (SCAD) [12] and minimax concave penalty
(MCP) [29] \- which were originally defined for regular LASSO, but can be
naturally extended to group LASSO setting. For details on the exact functional
form $P_{c}(\cdot)$ for those penalties see [17]. Three aforementioned
penalties - adaptive, SCAD, MCP - have an extra tuning parameter (generally
denoted as ”$\gamma$” in the original papers), that in each case we picked in
a data-driven way according to Bayesian Information Criterion (BIC) [21].

To implement all of the aforementioned group LASSO penalties for optimization
task (1) we utilized $grpreg$ package [7], along with 20-fold cross-validation
(CV) [8] and BIC [21] to select several candidate values for the tuning
parameter $\lambda$. For CV, we opted to use a $\lambda$ value that achieved
CV test error within one standard deviation of the minimal CV test error,
while providing a sparser model (selecting less variables). Moreover, to
reduce the randomness effect of the 20-fold CV procedure, we ran it five times
for each season, only counting features that got selected at least four out of
five times. We used BIC and sparse CV $\lambda$-value because, amongst all of
the most ubiquitous criteria, these ones imposed strongest penalties on overly
complicated models, leading to a small set of the most important features.
Most importantly, if a complementary football feature gets selected by these
criteria, it is virtually guaranteed to also be selected by any other
ubiquitous criteria.

Each complementary statistic was evaluated on consistency of its selection
across 2009-2019 FBS seasons by each of four main variable selection
methodologies outlined in Section 2 (classic and adaptive group LASSO, group
SCAD and MCP penalties). Specifically, each method was given equal weight of
$\frac{1}{4}$ in the final consistency metric, while within any individual
method an equal weight of $\frac{1}{2}$ was granted to each of the two
considered tuning parameter selection approaches (BIC and sparse CV
$\lambda$-value).

Lastly, to calculate the final adjustment for each season, natural cubic
splines were fitted solely with the complementary football features that were
consistently picked across variable selection methods, while centering
variables representing complementary statistics around their respective means.
That way, team $i$’s offensive worth $\mu+\alpha_{i}$ and defensive worth
$\mu+\beta_{i}$ will be projected not only onto league-average opponent and
neutral site, but also onto league-average complementary side.

####  2.2.4 Modeling assumptions given response variable type

We considered six statistics to be adjusted, breaking it down by the type of a
statistic (points, touchdowns, yards) and whether ”counter-plays” are being
accounted for or not. ”Counter-plays” are plays that happen directly after a
turnover, before (if at all) the complementary side takes the field. For
example, we could solely count touchdowns that your offense scored, or we
could also take the margin between touchdowns they scored and touchdowns they
allowed the opposing defense to score (via a returned interception or fumble).
Same for points and yards. The primary reason for creating these margin-
statistics is to eventually use them as features to forecast game outcomes,
given that they account for totality of the impact whenever a particular side
took the field, rather than reflecting solely their positive impact. Lastly,
due to the discrete, low-count, nature of the touchdowns category (while
points and yards are well-approximated by continuous distributions), we also
considered tailored count response approaches (e.g. Poisson and negative
binomial generalized linear models), but classic linear regression yielded a
fit of either superior or similar quality to that of the alternatives, hence
was picked moving forward.

##  3 Results

This section is broken up into two main parts based on the complementary
football feature sets being considered. First part works with the entire set
of statistical categories originally obtained from www.sports-
reference.com/cfb/ and www.cfbstats.com/. Second part discusses using an
efficiency-based feature set of solely per-play and per-possession statistics
as a way to partially address the issue of reverse causality permeating the
first set. Both feature sets can be found in the Appendix.

###  3.1 Using entire feature set of complementary football statistics.

####  3.1.1 Collinearity.

Prioritizing interpretation and consistency of variable selection results, we
disposed of both perfect and near-perfect multicollinearity amongst our
features. For perfect multi-collinearity cases of one statistical category
representing a sum of sub-categories (e.g. passing attempts equals completions
plus incompletions) we chose to go with the sub-categories (e.g. retain
completions and incompletions, while dropping passing attempts). The only
exceptions were the cases where breakdown into sub-categories was not useful
in terms of its effects on the complementary side. For example, as long as it
is a turnover, distinguishing between fumbles and interceptions is
unnecessary, as the impact on the complementary side won’t change based solely
on that aspect. On the other hand, there is a considerable difference between
scoring turnovers (when defense instantly scores) and non-scoring ones
(complementary side gets the ball back in the latter case, doesn’t in the
former), or between pass completions and incompletions (the latter stops the
game clock, the former - doesn’t). For near-perfect multi-collinearity cases,
predominantly represented by strong pairwise correlations (over 0.90 in
magnitude), we picked the more ubiquitously used statistic to represent the
highly correlated group. For example, points per game being are strongly
correlated with touchdowns per game (consistently over 0.95 each season), and
we pick points per game as a more all-encompassing and recognizable statistic.
For the feature set resulting from this step, which ended up containing
roughly $50$ variables, see the Appendix.

####  3.1.2 Selected complementary statistics, issue of reverse causality.

Among statistical categories that were consistently selected as affecting
team’s complementary side in regards to points scored/allowed in a game, we
got rushing attempts, pass incompletions, non-scoring turnovers, number of
attempted special teams returns of the ball, and yards gained on those
returns. Some of these are perfectly reasonable, e.g. non-scoring turnovers
provide the complementary side with the ball and a potentially easy scoring
opportunity. For a statistic such as pass incompletions, while expectedly
having a positive impact on the points scored on complementary side (stopping
the clock, hence leaving more game time to complementary side), it rarely
features in football discussions of the most pivotal factors driving
complementary scoring. Lastly, several statistical categories simply don’t
make as much sense, e.g. number of attempted returns or rushes, regardless of
yards gained on such plays.

Why do we observe such unintuitive strong effects for some of the statistics?
Reverse causality is the most likely answer. For example, while pass
incompletions facilitate complementary side scoring, they can also be a
partial function of play selection by the team that’s behind in the score.
Such team would be forced to throw the ball more often in order to gain yards
quicker and catch up, while also stopping the clock in case of an unsuccessful
play (unlike an ineffective rushing attempt, an incompleted pass actually
stops the clock), trying to score while taking as little of the game clock as
possible. For rushing attempts, on the other hand, reverse logic applies in
terms of play selection: the leading team is likely to call more rushing plays
to keep the game clock running, while the trailing team is less likely to do
so, creating a reverse-causal component in estimating the effect of that
complementary football feature. Lastly, the special teams return play data
(attempts, yards) suffers from the same issue, but has an even clearer time-
sequential component to it: a return attempt tends to follow directly after
someone having just scored on the complementary side. Despite lack of access
to sequential play-by-play data, game-by-game data allowed to calculate that
68% of attempted returns followed after a kickoff, which in its turn is most
likely to come after a score on the complementary side (60% chance it was a
7-point touchdown, 18% chance - a 3-point field goal). This clearly shows a
reverse-causal pattern of a special teams return following a complementary
side score, not the other way around.

How can we potentially alleviate these reverse-causal effects? There are
several potential solutions, the most obvious one being to obtain play-by-play
data across all games. With careful consideration and proper domain
understanding, knowledge of play sequencing could provide us with further
insight into the causal flow of within-game dynamics. At the very least, it
could make clear that return attempts by one’s offense typically directly
follow after the score by the opposing offense, avoiding the aforementioned
reverse-causal trap. Unfortunately, in many cases one lacks the capability to
access and process such fine-grained play-by-play data, and therefore is left
looking for more feasible alternatives.

Second approach could be adjusting for endogeneity, which broadly refers to
cases of an explanatory variable being correlated with an error term, hence
violating one of the least squares modeling assumptions and subsequently
biasing the estimates. Reverse causality is a special case of endogeneity,
with explanatory variables suffering from it being called ”endogenous”, due to
being determined by factors ”within the system” (e.g. within the game
situation), as opposed to ”exogenous”, which are fully defined by external
factors (e.g. weather). Issue of endogeneity has been discussed across several
application domains [19, 9, 23], with instrumental variables [6] proposed as
one of the most ubiquitous remedial approaches. A strong instrumental is an
exogenous variable that strongly correlates with the endogenous explanatory
variable, while being uncorrelated with the error term. Such variables have
been notoriously difficult to find, resulting into cases where weak
instruments only exacerbated estimation bias [5]. Moreover, in our specific
case, virtually every complementary football feature can be considered
endogenous due to being a partial byproduct of within-game dynamics, making it
that much harder to find good exogenous instrumental variables, which left us
looking for other alternatives.

Lastly, with unavailability of play-by-play data and complexity of endogeneity
adjustment logistics, we decided to pursue an intuition-driven approach of
focusing solely on efficiency statistics. While such cumulative metrics as
pass incompletions or rushing attempts can be a strong function of strategic
play selection within the context of the game, efficiency-based statistics are
more reflective of team’s effectiveness once a certain play was called, rather
than heavily depending on the overall game context. For example, if a team
pursues a pass-heavy strategy in an effort to catch up with the opponent,
their cumulative passing statistics - attempts, completions, incompletions,
yards - will inevitably increase, while their efficiency-based alternatives -
completion percentage, yards per attempt - are not guaranteed to change one
way or the other. Hence, to reduce the chance of reverse causality, we opt for
using solely efficiency-based statistical categories, showing how successful a
team is once a certain play was called, as opposed to getting overly affected
by contextual strategic play choices.

###  3.2 Using only efficiency-based complementary football statistics

####  3.2.1 Efficiency-based features

Having disposed of collinearity issues and unnecessary statistical category
breakdowns (see Section 3.1.1 for analogous reasoning), in the row names of
Figure 1 one can witness the efficiency-based complementary football
statistics that ended up being utilized for variable selection. For each
statistical category we had to decide between using a per-play or per-
possession efficiency. We selected to look at plays that terminate a
possession (turnover, points scored, punt, etc) on per-possession basis,
making them into possession percentages. For example,
$Off\\_ST.NonScoring.TO.PossPct$ corresponds to the percentage of possessions
ending in the team’s offensive side turning the ball over to their defensive
side (as opposed to turning it over for an instant score by the opposing
defensive side, which is captured via $Def\\_ST.TD.PossPct$). For offensive
scoring, we decided against breaking it down scoring possessions into
touchdowns and field goals, in big part due to both outcomes resulting into
same direct impact on the complementary side (a kickoff). Nonetheless,
percentage of scoring (and non-scoring) possessions was heavily correlated
with a more recognizable and informative metric of points per possession,
hence was dropped in favor of the latter. Lastly, punts and safeties were
combined into a single category because, despite safeties also resulting into
two points scored (unlike zero points in case of punts), the biggest advantage
of either play is the possession turning over to the other side, along with
the similarity of the ball getting punted across the field in both cases.

On the other hand, virtually all other metrics were treated on per-play basis,
including yards or first downs gained, overall tackles, tackles for loss,
forced fumbles. For yards gained, there were several alternative
considerations to simply looking at average yards per play by the offense:
breaking it down into yards per pass and rush attempt, or also including
special teams returns as plays into the overall per-play yardage. The former
idea was tested via cross-validation, having shown that total yards per play
outperformed the pass/rush yardage breakdown approach. The latter idea
resulted into special teams return yardage overly affecting the per-play
efficiency values, leading to this statistic getting selected every single
time, mostly due to the aforementioned reverse causality issue pertaining to
return play data. Meanwhile, some other statistics would result into
percentages due to their binary nature, e.g. pass completion, third down
conversion, blocked kicks/punts. One may notice that we didn’t include
statistics such as field goal percentage, average yards per punt or per
attempted return, red zone conversion percentage, etc. That is due to the
issue of a low (sometimes zero) denominator, where efficiency is calculated
over a really small sample size within a game (e.g. only over two field goal
attempts), and hence is not necessarily representative of team’s true ability
in that category, or the statistic’s impact on the game. Moreover, in
relatively frequent cases when sample size is zero, the value has to be set to
zero or left undefined, neither option being optimal.

Lastly, such statistical categories as ”sacks to pass attempts ratio” were
defined in an effort to adjust for effects of pass-heavy play calling being
conducive to such defensive plays as sacks or quarterback hurries. While a
sack is not recorded as a passing attempt, and ideally one needs the number of
quarterback dropbacks to create a proper efficiency metric, due to data
limitations this was a reasonable alternative. In similar vein, percentage of
forced fumbles by defense was only counted for rushing attempts and pass
completions - the only plays where a fumble was possible.

Figure 1: A weighted proportion of times a natural cubic spline representing a
particular complementary football statistic (rows) was selected across
different group penalties and tuning parameter selection approaches. Figure 2:
Non-linear effects of complementary statistics on respective response
categories, controlling for strength of schedule and homefield factor. Left:
Effect of non-scoring turnovers forced by your defense on points scored by
your offense. Right: Effect of per-play yardage allowed by your defense on
yards gained by your offense. These particular effect graphics were obtained
for a Navy versus Army match-up at a neutral site. Changing the match-up would
only affect the intercepts, hence there’s no loss of generality in depicting
the nature of the effect.

####  3.2.2 Selected complementary statistics.

Figure 1 demonstrates the efficiency-based features of complementary football
(rows) that affected each respective statistical category to be adjusted
(columns). As described in Section 2.2.4, ”margin” accounts for counter-plays
in a respective category.

One can witness nonscoring turnover percentage forced by your defensive side
(or committed by opponent’s offensive side, hence the ”$Off\\_ST$” prefix) as
the sole statistic that is consistently picked ($\geq 74\%$ of the time) in
terms of impacting points and touchdowns scored by your (complementary)
offensive side. Meanwhile, for both total yardage response categories, the
only complementary football statistic picked with strong stability is the
total yards per play ($\geq 94\%$ of the time). To evaluate the nature of
these consistently selected complementary statistics, we fit the respective
natural cubic spline regressions, with Figure 2 depicting the effects of your
defense’s forced non-scoring turnovers on your offensive scoring (left), and
of per-play yardage allowed by your defense (or gained by opponent’s offense,
hence the ”$Off$” prefix) on total yardage gained by your offense (right).
Each non-scoring turnover forced by your defense mostly leads to higher
offensive scoring, with majority of years exhibiting an approximately linear
relationship. This is a reasonable finding due to each such turnover giving
your offense the ball with a high potential for a good scoring opportunity. On
the other hand, yards allowed per play by your defense exhibit a clear non-
linear impact on your offense’s ability to gain yards, which actually aligns
well with the intuition. It is reasonable to expect for defenses allowing
either extremely few ($<2.5$) or very many ($>5$) yards per play to take less
of the game clock themselves, albeit for polar opposite reasons. Good defenses
(low per-play yardage allowed) would quickly stop the opposing offense from
scoring, and get the ball over to their complementary offense (typically via a
punt), putting the latter in a position to gain yards of their own. Bad
defenses (high per-play yardage allowed) would let the opposing offense score
points quickly, which still leads to their complementary offense getting the
ball. Meanwhile, if your defense allows moderate per-play yardage ($2.5$ to
$5$), that typically makes for longer offensive possessions by the opponent,
taking away from your own offensive unit the game time that could be used to
gain more yards.

Note that, given the two-sided nature of game sports, the reasoning above
could also be observed from the opposing team’s viewpoint of their offensive
side affecting their complementary defensive side. Each gain on offense could
be considered as a loss on defense (e.g. a touchdown scored on offense is a
touchdown allowed on defense) and vice versa (e.g. a sack gained by the
defense is a sack allowed by the offense).

####  3.2.3 Ranking shifts in points scored per game.

To illustrate the mechanism behind adjustments for complementary football
features, we showcase its impact on a team’s offensive scoring ranking if one
accounts for turnovers forced by that team’s defense. Table 1 demonstrates the
ranking shifts in points per game scored by offense in 2012 season when, on
top of adjusting for strength of schedule and homefield factor, defense’s
ability to turn the ball over to the offense is accounted for. Oregon’s
scoring, while still remaining the best in FBS, dipped strongly (by $1$ point
per game) once projected onto an average complementary unit, vastly due to
their defense being Top-4 in forcing the non-scoring turnovers (same happened
to the Washington offense, complemented by a Top-5 defense). On the other
hand, Texas A&M’s offense climbed into Top-3 as it didn’t benefit from as
opportunistic of a defense (ranked outside Top-100). The biggest increases in
projected points per game happened to Texas Tech and South Florida offenses,
both working aside some of the worst defenses in terms of forcing turnovers
(ranked 124th and 123rd, respectively, out of 125 teams). The intuition behind
ranking shifts in other years was identical.

Team | Points Scored | Non-Scoring Turnovers Forced
---|---|---
| (Per Game) | (Per Possession)
| Oregon
---
Louisiana Tech
Texas A&M
Oklahoma State
Baylor
Clemson
Oklahoma
Alabama
West Virginia
Tennessee
Texas Tech
Georgia
….
South Florida
….
Washington
….
| Value | (Shift) | Rank | (Shift)
---|---|---|---
48.43 | (${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\Downarrow}$ 1.02) | 1 | (${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}$)
47.12 | (${\color[rgb]{0.0,0.5,0.0}\definecolor[named]{pgfstrokecolor}{rgb}{0.0,0.5,0.0}\Uparrow}$ 0.19) | 2 | (${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}$)
46.67 | (${\color[rgb]{0.0,0.5,0.0}\definecolor[named]{pgfstrokecolor}{rgb}{0.0,0.5,0.0}\Uparrow}$ 0.75) | 3 | (${\color[rgb]{0.0,0.5,0.0}\definecolor[named]{pgfstrokecolor}{rgb}{0.0,0.5,0.0}\Uparrow}$ 1)
46.15 | (${\color[rgb]{0.0,0.5,0.0}\definecolor[named]{pgfstrokecolor}{rgb}{0.0,0.5,0.0}\Uparrow}$ 0.11) | 4 | (${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\Downarrow}$ 1)
45.25 | (${\color[rgb]{0.0,0.5,0.0}\definecolor[named]{pgfstrokecolor}{rgb}{0.0,0.5,0.0}\Uparrow}$ 0.01) | 5 | (${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}$)
41.02 | (${\color[rgb]{0.0,0.5,0.0}\definecolor[named]{pgfstrokecolor}{rgb}{0.0,0.5,0.0}\Uparrow}$ 0.08) | 6 | (${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}$)
40.93 | (${\color[rgb]{0.0,0.5,0.0}\definecolor[named]{pgfstrokecolor}{rgb}{0.0,0.5,0.0}\Uparrow}$ 0.26) | 7 | (${\color[rgb]{0.0,0.5,0.0}\definecolor[named]{pgfstrokecolor}{rgb}{0.0,0.5,0.0}\Uparrow}$ 1)
40.57 | (${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\Downarrow}$ 0.17) | 8 | (${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\Downarrow}$ 1)
39.74 | (${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\Downarrow}$ 0.06) | 9 | (${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}$)
39.69 | (${\color[rgb]{0.0,0.5,0.0}\definecolor[named]{pgfstrokecolor}{rgb}{0.0,0.5,0.0}\Uparrow}$ 0.58) | 10 | (${\color[rgb]{0.0,0.5,0.0}\definecolor[named]{pgfstrokecolor}{rgb}{0.0,0.5,0.0}\Uparrow}$ 2)
39.44 | (${\color[rgb]{0.0,0.5,0.0}\definecolor[named]{pgfstrokecolor}{rgb}{0.0,0.5,0.0}\Uparrow}$ 1.22) | 11 | (${\color[rgb]{0.0,0.5,0.0}\definecolor[named]{pgfstrokecolor}{rgb}{0.0,0.5,0.0}\Uparrow}$ 4)
38.98 | (${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\Downarrow}$ 0.19) | 12 | (${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\Downarrow}$ 1)
… | … |  |
21.99 | (${\color[rgb]{0.0,0.5,0.0}\definecolor[named]{pgfstrokecolor}{rgb}{0.0,0.5,0.0}\Uparrow}$ 1.16) | 87 | (${\color[rgb]{0.0,0.5,0.0}\definecolor[named]{pgfstrokecolor}{rgb}{0.0,0.5,0.0}\Uparrow}$ 6)
… | … |  |
23.72 | (${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\Downarrow}$ 1.02) | 80 | (${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\Downarrow}$ 5)
… | … |  |
| Value | Rank
---|---
0.19 | 4
0.13 | 52
0.08 | 114
0.12 | 75
0.12 | 61
0.12 | 67
0.09 | 111
0.17 | 13
0.10 | 97
0.09 | 108
0.06 | 124
0.16 | 23
… | …
0.06 | 123
… | …
0.19 | 5
… | …
Table 1: Final 2012 season rankings (and ranking shifts) in points scored per game by team’s offense, adjusted for defense’s ability to turn the ball over to team’s offense. Team | Points Allowed | Non-Scoring Turnovers Committed
---|---|---
| (Per Game) | (Per Possession)
| Alabama
---
Notre Dame
Florida
Michigan State
Stanford
Brigham Young
South Carolina
Florida State
Rutgers
Kansas State
Louisiana State
Texas Christian
….
Pittsburgh
….
Illinois
….
| Value | (Shift) | Rank | (Shift)
---|---|---|---
8.70 | (${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\Uparrow}$ 0.42) | 1 | (${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}$)
11.85 | (${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\Uparrow}$ 0.42) | 2 | (${\color[rgb]{0.0,0.5,0.0}\definecolor[named]{pgfstrokecolor}{rgb}{0.0,0.5,0.0}\Uparrow}$ 1)
11.88 | (${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\Uparrow}$ 0.80) | 3 | (${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\Downarrow}$ 1)
14.03 | (${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\Uparrow}$ 0.46) | 4 | (${\color[rgb]{0.0,0.5,0.0}\definecolor[named]{pgfstrokecolor}{rgb}{0.0,0.5,0.0}\Uparrow}$ 1)
14.05 | (${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\Uparrow}$ 0.88) | 5 | (${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\Downarrow}$ 1)
14.26 | (${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\Uparrow}$ 0.10) | 6 | (${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0}$)
15.74 | (${\color[rgb]{0.0,0.5,0.0}\definecolor[named]{pgfstrokecolor}{rgb}{0.0,0.5,0.0}\Downarrow}$ 0.60) | 7 | (${\color[rgb]{0.0,0.5,0.0}\definecolor[named]{pgfstrokecolor}{rgb}{0.0,0.5,0.0}\Uparrow}$ 1)
15.97 | (${\color[rgb]{0.0,0.5,0.0}\definecolor[named]{pgfstrokecolor}{rgb}{0.0,0.5,0.0}\Downarrow}$ 0.44) | 8 | (${\color[rgb]{0.0,0.5,0.0}\definecolor[named]{pgfstrokecolor}{rgb}{0.0,0.5,0.0}\Uparrow}$ 1)
16.56 | (${\color[rgb]{0.0,0.5,0.0}\definecolor[named]{pgfstrokecolor}{rgb}{0.0,0.5,0.0}\Downarrow}$ 0.31) | 9 | (${\color[rgb]{0.0,0.5,0.0}\definecolor[named]{pgfstrokecolor}{rgb}{0.0,0.5,0.0}\Uparrow}$ 2)
16.84 | (${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\Uparrow}$ 0.88) | 10 | (${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\Downarrow}$ 3)
17.06 | (${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\Uparrow}$ 0.25) | 11 | (${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\Downarrow}$ 1)
17.31 | (${\color[rgb]{0.0,0.5,0.0}\definecolor[named]{pgfstrokecolor}{rgb}{0.0,0.5,0.0}\Downarrow}$ 0.60) | 12 | (${\color[rgb]{0.0,0.5,0.0}\definecolor[named]{pgfstrokecolor}{rgb}{0.0,0.5,0.0}\Uparrow}$ 2)
… | … |  |
25.07 | (${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\Uparrow}$ 1.26) | 51 | (${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\Downarrow}$ 8)
… | … |  |
27.41 | (${\color[rgb]{0.0,0.5,0.0}\definecolor[named]{pgfstrokecolor}{rgb}{0.0,0.5,0.0}\Downarrow}$ 1.27) | 64 | (${\color[rgb]{0.0,0.5,0.0}\definecolor[named]{pgfstrokecolor}{rgb}{0.0,0.5,0.0}\Uparrow}$ 5)
… | … |  |
| Value | Rank
---|---
0.09 | 21
0.10 | 32
0.09 | 19
0.09 | 14
0.10 | 29
0.13 | 70
0.13 | 77
0.15 | 92
0.13 | 69
0.08 | 9
0.10 | 33
0.16 | 104
… | …
0.06 | 2
… | …
0.20 | 124
… | …
Table 2: Final 2012 season rankings (and ranking shifts) in points allowed per
game by team’s defense, adjusted for offense’s ability to avoid turning the
ball over to opposing offense.

We additionally provide Table 2 to demonstrate the defensive ranking shifts
for the same season. Symmetrically, the strongest positive (negative) impacts
happened to the defensive units accompanied by offensive sides unable (able)
to consistently take care of the ball and avoid turning it over to the
opponent’s offense. For example, teams like Illinois, Florida State, South
Carolina, all sporting below-average offensive units in terms of ball
security, had their defensive units move up the ranking. On the other hand,
teams like Florida, Stanford, Pittsburgh, Kansas State, all experienced a
considerable drop defensive values due to being accompanied by offensive units
that took care of the ball (all in the Top-30).

####  3.2.4 Using adjusted statistics as features for game outcome prediction.

In addition to exploration of the most impactful complementary football
features, we compared predictive performances for statistics that were
adjusted in various ways. Specifically, we pitched the following four methods
against one another: no adjustment (classical averages), adjusting solely for
strength of schedule and homefield factor (SOS+HF), adjusting also for
efficiency-based complementary statistics (our main method), or for any
complementary statistic. The last approach is included for sanity check
purposes, to demonstrate how reverse causality hurts predictive performance
(see Section 3.1.2 for more details).

Table 3 shows the results of using logistic regression for binary game outcome
prediction (whether team won or lost), and classic linear regression for
modeling the score differential, respectively. As our five predictors we used
offensive and defensive points-per-game margins for both teams, along with the
homefield indicator. As prediction quality metrics, we used area under the
curve (AUROC) [13] for binary game outcome, and mean absolute error (MAE) for
score differential. To produce training errors, we utilized the full season
data for each year in the 2009-2019 span. For test errors, we used the first 9
weeks to train the model and subsequently yield predictions for the rest of
the season (using other cut-offs, e.g. 8, 10, 11 weeks, led to similar
results).

| Binary Game Outcome | Score Differential
---|---|---
Adjustment | (Won/Lost, AUROC) | (MAE)
| None
---
SOS+HF
SOS+HF & Efficiency Cmpl.
SOS+HF & Any Cmpl.
| Training | Test
---|---
0.790 (0.017) | 0.704 (0.028)
0.810 (0.017) | 0.715 (0.026)
0.808 (0.016) | 0.717 (0.033)
0.790 (0.019) | 0.707 (0.026)
| Training | Test
---|---
14.76 (0.57) | 17.39 (1.20)
13.78 (0.58) | 17.12 (1.18)
13.83 (0.57) | 17.13 (1.18)
14.58 (0.58) | 17.53 (0.93)

Table 3: Predicting binary game outcomes (won/lost) via logistic regression,
and numerical scoring differentials via linear regression, across 2009-2019
seasons by using variously adjusted points per game margin stats. Training
metric obtained from fitting the model to entire season, while test metric is
obtained from using only first 9 weeks to predict the rest of the season.

Performances are mostly similar across the methods, with all the average
metrics being within less than two standard deviations of one another.
Nonetheless, SOS+HF and our method consistently show better average
performances (higher AUROC, lower MAE values). When adjusting for any
complementary football features, rather than only efficiency-based ones, the
performance was as suboptimal as the classical averages, confirming the
dangers of reverse causality in case of plainly using the game totals.

##  4 Discussion and future work

We have applied several group penalization techniques in combination with
natural cubic splines to detect the most consistent features of American
college football in terms of the impact on the complementary side of your team
(e.g. your offense affecting your defense, or vice versa). Partially
constrained optimization was implemented to also guarantee the adjustment for
strength of the opponent and homefield factor. Issue of reverse causality for
internal dynamics of a football game was discussed, along with several
remedies, the most feasible being to focus solely on the efficiency-based
statistics (e.g. yards per play, points per possession, etc). Benefits of
using efficiency statistics were shown in evaluating game outcome predictions
based on various adjustments, where efficiency-based complementary statistics
outperformed arbitrary complementary statistics, with the latter suffering
from reverse causality more heavily.

Among the findings, we unsurprisingly showed the non-scoring turnovers
(possession percentage) forced by your defensive unit as the most consistent
feature in positively impacting the scoring by your offensive unit. This is
reasonable due to such turnovers giving the your offensive side an extra
possession along with the potential for a good starting field position. Hence,
the ranking shifts in points or touchdowns scored per game that resulted from
this adjustment tended to penalize offensive units complemented by a defensive
side that forced many non-scoring turnovers (conversely, in points or
touchdowns allowed, most penalized were defensive units complemented by
offensive sides that committed few turnovers).

On the other hand, team’s defensive statistic mostly affecting the total
yardage gained by the offensive unit was the yards allowed per play. This
effect was non-linear, with defenses allowing a moderate per-play yardage
(about 3 to 5) being less conducive to their complementary offensive units
gaining yards in that game, while defenses on more extreme ends of the
spectrum - allowing either tons (e.g. $\geq 6$) or very few (e.g. $\leq 2.5$)
yards per play - tend to give more opportunities for their offense to gain
yards. As thoroughly discussed in Section 3.2.2, this finding is intuitive due
to the moderate yards per play generally indicating the inability of your
defense to stop the opposing offense from controlling the game clock, hence
preventing your offense from gaining yards by keeping it on the sidelines.

In future work, we plan on gaining access to the play-by-play data and
incorporating it into the analysis for the purpose of alleviating the reverse
causality issue, resulting into more accurate evaluation of complementary
football feature effects. Although using solely efficiency-based statistics
leads to reasonable results, it constitutes only a partial solution for
scenarios of limited data access (e.g. game-by-game total statistics rather
than play-by-play) and could still suffer from endogeneity, which manifests in
lack of improvement over pure strength-of-schedule and homefield factor
(SOS+HF) adjustment for game outcome predictions. Lastly, one may consider
trying to find strong exogenous instrumental variables, which has been a
notoriously difficult task.

## Acknowledgments

The author is grateful to www.sports-reference.com/cfb/ and www.cfbstats.com/
for providing open access comprehensive game-by-game college football data.

## Declaration of interest

We confirm that there are no known conflicts of interest associated with this
publication and there has been no financial support that could have influenced
its outcome.

## Appendix

###  4.1 Feature sets.

####  4.1.1 All features (not only efficiency-based).

After disposing of strong pairwise correlations (e.g. $Off.Pts$, $Off.TD$,
$Off\\_w\\_ST.TD$ all have a correlation of $0.95+$ with the more all-
encompassing and ubiquitously used $Off\\_w\\_ST.Pts$, etc), perfect multi-
collinearities (e.g. instead of including passing attempts, completions and
incompletions, only retaining the individual breakdown stats of completions
and incompletions), unneeded statistic breakdowns (e.g. there is no difference
to the complementary side if it was an individual or assisted tackle, whether
the turnover was an interception or a fumble, etc), and all low-denominator
efficiency statistics (field goal, red zone, fourth down conversion
percentages, punt and return yards per attempt, etc), we were left with the
features in Table 4. Note that the overall turnovers by offense and offensive
special teams are not included because they cause multi-collinearity with non-
scoring turnovers and defensive (plus defensive special teams) touchdowns,
with that breakdown warranted due to difference of the impact on the
complementary side.

Name | Description
---|---
$Off.Tot$ | Total first downs
$Off.Cmp$ | Pass completions
$Off.Incmp$ | Pass incompletions
$Off.Pct$ | Pass completion percentage
$Off.Pass.Yds$ | Total passing yards
$Off.Pass.Avg$ | Average yards per pass attempt
$Off.Rush.Att$ | Rush attempts
$Off.Rush.Yds$ | Total rushing yards
$Off.Rush.Avg$ | Average yards per rush attempt
$Off.FGM$ | Field goals made
$Off.FGF$ | Field goals failed
$Off.X3rdDown..Conversions$ | Third downs converted
$Off.X3rdDown..Failed$ | Third downs failed
$Off.X3rdDown..Conversion..$ | Third down conversion percentage
$Off.X4thDown..Conversions$ | Fourth downs converted
$Off.X4thDown..Failed$ | Fourth downs failed
$Off.RedZone..TD$ | Red zone touchdowns
$Off.RedZone..FG$ | Red zone field goals
$Off.XPF$ | Extra points failed
$Off.ST.Punting.Punts$ | Punts
$Off.ST.Return.TD$ | Return touchdown (punt or kickoff)
$Off.ST.Return.Yds$ | Total return yards (punt or kickoff)
$Off.ST.Return.Ret$ | Total returns attempted (punt or kickoff)
$Def.Sfty$ | Safeties
$Def.QB.Hurries$ | Quarterback hurries
$Def.KicksPunts.Blocked$ | Kicks or punts blocked
$Def.Tackles.Tot$ | Total tackles
$Def.Tackles.Loss$ | Tackles for loss
$Def.Tackles.Sk$ | Sacks
$Def.Int.PD$ | Passes defended
$Def.Fumbles.FF$ | Fumbles forced
$Off\\_w\\_ST.Pts$ | Points scored
$Off\\_w\\_ST.Tot.TO.NonScoring$ | Non-scoring turnovers
$Def\\_w\\_ST.TD$ | Touchdowns scored (by defensive side )
$Def\\_w\\_ST.Yds$ | Yards gained (by defensive side)

Table 4: Glossary of all original complementary football features. $Off$
stands of ”offense-only”, $Def$ \- ”defense-only”, $Off.ST$ \- ”offensive
special teams”, $Off\\_w\\_ST$ \- ”offense with offensive special teams”,
$Def\\_w\\_ST$ \- ”defense with defensive special teams”.

####  4.1.2 Efficiency-based features.

Just like in case of the original feature set, we dispose of strong pairwise
correlations, perfect multi-collinearities, unneeded statistic breakdowns, and
all low-denominator efficiency statistics. As described in the main
manuscript, for each statistical category we had to decide on whether to use a
per-play or per-possession efficiency. We’ve selected to look at plays that
terminate a possession, e.g. turnover, points scored, punt, on per-possession
basis, making them into possession percentages. For example,
$Off\\_w\\_ST.NonScoring.TO.PossPct$ corresponds to the percentage of
possessions that ended in the offense turning the ball over to their defense
(as opposed to turning it over for an instant score by the opposing defense,
which is captured via $Def\\_ST.TD.PossPct$). On the other hand, virtually all
other metrics were treated on per-play basis, including yards or first downs
gained, overall tackles, tackles for loss, forced fumbles. Lastly, one can
notice such statistical categories as ”sacks to pass attempts ratio”, which
were defined in an effort to adjust for effects of pass-heavy play calling
being conducive to such defensive plays as sacks or quarterback hurries. While
it’s understood that sack (or quarterback hurry) is not recorded as a passing
attempt, and ideally one would need the number of quarterback dropbacks to
make it a proper efficiency percentage, but due to data limitations this was a
reasonable alternative. In similar vein, percentage of forced fumbles by
defense was only counted for rushing attempts and pass completions, meaning
the plays where fumble was possible. Table 5 demonstrates all the efficiency-
based features considered.

Name | Description
---|---
$Off.Pct$ | Pass completion percentage
$Off.Pass.Avg$ | Average yards per pass attempt
$Off.Rush.Avg$ | Average yards per rush attempt
$Off.Tot.Avg$ | Average yards per play
$Off.X1stDowns.Avg$ | Average first downs gained per play
$Off.X3rdDown..Conversion..$ | Third down conversion percentage
$Def.KicksPunts.Blocked.Pct$ | Percentage of kicks and punts blocked
$Def.Tackles.Tot.Avg$ | Tackles per play
$Def.Tackles.Loss.Avg$ | Tackles for loss per play
$Def.QB.Hurries.To.Pass.Att.Ratio$ | Ratio of QB hurries to pass attempts
$Def.Tackles.Sk.To.Pass.Att.Ratio$ | Ratio of sacks to pass attempts
$Def.Int.PD.Avg$ | Defended passes per pass attempt
$Def.Fumbles.FF.Avg$ | Forced fumbles per play
$Off\\_w\\_ST.Punt.Sfty.PossPct$ | Punt or safety
$Off\\_w\\_ST.Tot.TO.NonScoring.PossPct$ | Non-scoring turnovers
$Off\\_w\\_ST.Pts.PerPos$ | Points scored
$Def\\_w\\_ST.TD.PossPct$ | Touchdowns scored (by defensive side)

Table 5: Glossary of all efficiency-based complementary football features.
$Off$ stands of ”offense-only”, $Def$ \- ”defense-only”, $Off.ST$ \-
”offensive special teams”, $Off\\_w\\_ST$ \- ”offense with offensive special
teams”, $Def\\_w\\_ST$ \- ”defense with defensive special teams”, $PossPct$ \-
”possession percentage”, $Avg$ \- per-play average, $PerPos$ \- ”per
possession”.

## References

  * [1]
  * Albert et al. [2005] Albert, J., Bennett, J. and Cochran, J. J. [2005], Anthology of statistics in sports, SIAM.
  * Bartko [1966] Bartko, J. J. [1966], ‘The intraclass correlation coefficient as a measure of reliability’, Psychological reports 19(1), 3–11.
  * Boshnakov et al. [2017] Boshnakov, G., Kharrat, T. and McHale, I. G. [2017], ‘A bivariate weibull count model for forecasting association football scores’, International Journal of Forecasting 33(2), 458–466.
  * Bound et al. [1995] Bound, J., Jaeger, D. A. and Baker, R. M. [1995], ‘Problems with instrumental variables estimation when the correlation between the instruments and the endogenous explanatory variable is weak’, Journal of the American statistical association 90(430), 443–450.
  * Bowden and Turkington [1990] Bowden, R. J. and Turkington, D. A. [1990], Instrumental variables, number 8, Cambridge university press.
  * Breheny and Breheny [2021] Breheny, P. and Breheny, M. P. [2021], ‘Package ‘grpreg”.
  * Browne [2000] Browne, M. W. [2000], ‘Cross-validation methods’, Journal of mathematical psychology 44(1), 108–132.
  * Duncan et al. [2004] Duncan, G. J., Magnuson, K. A. and Ludwig, J. [2004], ‘The endogeneity problem in developmental studies’, Research in human development 1(1-2), 59–80.
  * Edwards and Archambault [1979] Edwards, J. and Archambault, D. [1979], ‘The home field advantage’, Sports, games, and play: Social and psychological viewpoints pp. 409–438.
  * Efron [1994] Efron, B. [1994], ‘Missing data, imputation, and the bootstrap’, Journal of the American Statistical Association 89(426), 463–475.
  * Fan and Li [2001] Fan, J. and Li, R. [2001], ‘Variable selection via nonconcave penalized likelihood and its oracle properties’, Journal of the American statistical Association 96(456), 1348–1360.
  * Hanley and McNeil [1982] Hanley, J. A. and McNeil, B. J. [1982], ‘The meaning and use of the area under a receiver operating characteristic (roc) curve.’, Radiology 143(1), 29–36.
  * Harrison [2020] Harrison, J. [2020], RSelenium: R Bindings for ’Selenium WebDriver’.  R package version 1.7.7.
https://CRAN.R-project.org/package=RSelenium

  * Harville [1977] Harville, D. [1977], ‘The use of linear-model methodology to rate high school or college football teams’, Journal of the American Statistical Association 72(358), 278–289.
  * Hastie et al. [2009] Hastie, T., Tibshirani, R. and Friedman, J. [2009], The elements of statistical learning: data mining, inference and prediction, 2 edn, Springer.
http://www-stat.stanford.edu/~tibs/ElemStatLearn/

  * Huang et al. [2012] Huang, J., Breheny, P. and Ma, S. [2012], ‘A selective review of group selection in high-dimensional models’, Statistical science: a review journal of the Institute of Mathematical Statistics 27(4).
  * Karlis and Ntzoufras [2003] Karlis, D. and Ntzoufras, I. [2003], ‘Analysis of sports data by using bivariate poisson models’, Journal of the Royal Statistical Society: Series D (The Statistician) 52(3), 381–393.
  * Kmenta [1986] Kmenta, J. [1986], ‘Elements of econometrics. 2d ed. new york city’.
  * Leng et al. [2006] Leng, C., Lin, Y. and Wahba, G. [2006], ‘A note on the lasso and related procedures in model selection’, Statistica Sinica pp. 1273–1284.
  * Neath and Cavanaugh [2012] Neath, A. A. and Cavanaugh, J. E. [2012], ‘The bayesian information criterion: background, derivation, and applications’, Wiley Interdisciplinary Reviews: Computational Statistics 4(2), 199–203.
  * R Core Team [2020] R Core Team [2020], ‘R: A language and environment for statistical computing. R Foundation for Statistical Computing. Version 4.0.3’.
https://www.R-project.org/

  * Rutz and Watson [2019] Rutz, O. J. and Watson, G. F. [2019], ‘Endogeneity and marketing strategy research: an overview’, Journal of the Academy of Marketing Science 47(3), 479–498.
  * Tibshirani [1996] Tibshirani, R. [1996], ‘Regression shrinkage and selection via the lasso’, Journal of the Royal Statistical Society. Series B (Methodological) pp. 267–288.
  * Vergin and Sosik [1999] Vergin, R. C. and Sosik, J. J. [1999], ‘No place like home: an examination of the home field advantage in gambling strategies in nfl football’, Journal of Economics and Business 51(1), 21–31.
  * Wang and Leng [2008] Wang, H. and Leng, C. [2008], ‘A note on adaptive group lasso’, Computational statistics & data analysis 52(12), 5277–5286.
  * Wickham [2019] Wickham, H. [2019], rvest: Easily Harvest (Scrape) Web Pages.  R package version 0.3.5.
https://CRAN.R-project.org/package=rvest

  * Yuan and Lin [2006] Yuan, M. and Lin, Y. [2006], ‘Model selection and estimation in regression with grouped variables’, Journal of the Royal Statistical Society: Series B (Statistical Methodology) 68(1), 49–67.
  * Zhang [2010] Zhang, C.-H. [2010], ‘Nearly unbiased variable selection under minimax concave penalty’, The Annals of statistics 38(2), 894–942.
  * Zou [2006] Zou, H. [2006], ‘The adaptive lasso and its oracle properties’, Journal of the American statistical association 101(476), 1418–1429.


solutionSolutionsolutionfile solutionfile[AoI_v2]

# Age of Information in Prioritized Random Access

Khac-Hoang Ngo, Giuseppe Durisi, and Alexandre Graell i Amat This project has
received funding from the European Union’s Horizon 2020 research and
innovation programme under the Marie Skłodowska-Curie grant agreement No
101022113. Department of Electrical Engineering, Chalmers University of
Technology, 41296 Gothenburg, Sweden

###### Abstract

Age of information (AoI) is a performance metric that captures the freshness
of status updates. While AoI has been studied thoroughly for point-to-point
links, the impact of modern random-access protocols on this metric is still
unclear. In this paper, we extend the recent results by Munari to prioritized
random access where devices are divided into different classes according to
different AoI requirements. We consider the irregular repetition slotted ALOHA
protocol and analyze the AoI evolution by means of a Markovian analysis
following similar lines as in Munari (2021). We aim to design the protocol to
satisfy the AoI requirements for each class while minimizing the power
consumption. To this end, we optimize the update probability and the degree
distributions of each class, such that the probability that their AoI exceeds
a given threshold lies below a given target and the average number of
transmitted packets is minimized.

##  I Introduction

The Internet of Things (IoT) foresees a very large number of devices, which we
will refer to as users, to be connected and exchange data in a sporadic and
uncoordinated manner. This has led to the development of modern random access
protocols [1]. In most of these protocols, the users transmit multiple copies
of their packets to create time diversity, and the receiver employs successive
interference cancellation (SIC) to decode. In particular, in the irregular
repetition slotted ALOHA (IRSA) protocol [2], the users draw the number of
copies from a degree distribution and transmit the copies in randomly chosen
slots of a fixed-length frame. A common design goal is to minimize the packet
loss rate (PLR), thus maximizing the chance to deliver packets to the receiver
successfully.

In many IoT applications, it is becoming increasingly important to deliver
packets successfully and to guarantee the timeliness of those packets
simultaneously. Examples include sensor networks, vehicular tracking, and
health monitoring. In these delay-sensitive applications, the packets carry
critical status updates that are required to be fresh. The age of information
(AoI) metric (see, e.g., [3] and references therein) has been introduced
precisely to account for the freshness of the status updates. It captures the
offset between the generation of a packet and its observation time. In [4],
the AoI in a system where independent devices send status updates through a
shared queue was analyzed. The AoI has been used as a performance metric to
design status update protocols in, e.g., [5, 6]. The first analytical
characterization of the AoI for a class of modern random access, namely IRSA,
has been recently reported in [7].

Since IoT devices are mostly battery-limited, their power consumption should
be minimized. By assuming that each packet transmission consumes a fixed
amount of energy, we can use the average number of transmitted packets per
slot as a proxy of the power consumption. When status updates are conveyed via
an IRSA protocol, the number of transmitted packets per user depends both on
the user activity, i.e., on how often the user has an update to transmit, and
on the degree distribution assigned to the user. This leads to a tension
between minimizing the AoI and minimizing the number of packets. Too sporadic
user activity leads to stale updates, but too frequent updates lead to channel
congestion and update failure. Furthermore, degree distributions with low
degrees lead to a low number of transmitted packets but high PLR, which
results in larger AoI, while degree distributions with high degrees achieve
low PLR at the cost of a larger number of transmitted packets. Therefore, the
user update probability and the degree distribution need to be carefully
selected.

In this paper, we consider an IoT monitoring system where users attempt to
deliver timely status updates to a receiver following the IRSA protocol. We
assume that the users are heterogeneous and their updates require different
levels of freshness. Accordingly, users are divided into different classes,
each with a different AoI requirement. Following similar lines as in [7], we
analyze the AoI evolution by means of a Markovian analysis and derive the age-
violation probability (AVP), i.e., the probability that the AoI exceeds a
certain threshold, for each class. We study the trade-off between the AVP and
the number of transmitted packets by investigating the impact of the update
probability and the degree distributions. Since the PLR of IRSA and, hence,
the AVP are not known in closed form, we propose an easy-to-compute PLR
approximation, which leads to an accurate approximation of the AVP. Our PLR
approximation is based on density evolution (DE) [2] and on existing PLR
approximations in the error-floor region [8] and the waterfall region [9]. We
jointly optimize the update probability and the degree distributions for each
class to minimize the number of transmitted packets while guaranteeing that
the AVP of each class lies below a given target. Our simulation results show
that the number of transmitted packets can be significant reduced with
optimized irregular degree distributions, compared to regular distributions.
Our experiments also suggest that using degrees up to $3$ is sufficient for a
setting where there are two classes containing respectively $800$ and $3200$
users, the framelength is $100$ slots, and the AoI of class-$1$ users and
class-$2$ users exceeds a threshold $7.5\\!\times\\!10^{4}$ and
$4.5\\!\times\\!10^{4}$ with probability as low as $10^{-5}$ and $10^{-3}$,
respectively.

##  II System Model and Problem Formulation

We consider a system with $U$ users attempting to deliver timestamped status
updates to a receiver through a wireless channel. Time is slotted and each
update is transmitted in a slot. We let the slot length be $1$ without loss of
generality. Each user belongs to one of $K$ classes with different AVP
requirements. Let $U_{k}$ be the number of users in class $k$, $k\in[K]$.111We
use $[m:n]$ to denote the set of integers from $m$ to $n$, and $[n]=[1:n]$. We
define the fraction of class-$k$ users as $\gamma_{k}=U_{k}/U$. We assume that
a class-$k$ user has a new update in each slot with probability $\mu_{k}$
independently of the other users. We further assume that slots containing a
single packet always lead to successful decoding, whereas slots containing
multiple packets (or, more specifically, unresolved collisions after SIC) lead
to decoding failures.

###  II-A Irregular Repetition Slotted ALOHA

We assume that the system operates according to the IRSA protocol. Time is
divided into frames of $M$ slots and users are frame- and slot-synchronous. A
user may generate more than one update during a frame, but only the latest
update is transmitted in the next frame. An active user in class $k$ sends
${\mathsf{L}}_{k}$ identical replicas of its latest update in
${\mathsf{L}}_{k}$ slots chosen uniformly without replacement from the $M$
available slots. The number ${\mathsf{L}}_{k}$ is called the degree of the
transmitted packet. It follows a class-dependent probability distribution
$\\{\Lambda^{(k)}_{\ell}\\}$ where
$\Lambda^{(k)}_{\ell}\\!=\Pr[{\mathsf{L}}_{k}\\!=\\!\ell]$. We write this
distribution using a polynomial notation as
$\Lambda^{(k)}(x)=\sum_{\ell=0}^{d}\Lambda^{(k)}_{\ell}x^{\ell}$ where $d$ is
the maximum degree. Note that $\\{\Lambda^{(k)}\\}$ may contain degree $0$.
When ${\mathsf{L}}_{k}=0$, the user discards the update. Upon successfully
receiving an update, the receiver is assumed to be able to determine the
position of its replicas. In practice, this can be done by including in the
header of the packet containing each update a pointer to the position of its
replicas. The receiver employs a SIC decoder. It seeks slots containing a
single packet, decodes the packet, then locates and removes the replicas.
These steps are repeated until no slots with a single packet can be found.

Note that a user in class $k$ has a new update in a frame with probability
$\sigma_{k}=1-(1\\!-\\!\mu_{k})^{M}$. Therefore, the number of class-$k$ users
transmitting over a frame is a binomial random variable of parameters
$(U_{k},\sigma_{k})$ with expected value $U_{k}\sigma_{k}$. The average
channel load of class $k$ is given by $G_{k}=U_{k}\sigma_{k}/M.$ The overall
average channel load is $G=\sum_{k=1}^{K}G_{k}$.

The average number of packets transmitted by a class-$k$ user per slot is
$\Phi_{k}=\sigma_{k}\dot{\Lambda}^{(k)}(1)/M,$ where $\dot{\Lambda}^{(k)}(x)$
denotes the first-order derivative of $\Lambda^{(k)}(x)$. The total average
number of transmitted packets per slot is given by

$\displaystyle\Phi=\sum_{k=1}^{K}U_{k}\Phi_{k}=\sum_{k=1}^{K}G_{k}\dot{\Lambda}^{(k)}(1).$
(1)

We use $\Phi$ as a proxy of the total power consumption.

###  II-B Age of Information

We define the AoI for user $i$ at slot $n$ as $\delta_{i}(n)=n-t_{i}(n)$,
where $t_{i}(n)$ denotes the timestamp of the last received update from user
$i$ as of slot $n$. Since the AoIs of users in the same class are
stochastically equivalent, we denote a representative of the AoIs of class-$k$
users as $\delta^{(k)}(n)$. The AoI grows linearly with time and is reset at
the end of a frame only when a new update is successfully decoded. We are
interested in the value of the AoI at the end of a generic frame
$j\in\mathbb{N}_{0}$. We will refer to this quantity simply as AoI hereafter.
For class $k$, this quantity is given by $\delta^{(k)}(jM)+M$. We define the
AVP as the probability that the AoI exceeds a certain threshold $\theta$ at
steady state. Specifically, the AVP for class $k$ is defined as

$\displaystyle\zeta^{(k)}(\theta)=\lim\limits_{j\to\infty}\Pr[\delta^{(k)}(jM)+M>\theta].$
(2)

We shall see in the next section that the AoI process is ergodic Markovian,
thus the limit in (2) exists. We consider the requirement that the AoI at
steady state of class $k$ exceeds a threshold $\theta_{k}$ with probability no
larger than $\epsilon_{k}$:

$\displaystyle\zeta^{(k)}(\theta_{k})\leq\epsilon_{k},\quad k\in[K].$ (3)

###  II-C Problem Formulation

Our goal is to design the update probabilities $\\{\mu_{k}\\}$ and the degree
distributions $\\{\Lambda^{(k)}\\}$ such that the AoI requirements in (3) are
satisfied, while the number of packets $\Phi$ is minimized:

$\displaystyle\operatorname*{minimize}_{\\{\mu_{k},\Lambda^{(k)}(x)\\}_{k=1}^{K}}~{}$
$\displaystyle\Phi\qquad\text{~{}subject to~{}}\eqref{eq:AoI_requirement}.$
(4)

##  III AoI Analysis

###  III-A Current AoI

Let $P^{(k)}$ denote the PLR of a class-$k$ user. It is convenient to denote
by $\xi_{k}=\sigma_{k}(1-P^{(k)})$ the probability that the AoI
$\delta^{(k)}(n)$ is reset. Also, let ${\mathsf{B}}_{k}\in[M]$ denote the
number of slots between the generation of a packet of a class-$k$ user and the
start of the subsequent frame, when the user can access the channel. It has
probability mass function
$\Pr[{\mathsf{B}}_{k}=b]=\mu_{k}(1-\mu_{k})^{b-1}/\sigma_{k},$ where the
numerator is the probability that the user has generated an update for the
last time $b$ slots before the end of a frame, and the denominator is the
probability that at least one update is generated during the frame. Whenever
an update is successfully decoded, the current AoI is reset to
${\mathsf{B}}_{k}+M\in[M+1:2M]$.

In what follows, it will be convenient to decompose an arbitrary integer $n$
as $n=\alpha_{n}+M\beta_{n}$, where $\alpha_{n}=n\\!\mod M$ and
$\beta_{n}=\lfloor n/M\rfloor$. Using this decomposition, we can write
$\delta^{(k)}(n)=\delta^{(k)}(M\beta_{n})+\alpha_{n},$ where the first term on
the right-hand side captures the age at the beginning of the current frame,
and the second is the offset from the start of the current frame up to the
observation time $n$. We set $n=0$ right after the reception of the first
update. Thus, the initial AoI is in $[M+1:2M]$, and $\delta^{(k)}(n)\geq
M+1,\forall n$. Therefore, the evolution of the AoI of a class-$k$ user is
fully characterized by the discrete-time, discrete-valued stochastic process

$\displaystyle\Omega^{(k)}_{j}=\delta^{(k)}(jM)-(M+1),\quad
j\in\mathbb{N}_{0},k\in[K],$ (5)

where $j$ is the frame index. Since each user operates independently over
successive frames, $\Omega^{(k)}_{j}$ is a Markovian process across $j$. The
one-step transition probabilities
$q^{(k)}_{n_{1},n_{2}}=\Pr\big{[}\Omega^{(k)}_{j+1}\\!=\\!n_{2}\,|\,\Omega^{(k)}_{j}\\!=\\!n_{1}\big{]}$
are given by

$\displaystyle
q^{(k)}_{n_{1},n_{2}}=\begin{cases}\xi_{k}\Pr[{\mathsf{B}}_{k}=n_{2}+1],&\text{for~{}}n_{2}\in[0:M-1],\\\
1-\xi_{k},&\text{for~{}}n_{2}=n_{1}+M,\\\ 0,&\text{otherwise.}\end{cases}$ (6)

To verify (6), note that $\Omega^{(k)}_{j}$ is reset with probability
$\xi_{k}$, and in this case, it is reset to a value $n_{2}+1$, where
$n_{2}\in[0:M-1]$, with probability $\Pr[{\mathsf{B}}_{k}=n_{2}+1]$. With
probability $1-\xi_{k}$, the variable $\Omega^{(k)}_{j}$ is simply incremented
by the framelength $M$.

We start with the following observation.

######  Proposition 1.

The stochastic process $\Omega^{(k)}_{j}$ is ergodic, and has steady-state
distribution

$\displaystyle\pi^{(k)}_{w}=\xi_{k}\left(1-\xi_{k}\right)^{\beta_{w}}\Pr[{\mathsf{B}}_{k}=\alpha_{w}+1],\quad
w\in\mathbb{N}_{0}.$ (7)

###### Proof.

The proof follows directly from the proof of the single-class case in [7,
Prop. 1]. ∎

It follows from Proposition 1 that the limit in (2) exists.

###  III-B Age-Violation Probability

It follows from (5) and (2) that
$\zeta^{(k)}(\theta)=\Pr[\Omega^{(k)}>\theta-2M-1],$ where the random variable
$\Omega^{(k)}$ has steady-state distribution $\\{\pi_{w}^{(k)}\\}$. The
following result holds.

######  Proposition 2.

The AVP is given by

$\\!\\!\\!\\!\\!\zeta^{(k)}(\theta)=\\\
\begin{cases}(1-\xi_{k})^{\beta_{\theta-2M}}\left[1-\frac{1-(1-\mu_{k})^{1+\alpha_{\theta-2M}}}{\sigma_{k}}\xi_{k}\right],&\text{for~{}}\theta>2M,\\\
1,&\text{otherwise}.\end{cases}$ (8)

###### Proof.

The proof follows similar steps as the proof of the single-class case in [7,
Prop. 3]. ∎

######  Example 1.

Consider a system with $U=4000$ users, framelength $M=100$, $K=2$ classes with
fractions $(\gamma_{1},\gamma_{2})=(0.2,0.8)$, AoI thresholds
$(\theta_{1},\theta_{2})=(7.5\times 10^{4},4.5\times 10^{4})$, and target AVPs
$(\epsilon_{1},\epsilon_{2})=(10^{-4},10^{-2})$. We further assume that
$\mu_{1}\\!=\\!\mu_{2}\\!=\\!\mu$. Thus,
$\Phi=\frac{U(1-(1-\mu)^{M})}{M}\sum_{u=1}^{U}\gamma_{k}\dot{\Lambda}^{(k)}(1)$,
which increases with $\mu$. We evaluate the AVPs for this scenario and plot
them as functions of $\Phi$ in Fig. 1 for three sets of regular degree
distributions, namely,
$\Lambda^{(1)}(x)=\Lambda^{(2)}(x)\in\\{x,x^{2},x^{3}\\}$. The PLR is computed
numerically. We vary $\Phi$ by varying $\mu$. Some remarks are in order.

  * •

For each class, the AVP first decreases and then increases with $\Phi$.
Indeed, when $\mu$ is low, collisions are unlikely. Although the updates are
successfully received with high probability, the sporadicity of the updates
entails a high AoI. When $\mu$ is high, users transmit frequently, and updates
fail with high probability due to collision, entailing a high AoI.

  * •

The target AVPs $(\epsilon_{1},\epsilon_{2})=(10^{-4},10^{-2})$ are not met
when $\Lambda^{(1)}(x)=\Lambda^{(2)}(x)=x$. The distributions
$\Lambda^{(1)}(x)\\!=\\!\Lambda^{(2)}(x)\\!=\\!x^{2}$ satisfy these
requirements with a minimum number of packets per slot $\Phi\approx 1.09$. The
distributions $\Lambda^{(1)}(x)\\!=\\!\Lambda^{(2)}(x)\\!=\\!x^{3}$ require a
higher $\Phi$ to achieve the same requirements, but can yield a reduction of
the AVP. For example, the more stringent requirements
$(\epsilon_{1},\epsilon_{2})=(10^{-5},10^{-3})$ can be met with about $1.92$
packets/slot. In general, distributions with low degrees can achieve mild AoI
requirements with a low $\Phi$, while higher degrees are needed to achieve
more stringent requirements.

The observations in Example 1 reveal the existence of a trade-off in the
choice of $\\{\mu_{k}\\}$ and $\\{\Lambda^{(k)}\\}$ to satisfy the AoI
requirements while minimizing $\Phi$.

$0$$0.5$$1$$1.5$$2$$2.5$$3$$10^{-6}$$10^{-5}$$10^{-4}$$10^{-3}$$10^{-2}$$10^{-1}$$10^{0}$$1.09$
packets/slot$1.92$packets/slotClass $2$Class
$1$$\epsilon_{1}$$\epsilon_{2}$Average number of packets per slot $\Phi$Age-
violation probability
$\zeta^{(k)}(\theta_{k})$$\Lambda^{(1)}(x)\\!=\\!\Lambda^{(2)}(x)\\!=\\!x$$\Lambda^{(1)}(x)\\!=\\!\Lambda^{(2)}(x)\\!=\\!x^{2}\\!$$\Lambda^{(1)}(x)\\!=\\!\Lambda^{(2)}(x)\\!=\\!x^{3}\\!$
Figure 1: The AVPs $\zeta^{(1)}(7.5\times 10^{4})$ and $\zeta^{(2)}(4.5\times
10^{4})$ vs. $\Phi$ for the scenario in Example 1 with
$(\epsilon_{1},\epsilon_{2})=(10^{-4},10^{-2})$. We consider three sets of
regular degree distributions, namely,
$\Lambda^{(1)}(x)=\Lambda^{(2)}(x)\in\\{x,x^{2},x^{3}\\}$.

##  IV Packet Loss Rate Approximation

The PLR for class-$k$ users can be derived as [10, Eq. (2)]

$P^{(k)}=\sum_{\ell=0}^{d}\Lambda_{\ell}^{(k)}P_{\ell},$ (9)

where $P_{\ell}$ is the probability that a degree-$\ell$ user (of any class)
is not resolved. The probability $P_{\ell}$ is determined by the overall
channel load $G$ and the average degree distribution
$\Lambda(x)=\sum_{\ell=0}^{d}\Lambda_{\ell}x^{\ell}$ with
$\Lambda_{\ell}=\sum_{k=1}^{K}\gamma_{k}\Lambda^{(k)}_{\ell}$, $\ell\in[0:d]$.
If $\Lambda_{0}>0$, then $P_{0}=1$ and $P_{\ell}$, $\ell\geq 1$, is the
probability that a degree-$\ell$ user is not resolved in a system with channel
load $\bar{G}=G(1-\Lambda_{0})$ and degree distribution
$\bar{\Lambda}(x)\\!=\\!\sum_{\ell=1}^{d}\bar{\Lambda}_{\ell}x^{\ell}$ with
$\bar{\Lambda}_{\ell}\\!=\\!\frac{1}{1-\Lambda_{0}}\sum_{k=1}^{K}\gamma_{k}\Lambda_{\ell}^{(k)}$.
Therefore, we assume without loss of generality that $\Lambda_{0}=0$ in the
remainder of the section. The PLR is not known in closed form in general, but
can be computed numerically. However, since the optimization (4) requires
repeated evaluation of the AVP, and thus of the PLR, simulation-based PLR
computation becomes inefficient. Therefore, we seek an easy-to-compute
approximation of the PLR that leads to an accurate approximation of the AVP.

The SIC process of IRSA is equivalent to graph-based iterative erasure
decoding of low-density parity-check (LDPC) codes. In the asymptotic regime
where $M\to\infty$, $P_{\ell}$ can be evaluated using DE as [2]

$P_{\ell,\text{DE}}=\lim\limits_{i\to\infty}(\eta_{i})^{\ell},$ (10)

where $\eta_{i}$ is the probability that an edge connected to a degree-$\ell$
user remains unknown in the decoding process. It can be computed in an
iterative manner as $\eta_{0}=1$,
$\eta_{i}=1-\exp(-G\dot{\Lambda}(\eta_{i-1}))$ where
$\dot{\Lambda}(x)=\mathop{}\mathrm{d}\Lambda(x)/\mathop{}\mathrm{d}x$. As
$M\to\infty$, $P_{\ell}=P_{\ell,\text{DE}}$ and the PLR $P^{(k)}$ drops at a
certain threshold value as the channel load $G$ decreases. That is, all but a
vanishing fraction of the class-$k$ users are resolved if the channel load is
below the decoding threshold. According to [10, Prop. 1], the thresholds for
all classes coincide and can be obtained by means of DE as the largest value
$G^{*}$ of $g$ such that $\nu>1-\exp(-g\dot{\Lambda}(\nu))$ for all
$\nu\in(0,1]$.

In the finite-framelength regime, the PLR is typically characterized by two
regions: a waterfall (WF) region near the DE threshold where the PLR decreases
sharply, and an error-floor (EF) region where the PLR flattens. In the WF
region, according to [9], the PLR can be approximated based on the finite-
length scaling of the frame-error rate of LDPC codes [11]. Specifically, the
overall PLR $\sum_{\ell=0}^{d}\Lambda_{\ell}P_{\ell}$ (averaged over classes)
can be approximated by

$\displaystyle P_{\text{WF}}=P_{G\to
1}Q\bigg{(}\frac{\sqrt{M}(G^{*}-\beta(\Lambda)M^{-2/3}-G)}{\sqrt{\alpha^{2}(\Lambda)+G(1-MG/U)}}\bigg{)}$
(11)

where $P_{G\to 1}$ is the PLR in the limit $G\to 1$ computed via DE,
$Q(\cdot)$ is the Gaussian Q-function, and
$\\{\alpha(\Lambda),\beta(\Lambda)\\}$ are scaling parameters computed as
specified in [12]. In the EF region, the PLR can be approximated using the
method proposed in [8]. In this region, decoding failures are mainly caused by
harmful structures in the corresponding bipartite graph, referred to as
stopping sets. A connected bipartite graph ${\mathcal{S}}$ is a stopping set
if all check nodes in ${\mathcal{S}}$ have a degree larger than one. By
enumerating the stopping sets, we can approximate $P_{\ell}$ by

$P_{\ell,\text{EF}}=\frac{(U\\!-\\!1)!}{\Lambda_{\ell}}\sum_{{\mathcal{S}}\in{\mathcal{A}}}\\!\frac{v_{\ell}({\mathcal{S}})c({\mathcal{S}})\binom{M}{\psi({\mathcal{S}})}\\!}{(U\\!-\\!v({\mathcal{S}}))!}\prod_{j=1}^{d}\\!\binom{M}{j}^{\\!\\!-v_{j}({\mathcal{S}})}\\!\frac{\Lambda_{j}^{v_{j}({\mathcal{S}})}}{v_{j}({\mathcal{S}})!\\!},$
(12)

where ${\mathcal{A}}$ is the set of considered stopping sets,
$v({\mathcal{S}})$ and $\psi({\mathcal{S}})$ are the number of variable nodes
and check nodes in ${\mathcal{S}}$, respectively, $v_{j}({\mathcal{S}})$ is
the number of degree-$j$ variable nodes in ${\mathcal{S}}$, and
$c({\mathcal{S}})$ is the number of graphs isomorphic with ${\mathcal{S}}$.

In [7], the PLR is approximated as

$\displaystyle P_{\ell}\approx P_{\text{WF}}+P_{\ell,\text{EF}}$ (13)

for the single-class case.222In the single-class case, this means that the
overall PLR is approximated by $P_{\rm
WF}+\sum_{\ell=0}^{d}\Lambda_{\ell}P_{\ell,{\rm EF}}$. In our paper, it is
more convenient to write the approximation in terms of $P_{\ell}$. It was
shown to be accurate for degree distributions with degrees at least $3$ (see,
e.g., [7, Fig. 4]). These degree distributions are typically considered when
the design goal is to minimize the EF or maximize the decoding threshold. In
our setting, however, it is of interest to consider degree distributions with
lower degrees to reduce $\Phi$.

In Figs. 2(a) and 2(b), we investigate the tightness of the approximations
(13) and $P_{\ell}\approx P_{\ell,\text{DE}}$ for the setup in Example 1 and
some degree distributions with degrees $1$ and $2$. Note that an accurate PLR
approximation in the WF region is crucial for the computation of the AVP,
whereas the AVP is insentitive to low values of the PLR in the EF region where
update sporadicity is the dominating factor. As shown in Fig. 2(a) for the
degree distributions $\Lambda^{(1)}(x)=\Lambda^{(2)}(x)=0.5x^{2}+0.5x^{3}$,
the approximation (13) is loose in the WF region. The reason is that the
finite-length scaling leading to (11) is not guaranteed to hold when the
bipartite graph contains degree-$2$ variable nodes [11]. The situation is even
worse when degree-$1$ users are present: $P_{\text{WF}}$ is near $P_{G\to 1}$
for all channel load. This makes the approximation (13) inaccurate, as shown
for the distributions $\\{\Lambda^{(1)}(x)=\ 0.7x+0.3x^{3}$,
$\Lambda^{(2)}(x)=0.7x^{2}+0.3x^{3}\\}$ in Fig. 2(b). On the other hand,
$P_{\ell,\text{DE}}$ is an accurate approximation of $P_{\ell}$ for large
values of $\Phi$ corresponding to $G>G^{*}$, although $P_{\ell,\text{DE}}$ is
much lower than $P_{\ell}$ when $G<G^{*}$.

In Figs. 2(c) and 2(d), we show the AVP computed with different approximations
of the PLRs. For both sets of degree distributions, setting
$P_{\ell}\\!\approx\\!P_{\ell,\text{DE}}$ yields an accurate approximation of
the AVP when $G>G^{*}$. When class-$1$ users are present as in Fig. 2(d), we
have $G^{*}=0$ and the AVP approximation obtained by setting
$P_{\ell}\\!\approx\\!P_{\ell,\text{DE}}$ is accurate for all $G>0$
(equivalently $\Phi>0$). For moderate values of $G$ corresponding to the WF
region, if there are both degree $2$ and higher degrees, the approximation
$P_{\ell}\\!\approx\\!P_{\ell,\text{DE}}$ leads to an optimistic approximation
of the AVP while the approximation (13) is pessimistic in the WF region, as
shown in Fig. 2(c). In this case, one needs to balance between these two
approximations. Our experiments suggest that using
$P_{\ell}\\!\approx\\!P_{\ell,\text{DE}}$ for $\ell=2$ and (13) for $\ell>2$,
which leads to the dashed-dotted green line in Fig. 2(c), results in an
accurate approximation of the AVP.

From the above observations, we propose the following heuristic PLR
approximation. For $G>G^{*}$, we set $P_{\ell}\approx P_{\ell,\text{DE}}$ for
all $\ell$. For $G\leq G^{*}$, we set $P_{\ell}\approx P_{\ell,\text{DE}}$ for
$\ell\leq 2$ and $P_{\ell}\approx P_{\text{WF}}+P_{\ell,\text{EF}}$ for
$\ell>2$. We summarize the proposed approximation as

$\displaystyle
P_{\ell}\approx\begin{cases}P_{\ell,\text{DE}},&\text{if~{}$\ell\leq 2$ or
$G>G^{*}$},\\\ P_{\text{WF}}+P_{\ell,\text{EF}},&\text{if~{}$\ell>2$ and
$G\leq G^{*}$}.\end{cases}$ (14)

$0$$1$$2$$3$$10^{-5}$$10^{-4}$$10^{-3}$$10^{-2}$$10^{-1}$$10^{0}$$\Phi$
(packets/slot)Packet loss
rateSimulation$P_{\ell}\\!\approx\\!P_{\text{WF}}\\!+\\!P_{\ell,\text{EF}}\\!\\!$$P_{\ell}\approx
P_{\ell,\text{DE}}$With (14)$0$$0.4$$0.8$$1.2$Overall channel load $G$ (a) PLR
for $\Lambda^{(1)}(x)=\Lambda^{(2)}(x)=0.5x^{2}+0.5x^{3}$ ($G^{*}\approx
0.792$)

$0$$1$$2$$3$$10^{-5}$$10^{-4}$$10^{-3}$$10^{-2}$$10^{-1}$$10^{0}$Class
$2$Class $1$$\Phi$ (packets/slot)Packet loss
rateSimulation$P_{\ell}\\!\approx\\!P_{\text{WF}}\\!+\\!P_{\ell,\text{EF}}\\!\\!$$P_{\ell}\approx
P_{\ell,\text{DE}}$With (14)$0$$0.46$$0.93$$1.39$Overall channel load $G$ (b)
PLR for $\Lambda^{(1)}(x)\\!=0.7x\\!+0.3x^{3}$,
$\Lambda^{(2)}(x)=0.7x^{2}\\!+\\!0.3x^{3}$ ($G^{*}=0$)

$0$$1$$2$$3$$10^{-6}$$10^{-5}$$10^{-4}$$10^{-3}$$10^{-2}$$10^{-1}$$10^{0}$Class
$1$Class $2$$\Phi$ (packets/slot)Age-violation
probability$0$$0.4$$0.8$$1.2$Overall channel load $G$ (c) AVP for
$\Lambda^{(1)}(x)=\Lambda^{(2)}(x)=0.5x^{2}+0.5x^{3}$ ($G^{*}\approx 0.792$)

$0$$1$$2$$3$$10^{-6}$$10^{-5}$$10^{-4}$$10^{-3}$$10^{-2}$$10^{-1}$$10^{0}$Class
$1$Class $2$$\Phi$ (packets/slot)Age-violation
probability$0$$0.46$$0.93$$1.39$Overall channel load $G$ (d) AVP for
$\Lambda^{(1)}(x)\\!=0.7x\\!+0.3x^{3}$,
$\Lambda^{(2)}(x)=0.7x^{2}\\!+\\!0.3x^{3}$ ($G^{*}=0$)

Figure 2: The PLR obtained from simulation or approximation and the
corresponding AVP vs. $\Phi$ and $G$ for the scenario in Example 1. We
consider two sets of degree distributions, namely,
$\Lambda^{(1)}(x)=\Lambda^{(2)}(x)=0.5x^{2}+0.5x^{3}$ and
$\\{\Lambda^{(1)}(x)=0.7x+0.3x^{3}$, $\Lambda^{(2)}(x)=0.7x^{2}+0.3x^{3}\\}$.

##  V Numerical Results

We next solve the optimization (4) for the scenario in Example 1 with
different target AVPs. We let $d=3$, i.e.,
$\Lambda^{(k)}(x)=\sum_{\ell=0}^{3}\Lambda^{(k)}_{\ell}x^{\ell}$. We keep the
same update probability $\mu_{1}\\!=\\!\mu_{2}\\!=\\!\mu$ for both classes and
control the relative difference between the probabilities of activating users
in different classes via $(\Lambda^{(1)}_{0},\Lambda^{(2)}_{0})$. As
$\Lambda^{(k)}_{3}=1-\sum_{\ell=0}^{2}\Lambda^{(k)}_{\ell}$, the optimization
variables are
$\big{(}\mu,\\{\Lambda^{(k)}_{\ell}\\}_{k\in\\{1,2\\},\ell\in\\{0,1,2\\}}\big{)}$.
The AVP is computed as in Proposition 2 with the approximated/simulated PLR.
We numerically solve (4) by means of the Nelder-Mead simplex algorithm [13], a
commonly-used search method for multidimensional nonlinear optimization.
However, we note that this heuristic method can converge to nonstationary
points and is highly sensitive to the initial values of
$\\{\mu,\Lambda^{(k)}\\}$. We try multiple initializations by sampling the
search space with a step $0.1$, and by running the optimization multiple
times.

TABLE I: Optimized update probability and degree distributions for Example 1 with approximate PLR and different target AVPs $(\epsilon_{1},\epsilon_{2})$ $\epsilon_{1}$ | $\epsilon_{2}$ | $~{}U\mu~{}$ | $~{}\Lambda^{(1)}_{0}~{}$ | $~{}\Lambda^{(1)}_{1}~{}$ | $~{}\Lambda^{(1)}_{2}~{}$ | $~{}\Lambda^{(1)}_{3}~{}$ | $~{}\Lambda^{(2)}_{0}~{}$ | $~{}\Lambda^{(2)}_{1}~{}$ | $~{}\Lambda^{(2)}_{2}~{}$ | $~{}\Lambda^{(2)}_{3}~{}$ | $\Phi$ | $[1\\!-\\!(1\\!-\\!\mu)^{M}](1\\!-\\!\Lambda_{0})$
---|---|---|---|---|---|---|---|---|---|---|---|---
$10^{-5}$ | $10^{-3}$ | $0.71$ | $0.05$ | $0$ | $0.24$ | $0.71$ | $0.05$ | $0$ | $0.25$ | $0.7$ | $1.77$ | $0.017$
$10^{-4}$ | $10^{-3}$ | $0.65$ | $0.19$ | $0$ | $0.11$ | $0.7$ | $0.02$ | $0$ | $0.46$ | $0.52$ | $1.59$ | $0.015$
$10^{-4}$ | $10^{-2}$ | $0.69$ | $0.25$ | $0.19$ | $0.56$ | $0$ | $0.39$ | $0.01$ | $0.6$ | $0$ | $0.87$ | $0.011$
$10^{-3}$ | $10^{-2}$ | $0.56$ | $0.24$ | $0.36$ | $0.4$ | $0$ | $0.21$ | $0.12$ | $0.67$ | $0$ | $0.79$ | $0.011$
$10^{-3}$ | $10^{-1}$ | $0.59$ | $0.13$ | $0.87$ | $0$ | $0$ | $0.52$ | $0.48$ | $0$ | $0$ | $0.33$ | $0.008$

The optimization results with approximate PLR for some different target AVPs
$(\epsilon_{1},\epsilon_{2})$ are shown in Table I. We observe that for the
mild requirement $(\epsilon_{1},\epsilon_{2})=(10^{-3},10^{-1})$, using
degrees $0$ and $1$ is sufficient. As the requirement becomes more stringent,
one needs to use an increasing fraction of degree $2$ and eventually degree
$3$. Also, the users should be activated more frequently, indicated by the
probability $[1-(1-\mu)^{M}](1-\Lambda_{0})$ shown in the last column of Table
I. As compared to the regular distributions in Fig. 1, our optimized irregular
degree distributions reduce $\Phi$ by about $8\%$ and $20\%$ for the
requirements $(\epsilon_{1},\epsilon_{2})=(10^{-5},10^{-3})$ and
$(\epsilon_{1},\epsilon_{2})=(10^{-4},10^{-2})$, respectively.

In Fig. 3, we show the AVP of the optimized distributions using approximate
PLR (14) shown in Table I and compare it with the optimized distributions
using simulated PLR. The latter can further reduce $\Phi$ by no more than
$0.04$ packets/slot. This confirms that our proposed PLR approximation is
sufficiently accurate. Further experiments also show that using degrees higher
than $3$ is not beneficial for the considered parameters.

$0$$1$$2$$3$$10^{-6}$$10^{-5}$$10^{-4}$$10^{-3}$$10^{-2}$$10^{-1}$$10^{0}$Class
$1$Class $2$$\epsilon_{1}$$\epsilon_{2}$$\Phi$ (packets/slot)Age-violation
probability (a) $(\epsilon_{1},\epsilon_{2})=(10^{-5},10^{-3})$

$0$$1$$2$$3$$10^{-6}$$10^{-5}$$10^{-4}$$10^{-3}$$10^{-2}$$10^{-1}$$10^{0}$Class
$1$Class $2$$\epsilon_{1}$$\epsilon_{2}$$\Phi$ (packets/slot)Age-violation
probability (b) $(\epsilon_{1},\epsilon_{2})=(10^{-4},10^{-3})$

$0$$1$$2$$3$$10^{-6}$$10^{-5}$$10^{-4}$$10^{-3}$$10^{-2}$$10^{-1}$$10^{0}$Class
$1$Class $2$$\epsilon_{1}$$\epsilon_{2}$$\Phi$ (packets/slot)Age-violation
probability (c) $(\epsilon_{1},\epsilon_{2})=(10^{-4},10^{-2})$

$0$$1$$2$$3$$10^{-6}$$10^{-5}$$10^{-4}$$10^{-3}$$10^{-2}$$10^{-1}$$10^{0}$Class
$1$Class $2$$\epsilon_{1}$$\epsilon_{2}$$\Phi$ (packets/slot)Age-violation
probability (d) $(\epsilon_{1},\epsilon_{2})=(10^{-3},10^{-2})$

Figure 3: The AVP vs. $\Phi$ for the optimized distributions for the scenario
in Example 1. Red solid lines represent the optimized distributions with
approximate PLR shown in Table I. Blue dashed lines represent the optimized
distributions with simulated PLR.

##  VI Conclusion

We investigated the trade-off between the AVP and power consumption in a
status-update system with multiple classes of users operating according to the
IRSA protocol. Specifically, we illustrate the benefits of jointly optimizing
the update probability and the degree distributions of each class to minimize
the average number of transmitted packets per slot. To perform this
optimization efficiently, we proposed an easy-to-compute PLR approximation,
which yields an accurate approximation of the AVP. Our simulation results
suggest that irregular distributions are needed, and degrees up to $3$ are
sufficient for the considered setting.

## References

  * [1] M. Berioli, G. Cocco, G. Liva, and A. Munari, “Modern random access protocols,” _Foundations and Trends in Networking_ , vol. 10, no. 4, pp. 317–446, Nov. 2016.
  * [2] G. Liva, “Graph-based analysis and optimization of contention resolution diversity slotted ALOHA,” _IEEE Trans. Commun._ , vol. 59, no. 2, pp. 477–487, Feb. 2011.
  * [3] A. Kosta, N. Pappas, and V. Angelakis, “Age of information: A new concept, metric, and tool,” _Foundations and Trends in Networking_ , vol. 12, no. 3, pp. 162–259, Nov. 2017.
  * [4] R. D. Yates and S. K. Kaul, “The age of information: Real-time status updating by multiple sources,” _IEEE Trans. Inf. Theory_ , vol. 65, no. 3, pp. 1807–1827, Sep. 2019.
  * [5] Z. Jiang, B. Krishnamachari, X. Zheng, S. Zhou, and Z. Niu, “Timely status update in wireless uplinks: Analytical solutions with asymptotic optimality,” _IEEE Internet Things J._ , vol. 6, no. 2, pp. 3885–3898, Apr. 2019.
  * [6] Y. Gu, H. Chen, Y. Zhou, Y. Li, and B. Vucetic, “Timely status update in Internet of Things monitoring systems: An age-energy tradeoff,” _IEEE Internet Things J._ , vol. 6, no. 3, pp. 5324–5335, Jun. 2019.
  * [7] A. Munari, “Modern random access: An age of information perspective on irregular repetition slotted ALOHA,” _IEEE Trans. Commun._ , vol. 69, no. 6, pp. 3572–3585, Jun. 2021.
  * [8] M. Ivanov, F. Brannstrom, A. Graell i Amat, and P. Popovski, “Broadcast coded slotted ALOHA: A finite frame length analysis,” _IEEE Trans. Commun._ , vol. 65, no. 2, pp. 651–662, Feb. 2017.
  * [9] A. Graell i Amat and G. Liva, “Finite-length analysis of irregular repetition slotted ALOHA in the waterfall region,” _IEEE Commun. Lett._ , vol. 22, no. 5, pp. 886–889, May 2018.
  * [10] M. Ivanov, F. Brannstrom, A. Graell i Amat, and G. Liva, “Unequal error protection in coded slotted ALOHA,” _IEEE Wireless Commun. Lett._ , vol. 5, no. 5, pp. 536–539, Oct. 2016.
  * [11] A. Amraoui, A. Montanari, T. Richardson, and R. Urbanke, “Finite-length scaling for iteratively decoded LDPC ensembles,” _IEEE Trans. Inf. Theory_ , vol. 55, no. 2, pp. 473–498, Feb. 2009.
  * [12] A. Amraoui, A. Montanari, and R. Urbanke, “Analytic determination of scaling parameters,” in _Proc. IEEE Int. Symp. Inf. Theory (ISIT)_ , Seattle, WA, USA, Jul. 2006, pp. 562–566.
  * [13] J. A. Nelder and R. Mead, “A simplex method for function minimization,” _The computer journal_ , vol. 7, no. 4, pp. 308–313, Jan. 1965.

solutionfile