我該如何適應Skellam迴歸？

是否有一種簡單的方法來擬合R中的多變量回歸，其中因變量根據Skellam distribution（兩個泊松分佈計數之間的差異）分佈？類似於：我該如何適應Skellam迴歸？

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

這應該適應固定效果。但理想情況下，我更喜歡隨機效應，因爲我知道固定效應可能會引入測量偏差。因此，我想理想的解決方案應該使用lme4或glmmADMB包。

或者，有沒有辦法將數據轉換爲應用更常用的迴歸工具？

來源

2015-09-03 bdu

我想指出'Skellam'包現在迴歸功能有限。請參閱：http://stats.stackexchange.com/questions/264252/is-there-a-reason-why-the-regression-in-the-r-skellam-package-uses-three-optimis – Alex

不完整的答案，但似乎不僅僅是評論。

混合效應似乎很難;你可以用AD Model Builder或Template Model Builder來完成，兩者都有內置的拉普拉斯近似設施。爲固定效應，你可以使用像

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

重新參數化dskellam(x, lambda1, lambda2)到的形式本質上是位置（幾何平均波長= gmlambda = sqrt(lambda1*lambda2)）和形狀（lambda表達式差：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)))

...但它可能需要一些小事才能使其工作。

來源

2015-09-05 01:16:20

非常感謝。我會嘗試，儘管我有點害怕「fu」「部分。 – bdu

以下是@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中移植貝塞爾函數。

來源

2015-09-05 16:08:29

非常感謝，我會嘗試作爲替代方案，並保持張貼。 – bdu

我該如何適應Skellam迴歸？

回答

相關問題