# A hierarchical Bayesian approach for incorporating expert opinions into parametric survival models: a case study of female Ixodes ricinus ticks exposed to various temperature and relative humidity conditions

Phrutsamon Wongnak, Severine Bord, Sophie Donnet, Thierry Hoch, Frederic Beugnet, Karine Chalvet-Monfray

► To cite this version:

Phrutsamon Wongnak, Severine Bord, Sophie Donnet, Thierry Hoch, Frederic Beugnet, et al. A hierarchical Bayesian approach for incorporating expert opinions into parametric survival models: a case study of female Ixodes ricinus ticks exposed to various temperature and relative humidity conditions. Ecological Modelling, 2022, 464, pp.109821. 10.1016/j.ecolmodel.2021.109821. hal-03428117

HAL Id: hal-03428117

https://hal.inrae.fr/hal-03428117

Submitted on 15 Nov 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

# A hierarchical Bayesian approach for incorporating expert opinions into parametric survival models: A case study of female *Ixodes ricinus* ticks exposed to various temperature and relative humidity conditions

Phrutsamon Wongnak<sup>a, b, *</sup>, Séverine Bord<sup>c</sup>, Sophie Donnet<sup>c</sup>, Thierry Hoch<sup>d</sup>, Frederic Beugnet<sup>e</sup>, Karine Chalvet-Monfray<sup>a, b</sup>

## ARTICLE INFO

**Keywords:**

Expert elicitation

Hierarchical Bayesian model

*Ixodes ricinus* ticks

Parametric survival model

## ABSTRACT

The survival of ectothermic species is heavily dependent on environmental conditions, such as temperature and water balance. Understanding their survival responses to abiotic factors could help predict impacts of climate change on their population dynamics and human health. However, making a statistical inference and formulating a predictive model for mortality rates can be challenging when the observation numbers are limited. This study proposed an expert opinion elicitation framework that integrates expert opinions as prior distributions for the effects of continuous explanatory variables, through a Bayesian Parametric Survival Model (B-PSM). A historical survival dataset of female *Ixodes ricinus* ticks (Acari: Ixodidae) with small sample size was used. A total of 6 acarologists were recruited as experts for interactive online interview sessions to provide their opinions on average survival time under 4 different temperature and humidity scenarios. Most experts shared similar opinions on the effects of abiotic variables, and none of the experts was confident in the interaction effect. The variation of the opinions across multiple experts was handled by two approaches: 1) pooling and 2) averaging methods. The results showed that the pooling approach retains the variations of expert opinions, it may also disregard some irrelevant opinions to the observed data. While the averaging approach forms a numerical consensus across all the experts, but it may be less informative when the opinions distinctly diverge. The survival time of *I. ricinus* was found to be best described by the Weibull distribution, suggesting the mortality rate of ticks increases over time (aging effects). Also, the posterior predictions revealed that *I. ricinus* ticks were susceptible to desiccation conditions, with an interaction effect with the temperature. Therefore, our results suggested that relative humidity is an important factor in the survival of *I. ricinus* that should not be disregarded when evaluating the impacts of climate change on their population dynamics. Finally, this study provided a guideline for implementing the B-PSM framework to incorporate expert opinions and develop predictive survival models that can be applied in other ecological contexts.

## 1. Introduction

Survival is a fundamental ecological process that defines the demographics of a population. In the face of anthropogenic climate change, the survival of invertebrate ectothermic species, such as insects and acari, has been considerably impacted by long-term alteration in abiotic conditions. Climate change is deemed to have a greater impact

on these invertebrate ectotherms than other stressors, such as changes in land use (Halsch et al., 2021). In order to survive such abiotic stress, species are forced to relocate to their remaining climatic niches, adapt to accommodate altered climates, or otherwise face extinction (Bates et al., 2014; Berg et al., 2010; Román-Palacios and Wiens, 2020). The survival of ectotherms is heavily dependent on environmental conditions and their physiological needs, such as optimal body temperature and

<footnote><sup>a</sup> Université Clermont Auvergne, INRAE, VetAgro Sup, UMR EPIA, 63122 Saint-Genès-Chamanelle, France</footnote>

<footnote><sup>b</sup> Université de Lyon, INRAE, VetAgro Sup, UMR EPIA, 69280 Marcy l'Etoile, France</footnote>

<footnote><sup>c</sup> Université Paris-Saclay, AgroParisTech, INRAE, UMR MIA-Paris, 75005, Paris, France</footnote>

<footnote><sup>d</sup> INRAE, Oniris, BIOEPAR, 44300, Nantes, France</footnote>

<footnote><sup>e</sup> Global Technical Services, Boehringer-Ingelheim Animal Health, 69007, Lyon, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: phrutsamon.wongnak@vetagro-sup.fr (P. Wongnak).</footnote>
adequate water balance (Rozen-Rechels et al., 2019). At the same time, rising temperatures may also help accelerate their growth and reproduction rates, reducing generation time and favouring evolutionary adaptation to the changing climate (Schmalensee et al., 2021). As a result, climate change may alter the distribution and community composition of ectothermic species globally, potentially affecting human well-being (Pecl et al., 2017). Climate change-induced alteration in the population dynamics of pollinators, such as bumblebees, and arthropod vectors, such as mosquitoes and ticks, may have had a negative impact on food security (Giannini et al., 2017) and infectious disease distributions (Dumic and Severnini, 2018; Lee et al., 2018), respectively. Understanding survival responses of invertebrate ectotherms to abiotic factors, particularly the interaction between thermo-regulation and hydro-regulation (Rozen-Rechels et al., 2019), could help predict climate change impacts on their population dynamics and human health (Nadeau et al., 2017).

A hard tick species *Ixodes ricinus* (Acari: Ixodidae) plays an important role in transmitting Lyme borreliosis, the most prevalent vector-borne disease in Europe. They spend the majority of their lives off-host, exposing themselves to surrounding environmental conditions. Like other tick species and invertebrate ectotherms, temperature and relative humidity have been shown to affect their survival (Needham and Teel, 1991). However, experimental studies reporting the effects of both temperature and relative humidity on the survival time of *I. ricinus* are currently limited (Herrmann and Gern, 2010; Lees, 1946; Milne, 1950). The most detailed reported dataset could be dated back to 1950, in which the sample size per experimental group was small (5 ticks/condition) (Milne, 1950). A predictive survival model that includes both temperature and relative humidity is still needed to evaluate the impacts of climate change on the population of *I. ricinus*, and eventually the distribution of Lyme borreliosis.

A Parametric Survival Model (PSM) is a statistical technique for analysing time-to-death data where the survival time is assumed to follow a probability distribution time, e.g., exponential, Weibull, or log-logistic distributions (Kleinbaum and Klein, 2012). PSM can be applied to explore the effects of explanatory variables on the survival and hazard/mortality rate as well as biological hypotheses about mortality rate's behavior (constant or time-varying). In the population ecology field, PSM has rarely been incorporated to estimate the impacts of several environmental factors on the mortality rate (Ergon et al., 2018), despite being a well-established statistical tool. One constraint for PSM application in certain ecological studies is the quality of available survival data, e.g., limited observable ranges of environmental variables, uncertainty on unobserved survival states, or limited sample sizes.

When existing observed data are deficient, expert opinion has been recognised as an adjunct or supplementary information that could help develop a statistical inference in ecological studies (Krueger et al., 2012). An expert opinion is an information given by knowledgeable individuals with in-depth experience in the topic of interest, often referred to as experts (Fazey et al., 2006). The variability of the opinions is typically handled by including multiple experts in the analyses, while the uncertainty around each expert's opinion can be addressed as probabilistic statements through the expert elicitation process (Albert et al., 2012; Colson and Cooke, 2018). Bayesian inference is a statistical framework that integrates *a priori* knowledge/belief, called prior distributions, and observed data to estimate *posterior* distributions of unknown parameters. It treats all quantities as random variables and handles probability as a measure of uncertainty. Therefore, the Bayesian framework could inherently accept elicited expert opinions as prior distributions for estimating unknown parameters (Kuhnert et al., 2010).

Theoretical simulation studies suggested that incorporating informative prior distributions for Bayesian survival models could improve posterior distributions of unknown parameters (Omurlu et al., 2015, 2009). Although employing expert opinions as prior distributions in Bayesian survival analysis is not currently widely used, it has been applied in various domains. For example, an expert elicitation

framework for the scale and shape parameters of the Weibull distribution has been established and applied in the engineering domain to study material degradation and fatigue processes (Bousquet, 2010; Compare et al., 2017; Singpurwalla, 1988). In the clinical research domain, expert opinions on the proportion of patients who survive within a pre-determined period have been used to assist in the prediction of survival functions in clinical studies, particularly on cancer patients (Cope et al., 2019; Hiance et al., 2009). In addition, a guideline for expert elicitation on clinical studies, including survival analysis, has been recently published (Bojke et al., 2019). In the ecological domain, expert opinions were used to estimate age-specific survival rates (Johnson et al., 2017) and assess the uncertainty of a judgment on the cause-of-death-specific mortality rate of wildlife (Walsh et al., 2018). To date, there has never been a report on expert opinion elicitation framework for the effects of multiple explanatory variables for survival analysis, particularly in ecological studies.

Therefore, this study proposed a hierarchical Bayesian PSM (B-PSM) framework that incorporates expert opinion as supplementary information for estimating the effects of continuous explanatory variables on the mortality rate when the existing survival dataset is limited. We explored a historical survival dataset of *I. ricinus* ticks exposed to various combinations of controlled temperature and relative humidity with a limitation of small sample size. We conducted online interview sessions with multiple experts, eliciting the opinions of the experts, then combined them (by either pooling or averaging the opinions) as prior distributions of the environmental effects on survival time in the B-PSM. As a result, we established expert opinion-guided predictive models for survival probability and mortality rate of female *I. ricinus* based on the temperature and relative humidity.

## 2. Material and methods

### 2.1. Tick survival data

The published survival dataset of *I. ricinus* ticks reported by Milne (1950) was used. Briefly, survival time *T* (in days) of female adult *I. ricinus* ticks was observed under a variety of laboratory conditions, controlling for temperature *Q* (5, 11, 19, and 25 °C) and relative humidity *U* (0, 50, 70, 85, and 95%). The age of adult ticks averaged 12 weeks post-moulting, at the beginning of the experiment. Five ticks were assigned to each exposure condition, pairing the temperature and relative humidity (*Q*, *U*). In total, the survival time of 100 ticks was observed in 20 different controlled conditions (Table S1). Considering the expected difficulties in achieving and maintaining a perfect dehydrated condition (*U* = 0%) in the laboratory by the time of the study (~1950), we substituted the condition *U* = 0% with *U* = 10% in further analysis.

### 2.2. Weibull survival model

We first assessed the most suitable distribution for the tick survival data from the exponential, Weibull, log-logistic, log-normal, logistic, and normal distributions, using the Frequentist Parametric Survival Model (F-PSM) in the *survival* package (Therneau, 2020; Therneau and Grambsch, 2000). Here, temperature *Q* and relative humidity *U* were treated as categorical explanatory variables. The F-PSM models were compared using the Akaike information criterion (AIC). As a preliminary result, the Weibull distribution yielded the smallest AIC, suggesting that it is the most suitable distribution of the survival data. We subsequently validated the suitability of the Weibull distribution by exploring the Weibull property of linearity of the natural log of time ln(*t*) (Kleinbaum and Klein, 2012) where the relationship between ln(*t*) and the log negative log of the Kaplan-Meier survival estimates ln(−ln*S*(*t*)) should be linear.

The Weibull distribution is a generalized form of the exponential distribution with two parameters: the scale parameter $\lambda$ and the shape
parameter p, denoted as $T \sim \mathcal{W}(\lambda, p)$. The probability density function of the survival time $f(t)$ can be written as $f(t) = S(t) \cdot h(t)$, where the survival $S(t)$, and hazard $h(t)$ functions are defined as $S(t) = \exp(-\lambda t^p)$, and $h(t) = \lambda p t^{p-1}$, respectively (Kleinbaum and Klein, 2012). The scale parameter $\lambda$ reflects the baseline hazard rate at $t = 1$; a higher $\lambda$ value indicates a higher mortality risk. While the shape parameter $p$ determines whether the hazard rate is constant ($p = 1$), increasing ($p > 1$) or decreasing ($p < 1$) over time. In other words, the shape parameter $p$ reflects the effects of age on the hazard rate. In our context, the survival time $T$ represents time-to-death; therefore, the hazard function $h(t)$ is equivalent to the mortality rate $\mu(t)$. Here, we described the tick survival time $T$ exposed to different conditions as:

$$T_{ij} \sim \mathcal{W}(\lambda_j, p) \qquad (1)$$

$T_{ij}$ is the survival time of tick $i \in \{1, ..., 5\}$, in an experimental condition $j \in \{1, ..., 20\}$, corresponding to an abiotic variable pair ($Q_j, U_j$). We assumed that temperature and relative humidity influenced the baseline hazard rate while not affecting the aging effects. Therefore, we assumed that $Q_j$ and $U_j$ have an effect on the scale parameter $\lambda_j$, denoted as $\lambda_j$, while $p$ remains constant across all experimental conditions.

## 2.3. Relationship between abiotic variables and survival time

To establish the form of a predictive model, the relationship between the tick survival time $T$ and abiotic variables ($Q, U$) were explored using the F-PSM approach. Here, the abiotic variables ($Q, U$) were treated as continuous explanatory variables. According to the preliminary results from Section 2.2, the survival time $T$ was assumed to follow a Weibull distribution. The AIC value was used to compare the goodness-of-fit of different tested F-PSM models.

According to the characteristics of the Weibull distribution, the survival regression model can be formulated in two ways: the proportional hazards and the acceleration failure time approaches (Kleinbaum and Klein, 2012). Here, we constructed the model as a proportional hazards regression model that estimates temperature and relative humidity effects through a log link function for the subsequent hierarchical B-PSM analysis. The proportional hazards approach was chosen for its simplicity of interpretation. A value of the survival regression coefficients greater than 0 indicates a higher risk, while a value less than 0 indicates a protective effect. The model formulation was modified from the F-PSM model with the lowest AIC value as:

$$\ln \lambda_j = \beta_0 + \beta_1 U_j^k + \beta_2 Q_j + \beta_3 U_j^k Q_j \qquad (2)$$

Given a condition $j$, the log-transformed of the scale parameter $\ln\lambda_j$ is described by a combination of non-linear effects of relative humidity $U_j$, and linear effects of temperature $Q_j$. Let $B = [\beta_0, \beta_1, \beta_2, \beta_3]$ be a vector of the survival regression coefficients: $\beta_0$ indicates the value of $\ln\lambda_j$ at a reference condition ($Q_j = 0$ °C, $U_j = 0\%$); $\beta_1, \beta_2$, and $\beta_3$ indicate the effects of $U_j, Q_j$, and the interactions between $U_j$ and $Q_j$ on $\ln\lambda_j$, respectively. Also, $k$ is a parameter describing the degree of non-linear effects of $U_j$ as a continuous variable. Accordingly, we could describe the survival time $T_{ij}$ of tick $i$ exposed to the constant temperature $Q_j$ and relative humidity $U_j$ as:

$$T_{ij} | Q_j, U_j \sim \mathcal{W} \left( \exp \left( \beta_0 + \beta_1 U_j^k + \beta_2 Q_j + \beta_3 U_j^k Q_j \right), p \right) \qquad (3)$$

## 2.4. Expert opinion

A total of 6 acarologists (N) experienced in handling/breeding *I. ricinus* ticks under laboratory conditions were recruited. The objective was to gain a priori expert knowledge on the effects of abiotic variables ($Q, U$) on tick survival time, specifically on the parameter $B$ ($\beta_0, \beta_1, \beta_2, \beta_3$). The experts were requested for their opinions on the average survival time $\bar{T}$ of ticks in four different conditions, and they were subsequently transformed into $B$ through the elicitation process. Each expert

was interviewed separately in a 1-hour online session with a Shiny-based interactive Web application, developed using the *shiny* package (Chang et al., 2020).

### 2.4.1. Expert opinions on the average survival time

Upon starting the interview session, the interviewers delivered the background and objectives of this study. Subsequently, each expert $e \in \{1, ..., N\}$ was requested to provide their opinions on the average survival time $\bar{T}_{e,c}$ of 12-week-old unfed female adult *I. ricinus* ticks exposed to 4 controlled constant conditions $c \in \{1, ..., 4\}$, described by a couple “temperature $Q_{e,c}$; relative humidity $U_{e,c}$”. As an example, we initially proposed default values for “$Q_{e,c}$; $U_{e,c}$” corresponding to each condition $c$ and their brief descriptions as: 1) “$Q_{e,1}$; $U_{e,1}$” = “5 °C; 10%” (cold; dry); 2) “$Q_{e,2}$; $U_{e,2}$” = “25 °C; 10%” (warm; dry); 3) “$Q_{e,3}$; $U_{e,3}$” = “5 °C; 95%” (cold; humid); 4) “$Q_{e,4}$; $U_{e,4}$” = “25 °C; 95%” (warm; humid). We allowed the experts to adjust the default values and give their opinions on conditions most compatible with their prior experience, within the range of observed temperature (5 – 25 °C) and relative humidity (10 – 95%).

In each condition $c$, the experts were asked for the following parameters: 1) The mean of $\bar{T}_{e,c}$, denoted as $\bar{T}_{e,c}^m$; 2) The high, and 3) The low values of $\bar{T}_{e,c}$, denoted as $\bar{T}_{e,c}^h$ and $\bar{T}_{e,c}^d$, respectively; 4) A confidence level $C_{e,c}$, ranging from 0 to 1, corresponding to the degree of confidence on their opinions.

### 2.4.2. Optimization of the expert opinions into distribution

During the interview session, the parameters $\bar{T}_{e,c}^m, \bar{T}_{e,c}^h, \bar{T}_{e,c}^d$, and $C_{e,c}$ given by the experts were simultaneously optimized to find a corresponding log-normal distribution describing the uncertainty of the average survival time; $\bar{T}_{e,c} \sim \mathcal{LN}(\mu_{e,c}, \sigma_{e,c})$. The distributions were optimized with the following properties: 1) An average of the log-normal distribution equals $\bar{T}_{e,c}^m$; 2) The probability of having $\bar{T}_{e,c}$ between $\bar{T}_{e,c}^m$ and $\bar{T}_{e,c}^h$ equals $C_{e,c}$. Then, the Shiny-based Web application instantaneously visualized the probability density curves of $\bar{T}_{e,c}$ corresponding to expert's answers (Fig. 1). The experts could adjust the answers until the density curves agree with their opinions.

Upon finishing the interview, the Web application displayed the density curves of all four conditions side-by-side. The interviewers asked the experts to revise their answers (if necessary) before final submission. Then, the expert data $Q_{e,c}, U_{e,c}, \bar{T}_{e,c}^m, \bar{T}_{e,c}^h, \bar{T}_{e,c}^d$ and $C_{e,c}$, and their optimized parameters $\mu_{e,c}$, and $\sigma_{e,c}$ were recorded.

## 2.5. Hierarchical Bayesian model without expert opinions

Initially, a hierarchical B-PSM was employed to estimate the parameters $B, p$, and $k$, shown in (1) – (3), explicitly from the observed data without the information from the experts, referred as Model 1. All the parameters were provided with uniform prior distributions (Fig. 2A): each element of $B \sim \mathcal{N}[-50, 50]$; $p \sim \mathcal{N}[0, 5]$; $k \sim \mathcal{N}[1, 6]$. Without expert opinions, all values within the given ranges were assigned an equal probability to be included in the model. Ranges of uniform distributions were guided by the preliminary estimates of the F-PSM. The value of $k$ should lie between 3 and 4, therefore the lower and upper bounds of the uniform distribution were extended to 1 and 6, respectively. Besides, the intrinsic characteristic of the shape parameter $p$ defined the lower boundaries of the prior as it cannot be a negative value. Furthermore, statistical hypotheses can be evaluated by including critical values into the prior distributions: 1) A value of $B = 0$ indicates that the corresponding covariate does not affect the survival; 2) a value of $p = 1$ implies that the mortality rate is constant, whereas values of $p$ greater and less than 1 indicate that the mortality rate is increasing and decreasing over time, respectively.
Expert opinions:

A questionnaire on an average survival time of female adult *Ixodes ricinus* exposed
to constant temperature and relative humidity

Name and Last name :

Anonymous

Condition 1: Low humidity (dry) and low temperature (cold)

<figure><img src="image_4.png" /></figure>

**Fig. 1.** An interface of Shiny-based Web application used for the expert opinion interview. The experts were asked to provide their opinion on the average survival time of female *I. ricinus* ticks exposed to 4 constant laboratory conditions. The user-input panel (left) allows the experts to give the average survival time corresponding to conditions (predetermined or user-adjusted). Expert's inputs were simultaneously optimized for a log-normal distribution, with a density curve (upper right) and its parameters (lower right) displayed. The shaded area represents the probability under the curve between the lower and the higher value of *c*, which is equal to the confidence level *c*. The red vertical dashed line indicates the expected value of the optimized probability distribution, which is equal to the expert's average survival time.

<figure><img src="image_5.png" /></figure>

**Fig. 2.** Directed acyclic graphs (DAGs) of the hierarchical Bayesian models. The Bayesian models estimate the parameters *B* (β₀, β₁, β₂, and β₃), *p*, and *k*; (A) Model 1: without expert opinions; (B) Model 2: pooling expert opinions; (C) Model 3: averaging expert opinions. Arrows indicate the relationship between parameters (eclipse), observed data (double-bordered eclipse), covariates (small rectangle), and prior distributions (rounded rectangle): stochastic relationship (solid arrow); deterministic relationship (dashed arrow). *T*ᵢⱼ denoted the survival time of tick *i* exposed to experimental condition *j* (*Q*ᵢ, *U*ᵢ); λⱼ and *p* denoted the Weibull distribution's scale and shape parameters, respectively; *k* denoted a parameter describing the degree of non-linear effects of *U*ᵢ. Expert opinion *e* on condition *c* was optimized for hyperparameters μ<sub>e,c</sub><sup>l</sup> and σ<sub>e,c</sub><sup>l</sup>, describing the average survival time T̄<sub>e,c</sub> at temperature *Q*<sub>e,c</sub> and relative humidity *U*<sub>e,c</sub>. The relationship between the elicited expert data *B*<sub>e</sub> and *B* could be expressed as: 1) *B* = *B*<sub>e=e</sub>, where *ε* represents an expert, whose opinions were accounted in the model (Model 2); 2) *B* = 1/N Σ<sub>e=1</sub><sup>N</sup> *B*<sub>e</sub>, where *N* is the number of experts (Model 3).
## 2.6. Hierarchical Bayesian model including expert opinions

In this section, the parameters **B**, **p**, and **k**, shown in (1) – (3) were estimated using expert opinions incorporated in the hierarchical B-PSM framework.

### 2.6.1. Expert opinion elicitation

The expert opinion elicitation process transformed concrete biological quantities provided by experts, such as the average tick survival time $\bar{T}_{e,c}$, into theoretical quantities, such as the model parameters **B**. Here, the average tick survival time $\bar{T}_{e,c}$ provided by the experts were linked to the parameters of a Weibull distribution describing the tick survival time $T_{e,c} \sim \mathcal{W}(\lambda_{e,c}, p)$ through the following relationship:

$$E[T_{e,c}] = \bar{T}_{e,c} = \lambda_{e,c}^{-1/p} \cdot \Gamma(1 + 1/p) \quad (4)$$

Which is equivalent to:

$$\ln \lambda_{e,c} = -p \cdot \ln \frac{\bar{T}_{e,c}}{\Gamma(1 + 1/p)} \quad (5)$$

Also, $\ln\lambda_{e,c}$ is linked to the regression parameters through Eq. (2) as follows:

$$\ln \lambda_{e,c} = \beta_{0,e} + \beta_{1,e} U_{e,c}^k + \beta_{2,e} Q_{e,c} + \beta_{3,e} U_{e,c}^k Q_{e,c} \quad (6)$$

The relationship in Eq. (6) across all conditions $c \in \{1, ..., 4\}$ given by expert **e** can be expressed in a form of matrix multiplication as:

$$\begin{bmatrix} \ln \lambda_{e,1} \\ \ln \lambda_{e,2} \\ \ln \lambda_{e,3} \\ \ln \lambda_{e,4} \end{bmatrix} = \begin{bmatrix} 1 & U_{e,1}^k & Q_{e,1} & U_{e,1}^k Q_{e,1} \\ 1 & U_{e,2}^k & Q_{e,2} & U_{e,2}^k Q_{e,2} \\ 1 & U_{e,3}^k & Q_{e,3} & U_{e,3}^k Q_{e,3} \\ 1 & U_{e,4}^k & Q_{e,4} & U_{e,4}^k Q_{e,4} \end{bmatrix} \times \begin{bmatrix} \beta_{0,e} \\ \beta_{1,e} \\ \beta_{2,e} \\ \beta_{3,e} \end{bmatrix} \quad (7)$$

Therefore, we can express $\beta_{0,e}$, $\beta_{1,e}$, $\beta_{2,e}$, and $\beta_{3,e}$ by as a function of $\ln\lambda_{e,c}$ through Eq. (8).

$$\begin{bmatrix} \beta_{0,e} \\ \beta_{1,e} \\ \beta_{2,e} \\ \beta_{3,e} \end{bmatrix} = \begin{bmatrix} 1 & U_{e,1}^k & Q_{e,1} & U_{e,1}^k Q_{e,1} \\ 1 & U_{e,2}^k & Q_{e,2} & U_{e,2}^k Q_{e,2} \\ 1 & U_{e,3}^k & Q_{e,3} & U_{e,3}^k Q_{e,3} \\ 1 & U_{e,4}^k & Q_{e,4} & U_{e,4}^k Q_{e,4} \end{bmatrix}^{-1} \times \begin{bmatrix} \ln \lambda_{e,1} \\ \ln \lambda_{e,2} \\ \ln \lambda_{e,3} \\ \ln \lambda_{e,4} \end{bmatrix} \quad (8)$$

$$\text{Finally, let } \mathbf{B}_e = \begin{bmatrix} \beta_{0,e} \\ \beta_{1,e} \\ \beta_{2,e} \\ \beta_{3,e} \end{bmatrix}, \quad \mathbf{X}_e = \begin{bmatrix} 1 & U_{e,1}^k & Q_{e,1} & U_{e,1}^k Q_{e,1} \\ 1 & U_{e,2}^k & Q_{e,2} & U_{e,2}^k Q_{e,2} \\ 1 & U_{e,3}^k & Q_{e,3} & U_{e,3}^k Q_{e,3} \\ 1 & U_{e,4}^k & Q_{e,4} & U_{e,4}^k Q_{e,4} \end{bmatrix} \quad \text{and } \Lambda_e = \begin{bmatrix} \ln \lambda_{e,1} \\ \ln \lambda_{e,2} \\ \ln \lambda_{e,3} \\ \ln \lambda_{e,4} \end{bmatrix}, \quad \text{Eqs. (5) and (9) provided a simplified expression of Eq. (8)}$$

with uncertain components ($\bar{T}_{e,c}$, $p$, and $k$) indicated as:

$$\mathbf{B}_e = [\mathbf{X}_e | k]^{-1} \times \left[ \mathbf{\Lambda}_e \middle| \bar{T}_{e,c}, p \right] \quad (9)$$

Consequently, the elicited expert opinion $\mathbf{B}_e$ provides the a priori information for the model parameter **B**, where $\bar{T}_{e,c} \sim \mathcal{W}(\mu_{e,c}, \sigma_{e,c})$; $p \sim \mathcal{W}[0, 5]$; $k \sim \mathcal{W}[1, 6]$ (Figs. 2B and 2C).

### 2.6.2. Pooling expert opinion

The Bayesian Model 2 accepted the variation of expert opinions by pooling different experts as prior distributions. Elicited opinions $\mathbf{B}_e$ from expert $\varepsilon \in \{1, ..., N\}$ were chosen as prior distribution with equal initial probability $\pi = \{\pi_1, ..., \pi_N\} = \{\frac{1}{N}, ..., \frac{1}{N}\}$, where $N$ is the total

number of experts and experts were chosen using a categorical distribution as $\varepsilon \sim \text{Cat}(\pi)$. Therefore, the relationship between the unknown parameter **B** and $\mathbf{B}_e$ could be expressed as $\mathbf{B} = \mathbf{B}_{e=\varepsilon}$ (Fig. 2B). The final results of this model may take more into account the elicited opinions of some experts, while others may be disregarded.

### 2.6.3. Averaging expert opinion

The Bayesian Model 3 consolidated expert opinions by averaging the elicited expert data $\mathbf{B}_e$ across all $N$ experts. The relationship between the unknown parameter **B** and $\mathbf{B}_e$ could be expressed as: $\mathbf{B} = \frac{1}{N} \sum_{e=1}^{N} \mathbf{B}_e$ (Fig. 2C).

## 2.7. Implementation of Markov chain Monte Carlo algorithm

All the analyses in our study were performed using R programming language version 3.6.0 (R Core Team, 2019). The hierarchical Bayesian models were run by the Markov chain Monte Carlo (MCMC) algorithms in JAGS (Plummer, 2003) using rjags package (Plummer, 2019). For each Bayesian model, we run three independent MCMC chains consisting of 10,500,000 iterations in total, with a burn-in period of 50,000. The autocorrelation was controlled by thinning out the MCMC samples, keeping every 700 values. The convergence of MCMC chains was inspected by trace plots and the Gelman-Rubin convergence test. Probability of direction (p-direction) was used as an index of effect existence for each parameter using bayestestR package (Makowski et al., 2019).

## 2.8. Model evaluation and validation

The posterior distributions were evaluated by simulating survival time of 100 ticks in 20 conditions, then compared with the observed data. Subsequently, the validity of the Bayesian estimations was evaluated using posterior predictions. For each Bayesian model, a total of 15,000 values of each posterior estimates ($\mathbf{B}$, $\mathbf{p}$, and $\mathbf{k}$) were used to simulate the replication of other previously published data: 1) Survival time of 30 female *I. ricinus* ticks in 6 conditions (5 ticks/condition), compared to survival time ranges by Lees (1946); 2) Survival proportions of female *I. ricinus* ticks after 2- and 3-days post-exposure to 5 conditions (100 ticks/condition), compared against the survival proportion after 2 days reported by Herrmann and Gern (2010). The 95% confidence intervals for the Herrmann and Gern data were calculated by assuming a binomial distribution.

## 2.9. Posterior predictions

Median survival time and the log of scale parameter $\ln\lambda$ were predicted across the observed range of temperature (5 to 25 °C) and relative humidity (10%–95%), using the median of posterior distributions from all 3 models. Additionally, behaviours of the survival probability $S(t)$ and the mortality rate $\mu(t)$ of female *I. ricinus* ticks exposed to unobserved conditions (relative humidity of 55%, 65%, 75%, 85% and 95% at temperature of 5, 15 and 25 °C) were predicted following Eqs. (10) and (11), respectively.

$$S(t|\mathbf{Q}, \mathbf{U}) = \exp(-\exp(\beta_0 + \beta_1 \mathbf{U}^k + \beta_2 \mathbf{Q} + \beta_3 \mathbf{U}^k \mathbf{Q}) \cdot t^{\mu}) \quad (10)$$

$$\mu(t|\mathbf{Q}, \mathbf{U}) = \exp(\beta_0 + \beta_1 \mathbf{U}^k + \beta_2 \mathbf{Q} + \beta_3 \mathbf{U}^k \mathbf{Q}) \cdot t^{\mu-1} \quad (11)$$

# 3. Results

## 3.1. Descriptive analysis

The survival data of female *I. ricinus* ticks reported by Milne (1950) revealed a clear positive relationship between relative humidity and survival time (Fig. 3), while the negative effect of temperature on the survival time was less pronounced (Figure S1). The Kaplan-Meier analysis (Figure S2) showed that the tick survival time was significantly
<figure><img src="image_6.png" /><figcaption>Fig. 3. Relationship between relative humidity and survival time of female adult *I. ricinus* ticks exposed 4 temperature conditions (Q = 5 °C, 11 °C, 19 °C, and 25 °C). The vertical axes displayed the survival time *T* on a logarithmic scale. The data were originally reported by Milne (1950).</figcaption></figure>

longer in the conditions with higher relative humidity (log-rank test: *p*-value < 0.001, degree of freedom (df) = 4), while the effects of temperature were not statistically significant (log-rank test: *p*-value = 0.6, df = 3).

## 3.2. Exploring the distributions for survival time

The AIC values of survival regression models with categorical abiotic variables showed that the survival time *T* was best described by the Weibull distribution (Model A2; AIC = 822.91), followed by the logistic (Model A3; AIC = 843.30), and the log-normal (Model A4; AIC = 852.24) distributions (Table S2). The suitability for assuming the Weibull distribution to describe the survival time was supported by the property of linearity of ln(*t*). The relationship between ln(*t*) and ln(–

ln̂S(*t*)), treating relative humidity as categorical variables, was approximately linear (Figure S3).

## 3.3. Exploring the effects of abiotic variables

The AIC values of survival regression models assuming a Weibull distribution suggested that the survival time *T* could be explained by a combination of non-linear effect of relative humidity (raised to the power of 4), liner effect of temperature, and their interaction (Model B11; AIC 815.54; Table S3).

## 3.4. Expert opinions

The expert opinions on the average survival time $\bar{T}_{e,c}$ of female *I. ricinus* ticks exposed various experimental conditions (Table 1) were elicited for the distributions of *B* ($\beta_0$, $\beta_1$, $\beta_2$, and $\beta_3$) as shown in Fig. 4. Most experts shared similar opinions on the abiotic effects on tick survival. Experts 1, 2, 4, and 5 suggested a protective effect of relative humidity with $\beta_1$ significantly less than 0 (*p* < 0.05). While none of the experts was confident in the effect of temperature and the interactions with 0 included in 95% confidence intervals. In contrast, the average opinion across all 6 experts suggested a significant protective effect of relative humidity and a significant negative effect of temperature on tick survival (Table S4).

## 3.5. Hierarchical Bayesian parametric survival models

Fig. 5 shows the posterior distributions of *B*, *k*, and *p* and their summary statistics estimated from the Bayesian models. All the parameters from all 3 models converged to similar values. Providing Model 2 with mixed experts opinion as prior distributions reduced the size of 95% credible intervals for $\beta_0$ and $\beta_1$, while Model 3 reduced the size of 95% credible intervals for $\beta_0$, $\beta_1$, and $\beta_3$. The credible intervals of parameters *p* and *k* of Models 2 and 3 were slightly bigger than Model 1. All models estimated the negative value of $\beta_1$ (*p*-direction = 1.00), positive value of $\beta_2$ (*p*-direction = 1.00), and the shape parameter of the Weibull distribution *p* greater than 1 (*p*-direction = 1.00). The interaction effect

<caption>**Table 1**</caption>

The opinions of experts *e* on the average survival time $\bar{T}_{e,c}$ of female *I. ricinus* ticks exposed to experimental conditions *c* of temperature $Q_{e,c}$ and relative humidity $U_{e,c}$. $\bar{T}_{e,c}^m$ denoted the mean of average survival time; $\bar{T}_{e,c}^h$ and $\bar{T}_{e,c}^d$ denoted the high and low range of $\bar{T}_{e,c}$; $C_{e,c}$ denoted the confidence level of the experts on their opinions. The expert data were optimized to log-normal distribution describing the uncertainty of $\bar{T}_{e,c} \sim LN(\mu_{e,c}, \sigma_{e,c})$.

<table><thead><tr><th>Expert<br>e</th><th>Condition<br>c</th><th>Temperature (°C)<br>Q<sub>e,c</sub></th><th>Relative humidity<br>U<sub>e,c</sub></th><th colspan="3">Average survival time (days)</th><th>Confidence level<br>C<sub>e,c</sub></th><th colspan="2">Optimized parameter<br>μ<sub>e,c</sub> σ<sub>e,c</sub></th></tr><tr><th></th><th></th><th></th><th></th><th>T&#773;<sup>m</sup><sub>e,c</sub></th><th>T&#773;<sup>h</sup><sub>e,c</sub></th><th>T&#773;<sup>d</sup><sub>e,c</sub></th><th></th><th></th><th></th></tr></thead><tbody><tr><td rowspan="4">1</td><td>1</td><td>5</td><td>0.30</td><td>42</td><td>14</td><td>50</td><td>0.60</td><td>3.73</td><td>0.73</td></tr><tr><td>2</td><td>25</td><td>0.30</td><td>10</td><td>5</td><td>20</td><td>0.80</td><td>2.29</td><td>0.52</td></tr><tr><td>3</td><td>5</td><td>0.95</td><td>90</td><td>60</td><td>120</td><td>0.90</td><td>4.50</td><td>0.21</td></tr><tr><td>4</td><td>25</td><td>0.95</td><td>150</td><td>120</td><td>180</td><td>0.95</td><td>5.01</td><td>0.10</td></tr><tr><td rowspan="4">2</td><td>1</td><td>5</td><td>0.10</td><td>10</td><td>5</td><td>15</td><td>0.90</td><td>2.30</td><td>0.32</td></tr><tr><td>2</td><td>25</td><td>0.10</td><td>5</td><td>3</td><td>8</td><td>0.90</td><td>1.60</td><td>0.30</td></tr><tr><td>3</td><td>7</td><td>0.90</td><td>250</td><td>120</td><td>360</td><td>0.90</td><td>5.52</td><td>0.30</td></tr><tr><td>4</td><td>20</td><td>0.90</td><td>60</td><td>30</td><td>120</td><td>0.80</td><td>4.09</td><td>0.52</td></tr><tr><td rowspan="4">3</td><td>1</td><td>5</td><td>0.10</td><td>15</td><td>7</td><td>30</td><td>0.60</td><td>2.69</td><td>0.81</td></tr><tr><td>2</td><td>25</td><td>0.10</td><td>7</td><td>4</td><td>15</td><td>0.90</td><td>1.94</td><td>0.37</td></tr><tr><td>3</td><td>5</td><td>0.95</td><td>305</td><td>244</td><td>365</td><td>0.80</td><td>5.72</td><td>0.16</td></tr><tr><td>4</td><td>20</td><td>0.95</td><td>365</td><td>275</td><td>397</td><td>0.95</td><td>5.90</td><td>0.05</td></tr><tr><td rowspan="4">4</td><td>1</td><td>5</td><td>0.30</td><td>180</td><td>30</td><td>365</td><td>0.50</td><td>5.19</td><td>1.00</td></tr><tr><td>2</td><td>20</td><td>0.30</td><td>30</td><td>7</td><td>60</td><td>0.70</td><td>3.38</td><td>1.00</td></tr><tr><td>3</td><td>5</td><td>0.95</td><td>1277.5</td><td>700</td><td>2190</td><td>0.95</td><td>7.15</td><td>0.29</td></tr><tr><td>4</td><td>20</td><td>0.80</td><td>547.5</td><td>270</td><td>912.5</td><td>0.80</td><td>6.31</td><td>0.47</td></tr><tr><td rowspan="4">5</td><td>1</td><td>5</td><td>0.10</td><td>2</td><td>1</td><td>3</td><td>0.95</td><td>0.68</td><td>0.26</td></tr><tr><td>2</td><td>25</td><td>0.10</td><td>2</td><td>1</td><td>3</td><td>0.95</td><td>0.68</td><td>0.26</td></tr><tr><td>3</td><td>5</td><td>0.95</td><td>365</td><td>270</td><td>545</td><td>0.95</td><td>5.90</td><td>0.17</td></tr><tr><td>4</td><td>25</td><td>0.95</td><td>270</td><td>180</td><td>365</td><td>0.95</td><td>5.60</td><td>0.18</td></tr><tr><td rowspan="4">6</td><td>1</td><td>5</td><td>0.10</td><td>8</td><td>5</td><td>15</td><td>0.50</td><td>2.05</td><td>0.73</td></tr><tr><td>2</td><td>25</td><td>0.10</td><td>2</td><td>1</td><td>5</td><td>0.70</td><td>0.57</td><td>0.69</td></tr><tr><td>3</td><td>8</td><td>0.95</td><td>120</td><td>90</td><td>180</td><td>0.80</td><td>4.79</td><td>0.25</td></tr><tr><td>4</td><td>25</td><td>0.95</td><td>18</td><td>15</td><td>30</td><td>0.50</td><td>2.89</td><td>0.43</td></tr></tbody></table>
<figure><img src="image_7.png" /><figcaption>**Fig. 4.** Probability distributions of elicited expert opinions on each element of B. The distributions from Experts 1 to 6 were used as prior distributions for Bayesian Model 2, while the averaged distributions were used for Bayesian Model 3. Given a condition *j* with temperature *Q*<sub>*j*</sub> and relative humidity *U*<sub>*j*</sub>, the logarithm of scale parameter of the Weibull distribution lnλ<sub>*j*</sub> is expressed as: lnλ<sub>*j*</sub> = β<sub>0</sub> + β<sub>1</sub>U<sub>*j*</sub><sup>k</sup> + β<sub>2</sub>Q<sub>*j*</sub> + β<sub>3</sub>U<sub>*j*</sub><sup>k</sup>Q<sub>*j*</sub>; A) β<sub>0</sub> indicates the value of log scale parameter lnλ<sub>*j*</sub> at a reference condition; B) β<sub>1</sub>, C) β<sub>2</sub>, and D) β<sub>3</sub> indicate the effects of U<sub>*j*</sub>, Q<sub>*j*</sub>, and the interactions between U<sub>*j*</sub> and Q<sub>*j*</sub> on lnλ<sub>*j*</sub>, respectively. Red vertical dashed lines mark the value of 0.</figcaption></figure>

between relative humidity and temperature (β₃) was estimated to be
significant by Model 1 (p-direction = 0.986), Model 2 (p-direction =
0.973), and Model 3 (p-direction = 0.998). Additionally, most of the
posterior estimates from Model 2 were taken from the prior distributions
of Expert 5 (51.86%), followed by Expert 4 (48.11%), and Expert 2
(0.03%), respectively.

3.6. *Model evaluation and validation*

The 95% credible interval of posterior predictions by all 3 models
captured the observed survival time of female *I. ricinus* ticks in all 20
conditions (Figure S4). The posterior predictions were similar across all
the Bayesian models with some remarkable exceptions: the predicted
survival times from Models 2 and 3 were slightly longer at high humidity
and high temperature than those of Models 1. Besides, the posterior
estimates were validated against the survival time ranges in all 6 con-
ditions reported by Lees (1946) (Figure S5), and they adeptly repro-
duced the proportion of tick survival after 2-days exposure to 4 of 5
conditions reported by Herrmann and Gern (2010), except for one
condition at Q = 25 °C and U = 89% (Figure S6).

3.7. *Predicted survival probability and mortality rate*

Predicted median survival time and log scale parameter of the Weibull distribution lnλ from all 3 Bayesian models were similar, except for conditions with the relative humidity close to 95% (Fig. 6). At high humidity, the predictions from Models 2 and 3 were less sensitive to the negative effects of high temperature. In general, the predicted lnλ and μ(t) were considerably sensitive to the temperature in dry conditions (Fig. 6 and S7). However, while the excessively high μ(t) induces rapid death, it does not contribute to notable differences in predicted S(t) at low humidity (Figure S7). In contrast, S(t) is susceptible to minor changes in μ(t) at high humidity. Additionally, the predicted tick survivals at a relative humidity of 95% were noticeably longer than at 85%.

4. Discussion

Disentangling the impacts of temperature and water regulations on
the physiological performance and survival of ectotherms, both in the
laboratory and in the field, is an important first step toward predicting
the effects of climate change on their populations (Rozen-Rechels et al.,

2019). However, formulating a predictive model for the survival rate
involving several external factors can be challenging when the observed
data is limited. The present study proposed a hierarchical Bayesian
parametric survival modeling (B-PSM) framework that incorporates
expert opinions as supplementary information on temperature and
relative humidity effects on the survival time of female *I. ricinus* ticks. In
previous expert elicitation studies across various domains, the quantities
to elicit were primarily the parameters describing the probability dis-
tribution of the survival time, specifically the shape and scale parame-
ters for the Weibull distribution, without considering the effects of
covariates (Bousquet, 2010; Compare et al., 2017; Cope et al., 2019;
Singpurwalla, 1988). To our knowledge, this study was the first to
demonstrate an expert opinion elicitation framework on the survival
regression coefficients for the effects of multiple continuous covariates
and their interaction, explicitly in the ecological domain.

Despite having a small sample size, tick survival data were initially
assessed with frequentist parametric survival models (F-PSMs) to un-
cover the most suitable probability distribution and model formulation.
This preliminary step was important for guiding the formulation of
biological hypothesis, questions for expert opinion interviews, expert
opinion elicitation, and hierarchical Bayesian framework. As suggested
by the best-fitted F-PSM model, we designed the interview questions and
elicitation scheme to capture the expert opinions on the effects of rela-
tive humidity, temperature, and their ambiguous interaction. Therefore,
experts were asked to provide their opinions on 4 different conditions,
corresponding to the number of equations needed to solve a linear
algebraic system with 4 variables (β₀, β₁, β₂ and β₃) as in Eq. (8).

Experts are often defined as individuals with the relevance and extent of their experience in a topic of interest (Fazey et al., 2006). This definition can be rather subjective. In some domains, expert status can also be officially certified through a specified training program. In our context, however, a certified expert in the biology of *I. ricinus* ticks exposed to different environmental conditions does not exist. Therefore, we invited acarologists who have been handling *I. ricinus* colonies for years to provide their opinions principally based on their experience together with existing publications.

During the interview, the uncertainty that might arise was managed
by 1) Conducting the expert interview as an interactive online session.
The interactive conversation between experts and the interviewers was
made to reduce any linguistic uncertainty that would result in misun-
derstanding the questions; 2) Requesting the opinions on parameters
<figure><img src="image_8.png" /><figcaption>Fig. 5. Posterior distributions of the parameters *B* (β₀, β₁, β₂, and β₃), *k*, and *p*: Model 1, without expert opinions (A to F); Model 2, pooling expert opinions (G to L); Model 3, averaging expert opinions (M to R). Histograms represent posterior distributions. Solid lines indicate prior distributions. Vertical solid and dashed lines show median and the 95% credible interval, respectively.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Model 2: Pooling experts</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Model 3: Averaging experts</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>B Model 1: Without expert</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Model 2: Pooling experts</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Model 3: Averaging experts</figcaption></figure>

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /><figcaption>**Fig. 6.** Predictions for effects of relative humidity and temperature from 3 Bayesian models: A) Median survival time (days); B) Log of scale parameter of the Weibull distribution, lnλ. The predictions were calculated using the median values from posterior distributions. Values of relative humidity were truncated between 50% and 95%.</figcaption></figure>
that are easily understandable and do not require an advance statistical background. This practice has been addressed in most elicitation frameworks for survival modeling, as model parameters do not always have straightforward physical or biological interpretations (Bousquet, 2010; Compare et al., 2017; Cope et al., 2019; Singpurwalla, 1988). For example, we asked for the average survival time $\bar{T}_{e,c}$ instead of the parameter reflecting the effects of relative humidity on the survival time $\beta_1$. The expert opinions were later transformed to the desired unknown parameters during the elicitation process; 3) Allowing the experts to give their opinions on the conditions ($Q_{e,c}$; $U_{e,c}$) that are compatible with their experience. Forcing the experts to opine on unfamiliar pre-defined conditions would have created uncertainty and reduced the validity of their opinions; 4) Providing graphical representations of their answers, in real-time. Here, we displayed the probability density curves for the average survival time $\bar{T}_{e,c}$ and allowed the experts to adjust their answers until the final curves agreed with their opinions. This allows the experts to validate and avoid over/underestimating their opinions. An interactive web-based application displaying the Kaplan-Meier curve was also demonstrated to help with the expert opinion interview on the expected survival time of leukaemia patients (Cope et al., 2019); 5) Allowing the experts to recheck their opinions before submitting the final answers.

Several experts were recruited to avoid uncertainty and ultimately minimize the variability and identify potential outliers. We did not strive to form a consensus among the experts as reviewed by Kuhnert et al. (2010), neither through a consensus meeting (Cope et al., 2019) nor the Delphi process (MacMillan and Marshall, 2006) for the following reasons: 1) Each expert delivered their opinions based on different conditions and backgrounds 2) We avoided influencing the opinions across all experts. As a result, the variability among experts was addressed by pooling their opinions (Model 2) or forming a numerical consensus by averaging across all experts (Model 3). Medians of posterior distributions of each parameter were similar across all 3 Bayesian models. The pooling approach (Model 2) retains the diversity among the experts and allows them to be incorporated in the model with different weights. This approach could disregard the opinions of some experts that are irrelevant to the observed data (Experts 1, 3, 6). Nonetheless, the expert opinions chosen by the model (Experts 2, 4, and 5) did not agree entirely, particularly on the interaction between relative humidity and temperature $\beta_3$. The posterior distributions of Model 2 converged by compromising the different opinions on the parameters, resulting in the larger posterior distribution of $\beta_3$ than Models 1 and 3. While the averaging approach (Model 3) ignores the variability and the uncertainty that may arise from different experts (Albert et al., 2012). Among the 3 models, Model 3 generated the smallest posterior distributions of the parameter set indicating the effects of abiotic factors on tick survival $B$ ($\beta_0$, $\beta_1$, $\beta_2$, and $\beta_3$). However, the averaging approach may not be appropriate when the elicited opinions distinctly diverge across all the experts. Averaging multiple distributions with wide variations in modes could result in a consensus distribution with a large variance and a non-representative expected value.

The posterior distribution of the shape parameter $p$ converged to values greater than 1, indicating the mortality rate $\mu(t)$ increases over time (aging effects). As ixodid ticks only feed once per life stage, the energy stored in their bodies as lipid contents is limited during the off-host period and diminishes over time (Herrmann et al., 2013; Pool et al., 2017). The depletion of lipid resources is a critical limiting factor in their survival (Alasmari and Wall, 2021). As a result, aging ticks with lower energy reserves are more susceptible to death.

Posterior distributions of parameters $B$ reflect how abiotic factors (relative humidity, temperature) and their interaction affect the scale parameter $\lambda$ and the mortality rate $\mu(t)$. The parameter $\beta_1$ less than 0 indicated a protective effect of relative humidity on tick survival, while the parameter $\beta_2$ greater than 0 showed that temperature confers an increased risk of mortality. Interestingly, the parameter $\beta_3$ (interac-

tion effect) converged to $-\beta_2$, suggesting that higher temperature contributes higher scale parameter $\lambda$ and mortality rate $\mu(t)$ only when the relative humidity is low (Fig. 6B and Figure S7D to F). The predicted survival time, however, is predominately influenced by the relative humidity (Fig. 6A and Figure S7A to C). For example, despite the mortality rate being markedly different in dry conditions across temperature ranges, ticks are predicted to rapidly die out within a few days in all temperature conditions. In addition, the predicted tick survival at high humidity was greatly sensitive to a small change in the mortality rate $\mu(t)$ (Figure S7).

The models predicted that survival time would become notably shorter when the relative humidity dropped from 95% to 85% (Fig. 6A and Figure S7). This reflects the critical equilibrium humidity for adults *I. ricinus* of 86% – 96% (Lees, 1946), at which the relative humidity drops below ticks' ability to adsorb water vapor from unsaturated air (Needham and Teel, 1991). Posterior predictions of all 3 models reproduced well the survival data of female adult *I. ricinus* in most conditions of previous studies (Herrmann and Gern, 2010; Lees, 1946). Disagreement on the predictions was found in one condition reported briefly by Herrmann and Gern (2010), where our model overestimated the tick survival at $Q = 25$ °C and $U = 89\%$. In general, the relationship between temperature and performance of an organism, such as survival or locomotor, should exhibit an optimal condition (Rozen-Rechels et al., 2019). However, the historical data used in this study did not capture the optimal and lethal temperatures for *I. ricinus*. Therefore, the survival at high-temperature conditions could be overestimated by our model. In addition, the model prediction concisely captures the high sensitivity to the desiccation trait of *I. ricinus*. In comparison to other ixodid tick species, *I. ricinus* has a higher water loss rate. Based on the classification system of Hadley (1994), *I. ricinus*, along with *Ixodes reticulatus* and *Ixodes uriae* are classified as “mesic” species having a moderate water loss rate (0.8 – 2.0%/h). While *Amblyomma cajennense*, *Amblyomma maculatum*, *Amblyomma americanum*, *Hyalomma dromedarii*, *Dermacentor andersoni*, *Dermacentor variabilis*, *Dermacentor albipictus*, and *Rhipicephalus sanguineus* are “xeric” species with a low water loss rate (<0.8%/h) (Benoit and Denlinger, 2010). Besides, the significant interaction effects in our model suggested that the thermo-regulatory and hydro-regulatory systems of *I. ricinus* are not independent. The hydration status of ectotherms could modify the thermal sensitivity of cell and tissue metabolism, protect against thermal stress (Rozen-Rechels et al., 2019). Therefore, future studies on the impacts of abiotic factors and climate change on the survival and population dynamics of *I. ricinus* should consider the non-additive effects of temperature and relative humidity.

Finally, before applying our predictive model for a population dynamics study of *I. ricinus*, one should consider the following limitations: 1) The model is valid within the observed range of abiotic conditions ($Q$ between 5 °C and 25 °C, $U$ between 10% and 95%). Freezing or high temperatures could induce instant mortality (MacLeod, 1935); 2) The model parameters were based on limited historical data. The available methodology for observation in 1950 may have affected the data accuracy. Ticks in the present day could have adapted to survive in the current environment, which may be different from ticks in 1950. Moreover, ticks from different locations could also respond differently to the environment; 3) The mortality rates predicted from our model were based on ticks that were constantly exposed to the controlled environment. In the natural setting, however, they could change their behavior and avoid unfavorable conditions to improve their survival, e.g., moving toward humid air.

## 5. Conclusion

We demonstrated an expert opinion elicitation framework that integrates expert opinions as prior distributions for the effects of continuous explanatory variables on the survival process of *I. ricinus* ticks. Here, we summarized the key processes and considerations for applying
our framework to other contexts:

1) Defining the most suitable probability distribution for the survival data, such as the exponential, Weibull, or log-logistic distribution, is an important initial step. It not only defines the subsequent mathematical framework in the elicitation process, but also allows us to formulate biological hypotheses on the nature of mortality/hazard rate.

2) Subsequently, the mathematical relationship between the parameters to elicit (survival regression coefficients of the covariates) and the quantity to interview the experts should be well defined. To help the expert accurately provide their opinions, the quantity to interview should have a straightforward biological meaning that does not require high statistical background knowledge to understand. Also, the number of the parameters to elicit determines the number of questions required to ask the experts.

3) A flexible elicitation framework could assist the expert in providing their opinions on the survival time that corresponds to their experience while avoiding their opinions on unfamiliar conditions.

4) Using an interactive web-based application during the interview could visually assist the experts and allow them to revise their opinions.

5) With multiple experts, we demonstrated two approaches to combining opinions: pooling (Model 2) and averaging (Model 3). The pooling approach considers the variations in expert opinions, but it may also disregard some opinions that are irrelevant to the observed data. The averaging approach, on the other hand, simplifies the calculation by achieving a numerical consensus of the opinions, but it may be less informative when the opinions distinctly diverge.

Our model predictions also highlighted the importance of the combined effects of relative humidity and temperature on the survival of *I. ricinus* ticks. Although the survival of *I. ricinus* is deemed to be more dependent on relative humidity, the historical data used in our study did not include the upper and lower lethal temperature ranges. Therefore, additional studies on the survival of *I. ricinus* ticks, involving more sample size and a wider range of experimental conditions, is still required to improve our understanding of the impacts of climate change on the dynamics of their populations and Lyme borreliosis. Finally, when the available survival data of ectoherms in several abiotic conditions are limited, the elicitation framework proposed in this study could be applied to help acquire and incorporate expert opinions to develop predictive survival models in other ecological contexts.

## Authors' contributions

Conception and design of study: PW, SB, SD, KCM; Acquisition of data: PW, SB, KCM; Analysis and/or interpretation of data: PW, SB, SD, KCM; Drafting the manuscript: PW, SB, KCM; Revising the manuscript critically for important intellectual content: PW, SB, SD, TH, FB, KCM; Approval of the version of the manuscript to be published: PW, SB, SD, TH, FB, KCM.

## Data availability

R codes related to the Shiny-based web application and hierarchical Bayesian analysis are available at http://doi.org/10.5281/zenodo.4569493.

## Declaration of Competing Interest

The authors whose names are listed immediately below certify that they have NO affiliations with or involvement in any organization or entity with any financial interest (such as honoraria; educational grants; participation in speakers' bureaus; membership, employment, consulancies, stock ownership, or other equity interest; and expert testimony

or patent-licensing arrangements), or non-financial interest (such as personal or professional relationships, affiliations, knowledge or beliefs) in the subject matter or materials discussed in this manuscript.

## Acknowledgement

The authors would like to express our gratitude to all the experts who participated in our study. Also, PW is grateful for a grant from the Association Nationale de la Recherche et de la Technologie, bourse CIFRE (n° 2019/0393).

## Supplementary materials

Supplementary material associated with this article can be found, in the online version, at doi:10.1016/j.ecolmodel.2021.109821.

## Reference

Alasmari, S., Wall, R., 2021. Metabolic rate and resource depletion in the tick *Ixodes ricinus* in response to temperature. Exp. Appl. Acarol. 83, 81–93. https://doi.org/10.1007/s10493-020-00568-1.

Albert, I., Donnet, S., Guihenneuc-Jouyaux, C., Low-Choy, S., Mengersen, K., Rousseau, J., 2012. Combining expert opinions in prior elicitation. Bayesian Anal 7, 503–532. https://doi.org/10.1214/12-BA717.

Bates, A.E., Pecl, G.T., Frusher, S., Hobday, A.J., Wernberg, T., Smale, D.A., Sunday, J. M., Hill, N.A., Dulvy, N.K., Colwell, R.K., Holbrook, N.J., Fulton, E.A., Slawinski, D., Feng, M., Edgar, G.J., Radford, B.T., Thompson, P.A., Watson, R.A., 2014. Defining and observing stages of climate-mediated range shifts in marine systems. Glob. Environ. Chang. 26, 27–38. https://doi.org/10.1016/j.gloenvcha.2014.03.009.

Benoit, J.B., Denlinger, D.L., 2010. Meeting the challenges of on-host and off-host water balance in blood-feeding arthropods. J. Insect Physiol. 56, 1366–1376. https://doi.org/10.1016/j.jinsphys.2010.02.014.

Berg, M.P., Toby Kiers, E., Driessen, G., van der Heijden, M., Kooi, B.W., Kuenen, F., Liefting, M., Verhoef, H.A., Ellers, J., 2010. Adapt or disperse: understanding species persistence in a changing world. Glob. Chang. Biol. 16, 587–598. https://doi.org/10.1111/j.1365-2486.2009.02014.x.

Bojke, L., Soares, M.F.O., Fox, A., Jankovic, D., Claxton, K.P., Morton, A., Sharples, L., Jackson, C.H., Taylor, A.C., Colson, A., 2019. Developing a reference protocol for expert elicitation in healthcare decision making. Health Technol. Assess. Rep 10.3310/hta25370.

Bousquet, N., 2010. Elicitation of weibull priors. arXiv Prepr. arXiv1007.4740 1–21.

Chang, W., Cheng, J., Allaire, J., Xie, Y., McPherson, J., 2020. Shiny: Web application Framework For R. R Packag. Version 1.5.0.

Colson, A.R., Cooke, R.M., 2018. Expert elicitation: using the classical model to validate experts' judgments. Rev. Environ. Econ. Policy 12, 113–132. https://doi.org/10.1093/reep/rex022.

Compare, M., Baraldi, P., Bani, I., Zio, E., Mc Donnell, D., 2017. Development of a Bayesian multi-state degradation model for up-to-date reliability estimations of working industrial components. Reliab. Eng. Syst. Saf. 166, 25–40. https://doi.org/10.1016/j.rres.2016.11.020.

Cope, S., Ayers, D., Zhang, J., Batt, K., Jansen, J.P., 2019. Integrating expert opinion with clinical trial data to extrapolate long-term survival: a case study of CAR-T therapy for children and young adults with relapsed or refractory acute lymphoblastic leukemia. BMC Med. Res. Methodol. 19, 1–11. https://doi.org/10.1186/s12874-019-0823-8.

Dumic, I., Severini, E., 2018. Ticking bomb": the impact of climate change on the incidence of lyme disease. Can. J. Infect. Dis. Med. Microbiol. 2018 https://doi.org/10.1155/2018/5719081.

Ergon, T., Borgan, Ø., Nater, C.R., Vindenes, Y., 2018. The utility of mortality hazard rates in population analyses. Methods Ecol. Evol 10.1111/2041-210X.13059.

Fazey, I., Fazey, J.A., Salisbury, J.G., Lindenmayer, D.B., Dovers, S., 2006. The nature and role of experiential knowledge for environmental conservation. Environ. Conserv. 33, 1–10. https://doi.org/10.1017/S037689290600275X.

Giannini, T.C., Costa, W.F., Cordeiro, G.D., Imperatriz-Fonseca, V.L., Saraiva, A.M., Biesmeijer, J., Garibaldi, L.A., 2017. Projected climate change threatens pollinators and crop production in Brazil. PLoS ONE 12, 1–13. https://doi.org/10.1371/journal.pone.0182274.

Hadley, N.F., 1994. Water Relations of Terrestrial Arthropods. Academic Press, New York.

Halsch, C.A., Shapiro, A.M., Fordyce, J.A., Nice, C.C., Thorne, J.H., Waetjen, D.P., Forister, M.L., 2021. Insects and recent climate change. Proc. Natl. Acad. Sci. U. S. A. 118, 1–9. https://doi.org/10.1073/PNAS.2002543117.

Herrmann, C., Gern, L., 2010. Survival of *Ixodes ricinus* (Acari: Ixodidae) under challenging conditions of temperature and humidity is influenced by *Borrelia burgdorferi* sensu lato infection. J. Med. Entomol. 47, 1196–1204. https://doi.org/10.1603/me10111.

Herrmann, C., Voordouw, M.J., Gern, L., 2013. *Ixodes ricinus* ticks infected with the causative agent of Lyme disease, *Borrelia burgdorferi* sensu lato, have higher energy reserves. Int. J. Parasitol. 43, 477–483. https://doi.org/10.1016/j.ijpara.2012.12.010.
Hiance, A., Chevret, S., Vincent, L., 2009. A practical approach for eliciting expert prior beliefs about cancer survival in phase III randomized trial. J. Clin. Epidemiol. 62, 431–437. https://doi.org/10.1016/j.jclinepi.2008.04.009 e2.

Johnson, F.A., Smith, B.J., Bonneau, M., Martin, J., Romagosa, C., Mazzotti, F., Waddle, H., Reed, R.N., Kettevrlin, J., Vitt, L.J., 2017. Expert elicitation, uncertainty, and the value of information in controlling invasive species. Ecol. Econ. 137, 83–90. https://doi.org/10.1016/j.ecolecon.2017.03.004.

Kleinbaum, D., Klein, M., 2012. Parametric Survival Models. Survival Analysis: A Self-Learning Text. Springer., pp. 289–362. https://doi.org/10.1007/978-1-4419-6646-9_7

Krueger, T., Page, T., Hubacek, K., Smith, L., Hiscock, K., 2012. The role of expert opinion in environmental modelling. Environ. Model. Softw. 36, 4–18. https://doi.org/10.1016/j.envsoft.2012.01.011.

Kuhnert, P.M., Martin, T.G., Griffiths, S.P., 2010. A guide to eliciting and using expert knowledge in Bayesian ecological models. Ecol. Lett. 13, 900–914. https://doi.org/10.1111/j.1461-0248.2010.01477.x.

Lee, H., Kim, J.E., Lee, S., Lee, C.H., 2018. Potential effects of climate change on dengue transmission dynamics in Korea. PLoS ONE 13, e0199205.

Lees, A.D., 1946. The water balance in *Ixodes ricinus* L. and certain other species of ticks. *Parasitology* 37, 1–20. https://doi.org/10.1017/S003118200013093.

MacLeod, J., 1935. *Ixodes ricinus* in relation to its physical environment: II. The factors governing survival and activity. *Parasitology* 27, 123–144.

MacMillan, D.C., Marshall, K., 2006. The Delphi process – an expert-based approach to ecological modelling in data-poor environments. Anim. Conserv. 9, 11–19. https://doi.org/10.1111/j.1469-1795.2005.00001.x.

Makowski, D., Ben-Shachar, M., Lüdecke, D., 2019. bayestestR: describing effects and their uncertainty, existence and significance within the Bayesian framework. J. Open Source Softw. 4, 1541. https://doi.org/10.21105/joss.01541.

Milne, A., 1950. The ecology of the sheep tick, *Ixodes ricinus* L.:microhabitat economy of the adult tick. *Parasitology* 40, 14–34. https://doi.org/10.1017/S0031182000017820.

Nadeau, C.P., Urban, M.C., Bridle, J.R., 2017. Climates past, present, and yet-to-come shape climate change vulnerabilities. *Trends Ecol. Evol.* 32, 786–800. https://doi.org/10.1016/j.tree.2017.07.012.

Needham, G.R., Teel, P.D., 1991. Off-host physiological ecology of ixodid ticks. *Annu. Rev. Entomol.* 36, 659–681.

Omurlu, I.K., Ozdamar, K., Ture, M., 2009. Comparison of Bayesian survival analysis and Cox regression analysis in simulated and breast cancer data sets. Expert Syst. Appl. 36, 11341–11346. https://doi.org/10.1016/j.eswa.2009.03.058.

Omurlu, I.K., Ture, M., Ozdamar, K., 2015. Bayesian analysis of parametric survival models: a computer simulation study based informative priors. J. Stat. Manag. Syst. 18, 405–423. https://doi.org/10.1080/09720510.2014.961763.

Pecl, G.T., Araújo, M.B., Bell, J.D., Blanchard, J., Bonebrake, T.C., Chen, I.-C., Clark, T. D., Colwell, R.K., Danielsen, F., Evengård, B., Falconi, L., Ferrier, S., Frusher, S., Garcia, R.A., Griffis, R.B., Hobday, A.J., Janion-Scheepers, C., Jarzyna, M.A., Jennings, S., Lenoir, J., Linnetved, H.I., Martin, V.Y., McCormack, P.C., McDonald, J., Mitchell, N.J., Mustonen, T., Pandolfi, J.M., Pettorelli, N., Popova, E., Robinson, S.A., Scheffers, B.R., Shaw, J.D., Sorte, C.J.B., Strugnell, J.M., Sunday, J. M., Tuanmu, M.-N., Vergés, A., Villanueva, C., Wernberg, T., Wapstra, E., Williams, S.E., 2017. Biodiversity redistribution under climate change: impacts on ecosystems and human well-being. Science (80..) 355, eaai9214. https://doi.org/10.1126/science.aai9214.

Plummer, M., 2019. *Rjags: Bayesian Graphical Models Using MCMC. R Packag*, pp. 4–10 version.

Plummer, M., 2003. JAGS: a program for analysis of Bayesian graphical models using Gibbs sampling, in: Proceedings of the 3rd International Workshop On Distributed Statistical Computing. pp. 1–10 10.1002/ana.1067.

Pool, J.R., Petronglo, J.R., Falco, R.C., Daniels, T.J., 2017. Neurobiology, physiology, biochemistry energy usage of known-age blacklegged ticks (acari : ixodidae): what is the best method for determining physiological age ? 54, 949–956 10.1093/jme/tjx042.

R Core Team, 2019. A language and environment for statistical computing. R. Found. Stat. Comput.

Román-Palacios, C., Wiens, J.J., 2020. Recent responses to climate change reveal the drivers of species extinction and survival. Proc. Natl. Acad. Sci. U. S. A. 117, 4211–4217. https://doi.org/10.1073/pnas.1913007117.

Rozen-Rechels, D., Dupoué, A., Lourdais, O., Chamaillé-Jammes, S., Meylan, S., Clobert, J., Le Galliard, J.F., 2019. When water interacts with temperature: ecological and evolutionary implications of thermo-hydroregulation in terrestrial ectotherms. Ecol. Evol. 9, 10029–10043. https://doi.org/10.1002/ece3.5440.

Schmalensee, L.Von, Hulda, K., Joacim, G., 2021. Thermal performance under constant temperatures can accurately predict insect development times across naturally variable microclimates. Ecol. Lett. 24, 1633–1645. https://doi.org/10.1111/ele.13779.

Singpurwalla, N.D., 1988. An interactive PC-based procedure for reliability assessment incorporating expert opinion and survival data. J. Am. Stat. Assoc. 83, 43–51. https://doi.org/10.1080/01621459.1988.10478563.

Therneau, T., 2020. A Package for Survival Analysis in R. R Packag. version 3.2-7.

Therneau, T.M., Grambsch, P.M., 2000. Modeling Survival Data: Extending the Cox Model. Springer, New York.

Walsh, D.P., Norton, A.S., Storm, D.J., Van Deelen, T.R., Heisey, D.M., 2018. Using expert knowledge to incorporate uncertainty in cause-of-death assignments for modeling of cause-specific mortality. Ecol. Evol. 8, 509–520. https://doi.org/10.1002/ece3.3701.