# How to estimate the parameters of data with greater tail (seems as negative binomial)

I have a data with only 70 observations. The data contain one discrete variable. This variable contains the number of patient with flu (count). So, we have 70 days and each day we count how many number of patient come to the clinic with flu. Then, I would like to model the relationship between the number of patients and their age, health and other variables. However, I must understand the shape for the count variable.

The mean of my data is 300 and the variance is 1826. Hence, I think, the distribution is negative binomial. I would like to have Density, distribution function, quantile function and random generation for the negative binomial distribution with parameters size and prob. in R using the following functions.

dnbinom(x, size, prob, mu, log = FALSE)
pnbinom(q, size, prob, mu, lower.tail = TRUE, log.p = FALSE)
qnbinom(p, size, prob, mu, lower.tail = TRUE, log.p = FALSE)
rnbinom(n, size, prob, mu)


I know that mu is the mean of my variable, which is easy to calculate. However, how do I know the size and prob?

Is the size can be taken as the number of the patient? for me, I think I cannot!

Any idea or help to figure this point?

• There are several versions (parameterizations) of the negative binomial distribution. Some count the number of failures before a particular number of successes and some count the trial on which the particular number of successes occurs. Please state the PDF of your version, say what the parameters mean, and show the formulas for the mean and variance.in terms of the parameters. Jul 2 '20 at 8:10

The other answers would be fine if you were simply estimating the marginal distribution, i.e. the overall distribution of the number of visits per day ignoring any other characteristics of the patient population. However, it sounds from the comments as though you want to do negative binomial regression, i.e. a model where you estimate the relationships between various covariates and the outcome. In this case, you want the conditional distribution to be negative binomial, and you can't estimate the distribution before you start - the model will do this as part of the fitting process.

Since you're using R, the standard tool for negative binomial regression would be the glm.nb() function in the MASS package (which is built-in; no need to install it). To fit the simplest possible model (i.e. intercept-only; the expected number of patients is the same independent of day or anything else you can measure) you would say

library(MASS)
model_fit <- glm.nb(count~1, data=your_data)


This would give you results that were actually very similar to fitdistr (which is also in the MASS package), because you are not using any covariates. If you wanted to look for trends over time, you could change the formula to count ~ day (an intercept term is still automatically included).

To evaluate whether there are problematic failures of model assumptions, I'd recommend the DHARMa package, starting with

## install.packages("DHARMa") ## if necessary
library(DHARMa)
vignette("DHARMa") ## read the documentation first!
plot(simulateResiduals(model_fit))


I foresee some difficulties/confusion with your goal of

model[ing] the relationship between the number of patients and their age, health and other variables

Your response variable (counts per day) is observed at the scale of the population. The "age, health and other variables" of the patients are measured at the individual level. I'm not sure how you're going to model this with a negative binomial regression. You could, for example

• test whether these characteristics vary over time (but then the response variables would be something like the mean age of the patients on a day, not a count)
• if you knew these characteristics for the general population, test whether their values in your patient population were significantly different from those in the general population (again, not a count response)
• Thank you so much your help is mazing. Could you please explain to me why do you use number here glm.nb(count~1,? Is that because we do not have any covariate here? Jul 3 '20 at 8:27
• yes. ~1 means an intercept-only model. Jul 3 '20 at 11:31

I think the negative binomial can describe the number of trials you need to attempt until you get a given number of successes. E.g., it can answer the question "how many people do I need to screen before I get 20 with flu?"

In your case you could find the parameters of the negative binomial using the R function fitdistr in package MASS:

library(MASS)

count <- rnbinom(70, size= 10, prob= 0.1) # Your data (you don't know size and prop)

params <- fitdistr(count, "negative binomial")
params
size           mu
11.048490   1009.657143
(   1.861756) (  36.503965)


Now you can use the estimated parameters for the *nbinom functions. Note that you need prob and size OR mu and size (not all the three parameters). If you want you can get the third parameter using the equation prob <- size/(size+mu) (from the R docs for rnbinom)

• Thank you so much for your help. Can I use fitdistr to select for me the best fit distribution? Jul 2 '20 at 8:30
• @Alice see if this helps stats.stackexchange.com/questions/132652/… Jul 2 '20 at 8:35
• I think fitdistribplus works only for continuous data (continuous distribution), my data is discrete. Jul 2 '20 at 8:45
• I have used your code with my data, and the estimated parameters are quite close (not exactly) to the real one. Jul 2 '20 at 9:44

As always, there are many different ways of estimating the parameters. Negative binomial is a distribution for the number of successes $$k$$ observed until total $$r$$ failures happen. Assuming $$r$$ is known, maximum likelihood estimator for the probability of success $$p$$ is

$$\hat p_\text{MLE} = \frac{r}{r+k}$$

it is biased, though there is minimum-variance, unbiased estimator

$$\hat p_\text{MVUE} = \frac{r-1}{r+k-1}$$

R enables you to use alternative parametrization, where instead of $$p$$, we can use $$\mu$$, defined as

An alternative parametrization (often used in ecology) is by the mean mu (see above), and size, the dispersion parameter, where prob = size/(size+mu). The variance is mu + mu^2/size in this parametrization.

• Thank you so much for your help. I just have a count of the number of patients each day. So, I really do not know what is the size. My data contain 70 observation (70 days). Can I use it as the size? That is size = 70 Jul 2 '20 at 8:16
• @Alice but what was the stopping criterion, i.e. "the number of failures $r$"?
– Tim
Jul 2 '20 at 8:19
• @Alice no, the number of days is different thing. Why do you assume that this data follows negative binomial distribution?
– Tim
Jul 2 '20 at 8:19
• I need to model this data using a regression model. The model I will use requires the distribution of the margins and requires (dnbinom pnbinom qnbinom rnbinom). Hence, I should find these functions before use the regression model. Jul 2 '20 at 8:24
• I assume the data is negative binomial because the variance is very large than the mean. Jul 2 '20 at 8:24