Generalize from one to two random effects, possibly correlated

swihart commented 2 years ago

`gnlrem()`: that is, `gnlrim()` that allows for two random effects, possibly correlated

the `i` changes to an `e`, but the `i` wasn't necessarily just an intercept since you had the freedom to put the random parameter anywhere. Just an aside.

✅ Create branch two_rand_effs ✅ Get int and int2 into gnlrim repo; Rd2roxygen; export, etc. rename to romberg_int1d() and romberg_int2d(). ✅ copy gnlrim and rename gnlrem and get it roxygen'd. 🔲 --~--> SHOULD I rename this gnlmer2rp ?

✅ Figure out how rand is handled; consider number of random parameters...

✅ Work up the lmer example (see gmail) for random intercept, slopes, etc. so we have a target example

So it appears I have correlated random intercept and slopes bivariate normal case implemented. See below for some tests which compare numbers with glmer and GLMMadaptive. In the spirit of the README where I said I wanted a data argument and random slopes, I'm reminded with my work with Rolf and will probably do an on.exit with attach approach using a LPN() wrapper of gnlrem. This way, the ugly internals (and the name gnlrem vs gnlmer2rp) don't matter as much and all the different arguments/customizations are prepackaged for an easier user experience and I can have a "pseudo" data argument that's a touch more responsible than plain attach() -- I fear getting a true data argument in there will mess up the finterp internals of defining functions.

swihart commented 2 years ago

This helped with the integrals: https://cran.r-project.org/web/packages/Rd2roxygen/vignettes/Rd2roxygen.html

install.packages("Rd2roxygen")
library(Rd2roxygen)
options(roxygen.comment = "##' ")
str(info <- parse_file("../rmutil/man/int2.Rd"))
cat(create_roxygen(info), sep = '\n')

Produces the following, which I copy and paste into romberg_int2d

##' Vectorized Two-dimensional Numerical Integration
##' 
##' \code{int} performs vectorized numerical integration of a given
##' two-dimensional function.
##' 
##' 
##' @param f The function (of two variables) to integrate, returning either a
##' scalar or a vector.
##' @param a A two-element vector or a two-column matrix giving the lower
##' bounds. It cannot contain both -Inf and finite values.
##' @param b A two-element vector or a two-column matrix giving the upper
##' bounds. It cannot contain both Inf and finite values.
##' @param eps Precision.
##' @param max The maximum number of steps, by default set to 16.
##' @param d The number of extrapolation points so that 2k is the order of
##' integration, by default set to 5; d=2 is Simpson's rule.
##' @return The vector of values of the integrals of the function supplied.
##' @author J.K. Lindsey
##' @keywords math
##' @examples
##' 
##' f <- function(x,y) sin(x)+cos(y)-x^2
##' int2(f, a=c(0,1), b=c(2,4))
##' #
##' fn1 <- function(x, y) x^2+y^2
##' fn2 <- function(x, y) (1:4)*x^2+(2:5)*y^2
##' int2(fn1, c(1,2), c(2,4))
##' int2(fn2, c(1,2), c(2,4))
##' int2(fn1, matrix(c(1:4,1:4),ncol=2), matrix(c(2:5,2:5),ncol=2))
##' int2(fn2, matrix(c(1:4,1:4),ncol=2), matrix(c(2:5,2:5),ncol=2))
##'

swihart commented 2 years ago

Made some headway on using pcubature -- looks faster and more capable of getting an answer...(?)

R version 4.1.2 (2021-11-01) -- "Bird Hippie"
Copyright (C) 2021 The R Foundation for Statistical Computing
Platform: x86_64-apple-darwin17.0 (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

Restarting R session...

> library(gnlrim)
> ########################################$$$$#########
> ## 2022-02-10                                      ##
> ## START two random parameters                     ##
> #####################################################
> # Derived expression for gamma
> g <- function(a) {
+   iu <- complex(real=0, imaginary=1)
+   return(abs(1 - iu * tan(pi * a / 2)) ^ (-1 / a))
+ }
> 
> bridgecloglog_rstable <- function(n, alpha, delta){
+   mult <- (delta/alpha)^(1/alpha)
+   X <- stabledist::rstable(n , alpha, beta=1, gamma=g(alpha), delta=0, pm=1)
+   Z <- log(mult * X)
+   Z
+ }
> 
> ## add in beta random effect
> sim_mrim_data <- function(n1, n2, J, a0, a1, v1=1,v2=2,rho=0.5, mrim="Logistic-BIVARIATE_NORM", alpha=1.89, gamma=1.2, delta=1){
+ 
+   if(mrim=="Logistic-BIVARIATE_NORM"){
+     print("HERE")
+     G <- function(n, v1, v2, rho){mvtnorm::rmvnorm(n=n, sigma=matrix(c(v1,rho*sqrt(v1*v2),rho*sqrt(v1*v2),v2),nrow=2))}
+     H <- function(x) plogis(x)
+   }
+ 
+   n <- n1 + n2
+   u <- round(apply(G(n,v1,v2,rho),2, function(W) rep(W,each=J)),2)
+ 
+   ## x <- c(rep(1, n1*J), rep(0, n2*J))
+   ##  x <-c(runif(n1*J, 0.5,1.5), runif(n2*J, 0,0.60))
+   x <-c(sample(c(1,2,3,4)/10,n1*J,TRUE), sample(c(1.2,2.2,3.2,4.2)/10,n2*J,TRUE))
+   eta <- round(a0 + a1*x,2)
+ 
+   eta_i <- round( (a0 + u[,1]) + (a1+u[,2])*x, 2)
+   py1 <- round(H(eta_i),2)
+   y <- rbinom(length(eta_i), 1, prob=py1 )
+ 
+   data.frame(id=rep(1:n, each=J),
+              j = rep(1:J),
+              x1 = x,
+              eta = eta,
+              u_i = u,
+              eta_i = eta_i,
+              py1 = py1,
+              y=y
+   )
+ 
+ }
> 
> 
> detach(summed_binom_dat)
Error in detach(summed_binom_dat) : invalid 'name' argument
> set.seed(6)
> binom_dat <-
+ #  sim_mrim_data(800,400, J=100, a0 = -2, a1 = 1)
+ #  sim_mrim_data(4000,2000, J=100, a0 = -2, a1 = 1)
+   sim_mrim_data(8,8, J=100, a0 = -2, a1 = 1)
[1] "HERE"
> data.table::setDT(binom_dat)
> 
> summed_binom_dat <-
+   binom_dat[, {j=list(r=sum(y), n_r=sum(y==0))}, by=c("id","x1")]
> data.table::setkey(summed_binom_dat,id, x1)
> summed_binom_dat
    id   x1  r n_r
 1:  1 0.10  3  17
 2:  1 0.20  3  19
 3:  1 0.30  4  25
 4:  1 0.40  5  24
 5:  2 0.10 15  16
 6:  2 0.20 18  11
 7:  2 0.30 11   9
 8:  2 0.40 13   7
 9:  3 0.10  7  24
10:  3 0.20  7  17
11:  3 0.30  5  15
12:  3 0.40  6  19
13:  4 0.10  1  20
14:  4 0.20  0  28
15:  4 0.30  2  20
16:  4 0.40  4  25
17:  5 0.10  2  16
18:  5 0.20  2  27
19:  5 0.30  4  26
20:  5 0.40  1  22
21:  6 0.10  8  20
22:  6 0.20  9  15
23:  6 0.30 10  17
24:  6 0.40  8  13
25:  7 0.10  6  16
26:  7 0.20  5  23
27:  7 0.30  2  16
28:  7 0.40  6  26
29:  8 0.10  1  17
30:  8 0.20  1  23
31:  8 0.30  1  25
32:  8 0.40  3  29
33:  9 0.12 13  17
34:  9 0.22  5  18
35:  9 0.32  8  19
36:  9 0.42  5  15
37: 10 0.12  8  16
38: 10 0.22  9  15
39: 10 0.32 11  14
40: 10 0.42 16  11
41: 11 0.12  2  26
42: 11 0.22  2  26
43: 11 0.32  2  15
44: 11 0.42  2  25
45: 12 0.12 14  13
46: 12 0.22  8  10
47: 12 0.32  9  17
48: 12 0.42 13  16
49: 13 0.12  2  26
50: 13 0.22  1  16
51: 13 0.32  2  26
52: 13 0.42  2  25
53: 14 0.12  5  16
54: 14 0.22  6  18
55: 14 0.32 14  14
56: 14 0.42 19   8
57: 15 0.12  5  18
58: 15 0.22  6  20
59: 15 0.32  4  15
60: 15 0.42 16  16
61: 16 0.12  4  21
62: 16 0.22  2  26
63: 16 0.32  3  22
64: 16 0.42  6  16
    id   x1  r n_r
> ## glmer with correlation between random intercept and random slope
> lme4::glmer(cbind(r, n_r) ~ x1 + (x1 | id), summed_binom_dat, binomial, nAGQ = 0)
Generalized linear mixed model fit by maximum likelihood
  (Adaptive Gauss-Hermite Quadrature, nAGQ = 0) [glmerMod]
 Family: binomial  ( logit )
Formula: cbind(r, n_r) ~ x1 + (x1 | id)
   Data: summed_binom_dat
      AIC       BIC    logLik  deviance  df.resid 
 302.5655  313.3599 -146.2828  292.5655        59 
Random effects:
 Groups Name        Std.Dev. Corr 
 id     (Intercept) 0.9496        
        x1          1.8024   -0.28
Number of obs: 64, groups:  id, 16
Fixed Effects:
(Intercept)           x1  
     -1.664        1.362  
> lme4::glmer(cbind(r, n_r) ~ x1 + (x1 | id), summed_binom_dat, binomial, nAGQ = 1)
Generalized linear mixed model fit by maximum likelihood
  (Laplace Approximation) [glmerMod]
 Family: binomial  ( logit )
Formula: cbind(r, n_r) ~ x1 + (x1 | id)
   Data: summed_binom_dat
      AIC       BIC    logLik  deviance  df.resid 
 302.5518  313.3463 -146.2759  292.5518        59 
Random effects:
 Groups Name        Std.Dev. Corr 
 id     (Intercept) 0.9495        
        x1          1.8039   -0.28
Number of obs: 64, groups:  id, 16
Fixed Effects:
(Intercept)           x1  
     -1.695        1.379  
> 
> 
> 
> attach(summed_binom_dat)
> 
> ybind <- cbind(r,n_r)
> period_numeric <- x1
> (rand.int.rand.slopes.nonzero.corr.CUBA <-
+     gnlrem(y=ybind,
+            mu = ~ plogis(Intercept + period_numeric*b_p + rand1 + rand2*b_p),
+            pmu = c(Intercept=-0.95, b_p=0.55),
+            pmix=c(var1=1, var2=1, corr12= 0.20),
+            p_uppb = c(  0,   2, 4.00, 4.00, 0.90),
+            p_lowb = c( -4,  -2, 0.05, 0.05,-0.90),
+            distribution="binomial",
+            nest=id,
+            random=c("rand1", "rand2"),
+            mixture="bivariate-normal-corr",
+            ooo=TRUE,
+            compute_hessian = FALSE,
+            compute_kkt = FALSE,
+            trace=1,
+            method='nlminb',
+            int2dmethod="cuba"
+     )
+ )
[1] "HERE-boy-0"
[1] 5
Intercept       b_p      var1      var2    corr12 
    -0.95      0.55      1.00      1.00      0.20 
[1] 151.7783
fn is  fn 
Looking for method =  nlminb 
Function has  5  arguments
par[ 1 ]:  -4   <? -0.95   <? 0     In Bounds  
par[ 2 ]:  -2   <? 0.55   <? 2     In Bounds  
par[ 3 ]:  0.05   <? 1   <? 4     In Bounds  
par[ 4 ]:  0.05   <? 1   <? 4     In Bounds  
par[ 5 ]:  -0.9   <? 0.2   <? 0.9     In Bounds  
[1] "HERE-boy-0"
Analytic gradient not made available.
Analytic Hessian not made available.
Scale check -- log parameter ratio= 0.69897   log bounds ratio= 0.3467875 
Method:  nlminb 
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
  0:     151.77828: -0.950000 0.550000  1.00000  1.00000 0.200000
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
  1:     150.85955: -1.88940 0.552229 0.764557 0.912504 -0.0333218
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
  2:     149.96193: -1.79929  1.51621 0.596072 0.842885 -0.204764
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
  3:     148.99004: -1.79163  1.43275 0.571740 0.764282 -0.287673
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
  4:     148.13210: -1.75699  1.28756 0.522186 0.603705 -0.466986
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
  5:     148.05020: -1.61151  1.50475 0.456672 0.538977 -0.543503
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
  6:     147.98670: -1.77023  1.61087 0.462606 0.511665 -0.531033
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
  7:     147.96894: -1.76072  1.58829 0.461122 0.485315 -0.554216
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
  8:     147.95548: -1.74020  1.56247 0.457002 0.507764 -0.539204
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
  9:     147.95298: -1.72509  1.52523 0.455023 0.520286 -0.531544
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
 10:     147.95247: -1.72365  1.52188 0.459325 0.521224 -0.523196
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
 11:     147.95247: -1.72293  1.52214 0.458531 0.521537 -0.523911
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
 12:     147.95247: -1.72304  1.52201 0.458659 0.521499 -0.523784
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
[1] "HERE-boy-0"
 13:     147.95247: -1.72304  1.52202 0.458659 0.521498 -0.523784
[1] "HERE-boy-0"
Post processing for method  nlminb 
Successful convergence! 
Save results from method  nlminb 
$par
 Intercept        b_p       var1       var2     corr12 
-1.7230417  1.5220164  0.4586590  0.5214982 -0.5237844 

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

$convcode
[1] 0

$value
[1] 147.9525

$fevals
function 
      16 

$gevals
gradient 
      86 

$nitns
[1] 13

$kkt1
[1] NA

$kkt2
[1] NA

$xtimes
user.self 
  154.771 

Assemble the answers
       Intercept      b_p     var1      var2     corr12    value
nlminb -1.723042 1.522016 0.458659 0.5214982 -0.5237844 147.9525
       fevals gevals niter convcode kkt1 kkt2   xtime
nlminb     16     86    13        0   NA   NA 154.771
> (rand.int.rand.slopes.nonzero.corr <-
+     gnlrem(y=ybind,
+            mu = ~ plogis(Intercept + period_numeric*b_p + rand1 + rand2*b_p),
+            pmu = c(Intercept=-0.95, b_p=0.55),
+            pmix=c(var1=1, var2=1, corr12= 0.20),
+            p_uppb = c(  0,   2, 4.00, 4.00, 0.90),
+            p_lowb = c( -4,  -2, 0.05, 0.05,-0.90),
+            distribution="binomial",
+            nest=id,
+            random=c("rand1", "rand2"),
+            mixture="bivariate-normal-corr",
+            ooo=TRUE,
+            compute_hessian = FALSE,
+            compute_kkt = FALSE,
+            trace=1,
+            method='nlminb'
+     )
+ )
[1] 5
Intercept       b_p      var1      var2    corr12 
    -0.95      0.55      1.00      1.00      0.20 
[1] 151.7783
fn is  fn 
Looking for method =  nlminb 
Function has  5  arguments
par[ 1 ]:  -4   <? -0.95   <? 0     In Bounds  
par[ 2 ]:  -2   <? 0.55   <? 2     In Bounds  
par[ 3 ]:  0.05   <? 1   <? 4     In Bounds  
par[ 4 ]:  0.05   <? 1   <? 4     In Bounds  
par[ 5 ]:  -0.9   <? 0.2   <? 0.9     In Bounds  
Analytic gradient not made available.
Analytic Hessian not made available.
Scale check -- log parameter ratio= 0.69897   log bounds ratio= 0.3467875 
Method:  nlminb 
  0:     151.77828: -0.950000 0.550000  1.00000  1.00000 0.200000
  1:     150.85953: -1.88940 0.552230 0.764558 0.912503 -0.0333215
  2:     149.96190: -1.79929  1.51621 0.596074 0.842880 -0.204767
  3:     148.99006: -1.79163  1.43275 0.571744 0.764281 -0.287671
  4:     148.13207: -1.75699  1.28757 0.522194 0.603713 -0.466973
  5:     148.05018: -1.61151  1.50474 0.456693 0.538992 -0.543493
  6:     147.98662: -1.77024  1.61092 0.462663 0.511586 -0.531111
  7:     147.96897: -1.76071  1.58835 0.461181 0.485267 -0.554269
nlminb function evaluation failure
Post processing for method  nlminb 
Save results from method  nlminb 
$fevals
[1] NA

$convcode
[1] 9999

$value
[1] 8.988466e+307

$par
[1] NA NA NA NA NA

$nitns
[1] NA

$gevals
[1] NA

$message
[1] "nlminb failure"

$kkt1
[1] NA

$kkt2
[1] NA

$xtimes
user.self 
 2251.421 

Assemble the answers
       Intercept b_p var1 var2 corr12         value fevals gevals
nlminb        NA  NA   NA   NA     NA 8.988466e+307     NA     NA
       niter convcode kkt1 kkt2    xtime
nlminb    NA     9999   NA   NA 2251.421

swihart commented 2 years ago

You made a simple mistake.

Which led you to discovering pcubature fdim which sped things up 16x. and also about GLMMadaptive and more about the nuts and bolts of how estimating glmms.

But it all came from a silly mistake. Observe a fit from a random intercept model:

Assemble the answers
       Intercept      b_p      var1    value fevals gevals niter convcode kkt1 kkt2 xtime
nlminb -1.723029 1.522013 0.8869422 147.9525     14     44    10        0   NA   NA  8.47

And when you were fitting random intercepts and slopes:

Assemble the answers
       Intercept      b_p     var1      var2     corr12    value fevals gevals niter convcode
nlminb -1.723042 1.522016 0.458659 0.5214982 -0.5237844 147.9525     16     86    13        0

Same likelihood? Wonky correlation estimates? The exact same likelihood made you think about how maybe the two random effects were combining into one. So you did this calculation on the random effects matrix:

> matrix(c(1,1.522013),nrow=1) %*%matrix(c(.458659,-0.5237844*sqrt(.458659*0.5214982) ,-0.5237844*sqrt(.458659*0.5214982), 0.5214982  ),ncol=2) %*% matrix(c(1,1.522013),ncol=1)
          [,1]
[1,] 0.8869415

So, when you linearly combine the random effects you get the random-intercept only variance ... and the non-one weight is the same as the fixed effect slope? You've just punked yourself! Look at the call:

(rand.int.rand.slopes.nonzero.corr.CUBA <-
    gnlrem(y=ybind,
           mu = ~ plogis(Intercept + period_numeric*b_p + rand1 + rand2*b_p),
           pmu = c(Intercept=-0.95, b_p=0.55),
           pmix=c(var1=1, var2=1, corr12= 0.20),
           p_uppb = c(  0,   2, 4.00, 4.00, 0.90),
           p_lowb = c( -4,  -2, 0.05, 0.05,-0.90),
           distribution="binomial",
           nest=id,
           random=c("rand1", "rand2"),
           mixture="bivariate-normal-corr",
           ooo=TRUE,
           compute_hessian = FALSE,
           compute_kkt = FALSE,
           trace=1,
           method='nlminb',
           int2dmethod="cuba"
    )
)

specifically

           mu = ~ plogis(Intercept + period_numeric*b_p + rand1 + rand2*b_p),

SHOULD BE

        mu = ~ plogis(Intercept + period_numeric*b_p + rand1 + period_numeric*rand2),

b_p is the fixed effect slope parameter. X1 is the DATA!!!!!

So now we run:

(rand.int.rand.slopes.nonzero.corr.CUBA <-
    gnlrem(y=ybind,
           mu = ~ plogis(Intercept + period_numeric*b_p + rand1 + period_numeric*rand2),
           pmu = c(Intercept=-0.95, b_p=0.55),
           pmix=c(var1=1, var2=1, corr12= 0.20),
           p_uppb = c(  0,   2, 4.00, 4.00, 0.90),
           p_lowb = c( -4,  -2, 0.05, 0.05,-0.90),
           distribution="binomial",
           nest=id,
           random=c("rand1", "rand2"),
           mixture="bivariate-normal-corr",
           ooo=TRUE,
           compute_hessian = FALSE,
           compute_kkt = FALSE,
           trace=1,
           method='nlminb',
           int2dmethod="cuba"
    )
)

and get

Assemble the answers
       Intercept      b_p      var1     var2     corr12    value 
nlminb -1.695897 1.379829 0.9103922 3.326855 -0.2841794 146.1776

> sqrt(0.9103922)
[1] 0.9541447
> sqrt(3.326855)
[1] 1.823967

which is more similar to:

> lme4::glmer(cbind(r, n_r) ~ x1 + (x1 | id), summed_binom_dat, binomial, nAGQ = 0)
Generalized linear mixed model fit by maximum likelihood (Adaptive Gauss-Hermite Quadrature, nAGQ = 0) [
glmerMod]
 Family: binomial  ( logit )
Formula: cbind(r, n_r) ~ x1 + (x1 | id)
   Data: summed_binom_dat
      AIC       BIC    logLik  deviance  df.resid 
 302.5655  313.3599 -146.2828  292.5655        59 
Random effects:
 Groups Name        Std.Dev. Corr 
 id     (Intercept) 0.9496        
        x1          1.8024   -0.28
Number of obs: 64, groups:  id, 16
Fixed Effects:
(Intercept)           x1  
     -1.664        1.362  
> lme4::glmer(cbind(r, n_r) ~ x1 + (x1 | id), summed_binom_dat, binomial, nAGQ = 1)
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: binomial  ( logit )
Formula: cbind(r, n_r) ~ x1 + (x1 | id)
   Data: summed_binom_dat
      AIC       BIC    logLik  deviance  df.resid 
 302.5518  313.3463 -146.2759  292.5518        59 
Random effects:
 Groups Name        Std.Dev. Corr 
 id     (Intercept) 0.9495        
        x1          1.8039   -0.28
Number of obs: 64, groups:  id, 16
Fixed Effects:
(Intercept)           x1  
     -1.695        1.379

swihart commented 2 years ago

Here's an example of uncorrelated/correlated and comparing to the glmer example (where I modified data to get a numeric variable and fit a slope/random slope just for didactic-proof-of-concept purposes):

########################################$$$$#########
## 2022-02-09                                      ##
## START two random parameters                     ##
#####################################################
library(lme4)
# from ?glmer
cbpp$period2 <- (cbpp$period == 2) + 0L
cbpp$period3 <- (cbpp$period == 3) + 0L
cbpp$period4 <- (cbpp$period == 4) + 0L
cbpp$period_numeric <- as.numeric(cbpp$period)

attach(cbpp)

## glmer with correlation between random intercept and random slope
glmer(cbind(incidence, size - incidence) ~ period_numeric + (period_numeric | herd), cbpp, binomial, nAGQ = 0)
# Generalized linear mixed model fit by maximum likelihood (Adaptive Gauss-Hermite
#                                                           Quadrature, nAGQ = 0) [glmerMod]
# Family: binomial  ( logit )
# Formula: cbind(incidence, size - incidence) ~ period_numeric + (period_numeric |      herd)
# Data: cbpp
# AIC      BIC   logLik deviance df.resid
# 194.3058 204.4325 -92.1529 184.3058       51
# Random effects:
#   Groups Name           Std.Dev. Corr
# herd   (Intercept)      1.0603
# period_numeric          0.3311   -0.86
# Number of obs: 56, groups:  herd, 15
# Fixed Effects:
#   (Intercept)  period_numeric
# -0.8765         -0.5843
##
##
#
#
## glmer wout correlation between random intercept and random slope
glmer(cbind(incidence, size - incidence) ~ period_numeric + (1 | herd) + (0+period_numeric|herd), cbpp, binomial, nAGQ = 0)
# Generalized linear mixed model fit by maximum likelihood (Adaptive Gauss-Hermite
#                                                           Quadrature, nAGQ = 0) [glmerMod]
# Family: binomial  ( logit )
# Formula: cbind(incidence, size - incidence) ~ period_numeric + (1 | herd) +
#   (0 + period_numeric | herd)
# Data: cbpp
# AIC      BIC   logLik deviance df.resid
# 194.6247 202.7261 -93.3124 186.6247       52
# Random effects:
#   Groups Name           Std.Dev.
# herd   (Intercept)    0.6634
# herd.1 period_numeric 0.0000
# Number of obs: 56, groups:  herd, 15
# Fixed Effects:
#   (Intercept)  period_numeric
# -0.8987         -0.5517

(rand.int.rand.slopes <-
    gnlrem(y=cbind(incidence, size - incidence),
           mu = ~ plogis(Intercept + period_numeric*b_p + rand1 + period_numeric*rand2),
           pmu = c(Intercept=-1, b_p=-1),
           pmix=c(var1=0.78, var2=0.40, corr12=0),
           p_uppb = c(  2,   2, 4.00, 4.00, 0),
           p_lowb = c( -2,  -2, 0.001, 0.00004, 0),
           distribution="binomial",
           nest=herd,
           random=c("rand1", "rand2"),
           mixture="bivariate-normal-corr",
           ooo=TRUE,
           compute_hessian = FALSE,
           compute_kkt = FALSE,
           trace=1,
           method='nlminb',
           int2dmethod="cuba"

    )
)
# [1] 5
# Intercept       b_p      var1      var2    corr12
# -1.00     -1.00      0.78      0.40      0.00
# [1] 100.0948
# fn is  fn
# Looking for method =  nlminb
# Function has  5  arguments
# par[ 1 ]:  -2   <? -1   <? 2     In Bounds
# par[ 2 ]:  -2   <? -1   <? 2     In Bounds
# par[ 3 ]:  0.001   <? 0.78   <? 4     In Bounds
# par[ 4 ]:  4e-05   <? 0.4   <? 4     In Bounds
# par[ 5 ]:  0   <? 0   <? 0     In Bounds
# Analytic gradient not made available.
# Analytic Hessian not made available.
# Scale check -- log parameter ratio= 0.39794   log bounds ratio= 0.0001085872
# Method:  nlminb
# 0:     100.09476: -1.00000 -1.00000 0.780000 0.400000  0.00000
# 1:     96.840530: -0.889769 -0.753339 0.781888 0.240449  0.00000
# 2:     93.411291: -0.921116 -0.611218 0.414012 4.00000e-05  0.00000
# 3:     93.361374: -0.915466 -0.592099 0.416240 0.0188371  0.00000
# 4:     93.269964: -0.909596 -0.586202 0.431892 4.00000e-05  0.00000
# 5:     93.268432: -0.902234 -0.569471 0.450896 0.00781833  0.00000
# 6:     93.245441: -0.904989 -0.573779 0.453331 4.00000e-05  0.00000
# 7:     93.244993: -0.906060 -0.560295 0.450888 4.00000e-05  0.00000
# 8:     93.240816: -0.909700 -0.566084 0.450205 4.00000e-05  0.00000
# 9:     93.238758: -0.922579 -0.565109 0.445501 4.00000e-05  0.00000
# 10:     93.238174: -0.922553 -0.562195 0.445604 4.00000e-05  0.00000
# 11:     93.237931: -0.925453 -0.562172 0.445910 4.00000e-05  0.00000
# 12:     93.237854: -0.927833 -0.560123 0.450825 4.00000e-05  0.00000
# 13:     93.237692: -0.931224 -0.558374 0.446415 4.00000e-05  0.00000
# 14:     93.237642: -0.932079 -0.560479 0.446628 4.00000e-05  0.00000
# 15:     93.237548: -0.933500 -0.559038 0.447681 4.00000e-05  0.00000
# 16:     93.237548: -0.933253 -0.559058 0.447741 4.00000e-05  0.00000
# 17:     93.237548: -0.933305 -0.559046 0.447725 4.00000e-05  0.00000
# 18:     93.237548: -0.933304 -0.559047 0.447725 4.00000e-05  0.00000
# Post processing for method  nlminb
# Successful convergence!
#   Save results from method  nlminb
# $par
# Intercept        b_p       var1       var2     corr12
# -0.9333043 -0.5590465  0.4477254  0.0000400  0.0000000
#
# $message
# [1] "relative convergence (4)"
#
# $convcode
# [1] 0
#
# $value
# [1] 93.23755
#
# $fevals
# function
# 25
#
# $gevals
# gradient
# 90
#
# $nitns
# [1] 18
#
# $kkt1
# [1] NA
#
# $kkt2
# [1] NA
#
# $xtimes
# user.self
# 2779.862
#
# Assemble the answers
# Intercept        b_p      var1  var2 corr12    value fevals gevals niter convcode kkt1 kkt2
# nlminb -0.9333043 -0.5590465 0.4477254 4e-05      0 93.23755     25     90    18        0   NA   NA
# xtime
# nlminb 2779.862

(rand.int.rand.slopes.nonzero.corr <-
    gnlrem(y=cbind(incidence, size - incidence),
           mu = ~ plogis(Intercept + period_numeric*b_p + rand1 + period_numeric*rand2),
           pmu = c(Intercept=-0.9333689, b_p=-0.558996),
           pmix=c(var1=0.3686549, var2=0.2531917, corr12= -0.10),
           p_uppb = c(  2,   2, 4.00, 4.00, 0.90),
           p_lowb = c( -2,  -2, 0.05, 0.04,-0.90),
           distribution="binomial",
           nest=herd,
           random=c("rand1", "rand2"),
           mixture="bivariate-normal-corr",
           ooo=TRUE,
           compute_hessian = FALSE,
           compute_kkt = FALSE,
           trace=1,
           method='nlminb',
           int2dmethod="cuba"
    )
)
# Intercept        b_p       var1       var2     corr12
# -0.9333689 -0.5589960  0.3686549  0.2531917 -0.1000000
# [1] 96.66891
# fn is  fn
# Looking for method =  nlminb
# Function has  5  arguments
# par[ 1 ]:  -2   <? -0.9333689   <? 2     In Bounds
# par[ 2 ]:  -2   <? -0.558996   <? 2     In Bounds
# par[ 3 ]:  0.05   <? 0.3686549   <? 4     In Bounds
# par[ 4 ]:  0.04   <? 0.2531917   <? 4     In Bounds
# par[ 5 ]:  -0.9   <? -0.1   <? 0.9     In Bounds
# Analytic gradient not made available.
# Analytic Hessian not made available.
# Scale check -- log parameter ratio= 0.9700533   log bounds ratio= 0.3467875
# Method:  nlminb
# 0:     96.668908: -0.933369 -0.558996 0.368655 0.253192 -0.100000
# 1:     93.764320: -0.925143 -0.600387 0.401532 0.0705280 -0.139629
# 2:     93.663329: -0.794393 -0.508991 0.587035 0.0400000 -0.403583
# 3:     93.368641: -0.929231 -0.476850 0.561841 0.0400000 -0.516573
# 4:     93.135918: -1.05063 -0.488291 0.606646 0.0400000 -0.391118
# 5:     93.061465: -0.987166 -0.544726 0.447698 0.0400000 -0.379357
# 6:     93.003421: -1.00495 -0.518751 0.510988 0.0400000 -0.413143
# 7:     92.980272: -1.01102 -0.512586 0.519505 0.0400000 -0.443200
# 8:     92.843109: -1.02957 -0.491499 0.572051 0.0400000 -0.613687
# 9:     92.690070: -1.03509 -0.488232 0.643033 0.0400000 -0.769081
# 10:     92.417551: -0.995270 -0.521391 0.751831 0.0400000 -0.878620
# 11:     92.301027: -0.921752 -0.570864 0.813047 0.0400000 -0.816547
# 12:     92.116208: -0.843831 -0.621220 0.952859 0.106774 -0.818105
# 13:     92.108905: -0.829141 -0.637993 0.988450 0.0980816 -0.819566
# 14:     92.098351: -0.851070 -0.621007 0.955995 0.0958164 -0.823257
# 15:     92.080683: -0.868195 -0.607672 0.990312 0.0975466 -0.837036
# 16:     92.078781: -0.876284 -0.610849 0.987859 0.0932515 -0.837657
# 17:     92.075977: -0.872079 -0.611615 0.995128 0.0986009 -0.836853
# 18:     92.075055: -0.867671 -0.619043  1.01318 0.0987198 -0.835824
# 19:     92.072255: -0.897652 -0.600761 0.995213 0.0929377 -0.839901
# 20:     92.071222: -0.898548 -0.599909 0.998049 0.0981050 -0.839632
# 21:     92.069745: -0.895347 -0.602814  1.00165 0.0959737 -0.839253
# 22:     92.068704: -0.888963 -0.603302  1.01083 0.100449 -0.839293
# 23:     92.067546: -0.912729 -0.599588  1.00949 0.100118 -0.840250
# 24:     92.065667: -0.911289 -0.597023  1.01329 0.0960900 -0.839790
# 25:     92.065181: -0.911061 -0.596072  1.01740 0.100429 -0.841468
# 26:     92.063686: -0.908216 -0.600430  1.02047 0.0988163 -0.842168
# 27:     92.062576: -0.904611 -0.599756  1.02507 0.100845 -0.843054
# 28:     92.061936: -0.908856 -0.597960  1.02725 0.0975579 -0.842813
# 29:     92.060866: -0.906395 -0.599946  1.03095 0.100915 -0.844201
# 30:     92.060197: -0.905387 -0.601117  1.03622 0.0986433 -0.844509
# 31:     92.059329: -0.908543 -0.596741  1.03836 0.0999210 -0.844653
# 32:     92.058422: -0.906826 -0.600893  1.04156 0.101522 -0.846202
# 33:     92.057351: -0.907661 -0.597961  1.04666 0.101453 -0.846465
# 34:     92.056886: -0.915429 -0.596897  1.05498 0.0984724 -0.847833
# 35:     92.053833: -0.912204 -0.599004  1.06449 0.103531 -0.851084
# 36:     92.053320: -0.910584 -0.595329  1.07546 0.105772 -0.851029
# 37:     92.050592: -0.909178 -0.598999  1.08645 0.106157 -0.853318
# 38:     92.047510: -0.902631 -0.602922  1.13649 0.114179 -0.862802
# 39:     92.046690: -0.917413 -0.599439  1.15437 0.114772 -0.861209
# 40:     92.046615: -0.918075 -0.597644  1.15287 0.113208 -0.864805
# 41:     92.046506: -0.914415 -0.599314  1.15109 0.114132 -0.862900
# 42:     92.046501: -0.915416 -0.599059  1.15193 0.114079 -0.862989
# 43:     92.046501: -0.915315 -0.599069  1.15184 0.114072 -0.863020
# 44:     92.046501: -0.915313 -0.599070  1.15185 0.114073 -0.863017
#
# Post processing for method  nlminb
# Successful convergence!
#   Save results from method  nlminb
# $par
# Intercept        b_p       var1       var2     corr12
# -0.9153133 -0.5990704  1.1518453  0.1140729 -0.8630170
#
# $message
# [1] "relative convergence (4)"
#
# $convcode
# [1] 0
#
# $value
# [1] 92.0465
#
# $fevals
# function
# 53
#
# $gevals
# gradient
# 264
#
# $nitns
# [1] 44
#
# $kkt1
# [1] NA
#
# $kkt2
# [1] NA
#
# $xtimes
# user.self
# 328.956
#
# Assemble the answers
# Intercept        b_p     var1      var2    corr12   value fevals gevals niter convcode kkt1 kkt2
# nlminb -0.9153133 -0.5990704 1.151845 0.1140729 -0.863017 92.0465     53    264    44        0   NA   NA
# xtime
# nlminb 328.956

########################################$$$$#########
## 2022-02-09                                      ##
## END two random parameters                       ##
#####################################################

Things (gnlrem vs glmer) look pretty similar!

swihart commented 2 years ago

I have a "big dataset" running on cluster; will post results here when it is done. I'm glad I figured out how to get pcubature working -- even with 16x speed up, it is going slowly. Can't imagine what romberg would have been!

swihart commented 2 years ago

Here it is. Put in some glmer and GLMMadaptive calls as well. Looks like we are in the neighborhood of those estimates. This is in the /Volumes/data/github/gnlrim/biowulf_cluster/cuba_romberg2d_tests_sameLikTest.out


R version 4.1.0 (2021-05-18) -- "Camp Pontanezen"
Copyright (C) 2021 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> ## withr::with_libpaths(new="./cuba_gnlrem",
> ##                      devtools::install_github("swihart/gnlrim")
> ##                      )
> library(gnlrim, lib.loc="./cuba_gnlrem")
> library(lme4)
> library(GLMMadaptive)
> ########################################$$$$#########
> ## 2022-02-10                                      ##
> ## START two random parameters                     ##
> #####################################################
> # Derived expression for gamma
> g <- function(a) {
+   iu <- complex(real=0, imaginary=1)
+   return(abs(1 - iu * tan(pi * a / 2)) ^ (-1 / a))
+ }
> 
> ## add in beta random effect
> sim_mrim_data <- function(n1, n2, J, a0, a1, v1=1,v2=2,rho=0.5, mrim="Logistic-BIVARIATE_NORM", alpha=1.89, gamma=1.2, delta=1){
+ 
+   if(mrim=="Logistic-BIVARIATE_NORM"){
+     print("HERE")
+     G <- function(n, v1, v2, rho){mvtnorm::rmvnorm(n=n, sigma=matrix(c(v1,rho*sqrt(v1*v2),rho*sqrt(v1*v2),v2),nrow=2))}
+     H <- function(x) plogis(x)
+   }
+ 
+   n <- n1 + n2
+   u <- round(apply(G(n,v1,v2,rho),2, function(W) rep(W,each=J)),2)
+ 
+   ## x <- c(rep(1, n1*J), rep(0, n2*J))
+   ##  x <-c(runif(n1*J, 0.5,1.5), runif(n2*J, 0,0.60))
+   x <-c(sample(c(1,2,3,4)/10,n1*J,TRUE), sample(c(1.2,2.2,3.2,4.2)/10,n2*J,TRUE))
+   eta <- round(a0 + a1*x,2)
+ 
+   eta_i <- round( (a0 + u[,1]) + (a1+u[,2])*x, 2)
+   py1 <- round(H(eta_i),2)
+   y <- rbinom(length(eta_i), 1, prob=py1 )
+ 
+   data.frame(id=rep(1:n, each=J),
+              j = rep(1:J),
+              x1 = x,
+              eta = eta,
+              u_i = u,
+              eta_i = eta_i,
+              py1 = py1,
+              y=y
+   )
+ 
+ }
> 
> 
> ##detach(summed_binom_dat)
> set.seed(5)
> binom_dat <-
+   sim_mrim_data(800,400, J=100, a0 = -2, a1 = 1)
[1] "HERE"
> #  sim_mrim_data(2,2, J=100, a0 = -2, a1 = 1)
> data.table::setDT(binom_dat)
> 
> summed_binom_dat <-
+   binom_dat[, {j=list(r=sum(y), n_r=sum(y==0))}, by=c("id","x1")]
> data.table::setkey(summed_binom_dat,id, x1)
> summed_binom_dat
        id   x1 r n_r
   1:    1 0.10 2  20
   2:    1 0.20 5  28
   3:    1 0.30 8  20
   4:    1 0.40 5  12
   5:    2 0.10 0  22
  ---                
4796: 1199 0.42 6  19
4797: 1200 0.12 1  23
4798: 1200 0.22 3  24
4799: 1200 0.32 1  22
4800: 1200 0.42 6  20
> 
> ## GLMMadaptive; first with bernoulli, then with binomial (increasing nAGQ)
> zeroes_ones <-
+ GLMMadaptive::mixed_model(fixed= y ~ x1,
+                           random = ~ x1 | id, ## ~ x1 || id,
+                           data=binom_dat,
+                           family=binomial())
> summary(zeroes_ones)

Call:
GLMMadaptive::mixed_model(fixed = y ~ x1, random = ~x1 | id, 
    data = binom_dat, family = binomial())

Data Descriptives:
Number of Observations: 120000
Number of Groups: 1200 

Model:
 family: binomial
 link: logit 

Fit statistics:
   log.Lik      AIC      BIC
 -52770.73 105551.5 105576.9

Random effects covariance matrix:
             StdDev    Corr
(Intercept)  0.9746        
x1           1.5622  0.5007

Fixed effects:
            Estimate Std.Err  z-value p-value
(Intercept)  -1.9210  0.0371 -51.7477 < 1e-04
x1            0.9779  0.0965  10.1345 < 1e-04

Integration:
method: adaptive Gauss-Hermite quadrature rule
quadrature points: 11

Optimization:
method: hybrid EM and quasi-Newton
converged: TRUE 
> summed_nAGQ_05 <-
+ GLMMadaptive::mixed_model(fixed= cbind(r,n_r) ~ x1,
+                           random = ~ x1 | id, ## ~ x1 || id,
+                           data=summed_binom_dat,
+                           family=binomial(),
+                           nAGQ=5)
> summary(summed_nAGQ_05)

Call:
GLMMadaptive::mixed_model(fixed = cbind(r, n_r) ~ x1, random = ~x1 | 
    id, data = summed_binom_dat, family = binomial(), nAGQ = 5)

Data Descriptives:
Number of Observations: 4800
Number of Groups: 1200 

Model:
 family: binomial
 link: logit 

Fit statistics:
   log.Lik      AIC      BIC
 -10816.32 21642.64 21668.09

Random effects covariance matrix:
             StdDev    Corr
(Intercept)  0.9741        
x1           1.5735  0.4955

Fixed effects:
            Estimate Std.Err  z-value p-value
(Intercept)  -1.9221  0.0372 -51.6953 < 1e-04
x1            0.9740  0.0966  10.0793 < 1e-04

Integration:
method: adaptive Gauss-Hermite quadrature rule
quadrature points: 5

Optimization:
method: hybrid EM and quasi-Newton
converged: TRUE 
> summed_nAGQ_11 <-
+ GLMMadaptive::mixed_model(fixed= cbind(r,n_r) ~ x1,
+                           random = ~ x1 | id, ## ~ x1 || id,
+                           data=summed_binom_dat,
+                           family=binomial())
> summary(summed_nAGQ_11)

Call:
GLMMadaptive::mixed_model(fixed = cbind(r, n_r) ~ x1, random = ~x1 | 
    id, data = summed_binom_dat, family = binomial())

Data Descriptives:
Number of Observations: 4800
Number of Groups: 1200 

Model:
 family: binomial
 link: logit 

Fit statistics:
   log.Lik     AIC      BIC
 -10816.15 21642.3 21667.75

Random effects covariance matrix:
             StdDev    Corr
(Intercept)  0.9747        
x1           1.5623  0.5006

Fixed effects:
            Estimate Std.Err  z-value p-value
(Intercept)  -1.9211  0.0371 -51.7450 < 1e-04
x1            0.9783  0.0965  10.1389 < 1e-04

Integration:
method: adaptive Gauss-Hermite quadrature rule
quadrature points: 11

Optimization:
method: hybrid EM and quasi-Newton
converged: TRUE 
> summed_nAGQ_21 <-
+ GLMMadaptive::mixed_model(fixed= cbind(r,n_r) ~ x1,
+                           random = ~ x1 | id, ## ~ x1 || id,
+                           data=summed_binom_dat,
+                           family=binomial(),
+                           nAGQ=21)
> summary(summed_nAGQ_21)

Call:
GLMMadaptive::mixed_model(fixed = cbind(r, n_r) ~ x1, random = ~x1 | 
    id, data = summed_binom_dat, family = binomial(), nAGQ = 21)

Data Descriptives:
Number of Observations: 4800
Number of Groups: 1200 

Model:
 family: binomial
 link: logit 

Fit statistics:
   log.Lik     AIC      BIC
 -10816.15 21642.3 21667.75

Random effects covariance matrix:
             StdDev    Corr
(Intercept)  0.9747        
x1           1.5623  0.5006

Fixed effects:
            Estimate Std.Err  z-value p-value
(Intercept)  -1.9211  0.0371 -51.7451 < 1e-04
x1            0.9783  0.0965  10.1388 < 1e-04

Integration:
method: adaptive Gauss-Hermite quadrature rule
quadrature points: 21

Optimization:
method: hybrid EM and quasi-Newton
converged: TRUE 
> ## summed_nAGQ_99 <-
> ## GLMMadaptive::mixed_model(fixed= cbind(r,n_r) ~ x1,
> ##                           random = ~ x1 | id, ## ~ x1 || id,
> ##                           data=summed_binom_dat,
> ##                           family=binomial(),
> ##                           nAGQ=99)
> ## summary(summed_nAGQ_99)
> 
> 
> ## glmer with correlation between random intercept and random slope
> lme4::glmer(cbind(r, n_r) ~ x1 + (x1 | id), summed_binom_dat, binomial, nAGQ = 0)
Generalized linear mixed model fit by maximum likelihood (Adaptive
  Gauss-Hermite Quadrature, nAGQ = 0) [glmerMod]
 Family: binomial  ( logit )
Formula: cbind(r, n_r) ~ x1 + (x1 | id)
   Data: summed_binom_dat
      AIC       BIC    logLik  deviance  df.resid 
 21658.18  21690.56 -10824.09  21648.18      4795 
Random effects:
 Groups Name        Std.Dev. Corr
 id     (Intercept) 0.9707       
        x1          1.5442   0.51
Number of obs: 4800, groups:  id, 1200
Fixed Effects:
(Intercept)           x1  
     -1.888        1.005  
> 
> 
> attach(summed_binom_dat)
> 
> ybind <- cbind(r,n_r)
> period_numeric <- x1
> 
> ## naive starting values (see other .out for starting values based on glmer)
> (rand.int.rand.slopes.nonzero.corr <-
+     gnlrem(y=ybind,
+            mu = ~ plogis(Intercept + period_numeric*b_p + rand1 + period_numeric*rand2),
+            pmu = c(Intercept = -0.01, b_p = 0.01),
+            pmix=c(var1=1, var2=1, corr12= 0.01),
+            p_uppb = c(  1,   2, 4.00, 4.00, 0.90),
+            p_lowb = c( -3,  -1, 0.05, 0.05,-0.90),
+            distribution="binomial",
+            nest=id,
+            random=c("rand1", "rand2"),
+            mixture="bivariate-normal-corr",
+            ooo=TRUE,
+            compute_hessian = FALSE,
+            compute_kkt = FALSE,
+            trace=1,
+            method='nlminb',
+            int2dmethod="cuba"
+     )
+ )
fn is  fn 
Looking for method =  nlminb 
Function has  5  arguments
par[ 1 ]:  -3   <? -0.01   <? 1     In Bounds  
par[ 2 ]:  -1   <? 0.01   <? 2     In Bounds  
par[ 3 ]:  0.05   <? 1   <? 4     In Bounds  
par[ 4 ]:  0.05   <? 1   <? 4     In Bounds  
par[ 5 ]:  -0.9   <? 0.01   <? 0.9     In Bounds  
Analytic gradient not made available.
Analytic Hessian not made available.
Scale check -- log parameter ratio= 2   log bounds ratio= 0.3467875 
Method:  nlminb 
  0:     12500.413: -0.0100000 0.0100000  1.00000  1.00000 0.0100000
  1:     11315.628: -0.796145 -0.0595843  1.57921  1.03912 0.210302
  2:     10953.431: -1.75783 0.187477  1.51454  1.13405 0.179955
  3:     10829.904: -1.90950 0.997787  1.22874  1.40413 0.587105
  4:     10824.880: -1.95599  1.02154  1.14586  1.41850 0.573017
  5:     10821.894: -1.94354  1.05721  1.05729  1.44528 0.569813
  6:     10820.624: -1.95558  1.08397 0.986604  1.50269 0.598915
  7:     10819.670: -1.92662  1.05815 0.988549  1.59418 0.609959
  8:     10818.453: -1.95844  1.02739 0.974054  1.78428 0.568812
  9:     10818.182: -1.92669  1.08054 0.999277  1.93889 0.460987
 10:     10817.684: -1.89037 0.969424 0.979647  2.05000 0.577625
 11:     10816.971: -1.92763 0.971929 0.961835  2.05359 0.568198
 12:     10816.879: -1.91212 0.992301 0.941792  2.07323 0.548949
 13:     10816.657: -1.92553 0.998135 0.952115  2.11140 0.542982
 14:     10816.583: -1.92343  1.00276 0.968238  2.15030 0.539123
 15:     10816.300: -1.91936 0.958633 0.961574  2.32331 0.510711
 16:     10816.287: -1.91712  1.00824 0.953534  2.49002 0.461464
 17:     10816.149: -1.92208 0.978332 0.948186  2.47106 0.497611
 18:     10816.145: -1.92130 0.978477 0.952787  2.46207 0.491957
 19:     10816.145: -1.92127 0.978949 0.951939  2.46442 0.492379
 20:     10816.145: -1.92128 0.978917 0.951943  2.46442 0.492391
Post processing for method  nlminb 
Successful convergence! 
Save results from method  nlminb 
$par
 Intercept        b_p       var1       var2     corr12 
-1.9212763  0.9789169  0.9519433  2.4644232  0.4923914 

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

$convcode
[1] 0

$value
[1] 10816.14

$fevals
function 
      26 

$gevals
gradient 
     123 

$nitns
[1] 20

$kkt1
[1] NA

$kkt2
[1] NA

$xtimes
user.self 
  45131.5 

Assemble the answers
       Intercept       b_p      var1     var2    corr12    value fevals gevals
nlminb -1.921276 0.9789169 0.9519433 2.464423 0.4923914 10816.14     26    123
       niter convcode kkt1 kkt2   xtime
nlminb    20        0   NA   NA 45131.5
> 
> ########################################$$$$#########
> ## 2022-02-16                                      ##
> ## END two random parameters                       ##
> #####################################################
>

swihart / gnlrim