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?
 A: 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)
A: 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.

A: 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)

