A fundamental pattern in ecology is that smaller organisms are more abundant than larger organisms. This pattern is known as the individual size distribution (ISD), which is the frequency of all individual sizes in an ecosystem, regardless of taxon. The ISD is described by power law distribution with the form f(x)=Cxλ, and a major goal of size spectra analyses is to estimate the ISD parameter λ. However, while numerous methods have been developed to do this, they have focused almost exclusively on estimating λ from single samples. Here, we develop an extension of the truncated Pareto distribution within the probabilistic modeling language Stan. We use it to estimate multiple ISD parameters simultaneously with a hierarchical modeling approach. The most important result is the ability to examine hypotheses related to size spectra, including the assessment of fixed and random effects, within a single Bayesian generalized (non)-linear mixed model.