[R] rzinb (VGAM) and dnbinom in optim - Solved.

[Katharine Mullen]

Did you notice that the residual function I gave flipped what your
original example problem used as starting values for munb and size? so you
can sometimes afford to give a bad starting value for these param.

in trying to calculate the gradient L-BFGS-B may try evaluate fn with
parameter values outside of the box constraints; to categorically avoid
problems you'd have to know the max. step size for each parameter.  it
seems better to use a try statement to catch problems and impose a big
penalty if the try fails, as in

library(VGAM)
phi <- 0
munb <- 72
size <- 0.7

x <- rzinb(200, phi=phi, size=size, munb=munb)
## x should have some non-zero elements

fit <- vglm(x~1, zinegbinomial,trace=TRUE)

fncZiNB <- function(xx, x){
   phi <- xx[1]
   munb <- xx[2]
   size <- xx[3]
   yy <- try(dzinb(x=x,phi=phi, size=size, munb=munb,log=TRUE))
   if(class(yy) == "try-error" || !all(is.finite(yy))) 1e10 else -sum(yy)
}

# size, the 3rd param must be strictly > 0
 solut <- optim(fn=fncZiNB,par=c(0, 72, .7),x=x,method =
"L-BFGS-B",  lower=c(0,0,1e-10),upper=c(1,Inf,Inf),
control=list(trace=TRUE),hessian=TRUE)

##optim
solut$par
##vglm
Coef(fit)

On Sat, 19 Apr 2008, Tim Howard wrote:

> Thank you very much for the help and the great suggestion on cleaning up
> the residual function.
>
> Setting the bounds and being very careful with the starting position (on
> my real data) do seem to do the trick, which I wouldn't have guessed
> given the error message that was getting thrown.
>
> Best,
> Tim
>
> >>> Thomas Yee <t.yee at auckland.ac.nz> 4/18/2008 6:02 PM >>>
> Hi,
>
> using something like
>
> lower=c(1e-5,1e-5,1e-5),upper=c(1-1e-5,Inf,Inf)
>
> is even more safer since it sometime evaluates the function slightly
> outside the lower and upper limits.
>
> cheers
>
> Thomas
>
>
>
> Katharine Mullen wrote:
>
> >I'm not going to look into what's happening in-depth but it looks like the
> >bounds for your parameters need to be set with care; the below (with
> >slight re-def. of your residual function) gives results that seem to match
> >vglm approximately, with 1e-10 standing in for a bound of 0 and 1-e-10
> >standing in for a bound of 1.
> >
> >library(VGAM)
> >phi <- 0.09
> >size <- 0.7
> >munb <- 72
> >x <- rzinb(200, phi=phi, size=size, munb=munb)
> >
> >fit <- vglm(x~1, zinegbinomial,trace=TRUE)
> >
> >fncZiNB <- function(xx, x){
> >   phi <- xx[1]
> >   size <- xx[2]
> >   munb <- xx[3]
> >   -sum(dzinb(x=x,phi=phi, size=size, munb=munb,log=TRUE))
> >}
> >
> > solut <- optim(fn=fncZiNB,par=c(0.09, 72, 0.7),x=x,method =
> >"L-BFGS-B",  lower=c(1e-10,1e-10,1e-10),upper=c(1-1e-10,Inf,Inf),
> >control=list(trace=TRUE),hessian=TRUE)
> >
> >##optim
> >solut$par
> >##vglm
> >Coef(fit)
> >
> >On Fri, 18 Apr 2008, Tim Howard wrote:
> >
> >
> >
> >>Dear R-help gurus (and T.Yee, the VGAM maintainer) -
> >>I've been banging my head against the keyboard for too long now, hopefully someone can pick up on the errors of my ways...
> >>
> >>I am trying to use optim to fit a zero-inflated negative binomial distribution.  No matter what I try I can't get optim to recognize my initial parameters. I think the problem is that dnbinom allows either 'prob' or 'mu' and I am trying to pass it 'mu'.  I've tried a wrapper function but it still hangs.  An example:
> >>
> >>#set up dummy data,load VGAM, which is where I am getting the zero-inflated neg binom
> >>#I get the same problem without VGAM, using dnbinom in the wrapper instead.
> >>library(VGAM)
> >>phi <- 0.09
> >>size <- 0.7
> >>munb <- 72
> >>x <- rzinb(200, phi=phi, size=size, munb=munb)
> >>
> >>#VGAM can fit using vglm... but I'd like to try optim
> >>
> >>
> >>>fit <- vglm(x~1, zinegbinomial,trace=TRUE)
> >>>
> >>>
> >>VGLM    linear loop  1 :  loglikelihood = -964.1915
> >>VGLM    linear loop  2 :  loglikelihood = -964.1392
> >>VGLM    linear loop  3 :  loglikelihood = -964.1392
> >>
> >>
> >>>Coef(fit)
> >>>
> >>>
> >>       phi       munb          k
> >> 0.1200317 62.8882874  0.8179183
> >>
> >>
> >>
> >>>#build my wrapper function for dzinb
> >>>fncZiNB <- function(phi, size, munb){
> >>>
> >>>
> >>+         -sum(dzinb(x=x,phi=phi, size=size, munb=munb,log=TRUE))}
> >>
> >>
> >>>#test the wrapper, seems to work.
> >>>fncZiNB(phi=0.09, size=0.7, munb=72)
> >>>
> >>>
> >>[1] 969.1435
> >>
> >>#try it in optim
> >>
> >>
> >>>solut <- optim(fn=fncZiNB,par=c(phi=0.09, munb=72, size=0.7),method = "L-BFGS-B", hessian=TRUE)
> >>>
> >>>
> >>Error in dnbinom(x = x, size = size, prob = prob, log = log.arg) :
> >>  argument "size" is missing, with no default
> >>
> >>I have the same problem with dnbinom.
> >>I'm sure I'm doing something obvious.... any suggestions?
> >>Session info below. Thanks in advance.
> >>Tim Howard
> >>New York Natural Heritage Program
> >>
> >>
> >>
> >>>sessionInfo()
> >>>
> >>>
> >>R version 2.6.2 (2008-02-08)
> >>i386-pc-mingw32
> >>
> >>locale:
> >>LC_COLLATE=English_United States.1252;LC_CTYPE=English_United States.1252;LC_MONETARY=English_United States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252
> >>
> >>attached base packages:
> >>[1] stats4    splines   stats     graphics  grDevices utils     datasets  methods
> >>[9] base
> >>
> >>other attached packages:
> >>[1] VGAM_0.7-6          randomForest_4.5-23
> >>
> >>______________________________________________
> >>R-help at r-project.org mailing list
> >>https://stat.ethz.ch/mailman/listinfo/r-help
> >>PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> >>and provide commented, minimal, self-contained, reproducible code.
> >>
> >>
> >>
>
>