r - 将 SPCM-ZIP 模型从 SAS 中的 PROC NLMIXED 转换为 R

Question

我正在尝试将平滑参数计数模型零膨胀泊松从在 SAS 中使用 PROC NLMIXED 转换为 R 中的东西。这个模型和 SAS 代码来自论文https://onlinelibrary.wiley.com/doi/abs/10.1111/乔斯.12280。

我拥有的 SAS 代码：

proc nlmixed data=a maxiter=10000 tech=dbldog gconv=0 HESS method=gauss subgrad=gradient;
paramters a0 0 a1 0 g0 0 g1 0 b0 0 b1 0 b2 0 b3 0 d0 0 d1 0 d2 0 d3 0 gamma 0 c 0;
bounds 0<=gamma<=100, &minwz<=c<=&maxwz;
G = transition function based on gamma, c, and spatial weight variable
linpinf = zero inflated regression
infprob = 1/(1+exp(-linpinfl))
lambda = exp(count regression)
if &yvar=0 then
  ll = log(infprob + (1-infprob*exp(-lambda));
else ll = log((1-infprob)) - lambda + &yvar*log(lambda) - lgamma(&yvar + 1);
model &yvar ~ general(ll);
run;

我没有写出所有的回归和转换函数以节省空间。

在 R 中，我已经能够使用 optim 来复制标准 ZIP 模型，但是一旦我添加了转换函数，我就没有得到类似的系数估计值。

我的 R 代码：

spcmzip <- function(beta,a0,a1,g0,g1,x0,x1,x2,d0,d1,d2,y) {  
gamma=beta[11]
c=beta[12]

G <- NA
if (gamma >= 0 && gamma <= 100 && c >= minwz && c<= maxwz) {
  G <- 1/(1+exp(-gamma*(wshale_tight-c)/stdwz))
}

lin      <- beta[5]*x0+beta[6]*x1+beta[7]*x2+beta[8]*d0*G+beta[9]*d1*G+beta[10]*d2*G
linpinfl <- beta[1]*a0+beta[2]*a1+beta[3]*g0*G+beta[4]*g1*G
infprob  <- 1/(1+exp(-linpinfl))
lambda   <- exp(lin)

yt <- c()
for (i in 1:length(y)) {
  if (y[i]==0) {
    yt[i] <- log(infprob[i] + (1-infprob[i])*exp(-lambda[i]))
  } else {
    yt[i] <- log((1-infprob[i])) - lambda[i] + y[i]*log(lambda[i]) - lgamma(y[i]+1)
  }
}

logl <- (sum(yt))
return (-logl) 
}

mle <- optim(c(0,0,0,0,0,0,0,0,0,0,0,0),spcmzip,a0=a0,a1=a1,g0=g0,g1=g1,x0=x0,x1=x1,
         x2=x2,d0=d0,d1=d1,d2=d2,
         y=y,method="Nelder-Mead",hessian=FALSE)

但是，使用“Nelder-Mead”方法并选择不使用粗麻布是我可以对其进行估计的唯一方法，并且系数与 SAS 相去甚远。我认为问题的一部分是 optim 没有进行双折线优化，这可能是为什么这不起作用的一个重要部分。

人们知道 R 中有没有可以解决这种优化问题的包？还是有人看到我的代码明显有问题？

标签： rsas