1 Introduction
Adrien-Marie Legendre was a French mathematician, who in the early 19th century took a particular interest in predicting the trajectories of comets. The understanding of the movements of celestial objects represented one of the main applications of mathematics at that time. There were great scientific challenges in astronomy, such as the stability of the solar system, but there were also more earthbound applications, such as celestial navigation. Legendre wanted to describe the trajectory of a comet by a parabolic path, and he wanted to use observations of points from a comet’s trajectory to infer the best path coefficients. When Legendre substituted the observed points into the path equation he was left with more linear equations in the unknown coefficients than the number of coefficients. There was no solution; the observations did not fall exactly on any parabolic path.
Legendre was not the first to deal with such a path fitting problem, and various techniques had been developed, none of which were completely satisfactory. One solution was to disregard some equations so that a unique solution could be found to the remaining. The resulting path would then match some observed points exactly. But which equations should be disregarded? And would the resulting path be a good approximation to the omitted observation points?
Legendre’s parabolic paths could only produce an approximation of the observations and the future trajectory of the comet. It would not be an exact match. The deviations—or errors—from predicted trajectories of celestial objects were primarily due to unsystematic inaccuracies in the measurements. For Legendre and his contemporaries there was no clear solution as to how deviations from a path should be described and handled. The statistical theory was in its infancy, and the numerical computations required to obtain a descriptive statistical understanding of their nature was very labor intensive. Today it is well known that such deviations can be described statistically, and the methods of path or curve fitting are beneficially studied using statistical theory.
Legendre found an elegant and applicable algebraic solution to his path fitting problem—as described in the section below—which was quickly adapted by astronomers and geodesists in France and Germany. He did, however, not realize the statistical nature of his method. This aspect was developed and expanded by others in the 19th and 20th century, and the theory evolved into the statistical subject we know today as regression.
Regression is a core pillar of statistics, and the statistical theory of regression is a theory about models and methods for relating observables, for example, the coordinates on the path of a comet. The general idea of relating observables through a combination of a systematic and an unsystematic component penetrates most areas of science, technology and business today. The observables are typically treated asymmetrically with one of them pointed out as the outcome variable and the others as predictor variables, and a regression model is a probability model of the outcome given the predictors. It is largely the same type of model that the much younger field of machine learning deals with under the term supervised learning.
In a regression model, the systematic component represents a predictable relation between the outcome and the predictors, while the unsystematic component represents an unpredictable residual variation. This unpredictable variation may be due to measurement errors, sampling inaccuracies, lack of information or other mechanisms, whose contribution cannot be captured by the systematic component.
In regression theory, the unsystematic component is treated statistically. That is, even if the residual variation of individual observations is unpredictable, it can be described and modeled in terms of probability distributions. From this perspective, the theory has to deliver regression models and methods for fitting such models to data that: i) adequately capture the systematic component; ii) adequately describe the unsystematic component; iii) adequately answer questions we might have about the relation between observables. Additionally, the theory investigates properties of such regression models and methods, e.g., how do we compare different methodologies, and how do we choose among models and methods?
1.1 A brief history of regression
We give in this section a historical perspective on regression. Tracking the origin of regression methods and its terminology provides a historical context for the theory developed in this book. It is a complex story, and the introduction will provide a bird’s-eye perspective on what this book is about.
1.1.1 Linear regression
Legendre’s book (Legendre 1805) on the determination of the orbit of comets was published in 1805 containing the appendix Sur la Méthode des moindres carrés. It outlined a method for fitting a linear equation with unknown coefficients to observations. The English translation of the title of the appendix is On the Method of least squares, and it is the earliest known appearance of the method in print. Independently of Legendre, Carl F. Gauss (Gauss 1809) discovered the Gaussian distribution and used it to derive the method of least squares from Laplace’s principle of inverse probability. Gauss was like Legendre motivated by problems in astronomy—specifically the description and prediction of the motions of planets. The invention of the method of least squares marks a turning point in the history of regression. Earlier work by Pierre-Simon Laplace, and related work by Roger J. Boscovich with applications to geodesy, using the method of least absolute deviation, had been less successful. The method of least absolute deviation leads to a nonlinear estimator that is hard to compute. By contrast, the least squares estimator is a linear estimator that solves the (linear) normal equation, and Gauss presented a useful method for its numerical computation. The generality of the linear model framework considered by Gauss is impressive and corresponds essentially to that treated in Chapter 2 in this book.
In the years after the invention of the method of least squares, Gauss and Laplace derived several theoretical results about the method including sampling properties of the least squares estimator and confidence intervals for the parameters. The logical foundation of the method of least squares was, however, confusing at the time. Gauss’s initial derivation was—as mentioned above— based on Laplace’s principle of inverse probability. Later, Laplace as well as Gauss sought justification of the method in the sampling properties of the resulting estimator. Specifically by investigating if the estimator has a form of minimal estimation error. Gauss showed what is today known as the Gauss-Markov theorem: the least squares estimator has minimal variance among all linear unbiased estimators.
The models that Gauss fitted to data using his method of least squares were not known to him as regression models. That terminology originates from Sir Francis Galton (Galton 1886). In studies of hereditary traits Galton observed an apparent regression (in the sense of a movement) towards a mediocre state in the offspring compared to the parent(s). Originally he considered the size of plant seeds and later the height of human individuals. What Galton had discovered was the phenomenon known today as regression toward the mean. He explained the phenomenon statistically. The trait is partly inherited from the parent and partly a consequence of additional variation. This results in a bivariate distribution showing a co-relation of the trait in parent and offspring. Though Galton correctly described the phenomenon statistically, he wrongly explained the regression phenomenon as being caused by inheritance from the ancestry of the parent in addition to the parent themselves.
To understand the co-relation Galton observed in his hereditary studies, he contributed to the theory of the bivariate normal distribution and introduced the general notion of correlation. The technical details of what is known today as the Pearson correlation coefficient were worked out by Karl Pearson, though. Galton, nevertheless, realized the usefulness of regression and correlation in various application areas, and the regression terminology was adapted by others, notably Karl Pearson. Interestingly, the initial developments of regression analysis in British mathematical statistics were independent of the continental literature, and the regression models were, in particular, not framed in the linear model terminology of Gauss.
In the first half of the 20th century, Sir Ronald A. Fisher had an enormous impact on mathematical statistics. He developed Analysis of Variance (ANOVA), introduced the likelihood concept and advocated the general principle of maximum likelihood estimation among other things. The primary application of ANOVA was to designed experiments (Fisher 1935), and though far from Gauss’s applications on the modeling of the motion of celestial objects, the ANOVA models fit into the framework of the linear model.
Fisher did not rely on or acknowledge the contributions of his predecessors to any great extent, and he rediscovered known results about the method of least squares. But he also contributed significantly to the theory, for instance by deriving sampling distributions of test statistics, notably the \(F\)-test statistic. His presentation and applications of ANOVA models have been highly influential on the current terminology, which is, for instance, reflected in the anova() function in R for testing linear hypotheses within the linear model. Subsequent developments with contributions from many people have turned the linear model into a unified and general statistical theory as outlined in one of the first unifying textbooks on the subject by Shayle R. Searle (Searle 1971).
The history of regression is intertwined with the general history of statistics, and for the history before 1935 the work by Anders Hald (Hald 2007) is a valuable source. For the more recent, and much richer, development of nonlinear regression we give a brief outline in the following sections that focus on a few selected original sources.
1.1.2 Beyond linear regression
Regression models, methods and applications proliferated in the second half of the 20th century. Computers played a pivotal role, first of all by easing the solution of nonlinear estimating equations by iterative methods, which paved the way for applications of nonlinear regression models and nonlinear estimation methods. The use of computer simulations and Markov Chain Monte Carlo (MCMC), in particular, also made Bayesian regression modeling a practical possibility. In this brief account of the history, it is impossible to do justice to all the different directions, and we will therefore narrow the focus to regression models treated in this book.
Even if linear regression models are flexible, some regression problems do not fit into the linear modeling framework. One of the early applications of a nonlinear regression model was to the modeling of mortality of an organism as a function of the dosage of a toxic agent. Chester I. Bliss (Bliss 1935) proposed a transformation method where dosage is log-transformed and death frequency is transformed by the inverse distribution function for the normal distribution with mean 5 and unit variance. The resulting transformed death frequencies were called probits. The model posited that the transformed data should fall on a straight line, and Bliss suggested using the method of weighted least squares for fitting such a line.
Already one year earlier, Bliss (1934) presented a similar idea with the probits computed from an S-shaped curve inspired by, but not identical to, the distribution function for the normal distribution. His 1935-paper had, however, an important appendix written by Fisher, which showed how to appropriately deal with dosages for which none or all of the organisms survived. In this appendix, Fisher outlined how the method of maximum-likelihood could be used for estimation. In doing so he came up with a correction step in the weighted least squares fit, which is what we today know as the Fisher scoring algorithm.
The probit model introduced by Bliss is one among a number of nonlinear regression models—notably models of counts—that share a similar structure. These regression models were unified by John Nelder and Robert Wedderburn (Nelder and Wedderburn 1972) in the framework they called generalized linear models, which includes the linear model as a special case but extends it in two notable ways: by including a nonlinear link function between the predictors and the outcome, and by allowing for a more flexible description of the outcome distribution.
For generalized linear models, the method of least squares is replaced by maximum likelihood estimation. Nelder and Wedderburn showed that the maximum likelihood estimator can be computed iteratively by the Fisher scoring algorithm, where each step in this algorithm is, in fact, the solution of a weighted least squares problem, and they called it the Iterative Weighted Least Squares (IWLS) algorithm. They also paraphrased some terminology and methodology from linear models, in particular the use of analysis of deviance as a generalization of ANOVA.
As a curious detail, Nelder and Wedderburn outlined in their original paper how generalized linear models could form a useful basis for courses in statistics, which would unify the treatment of a number of statistical models that had previously been treated separately. The theory of generalized linear models was further developed in the 70s and 80s, and an authoritative reference on this classical theory is the book by McCullagh and Nelder (1989). In addition to his theoretical contributions, Nelder initially chaired the GLIM Working Party, who developed the statistical software program GLIM (Generalized Linear Interactive Modeling). This program contributed to making generalized linear models applicable in practice, and it was a source of inspiration for the later implementations (Chambers and Hastie 1991) in S and R.
Another type of regression problem that does not fit directly into the linear modeling framework is the modeling of survival times. As for generalized linear models there is a question of modeling scale; at which scale is survival time most appropriately associated to other variables? It can be handled similarly as for generalized linear models via a link function, but there is an additional practical problem with observations being right censored.
Survival models and estimation of life tables have a long history, in particular in actuarial science and medical statistics, but a systematic approach to survival regression models was not made until the seminal contributions by Sir David R. Cox (Cox 1972). He introduced the semiparametric proportional hazards model, which has since become the most widely used survival regression model. Cox made two major contributions in his paper. He introduced the partial likelihood, which provided a means for deriving a sensible estimating equation, and he presented approximate sampling distributions of the resulting estimator.
The paper by Cox inspired a rapid development in the theory of survival regression models, which clarified some of his heuristic arguments and extended his model class considerably. Central to this development was the use of methods from the theory of stochastic processes, in particular continuous time martingale theory. Modern survival analysis theory is heavily based on the use of stochastic processes, and an authoritative account of this theory is the by now classical book Statistical Models Based on Counting Processes (Andersen et al. 1993).
On the practical side, the estimating equation derived from Cox’s partial likelihood in the semiparametric framework is nonlinear—and so are parametric likelihood based estimating equations. The estimator thus has to be computed by iterative methods. As a curious historical detail, practical computations could be carried out using the framework of generalized linear models, and the GLIM program could be used to fit the proportional hazards models in Cox’s semiparametric framework as well as for several parametric models. This made these models and methods applicable early on, but specialized software implementations were later developed, such as the survival package for Splus by Terry Therneau (Therneau and Grambsch 2000). The package was ported to R, and it is one of the classical core packages for survival analysis using R.
1.1.3 Regression and machine learning
There were notable developments outside of mainstream statistics related to regression models, but little interaction between the different scientific communities. Frank Rosenblatt (Rosenblatt 1962) explored perceptrons and artificial neural networks for binary classification from multiple input variables. Support vector machines (Cortes and Vapnik 1995) were developed for similar purposes, with the conceptual framework established in 1965 by Vladimir Vapnik.
Artificial neural networks and support vector machines are today core subjects in the field of machine learning. Both can be understood as regression models of a binary outcome variable, and both can be altered to model a continuous outcome variable as well. The models are typically presented as classification models rather than probabilistic models, and machine learning theory and practice focus on classification performance over interpretable relations between the predictors and the outcome. Statistics is historically connected to other scientific fields, where explanations, interpretations and answers to inferential questions are central. That is, questions about the nature and magnitudes of the relations between any specific predictor and the outcome. Inferential questions have therefore driven the theoretical and methodological development in statistics, while classification performance has driven the development in machine learning.
In 2001 Leo Breiman (Breiman 2001) divided the data modeling communities into “two cultures”: those using stochastic (or probabilistic) data models and those using algorithmic models. Breiman ascribed the developments of algorithmic modeling largely to the machine learning community, and he strongly criticized the conventional statistical community for their (simplistic) probabilistic data modeling. Undoubtedly, contemporary regression modeling has benefited greatly from developments in machine learning, and the “two cultures” are linked much closer today than ever before.
The development of lasso (Tibshirani 1996) is an interesting early example. It was first developed for the linear regression model, but the lasso estimator is nonlinear and biased—as opposed to the least squares estimator. It required specialized algorithms for its computation, but it turned out to have certain computational and statistical benefits that made it possible to tackle large scale data analysis in new and useful ways.
Not long after the introduction of lasso, in 2001, the first edition of The Elements of Statistical Learning (Hastie et al. 2009) was published. It is hardly a coincidence that this was the same year as Breiman’s paper was published. The book had a notable impact on the statistics community by bridging the gap between the machine learning literature and the statistics literature. Along with lasso, the book introduced a range of nonlinear regression models and optimization methods to a broad statistical audience. This was a timely publication and in combination with new applications to large scale data modeling it spurred a substantial development of new regression modeling methods and exciting challenges.
1.1.4 Computation
In the long history of regression, one particularly influential development on how we work with regression models today is the growth of computational resources.
Data analysis with regression models has been completely transformed by the increase of computer power and the developments of statistical programming environments such as R. The computer serves all data analysts, and interactive and exploratory data analysis in the spirit of Tukey (1977) is now routine. Nonparametric and resampling based statistical methods, that require substantial computations, have also made data analysis possible in situations where we lack theoretical methods or where we want to impose minimal modeling assumptions.
The computer has, in addition, become the lab bench of statistics, and it has brought an experimental component to statistics. Most methodologies and theoretical results are today scrutinized via simulation studies, which can inform us about the applicability and limitations of the methods and the theory.
When the history of regression is viewed from our vantage point, with the computational resources we have available today, some of the theory and methodologies may appear outdated and irrelevant. And some surely are! Nobody should use least squares methods because solving linear equation systems is possible by hand, and nobody should fit generalized linear models using the arcane syntax and model specification of GLIM. However, important residues of history are alive and central to regression today because they represent solid and well understood methodology, and these residues function as the foundation we build on.
1.1.5 The four components to regression modeling
The history of regression contains four recurrent components.
- Models
- Estimation methodology
- Sampling properties
- Model skepticism and critique
The first component—the model development—was, and is, driven by applications, and the second component was developed to compute estimates of model parameters from data, that is, to fit models to data. Over the course of history, a supply of models, such as the linear and generalized linear models, have been developed together with methods or principles for parameter estimation, such as the methods of least squares or maximum likelihood.
The third component—the sampling properties—was developed to understand theoretical properties of the various estimation methodologies and to provide (frequentistic) quantification of uncertainty.
The developments of components one through three are all good and well, but from Gauss to Cox there has been grave concerns about whether the model assumptions that supported the theory could be justified for particular applications. Techniques of more or less formal character have been developed and used throughout history for validating if model assumptions are plausible. Though such techniques will never provide certainty, they serve as important sanity checks that keep models and methodology aligned with data.
The four components: models; estimation methodology; sampling properties; and model critique, are recurring components in this book as they are in history. The intention of this book is to teach the reader these four cornerstones of statistics in the framework of regression modeling.
1.2 What this book is—and what it is not
This is a textbook for a graduate-level course in regression modeling and regression theory. The purpose of the book is to make the reader develop their own independent understanding of what regression models are, what they are capable of, and how we can fit them to data and work with them to answer scientific questions. There exists today a very broad range of such models and related methodology, and this book had to be selective.
This book is therefore not
- a survey of regression models and methods
- a cookbook for regression modeling in practice
- an encyclopedia for answering all questions related to regression modeling
- a guide to state-of-the-art regression modeling
In fact, most models treated in this book are somewhat classical and not state of the art, though they are widely used in practice. It is a deliberate choice to use relatively simple and classical models as illustrative examples of general theory and methodology—even if new applications may need to look beyond this book for state-of-the-art solutions.
This book is not a normative prescription either of how a data analyst should analyze data using regression modeling techniques. The book presents a development of concrete model classes and practically usable methodologies, and a description of properties and justifications of various regression analysis techniques. The author is concerned with accurately developing the theory and methodologies in sufficient detail for the reader to completely absorb them. The author is certainly also concerned with developing good habits around practical data analysis but less convinced about having the definitive say on what the reader should do.
To be concrete, the book primarily covers the following three regression model classes:
- linear models
- generalized linear models
- parametric and proportional hazards regression models
These are classical subjects and an obligatory part of an education in statistics. At the same time they serve as good examples when developing general statistical principles and methodologies such as estimation methodology, likelihood based methods as well as methods for interval estimation and statistical tests.
The practical value of the three classical model classes is boosted substantially when combined with two techniques the book also covers:
- basis (or feature) expansions
- regularized estimation
The three primary questions we seek to answer with regression models are:
- quantifications of conditional associations
- tests of conditional independencies
- predictions of outcomes given the value of the predictors
The first question deals with quantification of the associational strength between the outcome and one predictor given the other predictors. The second question deals with formally testing the hypothesis that there is no association between the outcome and one predictor given the other predictors. These two questions are inferential and ask about specific properties of an unknown conditional distribution. The last question is directly asking for a (good) model of the conditional distribution of the outcome given the predictors.
The book deals with the questions above exclusively in terms of associational regression modeling, where the regression model relates the distribution of the outcome variable to one or more predictor variables in terms of a full or partial specification of the conditional distribution of the outcome given the predictors. Importantly, associational regression models do not necessarily capture causal relations, that is, relations that persist under interventions on or manipulations of the observables.
Associational regression models do play a role when inferring causal relations, but we refer to the substantial literature on causal inference for a discussion of the untestable structural assumptions that we need to make to be able to draw causal inference. A thorough understanding of the subject matter field, or knowledge from other studies, might allow us to make such necessary structural assumptions, but they are typically untestable using only the dataset we have at hand for the regression analysis.
Associational models are, however, directly applicable for predicting outcomes given observations of the predictors, which is used for medical diagnosis and prognosis, in businesses to predict customer behavior, to predict risk in insurance companies, pension funds and banks, to make weather forecasts and in many other areas where it is of interest to know what we cannot (yet) observe.
It is fairly clear what makes a model good for prediction. Given a method for quantification of predictive accuracy, the best model is the most accurate model. The accuracy of a prediction is typically quantified by a loss function, which actually quantifies how inaccurate a prediction is. The smaller the loss is the more accurate is the prediction. The expected loss quantifies how accurate the predictive model is on average. Predictive performance is something we can measure, and the book develops general methodology for quantification of a model’s predictive strength.
It is slightly more tricky to answer the inferential questions, partly because methods for doing so generally rely on some level of correctness of the regression model, that is, does the model fit the data? The two final points, that we emphasize in this introduction as part of the book, are:
- model diagnostics
- model-robust inference
Diagnostics are post-model fitting procedures for discovering ways that data violate model assumptions. If the model does not fit we cannot trust inferential conclusions based on model assumptions. The two ways around such problems are: either alter the model so that it fits the data; or rely on inferential procedures that are robust to model misspecifications. Both perspectives will be developed in the book.
1.3 Organization
The book consists of a combination of theory and larger case studies. The theory chapters focus on presenting the mathematical framework for regression modeling together with the relevant mathematical theory such as derivations of likelihood functions, estimators, estimation algorithms and properties of estimators. Other chapters focus on more general methodological questions. This includes developing methods for dealing with the complex decision process that practical data analysis is.
While smaller data examples are presented in the theory chapters, the larger case studies are developed to better illustrate how theory supports applications. The hope is that this will ease the transition from theory to practice. The price to pay is that there are some distractions due to real world data challenges. Data rarely behaves well. There are missing values and outliers, the models do not fit the data perfectly, the data comes with a strange encoding of variables and many other issues. Issues that require decisions to be made and issues on which many textbooks on statistical theory are silent.
By working through the case studies in detail it is the hope that many relevant practical problems are illustrated and appropriate solutions are given in such a way that the reader is better prepared to turn the theory into applications on their own.
1.4 Using R
We use the programming language R1 throughout to illustrate how good modeling strategies can be carried out in practice on real data, but the book will not provide an introduction to the R language. Consult the R manuals (R Core Team 2025+) or the many textbooks on R programming and data analysis (Wickham 2019; Wickham et al. 2023). The classical and influential book Modern Applied Statistics with S (Venables and Ripley 2002) and the corresponding MASS package (comes with the R distribution) should also be mentioned. Several classical statistical models and methods are supported by the MASS package and documented in that book.
The case studies in this book are complete with R code that covers all aspects of the analysis. They represent an integration of data analysis with documentation and reporting. This is an adaptation to data analysis of literate programming (Knuth 1984). The main idea is that the writing of a report that documents the data analysis and the actual data analysis are merged into one document. This supports the creation of reproducible analysis, which is a prerequisite for reproducible research. To achieve this integration the book was written using the Quarto2 publishing system and the R package knitr3. Quarto and knitr allows you to integrate text with chunks of R code and produce an output document in various formats, e.g., a webpage or a pdf document, that include text, code and output from evaluating the code. The use of the RStudio4 or Positron5 integrated development environments (IDE) is also recommended. Both are developed by Posit PBC as open source projects and the desktop versions are free to use.
The book source is available from the RwR GitHub repository, which also contains information on how to recreate the R environment used for writing the book and how to obtain the RwR package containing the datasets used.