Comparando dos vectores de la distribución binomial negativa en R

Question

Estoy usando R y tienen dos vectores de valores discretos. Ellos no son estrictamente hablando categórica debido a que los valores propios son el número de puntos contados en la imagen de una celda (de todo el vector de todas las celdas de la imagen). Hay dos vectores: referencia y un vector con puntos de cuenta después de algunos perturbación

Lo que yo creo es que este tipo de datos debe seguir distribución binomial negativa y una especie de bondad de ajuste debe dar un valor de p y algunas estadísticas que describen si las dos distribuciones difieren significativamente.

Lo que la gente me aconsejó es que la prueba de la chi cuadrado que hacer el truco, pero en mi comprensión de la chi cuadrado considera todos los valores sólo como una categoría e ignora el hecho de que estos son los números y si digamos que el número de células con 5 puntos disminuido un poco, mientras que el número de células con 4 puntos mayor que no es el mismo como si a la misma situación ocurre con 0 puntos y 6 puntos de las categorías.

Sin embargo, lo que no encuentro es una prueba que podía hacer frente a la negativa de distribuciones binomiales. Espero que describe el problema con claridad. Así que si alguien sabe alguna prueba de que iba a tratar con este tipo de datos o si alguien piensa que mis suposiciones son incorrectas, son bienvenidos a compartir sus ideas.

Ejemplo 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)

Ejemplo 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

Su variable dependiente es un count ("número de puntos contados en la imagen de un celular"). Preguntar si la distribución de la cuenta es similar en los dos grupos es conceptualmente el mismo que preguntar si la pertenencia al grupo es importante para la distribución de cargos.

Sugiero una regresión de Poisson como un primer paso donde se modelo el punto de contar con el grupo de pertenencia. En un segundo paso, entonces se podría tratar de decidir si la distribución de Poisson supuesto de que la "varianza condicional = media condicional" es violado, lo que sugiere un movimiento a una cuasi-modelo de Poisson, a una Poisson con el modelo de heterocedasticidad-consistente (HC) error estándar de las estimaciones, o a un modelo binomial negativo.

Dado datos c.dots y w.dots como en el OP ejemplo 1: en primer lugar, crear un marco de datos con predichos de la variable Y= número de puntos y predictor X= factor de pertenencia al grupo. Luego hacemos un estándar de regresión de Poisson

> 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

Esto indica un predictor significativo miembro del grupo" = w" resultante de maniquí de la codificación de la agrupación de los factores (2 grupos => 1 dummy predictor, c es el nivel de referencia). Para la comparación, podemos ejecutar el cuasi-Poisson modelo que tiene una mayor dispersión de parámetros para la varianza condicional.

> 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

Las estimaciones de los parámetros son los mismos, pero los errores estándar de estas estimaciones es ligeramente más grande. La estimación de la dispersión de parámetros es ligeramente mayor que 1 (el valor de la distribución de Poisson-modelo), lo que indica una sobredispersión. Un enfoque alternativo es el uso de un modelo de Poisson con HC-consistente errores estándar:

> 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 ***

De nuevo, algo más grande de los errores estándar. Ahora el modelo binomial negativo:

> 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 

Usted puede probar el modelo binomial negativo contra el modelo de Poisson en una prueba de razón de verosimilitud para el modelo de comparación:

> 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

El resultado indica que los datos son poco probable que se de un modelo de Poisson.

Para el OP del ejemplo 2, todas estas pruebas no son significativas.

Tenga en cuenta que estoy ligeramente acortado la salida de glm() y glm.nb().

Answer 2

Primero de todo, +1 para @caracal la respuesta.

Otra posibilidad sería como sigue: ajuste de un negbin de distribución para ambos conjuntos de datos (por ejemplo, el uso de fitdistr en la MASS paquete). Las estimaciones de los parámetros son (asintóticamente) se distribuye normalmente, y podemos obtener la covarianza de las estimaciones de esta manera:

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

Por desgracia, no pude encontrar una buena manera de poner a prueba la hipótesis nula de "ambos vectores son generados por el mismo negbin vector de parámetros", que debe ser multivariante, la generalización de la norma uno-dimensional two-sample t-test con agruparon desviaciones. Esto puede ayudar, pero no tengo acceso. Actualización: me hizo pensar en algo, ver a continuación.

Sin embargo, lo que podemos hacer es arrancar el ajuste y la trama el bootstrap se ajusta:

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"))

Cada punto corresponde a una bootstrap equipada negbin vector de parámetros. Las líneas horizontales y verticales dar el conjunto de los parámetros para los datos originales. Resultados:

A mí me parece que los datos que en el ejemplo 1, bastante, obviamente, no provienen de la misma distribución, mientras que aquellos en el ejemplo 2 bien podría hacerlo.

Por supuesto, todo esto está sujeto a la negbin distribución de estar en lo correcto. Sin embargo, uno podría hacer un ejercicio similar con un modelo de Poisson u otras formas de contabilidad para la sobredispersión.

Actualización: ahora, queremos probar la cantidad de solapamiento entre las distribuciones conjuntas de la negbin estimaciones de los parámetros de la base de los dos originales de vectores. Vamos a empezar con un simple unidimensional ejemplo. Supongamos que tenemos estimación de medias y desviaciones estándar de las dos normales densidades. Bajo la hipótesis nula, estos dos densidades son idénticos. Así pues, una posible estadístico de prueba es la cantidad de superposición bajo la curva de densidad:

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")

Por lo tanto, calcular esta área mediante la simple integración numérica:

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

Fácilmente podemos generalizar este enfoque de dos dimensiones. En este caso, tenemos que considerar la estimación de los medios y su covarianza elipses y superposición de una cuadrícula bidimensional para la integración. Voy a utilizar el bootstrap de los valores anteriores para determinar las dimensiones de la red.

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

Por supuesto, tanto el uno y el dos dimensiones de la integración puede llevarse a cabo con lápiz y papel y/o un equipo con álgebra sistema. integrate() en R trabaja para el caso unidimensional, me imagino que hay herramientas en R para integrar dos dimensiones de funciones. Además, supongo que voy a tener un par de downvotes para el doble bucle de arriba, pero no pude conseguir un vectorizados función (outer() o mapply()) para el trabajo - cualquier comentario será bienvenido. Finalmente, igualmente espaciado de la cuadrícula no es probablemente la forma más precisa para hacer este tipo de integración.

Totalmente vectorizados aplicación es sencillo, 'integral' no tiene que ser una matriz 2D:

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

La presente respuesta se supone que ambas distribuciones son en realidad Binomial Negativa (NB), y también se supone que las dos muestras independientes.

En este marco, se tiene un modelo con cuatro parámetros, dos para cada distribución, y quiere probar una restricción con dos restricciones sobre los parámetros. Esto se puede hacer por la norma de Máxima Verosimilitud (ML) de la teoría. La prueba de razón de Verosimilitud (LR) es sencillo: el fitdistr función de la MASA utilizada anteriormente por Stephan Kolassa puede proporcionar el maximiza las probabilidades necesarias en la LR. Otra prueba que se conoce como la puntuación de la prueba o Multiplicador de Lagrange de la prueba (LM); tiene buenas características teóricas, sino que necesita de una matriz de información y una puntuación de vector de therestricted estimación ML. Una tercera posibilidad sería una prueba de Wald.

Deje que $ \boldsymbol{\psi}_{\textrm{c}} := [r_{\textrm{c}},\,\pi_{\textrm{c}}]'$ se el vector de parámetros para la "c" de la muestra, con la notación similar para la "w". Aquí $r$ representa size, mientras que $\pi$ es sinónimo de prob en R, es decir, la probabilidad de éxito. El pleno del vector de parámetros es $\boldsymbol{\theta} := [\boldsymbol{\psi}_{\textrm{c}}',\, \boldsymbol{\psi}_{\textrm{w}}']'$. Quieres poner a prueba la hipótesis $H_0:\,\boldsymbol{\psi}_{\textrm{c}} = \boldsymbol{\psi}_{\textrm{w}}$. La PELÍCULA de la prueba se basa en la distribución aproximada $$ \mathbf{U}'\mathbf{I}^{-1}\mathbf{U} \sim \chi^2(d) $$ donde el grado de libertad $d$ es el número de escalar restricciones, aquí $2$. El vector $\mathbf{U}$ y la matriz $\mathbf{I}$ son la puntuación de vectores y de la matriz de información en la estimación restringida, decir $\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$. La LR estadística de aquí puede ser obtenida por darse cuenta de que $\mathbf{I}$ es de bloque diagonal con dos $2\times 2$ bloques que conduce a dos contribuciones derivadas de las dos muestras. La puntuación y observa la información de la matriz puede ser utilizada en forma cerrada. Con el fitdistr función, la parametrisation difiere de la anterior y usa $\mu:= r\times (1-\pi)/\pi$. En lugar de un cambio de parámetro, podemos con un poco más de esfuerzo el uso de un concentrado de la log-verosimilitud.

Las dos pruebas que se tienen aquí muy pequeño $p$-valores por eso uno debe claramente rechazar la hipótesis nula.

El uso de simulaciones, parece que la información de la matriz de $\mathbf{I}$ puede en algunos casos ser mal condicionado, y la PELÍCULA de la prueba estadística puede entonces indebidamente a ser negativo. En este caso al menos, el LR test debe ser preferido para la PELÍCULA de la prueba.

## 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)