Как я могу соответствовать регрессии Скеллама?

Неполный ответ, но кажется, что это больше, чем комментарий.

Смешанные эффекты кажутся сложными; вы можете сделать это с помощью AD Model Builder или Конструктор моделей шаблонов, оба из которых имеют встроенные средства для аппроксимации Лапласа. Для фиксированных эффектов вы можете использовать что-то вроде

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

Перепараметризируйте dskellam(x, lambda1, lambda2) в форме, которая по существу является местоположением (среднее геометрическое лямбда = gmlambda= sqrt(lambda1*lambda2)) и формой (разница в лямбдах: ldiff=sqrt(lambda1/lambda2) (так что lambda1=gmlambda*ldiff, lambda2=gmlambda/ldiff).

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

Тогда что-то вроде этого должно работать:

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

... но может потребоваться некоторая возня, чтобы заставить его работать.

Вот вариант ответа @BenBolker, параметризующий с точки зрения среднего значения и дисперсии.

Вы могли бы получить что-то похожее на GLM, если бы вы записали логарифмическое правдоподобие как функцию среднего значения и дисперсии, выразили среднее значение как линейную функцию ковариат и использовали optim() для получения MLE и гессиана.

Среднее mu1-mu2, дисперсия mu1+mu2. Два параметра могут быть записаны как функции среднего значения и дисперсии, т.е.:

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

Ограничение состоит в том, что mu1 и mu2 положительны. Для этого среднее значение Скеллама должно быть меньше дисперсии. Это предлагает перепараметризовать дисперсию как:

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

И относитесь к va_st как к оцениваемому параметру.

Собираем все это вместе:

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

Оптимизировать:

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

Заботясь о том, чтобы optim действительно сходится. MLE и станд. ошибки коэффициентов регрессии:

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

И я добавлю, что для включения случайных эффектов может быть возможно использовать MCMC с параметризациями в любом ответе, с предварительно упакованным сэмплером, таким как STAN или JAGS (вам нужно будет использовать уловка в JAGS). Самой сложной частью этого может быть портирование функции Бесселя в PMF Skellam.