3

Is there an easy way to fit a multivariate regression in R in which the dependent variable is distributed in accordance with the Skellam distribution (difference between two Poisson-distributed counts)? Something like:

myskellam <- glm(A ~ B + C + D, data = mydata, family = "skellam")

This should accommodate fixed effects. But ideally, I would prefer random effects as I understand that fixed effects may introduce measurement biases. Therefore I guess the ideal solution should be using the lme4 or glmmADMB package.

Alternatively, is there a way to transform the data to apply more usual regression tools?

1

2 Answers 2

6

Incomplete answer, but seems a bit more than a comment.

Mixed effects seem hard; you could do it with AD Model Builder or Template Model Builder, both of which have built-in facilities for Laplace approximation. For fixed effects you can use something like

library("skellam")
library("bbmle")

Reparameterize dskellam(x, lambda1, lambda2) to a form that is essentially location (geometric mean lambda=gmlambda=sqrt(lambda1*lambda2)) and shape (difference in lambdas: ldiff=sqrt(lambda1/lambda2) (so that lambda1=gmlambda*ldiff, lambda2=gmlambda/ldiff).

 dskellam2 <- function(x, gmlambda, ldiff, log=FALSE) {
     dskellam(x,gmlambda*ldiff,gmlambda/ldiff,log=log)
 }

Then something like this should work:

 mle2(A~dskellam2(gmlambda=exp(logmu),ldiff=exp(logs), data=mydata,
      parameters=list(logmu~B+C+D),
      start=list(logmu=0,logs=0)))

... but it could take some fussing to make it work.

1
  • Thanks very much. I will try that although I am a bit afraid of the "fussing" part.
    – bdu
    Sep 7, 2015 at 9:37
3

Here's an variant on @BenBolker's answer, parameterizing in terms of the mean and variance.

You could get something GLM-like if you write the log-likelihood as a function of the mean and variance, express the mean as a linear function of covariates, and use optim() to get the MLE and Hessian.

The mean is mu1-mu2, the variance is mu1+mu2. The two parameters can be written as functions of the mean and variance, ie:

mu1 <- (mn+va)/2
mu2 <- (va-mn)/2

The constraint is that mu1 and mu2 are positive. To achieve this, the mean of the Skellam must be less than the variance. This suggests to reparameterize the variance as:

va <- max(abs(mn)) + va_st

And treat va_st as the parameter to estimate.

Putting that all together:

library(skellam)
logL_Skellam <- function(pars, X, Y){
    mn <- X %*% pars[1:ncol(X)]
    va_st <- exp(pars[ncol(X)+1]) # constrain to be positive
    va <- max(abs(mn)) + va_st
    # convert to parameters of skellam
    mu1 <- (mn+va)/2
    mu2 <- (va-mn)/2
    # optim minimizes so return negative LL
    -sum(dskellam(Y, mu1, mu2, log=T)) 
}

Optimise:

# simulated data
set.seed(103)
npars <- 3
nobs <- 300
X <- cbind(1, matrix(rnorm(nobs*(npars-1)), nrow=nobs))
beta <- rnorm(npars)
va <- max(abs(X%*%beta)) + 3.3
Y <- rskellam(nobs, (X%*%beta+va)/2, (va-X%*%beta)/2)

# fit
fit <- optim(c(0,0,0,5), logL_Skellam, X=X, Y=Y, hessian=T)

Taking care that optim actually converges. The MLE and std. errors of the regression coefficients:

fit$par[1:npars] # MLE
sqrt(diag(solve(fit$hessian)))[1:npars] # std error

And I'll add that to include random effects, it may be possible to use MCMC with the parameterizations in either answer, with a prepackaged sampler like STAN or JAGS (you'd need to use the ones trick in JAGS). The hardest part to this may be porting the Bessel function in the Skellam PMF.

1
  • Thanks very much too, I will try this as an alternative and keep you posted.
    – bdu
    Sep 7, 2015 at 9:38

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.