We describe techniques to fit linear and generalized linear mixed models with different types of response variables and their respective parameter estimation. We consider the traditional frequentist and Bayesian approaches to model development and parameter estimation. For both paradigms we develop Markov chain Monte Carlo (MCMC) algorithms to fit the models. We first show the development of the linchpin variable sampler and its application to Bayesian hierarchical linear mixed models. We start by looking at a joint density that can be factored into the product of a conditional density from which is easy to generate observations and a more complex marginal density. We propose an algorithm that samples sequentially from the marginal density using an MCMC step and directly from the conditional density. We show theoretical properties that are consequences of this sampling technique and empirical properties by studying computational results. Secondly, we turn our attention to generalized linear mixed models (GLMMs). In particular we work with Bernoulli, Poisson, negative binomial, and gamma responses with normal and t distributed random effects. For each of these models we implement a Monte Carlo Expectation-Maximization (MCEM) algorithm in the package mcemGLM. In addition to parameter and standard error estimation, the package performs analysis of variance tests, linear combination estimation, and other functionality.