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