# An Automated Profile-Likelihood-Based Algorithm for Fast Computation of the Maximum Likelihood Estimate in a Statistical Model for Crash Data

Authors: Issa Cherif Geraldo
Journal: Journal of Applied Mathematics
 (2022)

---

## Abstract

Numerical computation of maximum likelihood estimates (MLE) is one of the most common problems encountered in applied statistics. Even if there exist many algorithms considered as performing, they can suffer in some cases for one or many of the following criteria: global convergence (capacity of an algorithm to converge to the true unknown solution from all starting guesses), numerical stability (ascent property), implementation feasibility (for example, algorithms requiring matrix inversion cannot be implemented when the involved matrices are not invertible), low computation time, low computational complexity, and capacity to handle high dimensional problems. The reality is that, in practice, no algorithm is perfect, and for each problem, it is necessary to find the most performing of all existing algorithms or even develop new ones. In this paper, we consider the computing of the maximum likelihood estimate of the vector parameter of a statistical model of crash frequencies. We split the parameter vector, and we develop a new estimation algorithm using the profile likelihood principle. We provide an automatic starting guess for which convergence and numerical stability are guaranteed. We study the performance of our new algorithm on simulated data by comparing it to some of the most famous and modern optimization algorithms. The results suggest that our proposed algorithm outperforms these algorithms.


---

## Body

## 1. Introduction

Letℓθ be a log-likelihood function where θ∈ℝd is a parameter vector whose structure will be specified later. In computing the maximum likelihood estimate (MLE) which is the point θ^ at which ℓθ attains its maximum, the very first algorithm one can try is the Newton-Raphson (NR) algorithm which is no longer to be presented. The success of this algorithm comes from an enviable and unequalled property: that of quadratic convergence when the starting guess is well chosen (i.e., close to the unknown solution). However, it also has drawbacks: it is not globally convergent (i.e., its success depends on the starting guess); it can be costly or impossible to implement in high-dimensional problems because of its need of inverting the Hessian matrix at each iteration. Many other algorithms can be used whenever NR cannot. We refer the reader to [1–4] for a comprehensive review of these algorithms. Of all these algorithms, we can mention quasi-Newton (which use approximations of the inverse of the Hessian matrix rather than inverting it), Fisher scoring (a purely statistical method which consists in replacing the Hessian matrix with its mathematical expectation in order to ensure the ascent property and therefore numerical stability), block optimization, and derivative-free optimization algorithms. We can also mention minorization-maximization (MM) algorithms [5, 6] which have made an extraordinary breakthrough in computational statistics and are increasingly used. The MM philosophy for maximizing a function ℓθ is to define, in the first M step, a minorizing function for the objective function, and to maximize, in the second M step, the minorizing function with respect to the parameter vector θ.Even if all these algorithms are considered as performing, they can suffer in some cases for one or many of the following criteria: global convergence (capacity of an algorithm to converge to the true unknown solution from all starting guesses), numerical stability (ascent property), implementation feasibility (for example, algorithms requiring matrix inversion cannot be implemented when the involved matrices are not invertible), low computation time, low computational complexity, and capacity to handle high dimensional problems. The reality is that, in practice, no algorithm is perfect, and for each problem, it is necessary to find the most performing of all existing algorithms or even develop new ones.In this paper, we consider the computing of the maximum likelihood estimate (MLE) of the vector parameter of a statistical model of crash frequencies proposed by [7]. The vector parameter has the form θ=α,βTT where α>0 is the parameter of interest and β is a vector of secondary parameters. Although secondary, the subvector β can contain an important number of components (up to several hundreds) depending on the structure of the data. Thus, the classical algorithms mentioned above may need a great number of iterations (and thus have a slow convergence) or simply fail to converge. For the model considered in this paper, Newton-Raphson and Minorization-Maximization (MM) algorithms have been implemented by [7, 8], but the numerical study of these algorithms has been limited to simple cases. The NR algorithm is known to be fast in simple cases and inefficient (or even impossible to implement) for high-dimensional problems. Moreover, its convergence depends on the starting guess (initial solution θ0 from which the algorithm starts). The MM algorithm, on the other hand, often requires a large number of iterations before converging to the solution.The main contribution of this paper is to build an estimation algorithm which converges faster than the other algorithms and whose convergence is guaranteed and is not affected by the dimension (the number of components ofθ). For this purpose, we exploit the splitting of the parameter vector into two blocks, and we apply the profile likelihood principle (see for example [9], p. 231) to reduce the search of θ to only the search of the parameter of interest α. Then, we develop a new estimation algorithm for which we provide an automatic starting guess from which convergence and numerical stability (ascent property) are guaranteed. This automatization of the starting guess reveals to be a true advantage for our algorithm over the others because it allows to circumvent the difficulty of finding an adequate starting guess. We show using simulated data that our algorithm outperforms Minorization-Maximization (MM) and Newton-Raphson algorithms which are two of the most famous and modern optimization algorithms.The remainder of this paper is organized as follows. In Section2, we describe the data and the estimation problem. The statistical model and the constrained maximum likelihood estimation of the parameters are also presented. In Section 3, we present the profile likelihood principle and then use it to design our new profile-likelihood-based algorithm (PLBA) for computing the MLE of θ. We also provide an automatic starting guess, and we prove that it guarantees convergence of the proposed algorithm. In Section 4, we prove that the proposed algorithm satisfies the ascent property. In Section 5, we report the results of the comparison of the PLBA with MM and NR algorithms. The paper finishes with some discussions and concluding remarks in Section 6.


## 2. Data, Statistical Model, and Problem Setup

The framework of this paper is the statistical analysis of crash data before and after the implementation of a road safety measure ats (s>0) experimental sites (called treated sites) where crashes are classified by severity in r (r>0) levels. This analysis is aimed at estimating the mean effect α>0 of the safety measure simultaneously on all the s sites. This mean effect must be understood in the multiplicative sense and therefore must be compared to 1. A value α<1 indicates a positive effect of the measure (an average reduction of 100×1−α% in the number of crashes) while a value α>1 indicates a negative effect of the measure (an average increase of 100×α−1% in the number of crashes) and α=1 indicates that the measure had no significant influence on the number of crashes.Let(1)x=x111x121⋯x1r1x211x221⋯x2r1⋮⋮⋯⋮⋮⋮⋯⋮x11kx12k⋯x1rkx21kx22k⋯x2rk⋮⋮⋯⋮⋮⋮⋯⋮x11sx12s⋯x1rsx21sx22s⋯x2rsbe a matrix of orders×2r where x1jk (respectively x2jk) is the number of accidents of severity level j (j=1,⋯,r) which occurred at site k in the period before (respectively, after) the implementation of the measure. Also let
(2)z=z11z21⋯zr1⋮⋮⋯⋮z1kz2k⋯zrk⋮⋮⋯⋮z1sz2s⋯zrsbe a matrix of order s×r where zjk (j=1,⋯,r, k=1,⋯,s) is the ratio of the number of accidents of severity j in the “after” period to the number of accidents of the same severity in the “before” period on a control site (a site where the measure has not been implemented) paired with treated site k.LetSr−1=p1,⋯,pr∈0,1r,∑j=1rpj=1 and nk=∑i=12∑j=1rxijk be the total number of accidents observed on the treated site k. N’Guessan et al. [7] proposed the following multinomial-based probability distribution of parameter vector θ for x:
(3)Pxθ=∏k=1snk!∏i=12∏j=1rxijk!∏i=12∏j=1rπijkθzxijk,where θ=α,βTT is the parameter vector, α>0 is the mean effect, β=β1T,⋯,βsTT∈Sr−1s is a vector of sr secondary parameters such that for all k=1,⋯,s,βk=β1k,⋯,βrkT∈Sr−1, and
(4)πijkθz=βjk1+αzk,βk,ifi=1,j=1,⋯,r,αzjkβjk1+αzk,βk,ifi=2,j=1,⋯,r,where zk,βk=∑j′=1rzj′kβj′k.The parameter of the model is the vectorθ=α,βTT∈ℝ+∗×Sr−1s, and the log-likelihood (logPxθ) of observed data x is given to one additive constant by
(5)ℓθ=∑k=1s∑j=1rx•jklogβjk+x2jklogα−x•jklog1+αzk,βk,wherex•jk=x1jk+x2jk (see [7] for more details).In this paper, we are interested in computing the maximum likelihood estimate (MLE)θ^ of the unknown vector parameter θ defined by
(6)θ^=argmaxθ∈ℝ+∗×Sr−1sℓθ.In the cases=1, [10] built an algorithm to solve Equation (6). In the next section, we present a profile-likelihood-based algorithm (PLBA) for computing the MLE in the general case s≥1. Our proposed PLBA is automated since we provide an automatic starting guess which guarantees convergence.


## 3. The Automated Profile-Likelihood-Based Algorithm (PLBA) for Computing the MLE

### 3.1. A Brief Reminder on Profile Likelihood

Let us rewrite the log-likelihood asℓθ=ℓα,β. If, for a given α, the MLE β^ of β may be written as a function β^α of α, that is,
(7)β^=β^α=argmaxβ∈Sr−1sℓα,β,then the profile likelihood function (see for example [9], p. 231) is
(8)ℓpα=ℓα,β^α,expressed as a function ofα only. The maximization of ℓpα is equivalent to that of ℓα,β [11].


### 3.2. Computation ofβ^ in Closed Form

Lemma 1.
Givenα>0, let β^=β∧1T,⋯,β∧sTT be the solution to Equation (7), where for all k=1,⋯,s, β^k=β∧1k,⋯,β∧rkT. The components of β^ are given by
(9)β^jk=x•jk/1+αzjk∑m=1rx•mk/1+αzmk,j=1,⋯,r;k=1,⋯,s.Proof.
Givenα>0, [7] proved that for all k=1,⋯,s, the components of β^k satisfy the following system of r equations:
(10)x•jk−nkβ^jk1+αzjk1+αzk,β^k=0,j=1,⋯,r.
Thus, for allk=1,⋯,s, we apply Theorem 3.5 of [12] to Equation (10) and get
(11)β^k=1nkMα,k/Δα,k−1x•1k1+αz1k,⋯,x•rk1+αzrkT,where
(12)Mα,k=Δα,kαX•kzkT1is a block-defined square matrix of order r+1, Δα,k=diag1+αz1k,⋯,1+αzrk is a diagonal matrix of order r, X•k=x•1k,⋯,x•rkT, and
(13)Mα,k/Δα,k=1nk∑m=1rx•mk1+αzmk
is the Schur complement ofΔα,k in Mα,k (see for example [13] (p. 34) for a reminder on the use of Schur complement for inverting a block matrix). The proof is thus completed.


### 3.3. Profile Likelihood

Theorem 2.
The profile likelihood function is defined, up to an additive constant independent ofα, by
(14)ℓpα=x2••logα−∑k=1s∑j=1rx•jklog1+αzjk,where x2••=∑k=1s∑m=1rx2mk.Proof.
Expression (5) is equivalent to
(15)ℓθ=∑k=1s∑j=1rx•jklogβjk+x2••logα−∑k=1snklog1+αzk,βk.
For allk=1,⋯,s, Equation (9) yields
(16)1+αzk,β^k=1+α∑j=1rzjkβ^jk=1+α∑j=1rzjkx•jk/1+αzjk∑m=1rx•mk/1+αzmk=∑j=1rx•jk/1+αzjk+∑j=1rαzjkx•jk/1+αzjk∑m=1rx•mk/1+αzmk=∑j=1rx•jk1+αzjk/1+αzjk∑m=1rx•mk/1+αzmk=nk∑m=1rx•mk/1+αzmk,and the relationships (9) and (16) enable us to write
(17)ℓpα=∑k=1s∑j=1rx•jklogx•jk1+αzjk−∑k=1s∑j=1rx•jklog∑m=1rx•mk1+αzmk+x2••logα−∑k=1snklognk∑m=1rx•mk/1+αzmk.
After some manipulations on the first, second and fourth terms, we get:(18)ℓpα=∑k=1s∑j=1rx•jklogx•jk−∑k=1s∑j=1rx•jklog1+αzjk−∑k=1snklog∑m=1rx•mk1+αzmk+x2••logα−∑k=1snklognk+∑k=1snklog∑m=1rx•mk1+αzmk.
Removing the third and sixth terms and the constants (first and fifth terms), we get (14).


### 3.4. Design of the PLBA

Lemma 3.
LetF be the mapping defined on ℝ+ by
(19)Fu=−x1••+∑k=1s∑j=1rx•jk1+uzjk,where x1••=∑k=1s∑j=1rx1jk. The MLE α^ of α is the unique root of F.Proof.
(i)
On the one hand, functionF is a one-to-one decreasing function (it is continuous and its derivative F′u is strictly negative for all u≥0) and Fℝ+=F∞,F0, where(20)F∞=limu⟶+∞Fu=−x1••,F0=−x1••+∑k=1s∑j=1rx•jk=−x1••+∑k=1s∑j=1rx1jk+x2jk=−x1••+x1••+x2••=x2••.
Since−x1••<0<x2••, Equation Fα=0 has a unique solution.
(ii)
On the other hand, the MLEα^, if it exists, is solution to the optimization problem(21)α^=argmaxα>0ℓpα.
The profile log-likelihoodℓpα being differentiable for every α>0, the MLE α^ is then solution to Equation ℓp′α=0, where
(22)ℓp′α=x2••α−∑k=1s∑j=1rx•jkzjk1+αzjk=1αx2••−∑k=1s∑j=1rx•jk−x•jk1+αzjk=1α−x1••+∑k=1s∑j=1rx•jk1+αzjk=Fααand ∑k=1s∑j=1rx•jk=x1••+x2••. Thus, Equation ℓp′α^=0 is equivalent to Fα^=0, and α^ is the unique root of F.EquationFu=0 seems fairly complicated to solve in closed form for any s>1 and must therefore be solved numerically. Obviously, there are many root-finding algorithms (see for example [14] (Chapter 3) or [15] (Chapter 3)). Here, we propose a numerical approximation of α^ using the following one-dimensional Newton-Raphson (NR) root-finding algorithm:
(23)αm+1=αm−FαmF′αm,m=0,1,2,⋯,where the starting guessα0 should be chosen in ℝ+ by the user. Our choice of NR algorithm is motivated by the fact that it converges quadratically to the solution if α0 is chosen near the unknown solution. To overcome the difficulty of the choice of α0, we prove that, if we set α0=0 as an automatic starting guess, then, the convergence of NR iterations (23) is always guaranteed.Theorem 4.
Ifα0=0, then, the sequence αm defined by NR iterations (23) converges to the MLE α^.To prove Theorem4, we need the following Lemma 5.Lemma 5.
Letφ be the function linked to NR iterations (23) such that αm+1=φαm, i.e., the function defined for all u∈ℝ+ by
(24)φu=u−FuF′u.(i)
Functionφ is increasing on 0,α^ and decreasing on α^,+∞.(ii)
The MLEα^ is the unique fixed point of φ.(iii)
For allu≥0,(25)φu≥u,ifu≤α^,φu≤u,ifu≥α^.(iv)
For allu≥0, φu≤α^.Proof of Lemma 5.
(i)
For allu≥0, we have(26)F′u=−∑k=1s∑j=1rx•jkzjk1+uzjk2<0.
Therefore,φ is differentiable and
(27)φ′u=1−F′u2−FuF′′uF′u2=FuF′′uF′u2,where for all u∈0,+∞,
(28)F′′u=∑k=1s∑j=1r2x•jkzjk21+uzjk3>0.
So the sign ofφ′u is the same as the one of Fu. As F is a decreasing function and Fα^=0, we have
(29)Fu≥0,ifu≤α^,Fu≤0,ifu≥α^.
It means thatφ is increasing on 0,α^ and decreasing on α^,+∞.
(ii)
Equationφu=u is equivalent to Fu=0, where F is defined by Formula (19). By Lemma 3, this equation has α^ as unique solution; hence, α^ is the unique fixed point of φ.(iii)
For allu≥0, φu−u=−Fu/F′u. The relation (25) is a simple consequence of (26) and (29).(iv)
Letu≥0. If u≤α^, then φu≤φα^=α^ because φ is increasing on 0,α^ (item (i)). If u≥α^, then φu≤φα^=α^ because φ is decreasing on α^,+∞. Thus, for all u≥0, φu≤α^.We may now prove Theorem4.Proof of Theorem 4.
Assume thatα0=0. Then, α0≤α^ and by item (iii) and (iv) of Lemma 5, α0≤φα0=α1≤α^. It can be easily proved by induction that
(30)α0≤α1≤⋯≤αm≤αm+1≤α^.
Thus, the sequenceαm is increasing and bounded; hence, it converges to the unique fixed point of φ which is α^.Remark 6.
Actually, from the proof of Theorem4, it is clear that any starting value α0∈0,α^ will guarantee convergence of the sequence αm generated by NR iterations (23). However, since α^ is unknown, it is difficult to find a value other than α0=0.We therefore propose an algorithm (see Algorithm1) starting from α0=0. The MLE α^ is computed using NR iterations (23); afterwards, β^ is computed from α^.Algorithm 1: PLBA for computing θ.^
Input:x, z and ε>0Output: MLE θ^1 Initializem=0 and α0=0 ;2repeat3αm+1=αm−Fαm/F′αm;4 Setm⟵m+1 ;5untilFαm<ε;6 Setα^=αm ;7 For everyk=1,⋯,s, compute β^k asβ^k=1/∑j=1rx•jk/1+α^zjkx•1k/1+α∧z1k,⋯,x•rk/1+α∧zrkT;8 Setθ^=α∧,β∧1T,⋯,β∧sTT.


## 3.1. A Brief Reminder on Profile Likelihood

Let us rewrite the log-likelihood asℓθ=ℓα,β. If, for a given α, the MLE β^ of β may be written as a function β^α of α, that is,
(7)β^=β^α=argmaxβ∈Sr−1sℓα,β,then the profile likelihood function (see for example [9], p. 231) is
(8)ℓpα=ℓα,β^α,expressed as a function ofα only. The maximization of ℓpα is equivalent to that of ℓα,β [11].


## 3.2. Computation ofβ^ in Closed Form

Lemma 1.
Givenα>0, let β^=β∧1T,⋯,β∧sTT be the solution to Equation (7), where for all k=1,⋯,s, β^k=β∧1k,⋯,β∧rkT. The components of β^ are given by
(9)β^jk=x•jk/1+αzjk∑m=1rx•mk/1+αzmk,j=1,⋯,r;k=1,⋯,s.Proof.
Givenα>0, [7] proved that for all k=1,⋯,s, the components of β^k satisfy the following system of r equations:
(10)x•jk−nkβ^jk1+αzjk1+αzk,β^k=0,j=1,⋯,r.
Thus, for allk=1,⋯,s, we apply Theorem 3.5 of [12] to Equation (10) and get
(11)β^k=1nkMα,k/Δα,k−1x•1k1+αz1k,⋯,x•rk1+αzrkT,where
(12)Mα,k=Δα,kαX•kzkT1is a block-defined square matrix of order r+1, Δα,k=diag1+αz1k,⋯,1+αzrk is a diagonal matrix of order r, X•k=x•1k,⋯,x•rkT, and
(13)Mα,k/Δα,k=1nk∑m=1rx•mk1+αzmk
is the Schur complement ofΔα,k in Mα,k (see for example [13] (p. 34) for a reminder on the use of Schur complement for inverting a block matrix). The proof is thus completed.


## 3.3. Profile Likelihood

Theorem 2.
The profile likelihood function is defined, up to an additive constant independent ofα, by
(14)ℓpα=x2••logα−∑k=1s∑j=1rx•jklog1+αzjk,where x2••=∑k=1s∑m=1rx2mk.Proof.
Expression (5) is equivalent to
(15)ℓθ=∑k=1s∑j=1rx•jklogβjk+x2••logα−∑k=1snklog1+αzk,βk.
For allk=1,⋯,s, Equation (9) yields
(16)1+αzk,β^k=1+α∑j=1rzjkβ^jk=1+α∑j=1rzjkx•jk/1+αzjk∑m=1rx•mk/1+αzmk=∑j=1rx•jk/1+αzjk+∑j=1rαzjkx•jk/1+αzjk∑m=1rx•mk/1+αzmk=∑j=1rx•jk1+αzjk/1+αzjk∑m=1rx•mk/1+αzmk=nk∑m=1rx•mk/1+αzmk,and the relationships (9) and (16) enable us to write
(17)ℓpα=∑k=1s∑j=1rx•jklogx•jk1+αzjk−∑k=1s∑j=1rx•jklog∑m=1rx•mk1+αzmk+x2••logα−∑k=1snklognk∑m=1rx•mk/1+αzmk.
After some manipulations on the first, second and fourth terms, we get:(18)ℓpα=∑k=1s∑j=1rx•jklogx•jk−∑k=1s∑j=1rx•jklog1+αzjk−∑k=1snklog∑m=1rx•mk1+αzmk+x2••logα−∑k=1snklognk+∑k=1snklog∑m=1rx•mk1+αzmk.
Removing the third and sixth terms and the constants (first and fifth terms), we get (14).


## 3.4. Design of the PLBA

Lemma 3.
LetF be the mapping defined on ℝ+ by
(19)Fu=−x1••+∑k=1s∑j=1rx•jk1+uzjk,where x1••=∑k=1s∑j=1rx1jk. The MLE α^ of α is the unique root of F.Proof.
(i)
On the one hand, functionF is a one-to-one decreasing function (it is continuous and its derivative F′u is strictly negative for all u≥0) and Fℝ+=F∞,F0, where(20)F∞=limu⟶+∞Fu=−x1••,F0=−x1••+∑k=1s∑j=1rx•jk=−x1••+∑k=1s∑j=1rx1jk+x2jk=−x1••+x1••+x2••=x2••.
Since−x1••<0<x2••, Equation Fα=0 has a unique solution.
(ii)
On the other hand, the MLEα^, if it exists, is solution to the optimization problem(21)α^=argmaxα>0ℓpα.
The profile log-likelihoodℓpα being differentiable for every α>0, the MLE α^ is then solution to Equation ℓp′α=0, where
(22)ℓp′α=x2••α−∑k=1s∑j=1rx•jkzjk1+αzjk=1αx2••−∑k=1s∑j=1rx•jk−x•jk1+αzjk=1α−x1••+∑k=1s∑j=1rx•jk1+αzjk=Fααand ∑k=1s∑j=1rx•jk=x1••+x2••. Thus, Equation ℓp′α^=0 is equivalent to Fα^=0, and α^ is the unique root of F.EquationFu=0 seems fairly complicated to solve in closed form for any s>1 and must therefore be solved numerically. Obviously, there are many root-finding algorithms (see for example [14] (Chapter 3) or [15] (Chapter 3)). Here, we propose a numerical approximation of α^ using the following one-dimensional Newton-Raphson (NR) root-finding algorithm:
(23)αm+1=αm−FαmF′αm,m=0,1,2,⋯,where the starting guessα0 should be chosen in ℝ+ by the user. Our choice of NR algorithm is motivated by the fact that it converges quadratically to the solution if α0 is chosen near the unknown solution. To overcome the difficulty of the choice of α0, we prove that, if we set α0=0 as an automatic starting guess, then, the convergence of NR iterations (23) is always guaranteed.Theorem 4.
Ifα0=0, then, the sequence αm defined by NR iterations (23) converges to the MLE α^.To prove Theorem4, we need the following Lemma 5.Lemma 5.
Letφ be the function linked to NR iterations (23) such that αm+1=φαm, i.e., the function defined for all u∈ℝ+ by
(24)φu=u−FuF′u.(i)
Functionφ is increasing on 0,α^ and decreasing on α^,+∞.(ii)
The MLEα^ is the unique fixed point of φ.(iii)
For allu≥0,(25)φu≥u,ifu≤α^,φu≤u,ifu≥α^.(iv)
For allu≥0, φu≤α^.Proof of Lemma 5.
(i)
For allu≥0, we have(26)F′u=−∑k=1s∑j=1rx•jkzjk1+uzjk2<0.
Therefore,φ is differentiable and
(27)φ′u=1−F′u2−FuF′′uF′u2=FuF′′uF′u2,where for all u∈0,+∞,
(28)F′′u=∑k=1s∑j=1r2x•jkzjk21+uzjk3>0.
So the sign ofφ′u is the same as the one of Fu. As F is a decreasing function and Fα^=0, we have
(29)Fu≥0,ifu≤α^,Fu≤0,ifu≥α^.
It means thatφ is increasing on 0,α^ and decreasing on α^,+∞.
(ii)
Equationφu=u is equivalent to Fu=0, where F is defined by Formula (19). By Lemma 3, this equation has α^ as unique solution; hence, α^ is the unique fixed point of φ.(iii)
For allu≥0, φu−u=−Fu/F′u. The relation (25) is a simple consequence of (26) and (29).(iv)
Letu≥0. If u≤α^, then φu≤φα^=α^ because φ is increasing on 0,α^ (item (i)). If u≥α^, then φu≤φα^=α^ because φ is decreasing on α^,+∞. Thus, for all u≥0, φu≤α^.We may now prove Theorem4.Proof of Theorem 4.
Assume thatα0=0. Then, α0≤α^ and by item (iii) and (iv) of Lemma 5, α0≤φα0=α1≤α^. It can be easily proved by induction that
(30)α0≤α1≤⋯≤αm≤αm+1≤α^.
Thus, the sequenceαm is increasing and bounded; hence, it converges to the unique fixed point of φ which is α^.Remark 6.
Actually, from the proof of Theorem4, it is clear that any starting value α0∈0,α^ will guarantee convergence of the sequence αm generated by NR iterations (23). However, since α^ is unknown, it is difficult to find a value other than α0=0.We therefore propose an algorithm (see Algorithm1) starting from α0=0. The MLE α^ is computed using NR iterations (23); afterwards, β^ is computed from α^.Algorithm 1: PLBA for computing θ.^
Input:x, z and ε>0Output: MLE θ^1 Initializem=0 and α0=0 ;2repeat3αm+1=αm−Fαm/F′αm;4 Setm⟵m+1 ;5untilFαm<ε;6 Setα^=αm ;7 For everyk=1,⋯,s, compute β^k asβ^k=1/∑j=1rx•jk/1+α^zjkx•1k/1+α∧z1k,⋯,x•rk/1+α∧zrkT;8 Setθ^=α∧,β∧1T,⋯,β∧sTT.


## 4. Ascent Property

The ascent property of Algorithm1 (the fact that the profile log-likelihood is increased monotonically by the algorithm) is given by Theorem 7.Theorem 7.
The sequenceαm generated by Algorithm 1 increases monotonically the profile log-likelihood ℓpα, that is
(31)ℓpαm+1≥ℓpαm,m=0,1,⋯Proof.
From (22) and (29), we deduce that, for all α>0, ℓp′α≥0 if α≤α^ and ℓp′α≤0 if α≥α^. Hence, the function ℓp is increasing on the interval 0,α^ and decreasing on α^,+∞. From (30), we know that the sequence αm is increasing and still belongs to the interval 0,α^. Then, for all iteration m, we have αm≤αm+1 and ℓpαm≤ℓpαm+1 because ℓp is increasing on 0,α^. The proof of Theorem 7 is then completed.


## 5. Simulation Study

We compare the performance of our PLBA with that of Newton-Raphson (NR) and MM algorithms in R software [16]. We implemented the NR algorithm using the R package pracma [17] and the MM algorithm using Theorem 3.3 of [8]. The choice of these two comparison algorithms is motivated by the following factors: (a) MM and NR algorithms are two of the most used algorithms in statistics for parameter estimation; (b) other algorithms such as quasi-Newton and derivative-free algorithms showed very low convergence proportions and their results are not reported here; (c) the results obtained for the particular case s=1 in [10] suggest that MM and NR algorithms are much more efficient than quasi-Newton and derivative-free algorithms.


### 5.1. Data Generation Principle

For given values ofs, r, n1, ..., ns, the components of matrix z are randomly generated from a uniform distribution on 0.5,2.5. Given the true parameter vector denoted θ0=α0,β10T,⋯,βs0TT, the true class probabilities πijkθ0z (see Equation (4)) are computed afterwards data matrix x (Equation (1)) is generated using the random generation function linked to the probability distribution function (3).The true parameter vectorθ0=α0,β10T,⋯,βs0TT is presented under the following six scenarios:Scenario 1:s=2, r=3,
(32)α0=0.85,β10=0.52,0.31,0.17T,β20=0.25,0.45,0.30T.Scenario 2:s=8, r=3,
(33)α0=0.85,βk0=0.50,0.15,0.35T,k∈1,3,5,βk0=0.43,0.32,0.25T,k∈2,4,7,βk0=0.35,0.35,0.30T,k∈6,8.Scenario 3:s=14, r=3,
(34)α0=1.02,βk0=0.50,0.15,0.35T,k∈1,3,5,9,11,βk0=0.43,0.32,0.25T,k∈2,4,7,12,13,βk0=0.35,0.35,0.30T,k∈6,8,10,14.Scenario 4:s=15, r=4,
(35)α0=1.02,βk0=0.40,0.10,0.30,0.20T,k∈1,3,5,9,13,14,βk0=0.45,0.10,0.25,0.20T,k∈2,4,7,12,βk0=0.25,⋯,0.25⏟4T,k∈6,8,10,11,15.Scenario 5:s=20, r=4,
(36)α0=1.25,βk0=0.62,0.16,0.07,0.15T,k∈1,3,5,9,11,13,15,19,βk0=0.50,0.15,0.10,0.25T,k∈2,4,7,12,14,17,βk0=0.25,⋯,0.25⏟4T,k∈6,8,10,16,18,20.Scenario 6:s=20, r=5,
(37)α0=1.25,βk0=0.63,0.16,0.06,0.05,0.10T,k∈1,3,5,9,11,13,15,19,βk0=0.32,0.18,0.25,0.10,0.15T,k∈2,4,7,12,14,17,βk0=0.20,⋯,0.20⏟5T,k∈6,8,10,16,18,20.For these different scenarios, the number of parameters (1+sr) is given by Table 1.Table 1
Number of parameters (1+sr) for the different scenarios.
Scenario123456Number of parameters725436181101For thenk’s (k=1,⋯,s), we have chosen two common values: a low value (n=50) and a great value (n=5000). Except for the proposed algorithm whose starting guess is automated, the other algorithms were given randomly generated starting guesses.


### 5.2. Results

Tables2–7 present the average results obtained for the different scenarios over 1000 replications (i.e., 1000 repetitions of the data generation and computation of θ^). In these tables, CPU times are given in seconds and CPU time ratios are calculated as the ratio of the mean CPU time of a given algorithm to the CPU time of the PLBA. Thus, the CPU time ratio of the PLBA is always equal to 1. The Mean Square Error (MSE) is defined as
(38)MSEθ^,θ0=11+srα∧−α02+∑k=1s∑j=1rβ∧jk−βjk02.Table 2
Results for Scenario 1 (values in brackets are standard deviations).
PLBAMMNRn=50α0.864 (0.180)0.864 (0.180)0.878 (0.181)β110.515 (0.071)0.515 (0.071)0.513 (0.071)β210.310 (0.066)0.310 (0.066)0.309 (0.067)β310.175 (0.054)0.175 (0.054)0.179 (0.054)β120.248 (0.062)0.248 (0.062)0.249 (0.063)β220.453 (0.068)0.453 (0.068)0.448 (0.068)β320.299 (0.064)0.299 (0.064)0.302 (0.064)Convergence proportion (%)10010054.3Iterations5.7 (0.5)27.8 (5.6)7 (1.2)CPU time (secs)0.00050.00430.0029Time ratio186Log-likelihood−191.95−191.95−191.95MSE8.2e−038.2e−038.4e−03n=5000α0.851 (0.017)0.851 (0.017)0.852 (0.017)β110.520 (0.007)0.520 (0.007)0.520 (0.007)β210.310 (0.007)0.310 (0.007)0.310 (0.007)β310.170 (0.005)0.170 (0.005)0.170 (0.005)β120.250 (0.006)0.250 (0.006)0.250 (0.006)β220.450 (0.007)0.450 (0.007)0.450 (0.007)β320.300 (0.007)0.300 (0.007)0.300 (0.007)Convergence proportion (%)10010055.9Iterations6 (0.1)36.3 (4.5)7 (1.3)CPU time (secs)0.00040.00560.0030Time ratio1147Log-likelihood−19380−19380−19380MSE7.8e−057.8e−057.6e−05Table 3
Results for Scenario 2 (values in brackets are standard deviations).
PLBAMMNRn=50Convergence proportion (%)10010035.5Iterations6 (0.1)29.9 (3)7.7 (1.3)CPU time (secs)0.00050.00870.0068Time ratio11814Log-likelihood−773.57−773.57−773.57MSE4.4e−034.4e−034.4e−03n=5000Convergence proportion (%)10010038.7Iterations6 (0)38.6 (2.6)7.7 (1.4)CPU time (secs)0.00050.01080.0069Time ratio12114Log-likelihood−78015.32−78015.32−78015.32MSE4.5e−054.5e−054.5e−05Table 4
Results for Scenario 3 (values in brackets are standard deviations).
PLBAMMNRn=50Convergence proportion (%)10010035.5Iterations6 (0)35.2 (3)8.2 (1.6)CPU time (secs)0.00050.01340.0113Time ratio12521Log-likelihood−1356.98−1356.98−1356.98MSE4.4e−034.4e−034.5e−03n=5000Convergence proportion (%)10010037.1Iterations6 (0)45.4 (2.5)8.2 (1.6)CPU time (secs)0.00050.01600.0106Time ratio13523Log-likelihood−137298.2−137298.2−137298.2MSE4.5e−054.5e−054.6e−05Table 5
Results for Scenario 4 (values in brackets are standard deviations).
PLBAMMNRn=50Convergence proportion (%)10010029.6Iterations6 (0)34.9 (2.9)8.3 (1.3)CPU time (secs)0.00060.01350.0149Time ratio12426Log-likelihood−1647.49−1647.49−1647.49MSE3.6e−033.6e−033.6e−03n=5000Convergence proportion (%)10010029.4Iterations6 (0)45.4 (2.6)8.3 (1.4)CPU time (secs)0.00050.01660.0142Time ratio13127Log-likelihood−166241.8−166241.8−166241.8MSE3.8e−053.8e−053.7e−05Table 6
Results for Scenario 5 (Values in brackets are standard deviations).
PLBAMMNRn=50Convergence proportion (%)10010025.9Iterations6 (0)42.2 (3.2)8.7 (1.8)CPU time (secs)0.00050.01810.0213Time ratio13642Log-likelihood−2102.95−2102.95−2102.95MSE3.2e−033.2e−033.3e−03n=5000Convergence proportion (%)10010026.3Iterations6.6 (0.5)56.2 (2.9)8.5 (1.6)CPU time (secs)0.00060.02340.0196Time ratio13731Log-likelihood−210936.75−210936.75−210936.75MSE3.5e−053.5e−053.5e−05Table 7
Results for Scenario 6 (values in brackets are standard deviations).
PLBAMMNRn=50Convergence proportion (%)10010019.6Iterations6 (0)41.1 (3.1)8.8 (1.4)CPU time (secs)0.00040.01770.0307Time ratio14374Log-likelihood−2308.19−2308.19−2308.19MSE2.8e−032.8e−032.8e−03n=5000Convergence proportion (%)10010021.9Iterations6.6 (0.5)56.4 (3.3)8.7 (1.6)CPU time (secs)0.00050.02270.0278Time ratio15061Log-likelihood−230371.92−230371.92−230371.92MSE3e−053e−053e−05Due to the increasing number of parameters, the MLEθ^ has only been included for Scenario 1 (Table 2).From these tables, it can be seen that PLBA and MM algorithms have always converged while the convergence proportion of NR algorithm decreases (from 55.9% to 19.6%) with the number of parameters (see Figure1). As far as the MSE is concerned, the trend for all the algorithms is that the MSE decreases when the sample size n increases.Figure 1
Convergence proportions for NR.By taking a look at the CPU times ratios, we notice that the CPU time ratios of the MM and NR algorithms are all well above 1. This means that the PLBA is significantly faster than these two algorithms. As shown on Figures2 and 3, the CPU time ratios of MM and NR algorithms increase with the number of parameters. The PLBA is on average 8 to 50 times faster than MM and 6 to 74 times quicker than NR.Figure 2
CPU time ratios for NR and MM algorithms (n=50).Figure 3
CPU time ratios for NR and MM algorithms (n=5000).


### 5.3. Analysis of the Results

It appears that our proposed PLBA outperforms the Minorization-Maximization (MM) and the full Newton-Raphson (NR) algorithms. It always converges; it requires a low number of iterations and little computation time. The number of iterations varies slightly and seems to stabilize around six iterations whatever the starting guess and the number of parameters to be estimated. The MM algorithm also has a convergence rate of 100%, but its number of iterations is quite high, and this may indicate some sensitivity to the starting guess. The NR algorithm has a convergence proportion at most slightly higher than 50% and this proportion decreases when the number of parameters increases. This is not surprising because at each iteration of the NR algorithm, a square matrix of order1+sr is numerically inverted. This numerical inversion can be complicated or impossible when the matrix to be inverted is ill-conditioned or singular. The fact that NR does not converge for all the replications is also not surprising since it is well known that NR may converge to bad values or may not converge at all if the starting guess is far from the true parameter vector.The good results obtained by our proposed algorithm PLBA can be explained by several factors. First, the principle of profile likelihood enables to reduce the search for1+sr solutions to that of the single solution α^ from which the remaining sr estimates β^jk (j=1,⋯,r, k=1,⋯,s) can be obtained by a simple formula. This also enables to reduce the computation time required by our proposed PLBA to converge. Secondly, the fact that the estimation of α relies on a one-dimensional NR enables our algorithm to enjoy the quadratic convergence of the NR algorithm. Thirdly, the definition of the automated starting guess ensures the convergence of our proposed algorithm PLBA, and the ascension property ensures its numerical stability.


## 5.1. Data Generation Principle

For given values ofs, r, n1, ..., ns, the components of matrix z are randomly generated from a uniform distribution on 0.5,2.5. Given the true parameter vector denoted θ0=α0,β10T,⋯,βs0TT, the true class probabilities πijkθ0z (see Equation (4)) are computed afterwards data matrix x (Equation (1)) is generated using the random generation function linked to the probability distribution function (3).The true parameter vectorθ0=α0,β10T,⋯,βs0TT is presented under the following six scenarios:Scenario 1:s=2, r=3,
(32)α0=0.85,β10=0.52,0.31,0.17T,β20=0.25,0.45,0.30T.Scenario 2:s=8, r=3,
(33)α0=0.85,βk0=0.50,0.15,0.35T,k∈1,3,5,βk0=0.43,0.32,0.25T,k∈2,4,7,βk0=0.35,0.35,0.30T,k∈6,8.Scenario 3:s=14, r=3,
(34)α0=1.02,βk0=0.50,0.15,0.35T,k∈1,3,5,9,11,βk0=0.43,0.32,0.25T,k∈2,4,7,12,13,βk0=0.35,0.35,0.30T,k∈6,8,10,14.Scenario 4:s=15, r=4,
(35)α0=1.02,βk0=0.40,0.10,0.30,0.20T,k∈1,3,5,9,13,14,βk0=0.45,0.10,0.25,0.20T,k∈2,4,7,12,βk0=0.25,⋯,0.25⏟4T,k∈6,8,10,11,15.Scenario 5:s=20, r=4,
(36)α0=1.25,βk0=0.62,0.16,0.07,0.15T,k∈1,3,5,9,11,13,15,19,βk0=0.50,0.15,0.10,0.25T,k∈2,4,7,12,14,17,βk0=0.25,⋯,0.25⏟4T,k∈6,8,10,16,18,20.Scenario 6:s=20, r=5,
(37)α0=1.25,βk0=0.63,0.16,0.06,0.05,0.10T,k∈1,3,5,9,11,13,15,19,βk0=0.32,0.18,0.25,0.10,0.15T,k∈2,4,7,12,14,17,βk0=0.20,⋯,0.20⏟5T,k∈6,8,10,16,18,20.For these different scenarios, the number of parameters (1+sr) is given by Table 1.Table 1
Number of parameters (1+sr) for the different scenarios.
Scenario123456Number of parameters725436181101For thenk’s (k=1,⋯,s), we have chosen two common values: a low value (n=50) and a great value (n=5000). Except for the proposed algorithm whose starting guess is automated, the other algorithms were given randomly generated starting guesses.


## 5.2. Results

Tables2–7 present the average results obtained for the different scenarios over 1000 replications (i.e., 1000 repetitions of the data generation and computation of θ^). In these tables, CPU times are given in seconds and CPU time ratios are calculated as the ratio of the mean CPU time of a given algorithm to the CPU time of the PLBA. Thus, the CPU time ratio of the PLBA is always equal to 1. The Mean Square Error (MSE) is defined as
(38)MSEθ^,θ0=11+srα∧−α02+∑k=1s∑j=1rβ∧jk−βjk02.Table 2
Results for Scenario 1 (values in brackets are standard deviations).
PLBAMMNRn=50α0.864 (0.180)0.864 (0.180)0.878 (0.181)β110.515 (0.071)0.515 (0.071)0.513 (0.071)β210.310 (0.066)0.310 (0.066)0.309 (0.067)β310.175 (0.054)0.175 (0.054)0.179 (0.054)β120.248 (0.062)0.248 (0.062)0.249 (0.063)β220.453 (0.068)0.453 (0.068)0.448 (0.068)β320.299 (0.064)0.299 (0.064)0.302 (0.064)Convergence proportion (%)10010054.3Iterations5.7 (0.5)27.8 (5.6)7 (1.2)CPU time (secs)0.00050.00430.0029Time ratio186Log-likelihood−191.95−191.95−191.95MSE8.2e−038.2e−038.4e−03n=5000α0.851 (0.017)0.851 (0.017)0.852 (0.017)β110.520 (0.007)0.520 (0.007)0.520 (0.007)β210.310 (0.007)0.310 (0.007)0.310 (0.007)β310.170 (0.005)0.170 (0.005)0.170 (0.005)β120.250 (0.006)0.250 (0.006)0.250 (0.006)β220.450 (0.007)0.450 (0.007)0.450 (0.007)β320.300 (0.007)0.300 (0.007)0.300 (0.007)Convergence proportion (%)10010055.9Iterations6 (0.1)36.3 (4.5)7 (1.3)CPU time (secs)0.00040.00560.0030Time ratio1147Log-likelihood−19380−19380−19380MSE7.8e−057.8e−057.6e−05Table 3
Results for Scenario 2 (values in brackets are standard deviations).
PLBAMMNRn=50Convergence proportion (%)10010035.5Iterations6 (0.1)29.9 (3)7.7 (1.3)CPU time (secs)0.00050.00870.0068Time ratio11814Log-likelihood−773.57−773.57−773.57MSE4.4e−034.4e−034.4e−03n=5000Convergence proportion (%)10010038.7Iterations6 (0)38.6 (2.6)7.7 (1.4)CPU time (secs)0.00050.01080.0069Time ratio12114Log-likelihood−78015.32−78015.32−78015.32MSE4.5e−054.5e−054.5e−05Table 4
Results for Scenario 3 (values in brackets are standard deviations).
PLBAMMNRn=50Convergence proportion (%)10010035.5Iterations6 (0)35.2 (3)8.2 (1.6)CPU time (secs)0.00050.01340.0113Time ratio12521Log-likelihood−1356.98−1356.98−1356.98MSE4.4e−034.4e−034.5e−03n=5000Convergence proportion (%)10010037.1Iterations6 (0)45.4 (2.5)8.2 (1.6)CPU time (secs)0.00050.01600.0106Time ratio13523Log-likelihood−137298.2−137298.2−137298.2MSE4.5e−054.5e−054.6e−05Table 5
Results for Scenario 4 (values in brackets are standard deviations).
PLBAMMNRn=50Convergence proportion (%)10010029.6Iterations6 (0)34.9 (2.9)8.3 (1.3)CPU time (secs)0.00060.01350.0149Time ratio12426Log-likelihood−1647.49−1647.49−1647.49MSE3.6e−033.6e−033.6e−03n=5000Convergence proportion (%)10010029.4Iterations6 (0)45.4 (2.6)8.3 (1.4)CPU time (secs)0.00050.01660.0142Time ratio13127Log-likelihood−166241.8−166241.8−166241.8MSE3.8e−053.8e−053.7e−05Table 6
Results for Scenario 5 (Values in brackets are standard deviations).
PLBAMMNRn=50Convergence proportion (%)10010025.9Iterations6 (0)42.2 (3.2)8.7 (1.8)CPU time (secs)0.00050.01810.0213Time ratio13642Log-likelihood−2102.95−2102.95−2102.95MSE3.2e−033.2e−033.3e−03n=5000Convergence proportion (%)10010026.3Iterations6.6 (0.5)56.2 (2.9)8.5 (1.6)CPU time (secs)0.00060.02340.0196Time ratio13731Log-likelihood−210936.75−210936.75−210936.75MSE3.5e−053.5e−053.5e−05Table 7
Results for Scenario 6 (values in brackets are standard deviations).
PLBAMMNRn=50Convergence proportion (%)10010019.6Iterations6 (0)41.1 (3.1)8.8 (1.4)CPU time (secs)0.00040.01770.0307Time ratio14374Log-likelihood−2308.19−2308.19−2308.19MSE2.8e−032.8e−032.8e−03n=5000Convergence proportion (%)10010021.9Iterations6.6 (0.5)56.4 (3.3)8.7 (1.6)CPU time (secs)0.00050.02270.0278Time ratio15061Log-likelihood−230371.92−230371.92−230371.92MSE3e−053e−053e−05Due to the increasing number of parameters, the MLEθ^ has only been included for Scenario 1 (Table 2).From these tables, it can be seen that PLBA and MM algorithms have always converged while the convergence proportion of NR algorithm decreases (from 55.9% to 19.6%) with the number of parameters (see Figure1). As far as the MSE is concerned, the trend for all the algorithms is that the MSE decreases when the sample size n increases.Figure 1
Convergence proportions for NR.By taking a look at the CPU times ratios, we notice that the CPU time ratios of the MM and NR algorithms are all well above 1. This means that the PLBA is significantly faster than these two algorithms. As shown on Figures2 and 3, the CPU time ratios of MM and NR algorithms increase with the number of parameters. The PLBA is on average 8 to 50 times faster than MM and 6 to 74 times quicker than NR.Figure 2
CPU time ratios for NR and MM algorithms (n=50).Figure 3
CPU time ratios for NR and MM algorithms (n=5000).


## 5.3. Analysis of the Results

It appears that our proposed PLBA outperforms the Minorization-Maximization (MM) and the full Newton-Raphson (NR) algorithms. It always converges; it requires a low number of iterations and little computation time. The number of iterations varies slightly and seems to stabilize around six iterations whatever the starting guess and the number of parameters to be estimated. The MM algorithm also has a convergence rate of 100%, but its number of iterations is quite high, and this may indicate some sensitivity to the starting guess. The NR algorithm has a convergence proportion at most slightly higher than 50% and this proportion decreases when the number of parameters increases. This is not surprising because at each iteration of the NR algorithm, a square matrix of order1+sr is numerically inverted. This numerical inversion can be complicated or impossible when the matrix to be inverted is ill-conditioned or singular. The fact that NR does not converge for all the replications is also not surprising since it is well known that NR may converge to bad values or may not converge at all if the starting guess is far from the true parameter vector.The good results obtained by our proposed algorithm PLBA can be explained by several factors. First, the principle of profile likelihood enables to reduce the search for1+sr solutions to that of the single solution α^ from which the remaining sr estimates β^jk (j=1,⋯,r, k=1,⋯,s) can be obtained by a simple formula. This also enables to reduce the computation time required by our proposed PLBA to converge. Secondly, the fact that the estimation of α relies on a one-dimensional NR enables our algorithm to enjoy the quadratic convergence of the NR algorithm. Thirdly, the definition of the automated starting guess ensures the convergence of our proposed algorithm PLBA, and the ascension property ensures its numerical stability.


## 6. Conclusion

In this article, we have built a profile-likelihood-based algorithm (PLBA) to compute, under constraints, the maximum likelihood estimate (MLE) of the parameter vector of a statistical model used in the analysis of a road safety measure applied tos sites presenting in total r accident severity levels. The parameter vector of the model is of the form θ=α,βTT, where α is the parameter of interest and β is a vector of sr secondary parameters. Using the likelihood equations, we obtained the closed-form expression of the components of β as a function of the main parameter α and then used the principle of profile likelihood to express the log-likelihood only as a function of α. We then built an algorithm mixing a one-dimensional Newton method to compute the estimate of the parameter of interest α and the computation of the secondary parameters from their closed-form expressions. The starting guess of our proposed algorithm is automated in such a way as to guarantee its convergence towards the MLE θ^. The numerical studies suggest that our PLBA outperforms the Minorization-Maximization (MM) and the full Newton-Raphson (NR) algorithms in terms of computation time. Our PLBA converges to estimates close to the true values of the parameter vector even for small sample sizes. They also suggest that the problem addressed in this paper is difficult to tackle for the full NR algorithm (the latter having a convergence proportion at most slightly higher than 50%).



---

*Source: 6974166-2022-10-26.xml*