General formulas are derived for the asymptotic bias in regression coefficients and variance components estimated by penalized quasi-likelihood (PQL) in generalized linear mixed models with canonical link function and multiple sets of independent random effects. Easily computed correction matrices result in variance component estimates that have satisfactory asymptotic behavior for small values of the variance components and significantly reduce bias for larger values. Both first-order and second-order correction procedures are developed for regression coefficients estimated by PQL. The methods are illustrated through an analysis of an experiment on salamander matings involving crossed male and female random effects, and their properties are evaluated in a simulation study.