通过 IWLS 实现具有两个过离散参数的负二项式答案

【问题标题】：Implementation of a Negative binomial with two overdispersion parameter through IWLS通过 IWLS 实现具有两个过离散参数的负二项式
【发布时间】：2020-11-16 12:16:41
【问题描述】：

理论

受Lindén and Mäntyniemi 的论文启发，使用两个过度离散参数来模拟线性和！二次均方差关系。请参见以下等式。 theta = 0 的情况对应于准泊松，omega = 1 的情况对应于正则负二项式。

sigma^2 = omega * mu + theta * mu^2

我主要跟随Ver Hoef and Boveng通过迭代加权最小二乘法（IWLS）来实现它。

实施 IWLS

我可以测试以下代码的准泊松和负二项式（欧米茄固定为 1），但无法验证它是否适用于欧米茄不等于 1 的负二项式。

IWLS <- function(Omega1,Theta1,tollvl=1,X=X){
X <- as.matrix(X) # model.matrix(), factor variables only
mu <- rep(mean(Y),length(X[,1]))
Omega <- Omega1
Theta <- Theta1
Y.head<-Y # 152 datapoints
tol=10
i=1

while(tol>tollvl){ # Iterative
W<-diag(c(mu/(Omega[i] + Theta[i]*mu)) ) 
# Theta = 0 -> Quasipoisson, Omega = 1 -> negbin, both unfixed (reason wrote this code)
Beta.head<-solve(t(X)%*%W%*%X,t(X), tol=9^-45)%*%W%*%Y.head
nu <- X%*%Beta.head
mu <- exp(nu)
Y.head <- nu + (Y-mu)/mu # derivative of log(nu) -> 1/nu
if(Theta1==0){
    Theta <- c(Theta,0)
    Omega <- c(Omega,mean( (Y-mu)^2/mu) )
} else{ 
    if(Omega1==1){
    Theta<- c(Theta,mean( (( (Y-mu)^2) - mu )/mean(mu)^2))  
    Omega <- c(Omega,1)
    } else{ # NB2
    ThetaF  <- mean( (Y-mu)^2 - mu )/(mean(mu)^2) 
    Omega <- c(Omega,mean( ((Y-mu)^2)/mu - mu*ThetaF) )
    Theta  <- c(Theta,mean( ((Y-mu)^2 - mu*Omega[i+1] )/(mean(mu)^2 )) )
    }
}
tol <- abs(Omega[i+1]-Omega[i])+abs(Theta[i+1]-Theta[i])
i <- i+1
}
list(Omega=Omega,
     Theta=Theta,
     Beta.head=Beta.head,
     SE = sqrt(diag(solve(t(X)%*%W%*%X))),
     FittedValues = mu)
}

问题/问题

主要关注的是以下行，它指定了双重过度分散的负二项式。我没有找到支持该代码的文献。我在这里尝试了一个 for() 循环以及关于平均值的不同方法，因为对每个 mu 取平均值或在最后取平均值对结果有很大影响。

else{ # NB2
    ThetaF  <- mean( (Y-mu)^2 - mu )/(mean(mu)^2) 
    Omega <- c(Omega,mean( ((Y-mu)^2)/mu - mu*ThetaF) )
    Theta  <- c(Theta,mean( ((Y-mu)^2 - mu*Omega[i+1] )/(mean(mu)^2 )) )
    }

此版本的代码适用于相对少量的协变量（与其他两个模型相比，Beta.head 相同，只是 SE 不同），其中模型矩阵的维度约为 152:5。对于更多列/多个因子变量，该算法无法使用 {In sqrt(diag(solve(t(X) %% W %% X)) 计算每个属性的标准误差) : NaN 产生}，表示多重线性。 Theta 和 omega 有时不会收敛，但即使它们收敛，也会出现不止一个变量的此消息。

要运行的测试代码：

set.seed(123)
Y.Test<-sample(1:100,152,replace = T)*rpois(152,1)+1 # Try to model true data
XXX<-model.matrix(Y.Test~as.factor(sample(1:4,152,replace = T))+as.factor(sample(1:3,152,replace = T)))
Neg.Bin2.Test<-IWLS(Omega1 =var(Y)/mean(Y),Theta1 = var(Y)/mean(Y)^2 - 1/mean(Y), tollvl =  10^-10,X=XXX)

如果多线性是一个问题，为什么它只对双重过度分散模型有问题。 omega=1 修复解决了这个问题。

【问题讨论】：

标签： r

【解决方案1】：

我想通了。有多个问题。

忘记了列表输出中的容差水平，因此：

SE = sqrt(diag(solve(t(X)%*%W%*%X, tol=9^-45)))

theta 的计算还有其他问题

去掉一个均值(mu)

else{ 
    if(Omega1==1){ # Negative binomial Case
    Theta<- c(Theta,mean( ((Y-mu)^2 - mu )/mu^2 ) ) # no mean(mu)
    Omega <- c(Omega,1) 
    }

这里我不得不添加一个额外的步骤

else{ # Negative binomial with 2 overdispersion parameters
    OmegaF <- mean( ((Y-mu)^2)/mu )
    ThetaF  <- mean( ((Y-mu)^2 - mu )/mu^2 ) 
    Omega <- c(Omega,mean( ((Y-mu)^2)/mu - mu*ThetaF) )
    Theta  <- c(Theta,mean( ((Y-mu)^2 - mu*Omega[i+1] )/mu^2 )) 
    }
}

【讨论】：