search

swihart / gnlrim

Generalized Nonlinear Random Intercept Models in R aka a freestanding `repeated::gnlmix()`
0 stars 0 forks source link

make method argument for "L-BGFS-B" and "bobyqa" with passable specific options from parent gnlrim call as well as add options for nlminb #19

Open swihart opened 5 years ago

swihart commented 5 years ago

I'm working up these cases in scratch_testing.R. Here's how it looks before any edits:

sim_data_binomial <- readRDS("TSJYO_sim_data_binomial.RDS")
> attach(sim_data_binomial)
> y_cbind <- cbind(num_y1, num_y0)
> start.time<-Sys.time()
> gapr_marg_lbfgs <-
+   gnlrim(y=y_cbind,
+          mu=~ stable_cdf2((a0 + a1*group)*(1^alpha)/phi + rand, c(alpha, 0, 1, 0)),
+          ##pmu=c( a0=bm0_start, a1=bm1_start, alpha=1.5, phi=0.50),
+          pmu=c( a0=0, a1=0, alpha=1.5, phi=0.50),
+          pmix=c(alpha=1.5, phi=0.50),
+          p_uppb = c(  Inf ,  Inf, 2-1e-5,   1-1e-5),
+          p_lowb = c( -Inf, -Inf,  0+1e-5,  0+1e-5),
+          distribution="binomial",
+          nest=id,
+          random="rand",
+          mixture="libstableR-subgauss-phi",
+          ooo=TRUE,
+          compute_hessian = FALSE,
+          compute_kkt = FALSE,
+          trace=1,
+          method="L-BFGS-B"
+   )
fn is  fn 
Looking for method =  L-BFGS-B 
Function has  4  arguments
par[ 1 ]:  -Inf   <? 0   <? Inf     In Bounds  
par[ 2 ]:  -Inf   <? 0   <? Inf     In Bounds  
par[ 3 ]:  1e-05   <? 1.5   <? 1.99999     In Bounds  
par[ 4 ]:  1e-05   <? 0.5   <? 0.99999     In Bounds  
Analytic gradient not made available.
Analytic Hessian not made available.
Scale check -- log parameter ratio= 0.4771213   log bounds ratio= 0.3010343 
Method:  L-BFGS-B 
iter   10 value 44.261411
final  value 43.861168 
converged
Post processing for method  L-BFGS-B 
Successful convergence! 
Save results from method  L-BFGS-B 
$par
         a0          a1       alpha         phi 
  7.5716369 -12.7028823   1.1503491   0.4495364 

$value
[1] 43.86117

$message
[1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"

$convcode
[1] 0

$fevals
function 
      21 

$gevals
gradient 
      21 

$nitns
[1] NA

$kkt1
[1] NA

$kkt2
[1] NA

$xtimes
user.self 
   80.313 

Assemble the answers

Note the warning message for "L-BGFS-B":

Warning message:
In optim(par = par, fn = ufn, gr = ugr, lower = lower, upper = upper,  :
  unknown names in control: print.level, typsize, ndigit, gradtol, stepmax, steptol, fscale, step.max
> end.time <- Sys.time()
> time.taken <- end.time - start.time
> lbfgs_time.taken <- time.taken
> time.taken
Time difference of 41.9453 secs

Now start bobyqa -- same model same data...

> start.time<-Sys.time()
> gapr_marg_bobyqa <-
+   gnlrim(y=y_cbind,
+          mu=~ stable_cdf2((a0 + a1*group)*(1^alpha)/phi + rand, c(alpha, 0, 1, 0)),
+          ##pmu=c( a0=bm0_start, a1=bm1_start, alpha=1.5, phi=0.50),
+          pmu=c( a0=0, a1=0, alpha=1.5, phi=0.50),
+          pmix=c(alpha=1.5, phi=0.50),
+          p_uppb = c(  Inf ,  Inf, 2-1e-5,   1-1e-5),
+          p_lowb = c( -Inf, -Inf,  0+1e-5,  0+1e-5),
+          distribution="binomial",
+          nest=id,
+          random="rand",
+          mixture="libstableR-subgauss-phi",
+          ooo=TRUE,
+          compute_hessian = FALSE,
+          compute_kkt = FALSE,
+          trace=1,
+          method="bobyqa"
+   )
fn is  fn 
Looking for method =  bobyqa 
Function has  4  arguments
par[ 1 ]:  -Inf   <? 0   <? Inf     In Bounds  
par[ 2 ]:  -Inf   <? 0   <? Inf     In Bounds  
par[ 3 ]:  1e-05   <? 1.5   <? 1.99999     In Bounds  
par[ 4 ]:  1e-05   <? 0.5   <? 0.99999     In Bounds  
Analytic gradient not made available.
Analytic Hessian not made available.
Scale check -- log parameter ratio= 0.4771213   log bounds ratio= 0.3010343 
Method:  bobyqa 
start par. =  0 0 1.5 0.5 fn =  71.81986 
At return
eval: 100 fn:      56.366656 par: 0.243883 -0.300255  1.13286 0.0955449
Post processing for method  bobyqa 
Successful convergence! 
Save results from method  bobyqa 
parameter estimates: 0.243882748096423, -0.300254923731886, 1.1328556034127, 0.0955449037422799 
objective: 56.3666561292822 
number of function evaluations: 100 
Assemble the answers

Note the warning messages involving maxfun:

Warning messages:
1: In (function (npt = min(n + 2L, 2L * n), rhobeg = NA, rhoend = NA,  :
  unused control arguments ignored
2: In commonArgs(par, fn, control, environment()) :
  maxfun < 10 * length(par)^2 is not recommended.
> end.time <- Sys.time()
> time.taken <- end.time - start.time
> bobyqa_time.taken <- time.taken
> time.taken
Time difference of 17.0707 secs

bobyqa is so quick, but convereges to something weird.

start nlminb -- the default. are tolerances too tight?

> start.time<-Sys.time()
> gapr_marg_nlminb <-
+   gnlrim(y=y_cbind,
+          mu=~ stable_cdf2((a0 + a1*group)*(1^alpha)/phi + rand, c(alpha, 0, 1, 0)),
+          ##pmu=c( a0=bm0_start, a1=bm1_start, alpha=1.5, phi=0.50),
+          pmu=c( a0=0, a1=0, alpha=1.5, phi=0.50),
+          pmix=c(alpha=1.5, phi=0.50),
+          p_uppb = c(  Inf ,  Inf, 2-1e-5,   1-1e-5),
+          p_lowb = c( -Inf, -Inf,  0+1e-5,  0+1e-5),
+          distribution="binomial",
+          nest=id,
+          random="rand",
+          mixture="libstableR-subgauss-phi",
+          ooo=TRUE,
+          compute_hessian = FALSE,
+          compute_kkt = FALSE,
+          trace=1,
+          method="nlminb"
+   )
fn is  fn 
Looking for method =  nlminb 
Function has  4  arguments
par[ 1 ]:  -Inf   <? 0   <? Inf     In Bounds  
par[ 2 ]:  -Inf   <? 0   <? Inf     In Bounds  
par[ 3 ]:  1e-05   <? 1.5   <? 1.99999     In Bounds  
par[ 4 ]:  1e-05   <? 0.5   <? 0.99999     In Bounds  
Analytic gradient not made available.
Analytic Hessian not made available.
Scale check -- log parameter ratio= 0.4771213   log bounds ratio= 0.3010343 
Method:  nlminb 
  0:     71.819860:  0.00000  0.00000  1.50000 0.500000
  1:     66.496600: 0.00326712 -0.00469494  1.52102 0.420905
  2:     59.810823: 0.0124732 -0.0157064  1.56418 0.263259
  3:     58.541368: 0.147661 -0.145900  1.35375 0.0953913
  4:     56.504180: 0.349478 -0.391080  1.31257 0.167135
  5:     47.982940:  2.14036 -3.57406 0.972066 0.126653
  6:     47.227792:  2.60489 -4.19355  1.15703 0.149198
  7:     46.806567:  2.88977 -4.90150 0.934450 0.100404
  8:     46.714417:  2.89010 -4.90138 0.907603 0.103759
  9:     46.643799:  2.89426 -4.90027 0.908689 0.116539
 10:     45.757018:  4.74674 -8.18016 0.903677 0.124912
 11:     45.006723:  7.22065 -11.0193 0.905033 0.216867
 12:     44.045505:  7.78897 -12.6115 0.906680 0.249119
 13:     43.962659:  8.41807 -13.7239 0.906992 0.262863
 14:     43.905812:  9.20864 -14.9814 0.907918 0.282674
 15:     43.889650:  10.0922 -16.3553 0.909293 0.303567
 16:     43.869667:  10.9824 -17.7519 0.910845 0.321277
 17:     43.864211:  10.7292 -17.3406 0.910824 0.313492
 18:     43.860172:  10.4114 -16.8245 0.911102 0.302526
 19:     43.858984:  10.4086 -16.8215 0.911393 0.301259
 20:     43.848049:  10.4282 -16.8557 0.911845 0.296819
 21:     43.837266:  10.4281 -16.8558 0.911414 0.294799
 22:     43.836871:  10.4281 -16.8558 0.911403 0.294736
 23:     43.836070:  10.4281 -16.8558 0.911381 0.294611
 24:     43.835746:  10.4281 -16.8558 0.911373 0.294562
 25:     43.835091:  10.4281 -16.8559 0.911356 0.294462
 26:     43.834826:  10.4281 -16.8559 0.911349 0.294422
 27:     43.834293:  10.4281 -16.8559 0.911335 0.294342
 28:     43.834078:  10.4281 -16.8559 0.911330 0.294310
 29:     43.834035:  10.4281 -16.8559 0.911329 0.294303
 30:     43.833683:  10.4281 -16.8559 0.911320 0.294252
 31:     43.833669:  10.4281 -16.8559 0.911320 0.294250
 32:     43.833666:  10.4281 -16.8559 0.911320 0.294250
 33:     43.833665:  10.4281 -16.8559 0.911320 0.294249
 34:     43.833664:  10.4281 -16.8559 0.911320 0.294249
 35:     43.833664:  10.4281 -16.8559 0.911320 0.294249
 36:     43.833664:  10.4281 -16.8559 0.911320 0.294249
 37:     43.833664:  10.4281 -16.8559 0.911320 0.294249
 38:     43.833664:  10.4281 -16.8559 0.911320 0.294249
 39:     43.833663:  10.4281 -16.8559 0.911320 0.294249
 40:     43.833663:  10.4281 -16.8559 0.911320 0.294249
 41:     43.833663:  10.4281 -16.8559 0.911320 0.294249
 42:     43.830784:  10.4281 -16.8559 0.907361 0.294249
 43:     43.830437:  10.4281 -16.8559 0.907360 0.294200
 44:     43.830367:  10.4281 -16.8559 0.907360 0.294190
 45:     43.830228:  10.4281 -16.8559 0.907359 0.294170
 46:     43.830171:  10.4281 -16.8559 0.907358 0.294162
 47:     43.830058:  10.4281 -16.8559 0.907358 0.294147
 48:     43.830013:  10.4281 -16.8559 0.907357 0.294140
 49:     43.829923:  10.4281 -16.8559 0.907357 0.294128
 50:     43.829887:  10.4281 -16.8559 0.907356 0.294123
 51:     43.829880:  10.4281 -16.8559 0.907356 0.294122
 52:     43.829822:  10.4281 -16.8559 0.907356 0.294113
 53:     43.829822:  10.4281 -16.8559 0.907356 0.294113
 54:     43.829822:  10.4281 -16.8559 0.907356 0.294113
 55:     43.829822:  10.4281 -16.8559 0.907356 0.294113
 56:     43.829822:  10.4281 -16.8559 0.907356 0.294113
 57:     43.829822:  10.4281 -16.8559 0.907356 0.294113
 58:     43.829822:  10.4281 -16.8559 0.907356 0.294113
 59:     43.829822:  10.4281 -16.8559 0.907356 0.294113
 60:     43.829455:  10.4281 -16.8559 0.905390 0.294113
 61:     43.829271:  10.4281 -16.8559 0.905390 0.294089
 62:     43.828910:  10.4281 -16.8559 0.905390 0.294040
 63:     43.828908:  10.4281 -16.8559 0.905390 0.294040
 64:     43.828905:  10.4281 -16.8559 0.905390 0.294039
 65:     43.828904:  10.4281 -16.8559 0.905390 0.294039
 66:     43.828902:  10.4281 -16.8559 0.905390 0.294039
 67:     43.828901:  10.4281 -16.8559 0.905390 0.294039
 68:     43.828901:  10.4281 -16.8559 0.905390 0.294039
 69:     43.828900:  10.4281 -16.8559 0.905390 0.294038
 70:     43.828900:  10.4281 -16.8559 0.905390 0.294038
 71:     43.828900:  10.4281 -16.8559 0.905390 0.294038
 72:     43.828900:  10.4281 -16.8559 0.905390 0.294038
 73:     43.828900:  10.4281 -16.8559 0.905390 0.294038
 74:     43.828900:  10.4281 -16.8559 0.905390 0.294038
 75:     43.828900:  10.4281 -16.8559 0.905390 0.294038
Post processing for method  nlminb 
Save results from method  nlminb 
$par
         a0          a1       alpha         phi 
 10.4280946 -16.8558713   0.9053898   0.2940384 

$message
[1] "function evaluation limit reached without convergence (9)"

$convcode
[1] 1

$value
[1] 43.8289

$fevals
function 
     200 

$gevals
gradient 
     433 

$nitns
[1] 75

$kkt1
[1] NA

$kkt2
[1] NA

$xtimes
user.self 
   866.51 

Assemble the answers

Note the warnings:

Warning messages:
1: In nlminb(start = par, objective = ufn, gradient = ugr, lower = lower,  :
  NA/NaN function evaluation
2: In nlminb(start = par, objective = ufn, gradient = ugr, lower = lower,  :
  NA/NaN function evaluation
> end.time <- Sys.time()
> time.taken <- end.time - start.time
> nlminb_time.taken <- time.taken
> time.taken
Time difference of 7.84591 mins

Slow and steady. Now, the nlminb marg model matches the conditional model in the TSJYO summary table, and the conditional model converged. This makes me think it might be a false alarm about not converging...

> rbind(gapr_marg_nlminb,gapr_marg_bobyqa,gapr_marg_lbfgs)
                 a0          a1     alpha       phi    value fevals gevals niter
nlminb   10.4280946 -16.8558713 0.9053898 0.2940384 43.82890    200    433    75
bobyqa    0.2438827  -0.3002549 1.1328556 0.0955449 56.36666    100     NA    NA
L-BFGS-B  7.5716369 -12.7028823 1.1503491 0.4495364 43.86117     21     21    NA
         convcode kkt1 kkt2   xtime
nlminb          1   NA   NA 866.510
bobyqa          0   NA   NA  31.702
L-BFGS-B        0   NA   NA  80.313

Compare the L-BFGS-B above to the conditional model below. Looks like a winner!

> readRDS("TSJYO_summary_table1.RDS")[5:7,]
     Approach ConvCode   LogLik      Bc0       Bc1     alpha    gam_g       phi
5 GAPR (Cond)        0 43.86117 16.84303 -28.25718 1.1503506 1.429786 0.4495233
6 GAPR (Marg)        1 43.83112 35.44501 -57.29514 0.9090899 2.190065 0.2943947
7       Truth       NA       NA 11.00000 -22.00000 1.2000000 2.000000 0.3699906
        Bm0       Bm1 CondProb0.00 CondProb0.33 CondProb0.67 CondProb1.00 MargProb0.00
5  7.571332 -12.70226    0.9886690     0.970962   0.12671626   0.01782203    0.9712008
6 10.434826 -16.86739    0.9870186     0.974194   0.11612771   0.02008772    0.9610480
7  4.000000  -8.00000    0.9840123     0.939380   0.06061995   0.01598768    0.9442213
  MargProb0.33 MargProb0.67 MargProb1.00 Ji N
5    0.9260358    0.2583170   0.04545054 60 5
6    0.9239563    0.2734200   0.05971113 60 5
7    0.8090422    0.1909578   0.05577873 60 5
swihart commented 5 years ago
swihart commented 5 years ago

Progress! But the estimates are still way off for nlminb:

> gapr_marg_nlminb <-
+   gnlrim(y=y_cbind,
+          mu=~ stable_cdf2((a0 + a1*group)*(1^alpha)/phi + rand, c(alpha, 0, 1, 0)),
+          ##pmu=c( a0=bm0_start, a1=bm1_start, alpha=1.5, phi=0.50),
+          pmu=c( a0=0, a1=0, alpha=1.5, phi=0.50),
+          pmix=c(alpha=1.5, phi=0.50),
+          p_uppb = c(  Inf ,  Inf, 2-1e-5,   1-1e-5),
+          p_lowb = c( -Inf, -Inf,  0+1e-5,  0+1e-5),
+          distribution="binomial",
+          nest=id,
+          random="rand",
+          mixture="libstableR-subgauss-phi",
+          ooo=TRUE,
+          compute_hessian = FALSE,
+          compute_kkt = FALSE,
+          trace=1,
+          method="nlminb",
+          abs.tol.nlminb=1e-5, ## 1e-20,
+          xf.tol.nlminb=2.2e-6, ##2.2e-14,
+          x.tol.nlminb=1.5e-3, ##1.5e-8,
+          rel.tol.nlminb=1e-4, ##1e-10,
+   )
fn is  fn 
Looking for method =  nlminb 
Function has  4  arguments
par[ 1 ]:  -Inf   <? 0   <? Inf     In Bounds  
par[ 2 ]:  -Inf   <? 0   <? Inf     In Bounds  
par[ 3 ]:  1e-05   <? 1.5   <? 1.99999     In Bounds  
par[ 4 ]:  1e-05   <? 0.5   <? 0.99999     In Bounds  
Analytic gradient not made available.
Analytic Hessian not made available.
Scale check -- log parameter ratio= 0.4771213   log bounds ratio= 0.3010343 
Method:  nlminb 
  0:     71.819860:  0.00000  0.00000  1.50000 0.500000
  1:     66.496600: 0.00326712 -0.00469494  1.52102 0.420905
  2:     59.810823: 0.0124732 -0.0157064  1.56418 0.263259
  3:     58.541368: 0.147661 -0.145900  1.35375 0.0953913
  4:     56.504180: 0.349478 -0.391080  1.31257 0.167135
  5:     47.982940:  2.14036 -3.57406 0.972066 0.126653
  6:     47.227792:  2.60489 -4.19355  1.15703 0.149198
  7:     46.806567:  2.88977 -4.90150 0.934450 0.100404
  8:     46.714417:  2.89010 -4.90138 0.907603 0.103759
  9:     46.643799:  2.89426 -4.90027 0.908689 0.116539
 10:     45.757018:  4.74674 -8.18016 0.903677 0.124912
 11:     45.006723:  7.22065 -11.0193 0.905033 0.216867
 12:     44.045505:  7.78897 -12.6115 0.906680 0.249119
 13:     43.962659:  8.41807 -13.7239 0.906992 0.262863
 14:     43.905812:  9.20864 -14.9814 0.907918 0.282674
 15:     43.889650:  10.0922 -16.3553 0.909293 0.303567
 16:     43.869667:  10.9824 -17.7519 0.910845 0.321277
 17:     43.864211:  10.7292 -17.3406 0.910824 0.313492
Post processing for method  nlminb 
Successful convergence! 
Save results from method  nlminb 
$par
         a0          a1       alpha         phi 
 10.7292082 -17.3405664   0.9108244   0.3134918 

$message
[1] "relative convergence (4)"

$convcode
[1] 0

$value
[1] 43.86421

$fevals
function 
      31 

$gevals
gradient 
      72 

$nitns
[1] 17

$kkt1
[1] NA

$kkt2
[1] NA

$xtimes
user.self 
   91.086
swihart commented 5 years ago

goal! But had to put lower bound on alpha to do it. Interesting.

> gapr_marg_nlminb2 <-
+   gnlrim(y=y_cbind,
+          mu=~ stable_cdf2((a0 + a1*group)*(1^alpha)/phi + rand, c(alpha, 0, 1, 0)),
+          ##pmu=c( a0=bm0_start, a1=bm1_start, alpha=1.5, phi=0.50),
+          pmu=c( a0=0, a1=0, alpha=1.5, phi=0.50),
+          pmix=c(alpha=1.5, phi=0.50),
+          p_uppb = c(  Inf ,  Inf, 2-1e-5,   1-1e-5),
+          p_lowb = c( -Inf, -Inf,  1.15,  0+1e-5),
+          distribution="binomial",
+          nest=id,
+          random="rand",
+          mixture="libstableR-subgauss-phi",
+          ooo=TRUE,
+          compute_hessian = FALSE,
+          compute_kkt = FALSE,
+          trace=1,
+          method="nlminb",
+          abs.tol.nlminb=1e-20, ## 1e-20,
+          xf.tol.nlminb=2.2e-10, ##2.2e-14,
+          x.tol.nlminb=1.5e-8, ##1.5e-8,
+          rel.tol.nlminb=1e-10, ##1e-10,
+   )
fn is  fn 
Looking for method =  nlminb 
Function has  4  arguments
par[ 1 ]:  -Inf   <? 0   <? Inf     In Bounds  
par[ 2 ]:  -Inf   <? 0   <? Inf     In Bounds  
par[ 3 ]:  1.15   <? 1.5   <? 1.99999     In Bounds  
par[ 4 ]:  1e-05   <? 0.5   <? 0.99999     In Bounds  
Analytic gradient not made available.
Analytic Hessian not made available.
Scale check -- log parameter ratio= 0.4771213   log bounds ratio= 0.0705775 
Method:  nlminb 
  0:     71.819860:  0.00000  0.00000  1.50000 0.500000
  1:     66.496600: 0.00326712 -0.00469494  1.52102 0.420905
  2:     59.810823: 0.0124732 -0.0157064  1.56418 0.263259
  3:     58.541368: 0.147661 -0.145900  1.35375 0.0953913
  4:     56.504180: 0.349478 -0.391080  1.31257 0.167135
  5:     46.515541:  2.84655 -4.73472  1.15000 0.174623
  6:     46.198757:  3.38396 -5.19292  1.15000 0.187130
  7:     44.753816:  4.55416 -7.76376  1.15000 0.250539
  8:     44.038791:  6.07374 -10.1432  1.15000 0.358069
  9:     43.916601:  6.70881 -11.2583  1.15000 0.397498
 10:     43.864832:  7.40337 -12.4510  1.15000 0.439471
 11:     43.862070:  7.55060 -12.6867  1.15000 0.447611
 12:     43.861449:  7.59436 -12.7462  1.15000 0.449680
 13:     43.861295:  7.59347 -12.7370  1.15000 0.449507
 14:     43.861240:  7.58402 -12.7191  1.15000 0.449129
 15:     43.861193:  7.57288 -12.7007  1.15000 0.448908
 16:     43.861174:  7.56972 -12.6973  1.15000 0.449090
 17:     43.861173:  7.56813 -12.6949  1.15021 0.449208
 18:     43.861171:  7.56942 -12.6975  1.15016 0.449278
 19:     43.861171:  7.56940 -12.6975  1.15019 0.449238
 20:     43.861170:  7.56931 -12.6975  1.15020 0.449298
 21:     43.861170:  7.56935 -12.6977  1.15023 0.449288
 22:     43.861168:  7.56897 -12.6980  1.15033 0.449343
 23:     43.861168:  7.57026 -12.7004  1.15033 0.449439
 24:     43.861168:  7.57128 -12.7022  1.15035 0.449521
 25:     43.861168:  7.57133 -12.7023  1.15035 0.449523
Post processing for method  nlminb 
Successful convergence! 
Save results from method  nlminb 
$par
         a0          a1       alpha         phi 
  7.5713326 -12.7022632   1.1503501   0.4495231 

$message
[1] "relative convergence (4)"

$convcode
[1] 0

$value
[1] 43.86117

$fevals
function 
      39 

$gevals
gradient 
     146 

$nitns
[1] 25

$kkt1
[1] NA

$kkt2
[1] NA

$xtimes
user.self 
   90.121 

Assemble the answers