以下是@BenBolker答案中的一個變體,參數用均值和方差來表示。
如果您將對數似然函數作爲均值和方差的函數,將平均值表示爲協變量的線性函數,並使用optim()
來獲得MLE和Hessian,那麼您可以得到類似GLM的東西。
均值爲mu1-mu2
,方差爲mu1+mu2
。的兩個參數可以被寫爲平均值和方差,即功能:
mu1 <- (mn+va)/2
mu2 <- (va-mn)/2
的約束是mu1
和mu2
都是正的。爲了達到這個目的,Skellam的平均值必須小於方差。這表明重新參數方差爲:
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和std。迴歸係數的錯誤:
fit$par[1:npars] # MLE
sqrt(diag(solve(fit$hessian)))[1:npars] # std error
我會補充說,包括隨機效應,有可能使用MCMC與任何答案參數化,與像STAN或JAGS預先包裝的採樣器(你會需要在JAGS中使用ones trick)。最難的部分可能是在Skellam PMF中移植貝塞爾函數。
我想指出'Skellam'包現在迴歸功能有限。請參閱:http://stats.stackexchange.com/questions/264252/is-there-a-reason-why-the-regression-in-the-r-skellam-package-uses-three-optimis – Alex