  One of the most celebrated findings in complex systems in the last decade is
that different indexes y (e.g., patents) scale nonlinearly with the
population~x of the cities in which they appear, i.e., $y\sim x^\beta, \beta
\neq 1$. More recently, the generality of this finding has been questioned in
studies using new databases and different definitions of city boundaries. In
this paper we investigate the existence of nonlinear scaling using a
probabilistic framework in which fluctuations are accounted explicitly. In
particular, we show that this allows not only to (a) estimate $\beta$ and
confidence intervals, but also to (b) quantify the evidence in favor of $\beta
\neq 1$ and (c) test the hypothesis that the observations are compatible with
the nonlinear scaling. We employ this framework to compare $5$ different models
to $15$ different datasets and we find that the answers to points (a)-(c)
crucially depend on the fluctuations contained in the data, on how they are
modeled, and on the fact that the city sizes are heavy-tailed distributed.
