Конвертировать модель SPCM-ZIP из PROC NLMIXED в SAS в R

Question

Я пытаюсь преобразовать сглаженный пуассон модели с гладким числом параметров из использования PROC NLMIXED в SAS во что-то в R. Эта модель и код SAS взяты из статьи https://onlinelibrary.wiley.com/doi/abs/10.1111/jors.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, которые могут решить такую ​​проблему оптимизации?Или кто-то видит что-то явно не так с моим кодом?