Skip to content
Prev 310599 / 398506 Next

maximum likelihood estimation in R

David Winsemius <dwinsemius <at> comcast.net> writes:
Copying and pasting the relevant code from Nabble:

library(GB2) 
library(bbmle)

lgb1=function(a,b0,b2,b3,b4,p,q){
 xb=b0+b2*fsex1[,2]+b3*fvtype1[,2]+b4*fvuse1[,2]
ll=sum(log(dgb2(loss1,a,exp(xb),p,q)))
 return(-ll)
 }

start=list(a=3.1,b0=2.5,
b2=.2,b3=1,b4=-.5,
p=7.2,q=.3)

mle2(lgb1,start)->fit1

summary(fit1)
Maximum likelihood estimation

Call:
mle2(minuslogl = lgb1, start = start)

Coefficients:
      Estimate  Std. Error     z value     Pr(z)    
a   3.0747e+00  6.4741e-01  4.7492e+00 2.042e-06 ***
b0  2.5327e+00  4.6887e-01  5.4016e+00 6.605e-08 ***
b2  2.0000e-01  3.9686e-11  5.0396e+09 < 2.2e-16 ***
b3  1.0000e+00  7.6565e-12  1.3061e+11 < 2.2e-16 ***
b4 -5.0000e-01  1.5312e-11 -3.2653e+10 < 2.2e-16 ***
p   7.1281e+00  8.0269e+00  8.8800e-01    0.3745    
q   3.5098e-01  8.6902e-02  4.0388e+00 5.372e-05 ***
---
Signif. codes:  0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1 

-2 log L: 10137.56


  I agree that it looks like the parameters got stuck at
their initial values, but without the data it's really hard to
know what went wrong.  I assume you started with reasonable
parameter values and did some sanity checking on the data ... ?

  You can use the 'trace' and 'browse_obj' arguments to 
get more detail about what's going on.

  For what it's worth, with a little more effort you could
use the formula interface to mle2: something like

## set up a "d*" function with a log argument
  dgb2B <- function(...,log=FALSE) {
     r <- dgb2(...)
     if (log) log(r) else r
  }

## set up a data frame to hold the data
dframe <- data.frame(loss1,sex=fsex1[,2],type=fvtype1[,2],use=fvuse1[,2])

mle2(loss1~dgb2B(shape1=a,scale=exp(logscale),shape2=p,shape3=q),
   parameters=list(logscale~sex+type+use),
   data=dframe,
   start=...)

[see the help page for the details of how "start" is specified
in this case]