<p>Performing Bayesian inference on large spatio-temporal models requires extracting inverse elements of large sparse precision matrices for marginal variances, as well as estimating model hyperparameters. Although direct matrix factorizations can be used for the inversion, such methods fail to scale well for distributed problems when run on large computing clusters. On the contrary, Krylov subspace methods for the selected inversion have been gaining traction. We propose a parallel hybrid approach based on domain decomposition, which extends the Rao-Blackwellized Monte Carlo estimator for distributed precision matrices. Our approach exploits the strength of Krylov subspace methods as global solvers and efficiency of direct factorizations as base case solvers to compute the marginal variances and the derivatives required for hyperparameter estimation using a divide-and-conquer strategy. By introducing subdomain overlaps, one can achieve greater accuracy at an increased computational effort with little to no additional communication. We demonstrate the speed improvements and efficient hyperparameter inference on both simulated models and a massive US daily temperature data.</p>