Daniela Florez

Antoine Stéphan
email: antoine.stephan@univ-lorraine.fr

Partrick Perré
email: patrick.perre@centralesupelec.fr

Romain Rémond
email: romain.remond@univ-lorraine.fr]

Probabilistic multi-objective optimization of wood torrefaction conditions using a validated mechanistic model

Keywords: Wood torrefaction, thermal modification, optimal conditions, treatment heterogeneity, overheating, computational model

This paper uses a comprehensive computational model to propose optimal wood torrefaction conditions by probabilistic optimization. Its main outcome is to propose tailor-made heat treatment conditions (temperature levels-duration of mild pyrolysis at temperature levels ranging from 200 to 300 °C) to meet users' expectations in terms of overall mass loss, duration and homogeneity of treatment. To this purpose, beech wood boards were torrefied with a usual 3-steps treatment schedule (drying, heating and cooling) under contrasting configurations in a well-instrumented device. The heterogeneity of the treatment within the wood sample was assessed through X-ray attenuation profile and water vapour sorption isotherm. These results allowed the model to be validated. In particular, it predicts the evolution of the mass loss and internal temperatures with good accuracy, including the temperature overshot. The results highlight the need to adjust the heat treatment schedule to each input parameter such as the wood piece dimensions and its initial moisture content or density, in order to limit the effect of exothermic reactions. The torrefaction model was then embedded in a probabilistic optimization process. A case study demonstrates the ability of the model to propose an alternative 3-steps treatment schedule able to reach the target mass loss while controlling the temperature overshot within the wood piece.

Introduction

Torrefaction of wood is a mild pyrolysis treatment at temperature levels ranging from 200 to 300 °C in the absence of oxygen. This induces the thermal alteration of the cell wall constituents of wood, namely hemicelluloses, cellulose and lignin, via a series of thermally activated chemical reactions coupled with heat and mass transfer phenomena. Chemical modification concerns mainly hemicelluloses degradation, lignins reticulation and modification of the crystalline structure of the cellulose [START_REF] Lv | TGA-FTIR Analysis of torrefaction of lignocellulosic components (cellulose, xylan, lignin) in isothermal conditions over a wide range of time durations[END_REF][START_REF] Qu | Experimental study of biomass pyrolysis based on three major components: hemicellulose, cellulose, and lignin[END_REF][START_REF] Rousset | Effect of severe thermal treatment on spruce and beech wood lignins[END_REF][START_REF] Shen | The pyrolytic degradation of wood-derived lignin from pulping process[END_REF][START_REF] Shen | The mechanism for thermal decomposition of cellulose and its main products[END_REF][START_REF] Wu | Low temperature pyrolysis characteristics of major components of biomass[END_REF][START_REF] Yang | Characteristics of hemicellulose, cellulose and lignin pyrolysis[END_REF].

Over the last 50 years, torrefaction has been used to effectively improve various wood properties such as durability and dimensional stability [START_REF] Bourgois | Characterization and analysis of torrefied wood[END_REF][START_REF] Stamm | Staybwood-Heat-stabilized wood[END_REF][START_REF] Viitaniemi | Decay-resistant wood created in a heating process -a heat-treatment process of wood development by VTT building Technology yields timber products with enhanced properties[END_REF]. These studies were intended for solid wood as a material.

Later, torrefaction was promoted for other purposes. Current perspectives of this technology aim at taking advantage of renewable resources, to valorize wood and biomass wastes, thus limiting the exploitation of natural resources. For example, in an attempt to replace fossil fuels, a pre-treatment by torrefaction of woody biomass is advised to obtain hydrophobic products with increased specific energy contents and improved grindability [START_REF] Almeida | Alterations in energy properties of eucalyptus wood and bark subjected to torrefaction: The potential of mass loss as a synthetic indicator[END_REF][START_REF] Ong | Variation of lignocellulosic biomass structure from torrefaction: A critical review[END_REF][START_REF] Pierre | Reduction of biomass resilience by torrefaction: apparent stiffness during failure (ASF) and specific failure energy (SFE) assessed by a custom impact device[END_REF][START_REF] Kota | Biomass torrefaction: An overview of process and technology assessment based on global readiness level[END_REF]. Besides, it also leads to increased anti-shrinkage efficiency, which reduces swelling during moisture recovery, which is an interesting feature for wood chips stored in a silo. Furthermore, torrefaction can be advantageous in processing heterogeneous cellulosic feedstock of different density, humidity content and size, such as those from forestry (logs, chips, sawdust, bark, etc.) or agricultural crops [START_REF] Kumar | Integrating torrefaction in the wood pellet industry: A critical review[END_REF][START_REF] Zulfiqar | Flow properties of biomass and coal blends[END_REF]. These characteristics make torrefaction suitable for integration into existing commercial regimes such as the white wood pellet value chain [START_REF] Kumar | Integrating torrefaction in the wood pellet industry: A critical review[END_REF][START_REF] Thrän | Moving torrefaction towards market introduction -Technical improvements and economic-environmental assessment along the overall torrefaction supply chain through the SECTOR project[END_REF].

Moreover, end-of-life wooden furniture may be recycled and decontaminated by removing the ureaformaldehyde and melamine-formaldehyde resins through a low temperature pyrolysis [START_REF] Girods | Thermal removal of nitrogen species from wood waste containing urea formaldehyde and melamine formaldehyde resins[END_REF][START_REF] Besserer | Cascading recycling of wood waste: A review[END_REF]. This opens the possibility to use recycled wooden furniture materials for energy production. Nonetheless, the introduction of the torrefaction technology in industrial manufacturing processes can pose serious safety issues with regards to fire and injury risks. In this respect, various authors have reported that overheating occurs within wood products induced by the exothermic reactions during torrefaction. Indeed, as wood has a low thermal conductivity, the internal source of heat due to exothermic reactions gives rise to an important temperature gradient needed to drive heat towards the exchange surfaces. The core temperature therefore increases with, and in turn, activates even more, the exothermic reactions [START_REF] Branca | Chemical characterization of volatile products of biomass pyrolysis under significant reaction-induced overheating[END_REF][START_REF] Blasi | On the experimental evidence of exothermicity in wood and biomass pyrolysis[END_REF][START_REF] Turner | An experimental and theoretical investigation of the thermal treatment of wood (Fagus sylvatica L.) in the range 200-260°C[END_REF]. A thermal runaway is then likely to occur at the scale of a boards stack or a packed bed of wood chips where a dual-scale interaction occurs. The wood piece supplies the energy to the airflow due to the exothermic reactions, the gas temperature increases along the stack or the packed bed which triggers even more the exothermic reactions along the flow direction [START_REF] Perré | A comprehensive dual-scale wood torrefaction model: Application to the analysis of thermal run-away in industrial heat treatment processes[END_REF]. The risk of thermal runaway depends on the size of the wood piece and process conditions (heating rate, airflow velocity, temperature, etc.).

In addition to safety issues, this produces also a heterogeneous treatment, both within and between the wood pieces, hence undesirable consequences on the mass yield and product quality. Nowadays, however, torrefaction programs are often empirical and customized according to the wood species (softwood, hardwood, density, etc.), the characteristics of the products to be torrefied (e.g., dimensions, geometry, initial moisture content, etc.) and the user expectations.

The motivation for this research arises from the difficulty of determining the optimal torrefaction schedule that could limit the overheating within the wood piece and ensure that products meet quality requirements, for example the hydrophobic character, energy content and grindability for energy purposes, or the anti-shrinkage efficiency, mechanical performance and color for solid wood. Since most of these properties are monotonic functions of the overall mass loss [START_REF] Almeida | Alterations in energy properties of eucalyptus wood and bark subjected to torrefaction: The potential of mass loss as a synthetic indicator[END_REF][START_REF] Almeida | Changes in wood-water relationship due to heat treatment assessed on micro-samples of three Eucalyptus species[END_REF][START_REF] Brischke | Quality control of thermally modified timber: Interrelationship between heat treatment intensities and CIE L*a*b* color data on homogenized wood samples[END_REF][START_REF] Candelier | Control of wood thermal treatment and its effects on decay resistance: a review[END_REF][START_REF] Candelier | Utilization of temperature kinetics as a method to predict treatment intensity and corresponding treated wood quality: Durability and mechanical properties of thermally modified wood[END_REF][START_REF] Candelier | Comparison of mechanical properties of heat treated beech wood cured under nitrogen or vacuum[END_REF], the latter could therefore be used as a target indicator to fine-tune the torrefaction parameters (duration-temperature levels) in order to meet quality specifications.

In the present work, the research approach consists of two complementary parts. In the first, in a series of experiments, beech (Fagus sylvatica) wood boards were torrefied using the usual 3-step program (drying, heating at high temperature and cooling) under different conditions in terms of the sample thickness and initial moisture content. The experimental data, including local values of internal sample temperatures and mass loss progression, were then compared with simulations of a torrefaction model [START_REF] Turner | An experimental and theoretical investigation of the thermal treatment of wood (Fagus sylvatica L.) in the range 200-260°C[END_REF][START_REF] Perré | A 3-D version of TransPore: a comprehensive heat and mass transfer computational model for simulating the drying of porous media[END_REF][START_REF] Rémond | Modeling the drying and heat treatment of lignocellulosic biomass: 2D effects due to the product anisotropy[END_REF][START_REF] Rousset | Choix et validation expérimentale d'un modèle de pyrolyse pour le bois traité par haute température : de la micro-particule au bois massif[END_REF]. The comparison between simulation and measurement in contrasting configurations allowed the robustness of the model to be assessed. In the second part, this model was integrated into an optimization algorithm with the purpose of proposing optimized heat treatment schedules (temperature levels and durations), meeting the user expectations in terms of mass loss, duration and homogeneity of treatment, while mitigating the undesirable overheating within the wood piece. A case study was conducted to experimentally validate one of the treatment schedules proposed by the optimization procedure.

2.

Measuring the torrefaction kinetics and treatment homogeneity under contrasting conditions

Materials

In order to focus on the effects of different factors such as the sample thickness and initial moisture content on the homogeneity of the heat treatment, the sampling was designed to limit the effect of the biological variability of wood (Fig. 1a). For that purpose, a green beech (Fagus sylvatica) log from the Vosges region in France, with a diameter of about 60 cm, was selected and sawn. Then, two flatsawn boards were selected at least 10 cm far from the bark to keep only the heartwood, and each board was cut into specimens of final dimensions 20 mm (R) x 150 mm (T) x 300 mm (L) and 40 mm (R) x 150 mm (T) x 280 mm (L). Note that the large volume/surface ratio of the solid wood samples used here are more demanding as they exacerbate the intra-particle temperature peak due to the importance of internal transfers. Such sizes also ensure the reliability and accuracy of the experimental measurements (mass evolution and internal temperature). Smaller volume/surface ratios encountered in energy applications could be easily predicted by this mechanistic model as these situations limit the coupling due to internal transfers. The boards were pre-conditioned at 12 % of moisture content on dry basis (MC) in a climatic chamber at 40 °C and 75 % of relative humidity.

Subsequently, the specimens were instrumented with thermocouples of 1 mm diameter with a gas-tight Swagelok connector to link the board to the probe. As illustrated in Figure 2, they were placed halflength and mid-width of the boards, 75 mm deep from the lateral surface (RL plan). One was inserted at the core (half thickness) and the second at about 2 mm from the surface (TL plan). Finally, the endpieces of the boards were coated with a silicone film maintained by two metallic parts strongly tied (see Fig. 2), to prevent mass transfer in the longitudinal direction and avoid any gas leakage due to overpressure during the treatment. This leads the boards to have a behavior similar to that of long boards with transfers mainly in the radial and tangential directions.

Methods

Experimental setup

The heat treatment experiments were carried out in an experimental reactor designed for this purpose (Fig. 2), inspired by Colin's thesis work [START_REF]Séchage en continu du bois énergie comme moyen de préconditionnement en vue de sa conservation thermochimique : approches expérimentale et numérique[END_REF]. It consists of an oven Memmert UF110 plus equipped with four polished stainless-steel heating walls (top, bottom, and sides) heated by electric resistances, two temperature probes placed in the middle and top of the oven chamber, and a tangential fan on the back wall. The maximum working temperature is 300 °C (setting accuracy of 0.5 °C) and the heating rate can be precisely tuned in the range of 1 to 10 °C min -1 . The fan ensures the circulation of the gases, and its velocity, adjustable in 10 % steps, was set at 40 % for all experiments. The complex air velocity inside the chamber is not known, but, as explained below, the heat transfer coefficient, the pertinent parameter, has been determined experimentally. Below the reactor, a digital balance is placed in a hermetic chamber, which communicates with the reactor chamber through a glass tube (inner diameter of 12 mm) embedded under the reactor bottom wall. Glass is chosen here because of its lower thermal conductivity compared to metal. The wood specimen is attached to a sample holder, placed in the reactor chamber, that stands directly on a tripod resting on the balance, both of which are clamped to the ends of a glass rod (8 mm diameter) passing through the glass tube. The system is continuously swept by a nitrogen flow (5 L min -1 ) at room temperature and fed into the weighing chamber. To prevent cooling of the reactor chamber by the flow of nitrogen, an aluminum plate (3 mm thick) was placed 3 cm above the bottom wall of the oven, thus creating a preheating chamber that retains and heats the entering gas before letting it flow through the main chamber. The volatile products are exhausted from the side of the reactor. The oxygen level is continuously measured at the outlet flow by using an oxygen analyzer with a zirconium oxide sensor (ZOA 100). The whole system was remotely controlled and monitored The overall external heat transfer coefficient of the reactor was determined by considering the contribution of the convective and radiative transfers. On the one hand, the convective heat transfer coefficient was assessed experimentally by the sponge method [START_REF] Challansonnex | Non-Fickian diffusion in biosourced materials: Experimental determination of the memory function using minute samples[END_REF] at a temperature of 50 °C, for which radiation is low, and under flow conditions similar to those used for the torrefaction experiments. A value of about 12 W m -2 K -1 was determined for this parameter. On the other hand, the contribution of the heat transfer by emitted radiation, from the oven surfaces to the boards, was estimated according to the Stefan-Boltzmann law. The radiative heat transfer at the torrefaction temperature (i.e., 230 °C) was estimated about 13 W m -2 K -1 . In this work, the overall external heat transfer coefficient of the torrefaction reactor was found to be of the order of 35 W m -2 K -1 .

Fig. 2.

Reactor device for wood torrefaction treatment and photo of one instrumented board.

Torrefaction conditions

The course of the torrefaction was monitored in real time through the following parameters:

temperatures of the reactor chamber (middle and top) and board (core and subsurface), oxygen level at the outlet flow (e.g., less than 2 %) and mass loss. The experimental conditions of the heat treatment (HT) schedules performed on the four single board specimens are detailed in Table 1. The usual 3-phase schedule was applied on these boards: (1) a drying period at 110 °C for at least 1 hour; (2) heating at a rate of 7 °C min -1 until the treatment temperature of 230 °C; then maintained at this level for 3 hours;

(3) and finally cooled down at a rate of 2 C min -1 until the temperature decreases to 150 °C in the reactor.

Samples "d" were oven-dried at 110 °C during 24 to 48 hours just before the torrefaction test. In total, four specimens were heat-treated. For the study of the effect of the thickness, two of the boards, namely HT20d and HT40d, were oven-dried at 103 °C to remove the residual water before torrefaction with the usual schedule. Then, the third sample, named HT40m and having an initial moisture content of 12 %, was directly heated at 230 °C by omitting the first drying step to study the effect of the residual water on the treatment. The fourth sample, denoted as HT40d-opt, was oven-dried at 103 °C and torrefied with an optimized treatment schedule proposed later in this work.

The overall mass loss (ML) due to thermal degradation should be determined as follows:

𝑀𝐿 = [ 𝑀 𝑜 -𝑀 𝑡 𝑀 𝑜 ] 100% ( 1
)
where Mo is the mass at the anhydrous state and Mt is the mass at a time t during thermal treatment.

However, in this work, the experimental measurements of mass loss do not distinguish between moisture evaporation and wood degradation. Thus, this mass loss is named apparent mass loss (MLapp), which is defined by the following equation:

𝑀𝐿 𝑎𝑝𝑝 = [ 𝑀 𝑖𝑛𝑖𝑡 -𝑀 𝑡 𝑀 𝑖𝑛𝑖𝑡 ] 100% (2)
where Minit is the mass of the board in either the moist or dry state, according to the initial conditions at the beginning of the experiment.

Assessing the final density profile and hygroscopicity

The density profiles of the heat-treated boards were scanned by using X-ray attenuation with a Dax5000 from Fagus-GreCon. Samples with dimensions of 50 mm x 50 mm were sawn from the geometric centre of the boards and subsequently dried at 103 °C until reaching constant mass before analysis (Fig. 1b).

Besides, the differences in water adsorption capacity between the core and the surface of the torrefied boards (Fig. 1b) were assessed by dynamic vapour sorption (DVS), using a DVS intrinsic apparatus from Surface Measurement Systems. Samples with 30 to 40 mg mass were subjected to an isothermal cycle at 23 °C of humidification by increasing the relative humidity of the surrounding atmosphere from 0 % up to 90 %, in steps of 15 %. The change in sample mass was followed at regular time intervals

(1 point/min). Equilibrium moisture content (EMC %) was considered as established when the change in the sample mass was less than 0.0005 %/min for a given relative humidity level, where this slope is evaluated over a time window of 10 min.

Optimization of torrefaction conditions

In this work, a program has been developed to determine tailor-made wood torrefaction conditions to meet the user expectations. The overall mass loss (ML) of the sample is a synthetic indicator of the treatment intensity [START_REF] Almeida | Changes in wood-water relationship due to heat treatment assessed on micro-samples of three Eucalyptus species[END_REF], and it can be used to determine some properties of the torrefied products such as the hygroscopicity, anti-swelling efficiency (ASE), color, grindability and energy content. Just as the final MC is the target of drying, the final ML can be defined as the target of wood torrefaction on the basis of the end user applications. The methodology implemented for the development of the program consists of three complementary parts:

• A torrefaction model computes the ML field over time and the temperature overshoot in the board for a given treatment schedule,

• An objective function is defined to describe the output targets that meet the user expectations,

• A metamodel calls the torrefaction model, starts with a given schedule and modifies it to minimize the objective function.

Each step is presented in detail in the following section.

Torrefaction model

The 2-D version of the computational model TransPore, developed in previous works [START_REF] Turner | An experimental and theoretical investigation of the thermal treatment of wood (Fagus sylvatica L.) in the range 200-260°C[END_REF][START_REF] Perré | A 3-D version of TransPore: a comprehensive heat and mass transfer computational model for simulating the drying of porous media[END_REF][START_REF] Rémond | Modeling the drying and heat treatment of lignocellulosic biomass: 2D effects due to the product anisotropy[END_REF][START_REF] Rousset | Choix et validation expérimentale d'un modèle de pyrolyse pour le bois traité par haute température : de la micro-particule au bois massif[END_REF], was used in this work. This computational model accounts for the thermal decomposition kinetics of the main components of the cell wall (cellulose, hemicelluloses and lignin) via a series of parallel, consecutive or competitive chemical reactions. Enthalpy values are associated with these reactions, which allows a two-way spatio-temporal coupling between heat and mass transfer and kinetics to be considered: activation of chemical reactions by the local temperature value and evolution of the temperature fields by the cumulative effect of the local heat source due to reactions. In this respect, the temperature computed by TransPore, for each control volume of the wood piece, is used to estimate the rate of the chemical reactions and their related heat sources (or sinks) and the production of volatiles. the literature [START_REF] Turner | An experimental and theoretical investigation of the thermal treatment of wood (Fagus sylvatica L.) in the range 200-260°C[END_REF][START_REF] Rémond | Modeling the drying and heat treatment of lignocellulosic biomass: 2D effects due to the product anisotropy[END_REF]. Finally, the board section is discretized by a set of 41 × 21 control volume, in the radial and tangential directions respectively.

Module for generating the torrefaction conditions to reach the target ML

A simple module written in Fortran has been developed to generate a wide variety of treatment schedules. This module is integrated in the TransPore model, in which the generated schedules are used as boundary conditions of the wood piece, which allows the ML to be computed versus time. In the following, the usual 3-step schedule is generated:

T max T 1 , t 1 R 2 R 1 R 1

Cooling

Heat treatment Drying T drying

Time → Temperature → Fig. 3. Example of the torrefaction schedule that can be proposed by the module. This includes a first drying step and two subsequent steps for wood torrefaction.

• Step 1 is a drying stage at a level ranging from 110 °C to 130 °C, maintained until the computed average MC reaches 1%.

• Step 2 is the heat treatment stage itself with two temperature plateaus, the intermediate temperature (T1) and the maximum torrefaction temperature (Tmax). The temperature change between plateaus is defined by the heating rate (R1) and the duration of the intermediate plateau is Δt1.

Step 2 ends as soon as the target ML is reached.

• Step 3 is the cooling stage, which lasts until the reactor temperature drops to 80 °C.

By defining the heating rate (R1), the parameters of the intermediate plateau (T1 , Δt1) and the maximal temperature (Tmax) as variables of the module, it is already possible to generate a wide variety of torrefaction schedules to reach the target ML (Fig. 3).

Objective function and variables

An objective function has been built up from three output targets defining the user expectations. The first can be defined by the heterogeneity of the ML field within the wood sample at the end of the treatment, since it is sought to perform a homogeneous treatment of the sample. The second concerns the treatment duration, as good homogeneity should not be sought at the expense of longer durations.

The third and last output target intends to limit the exothermic peak which is at the origin of the thermal runaway, hence heterogeneity, observed in industrial ovens. Accordingly, each of the objective variables can be defined as follows:

Firstly, the standard deviation of the ML at the end of the treatment, which is given as a fraction of its mean value 𝑀𝐿 ̅̅̅̅ ,

𝑍 1 = 𝜎/𝑀𝐿 ̅̅̅̅ (3)
with σ defined as:

𝜎 = √ 1 𝑛𝑐𝑣 ∑(𝑀𝐿(𝑗) -𝑀𝐿 ̅̅̅̅ ) 2 𝑛𝑐𝑣 𝑗=1 (4)
In the above expression, ncv is the total number of control volumes in the board section, and ML(j) is the local ML at the control volume j.

Secondly, the rate of torrefaction progression establishes the relationship between the duration of the treatment tend and the maximum duration ttarget accepted by the user:

𝑍 2 = max (0 ; 𝑡 𝑒𝑛𝑑 -𝑡 𝑡𝑎𝑟𝑔𝑒𝑡 𝑡 𝑡𝑎𝑟𝑔𝑒𝑡 ) (5)
Thirdly, the term considering the temperature overshoot, which can be considered as a safety factor to reduce the risk of the thermal runaway, is defined by:

𝑍 3 = max (0 ; 𝑇 𝑖𝑛 -𝑇 𝑠𝑝 𝑇 𝑠𝑝 ) (6)
where 𝑇 𝑖𝑛 = max (𝑇(𝑗) 𝑗=1,𝑛𝑐𝑣 ) during the treatment period, and Tsp is the gas temperature around the board surface.

The difficulty remains in finding the right compromise between a short treatment and a homogeneous treatment, while limiting the temperature overshoot. Finally, the global objective function is given as a weighting average of all sub-functions:

𝑔 = 𝑤 1 • 𝑍 1 + 𝑤 2 • 𝑍 2 + 𝑤 3 • 𝑍 3 (7)
where 𝑤 𝑖 is the weighting factor, such as 𝑤 1 + 𝑤 2 + 𝑤 3 = 1, whose value is chosen by the decision maker according to his own constraints.

Minimization of the objective function

For given torrefaction conditions, the TransPore code predicts the values of the objective variables (Z1, Z2 , Z3). A metamodel, developed here in the versatile language Python, calls TransPore code by using the F2PY package (TransPore being written in the very fast Fortran language) and estimates the value of the objective function (Eq. 7) [START_REF] Peterson | 2PY: a tool for connecting Fortran and Python programs[END_REF]. Simulated annealing (SA) has been adopted in the metamodel to minimize this function. We used the dual annealing function included in the SciPy Python Package. A compromise between the quality of the objective function minimum and the number of calls of such function was looked for to limit the time and resource consumption of the SA (Annex A). The program takes 5 minutes for the whole optimization procedure, with two Xenon W cores at 3.5 GHz. The parameters of the SA function were set at the default values [START_REF] Xiang | Generalized simulated annealing for global optimization: the GenSA package[END_REF]: the initial artificial temperature value was 5230, the value of the parameter for visiting distribution was 2.62, the value of the parameter for acceptance distribution was -5, the value of the ratio triggering the restart of the process was 2.10 -5 , and the local search strategy was deactivated.

Finally, by fixing a target ML, the torrefaction conditions defined by R1 , T1, Δt1, and Tmax can be tuned by the metamodel to find the best compromise between the three sub-criteria to minimize the objective function.

Results and discussion

Experimental and simulated results

Dynamics of torrefaction treatment

Figures 4a and 4b present the experimental results for the temporal evolution of the core and subsurface temperatures, as well as the MLapp for the HT20d and HT40d boards. Both cases were conducted with the usual 3-step schedule (i.e., drying, heating at 230 °C and cooling). It can be observed that at the end of the drying plateau, the core (blue triangle markers) and subsurface (red square markers) temperatures of both boards reach the dry bulb temperature, and the moisture is almost completely removed, as evidenced by the unchanging MLapp (grey circle markers). At this point, the remaining MC is almost zero, as shown in Table 2. Subsequently, during heating to the treatment plateau at 230 °C, the temperature in the core exceeds the one at the subsurface of the boards, and both exceed the reactor temperature (black cross markers) until reaching a maximum peak. As is well known, this phenomenon is attributed to the exothermicity of the chemical reactions activated by the rapid increase of the internal temperature of the boards. On the other hand, when the board thickness is doubled, the onset of the exothermic peak is delayed and its intensity is greater. This is due to a longer time for heat conduction towards the core of the board, and then, when exothermic reactions are activated, the thermal insulation ensured by the sample thickness limits the cooling of the product by the exchange surfaces. Nonetheless, it was also observed that the differences between the core and subsurface temperatures were the same, about 2 °C, for both specimens. This result is unexpected since a greater difference is likely to be observed by doubling the board thickness. This may be attributed to the internal location of the subsurface thermocouple, which was placed 6 mm under the surface of the thickest board, instead of 2 mm as for the thinnest one. Table 2 The subsequent comparison case, shown in Fig. 4c, concerns the evolution of the thermal treatment of 309 the 40 mm thick board with an initial MC at 12 %, which was directly heated at 230 °C for 200 min 310 (without previous drying). Due to the high initial MC, the internal temperature of the board is delayed 311 from the oven temperature, as much of the energy of the system is used to supply the latent heat of 312 evaporation required to remove the moisture. Thus, because of the massive and sudden water evaporation, the profile of the core temperature exhibits a plateau at 100 °C followed by a sharp increase up to 280 °C after 200 min of heating. Regarding the MLapp, as expected, it is higher for the moist board (30 %), but with a true value of ML of the order of 18 %, which is comparable to the dry board (as shown in Table 2), regardless of the differences in the initial MC. In addition to the results from the three torrefaction experiments presented here, additional experiments (from reference [START_REF] Turner | An experimental and theoretical investigation of the thermal treatment of wood (Fagus sylvatica L.) in the range 200-260°C[END_REF] and several unpublished data) using different temperature-levels, duration and heating rates, were tested for model validation. The model predictions for each board are compared to the respective experimental data and plotted as solid lines in Fig. 4. In all the presented cases, the actual temperature in the reactor chamber is used as a boundary condition in the computational model

TransPore. As can be seen, the model succeeds in predicting with good accuracy the internal board temperatures (core and subsurface) for the very contrasting experimental conditions of sample thickness and initial MC. Regarding the MLapp progression, the comparisons are correct for the thickest boards, but it is less good for the HT20d one, since the model seems to overestimate the MLapp with a more advanced degradation kinetics. For the three different experimental configurations, the occurrence of the exothermic peak is captured well by the model, but its intensity is under-estimated, notably for the 40 mm thick boards. In the case of the moist board (HT40m), the observed over-predictions on temperature evolution may be related to the local pressure at the internal position of the probes, which seems to be at atmospheric pressure with a plateau at 100°C despite the intense evaporation of the bound water, while the model predicts a local overpressure of about 2.5 bar. A possible cause would be a leakage around the tight connector that links the board to the probe.

Heterogeneity of treatment intensity across the sample thickness

The intensity of the different heat treatments was then experimentally compared by measuring the changes in the density profile and water sorption capacity of the beech boards before and after torrefaction. The ability of the Grecon X-ray beam to detect density changes after heat treatment was first verified using a series of homogeneously heat-treated thin beech specimens with different mass losses. It was found that the Grecon analyzer is capable of accurately measuring small changes in mean density, even as little as 10 kg/m 3 . Figure 5a compares the density profile of the HT40d specimen with the untreated reference, both of which have the same pattern of annual growth rings. Since the X-ray beam scans were performed in the radial direction, i.e., the X-ray passes through the sample in the tangential direction, the growth rings were highlighted, and the alternation between earlywood (EW, lower density) and latewood (LW, higher density) was distinguished by the lowest and the highest density values of each peak, respectively.

As can be seen, the case corresponding to the heat-treated specimen displays more pronounced and narrower peaks that are positioned at lower density values than the untreated reference, due to the degradation and shrinking of the cell wall during torrefaction. At the same time, the density reduction seems to be greater in earlywood than in latewood, as evidenced by the decrease of the respective peaks, notably those placed at the half-length of the board, where the highest temperature overshoot was recorded. These results are consistent with studies performed on different wood species [START_REF] Hamada | Thermal stability of Abies alba wood according to its radial position and forest management[END_REF][START_REF] Shchupakivskyy | Application of high-frequency densitometry to detect changes in early-and latewood density of oak (Quercus robur L.) due to thermal modification[END_REF]. The higher susceptibility of earlywood to thermal degradation can be attributed to a lower proportion of the cellulose content, since its S2 layer is thinner, thus increasing the proportion of the more thermally sensitive components, namely hemicelluloses and lignin, to heat degradation, while cellulose remains stable at the recorded temperatures. In addition, greater shrinkage of LW compared to EW due to cell wall component degradations could also explain this observation, since shrinkage increases density, hence counterbalancing the mass loss per unit volume.

These trends obviously cannot be captured by the model, which considers a homogeneous density in the board at the initial state (here 698 kg m -3 for the untreated 40 mm thick board). Fig. 5b shows that the model simulates a parabolic density profile at the end of the treatment, with a slightly lower density at the core than at the surface. For all other specimens, the changes in their density profile are not presented here because the pattern of the annual rings was very different from the untreated reference, making it difficult to draw conclusions about the extent of the local degradation. The mean density of all heat-treated specimens is shown in Table 3. Fig. 6 shows the changes in hygroscopicity after the heat treatment, as well as the differences in the adsorption isotherms between the core (C, solid marks) and the outer surface layer (S, empty marks) of the dry boards, namely HT20d and HT40d. In both cases, the EMC of the untreated beech decreases over the analyzed range of relative humidity (i.e., 0 % to 90 %). It is well known that heat treatment reduces the sorption isotherms. The measurement presented here perfectly confirms this: first, a slight reduction at the surface, in agreement with the density and temperature evolution; second, a larger reduction at the core due to the overshoot in temperature, which triggers kinetics; third, an effect of the thickness, which exacerbates the overshoot due to a thicker insulation.

Numerical analysis of the time-temperature couples for a single step treatment

From the above, it becomes evident that the intensity of the heat treatment is very sensitive to the wood properties and initial conditioning, hence the need to adapt the treatment schedule to each particular case. To further investigate the effect of board thickness on the time-temperature couple, when a single temperature step treatment is applied, a series of simulations were carried out for varying the board thickness, namely 10 mm, 18 mm, 34 mm and 54 mm (in the radial (R) direction of the wood), for a 150 mm wide board (in the tangential (T) direction) and for a target ML equal to 15 %. At the beginning of the treatment, the boards are assumed to be anhydrous and at a temperature of 100 °C. Subsequently, these are subjected to a single-step schedule at constant temperature, for an undefined period of time until the desired ML is reached. Fig. 7 shows the time-temperature couples proposed to reach a target ML of 15 % for each thickness (l, color of the circles), as well as the values (diameter of the circles) of two objective variables, Z1 (Fig. 7a) and Z3 (Fig. 7b). Remember that Z1 (Eq. 3) estimates the relative ML heterogeneity at the end of the thermal treatment and Z3 (Eq. 6) quantifies the temperature overshoot. These graphs highlight the effect of thermal activation of the exothermic reactions: a clear time-temperature equivalence is revealed by the time-temperature couples. For instance, for a 10 mm thick board, the treatment duration is 82 h at 180 °C compared to 0.2 h at 280 °C. In fact, at low temperature levels (e.g., below 200 °C), the duration of the heat treatment is independent of the board thickness: the time required to reach a ML of 15 % is much longer than the characteristic time of heat diffusion in the board thickness. Whatever the thickness, the temperature remains quasi-uniform and the overshoot is low (see Fig. 7a). In contrast, as the temperature rises above 220 °C, the residence time is of the same order of magnitude as the diffusion time constant of the wood piece. Whether during heating or in the exothermic phase, step temperature profiles are established during the process. Consequently, the treatment time and its heterogeneity depend on the board thickness. This is confirmed by the relative heterogeneity at the end of torrefaction, quantified by Z1 (Fig. 7a).

This criterion increases with temperature and, in a more intensive manner, with the sample thickness.

Note that the area where the temperature is maximum does not correspond necessarily to greater degradation rates since, as mentioned above, the degradation profile is a result of the history of the temperature profile in time (Fig. 8).

Finally, Fig. 7b represents the evolution of the temperature overshoot or Z3 as a function of temperature, where an inflection point around 220 °C to 240 °C is revealed for the thickest samples. In fact, given that the duration of the treatment is similar to or less than the diffusion constant time, the temperature increase at the core of the board is not high enough to activate the exothermic reactions. In this case, the temperature overshoot is attenuated and the treatment is larger at the surface than in the core of the board. It is seen here that the impact of the torrefaction conditions on the objective variables is quite subtle.

The empirical approach to find the adapted sequence of the temperature levels-duration couples is a very laborious task. Accordingly, the optimisation procedure proposed in this work makes sense to find the most suitable torrefaction conditions to meet the user expectation.

Case study: experimental validation of 4-step torrefaction program proposed by the model

In the torrefaction experiments presented above, it seems difficult to limit the magnitude of the exothermic peak with the usual 3-step program. In fact, a gradual temperature increase or the addition of an intermediate step in the treatment program causes the chemical reactions to proceed at a reasonable rate while limiting the thermal runaway. In an attempt to reduce these undesirable effects, which is such a difficult task, the optimization algorithm of the objective function was employed. More specifically, in this case study, the process was optimized only in terms of the duration of the treatment (Z 2 ) and reduction of the temperature overshoot (Z3) by fixing, in Eq. ( 7), w1=0, w2=0.3 and w3=0.7. This choice was made because these two variables can be easily evaluated during the process. In this sense, the optimization variables considered according to the 4-step treatment schedule, illustrated in Fig. 3, were T1 and Δt1, while the values of the other variables were fixed as: Tmax=230 °C; Tdrying=110 °C; R1= 5 °C min -1 ; ttarget =5 h (without the cooling step); for a target ML equal to 15 %. T1 can be tuned by the SA in the metamodel between 190 °C and 210 °C and Δt1 between 100 min and 240 min.

It is interesting to note that the optimal treatment schedule given by the SA method implemented in this work chooses the maximum temperature allowed for T1 to limit both the residence time at 230 °C and the duration of the treatment. Moreover, it uses the cooling phase to control the exothermic peak. 2 and3, respectively. The temperature overshoot at the core is now 8 °C on this optimized schedule, as against 41 °C for the twin board on the first schedule (compare HT40d-opt vs. HT40d, Table 2) for an identical 3), which is lower than the fixed target of 15 %.

From the above, the model proved its ability to propose optimized schedules able to fulfill the desired criteria, but also to correctly predict the experimental behavior for this schedule with 2 treatment plateaus, more demanding than the initial experiments.

Conclusions

An optimization program proposing tailor-made torrefaction schedules (temperature levels-duration)

has been derived to control the torrefaction process of wood. This opens up the possibility of controlling the ultimate product qualities according to the user expectations. In this regard, the main outcomes of this work can be summarised as follows:

• A torrefaction device was developed to thermally treat instrumented wood sample under contrasting conditions,

• Parameters such as sample thickness and initial moisture content influence the temperature overshoot and overall mass loss,

• The water sorption isotherms are accurate for evaluating the heterogeneity of this treatment,

• A mechanistic model was used to successfully simulate these experiments,

• The use of an optimization procedure allows the treatment schedule to be adapted while mitigating the temperature overshoot,

• The predictive potential of the model has been experimentally validated with this optimized schedule.

Discrepancies were evidenced in comparisons between measurements and simulations of the intensity of the exothermic peak and the magnitude of the mass loss. An adjustment of some parameters of the heat treatment model or the implementation of the new model proposed by [START_REF] Perré | A robust and frugal model of biomass pyrolysis in the range 100-800 °C: Inverse analysis of DAEM parameters, validation on static tests and determination of heats of reaction[END_REF] appears a promising development of this work. Furthermore, this optimization approach, validated at the scale of one board, Annex B: Formulation of the torrefaction model. Summarized from [START_REF] Turner | An experimental and theoretical investigation of the thermal treatment of wood (Fagus sylvatica L.) in the range 200-260°C[END_REF] The total behaviour of wood is obtained from the evolution of its principal components: Wood = % cellulose + % lignin + % hemicellulose

The following system of differential equations is used for the degradation of the wood components (the kinetics parameters are depicted in table A1).

-Hemicellulose [START_REF] Lv | TGA-FTIR Analysis of torrefaction of lignocellulosic components (cellulose, xylan, lignin) in isothermal conditions over a wide range of time durations[END_REF] where  H represents the mass of hemicellulose; k1 and k2 the reaction rates and  (2)

         = - =   = - +   =   =   1

Fig. 1 .

 1 Fig. 1. Overview of the experimental study: (a) sampling pre-treatment of the wood boards and (b) post-

For a better

  understanding of the intricate two-way coupling between the heat and mass transfer and chemical mechanisms, the formulation of the model at the board scale is summarized in Annex B, and it is presented in more details in the work of Rémond et al. (2010) and Turner et al. (2010) [22,30]. In this work, the initial MC, dimensions and wood density measured on the board samples from our experiments are the input data for the computational model. Additional model parameters come from

Fig. 4 .

 4 Fig. 4. Temporal evolution of the MLapp (grey circles and dot/dashed curves), core (blue triangles and 305

Fig. 5 .

 5 Fig. 5. Density profile across the thickness (radial direction on the RT plan) of the HT40d board (a)

Fig. 6 .

 6 Fig. 6. Effect of treatment on water adsorption capacity: untreated (black dots) and torrefied samples

Fig. 7 .

 7 Fig. 7. Comparison of the time-temperature couples and their respective values of: (a) Z1 and (b) Z3,

Fig. 8 .

 8 Fig. 8. (a) Field of the maximum temperature Tmax reached throughout the process for each location at

Fig. 9 .

 9 Fig. 9. Comparison between model predictions and experimental data of heat-treatment progression of

  4 h). The overall ML reached experimentally is 11.8 % (Table

Fig. A. 1 .

 1 Fig. A.1. Minimum value of the objective function called by the algorithm: N represents the number of

  coal or tar and volatile matters produced during reaction i = 1

Table 1 .

 1 Samples and conditions of the torrefaction experiments

	Specimen name	Board thickness	Initial MC	Heat treatment schedule
		(mm)	(%)
	HT20d	20	Oven-dry (d)	usual
	HT40d	40	Oven-dry (d)	usual

Table 2 .

 2 Experimental characteristics of the exothermic peaks (onset time, core and subsurface temperatures, T°), and overall mass loss ( MLapp and ML) of the heat-treated boards.

			Maximum of the exothermic peak
	Specimen name	MC (after drying step)	Time (after drying step)	Core_T°	Subsurface_T°	MLapp ML
		(%)	(min)	(°C)	(°C)	(%)	(%)
	HT20d	0	75	239	237	-	16.6
	HT40d	0.7	120	271	269	19.8	19.1
	HT40m	12.1 a	215	280	266	30.0	17.8
	HT40d-opt 1.5	245	238	238	13.3	11.8

Table 3 .

 3 Mean density of heat treated and untreated beech wood measured by tomography. Two specimens were analysed at each time.

	Specimen	Mean density	Standard deviation
		(kg m -3 )	(kg m -3 )
	Untreated, 20 mm	671.8	26.1
	HT20d	558.6	54.1
	Untreated, 40 mm	698.5	34.4
	HT40d	580.2	8.7
	HT40m	575.0	3.0
	HT40d_opt	604.1	7.0

Funding

This research did not receive any specific grant from funding agencies in the public, commercial, or notfor-profit sectors.

7.

could be developed in a double-scale modelling of a board stack [START_REF] Perré | A comprehensive dual-scale wood torrefaction model: Application to the analysis of thermal run-away in industrial heat treatment processes[END_REF] or a packed bed of chips to optimize the schedules and avoid the thermal runaway in the industrial heat treatment process.

6.

where  C is the mass of cellulose and k3 and k4 are the reaction rates for cellulose.

-Lignin

where  L is the mass of lignin and k5 is the reaction rate for lignin. In the above equations we model the reaction rates according to the Arrenhius Law The macroscopic conservation equations governing heat and mass transfer phenomena in porous media during drying [START_REF] Perré | A 3-D version of TransPore: a comprehensive heat and mass transfer computational model for simulating the drying of porous media[END_REF], [START_REF] Whitaker | Simultaneous Heat, Mass, and Momentum Transfer in Porous Media: A Theory of Drying[END_REF] were extended in [START_REF] Turner | An experimental and theoretical investigation of the thermal treatment of wood (Fagus sylvatica L.) in the range 200-260°C[END_REF] to account for the effects of chemical reactions (component degradation, heat source, heat sink and volatiles production). The control volume (CV) method is implemented to discretize the conservation laws. Thereafter, an efficient inexact Newton method is used to solve the complicated nonlinear system that describes the drying and heat treatment processes.


Konan Jean

Geoffroy Kouakou
email: geoffroy.kouakou@inphb.ci

Essoham Ali
email: essoham.ali@univ-ubs.fr

On IPW-based estimation for zero-inflated Modified Borel Tanner regression with missing covariates

Keywords: Count data, Asymptotic properties, Missing data, Inverse probability weighting, Mathematics Subject Classification: 62-08, 62F12, 62D10, 62P10

In this study, we consider the zero-inflated Modified Borel-Tanner (ZIMBT) model regression. This model has been widely used to study the effects of covariates in count data sets that have many zeros. Since, in practice, some of the covariates involved in regression modelling often have missing values, we propose methods for estimating the parameters of the ZIMBT regression model with missing random covariates. Assuming that the covariates can be of any nature, i.e. discrete, continuous or categorical, we estimate the unknown selection probability using a parametric method and then non-parametrically (using a kernel estimator). We propose two inverse probability weighting (IPW) methods for estimating the parameters of the ZIMBT regression model with missing random covariates. Under certain regularity conditions, the asymptotic properties of the proposed estimators are studied. Also, the performances of the proposed estimators are evaluated by Monte Carlo simulations. The practical use of the proposed methodology is illustrated with The RAND Health Insurance Experiment Data.

Introduction

In statistics, count data is a type of statistical data that takes non-negative integer values and whose values are derived from a count process. Count data are often modeled by a Poisson distribution. The Poisson regression model has been widely used as a benchmark for the analysis of discrete counting models. But one of the main drawbacks in practice is that count data often exhibit overdispersion, so the Poisson regression model is not well suited to the data. Modelling methods adapted to this type of data have been widely explored in the literature. Generalized linear models provide a powerful framework for the analysis of such data. As a result, [START_REF] Lambert | Zero-inflated Poisson regression, with an application to defects in manufacturing[END_REF] popularized the ZI Poisson (ZIP) among all ZI models.

In addition, other discrete probability distributions, which describe count phenomena, have been proposed in the literature. In this paper, we work on a regression model based on the use of the modified Borel-Tanner distribution. The Modified Borel Tanner family of distributions has probability mass function (PMF) of the form: f MBT (X; α) := P(X = x) = Γ(2x + 1) Γ(x + 2)Γ(x + 1)

α x (1 + α) 2x+1 , , x = 0, 1, • • •
where α ∈ (0, 1), and the notation used is X ∼ MBT(α). This distribution can also be re-written as

P(X = x) = C x α x (1 + α) 2x+1 , x = 0, 1, • • • (1.1)
where C x is defined by

C x = 1 x + 1 2x x
are the Catalan numbers denoted by C n (see [START_REF] Olver | NIST Handbook of Mathematical Functions Hardback and CD-ROM[END_REF]). The modified Borel Tanner distribution exhibits a number of interesting properties. For a full description of the MBT model, readers may refer to [START_REF] Gòmez-Déniz | The Modified Borel-Tanner (MBT) Regression Model[END_REF]. Obviously, the properties available to the MBT distribution give it notable scientific importance, and therefore open up a large number of new avenues of research. For example, we can mention the work of [START_REF] Ghosh | Some Structural Properties Related to the Borel-Taner Distribution and its Application[END_REF] who proposed to study certain structural properties that had not been explored before, notably the recursive relationship between successive probabilities. The bivariate Borel Tanner model developed by [START_REF] Gupta | A bivariate Borel-Tanner distribution and its approximation[END_REF] would prove as useful as a univariate distribution.

The Zero-inflated Modified Borel-Tanner (ZIMBT) regression model was proposed by Anwar [START_REF] Hassan | Zero-inflated Modified Borel-Tanner Regression Model for Count Data[END_REF] to account for excess zeros. Since then, the model has been applied in various fields, such as infected blood cells (IBCs) per mm 2 on microscope slides prepared from randomly selected individuals. In these applications, maximum likelihood estimation (MLE) is used to derive parameter estimates.

Quite often, observations are incomplete, i.e. the values of one or more variables are missing. These are known as missing data observations. For example, in the case of a survey, some people do not answer certain questions, or some answers are forgotten during data entry. Sometimes, it's simply a problem of data coding (for example, the presence of characters that cannot be interpreted by the data-reading function) that causes certain values to disappear.

When observations with missing data are part of the observations used to build a model, a simple solution is to ignore these incomplete observations. Although frequently used by default, this solution proves simplistic, as we shall see later, because its consequences can be not only a reduction in model performance, but potentially a strong modeling bias that can lead to an inoperative model. Before deciding to ignore incomplete observations, it is therefore necessary to better understand why data is missing. A solution that often proves better (from the point of view of the performance of the decision model built from the data) is to estimate (or impute) the missing data and treat the estimated values as measured values. Finally, the ideal solution would be to correct the errors responsible for the missing data, or to complete the data collection, but this solution is not always accessible. If the model has already been built, and it is the observations to which the model is to be applied that are missing data, it is necessary to examine carefully the impact on the decision of different estimation methods.

There are three commonly recognized mechanisms for missing data: missing completely at random (MCAR), missing at random (MAR), and missing not at random (MNAR) [START_REF] Rubin | Inference and missing data[END_REF]. In recent years, estimation of regression models with missing covariates has been the subject of a large literature on estimation in regression models with missing covariates, see [START_REF] Chen | Model selection for zero-inflated regression with missing covariates[END_REF], [START_REF] Lukusa | Semiparametric estimation of a zero-inflated Poisson regression model with missing covariates[END_REF], [START_REF] Diallo | Estimation in zero-inflated binomial regression with missing covariates[END_REF], [START_REF] Jean | Estimation of zero-inflated bivariate Poisson regression with missing covariates[END_REF].

Although many works provide a thorough review of ZI models with missing data, to our knowledge, no previous research has studied how to estimate the parameters of a ZIMBT regression model when covariates are missing. Motivated by this work, we study estimation in ZIMBT regression when certain covariate values are missing for some of the individuals in the sample with application to a dataset from the RAND health insurance experiment.

In this work, we propose inverse probability weighting (IPW) methods to estimate the parameters of a ZIMBT regression model with MAR covariates. We establish the consistency and asymptotic normality of these estimates.

The progression of this work is organized as follows. In section 2, we describe the data, model and notations that will be used throughout the paper. In section 3, we present our estimates, regularity conditions and large sample properties. All proofs are reported in the appendix. In section 4, we conduct a simulation study to evaluate the finite-sample performance of the proposed estimates. The applicability of the proposed methods is demonstrated in Section 5. Section 6 presents the discussion.

Data and model

The pmf of the response variable is given by

P(Y = y) =      π + (1 -π) 1+λ 1+2λ if y = 0; (1 -π)C y λ y (1+λ) 1+y (1+2λ) 1+2y if y > 0.
(2.2) where P(Y = y) corresponds to the expression (1.1), λ > 0 and 0 ≤ p ≤ 1. We will say that the response variable follows a Zero inflated modified Borel-Tanner regression model and we will denote it as Y ∼ ZIMBT(λ, π). From this, it is straightforward to see that the mean and variance of Y i are given by

E(Y ) = (1 -π)λ and var(Y ) = (1 -π)λ [λ(1 + λ)(1 + 2λ) + πλ]
In ZIMBT regression, the mixing probabilities π i is usually modeled by a logistic regression : logit(π i (γ) = γ ⊤ W i and λ i is classically modeled as λ i = exp(β ⊤ X i ), where

X i = (1, X i2 , • • • , X ip ) ⊤ and W i = (1, W i2 , • • • , W iq ) ⊤ are
random vectors of predictors or covariates (both categorical and continuous covariates are allowed) and ⊤ denotes the transpose operator. Here, β and γ are respectively p and q-dimensional vectors of unknown regression parameters to be estimated. Let {(Y i , X i , W i ), i = 1, • • • , n} be a sample of independent observations and ψ = (β ⊤ , γ ⊤ ) ⊤ denote the whole unknown k-dimensional (k := p + q) parameter. The loglikelihood function ℓℓ n based on the observed sample is :

ℓℓ n (ψ) = n i=1 J i log e γ ⊤ W i + 1 + e β ⊤ X i 1 + 2e β ⊤ X i -log 1 + e γ ⊤ W i +(1 -J i ) y i βX i + (y i + 1) log(1 + e β ⊤ X i ) -(2y i + 1) log(1 + 2e β ⊤ X i ) , = n i=1 ℓ i (ψ).
(2.3)

where

J i = 1 {Y i =0} . The maximum likelihood estimator (MLE) ψ F,n = ( β ⊤ , γ ⊤ ) ⊤ of ψ is obtained by solving the score equation U F,n (ψ) = 0, where U F,n (ψ) = 1 √ n ∂ℓℓ n (ψ) ∂ψ = 1 √ n n i=1 ∂ℓ i (ψ) ∂ψ .
(2.4)

This estimating equation can be solved by using the expectation-maximization (EM) algorithm [START_REF] Hall | Zero-inflated Poisson and binomial regression with random effects: a case study[END_REF] or by a direct maximization of ℓℓ n . The next section describes the problem and the proposed estimator.

The proposed estimates

Problem formulation, notation and target parameter

Our first estimate is based on the inverse probability weighted parametric method idea. Let X obs and X mis denote the observed and missing components of X for a given subject, W obs and W mis denote the observed and missing components of W for a given subject. Let D = (Y, X (obs),⊤ , W (obs),⊤ ) ⊤ denote the vector of variables that are always observed on each individual. Let δ be the variable of observation indicators whose elements are zero or one depending on whether the corresponding element of X, W is missing or observed. We assume throughout that the data X (mis) , W (mis) are missing at random (MAR).

The IPW method was proposed by [START_REF] Horvitz | A generalization of sampling without replacement from a finite universe[END_REF] to estimate regression models with monotonic missing covariates. The IPW consists of inverse weighting of the observed data with the weight P(

δ i = 1|D i , X (mis) i , W (mis) i , W (mis) i ) = P(δ i = 1|D i ).
(3.5) Therefore, IPW estimator of ψ in model (2.3) is the solution of the weighted equation

1 √ n n i=1 δ i P(δ i = 1|D i ) Ψ i (ψ) = 0. (3.6)
Based parametric and nonparametric estimators of π(O i )'s, we propose two estimators of ψ.

In the following, two IPW approaches are proposed according to whether the selection probability π(D) is estimated by a parametric or non-parametric method.

Inverse Probability Weighted Parametric

Assuming that the selection probability π(D i ) can be specified by a parametric model a(ω, D i ), where ω is a vector of unknown regression parameters with true value ω 0 . Let ω n be the maximum likelihood estimates of ω 0 . Under some regularity conditions, which we define in the following section, ω n can be obtained as

ω n = arg max ω n i=1 a(ω, D i ) δ i 1 -a(ω, D i ) 1-δ i .
By straightforward calculations, it can be shown that

√ n(ω n -ω 0 ) = n -1/2 n i=1 δ i -a(ω 0 , D i ) a(ω 0 , D i ) 1 -a(ω 0 , D i ) Ω -1 (ω 0 ) ȧ(ω, D i ) + o P (1), where ȧ(ω, D i ) = ∂a(ω, D i ) ∂ω and Ω(ω) = E ȧ(ω, D) ȧ(ω, D) ⊤ a(ω, D)(1 -a(ω, D)) .
Having estimated ω, we can now derive the estimator ψ W,n of ψ by solving the following equation

U W,n (ψ, ωn ) = 1 √ n n i=1 δ i a(ω n , D i ) Ψ i (ψ) = 0. (3.7)
This first approach yields consistent estimates of the regression coefficients in the model, provided that π is correctly specified. However, misspecification may lead to biased estimated regression coefficients. To overcome this difficulty, we consider nonparametric kernel smoothers for the selection probabilities. Hence, we consider a nonparametric model for the selection probability.

Regularity conditions and large sample properties

The following regularity conditions are needed to establish the large sample properties of the estimators ψ W,n .

A.1 The true ψ 0 , ω 0 belongs to the interior of some known compact set respectively N , B of R d and R q .

A.2 The function a(ω, D) is strictly greater than 0 for all value in the support of D and all ω ∈ K.

A.3

The function a(ω, D) is differentiable with respect to ω. For every ω, ω ∈ B, a(ω, D) -a( ω, D) ≤ κ(D) ω -ω for some bounded function κ with E[κ(D)] = ρ.

A.4 Ψ 1 (ψ) has a finite second moment definite in a neighborhood of the true ψ. Moreover,

E Ψ ⊗2 1 (ψ)
is positive definite in a neighborhood of the true ψ. Here for any column vector a, a ⊗2 = aa ⊤ .

A.5 In a neighborhood of the true ψ, the first and second derivatives of U W,n (ψ, ω) with respect to ψ are uniformly bounded above by a function of (Y, X, W), whose expectations exist. Furthermore, -1

√ n ∂U W,n (ψ,ω) ∂ψ ⊤
converges to some a positive definite matrix Λ(ψ, ω) as n tends to infinity, in a neighborhood of the true ψ.

We first state the large sample properties of the inverse probability weighting parametric estimate Theorem 3.1. Assume that regularity assumptions (A.1)-(A.5) hold, it follows that ψ W,n is a consistent estimator of ψ, and √ n( ψW,n -ψ 0 ) is asymptotically normally distributed with mean zero and covariance matrix ∆ W , where

∆ W := Λ(ψ 0 ) -1 Φ(ψ 0 , ω 0 ) -Γ(ψ 0 , ω 0 )Ω(ω 0 ) -1 Γ(ψ 0 , ω 0 ) ⊤ [Λ(ψ 0 ) -1 ] ⊤ , with Λ(ψ) = E [Ψ(ψ)] ⊗2 , Φ(ψ, ω) = E δ a(ω, D) 2 Ψ(ψ)Ψ(ψ) ⊤ , Γ(ψ, ω) = -E δ a(ω, D) 2 Ψ(ψ) ȧ(ω, D) ⊤ .
In practice, we use a consistent estimator of the covariance matrix ∆ W . A consistent estimator of ∆ W is given by

∆ W,n := Λ n ( ψ n , ω n ) -1 Φ n ( ψ n , ω n )-Γ n ( ψ n , ω n )Ω n ( ψ n , ω n ) -1 Γ n ( ψ n , ω n ) ⊤ Λ n ( ψ n , ω n ) -1 ⊤ , where Λ n (ψ, ω) = - 1 n n i=1 δ i a(ω, D i ) ∂ 2 ℓ i (ψ) ∂ψ∂ψ ⊤ , Φ n (ψ, ω) = 1 n n i=1 δ i a(ω, D i ) 2 Ψ i (ψ) Ψ i (ψ) ⊤ , Γ n (ψ, ω) = - 1 n n i=1 δ i a(ω, D i ) 2 Ψ i (ψ) ȧ(ω, D i ) ⊤ , Ω n (ψ, ω) = 1 n n i=1 ȧ(ω, D i ) ȧ(ω, D i ) ⊤ a(ω, D i )(1 -a(ω, D i ))
.

Inverse Probability Weighted Semiparametric

Our second estimate is based on the multiple inverse probability weighting semi parametric approach. It is assumed that D is a vector containing a mix of discrete, categorical and continuous variables. Let D = D b , D c , where D c ∈ R d is a d-dimensional continuous random vector, and where D b is an b-dimensional discrete and categorical random vector. Let Z := X (mis) , W (mis) denote the covariate vector MAR. For simplicity's sake, let Z be a univariate covariate. Assume that δ is a dummy variable that is 1 when Z is observed; 0 otherwise.

Let K be a rth-order (r > d) kernel function of d variables with a finite support, and let K h (.) = K(./h), where h a bandwidth that satisfies the conditions we give in the following. A kernel estimator of π(D) is given by

π(o) = π(o b , o c ) = n k=1 δ k K h (D b k = o b , D c k -o c ) n j=1 K h (D b j = o b , D c j -o c )
.

Following [START_REF] Wang | Weighted semiparametric estimation in regression analysis with missing covariate data[END_REF], [START_REF] Lukusa | Semiparametric estimation of a zero-inflated Poisson regression model with missing covariates[END_REF]

, | π(D) -π(D)| = η n , where η n = {nh 2r + 1/(nh 2d )} 1/2 . Let f (o) denotes the probability density function of D at D = o. Also, let f (o) = 1 h d n k=1 K h (D b k = o b , D c k -o c ) an estimator of f (o).
Hence, the semiparametric IPW estimator of ψ is the solution of the weighted equation

U W s,n (ψ, π) = 1 √ n n i=1 δ i π(D i ) Ψ i (ψ) = 0.
(3.8)

Regularity conditions and large sample properties

The following regularity conditions are needed to establish the large sample properties of the estimators ψ Ws . B.2 K(u) is a kernel function of order r, that is, K(u)du = 1, u m K(u)du = 0 for m = 1, ..., (r-1), u r K(u)du ̸ = 0 and K 2 (u)du < ∞. Moreover, h satisfies nh 2r → 0 and nh 2d → ∞, as n → ∞.

B.3

The density function f (o) of D is bounded away from zero, and has r continuous and bounded partial derivatives with respect to the continuous components of (D b , D c ) almost everywhere.

B.4

In a neighborhood of the true ψ, the first and second derivatives of Ψ(ψ) with respect to ψ are uniformly bounded above by a function of (Y, X, W), whose expectations exist.

B.5 E 1 π(D) [Ψ 1 (ψ)] ⊗2 is
Σ W s := Σ 1 (ψ) -1 Σ 2 (ψ, π) -Σ * 2 (ψ, π) -Σ * 0 (ψ, π) [Σ 1 (ψ) -1 ] ⊤ , with Σ 1 (ψ) = E [Ψ(ψ)] ⊗2 , Ψ * 1 (ψ) = E Ψ 1 (ψ)|D b 1 , D c 1 , Σ * 0 (ψ, π) = E [Ψ * 1 (ψ)] ⊗2 Σ 2 (ψ, π) = E 1 π(D b , D c ) [Ψ 1 (ψ)] ⊗2 , Σ * 2 (ψ, π) = E 1 π(D b , D c ) [Ψ * 1 (ψ)] ⊗2
A consistent estimator of the asymptotic variance Σ W s is defined by

Σ W s,n := Σ 1,n ( ψ W s , π) -1 Σ 2,n ( ψ W s , π) -Σ * 2,n ( ψ W s , π) -Σ * 0,n ( ψ W s , π) [Σ 1,n ( ψ W s , π) -1 ] ⊤ , where Σ 1,n (ψ, π) = 1 n n i=1 δ i π(D b i , D c i ) Ψ i (ψ) ⊗2 , Σ * 2,n (ψ, π) = 1 n n i=1 δ i π 2 (D b i , D c i ) Ψ * i (ψ) ⊗2 , Σ * 0,n (ψ, π) = 1 n n i=1 Ψ * i (ψ) ⊗2 , with Ψ * i (ψ) = n k=1 δ k Ψ k (ψ)K h (D b k = D b i , D c k -D c i n ℓ=1 δ ℓ K h (D b ℓ = D b i , D c ℓ -D c i .
The proofs of Theorems 3.1-3.2 are given in Appendices B.

Numerical study

In this section, we carry out a simulation experiment to assess and compare the finite sample performance of the proposed estimates. The simulation setup is as follows. The outcome Y is generated from the ZIMBT regression model (2.2) defined by :

log(µ i ) = β 1 + β 2 X i2 + β 3 X i3 + β 4 X i4 + β 5 X i5 and logit(π i ) = γ 1 + γ 2 W i2 + γ 3 W i3 + γ 4 W i4 .
Covariates X i2 to X i5 are generated from the standard normal distribution, uniform distribution on [-1, 1], exponential distribution with rate 1 , Bernoulli distribution with probability 0.3. W i3 and W i4 be independently drawn from normal distribution with mean 1 and standard deviation 1.5 and Bernoulli distribution with probability 0.8 respectively.respectively. The linear predictors in logit(π i ) and log(µ i ) are allowed to share some common terms by letting W i2 = X i2 . Here, we take β = (-0.6, 0.75, -1, 0.5, -0.25) ⊤ . We consider successively two values for γ, namely : γ = (0.2, -0.5, 0.8, -1.5) ⊤ , γ = (0.25, -0.4, 0.7, 0) ⊤ . With these values, the average proportion c of zero-inflated data in the simulated data sets is 0.20, and 0.40 respectively. Under the MAR assumption, missingness indicators δ are simulated from a logistic regression model with selection probability P(δ = 1|D) = logit -1 (ω ⊤ D).

To establish a benchmark for comparisons, we also include an estimator based on the full data (FD) and the complete cases (CC) estimate, on complete cases only. Now, we compare the finite-sample performances of the following estimators: ψ F D : ML estimator that is the solution of the estimating equations in Equation (2.4) For each experimental configuration and each estimation method, the empirical bias, standard error (ASE), standard deviation (SD) and Root mean square deviation (RMSE) were computed . All calculations are made using the statistical software R (R Core Team, 2018).

ψ CC : the CC estimator that is the solution of U CC,n (ψ) = 1 n n i=1 δ i Ψ i (ψ) = 0 ψ W,p : IPW estimator
The results are presented in Tables 1 and2 . As expected, the performance of the estimators improves with increasing sample size. Parameter bias for all four estimators, whatever the sample size, is low. But the lowest bias is generally obtained by the FD estimate (i.e. in the ideal situation where there is no missing data). It should be noted that for the parameter γ the bias is not stable whatever the proportion of zero-inflation and missing data. These observations illustrate the general fact that accurate estimation in a zero-inflation regression model requires a balance between susceptible and non-susceptible subpopulations (i.e., a sufficient number of zero-and non-zero observations must be available to accurately estimate zero-inflation probabilities and the count submodel).

The standard deviation is the smallest when there are no missing data and the largest for the CC estimate. The IPW parametric and IPW semi-parametric estimates have standard deviations of similar magnitude. The semi-parametric IPW estimate has slightly larger variability. Finally, the IPW parametric and semi-parametric IPW estimate have smaller RMSE than the CC estimate in every scenario, and its RMSE is usually close to that of the FD estimate (see Table 2).

In summary, the IPW parametric and semi-parametric IPW estimation methods had comparable performances in terms of SDs, SEs and RMSEs. Additionally, for each estimation method, except for the ML estimation method, the agreement between SDs and SEs tended to worsen as the missing rate increased. Moreover, the SDs and SEs of all these estimators were decreased as the sample size was increased.

A Practical Application

RAND Health Insurance Experiment Data

The RAND Health Insurance Experiment (HIE), the most important health insurance study ever conducted, addressed two key questions in health care financing: How much more medical care will people use if it is provided free of charge? What are the consequences for their health? For illustration we apply the proposed estimates to a data set data from the RAND Health Insurance Experiment previously used by [START_REF] Cameron | Microeconometrics: Methodsand Applications[END_REF]. They estimated a more complete set of models and carried out a deeper analysis of the data than is possible or desirable here. The experiment, conducted by the RAND Corporation from 1974 to 1982, has been the longest running and largest controlled social experiment in medical care research. The main goal of the experiment was to assess how the patient's use of health services is affected by types of randomly assigned health insurance, including both fee-for-service and health maintenance organizations (HMOs). In the experiment the data were collected from about 8,000 enrollees in 2,823 families, from six sites across the country. Each family was enrolled in one of 14 different health insurance plans for either three or five years. The plans ranged from free care to 95% coinsurance below a maximum dollar expenditure (MDE), and also included assignment in a prepaid group practice. The data file consists of utilization, expenditures, demographic characteristics and health status.The final sample consists of 20,190 observations; each observation represents data for an experimental subject in a given year.

In this illustration, the response variable analyzed is the number of contacts with a Medical Doctor (MDU). Note that the ratio of zeros in the data set is 6308/20190=0.31 (see Figure ( 1)) suggests that the response is zero-inflated. Other variables were collected for each individual in the sample, such as (i) socio-economic variables: age (in years, age), education level (number of years of education, denoted by school), family income (in ten thousand dollars, denoted by income); (ii) indicating health status, ghindx represents an index of the general state of health of the people surveyed. This index is constructed from various physical and mental health indicators, and provides a synthetic measure of participants overall health.

Among the data from the RAND experiment on health insurance, the covariates ghindx and school contain missing data. The rate of missing data for ghindx is 25.87%. This high proportion of non-responses may be due to the fact that this is the most sensitive and confidential information that respondents did not wish to declare. Furthermore, the level of education of four individuals in the sample is missing. In total, 25.89% of the individuals in our sample have at least one missing item.

Table 4 presents the results of the analysis using three different estimation methods: CC, parametric IPW and semi-parametric IPW. The model parameter estimates obtained from these methods are very similar. The reasons for this are undoubtedly the large sample size and the average proportion of missing data. It can therefore be said that the proposed methods are good alternatives to CC analysis. All methods conclude that the estimates of age, income and education are statistically significantly different from zero and positive, indicating that older respondents are more likely to consult a doctor than younger ones.

In addition, the variables of income and school were identified as the most decisive factors in the decision never to consult an outpatient doctor. The probability of never consulting an outpatient doctor decreased with the number of years of education. In particular, education can make individuals more informed consumers of medical care services, so they may tend to diversify their use of medical care services. For example, focusing on the income variable, our point estimate at indicates that spending under the 95% coinsurance plan is lower than under the free care plan. We can reject the null hypothesis that income in the positive cost-sharing plans is equal to that in the free care plan.

Finally, it should be pointed out that we can continue to reject the nullity of non-response on the ghindx variable, which provides a synthetic measure of respondents' general health status. This important and statistically significant variable is (for us) one of the most surprising results of the RAND experiment. This indicates that the response of hospitalized patients on various indicators of physical and mental health brief health status in general is not robust to plausible adjustments for non-participation bias, and so the RAND data do not necessarily reject (though neither do they confirm) the hypothesis of lack of responsiveness on inpatient health status.

Discussion

In this work, we studied a problem of zero-inflated Modified Borel Tanner regression model with some covariates MAR. We proposed a parametric weighting estimator and a semiparametric kernel weighting estimator by estimating of the selection probability via a regression binary model and the kernel estimator. Then, we derived their limit theories and study their performances by a simulation study and a real example. In both situations, the proposed methods outperformed the true weight and CC estimator. An application of the proposed methods on a dataset from the RAND experiment on health insurance confirms the robustness of the proposed methods in the presence of missing data.In article, we consider monotone missing data. However, in many situations, the missing data pattern is not monotone (i.e. covariates are arbitrarily missing between individuals). The IPWs are mainly limited to monotone missing data because estimating the missingness mechanism in the case of non-monotone missing data is difficult. The adaptation of the IPW approach to non-monotone missing data in zero-inflated regression is a non-trivial task that also deserves further research. The approach proposed by [START_REF] Sun | On Inverse Probability Weighting for nonmonotone missing at random data[END_REF]) for estimating when the absence mechanism in the case of non-monotonic missing data could be useful.

The code below determines the initial estimates of β 1 , • • • , β 5 , γ 1 , • • • , γ 4 and calculates the MLE mleFD=maxLik (logLik=loglikfunFD,start=c(rep(0,length(phi))))

Estimates, standard errors and several other summaries can be obtained using the R function summary.

Appendix B : proofs of asymptotic results

Proof of Theorem 3.1. The consistency of ψ W,n , can be proved by checking that the conditions of Foutz's inverse function theorem [START_REF] Foutz | On the unique consistent solution to the likelihood equations[END_REF] are satisfied. Firstly, let us show that

∂U W,n (ψ,ωn) ∂ψ ⊤
exists and is continuous in an open neighborhood of ψ 0 . To justify this, it can be noted that U W,n (ψ, ωn ) is twice differentiable with respect to ψ and its second derivatives are continuous.

Secondly, let us show that n -1/2 U W,n (ψ 0 , ωn ) converges in probability to 0 when n → ∞.

To justify this, we note

n -1/2 U W,n (ψ 0 , ωn ) = n -1/2 U W,n (ψ 0 , ωn ) -n -1/2 U W,n (ψ 0 , ω 0 ) + n -1/2 U W,n (ψ 0 , ω 0 ). Then n -1/2 U W,n (ψ 0 , ωn ) -n -1/2 U W,n (ψ 0 , ω 0 ) = 1 n n i=1 δ i 1 a( ω n , D i ) - 1 a(ω 0 , D i ) Ψ i (ψ 0 ) , ≤ 1 n n i=1 a(ω 0 , D i ) -a(ω n , D i ) a(ω 0 , D i )a(ω n , D i ) Ψ i (ψ 0 ) , ≤ 1 n n i=1 a(ω 0 , D i ) -a(ω n , D i ) Ψ i (ψ 0 ) × a(ω 0 , D i )a(ω n , D i ) -1
.

Conditions (A.2) and (A.4) guarantee the existence of a positive constant k 1 such that

n -1/2 U W,n (ψ 0 , ωn ) -n -1/2 U W,n (ψ 0 , ω 0 ) ≤ k 1 n n i=1 π(ω 0 , D i ) -π(ω n , D i ) .
Moreover, under (A.3), it follows that

n -1/2 U W,n (ψ 0 , ωn ) -n -1/2 U w,n (ψ 0 , ω 0 ) ≤ k 1 n n i=1 κ(O i ) ωn -ω 0 ) , ≤ k 1 ρ + o P (1) ωn -ω 0 .
Hence, consistency of ωn to ω 0 implies that n -1/2 U W,n (ψ 0 , ωn ) -n -1/2 U W,n (ψ 0 , ω 0 ) converges in probability to 0. Then, n -1/2 U W,n (ψ 0 , ω 0 ) converges in probability to 0. Finally, n -1/2 U W,n (ψ 0 , ω 0 ) converges in probability to 0 as n → ∞. Therefore, it follows that n -1/2 U W,n (ψ 0 , ωn ) converges to 0 when n tends to infinity.

Thirdly, we show that n -1/2 ∂U W,n (ψ,ωn) ∂ψ ⊤ converges in probability to a fixed function Λ(ψ, ω 0 ), uniformly in an open neighbourhood of ψ 0 when n → ∞.

To see this, let's take ψ ∈ V ψ 0 , we have

n -1/2 ∂U W,n (ψ, ωn ) ∂ψ ⊤ = n -1/2 ∂U W,n (ψ, ωn ) ∂ψ ⊤ -n -1/2 ∂U W,n (ψ, ω 0 ) ∂ψ ⊤ + n -1/2 ∂U W,n (ψ, ω 0 ) ∂ψ ⊤ . Convergence to 0 of n -1/2 ∂U W,n (ψ,ωn) ∂ψ ⊤ -n -1/2 ∂U W,n (ψ,ω 0 ) ∂ψ ⊤
is a consequence of the consistency of ω n to ω 0 , under conditions (A.2) to (A.5). This can be shown using the same arguments as in the proof of the second condition. Other details of the calculations are omitted.

Next, by (A.5), we have n -1/2 ∂U W,n (ψ,ωn) ∂ψ ⊤ converges in probability to the matrix Λ(ψ, ω 0 ) when n → ∞. Hence, it follows that n -1/2 ∂U W,n (ψ,ωn) ∂ψ ⊤ converges in probability to the matrix Λ(ψ, ω 0 ). Next, under conditions (A.2) and (A.5), n → ∞, n -1/2 ∂U W,n (ψ,ωn) ∂ψ ⊤ converges in probability to Λ(ψ, ω 0 ), uniformly on V ψ 0 . Condition 3 is therefore checked. The three conditions of [START_REF] Foutz | On the unique consistent solution to the likelihood equations[END_REF] inverse function theorem are verified. Thus, we conclude that ψW,n converges in probability to ψ 0 . Now let's prove that ψW,n is asymptotically Gaussian. First, we do a Taylor series expansion of U W,n ( ψn , ωn ) at (ψ 0 , ω 0 ). We get,

0 = U W,n ( ψn , ωn ) = U W,n (ψ 0 , ω 0 ) + ∂U W,n (ψ 0 , ω 0 ) ∂ψ ⊤ ( ψn -ψ 0 ) + ∂U W,n (ψ 0 , ω 0 ) ∂ω ⊤ (ω n -ω 0 ) + o P (1). It follows that 0 = U W,n (ψ 0 , ω 0 )+n -1/2 ∂U W,n (ψ 0 , ω 0 ) ∂ψ ⊤ √ n( ψn -ψ 0 )+n -1/2 ∂U W,n (ψ 0 , ω 0 ) ∂ω ⊤ √ n(ω n -ω 0 )+o P (1).
Furthermore, it can be shown that

n -1/2 ∂U W,n (ψ 0 ,ω 0 ) ∂ω ⊤ converges in probability to Γ(ψ 0 , ω 0 ). It follows that √ n( ψn -ψ 0 ) = - 1 √ n ∂U W,n (ψ 0 , ω 0 ) ∂ψ ⊤ -1 U W,n (ψ 0 , ω 0 ) + 1 √ n ∂U W,n (ψ 0 , ω 0 ) ∂ω ⊤ √ n(ω n -ω 0 ) + o P (1). = - 1 √ n ∂U W,n (ψ 0 , ω 0 ) ∂ψ ⊤ -1 U W,n (ψ 0 , ω 0 ) + n -1/2 n i=1 δ i -a(ω, D i ) a(ω, D i ) 1 -a(ω, D i ) Γ(ψ 0 , ω 0 )Ω -1 (ω 0 ) ȧ(ω, D i ) + o P (1).
We have already shown that -1

√ n ∂U W,n (ψ 0 ,ω 0 ) ∂ψ ⊤
converges in probability to Λ(ψ 0 ) (in the proof of condition 3).Then, we have

E δ -a(ω 0 , D) a(ω 0 , D) 1 -a(ω 0 , D) 2 ȧ(ω 0 , D) ȧ(ω 0 , D) ⊤ = Ω(ω 0 ). Therefore, we have var δ -a(ω 0 , D) a(ω, D) 1 -a(ω 0 , D) Γ(ψ 0 , ω 0 )Ω -1 (ω 0 ) ȧ(ω 0 , D) = Γ(ψ 0 , ω 0 )Ω -1 (ω 0 ) Γ(ψ 0 , ω 0 ) ⊤ . Besides, cov δ i a(ψ 0 , ω 0 ) Ψ i (ψ 0 ), δ i -a(ω 0 , D i ) a(ω 0 , D i ) 1 -a(ω 0 , D i ) Γ(ψ 0 , ω 0 )Ω -1 (ω 0 ) ȧ(ω 0 , D i ) = E δ i a(ω 0 , D i ) 2 Ψ i (ψ 0 ) ȧ(ω 0 , D i ) ⊤ δ i -a(ω 0 , D i ) 1 -a(ω, D i ) Ω -1 (ω 0 ) Γ(ψ 0 , ω 0 ) ⊤ , = E δ i a(ω 0 , D i ) 2 Ψ i (ψ 0 ) ȧ(ω 0 , D i ) ⊤ Ω -1 (ω 0 ) Γ(ψ 0 , ω 0 ) ⊤ , = -Γ(ψ 0 , ω 0 )Ω -1 (ω 0 ) Γ(ψ 0 , ω 0 ) ⊤ . Since var U W,n (ψ 0 , ω 0 ) = Φ(ψ 0 , ω 0 ). Hence var U W,n (ψ 0 , ω 0 ) + n -1/2 n i=1 δ i -a(ω, D i ) a(ω, D i ) 1 -a(ω, D i ) Γ(ψ 0 , ω 0 )Ω -1 (ω 0 ) ȧ(ω, D i ) = Φ(ψ 0 , ω 0 ) -Γ(ψ 0 , ω 0 )Ω(ω 0 ) -1 Γ(ψ 0 , ω 0 ) ⊤ .
Therefore, using Central Limit Theorem and Slutsky's Theorem, we have √ n( ψn -ψ 0 ) converges in distribution to a Gaussian vector with mean zero and covariance matrix ∆ W .

Proof of Theorem 3.2. The consistency of ψ W s can be established by verifying the conditions of inverse function theorem of [START_REF] Foutz | On the unique consistent solution to the likelihood equations[END_REF]. We first show that ψ W s is a consistent estimator of ψ. To do so, we begin by proving that U W s,n (ψ, π) -U W,n (ψ, π) converges in probability to 0 when n tends to infinity. Noting that, |π -π| = η n , where η n = {nh 2r + 1/(nh 2d )} 1/2 . Using some standard techniques,

U W s,n (ψ, π) -U W,n (ψ, π) = 1 √ n n i=1 1 π(D b i , D c i ) - 1 π(D b i , D c i ) δ i Ψ i (ψ), = - 1 √ n n i=1 π(D b i , D c i ) -π(D b i , D c i ) π 2 (D b i , D c i ) + O P π(D b i , D c i ) -π(D b i , D c i ) 2 δ i Ψ i (ψ), = - 1 √ n n i=1 π(D b i , D c i ) -π(D b i , D c i ) π 2 (D b i , D c i ) + O P (η 2 n ) δ i Ψ i (ψ), = - 1 √ n n i=1 n j=1 δ j K h (D b j = D b i , D c j -D c i )/ n j=1 K h (D b j = D b i , D c j -D c i ) -π(D b i , D c i ) π 2 (D b i , D c i ) + O P (η 2 n ) δ i Ψ i (ψ), = - 1 n 3/2 n j=1 n i=1 [δ j -π(D b i , D c i )]K h (D d j = D d i , D c j -D c i ) π 2 (D b i , D c i )h d f (D d = D d i , D c -D c i ) δ i Ψ i (ψ) + O P (η n ), = - 1 n 3/2 n j=1 n i=1 δ j -π(D b i , D c i ) δ i -π(D b i , D c i ) K h (D b j = D b i , D c j -D c i ) π 2 (D b i , D c i )h d f (D b = D b i , D c -D c i ) Ψ i (ψ) - 1 n 3/2 n j=1 n i=1 π(D b i , D c i ) δ j -π(D b i , D c i ) K h (D b j = D b i , D c j -D c i ) π 2 (D b i , D c i )h d f (D b = D b i , D c -D c i ) Ψ i (ψ) + O P (η n ), = -T 1n -T 2n + O P (η n ),
where

T 1n = 1 √ n 3 n j=1 n i=1 δ j -π(D d j , D c j ) δ i -π(D b i , D c i ) K h (D b j = D b i , D c j -D c i ) π 2 (D b i , D c i )h d f (D b = D b i , D c -D c i ) Ψ j (ψ), T 2n = 1 √ n 3 n j=1 n i=1 π(D b i , D c i ) δ j -π i (D b i , D c i ) K h (D b j = D b i , D c j -D c i ) π 2 (D b i , D c i )h d f (D b = D b i , D c -D c i ) Ψ j (ψ).
Let Ψ * i (ψ) = E Ψ(ψ|D i ) . T 2n can be expressed as follows :

T 2n = 1 √ n 3 n j=1 n i=1 π(D b i , D c i ) δ j -π(D b i , D c i ) K h (D b i = D b j , D c i -D c j ) π 2 (D b i , D c i )h d f (D b = D b i , D c -D c i ) Ψ i (ψ) = 1 √ n n j=1 δ j -π(D b j , D c j ) π(D b j , D c j ) n i=1 K h (D b i = D b j , D c i -D c j )Ψ i (ψ) h d f (D b = D b i , D c -D c i ) = 1 √ n n j=1 [δ j -π(D b j , D c j )] π(D b j , D c j ) Ψ * i (ψ) + O P (η 2 n ) T 2n = 1 √ n n j=1 δ j -π(D b j , D c j ) π(D b j , D c j ) Ψ * j (ψ) + O P (η n ).
Based on T 1n , let

τ ij = δ j -π(D b i , D c i ) δ i -π(D b i , D c i ) K h (D b j = D b i , D c j -D c i ) π 2 (D b i , D c i )h d f (D b = D b i , D c -D c i ) Ψ i (ψ).
Assessing, E(T 1n ) (see (?),and [START_REF] Lee | Estimation of a zero-inflated Poisson regression model with missing covariates via nonparametric multiple imputation methods[END_REF]), it can be shown that var

(T 1n ) = O(η n ) and E(T 1n ) = O(η n ) and var(T 1n ) = O(η n ). Finally, U W s,n (ψ, π) -U W,n (ψ, π) becomes U W s,n (ψ, π) -U W,n (ψ, π) = 1 √ n n i=1 δ i -π(D b i , D c i ) π(D b i , D c i ) Ψ * i (ψ) + O P (η n ).
(6.9)

Hence, we have

U W s,n (ψ, π) = 1 √ n n i=1 δ i Ψ i (ψ) π(D b i , D c i ) + δ i -π(D b i , D c i ) π(D b i , D c i ) Ψ * i (ψ) + O P (η n ).
By mimicking the proof of (6.9), Σ 1,n (ψ, π)

-Σ 1,n (ψ, π) = D P (η n ), where Σ 1,n (ψ, π) = -1 √ n ∂U W,n (ψ,π) ∂ψ ⊤ . Because E Σ 1,n (ψ, π) = 1 n n i=1 E E δ i π(D i ) -∂Ψ i (ψ) ∂ψ ⊤ D i = E -∂Ψ i (ψ) ∂ψ ⊤
= Σ 1 (ψ), it can be shown with the weak law of large numbers that Σ 1,n (ψ, π) converges in probability to Σ 1 (ψ). So, using Slutsky's theorem, it can be shown that Σ 1,n (ψ, π) converges in probability to Σ 1 (ψ). It can be shown by assumption (B.6) that the convergence of Σ 1,n (ψ, π) to Σ 1 (ψ) is uniform in a neighborhood of the true ψ. Therefore, by the inverse function theorem of [START_REF] Foutz | On the unique consistent solution to the likelihood equations[END_REF], there exists a unique consistent solution of the estimating equation U W s,n (ψ, π) = 0 in a neighborhood of the true ψ. Thus, it follows that ψ W s is a consistent estimator of ψ. Now, we derive the asymptotic distribution of √ n( ψ W s -ψ). By a Taylor expansion of U W s,n ( ψ W s , π) around (ψ, π), we have

0 = U W s,n ( ψ W s , π) = U W s,n (ψ, π) + 1 √ n ∂U W s,n (ψ, π) ∂ψ ⊤ √ n( ψ W s -ψ) + O P (η n ). Further √ n( ψ W s -ψ) = 1 √ n ∂U W s,n (ψ, π) ∂ψ ⊤ -1 U W s,n (ψ, π) + O P (η n ) = 1 √ n ∂U W s,n (ψ, π) ∂ψ ⊤ -1 U W,n (ψ, π) + U W s,n (ψ, π) -U W,n (ψ, π) + O P (η n ).
We have

var U W s,n (ψ, π)] = var U W,n (ψ, π) + (U W s,n (ψ, π) -U W,n (ψ, π)) = var U W,n (ψ, π) + var U W s,n (ψ, π) -U W,n (ψ, π) + 2cov U W,n (ψ, π), (U W s,n (ψ, π) -U W,n (ψ, π)) .
On the other hand, we have

var U W,n (ψ, π) = E 1 π(D b 1 , D c 1 ) Ψ 1 (ψ) ⊗2 . Then, it follows cov U W,n (ψ, π), U W s,n (ψ, π) -U W,n (ψ, π) = -E δ 1 δ 1 -π(D b 1 , D c 1 ) π 2 (D b 1 , D c 1 ) Ψ 1 (ψ) ⊗2 + O(η n ), = -E E δ 1 δ 1 -π(D b 1 , D c 1 ) π 2 (D b 1 , D c 1 ) Ψ * 1 (ψ) ⊗2 D b 1 , D c 1 + O(η n ), = -E 1 -π(D b 1 , D c 1 ) π(D b 1 , D c 1 ) Ψ * 1 (ψ) ⊗2 + O(η n ),

B. 1

 1 Let supp(D) denote the support of D = (D b , D c ). For any o ∈ supp(D), there exists c 1 such that the selection probability π(o) > c 1 > 0. In addition, π(o) has r continuous and bounded partial derivatives with respect to the continuous components of D almost everywhere.

  that is the solution of the estimating equations in Equation(3.7) ψ W,s : semiparametric IPW estimator that is the solution of the estimating equations in Equation (3.8) Sample sizes of n = 500, and 2000 were considered. Results are based on 1000 replications.

  finite and positive definite in a neighborhood of the true ψ.

B.6 The first derivatives of U W s,n (ψ, π) with respect to ψ exist almost surely in a neighborhood of the true ψ. Moreover, in such a neighborhood, the first derivatives are uniformly bounded above by a function of D, whose expectations exist. Theorem 3.2. Assume that assumptions (A.1) and (B.1) -(B.6) hold. Then, ψ W s is a consistent estimator of ψ and √ n( ψ W s -ψ) is asymptotically normally distributed with mean zero and covariance matrix Σ W s , where

Table 1 :

 1 Simulation results for ψ = (-0.6, 0.75, -1, 0.5, -0.25, 0.25, -0.4, 0.7, 0) ⊤ . The case where average fraction of missing data (25%) and average zero-inflated data (40%).

					βn				γn
			β1,n	β2,n	β3,n	β4,n	β5,n	γ1,n	γ2,n	γ3,n	γ4,n
	n = 500
	ψ FD	bias SD	-0.0450 -0.0202 0.0225 -0.0249 -0.0201 0.2978 0.1953 0.3473 0.2106 0.3121	0.0115 0.6065	-0.0292 0.6278	0.1250 0.2748	-0.0697 0.6469
		SE	0.2956	0.1985	0.3644	0.2206	0.3711	0.6845	0.6431	0.3014	0.7323
		RMSE	0.2997	0.1954	0.3463	0.2110	0.3112	0.6036	0.6253	0.3007	0.6474
	ψ CC	bias SD	0.1356 0.3060	-0.1333 0.0494 -0.0200 -0.2227 0.2134 0.3602 0.2173 0.3213	-0.2380 -0.0136 0.7527 0.8928	0.3341 0.4478	-0.2304 0.8598
		SE	0.3273	0.2197	0.3905	0.2465	0.3920	0.9377	0.9350	0.5456	1.0066
		RMSE	0.3333	0.2507	0.3618	0.2171	0.3896	0.7858	0.8884	0.5569	0.8859
	ψ W,P	bias	-0.0435 -0.0133 0.0227 -0.0174 -0.0204	0.0150	-0.0250	0.1263	-0.0675
		SD	0.2994	0.2059	0.3695	0.2116	0.3196	0.6526	0.7303	0.3157	0.7136
		SE	0.2998	0.2036	0.3690	0.2269	0.3756	0.6918	0.6545	0.3140	0.7411
		RMSE	0.3011	0.2053	0.3684	0.2112	0.3186	0.6495	0.7270	0.3385	0.7132
	ψ W,S	bias	-0.0313 -0.1421 0.0343 -0.0019 -0.0090	0.1280	0.0148	0.0768	-0.0694
		SD	0.3060	0.2138	0.3738	0.2195	0.3168	0.5943	0.6922	0.2775	0.6563
		SE	0.3060	0.2052	0.3792	0.2348	0.3812	0.6447	0.6080	0.2830	0.6932
		RMSE	0.3061	0.2558	0.3735	0.2184	0.3153	0.6049	0.6888	0.2866	0.6567
	n = 2000
	ψ FD	bias SD	-0.0149 -0.0269 0.0277 -0.0295 0.1452 0.0891 0.1687 0.0979	0.0048 0.1855	-0.0167 -0.0080 0.3204 0.2848	0.0295 0.1208	-0.0036 0.3383
		SE	0.1384	0.0915	0.1672	0.0959	0.1752	0.3109	0.2815	0.1212	0.3302
		RMSE	0.1459	0.0931	0.1709	0.1022	0.1855	0.3206	0.2848	0.1243	0.3382
	ψ CC	bias SD	0.1528 0.1499	-0.1201 0.0450 -0.0245 -0.1694 0.0971 0.1749 0.1028 0.1903	-0.2359 0.3960	0.0009 0.3574	0.1633 0.1586	-0.1367 0.4236
		SE	0.1518	0.0998	0.1777	0.1067	0.1827	0.3772	0.3509	0.1606	0.4083
		RMSE	0.2140	0.1544	0.1805	0.1057	0.2547	0.4607	0.3572	0.2276	0.4450
	ψ W,P	bias	-0.0168 -0.0254 0.0272 -0.0274	0.0067	-0.0137 -0.0076	0.0331	-0.0052
		SD	0.1502	0.0953	0.1771	0.1004	0.1880	0.3423	0.3060	0.1288	0.3617
		SE	0.1385	0.0917	0.1677	0.0963	0.1755	0.3121	0.2830	0.1220	0.3318
		RMSE	0.1511	0.0986	0.1791	0.1040	0.1881	0.3424	0.3060	0.1329	0.3616
	ψ W,s	bias	0.0117	-0.1330 0.0331 -0.0246	0.0137	0.0808	0.0393	-0.0066	0.0083
		SD	0.1495	0.0957	0.1774	0.0998	0.1875	0.3152	0.2886	0.1190	0.3292
		SE	0.1419	0.0930	0.1706	0.0997	0.1775	0.2962	0.2652	0.1131	0.3142
		RMSE	0.1498	0.1639	0.1804	0.1028	0.1879	0.3252	0.2911	0.1191	0.3291

Table 2 :

 2 Simulation results for ψ = (-0.6, 0.75, -1, 0.5, 0.2, -0.5, 0.8, -1.5) ⊤ . The case where average fraction of missing data (14 %) and average zero-inflated data (20%).

					βn				γn
			β1,n	β2,n	β3,n	β4,n	β5,n	γ1,n	γ2,n	γ3,n	γ4,n
	n = 500
	ψ FD	bias SD	-0.0280 -0.0229 0.0460 -0.0337 -0.0027 0.1976 0.1331 0.2387 0.1432 0.2506	-0.0413 -0.0177 0.0533 -0.0098 0.4714 0.4058 0.1838 0.4841
		SE	0.1981	0.1317	0.2413	0.1392	0.2512	0.4632	0.4194	0.1832	0.4897
		RMSE	0.1995	0.1349	0.2430	0.1471	0.2505	0.4730	0.4060	0.1913	0.4839
	ψ CC	bias SD	0.1391 0.2056	-0.1165 0.0593 -0.0277 -0.1807 0.1430 0.2475 0.1523 0.2557	-0.2743 -0.0191 0.2050 -0.1476 0.5958 0.5541 0.2563 0.6253
		SE	0.2180	0.1445	0.2577	0.1556	0.2626	0.5901	0.5674	0.2730	0.6318
		RMSE	0.2482	0.1844	0.2544	0.1547	0.3130	0.6557	0.5541	0.3281	0.6422
	ψ W,p	bias	-0.0296 -0.0189 0.0434 -0.0316 -0.0019	-0.0405 -0.0164 0.0604 -0.0108
		SD	0.2060	0.1432	0.2514	0.1490	0.2552	0.4945	0.4500	0.2049	0.5178
		SE	0.1986	0.1328	0.2427	0.1400	0.2517	0.4665	0.4238	0.1872	0.4942
		RMSE	0.2080	0.1444	0.2550	0.1522	0.2551	0.4959	0.4501	0.2135	0.5177
	ψ W,s	bias	-0.0063 -0.1294 0.0471 -0.0240	0.0069	0.0613	0.0304	0.0151	0.0018
		SD	0.2048	0.1411	0.2506	0.1475	0.2529	0.4560	0.4175	0.1877	0.4702
		SE	0.2036	0.1348	0.2478	0.1454	0.2551	0.4410	0.3952	0.1716	0.4656
		RMSE	0.2048	0.1914	0.2549	0.1494	0.2529	0.4598	0.4184	0.1882	0.4699
	n = 2000
	ψ FD	bias SD	0.0001 0.1134	-0.0272 0.0445 -0.0360 -0.0064 0.0763 0.1316 0.0815 0.1483	-0.0271 -0.0282 0.0640 -0.0982 0.3422 0.4427 0.2040 0.4587
		SE	0.1156	0.0775	0.1349	0.0809	0.1484	0.3371	0.4287	0.2031	0.4670
		RMSE	0.1133	0.0809	0.1389	0.0890	0.1484	0.3431	0.4434	0.2137	0.4689
	ψ CC	bias SD	0.0639 0.1185	-0.0472 0.0325 -0.0259 -0.0631 0.0825 0.1380 0.0851 0.1537	-0.1989 -0.0278 0.1468 -0.1089 0.4089 0.5194 0.2464 0.5474
		SE	0.1234	0.0832	0.1424	0.0872	0.1548	0.4066	0.5014	0.2479	0.5593
		RMSE	0.1346	0.0950	0.1417	0.0889	0.1661	0.4546	0.5199	0.2867	0.5578
	ψ W,p	bias	-0.0007 -0.0252 0.0416 -0.0347 -0.0064	-0.0182 -0.0304 0.0705 -0.1133
		SD	0.1168	0.0812	0.1374	0.0834	0.1510	0.3701	0.4814	0.2292	0.4992
		SE	0.1159	0.0777	0.1353	0.0811	0.1486	0.3466	0.4450	0.2107	0.4770
		RMSE	0.1167	0.0850	0.1435	0.0903	0.1510	0.3704	0.4821	0.2397	0.5117
	ψ W,s	bias	-0.0043 -0.0600 0.0384 -0.0321	0.0018	0.0865	0.0008	0.0223 -0.1032
		SD	0.1177	0.0811	0.1386	0.0833	0.1512	0.3375	0.4473	0.2132	0.4611
		SE	0.1172	0.0787	0.1368	0.0824	0.1495	0.3231	0.4050	0.1936	0.4491
		RMSE	0.1177	0.1009	0.1437	0.0892	0.1511	0.3482	0.4471	0.2143	0.4723

Table 3 :

 3 Descriptive summary (RAND Health Insurance Experiment Data).

	Variable min 1st Quartile Median 3rd Quartile max Mean Std dev
	mdvis	0.00	0.00	1.00	4.00	77.00	2.86	4.50
	age	0.00	11.46 24.20	37.42	64.28 25.72	16.77
	income	0	5334	7976	10475	29238 8037 4058.37
	school 0.00	11.00 12.00	13.00	25.00 11.96	2.81
	ghindx 0.00	0.00 67.00	81.50	100.00 54.18	34.84

Appendix A. R code for fitting the ZIMBT model

The code below fits the ZIMBT model to a data set simulated as in Section 4. In this code, b, g represent β, γ respectively. Before running this code, the user needs to specify the design matrices X and W. Using the distribution equation 1.1. X ∼ MBT(α), we have:

The following code builds the ZIMBT loglikelihood :

As a result, it follows that

Therefore, using the Central Limit Theorem, we have

Hence, using Slutsky's theorem, we have √ n( ψ W s -ψ) converges in distribution to a centered multivariate normal distribution with variance-covariance matrix Σ W s , where

Since ψ W s P → ψ and π P → π, it can be shown with the weak law of large numbers and Slutsky's theorem that Σ * 0,n ( ψ W s , π)

, where P → denotes convergence in probability. Hence, Σ W s,n is a consistent estimator of Σ W s .

Declarations

Conflict of interest. The authors have no conflicts of interest to declare that are relevant to the content of this article.