Individual behavioral type captured by a Bayesian model comparison of cap making by sponge crabs

Animal collection

From December 2015 to April 2017, 40 individuals (21 males, 19 females) of Lauridromia dehaani (Brachyura: Dromiidae) (Fig. 1A) were captured using a gill net at the Sakai fishing port, Minabe town, Wakayama, Japan (33°44′N, 135°20′E). We conducted behavioral experiments of cap making on 38 individuals (20 males, 18 females) and video recorded the behaviors of 2 individuals (4.30 cm and 7.19 cm of the carapace widths for each) in a tank filled with filtered natural seawater (about 3.4 % of the salinity) at Shirahama Aquarium, Seto Marine Biological Laboratory, Kyoto University (33°41′N, 135°20′E), from December 2015 to June 2017. For the behavioral experiments, we successfully sampled 8 individuals repeatedly (3 or 5 times for each). Thus, we only sampled one observation from one of the other 30 individuals. Note that although the sample sizes of the behavioral acts for each individual are different our method is still applicable. Before the experiments, all individuals were retained in the tanks (19.5–23.8 °C, light on: 8–17, light off: 17–8) of the aquarium for more than two days for acclimation. We measured their carapace width (cm) (Fig. 1B) as a proxy for the body size, and divided them into three levels depending on whether they lacked any of the fourth and fifth pereiopods: (O) none of the fourth and fifth pereiopods were absent, (1) one of them was absent, (2) both fourth and fifth of each side were absent.

Experimental setup and procedure

We prepared three sizes of white melamine foam that is found commonly worldwide (most notably manufactured by BASF of Germany) and often used in general households (in Japan, it is called Gekiochi-kun, LEC, Inc.) (S: 20 mm × 30 mm × 40 mm; M: 30 mm × 60 mm × 85 mm; L: 30 mm × 140 mm × 150 mm). We used this sponge because it is easy to sink.

First, to confirm that the cap making behavior of the crab Lauridromia dehaani is similar to the behavior in the reports (Dembowska, 1926; McLay, 1983), we video-recorded the behavior from two crabs. The two crabs were used only for video recording in the aquarium (310 mm × 180 mm × 240 mm, W × L × H). We started the recording from 9 am to 10 am in the morning, and stopped for 2 h after the crabs stopped cap making. We used a red light transmitted through a polyvinyl chloride board and excluded another light source by enclosing the aquarium. We obtained 5 recordings for each crab.

Secondly, we performed experiments on choice of cap size, trimming, and excavating behaviors. S size sponge was smaller than all crabs, whereas the L size was larger than all crabs. Each sponge was put pseudo-randomly on either side and at the back center of the cage (700 mm × 470 mm × 190 mm, W × L × H, Fig. 1C), which floated in the tank. Then, crabs were introduced to the front center of a cage floating in the tank, thereby the distance between each sponge and the crab was equal. We started a trial from 9 am to 10 am and checked whether the crab carried any sponge once a day. We counted the days when the crab carried a sponge. The latency, measured as the number of days to produce a cap, is modeled as a random variable. Note that the crab was assumed to make a cap at night, because it is considered nocturnal (McLay, 1983). If it did, we collected the sponge; otherwise, the crab and the three sponges remained in the cage. When the crab did not carry any sponge for five days, we stopped the trial. We destroyed all sponges that the crabs processed and measured their whole area (cm²), and area of the concave part (cm²) excavated by a crab from the pictures taken 46 cm above the sponges. The trimmed area and concave area are modeled as random variables. Although in the beginning we only performed one trial for one individual (N_animal = 30), we obtained five trials for one individual after February 2017 (N_animal = 8) to examine the behavior type. Our hypothesis that a behavioral type would be formed was conceived after the day, but we did not change the experimental condition.

Statistical modeling

To quantify the behavioral type in the experiment, we constructed 26 statistical models (Table 1) for the four different aspects: (1) choice of sponge size (6 models), (2) amount of sponge trimmed by cutting (8 models), (3) size of cavity (6 models), and (4) latency to produce a cap (6 models). In each case, we built the statistical models specifying individual behavioral types as hierarchical structures with parameters and performed MCMC samplings from the posterior distribution. Also, we conditioned the variables with the carapace width, levels of leg absence, and gender. We specified the models in the probabilistic programming language Stan (Stan Development Team, 2018). We used non-informative uniform priors for some parameters unless otherwise explicitly described. The performed samplings from the posterior distributions using No-U-Turn Sampler (NUTS) implemented as a Hamiltonian Monte Carlo (HMC) sampler in Stan. Sampling convergence was visually diagnosed by trace plots and quantitatively via the Gelman-Rubin convergence statistic, R_hat (Gelman & Rubin, 1992). All sampled draws were judged to be converged when R_hat < 1.10, were used to construct predictive distributions with WAIC on each model. All computations were performed in the R statistical environment (R Core Team, 2018), and the Stan codes for each model were compiled and executed through the R package rstan (Stan Development Team, 2018).

We compared the predictive performances of the models using WAIC (Watanabe, 2018; Watanabe, 2010b). It should be emphasized that the WAIC of a hierarchical model can be defined in several ways depending on how a predictive distribution is defined. In our case, as we would like to construct a new distribution regarding a new act of a new individual, we have to marginalize the intermediate parameters assigned to each individual in the statistical model (Watanabe, 2018). This is because we are interested in the prediction of a new behavioral act when we get a new individual and a new behavioral act instead of the prediction of a new behavioral act from the individuals sampled in this study. By performing this procedure, we can equally compare a hierarchical model with a non-hierarchical model, because the focus of the prediction in a non-hierarchical model is on a new behavioral act of a new individual.

Here we briefly describe the basic procedure based on Watanabe (2018). Let Xⁿ = (X₁, …, X_n) be an sample from the unknown true distribution and p(x|w) a statistical model with w assigned to each individual. Furthermore, w is assumed to be taken from φ(w|w₀) and w_0 from φ (w_0) to form a hierarchical structure. In learning step, w₁, …, w_Nanimal is prepared. In prediction step, our statistical model is built like p_model(x|w₀) by marginalizing w out: (1)${\displaystyle p_{model}\left(x|w_0\right)=\int p\left(x|w\right)\phi \left(w|w_0\right)dw.}$

WAIC is a measure for the degree of accuracy of an approximation of a predictive distribution to the true distribution generating data. For our hierarchical model, the predictive distribution is defined as 𝔼_w0[p_model(x|w₀)]. Then, the marginal-level WAIC for a hierarchical model is defined as: (2)${\displaystyle WAI{C}_h=-\frac{1}{N}\sum _{i=1}^{N}log{\mathbb{E}}_{w_0}\left[p_{model}\left(X_i|w_0\right)\right]+\frac{1}{N}\sum _{i=1}^N{\mathbb{V}}_{w_0}\left[log{p}_{model}\left(X_i|w_0\right)\right]}$

where 𝔼_w0[] and 𝕍_w0[] are the average and variance operators, respectively, of the posterior distribution of w₀. w₀ is estimated so that practically the MCMC sample is used, thus numerical integration is required. In this study, the computation is implemented in the ‘function’ block in the Stan codes using the Simpson’s rule for the numerical integration and the log_sum_exp function provided in Stan (see Supplemental Information 3).

On the contrary, the WAIC for a non-hierarchical model is defined for a statistical model p_model(x|w): (3)${\displaystyle WAI{C}_{non-h}=-\frac{1}{N}\sum _{i=1}^{N}log{\mathbb{E}}_w\left[p_{model}\left(X_i|w\right)\right]+\frac{1}{N}\sum _{i=1}^N{\mathbb{V}}_w\left[log{p}_{model}\left(X_i|w\right)\right]}$

where 𝔼_w[] and 𝕍_w[] are the average and variance operators, respectively, of the posterior distribution of w. Note that the often used conditional-level WAIC is described in the Discussion.

Choice of material size (model 1_1)

To provide an overview of the specified models, we describe only the best-performing models in terms of WAIC here. The other models are summarized in Table 1.

We formulate a tendency toward a choice as μ[n, m] ( m = 1, 2, 3 for M, L, skip, respectively):

(4)${\displaystyle \mu \left[n,1\right]=0,}$ (5)${\displaystyle \mu \left[n,2\right]=a_{choic{e}_L}\left[ID\left[n\right]\right]+b_{choic{e}_L}\cdot CarapaceWidth\left[n\right]+c_{choic{e}_L}\cdot LegLack\left[n\right],}$ (6)${\displaystyle \mu \left[n,3\right]=d_{choic{e}_0}+e_{choic{e}_0}\cdot CarapaceWidth\left[n\right]+f_{choic{e}_0}\cdot LegLack\left[n\right],n=1,\dots ,N_{act}}$

where N_act = 68 is the total number of behavioral acts, and ID represents animal identity (from 1 to N_animal =38). μ is linked to the linear predictor in terms of the carapace width, CarapaceWidth and the level of absence of leg, LegLack. The choice of an M size is fixed to zero. a_choiceL is for each individual, thus it is hierarchized. d_choice0 is not hierarchized. The distribution of a_choiceL is defined as the normal distribution with the mean a_choiceL0 and standard deviation a_choiceLs:

(7)${\displaystyle a_{choic{e}_L}\left[k\right]\sim Normal\left(a_{choic{e}_{L0}},a_{choic{e}_{Ls}}\right),k=1,\dots ,N_{animal}.}$

The actual choice Choice is defined as the categorical distribution with the softmax function:

(8)${\displaystyle Choice\left[n\right]\sim Categorical\left(softmax\left(\mu \left[n,\right]\right)\right),n=1,\dots ,N_{act}.}$

Thus, in this case, a statistical model p_model(x|w₀) is set using the parameters: (9)${\displaystyle w_0=\left(a_{choic{e}_{L0}},a_{choic{e}_{Ls}},b_{choic{e}_L},c_{choic{e}_L},d_{choic{e}_0},e_{choic{e}_0},f_{choic{e}_0}\right).}$

Note that w = a_choiceL is marginalized out when we build the predictive distribution, so that it is not included in w₀. The choice Choice[n] is modeled as a random variable X_n. CarapaceWidth[n] and LegLack[n] are the conditioning variables.

Cap making using an artificial sponge

The behaviors of the two crabs were video recorded to confirm the Cap making behavioral sequence when using an artificial sponge (Fig. S1). The crabs grasped either side of the sponge using their second and third pereiopods, and trimmed small pieces of the sponge using their chelae (Fig. 2A upper left, upper right, Video S1). They first visited the two sides of the sponge. To make a cavity, the crabs rotated their bodies backward and grasped the sponge by the fourth and fifth pereiopods. By repeating these behaviors, the crabs made a groove to cut off a portion of a sponge. In 9 of the 10 trials conducted, it took about 50 min to cut the portion, and the crabs started excavating as soon as they finished the trimming behavior. In the other trial, it took 19 min.

Next, the crabs made cavities by tearing off small pieces of a sponge (Fig. 2A bottom, Video S2). It took 11 min on average to excavate the cavity. Then, the crabs rotated their bodies backward in order to catch the excavated sponges with the fourth and fifth pereiopods while they kept the portion grasped by the second and third pereiopods. Finally, the crabs released the second and third pereiopods from the cap and carried it off (Figs. 2B, 2C). In terms of behavior, the crabs often rotated their bodies forward, dorso-ventrally, to enlarge the cavity. It is rare for them to move laterally. They repeated the excavation activities up to eleven times per night and it took up to 4 h. When the crabs rotated their bodies, the direction of rotation was maintained along with the sponge. While the crabs cut the sponge, they actively moved around the sponge. In contrast, they persistently stayed under the sponge during excavation.

Sponge choice

None of the 38 animals chose the S size sponge, and 7 animals skipped the cap making behavior (Fig. 3A). Therefore, we defined the choice as a random variable taking either of the three values: M, L, or skipping. The hierarchical model assuming behavioral types 1_1 (Figs. 3A, 3B) outperformed the non-hierarchical one in terms of WAIC (2.99 nat/sample in the difference, Figs. 3A–3D, Table 1). The posterior probability of the behavioral choices was more widely variable on model 1_1 than in model 1_6 depending on the individual difference specified as a_choiceL (Fig. 3B). To show the variability, the probability of choice sampled from the posterior distribution from the highest density is visualized in white lines (Figs. 3A, 3C). Note that the variability of the choice probability in the white curved lines is smaller than the model 1_1 even if the number of lines are the same. Although the body size of the animal indicated with the white arrowhead (Fig. 3A) is small, it preferably selected the size L. This indicates a large inter-individual difference. In the case of either the hierarchical or the non-hierarchical model, the behavioral choice of the sponges was better explained by the carapace width (Figs. 3A, 3C; Table 1). The estimated information gained by the model 1_1 against model 1_3 is 1.35 nat/sample (Table 1). This suggests that larger crabs tend to choose L size sponge rather than M size. However, the crabs larger than about nine cm carapace width did not choose the sponges.

Trimming

After a choice of either the M or L size sponges, the crabs decided whether to trim the extra parts of the sponges (Figs. 4A–4D). Here, we modeled the size of an area in a sponge that was trimmed (N_animal = 30). The trimming behavior for the sponge followed three patterns (Fig. 4D). They cut off (1) all four corners of a sponge, (2) one corner of the sponges elliptically, or (3) two corners of the sponge linearly. The crabs trimmed the white area (Fig. 4D) and started excavating a cavity (Fig. 4C Twenty-three crabs skipped the trimming behavior in 33 trials.

After trimming or skipping, they started excavating. For the behavioral act of trimming, a non-zero data point indicating a trimmed size of the sponge was recorded (Figs. 4A–4C). The size decreased with the increase of the carapace width. If a crab skipped trimming, a data point was recorded at zero, meaning no trimming. Almost all the crabs that chose the M size sponges decided not to trim the sponge except for one individual. Meanwhile, they trimmed less amounts of the sponges relative to the increase of their body sizes when they chose the L size sponges.

The WAIC of the hierarchical model 2_1 was −2.08 and that of comparable non-hierarchical model 2_6 was 7.40 (9.48 nat/sample in difference, Figs. 4A, 4D, Table 1), indicating that the hierarchical model is significantly better than the non-hierarchical one.

Cavity size

Six crabs just cut the sponge and did not excavate the sponge. We modeled the size of a cavity in a cap (N_animal = 30) as a random variable taken from the gamma distribution with the log link function (Fig. 5). The size increased with the carapace width, and the model considering individual behavioral types performed best (Table 1). The WAIC of the hierarchical model 3_1 is slightly smaller than that of the comparable non-hierarchical model 3_2 (0.08 nat/sample in the difference; Figs. 5A, 5B, Table 1). The individual with the arrowhead (Fig. 5A) made relatively large cap cavities, indicating large inter-individual differences. As expected, larger crabs made larger cavities. The difference of the WAIC was about 0.1 (Fig. 5B). The predictability improvement is relatively small against that of sponge choice, suggesting that the individual behavioral type would have a lower effect in the determination of cavity size.

Latency

We modeled the latency for cap making (from the choice of sponge to carrying) by 32 crabs as a random variable taken from the zero-inflated distribution (Fig. 6). No obvious relation was found between the carapace width and the latency, and the number of crabs that had carried the cap by the next day. However, the hierarchical model 4_1 outperformed the non-hierarchical model 4_2 (WAIC values, 1.10 and 1.28 respectively, thus 0.18 nat/sample in the difference).

Statistical modeling from the predictive point of view

First, as a general theory, we state that the model construction for prediction is appropriate in our case because of the difference between prediction and discovery of true distribution. Second, from a mathematical point of view, we explain the maximum likelihood (ML) estimation and Bayesian inference, which are typical methods of statistical model construction, and discuss the validity of Bayesian inference and WAIC. Lastly, we consider the novelty of our statistical modeling in terms of the type of prediction and the difference between conditioned-level WAIC and marginal-level WAIC.

Statistical models are broadly divided into those for prediction and those for discovering the true distribution. The two situations are as follows; (a) when there is no distribution that generates data in the finite set of models under consideration, (b) when there is a true distribution generating data in the finite set of models considered. Each situation is formulated as follows. Let n be the sample size.

In case (a), the model is constructed by minimizing the generalization loss G_n:

(25)${\displaystyle G_n=-log\int q\left(x\right)log{p}^{\ast }\left(x\right)dx,}$ where q(x) is the true distribution of the data and p^∗(x) is the predictive distribution. Predictive distribution is a probability distribution of the unknown new data based on the training data. For instance, the predictive distribution by ML estimation $p_{ML}^{\ast }\left(x\right)$ is defined by

(26)${\displaystyle p^{\ast }\left(x\right)=p\left(x|\hat{w}\left(Data\right)\right),}$ the likelihood model p(x|w) whose parameter w is plugged in the ML estimator (MLE) $\hat{w}\left(Data\right)$. Bayesian predictive distribution is defined by the expectation of p(x|w) overall a posterior distribution p_posterior(w|Data): (27)${\displaystyle p^{\ast }\left(x\right)=\int p\left(x|w\right)p_{posterior}\left(w|Data\right)dw={\mathbb{E}}_w\left[p\left(x|w\right)\right].}$

The generalization loss G_n is decomposed into two terms as below. The first term is independent from the model and the second one is the difference from the true distribution and the predictive one, called Kullback-Leibler divergence: (28)${\displaystyle G_n=-log\int q\left(x\right)log{p}^{\ast }\left(x\right)dx}$ (29)${\displaystyle =-log\int q\left(x\right)logq\left(x\right)dx+\int q\left(x\right)log\frac{q\left(x\right)}{p^{\ast }\left(x\right)}dx.}$

Thus, the generalization loss quantifies how far a predicted distribution is from the true distribution.

In the case (b), the model is constructed by minimizing the negative log marginal likelihood F_n: (30)${\displaystyle F_n=-log{Z}_n,}$

where p_likelihood(Data|w) is a likelihood, φ(w) is a prior distribution, and Z_n = ∫p_likelihood(Data|w)φ(w)dw is the marginal likelihood. Note that Z_n is equal to the normalizing constant of the posterior distribution since (31)${\displaystyle p_{posterior}\left(w|Data\right)=\frac{p_{likelihood}\left(Data|w\right)\phi \left(w\right)}{Z_n}}$

holds by the definition of conditional probability. In statistical mechanics, F_n is called ‘free energy’. Because of Data = (x₁, …, x_n), the true distribution of a dataset is $Q\left(x_1,\dots ,x_n\right)={\prod }_{i=1}^{n}q\left(x_i\right)$. For simplicity, we put xⁿ = (x₁, …, x_n). We rewrite the marginal likelihood Z_n = Z(xⁿ) to emphasize that it is a probability distribution of a dataset. We consider the expectation of F_n overall Q(xⁿ). Just like the generalization loss, it is decomposed into the model independent term and the Kullback-Leibler divergence from Q(xⁿ) to Z(xⁿ): (32)${\displaystyle \int Q\left(x^n\right)F_{n}d{x}^n=-\int Q\left(x^n\right)logZ\left(x^n\right)d{x}^n}$ (33)${\displaystyle =-\int Q\left(x^n\right)logQ\left(x^n\right)d{x}^n+\int Q\left(x^n\right)log\frac{Q\left(x^n\right)}{Z\left(x^n\right)}d{x}^n.}$

Thus, the expected free energy quantifies how different the marginal likelihood is from the true distribution. By definition, Z(xⁿ) is the distribution of a dataset, constructed by the model. In a practical application, a dataset is only an obtained set, hence F_n is to be minimized.

For the above formulation, the following facts are known in statistics. In case (a), the ‘effective method’ would be a model selection method which can select the model while minimizing the generalization loss if n → ∞. Obviously, minimizing G_n is the ‘effective method’. Moreover, minimizing WAIC is also ‘effective’.

On the other hand, in case (b), a model selection method that can select the true distribution if n → ∞, is called ‘consistent method’. Although statistical modeling in case (b) is formulated by minimizing the expected free energy ∫Q(xⁿ)F_ndxⁿ, minimizing the free energy F_n is ‘consistent’.

Our statistical analysis considers (a). We argue that no distribution generates data in the finite set of models, because our models are descriptive; they are not mechanistic models which represent individual behaviors of the crabs. Hence, it is appropriate to construct predictive models.

Next, we explain estimation methods. ML estimation and Bayesian inference are the typical methods used. However, they can be understood in a unified framework.

An analyst arbitrarily designs the simultaneous distribution of a pair (x, w) (34)${\displaystyle p\left(x,w\right)=p\left(x|w\right)\phi \left(w\right),}$

where x is an observable variable and w is a latent variable. In the ML estimation, the existence of the MLE $\hat{w}$ is assumed and φ(w) is set to $\delta \left(w-\hat{w}\right)$ whose realization is limited to the MLE $\hat{w}$. As a formality, this can be interpreted as the parameter w, which is the MLE $\hat{w}$.

Let w be a real number. For an arbitrary real number a, the function δ(w − a) satisfies

(35)${\displaystyle \delta \left(w-a\right)=\left\{\begin{array}{c}\infty w=a\hfill \\ 0w\ne a\hfill \end{array}\right.\text{and}\int \delta \left(w-a\right)dw=1.}$

Thus, it is clear that $p\left(Data|w\right)\delta \left(w-\hat{w}\right)∕Z_n=\delta \left(w-\hat{w}\right)$ holds.

In Bayesian inference, φ(w) is a prior distribution. ML estimation can be described in the same way as the Bayesian inference; ML estimation is the case when the prior distribution is fixed to $\delta \left(w-\hat{w}\right)$. Accordingly, ML estimation and Bayesian inference can be understood in a unified way for the parameter estimation. The method to be used depends on the purpose of statistical modeling.

In the construction of a predictive distribution, the following theorem has been proven (Watanabe, 2018). Let $G_n^{ML}$ be the generalization loss of ML estimation and $G_n^{Bayes}$ be the one of Bayesian inference. The overall expectation symbol across the dataset is denoted by 𝔼[⋅] = ∫Q(xⁿ)[⋅]dxⁿ. There are the constants μ, λ (0 < λ < μ), which are dependent on a model and the true distribution such that

(36)${\displaystyle \mathbb{E}\left[G_n^{ML}\right]=S+\frac{\mu }{n}+o\left(\frac{1}{n}\right),}$ (37)${\displaystyle \mathbb{E}\left[G_n^{Bayes}\right]=S+\frac{\lambda }{n}+o\left(\frac{1}{n}\right)}$

hold, where S =  − ∫q(x)logq(x)dx is the entropy. Especially, when a prior distribution in Bayesian inference is strictly positive and bounded (0 < φ(w) < ∞),

(38)${\displaystyle 0<\lambda \leqq d∕2\leqq \mu }$

holds, where d is the dimension of the parameter. If the likelihood and the posterior distribution can be approximated by a normal distribution, the equal signs hold.In this study, we use hierarchical models. For them, the likelihood and posterior distribution cannot be approximated by any normal distribution. Therefore, the Bayesian inference can make the generalization loss smaller than that of ML estimation, i.e., it is appropriate for constructing a predictive model:

(39)${\displaystyle \mathbb{E}\left[G_n^{Bayes}\right]+o\left(1∕n\right)<\mathbb{E}\left[G_n^{ML}\right].}$

Third, we discuss the evaluation criterion used in our analysis. We consider the case (a) as appropriate, so we are to minimize the generalization loss. In our case, the model selection should be ‘effective’, thus neither using ML estimation nor the coefficient of determination R² is appropriate. However, unfortunately, the generalization loss cannot be computed since the true distribution q(x) is unknown.

When a model is set in which the likelihood can be approximated by a normal distribution, the Akaike information criterion (AIC) can estimate the generalization loss $G_n^{ML}$ with the theoretical proof (Akaike, 1974). Moreover, minimizing AIC is the ‘effective method’ if the above assumption is satisfied. Hence, let A_n be AIC

(40)${\displaystyle A_n=-\frac{1}{n}\sum _{i=1}^{n}logp\left(x_i|\hat{w}\right)+\frac{d}{n},}$

we have

(41)${\displaystyle \mathbb{E}\left[G_n^{ML}\right]=\mathbb{E}\left[A_n\right]+o\left(\frac{1}{n}\right).}$

Note that we write the AIC in the scale of the generalization loss, not the deviance scale (2nA_n).

As mentioned above, our model set includes a hierarchical model, thus it is not appropriate to minimize AIC in the ML estimation. We chose inference in order to decrease the generalization loss $G_n^{Bayes}$. In addition, there are other advantages of the inference. The widely applicable information criterion (WAIC) Watanabe (2010a) can estimate $G_n^{Bayes}$ with mathematical proof even if the likelihood and the posterior cannot be approximated by any normal distribution. Moreover, minimizing WAIC is “effective”. Let W_n be WAIC

(42)${\displaystyle W_n=-\frac{1}{n}\sum _{i=1}^{n}log{\mathbb{E}}_w\left[p\left(x_i|w\right)\right]+\frac{1}{n}{\mathbb{V}}_w\left[logp\left(x_i|w\right)\right],}$

we have

(43)${\displaystyle \mathbb{E}\left[G_n^{Bayes}\right]=\mathbb{E}\left[W_n\right]+o\left(\frac{1}{n^2}\right),}$

without the assumption of normality for the posterior distribution (Watanabe, 2010a).

Another model evaluation criterion in Bayesian inference is well-known: widely applicable Bayesian information criterion (WBIC) Watanabe (2013). Let $W_n^B$ be WBIC. WBIC approximates the free energy: $W_n^B=F_n+O_p\left(\sqrt{logn}\right)$ Watanabe (2013). It is useful for statistical modeling of case (b). However, in our case, we do not suggest that the mechanistic model is included in the considered model set; instead, we perform for the case (a). Therefore, we concluded that the model evaluation using WAIC was more appropriate than using WBIC. Indeed, WAIC is useful for real data even if we are limited to behavioral data (Wakita, Kagaya & Aonuma, 2020; Barrett, McElreath & Perry, 2017).

Lastly, we discuss the difference between conditioned-level WAIC and marginal-level WAIC. Although WAIC is beginning to be used for evaluating models with empirical data, we should be careful to compute the value of a hierarchical model. Watanabe (2018) introduces two different definitions of WAIC depending on two different predictions. The often-used definition of WAIC for a hierarchical model is the first case in the book:

(44)${\displaystyle p_{posterior}\left(w_k|{\left(x_k\right)}^l\right)\propto \phi \left(w_k|w_0\right)\prod _{l=1}^{N_{act}}{p}_{model}\left(x_{lk}|w_k\right)}$ (45)${\displaystyle WAI{C}_k=-\frac{1}{N_{act}}\sum _{l=1}^{N_{act}}log{\mathbb{E}}_{w_k}\left[p_{model}\left(x_{lk}|w_k\right)\right]+\frac{1}{N_{act}}\sum _{l=1}^{N_{act}}{\mathbb{V}}_{w_k}\left[log{p}_{model}\left(x_{lk}|w_k\right)\right]}$ (46)${\displaystyle WAI{C}_{conditioned}=\sum _{k=1}^{N_{animal}}WAI{C}_k}$ (47)${\displaystyle k=1,\dots ,N_{animal};l=1,\dots ,N_{act}}$

where ${\left(x_k\right)}^l=\left(x_k^1,\dots ,x_k^{N_{act}}\right)$ represents all given data for an individual. Note that the number of behavioral acts of the k-th animal is the same (balanced) for simplicity here (The number is unbalanced in our data). It should be noted that the statistical model p_model(x_lk|w_k) is conditioned upon the w_k assigned to each individual. In other words, in prediction, this model focuses on new acts of the already obtained individuals, whereas our focus is usually on a new act of another individual in order to compare models. In many cases, when studying ‘animal personality’, we are not usually interested in how our models specifically explain the sampled animals personalities. Instead, we are interested in the distribution of a focused species. This is the reason why we did not use this conditioned-level WAIC. However, it is generally more convenient to use conditional-level likelihood within most Bayesian software, with the unfortunate consequence that conditional-level WAIC is often used. For example, Mitchell et al. (2016) uses conditioned-level WAIC to compare a hierarchical model with a non-hierarchical model. Furthermore, the use of the measure R² for the evaluation of the model in terms of ‘variance explained’ is based on how we can minimize the variation in the sample obtained, and not focused on the prediction of the true distribution (Nakagawa & Schielzeth, 2010).

In summary, we took ‘a’ best approach from the predictive point of view and explored ‘a’ best model rather than a ‘correct’ model, because it is a natural assumption that the distribution we build would never be correct in the empirical modeling of a natural behavior of the sponge crab.

Making cost and size choice: why did the crabs skip carrying the sponge?

The crabs in our experiments that did not carry caps were larger than those that carried caps. One possibility for the reason would be that when they grow up to some extent, the predators might avoid the crab and the relative energetic cost to make caps might increase. We speculate that this might be a reason why the large crabs did not make nor carry the caps.

Another possibility is that the sponges used in this experiment were smaller than those of the necessary size for the crabs. Dembowska (1926) reported that the proportion of caps that fit the size of D. personata tended to decrease with the increase of the size of the crabs, and we considered that the decision to skip carrying the caps was because there were few sponges that fit the large crabs. Similarly, the large crabs that skipped cap making and carrying, would carry a cap if a sponge size would be larger than the L size sponge. In addition, no individuals carried the S sponge, because it was too small for all of the crabs to carry. It is likely that a younger and smaller crab than those used in this experiment would carry the S sponge. However, we can not exclude the possibility that the variations of the morphology of dactyls or the molt stage affect the behavior.