Biostatistics 16, 1–9
DOI: 000
October 2026
Neural Networks as Linear Regression: An Introduction for Statisticians
Abigail Loe1,∗ , Susan Murray2 , and Zhenke Wu2 1
Department of Mathematics, Statistics and Computer Science, Macalester College,
arXiv:2606.23601v1 [stat.ML] 22 Jun 2026
Saint Paul, Minnesota, U.S.A. 2
Department of Biostatistics, University of Michigan, Ann Arbor, Michigan, U.S.A. *email: [email protected]
Summary:
Neural networks are a commonly used prediction tool in computer science and statistics. However, the
barrier to entry of this interesting field remains high, particularly for classical statisticians trained in a frequentist perspective. In this letter, we demystify neural networks by describing networks that approximate a linear regression and describe common customizations that provide a foundation for further study. Key words:
Approachable neural networks; Linear regression.
This paper has been submitted for consideration for publication in The American Statistician
Neural Networks as Linear Regression: An Introduction for Statisticians
1
1. Statistical Review of Neural Networks Neural networks are a popular way to model potentially complex functions with many variations and adaptations (see, for example Zazo et al., 2016; Radha et al., 2019; Peng et al., 2023). Although many features of neural networks have a clear statistical basis, jargon for describing these features are often quite different for computer science and statistical researchers, causing much unnecessary confusion. To familiarize statistical audiences with the neural network algorithm used in our research, we first introduce a neural network algorithm that corresponds to linear regression in Section 1.1. Then in Section ?? we describe a neural network algorithm appropriate for use with the longitudinal data. Additional practical considerations of applying this neural network algorithm in our setting are given in Section 1.2.
1.1 Linear Regression as a Neural Network In the language of computer scientists, neural networks are a user-specified architecture (model) with hidden layers of nodes, consisting of weights on the edges and biases (similar to intercept; different from the bias typically used in statistics) applied to model inputs for the purpose of predicting an outcome. Consider the regression model g(E[Y |Z 1 , . . . , Z p ]) = Z⊤ β, where Y , Z j , j = 1 . . . p are n-dimensional vectors, β is of dimension p × 1, and Z is the (p+1)×n-dimensional design matrix, which includes the intercept term. Here and throughout this manuscript, we use bold face type to denote matrices or vectors. The g(·) function, or link function in statistical language, is typically the identity function for linear regression settings. Generally, we will use the term g −1 (η) to denote the term that is used to minimize squared-error loss, i.e. {Y − g −1 (η)}⊤ {Y − g −1(η)}, where η = Z⊤ β and η is of dimension n × 1. Translating neural network to statistical terminology: (1) model inputs are covariates, Z j , j = 1, . . . , p; (2) weights are linear regression slope parameters βj , j = 1, . . . , p; (3)
2
Biostatistics, October 2026
bias corresponds to the linear regression model intercept, β0 ; (4) a node is the functional g −1(Z⊤ β), which reduces to Z⊤ β in the setting mirroring linear regression with an identity link. Individual nodes make up the hidden layers of the neural network, termed hidden b just the vector of outcome because neural network algorithms do not report estimates, β, b Computer scientists refer to g −1 (·) as an activation function predictions, Yb = g −1 (Z⊤ β). that estimates the desired outcome, often in several different layers within the algorithm. Hereafter, we use the terms “slope” to refer to a weight, “intercept” to refer to bias, and “parameter estimates” to refer collectively to weights and biases used in the neural net contexts. Figure 1 displays three increasingly more sophisticated neural networks in Panels (A), (B) and (C), where we include an additional subscript on η to reflect changes across panels. Panel (A) displays a linear regression model in a manner familiar to neural network researchers with a single hidden layer composed of one node; the corresponding statistical optimization function is displayed below the network. Covariates are represented by circles on the left of the diagram, while the linear predictor is displayed in the rectangular box. Researchers have found that an increase in the number of nodes in the hidden layer may improve prediction, particularly in the case where the relationship between covariates and outcome is more complex than a linear model. Figure 1, Panel (B) depicts a neural network with a single hidden layer composed of three nodes stacked on top of each other. The term, ηB , from Panel (B) is no longer a simple linear predictor as in panel (A), and involves parameter vectors α1 , . . . , α3 in addition to potentially different parameters, β1 , . . . , β3 applied to the same set of predictors, Z. The form of gℓ−1 , ℓ = 1, 2, 3 may also vary across nodes. A network with two hidden layers is depicted in Figure 1, Panel (C); hidden layers 1 and 2 have three and two nodes, respectively. The term, ηC , in Panel (C) has increased in complexity from Panel (B). Within ηC , readers can identify terms resembling ηB fed into two
Neural Networks as Linear Regression: An Introduction for Statisticians
3
−1 different activation functions h−1 1 (·) and h2 (·). Two issues emerge at this stage of neural
network complexity that require comment. First, the number of parameters involved in minimizing squared error loss based on the η in panel C has more than doubled. While the networks in Panels (B) and (C) have been found to potentially improve prediction compared to the network in Panel (A), parameter estimates are no longer guaranteed to be uniquely identifiable; it may be the case that there are multiple local minima of the loss function that are not the global minimum (Goodfellow et al., 2016), as well as multiple different matrices of parameter estimates that may provide the same predictions. This is an open area of research, with mathematicians, statisticians and computer scientists attempting to find conditions that guarantee network architecture and parameter estimate uniqueness (see, for example, Sussmann, 1992; Fefferman et al., 1994; Vlačić and Bölcskei, 2022). In practice, however, most theorists recognize that it is generally sufficient if a local minimum is found; i.e. even local minima have sufficiently good predictive performance (Goodfellow et al., 2016). The second issue that arises in Panel C is the potential for scaling differences to emerge between predicted outcomes generated within the hidden layers as the result of using so many different activation functions, gℓ−1 and h−1 ℓ2 , ℓ1 = 1, 2, 3 and ℓ2 = 1, 2, within the neural 1 network architecture. Computer scientists have approached this issue pragmatically by using the tanh function to rescale elements within the neural net architecture to be between a common -1 to 1 scale to reduce problematic differences. Some practitioners refer to the tanh operation as “normalizing,” though we should note that this is not normalizing in a statistical sense. Under certain assumptions, including non-colinearity of covariates, a linear regression squared error loss function has a closed form solution for parameter estimates. This is not typically the case for neural networks. Instead, iterative algorithms moving in smaller steps
4
Biostatistics, October 2026
towards loss function (local) minima, via a hyperparameter known as the learning rate; each iteration is referred to as an epoch in neural network literature. Popular algorithms are stochastic gradient descent, batch gradient descent, “Adam” (Kingma and Lei, 2014), and adaptive gradient descent (Duchi et al., 2011), each with mini-batch variations that fall under the framework proposed by Robbins and Monro (1951). For a comprehensive review, see Bottou et al. (2018). Because of the large number of parameters built into neural network architectures, these algorithms are particularly susceptible to overfitting. A popular method used by both statisticians and computer scientists to prevent predictions being overfit to the data is to use different data cohorts for model training, validating, and testing. These terms have been used differently depending on the background of an individual researcher, so a quick summary of terms and techniques follows. First, data are split into model building and model testing cohorts. The model building cohort determines the final architecture of the neural network using a k-fold cross-validation algorithm to evaluate architecture hyperparameters including, for instance, the number of hidden layers and the number of iterations used in the optimization algorithm. For a particular architecture defined using a set of hyperparameters, the k-fold cross-validation algorithm (1) splits the model building cohort into training and validation cohorts along a (k − 1) : 1 ratio, (2) minimizes the loss function (squared-error loss for linear regression) in the training cohort to establish a working neural network, (3) uses the working neural network applied to the validation cohort to obtain predictions, Yb , and estimates of the loss function. Steps 1-3 are repeated for each of the k data splits with loss function estimates across the k validation cohorts averaged. This process is repeated for different sets of hyperparameters with a final architecture selection made based on smallest average estimated loss function across the validation cohorts. Once the architecture hyperparameters are determined via k-fold cross-validation, the full
Neural Networks as Linear Regression: An Introduction for Statisticians
model building cohort is used to rebuild the neural network (via re-minimizing the loss function in this cohort). The testing cohort is then applied to the final neural network to determine overall performance. For a detailed description of hyperparameter tuning and neural networks, see Goodfellow et al. (2016), Chapter 11, where decisions on optimization algorithm learning rate and number of epochs are discussed at length. We will later summarize how these three groups may be more flexibly defined in the context of longitudinal outcomes in Section 1.2. To illustrate a typical process, we now describe the process of fitting the linear regression neural network depicted in Figure 1, Panel (A) on training data, though the process generalizes to the architectures displayed in Panels (B) and (C). For each epoch, a = 1, 2, . . . na (where na is often between 10 and 200):
b (a−1) for the linear regression loss function Y b (a−1) = Z ⊤ βe(a−1) , (1) Calculate fitted values of Y where (a) tracks the epoch. For the first epoch, traditionally βe(0) = 0, though some authors note that sufficient convergence to a local minima can depend on the initial starting values (Wang et al., 1993). (2) Calculate the updated loss function values. For linear regression, this would be squared error b = (Y − Y b (a−1) )⊤ (Y − Y b (a−1) ) = {Y − g −1(Z ⊤ βe(a−1) }⊤ {Y − loss with the formula L(Y, Y) g −1(Z ⊤ βe(a−1) }, where g −1 is typically taken to be the identity function. More generally, loss functions mimic those seen in generalized linear models, reflecting the assumed distribution of Y. (3) Take partial derivatives of the loss function with respect to each parameter β0 , . . . , βp to determine the gradient of the loss function. Update βe(a−1) to βe(a) in a step size corresponding to the learning rate and optimization algorithm of your choice. For example, with stochastic gradient descent, βe(a−1) = βe(a) − ▽ ∂ βe∂(a) L, where ▽ is the learning rate, and ∂ βe∂(a) L is the partial derivatives of the loss function with respect to each parameter.
5
6
Biostatistics, October 2026
(4) Repeat steps 1-3 until na is reached.
1.2 Useful Extensions to Bare-bones Networks Neural networks have many proposed extensions. In this section, we selectively describe (1) common parameters available from the off-the-shelf Python module torch.nn, and (2) a few additional neural network design choices in common use for tailoring analyses. Options in torch.nn include: η (computer scientists call this hidden state), parameter pruning (computer scientists call this the drop-out proportion), number of network layers, and the encoding of similarities within individuals via embeddings. A hyperparameter favored by computer scientists is allowing ηt to be extended from a scalar (1)
(2)
(h)
to a matrix in Rn(t)×h , with columns ηt , ηt , . . . ηt , and similarly Z i (t) to be extended from a vector in p × 1 to a matrix of dimension Rp×n(t) . Computer scientists call h the hidden dimension. For dimensions h > 1, the inverse link is applied to a linear combination of the (1)
(h)
h columns of ηt , i.e. g −1(γ0 + γ1 ηt + . . . ηt ). Beyond increasing h, computer scientists have also found that vertically stacking additional architecture at time t increases prediction performance (Graves, 2014). General recommendations are for at most 4 layers of architecture to be “stacked” on top of each other (as used in Sutskever et al., 2014). For instance, the entire architecture at time t seen in Figure 1C might be stacked on top of a similar architecture, where the ηC elements taken from the right side of the top stacked architecture play the role of the covariate inputs (Zp terms) on the left of the bottom stacked architecture. Additional subscripts on parameters and link functions in the bottom stacked architecture would complete this version of stacking architecture. The final two off-the-shelf options include (1) regularization via drop-out and (2) embeddings. Regularization via drop-out refers to randomly forcing a proportion of parameters, pd , to 0 during each epoch of training (Srivastava et al., 2014; Zaremba et al., 4 09). Subject (and potentially time point) embeddings are a concept borrowed from text mining
Neural Networks as Linear Regression: An Introduction for Statisticians
7
(Sutskever et al., 2014; Mikolov and Zweig, 2012; Mikolov et al., 2013). For each individual, a dimension b vector of embedding values may be allocated to quantify similarity of outcomes experienced by individuals. These embedding values are treated as parameters to be estimated in addition to other neural network outputs, and are available after training is complete. Training, validation and testing cohorts described in Section 1.1 may be used in a gridsearch to optimize the number of training iterations selected from the set, {1, . . . , na }, and the learning rate, ▽ selected from a size m grid of candidate learning rates, {▽1 , ▽2 . . . , ▽m },. All hyperparameters may be selected via a grid search algorithm, though a large number of points in a grid increases significantly increases computation time. As such, we chose to focus on learning rate and number of training epochs, two hyperparameters generally recognized to have an outsized effect on network reliability (Goodfellow et al., 2016) Loss function decreases for each additional training iteration are expected in the training set, but not necessarily in the validation set. When increasing the number of training iterations from 1 to na in the validation set for a particular learning rate, ▽ℓ , ℓ = 1, . . . , m if maxp consecutive increases in the loss function are seen, nopt is set to the number of training iterations, prior to the start of when these increases were seen (called early stopping in computer science literature); otherwise nopt = na . The value of maxp is a user chosen value traditionally between 5 and 10 that is called the tolerance or patience in computer science literature. Further details on early stopping may be found in Prechelt (1998). Validation loss and recommended number of epochs are recorded for each candidate learning rate, with final learning rate and recommended epochs, ▽opt and nopt respectively, selected based on the pairing that yields minimal loss in the validation set.
[Figure 1 about here.]
8
Biostatistics, October 2026
Acknowledgments
Conflict of Interest: None declared.
References
Bottou, L., Curtis, F. E., and Nocedal, J. (2018). Optimization methods for large-scale machine learning. Duchi, J., Hazan, E., and Singer, Y. (2011). Adaptive subgradient methods for online learning and stochastic optimization. Journal of Machine Learning Research 12,. Fefferman, C. et al. (1994). Reconstructing a neural net from its output. Revista Matemática Iberoamericana 10, 507–556. Goodfellow, I., Bengio, Y., Courville, A., and Bengio, Y. (2016). Deep learning, volume 1. MIT press Cambridge. Graves, A. (2014). Generating sequences with recurrent neural networks. Kingma, D. P. and Lei, J. L. (2014). Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980 . Mikolov, T., Chen, K., Corrado, G., and Dean, J. (2013). Efficient estimation of word representations in vector space. Mikolov, T. and Zweig, G. (2012). Context dependent recurrent neural network language model. In 2012 IEEE Spoken Language Technology Workshop (SLT), pages 234–239. Peng, Y., Chen, T., Xiao, F., and Zhang, S. (2023). Remaining useful lifetime prediction methods of proton exchange membrane fuel cell based on convolutional neural networklong short-term memory and convolutional neural network-bidirectional long short-term memory. Fuel Cells 23, 75–87. Prechelt, L. (1998). Automatic early stopping using cross validation: quantifying the criteria. Neural networks 11, 761–767.
Neural Networks as Linear Regression: An Introduction for Statisticians
9
Radha, M., Fonseca, P., Moreau, A., Ross, M., Cerny, A., Anderer, P., Long, X., and Aarts, R. M. (2019). Sleep stage classification from heart-rate variability using long short-term memory neural networks. Scientific reports 9, 14149. Robbins, H. and Monro, S. (1951). A stochastic approximation method. The Annals of Mathematical Statistics pages 400–407. Srivastava, N., Hinton, G., Krizhevsky, A., Sutskever, I., and Salakhutdinov, R. (2014). Dropout: a simple way to prevent neural networks from overfitting. The Journal of Machine Learning Research 15, 1929–1958. Sussmann, H. J. (1992). Uniqueness of the weights for minimal feedforward nets with a given input-output map. Neural Networks 5, 589–593. Sutskever, I., Vinyals, O., and Le, Q. V. (2014). Sequence to sequence learning with neural networks. Vlačić, V. and Bölcskei, H. (2022). Neural network identifiability for a family of sigmoidal nonlinearities. Constructive Approximation 55, 173–224. Wang, C., Venkatesh, S., and Judd, J. (1993). Optimal stopping and effective machine complexity in learning. Advances in Neural Information Processing Systems 6,. Zaremba, W., Sutskever, I., and Vinyals, O. (2014-09). Recurrent neural network regularization. Arxiv . Zazo, R., Lozano-Diez, A., Gonzalez-Dominguez, J., T. Toledano, D., and GonzalezRodriguez, J. (2016). Language identification in short utterances using long short-term memory (lstm) recurrent neural networks. PloS one 11, e0146917.
10
Biostatistics, October 2026
(A)
(B)
Hidden Layer
. . .
. . .
Hidden Layer 1
(C)
Hidden Layer 2
. . .
Figure 1: Three fully connected neural network architectures of increasing complexity, with squared error loss. Vectors of covariates are depicted in circles on the left, while nodes containing estimable parameters are displayed in rectangular boxes. Networks are called “fully connected” because each covariate is passed to each node in the first layer, and for networks with more than one layer, outputs from previous layers are passed to all subsequent layers. Panel (A) displays an architecture with one hidden layer containing one node, akin to linear regression in the special case where g is the identity link. Panel (B) displays a network architecture with three nodes in one hidden layer, that are ultimately combined to obtain the linear component, ηB . As Panel (C) contains multiple hidden layers, Panel (C) proposes a network that may also be termed a “deep learner.” Alt. text: The figure containing three different neural networks. Panel (A) contains the simplest, while Panels (B) and (C) increase in complexity. Under each network architecture, the corresponding loss function is displayed.