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

[Tim Howard]

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