How to compare linear vs negative binomial vs random fit for data? I'm trying to plot number of offspring (x-axis) vs size of offspring (y-axis). I want to check if the size of offspring has a linear, or negative binomial or randomly fit. I'm trying to statistically compare these models. I have little knowledge on how exactly to go about it. One way I've read about is using maximum likelihood estimations and comparing AIC values. Are there any good resources (books/ tutorials) showing how to do multiple comparisons? Are there other ways of comparing?
 A: The first question you must ask yourself is the nature of the $y$ variable you want to consider. Can the size of offspring be modelled as a count data variable ? Does it have a distribution that would make sense for a count data process? (Count data models can be used to model non-discrete $y$ variable). 
More importantly, you will not make the same assumptions regarding the relationship between your $y$ and $X$s in a OLS regression and in a count data regression. Before turning out to negative binomial, I would recommend you to have a look at Poisson models. Poisson regression is a glm model. In Poisson regression, the link function gives the relationship between endogenous and exogenous variables writes as
$$
    \lambda_i = \mathbb{E}_{f, \theta}(Y_{i}|X_{i}) = \exp(X_{i} \beta)
$$
With a Poisson data generating process, $\mathbb{E}(Y_i |X_{i}) = \mathbb{V}(Y_i|X_{i}) = \lambda_i$. Thus, a Poisson model yields the restriction that conditional mean and variance are equal. Note that, by definition, heteroskedasticity is introduced in a Poisson model. This is not the case in standard OLS. Does this relationship between number and size of offsprings seem consistent ?  With R, you can use glm with family to fit this kind of model.
If you think the mean = variance assumption is unsuitable, you can turn to negative binomial models (or other types of model that enable excess-variance). In that case, $\mathbb{E}(Y_i |X_{i}) = \mu_i$ and $\mathbb{V}(Y_i |X_{i}) = \mu_i + \mu_i^2/\theta$. Once again, link function writes $\mu_i = \mathbb{E}_{f, \theta}(Y_{i}|X_{i}) = \exp(X_{i} \beta)$. This model is sometimes called Negative Binomial I. You have an additional parameters $\theta$ (or $\alpha = 1/\theta$ in some parametrizations). To fit these models, you can use MASS::glm.nb function. 
There are several ways to estimate these models. Rather than handcoding your maximum likelihood estimator, I would recommend you to use the techniques available in glm and MASS::glm.nb routines. The choice of the best model will depend (i) on the suitability of your distributional assumptions if there are some arguments in favor of an approach (ii) on some performance criterion. Many choices possible, AIC, log-likelihood, BIC...
