When using a model-based approach to geostatistical problems, often, due to the complexity of the models, inference relies on Markov chain Monte Carlo methods. This article focuses on the generalized linear spatial models, and demonstrates that parameter estimation and model selection using Markov chain Monte Carlo maximum likelihood is a feasible and very useful technique. A dataset of radionuclide concentrations on Rongelap Island is used to illustrate the techniques. For this dataset we demonstrate that the log-link function is not a good choice, and that there exists additional nonspatial variation which cannot be attributed to the Poisson error distribution. We also show that the interpretation of this additional variation as either micro-scale variation or measurement error has a significant impact on predictions. The techniques presented in this article would also be useful for other types of geostatistical models.