Sequential change-point detection in a multinomial logistic regression model

Fuxiao Li; Zhanshou Chen; Yanting Xiao

doi:10.1515/math-2020-0037

Article Open Access

Sequential change-point detection in a multinomial logistic regression model

Fuxiao Li , Zhanshou Chen and Yanting Xiao

Published/Copyright: July 29, 2020

Published by

Become an author with De Gruyter Brill

Submit Manuscript Author Information Explore this Subject

$Open Mathematics$

From the journal Open Mathematics Volume 18 Issue 1

Abstract

Change-point detection in categorical time series has recently gained attention as statistical models incorporating change-points are common in practice, especially in the area of biomedicine. In this article, we propose a sequential change-point detection procedure based on the partial likelihood score process for the detection of changes in the coefficients of multinomial logistic regression model. The asymptotic results are presented under both the null of no change and the alternative of changes in coefficients. We carry out a Monte Carlo experiment to evaluate the empirical size of the proposed procedure as well as its average run length. We illustrate the method by using data on a DNA sequence. Monte Carlo experiments and real data analysis demonstrate the effectiveness of the proposed procedure.

Keywords: categorical time series; sequential change-point detection; multinomial logistic regression model; partial likelihood score process

MSC 2010: 60F15; 62G30

1 Introduction

Structural change is of central importance in many fields of research and data analysis. In the past few decades, studies of change-point detection have attracted a great deal of attention in fields such as statistics, engineering, economics, climatology and bioscience. For surveys, we refer to the studies by Csörgö and Horváth [1], Perron [2], Gombay [3], Chen and Tian [4], Na et al. [5], Zou et al. [6], Ross [7], Robbins et al. [8], Li et al. [9], Cao et al. [10] and Chen [11]. Most of these studies deal with the classical linear regression model either with independent identical distribution errors or assuming some type of dependence. However, the research of detection of change-point in generalized linear models has received increased attention in recent years. Antoch et al. [12] presented some results on testing for changes in generalized linear models based on overall maximum-type test statistic. Zhou and Liang [13] proposed an estimating procedure for estimating the change-point along with other regression coefficients under the generalized linear model framework.

Categorical time series data are frequently encountered in biomedicine, social science and genetics. As generalized linear regression models for categorical time series allow for parsimonious modeling and incorporation of random time-dependent covariates, Fokianos and Kedem [14] suggested the generalized linear model for categorical time series modeling. For change-point detection in categorical time series, there are retrospective procedures which deal with the detection of a structural change within an observed data set of fixed size, whereas sequential procedures check the stability hypothesis each time a new observation is available (see Horváth et al. [15]). Fokianos et al. [16] provided a new statistical procedure based on the partial likelihood score process for the retrospective detection of change in binary time series. Gombay et al. [17] discussed the retrospective detection of change in categorical time series using the same method. Hudecová [18] studied the change-point problem within the framework of autoregressive models for binary time series. Wang et al. [19] proposed a novel change-point detection procedure motivated by high-dimensional homogeneity tests to estimate the locations of multiple change-points in multinomial data with a large number of categories.

The above study of detection of change is retrospective change-point detection, in regard to sequential change-point detection (on-line monitoring, sequential test or priori test) whereby data are not observed at once, but arrive in a sequential manner – one by one, Xia et al. [20] introduced two procedures to sequentially detect structural change in generalized linear models with assuming independence. Höhle [21] proposed a CUSUM control chart method based on the generalized likelihood ratio statistic for sequential change-point detection in regression models for categorical time series.

The goal of this article is to propose a sequential change-point detection procedure based on the partial likelihood score process for the coefficients of multinomial logistic regression model for categorical time series. Score test for the sequential detection of changes in time series models has been studied by many authors. Related work can be found in the studies by Gombay and Serban [22] and Gombay et al. [23]. Compared with the CUSUM detection method, the alternative value of the parameter of the proposed method does not have to be estimated, and in our current model the parameter estimation is complex. Also, the probability of false alarms of the proposed score test is lower.

In this article, we propose a sequential test statistic based on the partial likelihood score process for the coefficients of multinomial logistic regression model. The asymptotic distribution of test statistic is derived under the null hypothesis and the consistency is proven under the alternative hypothesis. We perform Monte Carlo simulations to explore the finite sample performance of the proposed test statistics in terms of empirical size as well as average run lengths. The results show that the proposed test statistic can make the alarms reliable signals of genuine change, and the price is an increased average run length in detection. A real data example is also provided to examine the efficiency of the proposed procedure.

The rest of the article is organized as follows. Section 2 describes the multinomial logistic regression model. Section 3 contains the proposed procedure, necessary assumptions and their asymptotics. In Section 4, Monte Carlo simulations and real data analysis are conducted to investigate the performance of the proposed procedure. Section 5 concludes the article. All proofs of the theorems are gathered in Appendix.

2 Multinomial logistic regression model

Consider a categorical time series Y t with l categories, where Y t = ( Y t 1 , … , Y t q ) ′ , q = l − 1 ,

Y t j = 1, if the j th category is observed at time t , 0, otherwise

for t = 1 , 2 , … and j = 1 , … , q . Denote the vector of conditional probabilities given ℱ t − 1 by π t = ( π t 1 , … , π t q ) ′ , where

π t j = E [ Y t j ℱ t − 1 ] = P ( Y t j = 1 ℱ t − 1 ) , j = 1 , … , q ,

for every t , σ field ℱ t − 1 is generated by previous observations and covariates Y t − 1 , Y t − 2 , … , Z t − 1 , Z t − 2 , … , that is,

ℱ t − 1 = σ { Y t − 1 , Y t − 2 , … , Z t − 1 , Z t − 2 , … } ,

with p × q matrices { Z t − 1 } , representing the covariate process. Define

Y t l = 1 − ∑ j = 1 q Y t j , π t l = 1 − ∑ j = 1 q π t j .

Following the definition of generalized linear models, the vector of conditional probabilities is linked to the covariate process through the following equation:

π t ( β ) = π t 1 ( β ) π t 2 ( β ) ⋮ π t q ( β ) = h Z t − 1 ′ β

with β , a p-dimensional vector of parameters. The inverse link function h is defined on R q and takes values in R q as well.

In this article, we investigate the multinomial logistic regression model, which is frequently employed in the analysis of nominal time series (Agresti [24], Section 9.2). Let β = ( β 1 ′ , … , β q ′ ) ′ , β j ∈ R d , j = 1 , … , q , then we have

(1) π t j ( β ) = exp β j ′ z t − 1 1 + ∑ i = 1 q exp β i ′ z t − 1 , j = 1 , … , q ,

where z t − 1 is the corresponding d-dimensional vector of stochastic time-dependent covariates independent of j. Obviously,

π t l ( β ) = 1 1 + ∑ i = 1 q exp β i ′ z t − 1 .

3 Sequential change-point detection

In this section, we propose a sequential change-point detection procedure for the coefficients of multinomial logistic regression model. Following the general paradigm of Chu et al. [25], sequential change-point detection uses the initial time period of length m to estimate a model, and its goal is to verify that the probability of the detection approaches α under the null of no change and one under the alternative of a change in parameters after the initial time period.

We first assume that there is no change in the regression parameter during the first m observations, i.e.,

β = β 0 , 1 ≤ t ≤ m .

We are interested in testing the hypothesis

H 0 : β = β 0 , for all π t ( β ) , t = m + 1 , m + 2 , … , κ m , H A : β = β 0 , for all π t ( β ) , t = m + 1 , m + 2 , … , m + k ⁎ , and β ≠ β 0 for π t , m + k ∗ < t ≤ κ m ,

where β = ( β 1 ′ , … , β q ′ ) ′ , β 0 is the true value of the parameter vector β and unknown, k ∗ is the unknown time when change-point occurs in some of the regression coefficients, and κ is some fixed positive integer larger than 1. The sequential change-point detection procedure is based on a detecting statistic Γ m + k and a boundary g ( m , k ) , we stop and reject H 0 at

τ ( m ) = inf { 1 < k < N : Γ m + k ≥ g ( m , k ) } , ∞ , if Γ m + k < g ( m , k ) ,

where N = ( κ − 1 ) m , this is also called the “closed-end” procedure [26]. The detector and boundary must be chosen such that

(2) lim m → ∞ P { τ ( m ) < ∞ } = α , under H 0 ,

(3) lim m → ∞ P { τ ( m ) < ∞ } = 1 , under H A .

Condition (2) ensures that the probability of a false alarm is asymptotically bounded by α , while condition (3) means that a change-point is detected with probability approaching 1.

Next, we will state the assumptions on the covariate vector z t and parameter vector β . We denote the i th component of a vector X as X ( i ) .

Assumption 1

The process { z t } is ergodic and stationary in the sense that for all t and all l ≥ 0 , ( z t , z t + 1 , … , z l ) has the same distribution as ( z 0 , z 1 , … , z l − t ) .

Assumption 2

All components of the covariate vector have finite fourth moment, that is, E | z ( i ) | 4 < ∞ , n → ∞ .

Assumption 3

For the true value of the parameter vector β 0 , we have β 0 ∈ Ω , where Ω ⊂ R p is an open subset.

From Assumptions 1–3, we can obtain the stationarity of various powers of z t ( i ) (see Fokianos et al. [16]),

1 m ∑ t =1 m z t ( i ) z t ( j ) ⟶ a . s . E ( z t ( i ) z t ( j ) ) , n → ∞ , 1 m ∑ t = 1 m z t ( i ) z t ( j ) z t ( l ) ⟶ a . s . E ( z t ( i ) z t ( j ) z t ( l ) ) , n → ∞ ,

for all i , j , l ∈ 1 , 2 , … , d , and

1 m ∑ t =1 m ( Z t − 1 ( Y t − π t ( β ) ) ) ( i ) ( Z t − 1 ( Y t − π t ( β ) ) ) ( j ) ⟶ a . s . E ( Z t − 1 Σ ( β ) Z t − 1 T ) ) ( i , j )

for all i , j ∈ 1 , 2 , … , q d , where

Σ ( β ) = ( σ i j ) q × q , σ i j ( β ) = − π t i ( β ) π t j ( β ) , i ≠ j , π t i ( β ) ( 1 − π t i ( β ) ) , i = j .

The estimation of the parameter vector β in model (1) follows the partial likelihood methodology described in Kedem and Fokianos [27].

The partial likelihood function is

PL ( β ) = ∏ t =1 m ∏ j =1 l π t j Y t j ( β ) ,

so that the partial log-likelihood function is given by

l ( β ) = log PL ( β ) = ∑ t = 1 m ∑ j = 1 l Y t j log π t j ( β ) ,

the partial score process is defined through the partial sum

(4) S m ( β ) = ∇ l ( β ) = ∑ t = 1 m Z t − 1 ( Y t − π t ( β ) ) .

The maximum partial likelihood estimator denoted by β ˆ m is a solution of the score equations S m ( β ) = 0 , and its asymptotic properties have been studied by Kedem and Fokianos [27].

The following results show that the score vector behaves, approximately, as a vector of the Wiener process.

Lemma 3.1

Under Assumptions 1–3, there exists a Wiener process ( W ( t ) ) t ≥ 0 with covariate matrix T such that if β is the true vector of coefficients, then the score vector admits the following approximation:

∥ S m ( β ) − W ( m ) ∥ = O ( m 1/2 − λ ) a . s .,

for some constant λ > 0 , S m ( β ) is defined by (4), and T is the limit of the sample information matrix

(5) 1 m S m ( β ) S m ′ ( β ) ⟶ a .s . T .

Let β ˆ m be the maximum partial likelihood estimator of the true value of the parameter β .

Lemma 3.2

Under Assumptions 1–3,

β ˆ m ⟶ P β , m → ∞ ,

and

m 1/2 T 1/2 ( β ˆ m − β ) ⟶ d N ( 0 , I ) , m → ∞ ,

where I = I p × p is the identity matrix, p = q d , T is given by (5).

Next, we consider the asymptotics under the null hypothesis and alternative hypothesis.

Theorem 3.1

Suppose that Assumptions 1–3 and H 0 hold,

If g ( s ) = c g 1 ( s ) , s ∈ ( 1 , ∞ ) , where g 1 is a given continuous real-valued function with inf 1 < s < ∞ g 1 ( s ) > 0 , then
lim m → ∞ P { τ ( m ) < ∞ } = lim m → ∞ P max 1 ≤ k ≤ N Γ ˆ m , k ∞ ≥ g m + k m = lim m → ∞ P max 1 ≤ k ≤ N ( κ m ) − 1 / 2 T ˆ − 1 / 2 S m + k ( β ˆ m ) ∞ ≥ g m + k m = P sup 1 < s < κ ∥ κ − 1 / 2 ( W ( s ) − s W ( 1 ) ) ∥ ∞ g 1 ( s ) ≥ c .
If g is a constant, g ( s ) ≡ c , then

lim m → ∞ P { τ ( m ) < ∞ } = lim m → ∞ P max 1 ≤ k ≤ N Γ ˆ ( m , k ) ∞ ≥ c = lim m → ∞ P max 1 ≤ k ≤ N ( κ m ) − 1 / 2 T ˆ − 1 / 2 S m + k ( β ˆ m ) ∞ ≥ c = P sup 1 < s < κ κ − 1 / 2 ( W ( s ) − s W ( 1 ) ) ∞ ≥ c ,

where N = ( κ − 1 ) m , κ is some fixed positive integer larger than 1,

Γ ˆ ( m , k ) = ( κ m ) − 1/2 T ˆ − 1/2 S m + k ( β ˆ m ) , T ˆ = 1 m S m ( β ˆ m ) S ′ m ( β ˆ m ) ,

S m + k ( β ˆ m ) = ∑ s = 1 m + k Z s − 1 ( Y s − π s ( β ˆ m ) )

c = c ( α ) is the ( 1 − α ) quantile point of asymptotic distribution, α is the significance level, W ( s ) is a p = q d -dimensional vector of the independent Wiener process and ∥ ⋅ ∥ ∞ denotes the maximum norm.

Theorem 3.2

Suppose Assumptions 1–3 and H A hold, if the coefficient changes from β 0 to β 0 ∗ at m + k ∗ , β 0 ∗ ( j ) = β 0 ( j ) + δ , β 0 ( j ) is the jth component of β 0 , j ∈ { 1 , 2 , … , q d } , where δ is a constant, δ ≠ 0 , then we have

max 1 < k < N ( κ m ) − 1 / 2 T ˆ − 1 / 2 S m + k ( β ˆ m ) ∞ ⟶ P ∞

as m → ∞ , where N , T ˆ − 1/2 and S m + k ( β ˆ m ) are defined in Theorem 3.1.

Proofs of Theorems 3.1 and 3.2 are postponed to Appendix.

4 Simulations and real data analysis

4.1 Simulations

In this section, we considered a multinomial logistic regression model with l = 3 categories and length κ m , where the history data set m = 100, 200, 500 and a monitoring length N = ( κ − 1 ) m with κ = 3 , 5 , 7 . The data were generated according to the model

log π t 1 π t 3 = β 1 ′ z t − 1 = β 10 + β 11 cos ( 2 π t / 12 ) + β 12 Y ( t − 1 ) 1 , log π t 2 π t 3 = β 2 ′ z t − 1 = β 20 + β 21 cos ( 2 π t / 12 ) + β 22 Y ( t − 1 ) 1 ,

with β 1 = ( β 10 , β 11 , β 12 ) = ( 0.3 , 1.25 , 0.5 ) ′ , β 2 = ( β 20 , β 21 , β 22 ) = ( − 0.2, − 2, − 0.75 ) ′ and z t − 1 = ( 1, cos ( 2 π t /12 ) , Y ( t − 1 ) 1 ) . The null hypothesis is rejected at level α at the first time k such that Γ ˆ ( m , k ) > g m + k m , let g 1 ( s ) = s , critical values c are determined by

P sup 1 < s < κ ∥ κ − 1/2 ( W ( s ) − s W ( 1 ) ) ∥ ∞ s ≥ c = α

and

P sup 1 ≤ s ≤ κ ∥ W ( s ) − s W ( 1 ) ∥ ∞ ≥ c = α ,

which could be obtained by Monte Carlo simulation.

Theorems 3.1 and 3.2 are for the case when all parameters are of interest. In some applications this is not the case. It is possible that some parameters are of special concern, and the others are nuisance parameters. In this simulation, we focused on β 10 and β 20 , which means p = 2 . Table 1 reports the simulated empirical sizes for the proposed procedure based on detecting statistic Γ ˆ ( m , k ) , and the significance level is chosen as 5%. It can be seen that the results are close to the significance level 5% in most cases.

Table 1

Empirical sizes of tests based on Γ ˆ ( m , k ) for the 5% significance level

m	g ( s ) ≡ c			g ( s ) = c s
m	N = 2 m	N = 4 m	N = 6 m	N = 2 m	N = 4 m	N = 6 m
100	0.054	0.053	0.054	0.048	0.042	0.056
200	0.062	0.033	0.036	0.042	0.06	0.056
500	0.072	0.038	0.058	0.044	0.054	0.054

Next, we considered the following alternative hypotheses:

H A : β 10 changes from 0.3 to 1.5 and β 20 changes from − 0.2 to − 1.5 at m + k ∗ ,

where k ∗ = 0.01 N , 0.05 N , 0.1 N , 0.2 N , 0.3 N . H A means that there exists a change-point in the regression coefficients. To evaluate the performance of the proposed score test, we made a comparison with the CUSUM detection method (Höhle [21]). The threshold of the CUSUM detection method could be obtained by simulation of the run length. In the control state, suppose that the average run length is set to 370, and then the threshold is selected as 1.7. Monte Carlo simulations were performed in R and all simulations were based on 2,000 replications.

Tables 2–5 show the conditional power, the first and third quartiles ( Q 1 , Q 3 ) , the median Q 2 , the mean (ARL), the maximum of the distribution of the run lengths and the probability of false alarms ( P τ ) . The conditional power is defined as the proportion of signals occurring after m + k ∗ but before the truncation point κ m . The probability of false alarms is obtained as the proportion of signals occurring before the change point m + k ∗ , denoted by P τ . Let k ¯ = m + k ∗ , ARL ≥ k ¯ is the average run length for runs stopping on or after k ¯ .

Table 2

Conditional powers, average of run lengths and probability of false alarms for score test with g ( s ) ≡ c under H A when m = 100 , N = 400

k ∗	Power	Q 1	Q 2	ARL	Q 3	Max	ARL ≥ k ¯	P τ
0.01 N	0.922	219	263	271.25	321	396	271.25	0.000
0.05 N	0.891	218	269	268.93	320	380	268.93	0.000
0.10 N	0.855	214	265	265.98	322	360	265.98	0.000
0.20 N	0.694	212	268	257.71	320	320	257.71	0.000
0.30 N	0.569	206	264	237.63	280	280	239.23	0.001

Table 3

Conditional powers, average of run lengths and probability of false alarms for score test with g ( s ) = c s under H A when m = 100 , N = 400

k ∗	Power	Q 1	Q 2	ARL	Q 3	Max	ARL ≥ k ¯	P τ
0.01 N	0.858	191	250.5	258.65	337	399	258.65	0.000
0.05 N	0.808	204	276	277.5	359	395	277.5	0.000
0.10 N	0.756	211	283	285.28	377	390	285.28	0.000
0.20 N	0.598	270	344	317.76	380	380	318.83	0.004
0.30 N	0.440	300	370	330.75	370	370	330.24	0.002

Table 4

Conditional powers, average of run lengths and probability of false alarms for score test with g ( s ) ≡ c under H A when m = 200 , N = 800

k ∗	Power	Q 1	Q 2	ARL	Q 3	Max	ARL ≥ k ¯
0.01 N	1.000	306	351	358.33	400	735	358.33
0.05 N	1.000	304	349.5	357.6	409	728	357.6
0.10 N	1.000	300	350	356.04	408	668	356.04
0.20 N	0.993	294	352.5	360.56	424	640	360.56
0.30 N	0.958	283	348	356.70	424	560	356.7

Table 5

Conditional powers, average of run lengths and probability of false alarms for score test with g ( s ) = c s under H A when m = 200 , N = 800

k ∗	Power	Q 1	Q 2	ARL	Q 3	Max	ARL ≥ k ¯	P τ
0.01 N	1.000	177	214	225.25	259	617	225.25	0.000
0.05 N	1.000	206	250.5	259.74	301	628	259.74	0.000
0.10 N	1.000	242	291	301.05	348	605	301.05	0.000
0.20 N	0.996	309	386	382.35	440	760	382.35	0.000
0.30 N	0.976	385	451	463.37	525	740	464.63	0.004

From Tables 2–5, it can be seen that the conditional power of score test increases as the sample size increases and decreases as the change-point occurs late. The general relation between power and change-point location in hypothesis testing is that a change-point occurs after a longer period of monitoring and tends to have a lower power. For the boundaries g ( s ) ≡ c and g ( s ) = c s , when m = 100 , the conditional power of the former is higher than that of the latter, and when m = 200 , the conditional power of two boundaries is almost 1.

Next, we investigate the numerical characteristics of run lengths, Q 1 , Q 2 , ARL, Q 3 and maximum of the distribution of run lengths increase with the increase in sample size. When the boundary g ( s ) ≡ c , these numerical characteristics are slightly influenced by the change-point location. For instance, when m = 100 , 200 , the ARLs are around 265 and 357, respectively. However, when the boundary g ( s ) = c s , these numerical characteristics increase as the change-point occurs late. For instance, when m = 100 , k ∗ = 0.01 N , 0.3 N , the ARLs are 258.65 and 330.75, respectively. In addition, compared with the boundary g ( s ) ≡ c , the ARLs of test statistic with g ( s ) = c s are shorter when the change-point is at the beginning of the series being monitored, but they are longer when the change-point occurs late. Finally, the probability of false alarm P τ in two cases is almost zero in most cases.

Table 6 illustrates the conditional power, Q 1 , Q 2 , ARL, Q 3 , maximum of the distribution of run lengths and P τ for the CUSUM detection method. It can be seen from Table 6 that the conditional power and ARLs decrease as change-point occurs late. When the change-point is at the beginning of the series being monitored, the condition power is high, the ARLs are short, and the probability of false alarms is low. However, when the change-point occurs after a period of monitoring, the conditional power decreases rapidly, and the probability of false alarms grows rapidly. The monitoring may stop before the change-point m + k ⁎ , thus the ARLs are negative when change-point occurs late. In comparison with the proposed score test, the ARLs are shorter, but the condition power is lower, and the probability of false alarms is higher.

Table 6

Conditional powers, average of run lengths and probability of false alarms for CUSUM detection method under H A when m = 100 , N = 400

k ∗	Power	Q 1	Q 2	ARL	Q 3	Max	ARL ≥ k ¯	P τ
0.01 N	0.886	47.8	118	153.64	244	396	153.79	0.01
0.05 N	0.651	−3	59	109.21	184	389	150.71	0.261
0.10 N	0.476	−17.2	2	56.22	94.2	360	125.96	0.481
0.20 N	0.276	−56	−32	−1.215	16	320	103.81	0.704
0.30 N	0.194	−95	−68	−42.91	−21	280	77.3	0.795

To summarize, the score test has good control over its type I errors (empirical size), the price is longer ARL. The CUSUM detection method has shorter ARL, but the probability of false alarms grows rapidly when the change-point occurs late.

4.2 Application to real data

In this section, we applied the proposed procedure to a DNA sequence, which was obtained from the gene BNRF1 of the Epstein-Barr Virus (Shumway and Stoffer [28]), and contains n = 1 , 197 observations. Let A, G, C and T denote adenine, guanine, cytosine and thymine nucleotides, respectively, then a sequence of letters A, G, C and T can be viewed as a nominal categorical time series. Fokianos and Kedem (2002) [27] argued in detail that the most appropriate model for these data is the multinomial logistic model with covariate vector z t − 1 = ( 1 , Y ( t − 1 ) 1 , Y ( t − 1 ) 2 , Y ( t − 1 ) 3 , Y ( t − 3 ) 1 , Y ( t − 3 ) 2 , Y ( t − 3 ) 3 ) ′ and link function

log i t ( π t 1 ) = log π t 1 1 − π t 1 − π t 2 − π t 3 = β 10 + β 11 Y ( t − 1 ) 1 + β 12 Y ( t − 1 ) 2 + β 13 Y ( t − 1 ) 3 + β 14 Y ( t − 3 ) 1 + β 15 Y ( t − 3 ) 2 + β 16 Y ( t − 3 ) 3 ,

log i t ( π t 2 ) = log π t 2 1 − π t 1 − π t 2 − π t 3 = β 20 + β 21 Y ( t − 1 ) 1 + β 22 Y ( t − 1 ) 2 + β 23 Y ( t − 1 ) 3 + β 24 Y ( t − 3 ) 1 + β 25 Y ( t − 3 ) 2 + β 26 Y ( t − 3 ) 3 ,

log i t ( π t 3 ) = log π t 3 1 − π t 1 − π t 2 − π t 3 = β 30 + β 31 Y ( t − 1 ) 1 + β 32 Y ( t − 1 ) 2 + β 33 Y ( t − 1 ) 3 + β 34 Y ( t − 3 ) 1 + β 35 Y ( t − 3 ) 2 + β 36 Y ( t − 3 ) 3 ,

where β has 21 components. Let the first 400 samples be the historical sample, then monitor the remaining sample one by one (401–1,197). The proposed score test statistic was computed and graphically reported in Figure 1. The results show that there exists a structural change, which occurs at 738 for β 34 .

$Figure 1 Regression coefficient β 34 {\beta }_{34} is monitored at α = 0.1 \alpha =0.1 level of significance.$

Figure 1

Regression coefficient β 34 is monitored at α = 0.1 level of significance.

To evaluate the validity of the above results, we analyzed the data by the retrospective change point detection procedure (Gombay et al. [17]), and the results show that there exists a structural change at 701. It is also demonstrated that there exists a certain delay, but the proposed score test is effective.

5 Concluding remarks

In this article, we propose a sequential change-point detection procedure based on the partial likelihood score process to detect the coefficients of multinomial logistic regression model. The asymptotic properties of score test statistic are derived under both the null of no change and the alternative of changes in coefficients. To evaluate the finite sample performance of the proposed score test statistic, we conduct a Monte Carlo simulation; simulation results show that the score test method can signal a genuine change reliably, and the price is an increased average run length in detection.

In Monte Carlo simulation, we find that the proposed score test is suitable for the case that the amount of structural change is relatively large and is not sensitive for small deviations. Thus in our simulation, we assume that β 10 changes from 0.3 to 1.5 and β 20 changes from −0.2 to −1.5 under the alternative hypotheses. When the amount of structural change is small, the conditional power is lower, and the ARL is longer, so how to establish a more effective test statistic is the next research subject.

Acknowledgement(s)

The authors are grateful to the editors and referees for useful suggestions and helpful comments for improving the article. This work was supported by the National Nature Science Foundation of China (Program No. 11801438), the Natural Science Basic Research Plan in Shaanxi Province of China (Program No. 2018JQ1089), Scientific Research Program Funded by Shaanxi Provincial Education Department (Program No.18JK0552) and the Innovation Capability Support Program of Shaanxi (Program No. 2020PT-023).

Appendix

Proofs of Lemmas 3.1 and 3.2 are similar to Propositions 1 and 2 of Gombay et al. [17].

Proof of Theorem 3.1

The proof can be done along the lines of retrospective change point detection as in Gombay et al. [17].

First, since S m ( β ˆ m ) = 0 , we can write

( κ m ) − 1/2 S k β ˆ m = ( κ m ) − 1/2 S k ( β ˆ m ) − k m S m β ˆ m .

For the ith component of the score vector, S k ( i ) , i = 1, … , q d , with β 0 denoting the true value of the parameter, we have

( κ m ) − 1/2 S k ( i ) β ˆ m = ( κ m ) − 1/2 S k ( i ) β ˆ m − k m S n ( i ) β ˆ m = ( κ m ) − 1/2 S k ( i ) ( β 0 ) − k m S m ( i ) ( β 0 ) + ( κ m ) 1/2 ∑ j =1 q d β ˆ m ( j ) − β 0 ( j ) 1 κ m ∑ t =1 k ( Z t − 1 ( Y t − π t ( β 0 ) ) ) ( i ) ( Z t − 1 ( Y t − π t ( β 0 ) ) ) ( j ) − k m ∑ t =1 m ( Z t − 1 ( Y t − π t ( β 0 ) ) ) ( i ) ( Z t − 1 ( Y t − π t ( β 0 ) ) ) ( j ) + E k , m ( i ) ,

where E k , m ( i ) have two orders of products of β ˆ m ( j ) − β 0 ( j ) .

Next, we show that

(6) max m + 1 ≤ k ≤ κ m 1 κ m ∑ t =1 k ( Z t − 1 ( Y t − π t ( β 0 ) ) ) ( i ) ( Z t − 1 ( Y t − π t ( β 0 ) ) ) ( j ) − k m ∑ t =1 m ( Z t − 1 ( Y t − π t ( β 0 ) ) ) ( i ) ( Z t − 1 ( Y t − π t ( β 0 ) ) ) ( j ) = o P ( 1 ) .

By Assumptions 1 and 2, we have as k → ∞ ,

1 k ∑ t =1 k ( ( Z t − 1 ( Y t − π t ( β 0 ) ) ) ( i ) Z t − 1 ( Y t − π t ( β 0 ) ) ) ( j ) − E ( Z t − 1 Σ ( β 0 ) Z t − 1 T ) ( i , j ) ) → a . s . 0,

then for every ε > 0 ,

lim k → ∞ P 1 k ∑ t =1 k ( ( Z t − 1 ( Y t − π t ( β 0 ) ) ) ( i ) ( Z t − 1 ( Y t − π t ( β 0 ) ) ) ( j ) ) − E ( Z t − 1 Σ ( β 0 ) Z t − 1 T ) ( i , j ) ) > ε = 0,

hence, as m → ∞ ,

max 1 ≤ k ≤ κ m 1 k ∑ t = 1 k ( ( Z t − 1 ( Y t − π t ( β 0 ) ) ) ( i ) ( Z t − 1 ( Y t − π t ( β 0 ) ) ) ( j ) − E ( Z t − 1 Σ ( β 0 ) Z t − 1 T ) ( i , j ) ) = o P ( 1 ) .

Based on this, we have as m → ∞ ,

max m + 1 ≤ k ≤ κ m / log ( κ m ) k κ m 1 k ∑ t =1 k Z t − 1 ( Y t − π t ( β 0 ) ) ( i ) Z t − 1 ( Y t − π t ( β 0 ) ) ( j ) − E Z t − 1 Σ ( β 0 ) Z t − 1 T ( i , j ) = o P ( 1 ) ,

and

max κ m / log ( κ m ) ≤ k ≤ κ m k κ m 1 k ∑ t =1 k ( ( Z t − 1 ( Y t − π t ( β 0 ) ) ) ( i ) ( Z t − 1 ( Y t − π t ( β 0 ) ) ) ( j ) − E ( Z t − 1 Σ ( β 0 ) Z t − 1 T ) ( i , j ) ) = o P ( 1 ) .

Hence, we get that

max m + 1 ≤ k ≤ κ m 1 κ m ∑ t =1 k ( ( Z t − 1 ( Y t − π t ( β 0 ) ) ) ( i ) ( Z t − 1 ( Y t − π t ( β 0 ) ) ) ( j ) − E ( Z t − 1 Σ ( β 0 ) Z t − 1 T ) ( i , j ) ) = o P ( 1 ) .

Following the proof of Proposition 3 (see Fokianos et al. [16])

max 1 ≤ k ≤ m k m 1 k ∑ t =1 m ( ( Z t − 1 ( Y t − π t ( β 0 ) ) ) ( i ) ( Z t − 1 ( Y t − π t ( β 0 ) ) ) ( j ) − E ( Z t − 1 Σ ( β 0 ) Z t − 1 T ) ( i , j ) ) = o P ( 1 ) ,

then we have

max m + 1 ≤ k ≤ κ m 1 κ m k m ∑ t =1 m ( Z t − 1 ( Y t − π t ( β 0 ) ) ) ( i ) ( Z t − 1 ( Y t − π t ( β 0 ) ) ) ( j ) − E ( Z t − 1 Σ ( β 0 ) Z t − 1 T ) ( i , j ) = max m + 1 ≤ k ≤ κ m m κ m k m 1 m ∑ t =1 m ( Z t − 1 ( Y t − π t ( β 0 ) ) ) ( i ) ( Z t − 1 ( Y t − π t ( β 0 ) ) ) ( j ) − E ( Z t − 1 Σ ( β 0 ) Z t − 1 T ) ( i , j ) = o P ( 1 ) .

Combining the above with m 1/2 β ˆ m ( i ) − β ˆ 0 ( i ) = o P ( 1 ) from Lemma 3.2, we get (6). The error terms E k , m ( i ) have higher orders of products of β ˆ m ( j ) − β 0 ( j ) , and it can be shown that they are o P ( 1 ) by Lemma 3.2.

From Lemmas 3.1 and 3.2, we get as m → ∞

max m + 1 ≤ k ≤ κ m ( κ m ) − 1/2 T − 1/2 S k ( β 0 ) − k m S m ( β 0 ) ⟶ D sup 1 < s < κ κ − 1/2 ( W ( s ) − s W ( 1 ) ) .

Finally, we replace T with T ˆ , and the results hold true by using Slutsky’s theorem.□

Proof of Theorem 3.2

Under the alternative hypothesis, β = β 0 , for all π t ( β ) , t = m + 1 , m + 2 , … , m + k ∗ , and β ≠ β 0 for π t , m + k ∗ < t ≤ κ m . Suppose the coefficient changes from β 0 to β 0 ⁎ at m + k ⁎ , β 0 ∗ ( j ) = β 0 ( j ) + δ , β 0 ( j ) is the jth component of β 0 , j = 1 , … , q d , where δ is a constant, δ ≠ 0 .□

For m + k ∗ < k ≤ κ m , we can write

max m + 1 ≤ k ≤ κ m ( κ m ) − 1 / 2 S k ( β ˆ m ) ≥ ( κ m ) − 1/2 S k ( β ˆ m ) = ( κ m ) − 1/2 ( S k 1 ( β ˆ m ) + S k 2 ( β ˆ m ) ) ,

where

S k 1 ( β ˆ m ) = ∑ t = 1 m + k ∗ − 1 Z t − 1 ( Y t − π t ( β ˆ m ) ) , S k 2 ( β ˆ m ) = ∑ t = m + k ∗ k Z t − 1 ( Y t − π t ( β ˆ m ) ) .

For the ith component of S k 1 ( β ˆ m ) , i = 1 , … , q d , we have

( κ m ) − 1/2 S k 1 ( i ) β ˆ m = ( κ m ) − 1/2 S k 1 ( i ) β 0 + ( κ m ) 1/2 ∑ j =1 q d β ˆ m ( j ) − β 0 ( j ) ( κ m ) − 1 ∑ t =1 m + k ∗ − 1 ( Z t − 1 ( Y t − π t ( β 0 ) ) ) ( i ) ( Z t − 1 ( Y t − π t ( β 0 ) ) ) ( j ) + E k 1 , m ( i ) , ( κ m ) − 1/2 S k 2 ( i ) β ˆ m = ( κ m ) − 1/2 S k 2 ( i ) β 0 ⁎ + ( κ m ) 1/2 ∑ j =1 q d β ˆ m ( j ) − β 0 ⁎ ( j ) ( κ m ) − 1 ∑ t = m + k ⁎ k Z t − 1 Y t − π t β 0 ⁎ ( i ) Z t − 1 Y t − π t β 0 ⁎ ( j ) + E k 2 , m ( i ) ,

where E k 1 , m ( i ) has two orders of products of β ˆ m ( j ) − β 0 ( j ) , E k 2 , m ( i ) has two orders of products of β ˆ m ( j ) − β 0 ⁎ ( j ) .

As m → ∞ , by Theorem 3.1 we have

(7) ( κ m ) − 1/2 S k 1 ( i ) β ˆ m = o P ( 1 ) ,

(8) ( κ m ) − 1/2 S k 2 ( i ) β 0 ⁎ = o P ( 1 ) ,

following Assumptions 1–3, we conclude

(9) ( κ m ) − 1 ∑ t = m + k ∗ k Z t − 1 ( Y t − π t β 0 ∗ ( i ) Z t − 1 ( Y t − π t β 0 ∗ ( j ) = o P ( 1 ) .

Since δ ≠ 0 and δ is a constant, as m → ∞ we have

(10) ( κ m ) 1/2 β ˆ m ( j ) − β 0 ⁎ ( j ) = ( κ m ) 1/2 β ˆ m ( j ) − β 0 ( j ) − δ → P ∞

by Lemma 3.2. Putting together (7), (8), (9) and (10), the proof of Theorem 3.4 is complete.

Conflict of Interest: The authors report no potential conflict of interest.

References

[1] Miklós Csörgö and Lajos Horváth, Limit Theorems in Change-Point Analysis, John Wiley & Sons Inc., New York, 1997.Search in Google Scholar

[2] Pierre Perron, Dealing with structural breaks, Palgrave Handbook of Econometrics 1 (2006), 278–352, 10.1016/j.gfj.2006.04.004.Search in Google Scholar

[3] Edit Gombay, Change detection in autoregressive time series, J. Multivariate Anal. 99 (2008), no. 3, 451–464, 10.1016/j.jmva.2007.01.003.Search in Google Scholar

[4] Zhanshou Chen and Zheng Tian, Modified procedures for change point monitoring in linear models, Math. Comput. Simulat. 81 (2010), no. 1, 62–75, 10.1016/j.matcom.2010.06.021.Search in Google Scholar

[5] Okyoung Na, Youngmi Lee, and Sangyeol Lee, Monitoring parameter change in time series models, Stat. Methods Appl. 20 (2011), no. 2, 171–199, 10.1007/s10260-011-0162-3.Search in Google Scholar

[6] Changliang Zou, Guosheng Yin, Long Feng, and Zhaojun Wang, Nonparametric maximum likelihood approach to multiple change-point problems, Ann. Stat. 42 (2014), no. 3, 970–1002, 10.1214/14-aos1210.Search in Google Scholar

[7] J. Ross Gordon, Parametric and nonparametric sequential change detection in R: The cpm package, J. Stat. Softw. 66 (2015), no. 3, 1–20, 10.18637/jss.v066.i03.Search in Google Scholar

[8] Michael W. Robbins, Colin M. Gallagher, and Robert B. Lund, A general regression change point test for time series data, J. Am. Stat. Assoc. 111 (2016), no. 514, 670–683, 10.1080/01621459.2015.1029130.Search in Google Scholar

[9] Fuxiao Li, Zheng Tian, and Zhanshou Chen, Monitoring distributional changes of squared residuals in GARCH models, Comm. Statist. Theory Methods 46 (2017), no. 1, 354–372, 10.1080/03610926.2014.995819.Search in Google Scholar

[10] Yang Cao, Liyan Xie, Yao Xie, and Huan Xu, Sequential change-point detection via online convex optimization, Entropy 20 (2018), no. 2, 108, 10.3390/e20020108.Search in Google Scholar PubMed PubMed Central

[11] Hao Chen, Sequential change-point detection based on nearest neighbors, Ann. Stat. 47 (2019), no. 3, 1381–1407, 10.1214/18-aos1718.Search in Google Scholar

[12] Jaromír Antoch, Gérard Gregoire, and Daniela Jarušková, Detection of structural changes in generalized linear models, Stat. Probabil. Lett. 69 (2004), no. 3, 315–332, 10.1016/j.spl.2004.06.028.Search in Google Scholar

[13] Hongling Zhou and Kung-Yee Liang, On estimating the change point in generalized linear models, IMS Collections Beyond Parametrics in Interdisciplinary Research: Festschrift in Honor of Professor Pranab K. Sen 1 (2008), 305–320, 10.1214/193940307000000239.Search in Google Scholar

[14] Konstantinos Fokianos and Benjamin Kedem, Regression theory for categorical time series, Stat. Sci. 18 (2003), no. 3, 357–376, 10.1214/ss/1076102425.Search in Google Scholar

[15] Lajos Horváth, Marie Hušková, Piotr Kokoszka, and Josef Steinebach, Monitoring changes in linear models, J. Stat. Plan. Infer. 126 (2004), no. 1, 225–251, 10.1016/j.jspi.2003.07.014.Search in Google Scholar

[16] Konstantinos Fokianos, Edit Gombay, and Abdulkadir Hussein, Retrospective change detection for binary time series models, J. Stat. Plan. Infer. 145 (2014), 102–112, 10.1016/j.jspi.2013.08.017.Search in Google Scholar

[17] Edit Gombay, Fuxiao Li, and Hao Yu, Retrospective change detection in categorical time series, Comm. Statist. Simulation Comput. 46 (2017), no. 14, 6831–6845, 10.1080/03610926.2015.1137595.Search in Google Scholar

[18] Šárka Hudecová, Structural changes in autoregressive models for binary time series, J. Stat. Plan. Infer. 143 (2013), no. 10, 1744–1752, 10.1016/j.jspi.2013.05.009.Search in Google Scholar

[19] Guanghui Wang, Changliang Zou, and Guosheng Yin, Change-point detection in multinomial data with a large number of categories, Ann. Stat. 46 (2018), no. 5, 2020–2044, 10.1214/17-AOS1610.Search in Google Scholar

[20] Zhiming Xia, Pengjiang Guo, and Wenzhi Zhao, Monitoring structural changes in generalized linear models, Comm. Statist. Theory Methods 38 (2009), no. 11, 1927–1947, 10.1080/03610920802549910.Search in Google Scholar

[21] Michael Höhle, Online change-point detection in categorical time series, in: T. Kneib, G. Tutz (eds), Statistical Modelling and Regression Structures, Physica, Springer-Verlag Berlin Heidelberg, 2010, pp. 377–397, 10.1007/978-3-7908-2413-1_20.Search in Google Scholar

[22] Edit Gombay and Daniel Serban, Monitoring parameter change in AR(p) time series models, J. Multivariate Anal. 100 (2009), no. 4, 715–725, 10.1016/j.jmva.2008.08.005.Search in Google Scholar

[23] Edit Gombay, Abdulkadir A. Hussein, and Stefan H. Steiner, Monitoring binary outcomes using risk-adjusted charts: a comparative study, Stat. Med. 30 (2011), no. 23, 2815–2826, 10.1002/sim.4305.Search in Google Scholar PubMed

[24] Alan Agresti, Categorical Data Analysis, Wiley, New York, 1990.Search in Google Scholar

[25] Chia-Shang James Chu, Maxwell Stinchcombe, and Halbert White, Monitoring structural change, Econometrica 64 (1996), no. 5, 1045–1065, 10.2307/2171955.Search in Google Scholar

[26] Zdenĕk Hlávka, Marie Hušková, Claudia Kirch, and Simos G. Meitanis, Monitoring changes in the error distribution of autoregressive models based on Fourier methods, Test 21 (2012), 605–634, 10.1007/s11749-011-0265-z.Search in Google Scholar

[27] Konstantinos Fokianos and Benjamin Kedem, Regression model for time series analysis, John Wiley & Sons, Hoboken, 2002.10.1002/0471266981Search in Google Scholar

[28] Robert H. Shumway and David S. Stoffer, Time Series Analysis and its Applications, Springer, New York, 2002.Search in Google Scholar

Received: 2019-07-06

Revised: 2020-03-17

Accepted: 2020-05-15

Published Online: 2020-07-29

This work is licensed under the Creative Commons Attribution 4.0 International License.

Articles in the same Issue

https://doi.org/10.1515/math-2020-0037

Keywords for this article

categorical time series; sequential change-point detection; multinomial logistic regression model; partial likelihood score process

Creative Commons

BY 4.0