In recent years, interest in spatial statistics has increased significantly. However, for large data sets, statistical computations for spatial models are a challenge, as it is extremely difficult to store a large covariance or an inverse covariance matrix, and compute its inverse, determinant or Cholesky decomposition. This talk will focus on spatial mixed models and discuss scalable matrix-free conditional samplings for their inference. The role of shrinkage in the estimation will be considered. Both Bayesian computations and frequentist method of inference will be considered. The work arose in collaboration with Somak Dutta at Iowa State University.