Bastien Fernandez May

Discontinuous generalized synchronization of chaos

Keywords: 2010 Mathematical Subject Classification. 37A99, 26A45 Forced systems, synchronization graphs, bounded variation

We study synchronization functions in basic examples of discontinuous forced systems with contractive response and chaotic driving. The forcing is given by baker-type maps and the response is assumed to depend monotonically on the drive. The resulting synchronization functions have dense sets of discontinuities and their graphs appear to be extremely choppy. We show that these functions have bounded variation when the contraction is strong, and conversely, that their total variation is infinite when the contraction becomes weak. In the first case, we also analyze in detail smoothness properties of the corresponding continuous component.

Introduction

Analyzing the asymptotic response to a random or erratic stimulus is a ubiquitous problem in Nonlinear Dynamics. The archetypical example is given by synchronization phenomena in directionally coupled systems [START_REF] Pecora | Synchronization in chaotic systems[END_REF]. (For convenience, focus will be on discrete time dynamics throughout this paper). When an autonomous forcing x t+1 = f (x t ) compels the iterations of a dissipative factor z t+1 = g(x t , z t ), the dynamics is known to be attracted by the invariant graph z = h(x) of an associated synchronization function h [START_REF] Afraimovich | Stochastic synchronization of oscillations in dissipative systems[END_REF][START_REF] Kocarev | Generalized synchronization, predictability, and equivalence of unidirectionally coupled dynamical systems[END_REF][START_REF] Rulkov | Generalized synchronization of chaos in directionally coupled chaotic systems[END_REF]. In brief, the response z is asymptotically locked (viz. conjugated) to the drive x.

In this context, the regularity and smoothness of h -which depend on f and g -prescribe those drive features (such as e.g. Lyapunov exponents, fractal dimensions, etc) that are conveyed to the factor. For instance, Lipschitz regularity implies some estimates on the attractor's dimension. Applications range from modeling of (low-pass) filters in signal analysis [START_REF] Kaplan | The Lyapunov dimension of a nowhere differentiable attracting torus[END_REF][START_REF] Badii | Dimension function and phase transition-like behavior in strange attractors[END_REF][START_REF] Broomhead | Linear filters and nonlinear systems[END_REF][START_REF] Hunt | Differentiable generalized synchronization of chaos[END_REF][START_REF] Carroll | Detecting recursive and nonrecursive filters using chaos[END_REF] to damage detection in material science [START_REF] Nichols | Use of fiber-optic strain sensors and Holder exponents for detecting and localizing damage in an experimental plate structure[END_REF].

The mathematical theory of synchronization functions is part of the broader study of inertial manifolds in dynamical systems. It goes back to the seminal work of Hirsch, Pugh and Shub [START_REF] Hirsch | Stable Manifolds and Hyperbolic Sets in Global Analysis[END_REF][START_REF] Hirsch | Invariant Manifolds[END_REF] who proved existence and continuity under the assumption that f is a homeomorphism and g is contracting for z. They also evaluated the Hölder continuity and showed Lipschitz regularity when the contraction of g is stronger than the stronger contraction rate of f (and additional mild assumptions on f and g). These results were perfected later on [START_REF] Campbell | The existence of inertial functions in skew product systems[END_REF][START_REF] Afraimovich | Synchronization in unidirectionally coupled systems: some rigorous results[END_REF][START_REF] Urias | Filters display inverse limit spaces[END_REF], especially by Stark [START_REF] Stark | Invariant graphs for forced systems[END_REF][START_REF] Stark | Regularity of invariant graphs for forced systems Ergod[END_REF] who established smoothness and extended them to non-uniform contractive responses (see also [START_REF] Singh | Scenarios for generalized synchronization with chaotic driving[END_REF][START_REF] Hu | Holder continuity of three types of generalized synchronization manifolds of non-autonomous systems[END_REF] for recent developments).

Studies of continuity have been pursued beyond the homeomorphic case, not only when the response function g remains blind to drive discontinuities, but also when f is not invertible [START_REF] Rulkov | Detectability of nondifferentiable generalized synchrony[END_REF][START_REF] Barreto | The geometry of chaos synchronization[END_REF]. However, continuity may not always hold in applications [START_REF] Badii | Dimension function and phase transition-like behavior in strange attractors[END_REF][START_REF] Broomhead | Linear filters and nonlinear systems[END_REF][START_REF] Hunt | Differentiable generalized synchronization of chaos[END_REF] and sensitive response functions also have to be considered.

To that aim, this paper considers discontinuous synchronization graphs in basic examples of skew-products (f, g) with chaotic driving f and contractive response g i.e.

|g(x, z) -g(x, z )| ≤ λ|z -z | (1)
where 0 < λ < 1. As noticed in the literature [START_REF] Afraimovich | Synchronization in unidirectionally coupled systems: some rigorous results[END_REF][START_REF] Stark | Invariant graphs for forced systems[END_REF][START_REF] Stark | Regularity of invariant graphs for forced systems Ergod[END_REF], the response h lacks monotonicity and may have infinitely many discontinuity points that generate a dense subset in phase space of the forcing. The appropriate notion to investigate in this case is the overall graph "length" (i.e. total variation of h), together with properties of the continuous and discontinuous components. We proceed to a thorough analytic study in the case where f is an extension of the baker's map and g(x, z) = λz + (1 -λ)x is linear [START_REF] Badii | Dimension function and phase transition-like behavior in strange attractors[END_REF][START_REF] Hunt | Differentiable generalized synchronization of chaos[END_REF][START_REF] Afraimovich | Synchronization in unidirectionally coupled systems: some rigorous results[END_REF]. The total variation is shown to possibly diverge depending on the contraction parameter λ. When this quantity remains finite, an analysis of regularity and of the derivative of the continuous component is given. The specific form of h here allows for results well beyond the standard theory of real functions. The basic properties of the synchronization function however do not depend on the piecewise affine nature of the dynamical system under study. In the last section, we present results that hold for more general discontinuous skew-product systems (f, g).

Piecewise linear skew-products and synchronization graphs

As announced before, the autonomous forcing in this study materializes via the generalized baker's map (x, y) → f (x, y) = (T a,b (x, y), T b (y)) of the unit square [0, 1] 2 into itself [START_REF] Badii | Dimension function and phase transition-like behavior in strange attractors[END_REF][START_REF] Hunt | Differentiable generalized synchronization of chaos[END_REF]. Here 0 < a, b < 1 and the mapping components write

T a,b (x, y) = ax if 0 ≤ y < b (1 -a)x + a if b ≤ y ≤ 1 and T b (y) = y b if 0 ≤ y < b y-b 1-b if b ≤ y ≤ 1
Remark: These notations are generic and will often be used with other parameters/variables throughout the text. For instance, T a (x) denotes the interval map above with parameter a and variable x (instead of b and y respectively).

Due to the contraction 0 < λ < 1 the response system z t → z t+1 = g(x t , z t ) = λz t + (1 -λ)x t results to be asymptotically locked to the forcing term [START_REF] Afraimovich | Stochastic synchronization of oscillations in dissipative systems[END_REF][START_REF] Kocarev | Generalized synchronization, predictability, and equivalence of unidirectionally coupled dynamical systems[END_REF][START_REF] Rulkov | Generalized synchronization of chaos in directionally coupled chaotic systems[END_REF]. More precisely, the large time behavior of the sequence {z t } is independent of z 0 and approaches {h(x t )} where the synchronization function h is given by

h(x) = (1 -λ) ∞ t=0 λ t T t+1 a (x), ∀x ∈ [0, 1]. ( 2
)
This property is a consequence of the equality z t -h(x t ) = λ t (z 0 -h(x 0 )).

The response h only depends on the first forcing variable because f is invertible and the first coordinate of the inverse f -1 (x, y) = (T a (x), T b,a (y, x)) only depends on x. (Besides, as a function of λ, the map h is smooth and strictly increasing and uniformly converges to T a when λ → 0.) Furthermore, this function solves the conjugacy equation

g(x, h(x)) = h • T a,b (x, y), ∀x, y ∈ [0, 1] 2 .
thanks to the property T a • T a,b (x, y) = x. When f is not invertible, h a priori depends on backward histories {x t } t≤0 , see e.g. [START_REF] Rulkov | Detectability of nondifferentiable generalized synchrony[END_REF][START_REF] Barreto | The geometry of chaos synchronization[END_REF]. Exponential convergence of the series (2) guarantees that h is welldefined on [0, 1]. The function can be viewed as a uniform limit of piecewise affine maps obtained by truncating the series to finite order (see equation [START_REF] Rulkov | Generalized synchronization of chaos in directionally coupled chaotic systems[END_REF] in the proof of Proposition 3.2 below). This comment provides a convenient way to numerically compute the graph of h with arbitrary accuracy, see Figure 1.

The map T a is right continuous and piecewise increasing. The same properties hold for all maps T t a . By uniform convergence, right continuity 4)). In the left pictures (λ < 1/2), h is of bounded variation; in the right ones (λ > 1/2), it has infinite variation. Notice the symmetry h

h(x) h(x) h(x) h(x) x x x = 0.4
(x) = 1 -h(1 -x -0) for a = 0.5 (bottom pictures).
transfers to h, viz. we have h(x + 0) = h(x) for all x. In addition, the maps T t a are all piecewise increasing; thus the left limit h(x -0) exists everywhere. This limit coincides with h(x) unless x is a pre-image of the discontinuity point a. More precisely, we have

h(x -0) > h(x) iff x ∈ D a := ∞ t=0 T -t a (a).
Since T a is expanding, D a is a dense subset of [0, 1]. To prove that all jump discontinuities are negative, we start to notice that uniform convergence yields

h(x -0) = (1 -λ) ∞ t=0 λ t T t+1 a (x -0), ∀x ∈ [0, 1].
By definition, for every point x ∈ D a there is a unique t 0 ≥ 0 such that T t 0 a (x) = a and T t a (x) = a when 0 ≤ t < t 0 . The map T a is continuous everywhere but at the point a; hence we have

T t+1 a (x -0) = T t+1 a (x), ∀0 ≤ t < t 0 . Using again T t 0 a (x) = a, we get T t 0 +1 a (x -0) = 1 and T t 0 +1 a (x) = 0.
However, these points 0 and 1 are fixed points; hence the same values hold for subsequent iterates. Altogether we obtain explicit estimates for jump discontinuities at every point of

D a h(x -0) -h(x) = (1 -λ) ∞ t=t 0 λ t > 0.
On the other hand, if x lies outside D a , we have lim y→x T t (y) = T t (x) for all t ≥ 0. Uniform convergence then implies that lim y→x h(y) = h(x) as claimed.

Beside negative jump discontinuities, the function h has positive increments in arbitrary small left-neighborhoods of every point outside D a . More precisely, for every x ∈ D c a and > 0 there exists y ∈ (x -, x) such that h(y) < h(x). Therefore, h is 'nowhere monotonous' (i.e. there is no interval on which h is monotonous) and its graph must be extremely wrinkled and choppy as Figure 1 indicates.

Proof that for every x ∈ D c a , there exists y < x and arbitrarily close such that h(y) < h(x). Given x ∈ D c a and > 0, let

t 0 = min{t ≥ 0 : ∃ n : x -< x n t < x},
where x n t is defined by T t a (x n t ) = a. The number t 0 exists because D a is dense. By definition of t 0 , all iterates of x n t 0 and x up to t 0 -1 must lie on the same side of the discontinuity point. By strict monotonicity of T t a this implies T t a (x n t 0 ) < T t a (x) when 0 ≤ t ≤ t 0 . Moreover, we have T t a (x n t 0 ) = 0 when t > t 0 and T t a (x) > 0 for all t since x ∈ D c a . It easily follows that h(x n t 0 ) < h(y).

Finally, uniform convergence and right continuity allow one to prove that the range of h is an entire interval, the unit interval indeed because of normalization, i.e. we have h([0, 1]) = [0, 1], see Figure 1.

Proof that Ran(h) = [0, 1]. The crucial point is to show that for every n we have h n ([0, 1]) = [0, 1 -λ n+1 ]
where the approximations h n are defined in the relation (4) below.

Consider the cylinder sets [θ 0 • • • θ n ] that are associated with the usual symbolic dynamics of (T a , [0, 1]), i.e. θ t ≡ H(T t a (x) -a) where H is the Heaviside function. Cylinder sets are intervals and their union (with length fixed) covers [0, 1].

The iterates T t+1 a are piecewise bijections and we have The previous property immediately implies that Ran(h) is dense in [0, 1]. Indeed for any strict sub-interval I in this set, there must be at least one point of h([0, 1]). To see this, it suffices to take n such that |h(x) -h n (x)| is uniformly smaller than |I|/2. Assume also that n is sufficiently large so that I is contained in [0, 1 -λ n+1 ]. Choose x such that h n (x) is the middle of I. Then h(x) must belong to I. Now, by continuity we know that every h(x) for x ∈ D c a can be realized as a limit lim y→x h(y). For points in D a , we similarly take right limits. This shows that h([0, 1]) is the right closure of a dense subset in [0, 1], namely it consists of this entire interval excepted the right boundary 1. But h(1) = 1, thus the proof is complete. measures the integral of the modulus of the derivative. As formally claimed in the next statement, the finiteness of the total variation depends on the contraction factor λ. When the variation is finite, we also simultaneously provide the decomposition into the difference h c -h d of two increasing functions and the one into the sum h c + (-h d ) of a continuous function h c and a step function -h d (both decompositions are granted by standard theorems on BV-functions [START_REF] Kolmogorov | Elements of the theory of functions and functional analysis[END_REF]). Some illustrations are given in Figure 2.

T t+1 a ([θ 0 • • • θ n ]) = [θ t+1 • • • θ n ] for 0 ≤ t ≤ n (where [θ n+1 • • • θ n ] should be understood as [0, 1]). This yields h n ([θ 0 • • • θ n ]) = (1 -λ) n t=0 λ t [θ t+1 • • • θ n ] which is an interval. Moreover, given two adjacent cylinders [θ 0 • • • θ n ] ≤ [ θ0 • • • θn ], their images [θ t+1 • • • θ n ]
Theorem 3.1 h is of bounded variation iff λ < 1 2 . Under this condition, the function h d given by

h d (x) = y∈Da : y≤x h(y -0) -h(y), ∀x ∈ [0, 1]
is well-defined (and is increasing) and the function h c defined by h c = h + h d is continuous and strictly increasing on [0, 1] with range [0, 2-2λ 1-2λ ].

Proof. We begin to prove that the variation of h is infinite when λ ≥ 1 2 . Since all discontinuities of h are negative jumps, the total variation is at least x∈Da h(x -0) -h(x). There are 2 n points in D a for which T n a (x) = a and the corresponding difference

h(x -0) -h(x) is equal to (1 -λ) ∞ t=n λ t .
Therefore the variation of h is at least

x∈Da h(x -0) -h(x) = ∞ n=0 (2λ) n . (3)
It easily follows that the total variation is infinite when λ ≥ 1 2 . To continue, we assume that λ < 1 2 . Relation (3) implies that the function h d is well-defined (increasing and bounded) on [0, 1]. Moreover

h d (x-0) = y∈Da : y<x h(y-0)-h(y), hence h d (x)-h d (x-0) = h(x-0)-h(x).
All discontinuities of h are contained in h d and thus the map h c defined by

h c = h + h d is continuous.
We now show that h c is strictly increasing and bounded, and as a consequence, that h is of bounded variation. Assume that x < y and using the same notations as in the proof above that h is nowhere monotonous, let t 0 = min{t ≥ 0 : ∃ n : x < x n t ≤ y}, where x n t is defined by T t a (x n t ) = a. We have

T t a (x) < T t a (x n t 0 ) ≤ T t a (y), 0 ≤ t ≤ t 0 . Hence h(x) < (1 -λ) t 0 -1 t=0 λ t T t+1 a (x n t 0 ) + ∞ t=t 0 λ t = h(x n t 0 -0) and (1 -λ) t 0 -1 t=0 λ t T t+1 a (x n t 0 ) = h(x n t 0 ) ≤ h(y)
(The sum from 0 to t 0 -1 equals 0 if t 0 = 0.) Using also the definition of h d and h c , we obtain

h c (x) = h(x) + h d (x) < h(x n t 0 -0) + h d (x) ≤ h(x n t 0 -0) + h d (x n t 0 -0) = h c (x n t 0 ), and
h c (y) = h(y) + h d (y) ≥ h(x n t 0 ) + h d (y) ≥ h(x n t 0 ) + h d (x n t 0 ) = h c (x n t 0 ), i.e. h c (x) < h c (y).
Finally we have h c (0) = 0 and h c (1) = h(1)+h d (1) < +∞ and the proof is complete.

We now study in detail smoothness properties of the continuous component h c . In a general setting, Kolmogorov's theorem says that monotonicity implies the existence of a finite derivative almost everywhere. In our case, the derivative h c (x) of the continuous component h c turns out to exist everywhere outside the countable set D a , provided that λ is further restricted. Proposition 3.2 If λ < min{a, 1 -a}, then the map h c is differentiable with bounded derivative at any point in the complement D c a . The derivative can be written

h c (x) = (1 -λ) ∞ t=0 λ t (T t+1 a (x)) , ∀x ∈ D c a .
In addition for any λ ≥ min{a, 1 -a}, there are points in D c a where h c (x) diverges.

Proof. The map h can be regarded as the uniform limit of the sequence {h n } where

h n (x) = (1 -λ) n t=0 λ t T t+1 a (x), ∀x ∈ [0, 1]. ( 4
)
Each map h n is right continuous and with negative jump discontinuities in the set

D n a = n t=0
T -t a (a). As expected, the map h n,d defined by

h n,d (x) = y∈D n a : y≤x h n (y -0) -h n (y), ∀x ∈ [0, 1]
is an increasing step function that contains all discontinuities of h n . Consequently, the subsequent map h n,c = h n +h n,d is a piecewise affine continuous map with finitely many affine parts. As such, it is absolutely continuous and hence differentiable almost everywhere in [0, 1] [START_REF] Kolmogorov | Elements of the theory of functions and functional analysis[END_REF]. The derivative h n,c is a summable function. The fundamental theorem of calculus then yields

h n,c (x) -h n,c (0) = x 0 h n,c (y)dy, ∀x ∈ [0, 1] (5)
The map h n,c is actually differentiable on [0, 1] \ D n a . Letting T a (a) = 1 1-a , the derivative h n,c can be uniquely continued to the following step function on [0, 1]

h n,c (x) = (1 -λ) n t=0 λ t (T t+1 a (x))
Since T a (x) ∈ { 1 a , 1 1-a }, when λ < min{a, 1 -a}, the sequence {h n,c } uniformly converges to a bounded map, say h c . Applying Lebesgue's dominated convergence theorem, one can take the limit n → ∞ in (5) to obtain the following equality

h c (x) -h c (0) = x 0 h c (y)dy, ∀x ∈ [0, 1]
Now, the map h c is uniformly approximated by step functions and is continuous at every point of [0, 1] \ D a . A standard result [START_REF]Theorem 6.1[END_REF] states that the map x → x 0 h c (y)dy is differentiable at every point of [0, 1] \ D a and with derivative h c . We conclude from the previous equality that h c is differentiable on [0, 1] \ D a with derivative h c .

The function h c can not be differentiable in D a when a = 1 2 because the right and left derivatives of T a are unequal at x = a. However, by artificially setting T a (a) = 1 1-a for the derivative at the discontinuity point, the domain of h c extends to the entire interval [0, 1] when λ ≤ min{a, 1 -a}.

Alternatively, the map h c can be regarded as the (well-defined) right derivative on [0, 1]. The extended map h c intriguingly shares several properties with the original function h. Assuming a = 1 2 , this map is right continuous with jump discontinuities at every point in D a . (This map is constant when a = 1 2 .) Unlike for h however, the signs of h c jumps do depend on parameters. The analysis reveals that all quantities

h c (• -0) -h c (•) are positive if (a - 1
2 )(λ -a(1 -a)) > 0 and negative otherwise. Moreover we have the following statement analogous to Theorem 3.1.

Proposition 3.3 Assume a = 1 2 .
The function h c is of bounded variation iff λ < a(1 -a). Under this condition, the map (h c ) d defined by

(h c ) d (x) = y∈Da : y≤x h c (y -0) -h c (y)
is well-defined on [0, 1] (and non-increasing or non-decreasing depending on the sign of a -1 2 ). In addition the function defined by (h c ) c = h c + (h c ) d is continuous and strictly increasing if a < 1 2 (resp. strictly decreasing if a > 1 2 ). It is differentiable at any x ∈ D c a and the derivative identically vanishes.

Proof. To begin, each map x → (T t+1 a (x)) is piecewise constant and right continuous; hence h c is also right continuous. As before, the set D a collects all discontinuity points. A similar calculation to that in the proof of [START_REF] Afraimovich | Stochastic synchronization of oscillations in dissipative systems[END_REF] shows that if x ∈ D a is such that T t 0 a (x) = a for some t 0 ≥ 0, then we have

h c (x -0) -h c (x) = C a,λ λ t 0 (T t 0 a (x))
where

C a,λ = (1 -λ) 1 a(1 -λ 1-a ) - 1 (1 -a)(1 -λ a )
.

In particular, the sign of h c (x -0) -h c (x) is independent of x. However, it depends on parameters via the constant C a,λ . The latter is positive iff (a -1 2 )(λ -a(1 -a)) > 0. Now in order to estimate the total variation, as before, we have to sum up all contributions from discontinuities. Given x such that T t 0 a (x) = a for some t 0 ≥ 0, the derivative (T t 0 a (x)) can be written

a -k (1 -a) -(t 0 -k)
where k is the number of those iterates {T t a (x)} t 0 -1 t=0 that are smaller than a. Up to few sequences, the symbolic dynamics of T a is the full shift on two symbols. Thus, for every 0 ≤ k ≤ t 0 , there are t 0 k points with k iterates lying below a. Consequently, the total variation of h c /C a,λ is bounded below by

(1 -λ) ∞ t=0 λ t t k=0 t k 1 a k (1 -a) t-k = 1 -λ λ ∞ t=0 λ a(1 -a) t+1 which is infinite when λ ≥ a(1 -a).
To prove bounded variation when λ < a(1 -a), we again proceed as before. Assume that a < 1 2 and consider the non-decreasing step function (h c ) d associated with h c , viz.

(h c ) d (x) = y∈Da : y≤x h c (y -0) -h c (y), ∀x ∈ [0, 1]
and also the continuous component defined by (h c ) c = h c + (h c ) d . Similar arguments to those in the proof of Theorem 3.1 show that the map (h c ) c is strictly increasing and bounded on [0, 1]. It follows that the variation of h c is finite when λ < a(1 -a). An analogous construction applies when a > 1 2 . As in the proof of Proposition (3.2), the proof that (h c ) c is differentiable proceeds by considering a sequence of uniform approximations of h c by piecewise affine functions with finitely many branches.

Extensions to Nonlinear Systems

The basic results on synchronization functions in piecewise monotonous forced systems actually do not rely on the affine assumption. They extend to more general systems. To present extensions, we start by introducing nonlinear generalized baker's maps with finitely many pieces.

Let N > 1 be an arbitrary integer and consider the finite collections {I i } N i=1 and {J i } N i=1 of intervals defined by

I i = [x i , x i+1 ) where 0 = x 1 < x 2 < • • • < x N +1 = 1 and J i = [y i , y i+1 ) where 0 = y 1 < y 2 < • • • < y N +1 = 1 for i = 1, • • • , N -1, together with I N = [x N , 1] and J N = [y N , 1]
. Consider two mappings T and S on [0, 1] defined by

T | I i ≡ T i and S| J i ≡ S i , i = 1, • • • , N
where each T i : I i → [0, 1) (and T N :

I N → [0, 1]
) is a C 1 increasing, one-toone and onto function and similarly for S i . Now define the map f on [0, 1] 2 by f (x, y) = (T (x, y), S(y)) where

T (x, y) = T -1 i (x) if y ∈ J i
The map f is easily checked to be invertible with inverse given by f -1 (x, y) = (T (x), S(y, x)) where

S(y, x) = S -1 i (y) if x ∈ I i
As in the piecewise affine case, the dynamics of the skew-product (f, g) (where g(x, z) = λz + (1 -λ)x still remains unchanged) is attracted by the graph of a function h defined by equation ( 2) where T a is replaced by T . The new synchronization function h shares several properties with the original one. It is piecewise increasing, right continuous with negative jump discontinuities at every point of the set D defined by

D = ∞ t=0 T -t ({x 2 , • • • , x N }) (6)
and h([0, 1]) = [0, 1]. Moreover, the conclusions of Theorem 3.1 equally repeat in this case provided that 1 2 is replaced by 1 N in the condition on the contraction parameter λ. Indeed, the only novelty sits in the number of t-preimages of discontinuity point {x 2 , • • • , x N } (where h(x-0)-h(x) = λ t ) which is now given by (N -1)N t .

If, in addition the map T is piecewise affine, then analogous to statements of Propositions 3.2 and 3.3 hold for arbitrary N > 2.

Next, we consider the case where the response g is also nonlinear, i.e. we assume -relation (1) holds for some 0 < λ < 1, -the maps g(•, z) are all strictly increasing and the family is equicontinuous.

In such cases, the existence of a globally attracting synchronization function has previously been established [START_REF] Afraimovich | Synchronization in unidirectionally coupled systems: some rigorous results[END_REF][START_REF] Stark | Regularity of invariant graphs for forced systems Ergod[END_REF]. Here, we complete this result by specifying some properties of this function.

Proposition 4.1 For any f (x, y) = (T (x, y), S(y)) and g(x, z) as above, there exists a function h whose graph z = h(x) attracts all sequences {(x t , z t )} of the skew-product (f, g). The function h has the following properties:

• it is right continuous,

• it is continuous in D c and h(x -0) > h(x) if x ∈ D,

• it is nowhere monotonous,

• it has bounded variation if λ < 1 N and writes h c -h d in this case, where h c is strictly increasing and continuous and h d is an increasing step function.

• it has infinite variation if λ > 1 N where λ := inf • The estimates on the bounded variation follow directly from the fact that if T t 0 (x) ∈ D, then (λ ) t 0 (h(1) -h(0)) ≤ h(x -0) -h(x) ≤ λ t 0 (h(1) -h(0))

Finally the properties of the components h c and h d can be obtained in a similar way as in the proof of Theorem 3.1.

Figure 1 :

 1 Figure 1: Examples of graph of h (of the function h 14 indeed -see equation (4)). In the left pictures (λ < 1/2), h is of bounded variation; in the right ones (λ > 1/2), it has infinite variation. Notice the symmetry h(x) = 1 -h(1 -x -0) for a = 0.5 (bottom pictures).

  and [ θt+1 • • • θn ] are either equal or adjacent. As a result, the intervals h n ([θ 0 • • • θ n ]) and h n ([ θ0 • • • θn ]) must intersect; hence h n ([0, 1]) must be an interval. Computing its extrema gives h n ([0, 1]) = [0, 1 -λ n+1 ].

Figure 2 :

 2 Figure 2: Examples of graphs of h c (blue dots -upper continuous curve) and h d (red dots -lower discontinuous curve) when h is of bounded variation (λ = 0.4); Left a = 0.2; Right a = 0.5.

x

  ,z-z =0 |g(x,z)-g(x,z )| |z-z |.

Total variation and component regularity of the response functionWith such an irregular graph for the synchronization function, the appropriate characteristic to evaluate is the total variation[START_REF] Kolmogorov | Elements of the theory of functions and functional analysis[END_REF]; this quantity

Acknowledgments I am grateful both to G. Mantica for inspiring discussions and especially to A. Quas for pointing out crucial errors in preliminary versions of the manuscript. I thank the Courant Institute (New York University) for hospitality. Support was partly provided by CNRS and by the EU Marie Curie fellowship PIOF-GA-2009-235741.

Proof. Existence proofs in the nonlinear case have already been published in the literature [START_REF] Afraimovich | Synchronization in unidirectionally coupled systems: some rigorous results[END_REF][START_REF] Stark | Regularity of invariant graphs for forced systems Ergod[END_REF]. We provide another proof here that is more suitable to the present framework.

Let g x (y) := g(x, y) and M = sup

where M λ = M 1-λ . Given an arbitrary pair (x, z) and t ≥ 0, define the tth iterate h t (x, z) as follows

By refining the arguments above, one shows that the interval [-M λ , M λ ] is not only invariant for h t but it is also absorbing. In particular every sequence {h t (x, z)} must be bounded.

Next choose two integers t > s. We have

which implies that {h t (x, z)} is a Cauchy sequence; hence the following limit exists for all

and is independent of z. Moreover, the continuous dependence on z implies the following conjugacy equation, i.e.

g(x, h(x)) = g x ( lim

for any y ∈ [0, 1]. Global attraction easily follows by using contraction once again.

The arguments for properties of h are very similar to those in the linear case:

• Right continuity follows from both the fact that all h t (x, y) are right continuous and uniform convergence in the definition of h.

• The existence of left limit h(x-0) and the continuity in D c are obtained similarly by using also the monotonicity of the g(•, z). Now, if T t 0 (x) ∈ D, then by continuity outside D, we have

and thus h(x -0) -h(x) > 0 by strict monotonicity.


Aba Diop
email: aba.diop@univ-lr.fr

Aliou Diop
email: alioudiop52@yahoo.fr

Jean-Franc ¸ois Dupuy
email: jean-francois.dupuy@univ-lr.fr

Maximum likelihood estimation in the logistic regression model with a cure fraction

Keywords: Zero-inflation, Maximum likelihood estimation, Consistency, Asymptotic normality, Simulations

Logistic regression is widely used in medical studies to investigate the relationship between a binary response variable Y and a set of potential predictors X. The binary response may represent, for example, the occurrence of some outcome of interest (Y = 1 if the outcome occurred and Y = 0 otherwise). In this paper, we consider the problem of estimating the logistic regression model with a cure fraction. A sample of observations is said to contain a cure fraction when a proportion of the study subjects (the so-called cured individuals, as opposed to the susceptibles) cannot experience the outcome of interest. One problem arising then is that it is usually unknown who are the cured and the susceptible subjects, unless the outcome of interest has been observed. In this setting, a logistic regression analysis of the relationship between X and Y among the susceptibles is no more straightforward. We develop a maximum likelihood estimation procedure for this problem, based on the joint modeling of the binary response of interest and the cure status. We investigate the identifiability of the resulting model. Then, we establish the consistency and asymptotic normality of the proposed estimator, and we conduct a simulation study to investigate its finite-sample behavior.

Introduction

Logistic regression is widely used to model binary response data in medical studies. An example of a binary response variable is the infection status (infected vs uninfected) with respect to some disease. A logistic regression model can be used to investigate the relationship between the infection status and various potential predictors. If Y i denotes the infection status for the i-th individual in a sample of size n (Y i = 1 if the individual is †Corresponding author infected, and Y i = 0 otherwise), and X i denotes the corresponding (p-dimensional, say) predictor, the logistic regression model expresses the relationship between Y i and X i in term of the conditional probability P(Y i = 1|X i ) of infection, as:

log P(Y i = 1|X i ) 1 -P(Y i = 1|X i ) = β ′ X i ,
where β ∈ R p is an unknown parameter to be estimated. An extensive literature has been devoted so far to statistical inference in logistic regression models. Estimation and testing procedures for this class of models are now well established and are available in standard statistical softwares. In particular, the maximum likelihood estimator of β is obtained by solving the following score equation:

n i=1 X i Y i - e β ′ Xi 1 + e β ′ Xi = 0.
Asymptotic results (consistency and asymptotic normality) for this estimator were given by [START_REF] Gouriéroux | Asymptotic properties of the maximum likelihood estimator in dichotomous logit models[END_REF] and [START_REF] Fahrmeir | Consistency and asymptotic normality of the maximum likelihood estimator in generalized linear models[END_REF], among others. We refer the reader to [START_REF] Hosmer | Applied logistic regression[END_REF] and [START_REF] Hilbe | Logistic regression models[END_REF] for detailed treatments and numerous examples.

In this paper, we consider the problem of estimation in the logistic regression model with a cure fraction. In medical studies, it often arises that a proportion of the study subjects cannot experience the outcome of interest. Such individuals are said to be cured, or immune. The population under study can then be considered as a mixture of cured and susceptible subjects, where a subject is said to be susceptible if he would eventually experience the outcome of interest. One problem arising in this setting is that it is usually unknown who are the susceptible, and the cured subjects (unless the outcome of interest has been observed). Consider, for example, the occurrence of infection from some disease to be the outcome of interest. Then, if a subject is uninfected, the investigator usually does not know whether this subject is immune to the infection, or susceptible albeit still uninfected.

Estimating a regression model with a cure fraction can be viewed as a zero-inflated regression problem. Zero-inflation occurs in the analysis of count data when the observations contain more zeros than expected. Failure to account for these extra zeros is known to result in biased parameter estimates and inferences. The regression analysis of count data with excess zeros has attracted much attention so far. For example, [START_REF] Lambert | Zero-inflated Poisson regression, with an application to defects in manufacturing[END_REF] proposed the zero-inflated Poisson (ZIP) regression model for count data with many zeros. This was further extended to a semiparametric ZIP regression model by [START_REF] Lam | Semiparametric analysis of zero-inflated count data[END_REF]. We refer to [START_REF] Dietz | On estimation of the poisson parameter in zero-modified poisson models[END_REF] and Xiang et al. (2007) for a review of various other extensions of the ZIP model. Other popular models are the zero-inflated binomial (ZIB) regression model (see, for example, [START_REF] Hall | Zero-inflated Poisson and binomial regression with random effects: a case study[END_REF]), and the zero-inflated negative binomial (ZINB) regression model (see, for example, Ridout et al. (2001)). Recently, [START_REF] Kelley | Zero inflation in ordinal data: incorporating susceptibility to response through the use of a mixture model[END_REF] proposed a zero-inflated proportional odds model (ZIPO) for ordinal outcomes, when some individuals are not susceptible to the phenomenon being measured. Various other models and numerous references can be found in [START_REF] Famoye | Zero-inflated generalized Poisson regression model with an application to domestic violence data[END_REF] and [START_REF] Lee | Multi-level zero-inflated Poisson regression modelling of correlated count data with excess zeros[END_REF].

In our paper, we consider the problem of estimating a logistic regression model from binary response data with a cure fraction, when the cure probability is modeled by a logistic regression. This can be viewed as a zero-inflated Bernoulli regression problem, where logistic link functions are used for both the binary response of interest (the probability of infection, say) and the zero-inflation probability (the probability of being cured). The literature on zero-inflated models is extensive but to the best of our knowledge, the theoretical and numerical issues related to the statistical inference in this model have not been yet investigated. In this paper, we intend to fill this gap. We first investigate the identifiability question in this model. Then, we turn to the problem of estimation. The estimator we propose is obtained by maximizing the joint likelihood for the binary response of interest and the cure indicator. We prove the almost sure asymptotic existence, the consistency, and the asymptotic normality of this estimator. Then, we investigate its finite-sample properties via simulations.

The rest of this paper is organized as follows. In Section 2, we describe the problem of logistic regression with a cure fraction, and we propose an estimation method adapted to this setting. The proposed procedure is based on a joint regression model for the binary response of interest and the cure indicator. In Section 3, we investigate the identifiability of this model, and we state some regularity conditions. In Section 4, we derive the asymptotic properties of the resulting estimator. Section 5 describes a simulation study, where we numerically investigate the small to large sample properties of this estimator. A real data example illustrates the methodology. A discussion and some perspectives are given in Section 6.

Logistic regression with a cure fraction

Notations and the model set-up

Let (Y 1 , S 1 , X 1 , Z 1 ), . . . , (Y n , S n , X n , Z n ) be independent and identically distributed copies of the random vector (Y, S, X, Z) defined on the probability space (Ω, A, P). For every individual i = 1, . . . , n, Y i is a binary response variable indicating say, the infection status with respect to some disease (that is, Y i = 1 if the i-th individual is infected, and Y i = 0 otherwise), and S i is a binary variable indicating whether individual i is susceptible to the infection (S i = 1) or immune (S i = 0). If Y i = 0, then the value of S i is unknown. Let X i = (1, X i2 , . . . , X ip ) ′ and Z i = (1, Z i2 , . . . , Z iq ) ′ be random vectors of predictors or covariates (both categorical and continuous predictors are allowed). We shall assume in the following that the X i 's are related to the infection status, while the Z i 's are related to immunity. X i and Z i are allowed to share some components.

The logistic regression model for the infection status assumes that the conditional probability P(Y = 1|X i , S i ) of infection is given by log

P(Y = 1|X i , S i ) 1 -P(Y = 1|X i , S i ) = β 1 + β 2 X i2 + . . . + β p X ip := β ′ X i (1)
if {S i = 1}, and by

P(Y = 1|X i , S i ) = 0 (2) if {S i = 0}
, where β = (β 1 , . . . , β p ) ′ ∈ R p is an unknown regression parameter measuring the association between potential predictors and the risk of infection (for a susceptible individual).

The statistical analysis of infection data with model (1) includes estimation and testing for β. Without immunity (that is, if S i = 1 for every i = 1, . . . , n), inference on β from the sample (Y 1 , X 1 , Z 1 ), . . . , (Y n , X n , Z n ) can be based on the maximum likelihood principle. When immunity is present, deriving the maximum likelihood estimator of β is no longer straightforward: if Y i = 0, we do not know whether {S i = 1}, so that (1) applies, or whether {S i = 0}, so that (2) applies.

One solution is to consider every individual i such that {Y i = 0} as being susceptible that is, to ignore a possible immunity of this individual. We may however expect this method to produce biased estimates of the association of interest (such a method will be evaluated in the simulation study described in section 5). Therefore in this paper, we aim at providing an alternative estimation procedure for β. This can be achieved if a model for immunity is available, as is explained in the next section.

The proposed estimation procedure

A model for the immunity status is defined through the conditional probability P(S = 1|Z i ) of being susceptible to the infection. A common choice for this is the logistic model (see, for example, [START_REF] Fang | Maximum likelihood estimation in a semiparametric logistic/proportional-hazards mixture model[END_REF] and [START_REF] Lu | Maximum likelihood estimation in the proportional hazards cure model[END_REF][START_REF] Lu | Efficient estimation for an accelerated failure time model with a cure fraction[END_REF] who considered estimation in various survival regression models with a cure fraction):

log P(S = 1|Z i ) 1 -P(S = 1|Z i ) = θ 1 + θ 2 Z i2 + . . . + θ q Z iq := θ ′ Z i (3)
where θ = (θ 1 , . . . , θ q ) ′ ∈ R q is an unknown regression parameter.

Remark 1. We note that the model defined by ( 1)-( 2)-( 3) can be viewed as a zero-inflated Bernoulli regression model, with logit links for both the binary response of interest and the zero-inflation component. As far as we know, no theoretical investigation of this model has been undertaken yet. Such a work is carried out in the following.

From (1), (2), and (3), a straightforward calculation yields that

P(Y = 1|X i , Z i ) = e β ′ Xi+θ ′ Zi (1 + e β ′ Xi )(1 + e θ ′ Zi )
.

Let ψ := (β ′ , θ ′ ) ′ denote the unknown k-dimensional (k = p+q) parameter in the conditional distribution of Y given X i and Z i . ψ includes both β (considered as the parameter of interest) and θ (considered as a nuisance parameter). Now, the likelihood for ψ from the independent sample (Y i , S i , X i , Z i ) (i = 1, . . . , n) (where S i is unknown when Y i = 0) is as follows:

L n (ψ) = n i=1    e β ′ Xi+θ ′ Zi (1 + e β ′ Xi )(1 + e θ ′ Zi ) Yi 1 - e β ′ Xi+θ ′ Zi (1 + e β ′ Xi )(1 + e θ ′ Zi ) 1-Yi    .
We define the maximum likelihood estimator

ψ n := ( β ′ n , θ ′ n ) ′ of ψ as the solution (if it exists) of the k-dimensional score equation ln (ψ) = ∂l n (ψ) ∂ψ = 0, ( 4
)
where l n (ψ) := log L n (ψ) is the log-likelihood function. In the following, we shall be interested in the asymptotic properties of the maximum likelihood estimator β n of β, considered as a sub-component of ψ n . We will however obtain consistency and asymptotic normality results for the whole ψ n . Before proceeding, we need to set some further notations.

Some further notations

Define first the (p × n) and (q × n) matrices

X =      1 1 • • • 1 X 12 X 22 • • • X n2 . . . . . . . . . . . . X 1p X 2p • • • X np      and Z =      1 1 • • • 1 Z 12 Z 22 • • • Z n2 . . . . . . . . . . . . Z 1q Z 2q • • • Z nq      ,
and let W be the (k × 2n) block-matrix defined as

W = X 0 pn 0 qn Z ,
where 0 ab denotes the (a×b) matrix whose components are all equal to zero (for any positive integer values a, b). Let also C(ψ) be the 2n-dimensional column vector defined as

C(ψ) = (A β (ψ) -B β (ψ)) ′ , (A θ (ψ) -B θ (ψ)) ′ ′ , where A β (ψ) = (A β i (ψ)) 1≤i≤n , B β (ψ) = (B β i (ψ)) 1≤i≤n , A θ (ψ) = (A θ i (ψ)) 1≤i≤n
, and B θ (ψ) = (B θ i (ψ)) 1≤i≤n are n-dimensional column vectors with respective elements

A β i (ψ) = 1 + e θ ′ Zi 1 + e β ′ Xi + e θ ′ Zi Y i , B β i (ψ) = e β ′ Xi+θ ′ Zi (1 + e β ′ Xi )(1 + e β ′ Xi + e θ ′ Zi ) , A θ i (ψ) = 1 + e β ′ Xi 1 + e β ′ Xi + e θ ′ Zi Y i , B θ i (ψ) = e β ′ Xi+θ ′ Zi (1 + e θ ′ Zi )(1 + e β ′ Xi + e θ ′ Zi )
.

Then, simple algebra shows that the score equation can be rewritten as ln (ψ) = WC(ψ) = 0.

If M = (M ij ) 1≤i≤a,1≤j≤b denotes some (a×b) matrix, we will denote by M •j its j-th column (j = 1, . . . , b) that is, M •j = (M 1j , . . . , M aj ) ′ . Then, it will be useful to rewrite the score vector as

ln (ψ) = 2n j=1 W •j C j (ψ).
We shall further note ln (ψ) the (k × k) matrix of second derivatives of l n (ψ) that is, ln (ψ) = ∂ 2 l n (ψ)/∂ψ∂ψ ′ . Let D(ψ) = (D ij (ψ)) 1≤i,j≤2n be the (2n × 2n) block matrix defined as

D(ψ) = D 1 (ψ) D 3 (ψ) D 3 (ψ) D 2 (ψ) ,
where D 1 (ψ), D 2 (ψ), and D 3 (ψ) are (n × n) diagonal matrices, with i-th diagonal elements (i = 1, . . . , n) respectively given by

D 1,ii (ψ) = e β ′ Xi+θ ′ Zi (1 + e β ′ Xi ) 2 (1 + e β ′ Xi + e θ ′ Zi ) , D 2,ii (ψ) = e β ′ Xi+θ ′ Zi (1 + e θ ′ Zi ) 2 (1 + e β ′ Xi + e θ ′ Zi ) , D 3,ii (ψ) = e β ′ Xi+θ ′ Zi (1 + e β ′ Xi )(1 + e θ ′ Zi )(1 + e β ′ Xi + e θ ′ Zi )
.

Then, some algebra shows that ln (ψ) can be expressed as

ln (ψ) = -WD(ψ)W ′ .
Note that the size of C(ψ), W, and D(ψ) depends on n. However, in order to simplify notations, n will not be used as a lower index for these vector and matrices. In the next section, we investigate the question of parameter identifiability in model ( 1)-( 2)-(3).

Identifiability and regularity conditions

We first state some regularity conditions that will be needed to ensure identifiability and the asymptotic results in Section 4:

C1 The covariates are bounded that is, there exist compact sets F ⊂ R p and G ⊂ R q such that X i ∈ F and Z i ∈ G for every i = 1, 2, . . . For every i = 1, 2, . . ., j = 2, . . . , p, k = 2, . . . , q, var[X ij ] > 0 and var[Z ik ] > 0. For every i = 1, 2, . . ., the X ij (j = 1, . . . , p) are linearly independent, and the Z ik (k = 1, . . . , q) are linearly independent.

C2 Let ψ 0 = (β ′ 0 , θ ′ 0 ) ′ denote the true parameter value. β 0 and θ 0 lie in the interior of known compact sets B ⊂ R p and G ⊂ R q respectively. C3 The Hessian matrix ln (ψ) is negative definite and of full rank, for every n = 1, 2, . . . Let λ n and Λ n be respectively the smallest and largest eigenvalues of WD(ψ 0 )W ′ . There exists a finite positive constant c 2 such that Λ n /λ n < c 2 for every n = 1, 2, . . .

C4

There exists a continuous covariate V which is in X but not in Z that is, if β V and θ V denote the coefficients of V in the linear predictors (1) and (3) respectively, then β V = 0 and θ V = 0. At a model-building stage, it is known that V is in X.

The conditions C1, C2, C3 are classical conditions for identifiability and asymptotic results in standard logistic regression (see, for example, [START_REF] Gouriéroux | Asymptotic properties of the maximum likelihood estimator in dichotomous logit models[END_REF] and [START_REF] Guyon | Statistique et économétrie -Du modèle linéaire aux modèles non-linéaires[END_REF]). The condition C4, which imposes some restrictions on the covariates, is required for identifiability of ψ in the joint model ( 1)-( 2)-(3) (we may alternatively assume that the continuous covariate V is in Z but not in X). In the following, we will assume that V is in X but not in Z, with β V := β l for some l ∈ {2, . . . , p}, and for the i-th individual, we will denote V i by X il . The condition C4 is discussed in greater details in the following two remarks.

Remark 2. We may relate the identifiability issue in model ( 1)-( 2)-( 3) to the problem of identifiability of mixtures of logistic regression models, which was investigated by [START_REF] Follmann | Identifiability of finite mixtures of logistic regression models[END_REF]. [START_REF] Follmann | Identifiability of finite mixtures of logistic regression models[END_REF] considered the case where there is a finite number c of components in the mixture (we consider here the case where c = 2, with one degenerate component) and the mixing probabilities are constant (here, the mixing probabilities given by ( 3) are allowed to depend on covariates). The authors have shown that finite mixtures of logistic regressions are identifiable provided that the number of unique covariate combinations values is sufficiently large. C4 can be viewed as a sufficient condition for achieving the same kind of requirement. A similar condition appears in [START_REF] Kelley | Zero inflation in ordinal data: incorporating susceptibility to response through the use of a mixture model[END_REF].

To understand C4, note that if X i = Z i , then exchanging the parameters β and θ in (1) and

(3) yields the same likelihood value L n (ψ), which is a cause of model non-identifiability. A similar remark holds if we invert the linear predictors β ′ X i and θ ′ Z i . The condition C4 evacuates these problems.

First, by asking one of the covariates to be significant in one and only one linear predictor, C4 prevents β ′ X and θ ′ Z from being of the same form, and the parameters are thus not exchangeable. Secondly, by assuming that we know, prior to model fitting, that there exists a covariate V which is in X but not in Z, C4 will force each linear predictor to be attached to the correct corresponding model ( 1) or (3). These facts are illustrated in a web-based supplementary document available at the following address: http://perso.univ-lr.fr/jfdupuy/supplementary.pdf. There, we provide the results of a simulation study which investigates numerically the identifiability of model ( 1)-( 2)-( 3). For each of the models considered in this study, we assume that C4 is satisfied: the linear predictors β ′ X i and θ ′ Z i share three covariates (one is continuous, two are discrete), and an additional continuous covariate is included in X i . Using the procedure described in Section 2, maximum likelihood estimates are obtained for β and θ, and are averaged over N = 1000 samples (we considered several combinations of sample size, proportion of immunes, proportion of infected among the susceptibles). Both parameters β and θ appear to be identifiable (the averaged estimates appear to be close to the true parameters, including those corresponding to the three shared covariates).

Remark 3. The condition C4 does not appear to be too restrictive in practice. Consider the example of the transmission of some disease by breastfeeding. If every child in the sample is breastfeeded, it can be expected that the length (in days, say) of the breastfeeding period (a continuous covariate) will influence the probability of infection, while the susceptibility probability will rather depend on risk factors such as say, the mother's infection status. It is also worth noting that the consequences of C4, in terms of model-building, are rather mild. At a model-building stage, we may be tempted to incorporate all available covariates in both linear predictors (1) and (3), and to remove irrelevant factors by using backward elimination. The condition C4 slightly restricts this fitting strategy, by imposing that one relevant continuous covariate is incorporated in one (and only one) linear predictor. This should often be doable in practice, since the statistician often gets some prior knowledge (from the clinicians, epidemiologists,. . . ) about the dataset to be analyzed.

We are now in position to prove the following result:

Theorem 1 (Identifiability). Under the conditions C1-C4, the model ( 1)-( 2

)-(3) is iden- tifiable; that is, L 1 (ψ) = L 1 (ψ * ) almost surely implies ψ = ψ * .
Proof of Theorem 1. Suppose that L 1 (ψ) = L 1 (ψ * ) almost surely. Under C1 and C2, there exists a positive constant c 1 such that for every x ∈ F , z ∈ G, and

ψ ∈ B × G, c 1 < P(Y = 1|x, z) < 1 -c 1 . Thus we can find a ω ∈ Ω, outside the negligible set where L 1 (ψ) = L 1 (ψ * ), and such that Y (ω) = 1 when X = x and Z = z. For this ω, L 1 (ψ) = L 1 (ψ * ) becomes e β ′ x+θ ′ z (1 + e β ′ x )(1 + e θ ′ z ) = e β * ′ x+θ * ′ z (1 + e β * ′ x )(1 + e θ * ′ z ) .
This can be rewritten as

1 + e -β ′ x 1 + e -β * ′ x = 1 + e -θ * ′ z 1 + e -θ ′ z . (5)
Now, under the condition C4, taking the partial derivative of both sides of (5) with respect to the l-th component of x (X il is a continuous covariate) yields

-β l e -β ′ x (1 + e -β * ′ x ) + β * l e -β * ′ x (1 + e -β ′ x ) (1 + e -β * ′ x ) 2 = 0
since the right-hand-side of (5) does not depend on x. Thus, it follows that

β l β * l = 1 + e β ′ x 1 + e β * ′ x .
Differentiating both sides of this equality with respect to the l-th component of x further yields (ββ * ) ′ x = 0, which implies that β = β * under C1. It remains to show that θ = θ * , which reduces to the identifiability problem in the standard logistic regression model. We have that θ = θ * under C1 (see [START_REF] Guyon | Statistique et économétrie -Du modèle linéaire aux modèles non-linéaires[END_REF] for example), which concludes the proof.

2

We now turn to the asymptotic theory for the proposed estimator.

Asymptotic theory

In this section, we establish rigorously the existence, consistency and asymptotic normality of the maximum likelihood estimator β n of β in model ( 1), obtained from a sample of binary response data with a cure fraction. In the sequel, the space R k of k-dimensional (column) vectors will be provided with the Euclidean norm, and the space R k×k of (k × k) real matrices will be provided with the spectral norm (we will use the same notation • for both). We first prove the following result:

Theorem 2 (Existence and consistency). Under the conditions C1-C3, the maximum likelihood estimator ψ n exists almost surely as n → ∞, and converges almost surely to ψ 0 , if and only if λ n tends to infinity as n → ∞.

Proof of Theorem 2. The principle of the proof is similar to [START_REF] Gouriéroux | Asymptotic properties of the maximum likelihood estimator in dichotomous logit models[END_REF] but the technical details are different. Three lemmas are needed. The first lemma essentially provides an intermediate technical result. Its proof is postponed to the appendix.

Lemma 1. Let φ n : R k -→ R k be defined as

φ n (ψ) = ψ + (WD(ψ 0 )W ′ ) -1 ln (ψ).
Then there exists an open ball B(ψ 0 , r) (with r > 0) such that φ n satisfies the Lipschitz condition on B(ψ 0 , r) that is,

φ n (ψ 1 ) -φ n (ψ 2 ) ≤ c ψ 1 -ψ 2 for all ψ 1 , ψ 2 ∈ B(ψ 0 , r), (6)
and 0 < c < 1.

Lemma 2. The maximum likelihood estimator ψ n exists almost surely as n → ∞, and converges almost surely to ψ 0 , if and only if (WD(ψ 0 )W ′ ) -1 ln (ψ 0 ) converges almost surely to 0.

Proof of Lemma 2. We first prove that the condition is sufficient. Thus, we assume that (WD(ψ 0 )W ′ ) -1 ln (ψ 0 ) converges almost surely to 0. Define

η n (ψ) = ψ -φ n (ψ) = -(WD(ψ 0 )W ′ ) -1
ln (ψ) and let ǫ be an arbitrary positive value. Then for almost every ω ∈ Ω, there exists an integer value n(ǫ, ω) such that for any n ≥ n(ǫ, ω), η n (ψ 0 ) ≤ ǫ or equivalently, 0 ∈ B(η n (ψ 0 ), ǫ). In particular, let ǫ = (1c)s with 0 < c < 1 such as in Lemma 1. Since φ n satisfies the Lipschitz condition (6) (by Lemma 1), the lemma 2 of [START_REF] Gouriéroux | Asymptotic properties of the maximum likelihood estimator in dichotomous logit models[END_REF] ensures that there exists an element of B(ψ 0 , s) (let denote this element by

ψ n ) such that η n ( ψ n ) = 0 that is, (WD(ψ 0 )W ′ ) -1 ln ( ψ n ) = 0.
The condition C3 implies that ln ( ψ n ) = 0 and that ψ n is the unique maximizer of l n . To summarize, we have shown that for almost every ω ∈ Ω and for every s > 0, there exists an integer value n(s, ω) such that if n ≥ n(s, ω), then the maximum likelihood estimator ψ n exists, and ψ nψ 0 ≤ s (that is, ψ n converges almost surely to ψ 0 ). We now prove that the condition that η n (ψ 0 ) converges almost surely to 0 is necessary. We use a proof by contradiction.

Assume that as n → ∞, ψ n exists and converges almost surely to ψ 0 , but η n (ψ 0 ) does not converge almost surely to 0. Then there exists a set Ω ⊂ Ω with P( Ω) > 0, such that if ω ∈ Ω, there exists ǫ > 0 such that for every m ∈ N, there exists n ≥ m with η n (ψ 0 ) > ǫ. Now, let t = ǫ d(1+c) , with d > 1 sufficiently large so that t ≤ r, where r is such as in Lemma 1. Then for every ψ ∈ B(ψ 0 , t), the following holds:

η n (ψ 0 ) -η n (ψ) = ψ 0 -φ n (ψ 0 ) -ψ + φ n (ψ) ≤ ψ 0 -ψ + φ n (ψ) -φ n (ψ 0 ) ≤ t(1 + c) = ǫ d ,
where the second to third line follows by Lemma 1. Therefore, for every ψ ∈ B(ψ 0 , t),

ǫ < η n (ψ 0 ) ≤ η n (ψ 0 ) -η n (ψ) + η n (ψ) ≤ η n (ψ) + ǫ d
and we conclude that for every ψ ∈ B(ψ 0 , t), η n (ψ) > ǫ(1 -1 d ) > 0. Since η n ( ψ n ) = 0, ψ n cannot belong to B(ψ 0 , t) for large n, which implies that ψ n does not converge almost surely to ψ 0 . This is the desired contradiction.

2 Lemma 3. (WD(ψ 0 )W ′ ) -1
ln (ψ 0 ) converges almost surely to 0 if and only if λ n tends to infinity as n → ∞.

Proof of Lemma 3. We first prove that the condition is sufficient that is, we assume that λ n tends to infinity as n → ∞. Define the (2n × k) matrix V = (D(ψ 0 ))

1 2 W ′ and the 2n-dimensional vector U = (D(ψ 0 )) -1 2 C(ψ 0 ). Then E[U ] = 0 and var[U ] = I 2n , (7)
where I 2n denotes the identity matrix of order 2n. To see this, note that

E[U ] = E[E[(D(ψ 0 )) -1 2 C(ψ 0 )|X, Z]] = E[(D(ψ 0 )) -1 2 E[C(ψ 0 )|X, Z]] = E[(D(ψ 0 )) -1 2 E[ (A β (ψ 0 ) -B β (ψ 0 )) ′ , (A θ (ψ 0 ) -B θ (ψ 0 )) ′ ′ |X, Z]].
For every i = 1, . . . , n,

E[A β i (ψ 0 ) -B β i (ψ 0 )|X, Z]] = E[A β i (ψ 0 ) -B β i (ψ 0 )|X i , Z i ]
] by independence between the individuals, and By [START_REF] Gouriéroux | Asymptotic properties of the maximum likelihood estimator in dichotomous logit models[END_REF] (proof of Lemma 4), if (7) holds, Λ n /λ n < c 2 for every n = 1, 2, . . ., and λ n tends to infinity as n → ∞, then

E[A β i (ψ 0 ) -B β i (ψ 0 )|X i , Z i ]] = 1 + e θ ′ 0 Zi 1 + e β ′ 0 Xi + e θ ′ 0 Zi P(Y i = 1|X i , Z i ) -B β i (ψ 0 ) = B β i (ψ 0 ) -B β i (ψ 0 ) = 0. Similarly, E[A θ i (ψ 0 ) -B θ i (ψ 0 )|X, Z]] = 0 for every i = 1, . . . , n and thus, E[C(ψ 0 )|X, Z]] = 0 and E[U ] = 0. Next, var[U ] = E[var[U |X, Z]] since E[U |X, Z] = 0. Moreover, var[U |X, Z] = (D(ψ 0 )) -1 2 var[C(ψ 0 )|X, Z](D(ψ 0 )) -1 2 , with var[C(ψ 0 )|X, Z] = var[ A β (ψ 0 ) ′ , A θ (ψ 0 ) ′ ′ |X, Z] a (2n × 2n) block-matrix of the form V 1 V 3 V 3 V 2 where V 1 , V 2 , and V 3 are (n × n) matrices. The i-th diagonal elements (i = 1, . . . , n) of V 1 , V 2 ,
(V ′ V) -1 V ′ U a.s. -→ 0 as n → ∞ that is, (WD(ψ 0 )W ′ ) -1
ln (ψ 0 ) converges almost surely to 0. We now prove that the condition is necessary. Assume that λ n does not tend to infinity as n → ∞. By [START_REF] Gouriéroux | Asymptotic properties of the maximum likelihood estimator in dichotomous logit models[END_REF]

(proof of Lemma 4), (V ′ V) -1 V ′ U (and therefore (WD(ψ 0 )W ′ ) -1
ln (ψ 0 )) cannot converge to 0, which concludes the proof.

2

Finally, Theorem 2 follows by Lemma 2 and Lemma 3.

2

We now turn to the convergence in distribution of the proposed estimator, which is stated by the following theorem:

Theorem 3 (Asymptotic normality). Assume that the conditions C1-C3 hold and that ψ n converges almost surely to ψ 0 . Let Σ n = WD( ψ n )W ′ and I k denote the identity matrix of order k. Then Σ 1 2

n ( ψ nψ 0 ) converges in distribution to the Gaussian vector N (0, I k ).

Proof of Theorem 3. A Taylor expansion of the score function is as

0 = ln ( ψ n ) = ln (ψ 0 ) + ln ( ψ n )( ψ n -ψ 0 )
where ψ n lies between ψ n and ψ 0 , and thus ln (

ψ 0 ) = -ln ( ψ n )( ψ n -ψ 0 ). Let Σ n := -ln ( ψ n ) = WD( ψ n )W ′ and Σ n,0 := WD(ψ 0 )W ′ . Now, Σ 1 2 n ( ψ n -ψ 0 ) = Σ 1 2 n Σ -1 2 n Σ -1 2 n Σ 1 2 n,0 Σ -1 2 n,0 Σ n ( ψ n -ψ 0 ) . (8)
The two terms in brackets in (8) converge almost surely to I k . To see this, we show for example that Σ

-1 2 n Σ 1 2 n,0 -I k a.s. -→ 0 as n → ∞. First, note that Σ -1 2 n Σ 1 2 n,0 -I k ≤ Λ 1 2 n Σ -1 2 n Λ -1 2 n Σ 1 2 n,0 -Σ 1 2 n , (9)
and

Λ -1 n Σ n,0 -Σ n = Λ -1 n W(D(ψ 0 ) -D( ψ n ))W ′ .
Note also that ψ n converges almost surely to ψ 0 (that is, for every ω ∈ Ω, where Ω ⊂ Ω and P( Ω) = 1). Let ω ∈ Ω. By the same arguments as in the proof of Lemma 1, for every ǫ > 0, there exists a positive n(ǫ, ω) ∈ N such that if n ≥ n(ǫ, ω), then [START_REF] Eicker | A multivariate central limit theorem for random linear vector forms[END_REF], this convergence holds if we can check the following conditions:

Λ -1 n W(D(ψ 0 ) -D( ψ n ))W ′ ≤ ǫ. Hence Λ -1 n W(D(ψ 0 ) -D( ψ n ))W ′ converges

It remains for us to show that Σ

-1 2 n,0 ( Σ n ( ψ n -ψ 0 )) converges in distribution to N (0, I k ), or equivalently, that (V ′ V) -1 2 V ′ U converges in distribution to N (0, I k ). Following
i) max 1≤i≤2n V i• Σ -1 n,0 V ′ i• -→ 0 as n → ∞, ii) sup 1≤i≤2n E[U 2 i 1 {|Ui|>α} ] -→ 0 as α → ∞, iii) inf 1≤i≤2n E[U 2 i ] > 0, where V i• and U i respectively denote the i-th raw of V and the i-th component of U , i = 1, . . . , 2n. Condition i) follows by noting that 0 ≤ max 1≤i≤2n V i• Σ -1 n,0 V ′ i• ≤ max 1≤i≤2n V i• 2 Σ -1 n,0 = max 1≤i≤2n 1 λ n V i• 2 ,
and that V i• is bounded above, by C1 and C2. Moreover, 1 λn tends to 0 as n → ∞, since ψ n converges almost surely to ψ 0 . Condition ii) follows by noting that the components U i of U are bounded under C1 and C2. Finally, for every i = 1, . . . , 2n, E[U 2 i ] = var[U i ] since U is centered. We have proved (see Lemma 3) that var[U ] = I 2n , thus for every i = 1, . . . , 2n, var[U i ] = 1, and finally, inf 1≤i≤2n E[U 2 i ] = 1 > 0. To summarize, we have proved that Σ -1 2 n,0 ( Σ n ( ψ n -ψ 0 )) converges in distribution to N (0, I k ). This result, combined with Slutsky's theorem and equation ( 8), implies that Σ 1 2 n ( ψ nψ 0 ) converges in distribution to N (0, I k ).

2

A simulation study and real data example

Study design

In this section, we investigate the numerical properties of the maximum likelihood estimator β n , under various conditions. The simulation setting is as follows. We consider the following models for the infection status:

log P(Y =1|Xi,Si) 1-P(Y =1|Xi,Si) = β 1 + β 2 X i2 if S i = 1 P(Y = 1|X i , S i ) = 0 if S i = 0
and the immunity status:

log P(S = 1|Z i ) 1 -P(S = 1|Z i ) = θ 1 + θ 2 Z i2 ,
where X i2 is normally distributed with mean 0 and variance 1, and Z i2 is normally distributed with mean 1 and variance 1. An i.i.d. sample of size n of the vector (Y, S, X, Z) is generated from this model, and for each individual i, we get a realization (y i , s i , x i , z i ), where s i is considered as unknown if y i = 0. A maximum likelihood estimator β n of β = (β 1 , β 2 ) is obtained from this incomplete dataset by solving the score equation (4), using the optim function of the software R. An estimate is also obtained for θ = (θ 1 , θ 2 ), but θ is not the primary parameter of interest hence we only focus on the simulation results for β n . The finite-sample behavior of the maximum likelihood estimator β n was assessed for several sample sizes (n = 100, 500, 1000, 1500) and various values for the percentage of immunes in the sample, namely 25%, 50%, and 75%. The case where it is known that there are no immunes in the sample was also considered. In this case, there is no missing information about the infection status and therefore, this case provides a benchmark for evaluating the performance of the proposed estimation method. We also considered different values for the proportion of infected individuals among the susceptibles. The desired proportions of immunes and infected were obtained by choosing appropriate values for the parameters β (the parameter of interest) and θ (the nuisance parameter). The following values were considered for β: i) model M 1 : β = (-.8, 1) (using these values, approximately 30% of the susceptibles are infected), ii) model M 2 : β = (1, .7) (approximately 70% of the susceptibles are infected), iii) model M 3 : β = (-.8, 0) (approximately 30% of the susceptibles are infected), iv) model M 4 : β = (1, 0) (approximately 70% of the susceptibles are infected).

Results

For each configuration (sample size, percentage of immunes, percentage of infected among susceptibles) of the design parameters, N = 1500 samples were obtained. Based on these 1500 repetitions, we obtain averaged values for the estimates of β 1 and β 2 , which are calculated as

N -1 N j=1 β (j) 1,n and N -1 N j=1 β (j)
2,n , where β

(j) n = ( β (j) 1,n , β (j)
2,n ) is the estimate obtained from the j-th simulated sample. For each of the parameters β 1 and β 2 , we also obtain the empirical root mean square and mean absolute errors, based on the N samples. When β 2 = 0 (respectively β 2 = 0), we obtain the empirical power (respectively the empirical size) of the Wald test at the 5% level for testing H 0 : β 2 = 0 (models M 1 and M 2 , see Tables 1 and2) (respectively models M 3 and M 4 , see Tables 1 and2). The null hypothesis H 0 : β 2 = 0 is the hypothesis that the predictor X 2 does not influence the risk of infection of susceptible individuals. The results are summarized in Tables 1 and2.

Tables 1 to 2 about here From these tables, it appears that the proposed maximum likelihood estimator β n provides a reasonable approximation of the true parameter value, even when the percentage of immunes is high. While the bias of β n stays limited, its variability increases with the immune fraction, sometimes drastically when the sample size is small. Consequently, when the sample size is small (n = 100) and/or the immune proportion is very high (75%), the power of the Wald test for nullity of the regression coefficient β 2 can be low, compared to the case where there are no immunes. But we note that for moderately large to large sample sizes (n ≥ 500), the dispersion indicators and the power of the Wald test indicate good performance of the maximum likelihood estimate, even when the immune proportion is up to 50%. The level of the Wald test for nullity of β 2 is globally respected except, for every immune proportion, when the sample size is small (n = 100).

We compare these results to the ones obtained from a "naive" method where: i) we consider every individual i such that {Y i = 0} as being susceptible but uninfected, that is we ignore the eventual immunity of this individual, ii) we apply a usual logistic regression analysis to the resulting dataset. The results of such "naive" analysis for model M 1 are given in Table 3 (the results for models M 2 , M 3 , M 4 yield similar observations and thus, they are not given here. However, the complete simulation study is available from the web-based supplementary document mentioned above).

Table 3 about here

From this table, it appears that ignoring the immunity present in the sample results in strongly biased estimates of β. The bias of the intercept estimate increases with the immune proportion. At the same time, the estimate of the regression coefficient β 2 is biased towards 0 for all values of the immune percentage and sample size. This results in a very low power for the Wald test of nullity of β 2 , and in a wrong interpretation of the relationship between the covariate X 2 and the binary response Y .

The quality of the Gaussian approximation to the large-sample distribution of β 2,n was also investigated. For each configuration of the design parameters, histograms of the β (j) 2,n (j = 1, . . . , N ) are obtained, along with the corresponding QQ-plots. The plots for the model M 1 are pictured on Figures 1 to 4 (the plots for the models M 2 , M 3 , M 4 are given in the web-based file).

Figures 1 to 4 about here

From these figures, it appears that the normal approximation stated in Theorem 3 is reasonably satisfied when the proportion of immunes is moderate (25%), provided that the sample size is sufficiently large (n ≥ 500, say). Consider the case when β 2 = 0. When the immune fraction is large (50%), the normal approximation still appears reasonable, provided that the sample size is at least 1000, or eventually 1500. When the immune proportion is very large (75%), the distribution of β 2,n can be highly skewed, in particular when the sample size is small. Consider the case when β 2 = 0. Then the finite-sample distribution of β 2,n appears to be symmetric, with heavy tails however, especially when the sample size is small. When the immune fraction is about 50% and the sample size is greater than or equal to 500, the normal distribution appears to fit reasonably well the distribution of β 2,n . Overall, these results indicate that a reliable statistical inference on the regression effect in the model ( 1) with a cure fraction should be based on a sample having, at least, a moderately large size (n ≥ 500, say) when the immune fraction is moderate (25%), or a large size (n ≥ 1000, say) when the immune proportion is large (50%).

A real data example

In this application, we consider a study of dengue fever, which is a mosquito-borne viral human disease. A dengue infection confers a partial and transient immunity against a subsequent infection (see [START_REF] Dussart | Study of dengue cases and the members of their households: a familial cluster analysis in the multinational denframe project[END_REF]). We consider here a database of size n = 528, which was constituted with individuals recruited in Cambodia, Vietnam, French Guiana, and Brazil [START_REF] Dussart | Study of dengue cases and the members of their households: a familial cluster analysis in the multinational denframe project[END_REF]). Each individual i was diagnosed for dengue infection and coded as Y i = 1 if infection was present and 0 otherwise. Note that if Y i = 0, then the i-th individual may either be immune at the time of analysis (due to a temporary immunity acquired following a previous infection) or susceptible to dengue infection, albeit not infected. We aim at estimating the risk of infection for those individuals, based on this data set which also includes the following covariates: age (a continuous bounded covariate) and weight (coded 0 in case of underweight and 1 otherwise). We first ran a standard logistic regression analysis of the model logit P(Y = 1|age, weight) = β 1 + β 2 age + β 3 weight. The results are displayed in Table 4. Then, we estimated the parameters β 1 , β 2 and β 3 using the methodology described in Section 2. Note first that the eventual immunity imparted by a past infection is only transient, thus there is no reason why an older individual (who has therefore been exposed longer to the risk of dengue fever) would have a greater probability of being immune than a younger one. In fact, individual susceptibility to the dengue infection may rather depend on whether the individual benefits or not from some preventive and control measures (such as the application of insecticides to larval habitats in his area, or appropriate water storage and waste disposal practices). Such informations are not available in our dataset. Age was therefore taken as the variable V in condition C4, and we fitted to the data the model ( 1)-( 2)-(3) with logit P(Y = 1|age, weight, S = 1) = β 1 + β 2 age + β 3 weight and logit P(S = 1| weight) = θ 1 + θ 2 weight. Since the Wald-type test of "θ 2 = 0" was not significant, we removed the weight from the model for susceptibility, resulting in a constant proportion of immunes. The final results of this fitting procedure are given in Table 4.

Table 4 about here

The fitted model produced the following estimate: 1-exp(0.497)/(1+exp(0.497)) ≈ 0.38 for the probability of being immune. Then, as expected, the estimated probabilities of infection obtained from our approach are larger than the ones derived from a standard analysis that does not take account of the possible immunity. For example, the probabilities of infection for individuals respectively aged 30 and 10 years old, with "normal" weight, are estimated by 0.29 and 0.55 (standard logistic regression) and 0.31 and 0.86 (from our approach). It is expected that underweighted subjects (those considered to be under a healthy weight) will have higher risks of infection. The probabilities of infection for underweighted individuals respectively aged 30 and 10 years old are estimated by 0.48 and 0.73 (standard logistic regression) and 0.97 and 0.99 (our approach). While both approaches provide the same qualitative conclusions: the probability of dengue infection is higher for younger individuals and in case of underweight (caused by malnutrition for example), they differ on their estimations of the risk of infection. Our approach takes account of the eventual immunity imparted by a past infection and therefore, it is reasonable to think that the resulting estimations of the infection probabilities provide a more realistic picture of the infection risk for this data set. In particular, the estimates provided by our approach suggest that underweight constitutes a major risk factor for dengue infection, irrespectively of age.

Discussion and perspectives

In this paper, we have considered the problem of estimating the logistic regression model from a sample of binary response data with a cure fraction. The estimator we propose is obtained by maximizing a likelihood function, which is derived from a joint regression model for the binary response of interest and the cure indicator, considered as a random variable whose distribution is modeled by a logistic regression (the proposed joint model can thus be viewed as a zero-inflated Bernoulli regression model, with logit links for both the binary response of interest and the zero-inflation component). We have established the existence, consistency, and asymptotic normality of this estimator, and we have investigated its finite-sample properties via simulations. Several open questions now deserve attention. The estimation approach proposed here relies on our ability to correctly specify the model for the binary immunity status. It is therefore of interest to investigate the effect of a misspecification of this model (and in particular, of the link function). The techniques and results by [START_REF] Czado | The effect of link misspecification on binary regression inference[END_REF] may be useful for that purpose. Another issue of interest deals with the inference in the logistic regression model with a cure fraction, in a high-dimensional setting. We have established the theoretical properties of our estimator in a low-dimensional setting that is, when a small number of potential predictors are involved. Several recent contributions (see for example [START_REF] Huang | The iterated lasso for high-dimensional logistic regression[END_REF] and [START_REF] Meier | The group Lasso for logistic regression[END_REF]) have considered the problem of estimation in the logistic model (without cure fraction) when the predictor dimension is much larger than the sample size (this problem arises, for example, in genetic studies where high-dimensional data are generated using microarray technologies). Extending our methodology to this setting constitutes another topic for further research.

Next, define S = {(i, j) ∈ {1, 2, . . . , 2n} 2 |D ij (ψ 0 ) = 0}. Then the following holds:

W(D(ψ 0 ) -D(ψ))W ′ = 2n i=1 2n j=1 W •i W ′ •j (D ij (ψ) -D ij (ψ 0 )) ≤ (i,j)∈S W •i W ′ •j D ij (ψ 0 ) D ij (ψ) -D ij (ψ 0 ) D ij (ψ 0 ) .
From C1 and C2, there exists a real constant c 3 (c 3 > 0) such that D ij (ψ 0 ) > c 3 for every (i, j) ∈ S. Moreover, D ij (•) is uniformly continuous on B × G, thus for every ǫ > 0, there exists a positive r such that for all ψ ∈ B(ψ 0 , r),

|D ij (ψ) -D ij (ψ 0 )| < ǫ. It follows that W(D(ψ 0 ) -D(ψ))W ′ ≤ ǫ c 3 (i,j)∈S W •i W ′ •j D ij (ψ 0 ) ≤ ǫ c 3 tr   (i,j)∈S W •i W ′ •j D ij (ψ 0 )   = ǫ c 3 tr   2n i=1 2n j=1 W •i W ′ •j D ij (ψ 0 )   = ǫ c 3 tr (WD(ψ 0 )W ′ ) ≤ ǫ c 3 Λ n k.
This in turn implies that ∂φn(ψ) ∂ψ ′ ≤ ǫΛnk c3λn < ǫc2k c3 . Now, choosing ǫ = c c3 c2k with 0 < c < 1, we get that ∂φn(ψ) ∂ψ ′ ≤ c for all ψ ∈ B(ψ 0 , r), and the result follows.

2 Ridout, M., J. Hinde, and C. G. B. Demétrio (2001) Note: * : empirical power of the Wald test at the level 5% for testing H 0 : β 2 = 0. For each percentage of immunes, the percentage of infected among the susceptibles is 30%. In the "naive" analysis, every uninfected individual (i.e. Y i = 0) is considered as susceptible.

  and V 3 are var[A β i (ψ 0 )|X, Z], var[A θ i (ψ 0 )|X, Z], and cov[A β i (ψ 0 ), A θ i (ψ 0 )|X, Z] respectively. Similar calculations as above yield: var[A β i (ψ 0 )|X, Z] = D 1,ii (ψ 0 ), var[A θ i (ψ 0 )|X, Z] = D 2,ii (ψ 0 ), and cov[A β i (ψ 0 ), A θ i (ψ 0 )|X, Z] = D 3,ii (ψ 0 ). Note also that V 1 , V 2 ,and V 3 are diagonal matrices, by independence between the individuals. It follows that var[C(ψ 0 )|X, Z] = D(ψ 0 ) and thus, var[U |X, Z] = I 2n and var[U ] = I 2n .

  surely to 0. Moreover, for n sufficiently large, there exists a positive constant c 4 < ∞ such that almost surely, Σ -I k converges almost surely to 0. The almost sure convergence to k follows by similar arguments.

Note:Figure 1 .Figure 2 .Figure 3 .

 123 Figure 1. Histograms and Q-Q plots for β2,n in model M1, with no immunes in the sample (the percentage of immunes is given in brackets). n is the sample size. All results are based on 1500 simulated datasets.

Figure 4 .

 4 Figure 4. Histograms and Q-Q plots for β2,n in model M1, with 75% of immunes.

Table 1 .

 1 . A score test for testing a zeroinflated Poisson regression model against zero-inflated negative binomial alternatives. Simulation results for models M 1 : β = (-.8, 1) and M 3 : β = (-.8, 0) : empirical size) of the Wald test at the level 5% for testing H 0 : β 2 = 0. For each percentage of immunes, the percentage of infected among the susceptibles is 30%. All results are based on 1500 replicates.

	Biometrics 57 (1), 219-223.
	Xiang, L., A. H. Lee, K. K. W. Yau, and G. J. McLachlan (2007). A score test for overdis-
	persion in zero-inflated Poisson mixed regression model. Statistics in Medicine 26 (7),
	1608-1622.

Note: n: sample size. (•): root mean square error. [•]: mean absolute error. * : empirical power ( †

Table 2 .

 2 Simulation results for models M 2 : β = (1, .7) and M 4 : β = (1, 0) : empirical power ( † : empirical size) of the Wald test at the level 5% for testing H 0 : β 2 = 0. For each percentage of immunes, the percentage of infected among the susceptibles is 70%.

				percentage of immunes in the sample
		0%			25%	50%			75%
	n	β 1,n	β 2,n	β 1,n	β 2,n	β 1,n	β 2,n	β 1,n	β 2,n
					Model M 2
	100	1.026	0.720	0.945	0.723	0.949	0.740	0.834	0.647
		(0.246) (0.273)	(0.780) (0.534)	(0.988) (0.788)	(1.549) (1.455)
		[0.196] [0.215]	[0.655] [0.376]	[0.829] [0.555]	[1.326] [0.933]
			0.746 *		0.132 *	0.118 *		0.088 *
	500	1.003	0.712	1.098	0.717	1.112	0.721	0.840	0.652
		(0.107) (0.115)	(0.651) (0.247)	(0.672) (0.279)	(0.969) (0.534)
		[0.086] [0.091]	[0.518] [0.202]	[0.534] [0.230]	[0.802] [0.421]
			1 *		0.503 *	0.418 *		0.168 *
	1000	1.003	0.707	1.078	0.711	1.096	0.719	0.842	0.657
		(0.071) (0.082)	(0.590) (0.215)	(0.571) (0.224)	(0.796) (0.439)
		[0.057] [0.065]	[0.428] [0.168]	[0.441] [0.181]	[0.670] [0.352]
			1 *		0.779 *	0.675 *		0.205 *
	1500	1.001	0.701	1.035	0.705	1.069	0.709	0.887	0.655
		(0.064) (0.065)	(0.450) (0.163)	(0.466) (0.177)	(0.604) (0.312)
		[0.050] [0.052]	[0.344] [0.135]	[0.358] [0.144]	[0.502] [0.257]
			1 *		0.986 *	0.926 *		0.300 *
					Model M 4
	100	1.030	0.001	1.110	0.007	1.154	0.017	0.913	-0.003
		(0.233) (0.234)	(0.852) (0.969)	(1.211) (1.347)	(1.775) (1.640)
		[0.182] [0.187]	[0.684] [0.587]	[0.995] [0.792]	[1.450] [0.865]
			0.058 †		0.072 †	0.083 †		0.066 †
	500	1.007	-0.005	1.105	0.020	1.123	0.054	0.915	-0.009
		(0.103) (0.103)	(0.609) (0.293)	(0.690) (0.318)	(0.817) (0.370)
		[0.081] [0.082]	[0.492] [0.180]	[0.562] [0.208]	[0.614] [0.215]
			0.046 †		0.050 †	0.063 †		0.051 †
	1000	1.003	0.000	1.091	-0.003	1.101	0.033	0.934	-0.003
		(0.071) (0.070)	(0.521) (0.198)	(0.578) (0.210)	(0.757) (0.256)
		[0.057] [0.055]	[0.437] [0.125]	[0.455] [0.135]	[0.600] [0.142]
			0.051 †		0.045 †	0.042 †		0.039 †
	1500	1.003	0.001	1.073	0.009	1.115	0.015	0.934	0.002
		(0.057) (0.057)	(0.480) (0.132)	(0.501) (0.139)	(0.633) (0.175)
		[0.046] [0.046]	[0.392] [0.087]	[0.400] [0.104]	[0.521] [0.109]
			0.042 †		0.040 †	0.046 †		0.047 †

Note: *

Table 3 .

 3 "Naive" analysis of model M 1 : β = (-.8, 1)

			percentage of immunes in the sample
		25%	50%		75%
	n	β 1,n	β 2,n	β 1,n	β 2,n	β 1,n	β 2,n
	100	-1.154	0.023	-1.632	0.017	-2.410	0.001
		( 0.428) (1.011)	(0.879) (1.025)	(1.776) (1.071)
		[0.365]	[0.977]	[0.833] [0.983]	[1.610] [1.003]
			0.049 *	0.057 *	0.052 *
	500	-1.128	0.087	-1.594	0.042	-2.305	0.002
		(0.344) (0.915)	(0.803) (0.963)	(1.513) (1.010)
		[0.328]	[0.913]	[0.794] [0.958]	[1.505] [0.997]
			0.049 *	0.051 *	0.053 *
	1000	-1.131	0.059	-1.590	0.050	-2.297	0.033
		(0.338) (0.941)	(0.795) (0.952)	(1.501) (0.970)
		[0.330]	[0.940]	[0.790] [0.950]	[1.497] [0.966]
			0.053 *	0.051 *	0.054 *
	1500	-1.127	0.050	-1.591	0.046	-2.302	0.039
		(0.332) (0.953)	(0.794) (0.955)	(1.504) (0.962)
		[0.327]	[0.952]	[0.791] [0.954]	[1.502] [0.960]
			0.051 *	0.050 *	0.053 *

Table 4 .

 4 Dengue fever data analysis

		naive analysis	model (1)-(2)-(3)
	parameter estimate	sd	estimate	sd
	β 1	1.552	0.255	7.654	1.485
	β 2	-0.055	0.007	-0.131	0.020
	β 3	-0.813	0.207	-4.501	1.059
	θ 1			0.497	0.159

Appendix

Proof of Lemma 1. Recall that I k denotes the identity matrix of order k. Then we write: