## Appendix A Markov Chain Monte Carlo Algorithm  
In this appendix, we provide a Gibbs sampling algorithm for calculating the posterior distributions of the parameters of interest, \(\mathbf{\theta}=(\mathbf{\beta},\sigma^{2},\lambda_{f},\mathbf{\gamma},\omega)\). The algorithm is modified based on the Markov chain Monte Carlo (MCMC) sampler proposed in [PERSON] et al. (2015) [5]and [PERSON] et al. (2018) [6]. To give the detailed description of the algorithm, we first re-express the spike-and-slab prior in (4) hierarchically using binary latent variables as  
\[\gamma_{m}\mid r_{m},a_{\gamma},b_{\gamma} \sim r_{m}\partial_{0}(\gamma_{m})+(1-r_{m})_{\mathcal{B}}(\gamma_{m} \mid a_{\gamma},b_{\gamma}),\] \[r_{m}\mid\omega \sim Bernoulli(\omega),\] \[\omega \sim Beta(a_{\omega},b_{\omega}),\;a_{\omega}=b_{\omega}=1/2.\] (A1)  
Then, after integrating over \(\omega\) and \(\mathbf{f}\), the joint posterior distribution is proportional to  
\[\pi(\mathbf{\theta}\mid\bar{y}_{1},\ldots,\bar{y}_{T}) \propto \mathcal{N}(\bar{\mathbf{y}}\mid\mathbf{Z}\beta,\sigma^{2}\mathbf{\Sigma}) \Bigg{\{}\prod_{m=1}^{M}[r_{m}\delta_{0}(r_{m})+(1-r_{m})g(\gamma_{m}\mid a_{ \gamma},b_{\gamma})]\Bigg{\}}\] \[\times\Gamma(\sum_{m}r_{m}+a_{\omega})\Gamma(M-\sum_{m}r_{m}+b_{ \omega})\] \[\times\mathcal{Gamma}(\sigma^{-2}\mid a_{\sigma},b_{\sigma}) \mathcal{Gamma}(\lambda_{f}\mid a_{\lambda},b_{\lambda}),\]  
where \(\bar{\mathbf{y}}=(\bar{y}_{1},\ldots,\bar{y}_{T})^{\top}\), \(\mathbf{Z}=(\mathbf{z}_{1}^{\top},\ldots,\mathbf{z}_{T}^{\top})^{\top}\) is the design matrix for linear effects, and \(\mathbf{\Sigma}=\text{diag}(1/n_{1},\ldots,1/n_{T})+\lambda_{f}K\). The algorithm proceeds proceeds iteratively by sampling from the full conditional posterior distributions of each parameter. Below, we outline the main steps of the Gibbs sampling procedure:  
1. Set initial values for all parameters in \(\mathbf{\theta}=(\mathbf{\beta},\sigma^{2},\lambda_{f},\gamma,\omega)\).
2. For each iteration, sample from the following full conditional posterior distributions: 1. Step 1: Sample the coefficients for linear fixed effects, \(\mathbf{\beta}\), from the full conditional distribution \[\mathbf{\beta}\mid\sigma^{2},\lambda_{f},\mathbf{r},\bar{\mathbf{y}}\sim\mathcal{N}(\mathbf{ \beta}\mid\mathbf{\Sigma}_{n}\mathbf{Z}^{\top}\mathbf{\Sigma}^{-1}\bar{\mathbf{y}},\sigma^{2} \mathbf{\Sigma}_{n}),\] where \(\mathbf{\Sigma}_{n}=(\mathbf{Z}^{\top}\mathbf{\Sigma}^{-1}\mathbf{Z})^{-1}\). 2. Step 2: Sample the variance, \(\sigma^{2}\), from the full conditional distribution given by \[\sigma^{-2}\mid\mathbf{\beta},\lambda_{f},\mathbf{r},\bar{\mathbf{y}}\sim\mathcal{Gamma} (\sigma^{-2}\mid a_{\sigma}+n/2,b_{\sigma}+SSE/2),\] where \(SSE=(\bar{\mathbf{y}}-\mathbf{Z}\mathbf{\beta})^{\top}\mathbf{\Sigma}_{n}^{-1}(\bar{\mathbf{y}}- \mathbf{Z}\mathbf{\beta})\). 3. Step 3: Sample \(\lambda_{f}\) using a Metropolis-Hastings method because the full conditional distribution does not the closed form and is proportional to \[\pi(\lambda_{f}\mid\mathbf{\beta},\sigma^{2},\mathbf{r},\bar{\mathbf{y}})\propto|\mathbf{ \Sigma}|^{-1/2}\exp\Big{\{}-SSE/(2\sigma^{2})\Big{\}}\lambda_{f}^{a_{\lambda} -1}\exp\{-b_{\lambda}/\sigma^{2}\}.\] Specifically, we generate a candidate sample from a gamma distribution, where the mean is set to the value of \(\lambda_{f}\) from the previous iteration, and the variance is adjusted to achieve an appropriate acceptance rate. 4. Step 4: Sample \((\mathbf{\gamma},\mathbf{r})\) using an adaptive Metropolis-Hastings algorithm jointly from the following distribution \[\pi(\mathbf{\gamma},\mathbf{r}\mid\mathbf{\beta},\sigma^{2},\lambda_{f},\bar{\mathbf{y}}) \propto \Gamma(\sum_{m}r_{m}+a_{\omega})\Gamma(M-\sum_{m}r_{m}+b_{\omega})\] \[\Bigg{\{}\prod_{m=1}^{M}[r_{m}\delta_{0}(r_{m})+(1-r_{m})g( \gamma_{m}\mid a_{\gamma},b_{\gamma})]\Bigg{\}}.\]  
For more details, see [PERSON] et al. (2015) [5].