Comparing two vectors from negative binomial distribution in R

Question

I'm using R and have two vectors of discrete values. They are not strictly speaking categorical because the values themselves are number of dots counted on the image of a cell (whole vector is all the cells on the image). There are two vectors: reference and a vector with dot counts after some perturbation

What I believe is that such data should follow negative binomial distribution and some sort of goodness of fit should give a p-value and some statistic describing whether the two distributions differ significantly.

What people advised me is that chi square test would do the trick but in my understanding chi square considers all values only as a category and ignores the fact that these are numbers and if lets say number of cells with 5 dots decreased a bit while number of cells with 4 dots increased that's not the same as if same situation would happen with 0-dots and 6-dots categories.

However what I didn't find is a test which could deal with negative binomial distributions. I hope I described the problem clearly. So if somebody know any test which would deal with such kind of data or if anybody thinks that my assumptions are wrong you are welcome to share your ideas.

Example 1

library(ggplot2)
c.dots = c(0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 1, 3, 0, 1, 0, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0,
           0, 1, 1, 0, 1, 1, 0, 1, 0, 0, 0, 2, 2, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 2,
           0, 0, 0, 1, 0, 0, 1, 0, 0, 3, 0, 1, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 1, 1,
           0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 2, 0, 0, 0, 1, 0,
           0, 0, 0, 0, 0, 0, 0, 2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, 0, 2, 0,
           0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 2, 0, 0, 0, 0, 0,
           0, 0, 1, 1, 0, 0, 1, 0, 2, 0, 1, 0, 2, 0, 0, 1, 0, 0, 1, 1, 0, 0, 3, 0, 0,
           0, 1, 0, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 0, 0, 0, 0, 0, 0, 0, 0)

w.dots = c(0, 0, 0, 1, 3, 1, 1, 1, 1, 0, 0, 2, 0, 0, 2, 1, 0, 1, 3, 0, 1, 0, 0, 0, 2,
           0, 2, 2, 0, 3, 1, 2, 1, 0, 2, 1, 0, 2, 0, 1, 2, 1, 0, 0, 1, 0, 1, 1, 0, 0,
           0, 1, 1, 0, 2, 0, 0, 1, 3, 0, 0, 1, 0, 2, 1, 0, 1, 1, 1, 1, 1, 1, 2, 1, 1,
           2, 4, 1, 0, 0, 2, 2, 0, 1, 0, 1, 3, 0, 2, 1, 1, 2, 0, 0, 0, 0, 0, 0, 1, 0,
           1, 1, 0, 1, 0, 0, 2, 0, 1, 0, 2, 1, 0, 1, 2, 0, 4, 2, 0, 1, 0, 2, 0, 1, 2,
           1, 1, 2, 1, 1, 3, 1, 0, 1, 0, 1, 2, 0, 1, 2, 0, 1, 1, 2, 2, 0, 3, 0, 1, 1,
           0, 0, 2, 0, 1, 1, 0, 1, 2, 0, 0, 1, 0, 1, 2, 0, 0, 4, 3, 0, 1, 0, 0, 1, 0,
           4, 0, 1, 0, 0, 0, 0, 0, 2, 1, 0, 0, 0, 0, 1, 1, 0, 3, 1, 1, 0, 4, 1, 1, 3)

chisq.test(rbind(table(w.dots), table(c.dots)))

nbrand = rnbinom(length(c.dots), mu = 1, size = 1)
ggplot() + 
  geom_density(aes(x=x), data=data.frame(x=c.dots), fill="red", alpha=0.5) + 
  geom_density(aes(x=x), data=data.frame(x=w.dots), fill="blue", alpha=0.5) +
  geom_density(aes(x=x), data=data.frame(x=nbrand), colour="green", alpha=0, linetype=3)

Example 2

library(ggplot2)
c.dots = c(1, 0, 0, 1, 0, 0, 3, 0, 1, 0, 3, 0, 2, 0, 0, 2, 2, 0, 0, 1, 1, 0, 0, 1, 0,
           0, 0, 0, 1, 0, 1, 0, 2, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0,
           1, 1, 1, 2, 0, 4, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 1, 3, 4,         
           0, 1, 1, 0, 1, 0, 2, 1, 2, 2, 3, 1, 0, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 0, 1,
           1, 1, 1, 1, 1, 1, 0, 1, 1, 0, 2, 0, 0, 3, 2, 2, 1, 0, 2, 0, 2, 2, 0, 0, 2,
           1, 0, 2, 0, 0, 2, 2, 1, 0, 0, 0, 0, 0, 1, 0, 3, 0, 1, 0, 0, 1, 0, 0, 0, 0,
           2, 1, 1, 0, 1, 0, 1, 1, 0, 1, 3, 0, 0, 1, 1, 1, 0, 1, 1, 0, 0, 1, 2, 1, 0,
           1, 2, 0, 0, 3, 3, 0, 1, 2, 0, 0, 1, 1, 0, 1, 1, 3, 1, 3, 0, 2, 0, 0, 0, 0)

w.dots = c(1, 3, 1, 1, 0, 1, 0, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 3, 0, 0, 0, 1, 2, 0,
           1, 0, 0, 0, 1, 0, 0, 1, 1, 1, 1, 0, 5, 1, 1, 0, 1, 0, 1, 0, 0, 1, 1, 1, 2,
           0, 1, 0, 3, 0, 0, 1, 2, 3, 1, 0, 0, 0, 2, 1, 1, 2, 0, 2, 0, 3, 0, 2, 0, 0,
           0, 0, 2, 0, 1, 0, 2, 0, 0, 1, 1, 2, 3, 0, 2, 2, 1, 0, 1, 0, 0, 1, 0, 1, 0,
           0, 1, 0, 0, 1, 0, 0, 1, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 2, 0,
           0, 1, 2, 1, 1, 1, 2, 1, 2, 3, 2, 0, 0, 0, 0, 0, 1, 2, 0, 0, 1, 0, 0, 1, 1,
           0, 0, 1, 2, 1, 0, 1, 2, 1, 1, 1, 0, 1, 0, 0, 0, 2, 1, 1, 1, 0, 0, 0, 0, 0,
           1, 2, 1, 2, 0, 1, 2, 1, 0, 1, 3, 2, 1, 0, 0, 0, 0, 0, 2, 1, 1, 2, 2, 1, 2)

chisq.test(rbind(table(w.dots), table(c.dots)))

nbrand = rnbinom(length(c.dots), mu = 1, size = 1)
ggplot() + 
  geom_density(aes(x=x), data=data.frame(x=c.dots), fill="red", alpha=0.5) + 
  geom_density(aes(x=x), data=data.frame(x=w.dots), fill="blue", alpha=0.5) +
  geom_density(aes(x=x), data=data.frame(x=nbrand), colour="green", alpha=0, linetype=3)

Answer 1

Your dependent variable is a count ("number of dots counted on the image of a cell"). Asking whether the distribution of counts is similar in two groups is conceptually the same as asking whether group membership matters for the distribution of counts.

I suggest a Poisson regression as a first step where you model the dot count with group membership. In a second step, one might then try to decide whether the Poisson assumption of "conditional variance = conditional mean" is violated, suggesting a move to a quasi-Poisson model, to a Poisson-model with heteroscedasticity-consistent (HC) standard error estimates, or to a negative binomial model.

Given data c.dots and w.dots as in the OP's example 1: We first create a data frame with predicted variable Y= number of dots and predictor X= factor with group membership. Then we run a standard Poisson regression

> dotsDf <- data.frame(Y=c(c.dots, w.dots),
+                      X=factor(rep(c("c", "w"), c(length(c.dots), length(w.dots)))))

> glmFitP <- glm(Y ~ X, family=poisson(link="log"), data=dotsDf)
> summary(glmFitP)                      # Poisson model
Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -0.9289     0.1125  -8.256  < 2e-16 ***
Xw            0.8455     0.1345   6.286 3.26e-10 ***

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 509.01  on 399  degrees of freedom
Residual deviance: 465.90  on 398  degrees of freedom
AIC: 862.16

This indicates a significant predictor "group membership = w" resulting from dummy coding the grouping factor (2 groups => 1 dummy predictor, c is the reference level). For comparison, we can run the quasi-Poisson model that has an extra dispersion parameter for the conditional variance.

> glmFitQP <- glm(Y ~ X, family=quasipoisson(link="log"), data=dotsDf)
> summary(glmFitQP)                     # quasi-Poisson model
Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -0.9289     0.1298  -7.158 3.96e-12 ***
Xw            0.8455     0.1551   5.450 8.85e-08 ***

(Dispersion parameter for quasipoisson family taken to be 1.330164)

    Null deviance: 509.01  on 399  degrees of freedom
Residual deviance: 465.90  on 398  degrees of freedom
AIC: NA

The parameter estimates are the same, but the standard errors of these estimates is slightly larger. The estimated dispersion parameter is slightly larger than 1 (the value in the Poisson-model), indicating some overdispersion. An alternative approach is to use a Poisson model with HC-consistent standard errors:

> library(sandwich)                     # for vcovHC()
> library(lmtest)                       # for coeftest()
> hcSE <- vcovHC(glmFitP, type="HC0")   # HC-consistent standard errors
> coeftest(glmFitP, vcov=hcSE)
z test of coefficients:

            Estimate Std. Error z value  Pr(>|z|)    
(Intercept) -0.92887    0.14084 -6.5952 4.246e-11 ***
Xw           0.84549    0.16033  5.2735 1.339e-07 ***

Again, somewhat larger standard errors. Now the negative binomial model:

> library(MASS)                         # for glm.nb()
> glmFitNB <- glm.nb(Y ~ X, data=dotsDf)
> summary(glmFitNB)                     # negative binomial model
Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -0.9289     0.1192  -7.791 6.66e-15 ***
Xw            0.8455     0.1456   5.806 6.40e-09 ***

(Dispersion parameter for Negative Binomial(3.212) family taken to be 1)

    Null deviance: 427.19  on 399  degrees of freedom
Residual deviance: 391.20  on 398  degrees of freedom
AIC: 856.87

              Theta:  3.21 
          Std. Err.:  1.46 

 2 x log-likelihood:  -850.867 

You can test the negative binomial model against the Poisson model in a likelihood ratio test for the model comparison:

> library(pscl)                         # for odTest()
> odTest(glmFitNB)
Likelihood ratio test of H0: Poisson, as restricted NB model:
n.b., the distribution of the test-statistic under H0 is non-standard

Critical value of test statistic at the alpha= 0.05 level: 2.7055 
Chi-Square Test Statistic =  7.2978 p-value = 0.003452

The result here indicates that the data are unlikely to come from a Poisson model.

For the OP's example 2, all these tests are non-significant.

Note that I slightly shortened the output from glm() and glm.nb().

Answer 2

First of all, +1 for @caracal's answer.

One other possibility would be as follows: fit a negbin distribution to both datasets (e.g., using fitdistr in the MASS package). The parameter estimates are (asymptotically) normally distributed, and we can get the covariance estimates this way:

vcov(fitdistr(w.dots,densfun="negative binomial"))
vcov(fitdistr(c.dots,densfun="negative binomial"))

Unfortunately, I couldn't find a good way to test the null hypothesis "both vectors are generated by the same negbin parameter vector", which should be a multivariate generalization of the standard one-dimensional two-sample t-test with pooled variances. This may help, but I don't have access. Update: I did think of something, see below.

However, what we can do is bootstrap the fits and plot the bootstrapped fits:

library(boot)
library(MASS)
n.boot <- 1000
w.fit.boot <- boot(w.dots,R=n.boot,
  statistic=function(xx,index)fitdistr(xx[index],densfun="negative binomial")$estimate)
c.fit.boot <- boot(c.dots,R=n.boot,
  statistic=function(xx,index)fitdistr(xx[index],densfun="negative binomial")$estimate)

plot(c.fit.boot$t,pch=21,bg="black",xlab="size",ylab="mu",log="xy",cex=0.5,
  xlim=range(rbind(w.fit.boot$t,c.fit.boot$t)[,1]),
  ylim=range(rbind(w.fit.boot$t,c.fit.boot$t)[,2]))
abline(v=c.fit.boot$t0[1]); abline(h=c.fit.boot$t0[2])
points(w.fit.boot$t,pch=21,bg="red",col="red",cex=0.5)
abline(v=w.fit.boot$t0[1],col="red"); abline(h=w.fit.boot$t0[2],col="red")
legend(x="bottomright",inset=.01,pch=21,col=c("black","red"),pt.bg=c("black","red"),
  legend=c("c.dots","w.dots"))

Each dot corresponds to one bootstrapped fitted negbin parameter vector. The horizontal and vertical lines give the fitted parameters for the original data. Results:

It looks to me like the data in example 1 pretty obviously do not come from the same distribution, while those in example 2 might well do so.

Of course, all this is subject to the negbin distribution being correct. However, one could do a similar exercise with a Poisson model or other ways of accounting for overdispersion.

Update: now, we want to test the amount of overlap between the joint distributions of the negbin parameter estimates based on the two original vectors. Let's start with a simple one-dimensional example. Assume we have estimated means and standard deviations of two normal densities. Under the null hypothesis, these two densities are identical. So one possible test statistic is the amount of overlap under the density curves:

means <- c(1,2)
sds <- c(.2,.3)
xx <- seq(means[1]-3*sds[1],means[2]+3*sds[2],by=.001)
plot(xx,dnorm(xx,means[1],sds[1]),type="l",xlab="",ylab="")
lines(xx,dnorm(xx,means[2],sds[2]))
polygon(x=c(xx,rev(xx)),y=c(rep(0,length(xx)),
  rev(pmin(dnorm(xx,means[1],sds[1]),dnorm(xx,means[2],sds[2])))),col="grey")

We therefore calculate this area using simple numerical integration:

sum(pmin(dnorm(xx,means[1],sds[1]),dnorm(xx,means[2],sds[2])))*mean(diff(xx))
[1] 0.04443207

We can easily generalize this approach to two dimensions. In this case, we have to consider the estimated means and their covariance ellipses and overlay a two-dimensional grid for integration. I'll use the bootstrapped values above to determine the grid dimensions.

library(mvtnorm)
nn <- 100
xx <- seq(min(c(c.fit.boot$t[,1],w.fit.boot$t[,1])),
  max(c(c.fit.boot$t[,1],w.fit.boot$t[,1])),length.out=nn)
yy <- seq(min(c(c.fit.boot$t[,2],w.fit.boot$t[,2])),
  max(c(c.fit.boot$t[,2],w.fit.boot$t[,2])),length.out=nn)
integral <- matrix(NA,nrow=nn,ncol=nn)
for ( ii in 1:nn ) {
  for ( jj in 1:nn ) {
    integral[ii,jj] <-
      min(dmvnorm(c(xx[ii],yy[jj]),mean=c.fit$estimate,sigma=vcov(c.fit)),
        dmvnorm(c(xx[ii],yy[jj]),mean=w.fit$estimate,sigma=vcov(w.fit)))
  }
}
sum(integral)*mean(diff(xx))*mean(diff(yy))
[1] 6.673166e-05  # for the data in example 1

Of course, both the one- and the two-dimensional integration can be carried out with pen & paper and/or a computer algrebra system. integrate() in R works for the one-dimensional case, I imagine that there are tools in R to integrate two-dimensional functions. In addition, I assume I will get a few downvotes for the double loop above, but I couldn't get a vectorized function (outer() or mapply()) to work - any comments would be welcome. Finally, an equally spaced grid is probably not the most accurate way to do this kind of integration.

A fully vectorized implementation is straight forward, 'integral' doesn't have to be a 2D matrix:

library(mvtnorm)
nn <- 100
xx <- seq(min(c(c.fit.boot$t[,1],w.fit.boot$t[,1])),
  max(c(c.fit.boot$t[,1],w.fit.boot$t[,1])),length.out=nn)
yy <- seq(min(c(c.fit.boot$t[,2],w.fit.boot$t[,2])),
  max(c(c.fit.boot$t[,2],w.fit.boot$t[,2])),length.out=nn)
#vectorized:
xy_pair <- cbind(rep(xx, each=nn, times=1), rep(yy, each=1, times=nn))
integral <- min(dmvnorm(xy_pair,mean=c.fit$estimate,sigma=vcov(c.fit)),
        dmvnorm(xy_pair,mean=w.fit$estimate,sigma=vcov(w.fit)))
sum(integral)*mean(diff(xx))*mean(diff(yy))

Answer 3

The present answer assumes that both distributions are actually Negative Binomial (NB), and also assumes the two samples to be independent.

In this framework, you have a model with four parameters, two for each distribution, and want to test a restriction with two constraints on the parameters. This can be done by standard Maximum Likelihood (ML) theory. The Likelihood Ratio test (LR) is straightforward: the fitdistr function from MASS used above by Stephan Kolassa can provide the maximised likelihoods needed in the LR. Another test is known as score test or Lagrange Multiplier test (LM); it has nice theoretical properties but needs an information matrix and a score vector from therestricted ML estimation. A third possibility would be a Wald test.

Let $ \boldsymbol{\psi}_{\textrm{c}} := [r_{\textrm{c}},\,\pi_{\textrm{c}}]'$ be the vector of parameters for the "c" sample, with similar notation for the "w". Here $r$ stands for size, while $\pi$ stands for prob in R i.e., the probability of success. The full parameter vector is $\boldsymbol{\theta} := [\boldsymbol{\psi}_{\textrm{c}}',\, \boldsymbol{\psi}_{\textrm{w}}']'$. You want to test the hypothesis $H_0:\,\boldsymbol{\psi}_{\textrm{c}} = \boldsymbol{\psi}_{\textrm{w}}$. The LM test is based on the approximate distribution $$ \mathbf{U}'\mathbf{I}^{-1}\mathbf{U} \sim \chi^2(d) $$ where the degree of freedom $d$ is the number of scalar restrictions, here $2$. The vector $\mathbf{U}$ and the matrix $\mathbf{I}$ are the score vector and the information matrix at the restricted estimate, say $\widehat{\boldsymbol{\theta}}_{0}= [{\widehat{\boldsymbol{\psi}}_{0}}',\, {\widehat{\boldsymbol{\psi}}_{0}}']'$.

Here $\mathbf{U}$ is a vector of length $4$ and $\mathbf{I}$ is $4 \times 4$. The LR statistic can here be obtained by noticing that $\mathbf{I}$ is block diagonal with two $2\times 2$ blocks leading to two contributions arising from the two samples. The score and observed information matrix can be used in closed form. With the fitdistr function, the parametrisation differs from above and uses $\mu:= r\times (1-\pi)/\pi$. Rather than a parameter change, we can with a little more effort use a concentrated log-likelihood.

The two tests here have very small $p$-values hence one should clearly reject the null hypothesis.

Using simulations, it seems that the information matrix $\mathbf{I}$ can in some rare cases be ill-conditioned, and the LM test statistic can then unduly become negative. In this cases at least, the LR test must be preferred to the LM test.

## fit NB distr from a sample X using concentrated logLik
fitNB <- function(X) {
  n <- length(X)
  loglik.conc <- function(r) {
    prob <- n*r / (sum(X) + n*r)
    sum( lgamma(r + X) - lgamma(r) - lgamma(X + 1) +
        r * log(prob) + X * log(1 - prob) ) 
  }
  ## find 'r' with an 1D optim...
  res <- optimize(f = loglik.conc, interval = c(0.001, 1000),
                  maximum = TRUE)
  r <- res$maximum[1]
      params <- c(size = r, prob = n*r / (sum(X) + n*r))
      attr(params, "logLik") <- res$objective[1]
  params
}
## compute score vector and info matrix at params 'psi' using closed forms
scoreAndInfo <- function(psi, X) {
  size <- psi[1]; prob <- psi[2]
  n <- length(X)
  U <- c(sum(digamma(size + X) - digamma(size) + log(prob)),  
         sum(size / prob - X / (1-prob) ))
  I <- matrix(c(- sum(trigamma(size + X) - trigamma(size)),  
                -n / prob, -n / prob,  
                sum( size / prob^2  + X / (1-prob)^2)),
              nrow = 2, ncol = 2)
  names(U) <- rownames(I) <- colnames(I) <- c("size", "prob")
  LM <-  as.numeric(t(U) %*% solve(I) %*% U)
  list(score = U, info = I, LM = LM)
}
## continuing on the question code a is for "all" 
c.fit <- fitNB(X = c.dots)
w.fit <- fitNB(X = w.dots)
a.fit <- fitNB(X = c(c.dots, w.dots))
## LR test and p.value
D.LR <- 2 * ( attr(c.fit, "logLik") + attr(w.fit, "logLik") ) -
  2 * attr(a.fit, "logLik") 
p.LR <- pchisq(D.LR, df = 2, lower.tail = FALSE)
## use restricted parameter estimate to compute the LM contributions
c.sI <- scoreAndInfo(psi = a.fit, X = c.dots) 
w.sI <- scoreAndInfo(psi = a.fit, X = w.dots) 
D.LM <- c.sI$LM + w.sI$LM 
p.LM <- pchisq(D.LM, df = 2, lower.tail = FALSE)

Answer 4

You can use a Bayesian technique. http://www.indiana.edu/~kruschke/BEST/ This gives you the posterior distribution of the sample means for comparison.

x<-c(rep(1,200),rep(2,200))
y<-c(c.dots,w.dots)

require(R2jags)  #open jags console.
dataList<-list(x=x,y=y, Ntotal=length(y))

modelstring = "
model {
for ( i in 1:Ntotal ) {
y[i] ~ dnegbin( p[x[i]] , r[x[i]] )
}
for ( j in 1:2) {
p[j] <- r[j]/(r[j]+m[j])
m[j] ~ dgamma(0.01, 0.01)
r[j] ~ dgamma(0.01, 0.01)
v[j] <- r[j]*(1-p[j])/(p[j]*p[j])
}
}
"
writeLines(modelstring,con="model.txt")

parameters=c("m")
adaptSteps = 1000              
burnInSteps = 1000               
nChains = 1                   
numSavedSteps=10000           
thinSteps=1
nPerChain = ceiling( ( numSavedSteps * thinSteps ) / nChains ) 

JagsModel = jags.model( "model.txt" , data=dataList  , 
                        n.chains=nChains , n.adapt=adaptSteps )

codaSamples = coda.samples( JagsModel , variable.names=parameters ,
                            n.iter=nPerChain , thin=thinSteps )

m <- as.matrix(codaSamples)
head(m)
plot(density(m[,1]))
lines(density(m[,2]))