Lorenzo Cattarino
2010-Nov-09 05:41 UTC
[R] convergence message & SE calculation when using optim( )
Hi R-users,
I am trying to estimate function parameters using optim(). My count
observations follows a Poisson like distribution. The problem is that I
wanna express the lambda coefficient, in the passion likelihood
function, as a linear function of other covariates (and thus of other
coefficients). The codes that I am using (except data frame) are the
following (FYI the parameters need to be positive):
myfun <- function(coeff, H1, H2, p, Range)
{
(coeff[1]+coeff[2]*H1+coeff[3]*H2+coeff[4]*p+H1*Range*coeff[5]+H2*Range*
coeff[6]+H1*H2*coeff[7])*exp((coeff[8]+coeff[9]*H1+coeff[10]*H2+coeff[11
]*p+H1*Range*coeff[12]+H2*Range*coeff[13]+H1*H2*coeff[14])*(Range-1))+co
eff[15]+coeff[16]*H1+coeff[17]*H2+coeff[18]*p+H1*Range*coeff[19]+H2*Rang
e*coeff[20]+H1*H2*coeff[21]
}
SS <- function(coeff,Range,H1,H2,p,steps)
{
sum((steps - myfun(coeff,Range,H1,H2,p))^2)
}
coeff <-
c(0.1,0.1,0.1,0.1,0.1,1,5,5,5,1,1,1,1,1,1,1,0.1,0.1,0.1,0.1,0.1)
scale <-
c(0.1,0.1,0.1,0.1,0.1,1,5,5,5,1,1,1,1,1,1,1,0.1,0.1,0.1,0.1,0.1)
est_coeff <- optim(par=coeff, fn = SS, H1=org_results$H1,
H2=org_results$H2, p=org_results$p, Range=org_results$Range,
steps=org_results$no.steps, method= 'L-BFGS-B', lower = rep(0, 21),
upper = rep(Inf, 21), control = list(trace=FALSE, parscale=scale),
hessian=TRUE)
this is the output:
$par
[1] 0.099794607 0.099841098 0.099896127 0.099899549 0.099856776
0.991269412 4.807153280
[8] 25.115556187 55.961674737 1.519658195 1.378913148 2.800328223
1.448902455 2.280837645
[15] 1.594648898 0.011581676 0.040651369 0.000000000 0.000000000
0.000000000 0.002717246
$value
[1] 14535187
$counts
function gradient
54 54
$convergence
[1] 0
$message
[1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"
$hessian
[,1] [,2] [,3] [,4]
[,5] [,6]
[1,] 0.000000e+00 0.000000e+00 0.000000e+00 -0.0023283064
0.0000000000 2.328306e-03
[2,] 0.000000e+00 0.000000e+00 0.000000e+00 -0.0023283064
0.0000000000 2.328306e-03
[3,] 0.000000e+00 0.000000e+00 0.000000e+00 0.0023283064
0.0000000000 0.000000e+00
[4,] -2.328306e-03 -2.328306e-03 2.328306e-03 0.0000000000
-0.0023283064 2.328306e-04
[5,] 0.000000e+00 0.000000e+00 0.000000e+00 -0.0023283064
0.0046566129 -2.095476e-03
[6,] 2.328306e-03 2.328306e-03 0.000000e+00 0.0002328306
-0.0020954758 0.000000e+00
[7,] 9.313226e-05 9.313226e-05 4.656613e-05 0.0023748726
0.0001396984 4.656613e-05
[8,] -4.284550e-01 -4.284550e-01 -1.916196e-01 -0.2018641680
-0.3814231604 -1.689885e-01
[9,] -4.284550e-01 -4.284550e-01 -1.916196e-01 -0.2018641680
-0.3814231604 -1.689885e-01
[10,] -1.892913e-01 -1.892913e-01 -1.240987e-01 -0.0880099833
-0.1706648618 -1.098961e-01
[11,] -1.974404e-01 -1.974404e-01 -8.568168e-02 -0.1315493137
-0.1781154424 -7.823110e-02
[12,] -3.885943e-01 -3.885943e-01 -1.704320e-01 -0.1802109182
-0.3480818123 -1.513399e-01
[13,] -1.706649e-01 -1.706649e-01 -1.108274e-01 -0.0789295882
-0.1504085958 -9.872019e-02
[14,] -1.892913e-01 -1.892913e-01 -1.240987e-01 -0.0880099833
-0.1706648618 -1.098961e-01
[15,] 2.813991e+00 2.813991e+00 1.214910e+00 1.3138633221
2.5336630642 1.093373e+00
[16,] 2.814224e+00 2.814224e+00 1.212582e+00 1.3138633221
2.5336630642 1.093373e+00
[17,] 1.213048e+00 1.213048e+00 7.892959e-01 0.5634501576
1.0943040252 7.092021e-01
[18,] 1.317821e+00 1.317821e+00 5.704351e-01 0.8870847523
1.1851079762 5.112961e-01
[19,] 2.537854e+00 2.537854e+00 1.089647e+00 1.1827796698
2.2817403078 9.855721e-01
[20,] 1.094304e+00 1.094304e+00 7.101335e-01 0.5122274160
0.9802170098 6.388873e-01
[21,] 1.215376e+00 1.215376e+00 7.939525e-01 0.5657784641
1.0919757187 7.094350e-01
[,7] [,8] [,9] [,10] [,11]
[,12] [,13]
[1,] 9.313226e-05 -0.42845495 -0.42845495 -0.18929131 -0.19744039
-0.38859434 -0.17066486
[2,] 9.313226e-05 -0.42845495 -0.42845495 -0.18929131 -0.19744039
-0.38859434 -0.17066486
[3,] 4.656613e-05 -0.19161962 -0.19161962 -0.12409873 -0.08568168
-0.17043203 -0.11082739
[4,] 2.374873e-03 -0.20186417 -0.20186417 -0.08800998 -0.13154931
-0.18021092 -0.07892959
[5,] 1.396984e-04 -0.38142316 -0.38142316 -0.17066486 -0.17811544
-0.34808181 -0.15040860
[6,] 4.656613e-05 -0.16898848 -0.16898848 -0.10989606 -0.07823110
-0.15133992 -0.09872019
[7,] 9.313226e-05 -0.18775463 -0.18775463 -0.12246892 -0.08759089
-0.16880222 -0.11008233
[8,] -1.877546e-01 0.12330711 0.12330711 0.07711351 0.05806796
0.11092052 0.06952323
[9,] -1.877546e-01 0.12330711 0.12330711 0.07711351 0.05806796
0.11092052 0.06952323
[10,] -1.224689e-01 0.07711351 0.07711351 0.05681068 0.03585592
0.06938353 0.05122274
[11,] -8.759089e-02 0.05806796 0.05806796 0.03585592 0.03864989
0.05215406 0.03213063
[12,] -1.688022e-01 0.11092052 0.11092052 0.06938353 0.05215406
0.09965152 0.06286427
[13,] -1.100823e-01 0.06952323 0.06952323 0.05122274 0.03213063
0.06286427 0.04656613
[14,] -1.224689e-01 0.07711351 0.07711351 0.05681068 0.03585592
0.06938353 0.05122274
[15,] 1.214491e+00 -0.79902820 -0.79902820 -0.49918890 -0.37578866
-0.71944669 -0.44936314
[16,] 1.214491e+00 -0.79898164 -0.79898164 -0.49918890 -0.37602149
-0.71921386 -0.44936314
[17,] 7.884577e-01 -0.49867667 -0.49867667 -0.36740676 -0.23562461
-0.45029446 -0.32992102
[18,] 5.675014e-01 -0.37644058 -0.37644058 -0.23585744 -0.25262125
-0.33900142 -0.20791776
[19,] 1.091324e+00 -0.71818940 -0.71818940 -0.45052730 -0.33900142
-0.64703636 -0.40209852
[20,] 7.117167e-01 -0.45038760 -0.45038760 -0.32759272 -0.21001324
-0.40465966 -0.29732473
[21,] 7.907860e-01 -0.49867667 -0.49867667 -0.36740676 -0.23562461
-0.45029446 -0.33038668
[,14] [,15] [,16] [,17]
[,18] [,19]
[1,] -0.18929131 2.813991e+00 2.814224e+00 1.213048e+00
1.317821e+00 2.537854e+00
[2,] -0.18929131 2.813991e+00 2.814224e+00 1.213048e+00
1.317821e+00 2.537854e+00
[3,] -0.12409873 1.214910e+00 1.212582e+00 7.892959e-01
5.704351e-01 1.089647e+00
[4,] -0.08800998 1.313863e+00 1.313863e+00 5.634502e-01
8.870848e-01 1.182780e+00
[5,] -0.17066486 2.533663e+00 2.533663e+00 1.094304e+00
1.185108e+00 2.281740e+00
[6,] -0.10989606 1.093373e+00 1.093373e+00 7.092021e-01
5.112961e-01 9.855721e-01
[7,] -0.12246892 1.214491e+00 1.214491e+00 7.884577e-01
5.675014e-01 1.091324e+00
[8,] 0.07711351 -7.990282e-01 -7.989816e-01 -4.986767e-01
-3.764406e-01 -7.181894e-01
[9,] 0.07711351 -7.990282e-01 -7.989816e-01 -4.986767e-01
-3.764406e-01 -7.181894e-01
[10,] 0.05681068 -4.991889e-01 -4.991889e-01 -3.674068e-01
-2.358574e-01 -4.505273e-01
[11,] 0.03585592 -3.757887e-01 -3.760215e-01 -2.356246e-01
-2.526212e-01 -3.390014e-01
[12,] 0.06938353 -7.194467e-01 -7.192139e-01 -4.502945e-01
-3.390014e-01 -6.470364e-01
[13,] 0.05122274 -4.493631e-01 -4.493631e-01 -3.299210e-01
-2.079178e-01 -4.020985e-01
[14,] 0.05681068 -4.991889e-01 -4.991889e-01 -3.674068e-01
-2.358574e-01 -4.505273e-01
[15,] -0.49918890 1.312200e+06 2.638980e+07 6.298560e+05
6.298560e+05 1.319490e+07
[16,] -0.49918890 2.638980e+07 7.437258e+08 1.266710e+07
1.266710e+07 3.718629e+08
[17,] -0.36740676 6.298560e+05 1.266710e+07 4.283021e+05
3.023309e+05 6.333552e+06
[18,] -0.23585744 6.298560e+05 1.266710e+07 3.023309e+05
4.283021e+05 6.333552e+06
[19,] -0.45052730 1.319490e+07 3.718629e+08 6.333552e+06
6.333552e+06 2.355132e+08
[20,] -0.32759272 3.149280e+05 6.333552e+06 2.141510e+05
1.511654e+05 4.011250e+06
[21,] -0.36740676 1.266710e+07 3.569884e+08 8.613631e+06
6.080210e+06 1.784942e+08
[,20] [,21]
[1,] 1.094304e+00 1.215376e+00
[2,] 1.094304e+00 1.215376e+00
[3,] 7.101335e-01 7.939525e-01
[4,] 5.122274e-01 5.657785e-01
[5,] 9.802170e-01 1.091976e+00
[6,] 6.388873e-01 7.094350e-01
[7,] 7.117167e-01 7.907860e-01
[8,] -4.503876e-01 -4.986767e-01
[9,] -4.503876e-01 -4.986767e-01
[10,] -3.275927e-01 -3.674068e-01
[11,] -2.100132e-01 -2.356246e-01
[12,] -4.046597e-01 -4.502945e-01
[13,] -2.973247e-01 -3.303867e-01
[14,] -3.275927e-01 -3.674068e-01
[15,] 3.149280e+05 1.266710e+07
[16,] 6.333552e+06 3.569884e+08
[17,] 2.141510e+05 8.613631e+06
[18,] 1.511654e+05 6.080210e+06
[19,] 4.011250e+06 1.784942e+08
[20,] 1.356290e+05 4.306815e+06
[21,] 4.306815e+06 2.427521e+08
What does the message about convergence mean? (I thought it was
something related to parameter scaling issues. That is why I provided a
'parscale' argument)
Also when trying to calculate SE, I got the following messages:
> OI<-solve(est_coeff$hessian)
Error in solve.default(est_coeff$hessian) :
Lapack routine dgesv: system is exactly singular
> se<-sqrt(diag(OI))
Warning message:
In sqrt(diag(OI)) : NaNs produced
>
I would appreciate it very much any suggestion you might give
Thanks for your help
Lorenzo
[[alternative HTML version deleted]]