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

From: Tim Howard <tghoward_at_gw.dec.state.ny.us>
Date: Sat, 19 Apr 2008 10:25:55 -0400

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.
>>
>>
>>



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. Received on Sat 19 Apr 2008 - 14:28:56 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Sat 19 Apr 2008 - 17:30:30 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.

list of date sections of archive