0

Why fitting a correlated shared gamma frailty model in R, I obtained the negative estimate for the variance of parameters of interest. I used optim, mle and mle2 with L-BFGS-B and Nelder-Mead optimization methods (just to see if there are different outputs) with different starting values. However, all of them gave me some negative values for the variance. An example of outputs is:

    Coefficients:
       Estimate   Std. Error
       a1      0.008622576          NaN
       b1      0.074853595          NaN
       a2      0.002150580 0.0001168565
       b2     -0.007083631 0.0019890865
       alpha1  0.599927849          NaN
       alpha2  1.099990231          NaN

       -2 log L: 5684.66

My questions are: 1. Why do I get the negative variance? Is it common? 2. Are parameter estimates reliable when I have negative variance. Since I see that they do not move away from the starting values (for those parameters with negative one). 3. If it is a problem, if there is any strategy to deal with it? Thank you a lot. Here is my code:

    LoglikU  <- function(a1,b1,a2,b2,alpha1,alpha2){
   S1  <- (1+(1/alpha1)*a1/b1*(exp(b1*age)-1))^{-alpha1}
   S2  <- (1+(1/alpha2)*a2/b2*(exp(b2*age)-1))^{-alpha2}
   S12 <- S1*S2
   loglik <- -base::sum(nn*log(S12) + ny*log(S1-S12) +yn*log(S2-S12) + yy*log(1-S1-S2+S12))
   return(loglik)

}

   estimU <- mle2(LoglikU, start=list(a1=0.008, b1=0.077,a2=0.002,b2=-0.00002,alpha1=0.6,alpha2=1),
           skip.hessian = F, method = "L-BFGS-B", 
           lower =c(a1=1e-6,b1=-Inf,a2=1e-6,b2=-Inf, alpha1=1e-6,alpha2=1e-6),
           upper = c(a1=+Inf,b1= +Inf,a2= +Inf,b2= +Inf,alpha1=+Inf,alpha2=+Inf),
           control = list(trace=TRUE, REPORT=500))

My data are:

    list( age = c( 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 
   18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 
   38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 
   58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 
   78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 
   98, 99, 100) , 
   nn = c( 4, 11, 8, 27, 38, 61, 81, 50, 41, 35, 40, 31, 29, 27, 
   26, 29, 19, 21, 18, 18, 32, 39, 38, 40, 41, 26, 41, 49, 35, 35, 28,    29, 28, 32, 
   31, 31, 26, 18, 24, 22, 22, 22, 29, 19, 18, 20, 22, 21, 17, 19, 19, 14, 11, 7, 
   6, 16, 8, 12, 9, 8, 8, 4, 6, 6, 4, 4, 6, 4, 2, 1, 3, 3, 5, 2, 3, 3, 1, 1, 2, 0, 
   2, 3, 1, 2, 0, 1, 0, 1, 0, 0, 1, 1, 2, 0, 0, 0, 0, 0, 0, 0, 1) , 
   yn = c( 0, 
   0, 0, 2, 3, 1, 8, 7, 8, 4, 5, 1, 3, 8, 2, 3, 3, 6, 6, 2, 7, 15, 16, 10, 15, 13, 
   12, 25, 10, 15, 18, 27, 19, 29, 30, 28, 26, 27, 32, 27, 41, 24, 40, 43, 47, 31, 
   39, 57, 38, 37, 45, 42, 32, 51, 67, 57, 62, 62, 42, 54, 41, 28, 20, 28, 17, 26, 
   18, 27, 30, 21, 27, 26, 9, 19, 10, 12, 15, 12, 11, 17, 10, 16, 9, 9, 10, 10, 5, 
   5, 5, 3, 1, 1, 1, 1, 0, 0, 0, 1, 0, 0, 0) , 
   ny = c( 0, 0, 0, 0, 0, 0, 0, 0, 
   0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 3, 0, 2, 1, 0, 1, 0, 0, 0, 
  1, 0, 2, 3, 2, 1, 1, 0, 1, 2, 0, 0, 0, 0, 0, 2, 0, 0, 3, 1, 1, 0, 1, 1, 0, 1, 
  1, 0, 1, 0, 0, 1, 1, 0, 0, 2, 0, 1, 0, 0, 1, 0, 0, 0, 2, 1, 0, 1, 0, 0, 1, 0, 
  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0) , 
  yy = c( 0, 0, 0, 0, 0, 0, 1, 
  2, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 3, 2, 4, 1, 2, 3, 4, 2, 2, 2, 
  4, 3, 1, 9, 4, 3, 4, 3, 4, 3, 8, 1, 7, 4, 4, 4, 10, 0, 2, 3, 3, 3, 4, 5, 3, 4, 
 2, 2, 3, 2, 2, 4, 1, 2, 4, 1, 0, 2, 2, 1, 2, 3, 3, 3, 3, 3, 5, 3, 3, 2, 1, 2, 
 0, 0, 0, 0, 3, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0) )
kjetil b halvorsen
  • 63,378
  • 26
  • 142
  • 467
bienco88
  • 61
  • 1
  • 6
  • 2
    There seems to be no assurance that your objective function is even defined--it looks entirely possible for it to attempt to take logs of negative values. The parameterization looks unstable, too. Are you sure this is actually a log likelihood? Have you checked that the routine actually converged to what appears to be a maximum? – whuber Mar 23 '17 at 15:55
  • Hi Whuber, Yes I am sure this is the log-likelihood function for the type of univariate gamma frailty, a special case of correlated gamma frailty model with two gamma distributions for frailty terms. The data above was analysed using SAS. My parameter estimates are some what close to what obtained using SAS. Only the problem I get is that sometimes I got negative estimate for the variance. – bienco88 Mar 23 '17 at 19:19
  • OK--*but how do you know the routine has converged to the maximum likelihood?* – whuber Mar 23 '17 at 19:21
  • I had run the analysis on some simulation data set. Do you have any better suggestion how to see it? – bienco88 Mar 23 '17 at 19:54
  • Yes: read the help page for `mle2` and inspect its return value for the codes and other information put there that tell you about whether it converged and how well. Even then there are no guarantees: see http://stats.stackexchange.com/questions/7308 for a case in point. The answers illustrate some techniques for chasing down apparent anomalies in optimization procedures. – whuber Mar 23 '17 at 21:05
  • 1
    Hi @whuber, I tried to use profile to plot my likelihood but it gave me errors: Error in as.double(y) : cannot coerce type 'S4' to vector of type 'double'. I am still working on it to fix. In the meanwhile, I used optimizer = "nlmimb" and the estimates of variance are ok new. The parameter estimates do not depend on the initial values any more. – bienco88 Mar 24 '17 at 10:55

0 Answers0