Running 'Coxph-Ex.R' [12s]
  Comparing 'Coxph-Ex.Rout' to 'Coxph-Ex.Rout.save' ... OK
  Running 'KaplanMeier-Ex.R' [2s]
  Running 'PI_OVL-Ex.R' [161s]
  Running 'Polr-Ex.R' [161s]
  Running 'Survreg-Ex.R' [4s]
  Comparing 'Survreg-Ex.Rout' to 'Survreg-Ex.Rout.save' ... OK
  Running 'bugfixes.R' [196s]
  Running 'intercepts-Ex.R' [4s]
  Comparing 'intercepts-Ex.Rout' to 'intercepts-Ex.Rout.save' ... OK
  Running 'mmlt-Ex.R' [256s]
  Running 'mmlt_ll_sc.R' [183s]
  Running 'stram-Ex.R' [13s]
  Comparing 'stram-Ex.Rout' to 'stram-Ex.Rout.save' ... OK
Running the tests in 'tests/PI_OVL-Ex.R' failed.
Complete output:
  > 
  > library("tram")
  Loading required package: mlt
  Loading required package: basefun
  Loading required package: variables
  Loading required package: mvtnorm
  > library("multcomp")
  Loading required package: survival
  Loading required package: TH.data
  Loading required package: MASS
  
  Attaching package: 'TH.data'
  
  The following object is masked from 'package:MASS':
  
      geyser
  
  > 
  > m <- BoxCox(dist ~ speed, data = cars)
  > m$negative
  [1] TRUE
  > coef(m)
      speed 
  0.3062335 
  > nd <- data.frame(speed = as.double(1:5 * 5))
  > 
  > PI(m, newdata = nd)
           [,1]      [,2]      [,3]      [,4] [,5]
  [1,] 0.860529 0.9848214 0.9994191 0.9999926    1
  > lp15 <- c(predict(m, newdata = data.frame(speed = 15)))
  > PI(m, newdata = nd, reference = lp15)
             [,1]     [,2] [,3]     [,4]      [,5]
  [1,] 0.01517862 0.139471  0.5 0.860529 0.9848214
  > PI(m, newdata = nd, reference = nd[3,,drop = FALSE])
             [,1]     [,2] [,3]     [,4]      [,5]
  [1,] 0.01517862 0.139471  0.5 0.860529 0.9848214
  > PI(m, newdata = nd, reference = nd)
            1         2         3         4
  2 0.8605290                              
  3 0.9848214 0.8605290                    
  4 0.9994191 0.9848214 0.8605290          
  5 0.9999926 0.9994191 0.9848214 0.8605290
  > 
  > F <- m$model$todistr$p
  > Q <- m$model$todistr$q
  > f <- function(p, b) F(Q(p) - b)
  > lp <- predict(m, newdata = nd, type = "lp")
  > integrate(f, lower = 0, upper = 1, b = lp[2] - lp[1])
  0.139471 with absolute error < 0.00012
  > 
  > OVL(m, newdata = nd)
           [,1]     [,2]       [,3]        [,4]         [,5]
  [1,] 0.443924 0.125728 0.02163296 0.002196175 0.0001292314
  > lp15 <- c(predict(m, newdata = data.frame(speed = 15)))
  > OVL(m, newdata = nd, reference = lp15)
           [,1]     [,2] [,3]     [,4]     [,5]
  [1,] 0.125728 0.443924    1 0.443924 0.125728
  > OVL(m, newdata = nd, reference = nd[3,,drop = FALSE])
           [,1]     [,2] [,3]     [,4]     [,5]
  [1,] 0.125728 0.443924    1 0.443924 0.125728
  > OVL(m, newdata = nd, reference = nd)
              1           2           3           4
  2 0.443924015                                    
  3 0.125727999 0.443924015                        
  4 0.021632964 0.125727999 0.443924015            
  5 0.002196175 0.021632964 0.125727999 0.443924015
  > 
  > d <- m$model$todistr$d
  > Q <- m$model$todistr$q
  > f <- function(p, b) pmin(1, d(Q(p) - b) / d(Q(p)))
  > lp <- predict(m, newdata = nd, type = "lp")
  > integrate(f, lower = 0, upper = 1, b = lp[2] - lp[1])
  0.4439239 with absolute error < 7.3e-05
  > 
  > m <- Colr(dist ~ speed, data = cars)
  > m$negative
  [1] FALSE
  > coef(m)
       speed 
  -0.5293216 
  > nd <- data.frame(speed = as.double(1:5 * 5))
  > 
  > PI(m, newdata = nd)
            [,1]      [,2]      [,3]      [,4]      [,5]
  [1,] 0.8589555 0.9781801 0.9975257 0.9997579 0.9999781
  > lp15 <- c(predict(m, newdata = data.frame(speed = 15)))
  > PI(m, newdata = nd, reference = lp15)
             [,1]      [,2] [,3]      [,4]      [,5]
  [1,] 0.02181988 0.1410445  0.5 0.8589555 0.9781801
  > PI(m, newdata = nd, reference = nd[3,,drop = FALSE])
             [,1]      [,2] [,3]      [,4]      [,5]
  [1,] 0.02181988 0.1410445  0.5 0.8589555 0.9781801
  > PI(m, newdata = nd, reference = nd)
            1         2         3         4
  2 0.8589555                              
  3 0.9781801 0.8589555                    
  4 0.9975257 0.9781801 0.8589555          
  5 0.9997579 0.9975257 0.9781801 0.8589555
  > 
  > F <- m$model$todistr$p
  > Q <- m$model$todistr$q
  > f <- function(p, b) F(Q(p) + b)
  > lp <- predict(m, newdata = nd, type = "lp")
  > integrate(f, lower = 0, upper = 1, b = lp[2] - lp[1])
  0.1410445 with absolute error < 9.2e-05
  > 
  > 
  > OVL(m, newdata = nd)
            [,1]      [,2]       [,3]       [,4]        [,5]
  [1,] 0.4205382 0.1323967 0.03705084 0.01000088 0.002672582
  > lp15 <- c(predict(m, newdata = data.frame(speed = 15)))
  > OVL(m, newdata = nd, reference = lp15)
            [,1]      [,2] [,3]      [,4]      [,5]
  [1,] 0.1323967 0.4205382    1 0.4205382 0.1323967
  > OVL(m, newdata = nd, reference = nd[3,,drop = FALSE])
            [,1]      [,2] [,3]      [,4]      [,5]
  [1,] 0.1323967 0.4205382    1 0.4205382 0.1323967
  > OVL(m, newdata = nd, reference = nd)
             1          2          3          4
  2 0.42053821                                 
  3 0.13239674 0.42053821                      
  4 0.03705084 0.13239674 0.42053821           
  5 0.01000088 0.03705084 0.13239674 0.42053821
  > 
  > d <- m$model$todistr$d
  > Q <- m$model$todistr$q
  > f <- function(p, b) pmin(1, d(Q(p) + b) / d(Q(p)))
  > lp <- predict(m, newdata = nd, type = "lp")
  > integrate(f, lower = 0, upper = 1, b = lp[2] - lp[1])
  0.4205378 with absolute error < 3.4e-05
  > 
  > 
  > m <- Coxph(dist ~ speed, data = cars)
  > m$negative
  [1] FALSE
  > coef(m)
      speed 
  -0.280805 
  > nd <- data.frame(speed = as.double(1:5 * 5))
  > 
  > PI(m, newdata = nd)
            [,1]      [,2]      [,3]      [,4]      [,5]
  [1,] 0.8028218 0.9431093 0.9854007 0.9963744 0.9991071
  > lp15 <- c(predict(m, newdata = data.frame(speed = 15)))
  > PI(m, newdata = nd, reference = lp15)
             [,1]      [,2] [,3]      [,4]      [,5]
  [1,] 0.05689072 0.1971782  0.5 0.8028218 0.9431093
  > PI(m, newdata = nd, reference = nd[3,,drop = FALSE])
             [,1]      [,2] [,3]      [,4]      [,5]
  [1,] 0.05689072 0.1971782  0.5 0.8028218 0.9431093
  > PI(m, newdata = nd, reference = nd)
            1         2         3         4
  2 0.8028218                              
  3 0.9431093 0.8028218                    
  4 0.9854007 0.9431093 0.8028218          
  5 0.9963744 0.9854007 0.9431093 0.8028218
  > 
  > F <- m$model$todistr$p
  > Q <- m$model$todistr$q
  > f <- function(p, b) F(Q(p) + b)
  > lp <- predict(m, newdata = nd, type = "lp")
  > integrate(f, lower = 0, upper = 1, b = lp[2] - lp[1])
  0.1971782 with absolute error < 2.2e-16
  > 
  > 
  > OVL(m, newdata = nd)
            [,1]      [,2]       [,3]       [,4]        [,5]
  [1,] 0.5223833 0.2153214 0.07528464 0.02386655 0.007148041
  > lp15 <- c(predict(m, newdata = data.frame(speed = 15)))
  > OVL(m, newdata = nd, reference = lp15)
            [,1]      [,2] [,3]      [,4]      [,5]
  [1,] 0.2153214 0.5223833    1 0.5223833 0.2153214
  > OVL(m, newdata = nd, reference = nd[3,,drop = FALSE])
            [,1]      [,2] [,3]      [,4]      [,5]
  [1,] 0.2153214 0.5223833    1 0.5223833 0.2153214
  > OVL(m, newdata = nd, reference = nd)
             1          2          3          4
  2 0.52238327                                 
  3 0.21532141 0.52238327                      
  4 0.07528464 0.21532141 0.52238327           
  5 0.02386655 0.07528464 0.21532141 0.52238327
  > 
  > d <- m$model$todistr$d
  > Q <- m$model$todistr$q
  > f <- function(p, b) pmin(1, d(Q(p) + b) / d(Q(p)))
  > lp <- predict(m, newdata = nd, type = "lp")
  > integrate(f, lower = 0, upper = 1, b = lp[2] - lp[1])
  0.5223832 with absolute error < 2.8e-06
  > 
  > 
  > m <- Lehmann(dist ~ speed, data = cars)
  > m$negative
  [1] TRUE
  > coef(m)
      speed 
  0.3478323 
  > nd <- data.frame(speed = as.double(1:5 * 5))
  > 
  > PI(m, newdata = nd)
            [,1]      [,2]      [,3]      [,4]      [,5]
  [1,] 0.8505805 0.9700647 0.9946083 0.9990486 0.9998327
  > lp15 <- c(predict(m, newdata = data.frame(speed = 15)))
  > PI(m, newdata = nd, reference = lp15)
             [,1]      [,2] [,3]      [,4]      [,5]
  [1,] 0.02993533 0.1494195  0.5 0.8505805 0.9700647
  > PI(m, newdata = nd, reference = nd[3,,drop = FALSE])
             [,1]      [,2] [,3]      [,4]      [,5]
  [1,] 0.02993533 0.1494195  0.5 0.8505805 0.9700647
  > PI(m, newdata = nd, reference = nd)
            1         2         3         4
  2 0.8505805                              
  3 0.9700647 0.8505805                    
  4 0.9946083 0.9700647 0.8505805          
  5 0.9990486 0.9946083 0.9700647 0.8505805
  > 
  > F <- m$model$todistr$p
  > Q <- m$model$todistr$q
  > f <- function(p, b) F(Q(p) - b)
  > lp <- predict(m, newdata = nd, type = "lp")
  > integrate(f, lower = 0, upper = 1, b = lp[2] - lp[1])
  0.1494195 with absolute error < 6.7e-15
  > 
  > 
  > OVL(m, newdata = nd)
            [,1]      [,2]       [,3]        [,4]        [,5]
  [1,] 0.4309584 0.1324664 0.03330627 0.007555076 0.001620911
  > lp15 <- c(predict(m, newdata = data.frame(speed = 15)))
  > OVL(m, newdata = nd, reference = lp15)
            [,1]      [,2] [,3]      [,4]      [,5]
  [1,] 0.1324664 0.4309584    1 0.4309584 0.1324664
  > OVL(m, newdata = nd, reference = nd[3,,drop = FALSE])
            [,1]      [,2] [,3]      [,4]      [,5]
  [1,] 0.1324664 0.4309584    1 0.4309584 0.1324664
  > OVL(m, newdata = nd, reference = nd)
              1           2           3           4
  2 0.430958441                                    
  3 0.132466424 0.430958441                        
  4 0.033306271 0.132466424 0.430958441            
  5 0.007555076 0.033306271 0.132466424 0.430958441
  > 
  > d <- m$model$todistr$d
  > Q <- m$model$todistr$q
  > f <- function(p, b) pmin(1, d(Q(p) - b) / d(Q(p)))
  > lp <- predict(m, newdata = nd, type = "lp")
  > integrate(f, lower = 0, upper = 1, b = lp[2] - lp[1])
  0.4309575 with absolute error < 5.2e-05
  > 
  > OVL(m, newdata = nd)
            [,1]      [,2]       [,3]        [,4]        [,5]
  [1,] 0.4309584 0.1324664 0.03330627 0.007555076 0.001620911
  > OVL(m, newdata = nd, conf.level = .95)
       Estimate          lwr        upr
  1 0.430958441 0.3206993334 0.56433049
  2 0.132466424 0.0632140413 0.26158351
  3 0.033306271 0.0098805486 0.10482119
  4 0.007555076 0.0013940936 0.03827161
  5 0.001620911 0.0001862683 0.01320327
  attr(,"conf.level")
  [1] 0.95
  attr(,"calpha")
  [1] 1.960371
  > 
  > lp15 <- c(predict(m, newdata = data.frame(speed = 15)))
  > OVL(m, newdata = nd, reference = lp15)
            [,1]      [,2] [,3]      [,4]      [,5]
  [1,] 0.1324664 0.4309584    1 0.4309584 0.1324664
  > OVL(m, newdata = nd, reference = lp15, conf.level = .95)
        Estimate          lwr         upr
  1 0.1251517592 5.851127e-02 0.252039347
  2 0.0312767073 9.088354e-03 0.100488045
  3 0.0070724877 1.278423e-03 0.036576212
  4 0.0015146035 1.705034e-04 0.012594431
  5 0.0003135728 2.197845e-05 0.004187655
  attr(,"conf.level")
  [1] 0.95
  attr(,"calpha")
  [1] 1.960371
  > 
  > OVL(m, newdata = nd[-3,,drop = FALSE], reference = nd[3,,drop = FALSE])
            [,1]      [,2]      [,3]      [,4]
  [1,] 0.1324664 0.4309584 0.4309584 0.1324664
  > OVL(m, newdata = nd[-3,,drop = FALSE], reference = nd[3,,drop = FALSE], conf.level = .95)
       Estimate        lwr       upr
  3-1 0.1324664 0.06321515 0.2615796
  3-2 0.4309584 0.32070160 0.5643272
  3-4 0.4309584 0.32070160 0.5643272
  3-5 0.1324664 0.06321515 0.2615796
  attr(,"conf.level")
  [1] 0.95
  attr(,"calpha")
  [1] 1.960326
  > 
  > OVL(m, newdata = nd, reference = nd)
              1           2           3           4
  2 0.430958441                                    
  3 0.132466424 0.430958441                        
  4 0.033306271 0.132466424 0.430958441            
  5 0.007555076 0.033306271 0.132466424 0.430958441
  > OVL(m, newdata = nd, reference = nd, conf.level = .95)
         Estimate         lwr        upr
  1-2 0.430958441 0.320699820 0.56432978
  1-3 0.132466424 0.063214279 0.26158268
  2-3 0.430958441 0.320699820 0.56432978
  1-4 0.033306271 0.009880609 0.10482062
  2-4 0.132466424 0.063214279 0.26158268
  3-4 0.430958441 0.320699820 0.56432978
  1-5 0.007555076 0.001394105 0.03827131
  2-5 0.033306271 0.009880609 0.10482062
  3-5 0.132466424 0.063214279 0.26158268
  4-5 0.430958441 0.320699820 0.56432978
  attr(,"conf.level")
  [1] 0.95
  attr(,"calpha")
  [1] 1.960362
  > 
  > OVL(m, newdata = nd, reference = nd, conf.level = .95, 
  +     calpha = univariate_calpha())
         Estimate         lwr        upr
  1-2 0.430958441 0.320719822 0.56430050
  1-3 0.132466424 0.063224062 0.26154855
  2-3 0.430958441 0.320719822 0.56430050
  1-4 0.033306271 0.009883093 0.10479729
  2-4 0.132466424 0.063224062 0.26154855
  3-4 0.430958441 0.320719822 0.56430050
  1-5 0.007555076 0.001394590 0.03825911
  2-5 0.033306271 0.009883093 0.10479729
  3-5 0.132466424 0.063224062 0.26154855
  4-5 0.430958441 0.320719822 0.56430050
  attr(,"conf.level")
  [1] 0.95
  attr(,"calpha")
  [1] 1.959964
  > OVL(m, newdata = nd, reference = nd, conf.level = .95, 
  +     calpha = adjusted_calpha())
         Estimate         lwr        upr
  1-2 0.430958441 0.320699820 0.56432978
  1-3 0.132466424 0.063214279 0.26158268
  2-3 0.430958441 0.320699820 0.56432978
  1-4 0.033306271 0.009880609 0.10482062
  2-4 0.132466424 0.063214279 0.26158268
  3-4 0.430958441 0.320699820 0.56432978
  1-5 0.007555076 0.001394105 0.03827131
  2-5 0.033306271 0.009880609 0.10482062
  3-5 0.132466424 0.063214279 0.26158268
  4-5 0.430958441 0.320699820 0.56432978
  attr(,"conf.level")
  [1] 0.95
  attr(,"calpha")
  [1] 1.960362
  > 
  > proc.time()
     user  system elapsed 
     2.43    0.43    2.85 
Running the tests in 'tests/Polr-Ex.R' failed.
Complete output:
  > 
  > library("MASS")
  > library("tram")
  Loading required package: mlt
  Loading required package: basefun
  Loading required package: variables
  Loading required package: mvtnorm
  > 
  > ### Windows diffs...
  > options(digits = 3)
  > 
  > tol <- .Machine$double.eps^(1/4)
  > 
  > cmp <- function(x, y)
  +     stopifnot(isTRUE(all.equal(x, y, tolerance = tol)))
  > 
  > (house.plr <- polr(Sat ~ Infl + Type + Cont, weights = Freq, data = housing))
  Call:
  polr(formula = Sat ~ Infl + Type + Cont, data = housing, weights = Freq)
  
  Coefficients:
     InflMedium      InflHigh TypeApartment    TypeAtrium   TypeTerrace 
          0.566         1.289        -0.572        -0.366        -1.091 
       ContHigh 
          0.360 
  
  Intercepts:
   Low|Medium Medium|High 
       -0.496       0.691 
  
  Residual Deviance: 3479.15 
  AIC: 3495.15 
  > summary(house.plr)
  
  Re-fitting to get Hessian
  
  Call:
  polr(formula = Sat ~ Infl + Type + Cont, data = housing, weights = Freq)
  
  Coefficients:
                 Value Std. Error t value
  InflMedium     0.566     0.1047    5.41
  InflHigh       1.289     0.1272   10.14
  TypeApartment -0.572     0.1192   -4.80
  TypeAtrium    -0.366     0.1552   -2.36
  TypeTerrace   -1.091     0.1515   -7.20
  ContHigh       0.360     0.0955    3.77
  
  Intercepts:
              Value  Std. Error t value
  Low|Medium  -0.496  0.125     -3.974 
  Medium|High  0.691  0.125      5.505 
  
  Residual Deviance: 3479.15 
  AIC: 3495.15 
  > 
  > (house.plr2 <- Polr(Sat ~ Infl + Type + Cont, weights = Freq, data = housing))
  
    Proportional Odds Regression Model 
  
  Call:
  Polr(formula = Sat ~ Infl + Type + Cont, data = housing, weights = Freq)
  
  Coefficients:
     InflMedium      InflHigh TypeApartment    TypeAtrium   TypeTerrace 
          0.566         1.289        -0.572        -0.366        -1.091 
       ContHigh 
          0.360 
  
  Log-Likelihood:
   -1740 (df = 8)
  
  > summary(house.plr2)
  
    Proportional Odds Regression Model 
  
  Call:
  Polr(formula = Sat ~ Infl + Type + Cont, data = housing, weights = Freq)
  
  Coefficients:
                Estimate Std. Error z value Pr(>|z|)    
  InflMedium      0.5664     0.1047    5.41  6.2e-08 ***
  InflHigh        1.2888     0.1272   10.14  < 2e-16 ***
  TypeApartment  -0.5724     0.1192   -4.80  1.6e-06 ***
  TypeAtrium     -0.3662     0.1552   -2.36  0.01828 *  
  TypeTerrace    -1.0910     0.1515   -7.20  5.9e-13 ***
  ContHigh        0.3603     0.0955    3.77  0.00016 ***
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  
  Log-Likelihood:
   -1740 (df = 8)
  Likelihood-ratio Test: Chisq = 170 on 6 degrees of freedom; p = <2e-16
  
  > 
  > cmp(coef(house.plr), coef(house.plr2))
  > ll <- logLik(house.plr)
  > attr(ll, "nobs") <- NULL
  > cmp(ll, logLik(house.plr2))
  > 
  > proc.time()
     user  system elapsed 
     1.59    0.29    1.87 
Running the tests in 'tests/bugfixes.R' failed.
Complete output:
  > 
  > library("tram")
  Loading required package: mlt
  Loading required package: basefun
  Loading required package: variables
  Loading required package: mvtnorm
  > library("sandwich")
  > library("survival")
  > library("numDeriv")
  > 
  > ### scores for fixed parameters
  > cars$int <- 1
  > a0 <- Lm(dist ~ speed, data = cars)
  > a1 <- Lm(dist ~ speed + int, data = cars, fixed = c("int" = 0))
  > 
  > s0 <- estfun(a0)
  > s1 <- estfun(a1)
  > s2 <- estfun(a1, parm = coef(as.mlt(a1), fixed = TRUE))
  > 
  > stopifnot(all.equal(s0[,"(Intercept)"], -s2[,"int"]))
  > stopifnot(all.equal(s1[,"(Intercept)"], -s2[,"int"]))
  > stopifnot(all.equal(s2[,"(Intercept)"], -s2[,"int"]))
  > 
  > ### log_first for count data (by Sandra Siegfried)
  > ### use: y + 1; log_first = TRUE, support = c(1, ...), bounds = c(1, ...)
  > library("MASS")
  > set.seed(29)
  > Nsim <- 100
  > b1 <- - 4.5
  > b0 <- 5
  > theta <- 2
  > h <- qlogis(ppois(0:100, lambda = 5))
  > dgp <- function(n = 4000){
  +   x <- runif(n, min = 0, max = 1) 
  +   log.mu <- b0 + b1 * x
  +   h.m <- matrix(h, nrow = length(h), ncol = length(x)) 
  +   p <- (plogis(t(h.m) - b1 * x) - runif(n))^2
  +   y <- max.col(-p) - 1
  +   m <- Colr(y ~ x, data = data.frame(y = y, x = x), 
  +             bounds = c(0L, Inf), fixed = c("x" = -b1),
  +             support = c(0L, floor(quantile(y, .9))), order = 10)
  +   
  +   ret <- data.frame(x = x, y = as.integer(y))
  +   attr(ret, "mC") <- m
  +   ret
  + }
  > d <- dgp()
  > d$y.p1 <- d$y + 1L
  > 
  > m1 <- Colr(y ~ x, data = d,
  +                   support = c(0L, as.numeric(max(d$y))),
  +                   bounds = c(0L, Inf),
  +                   order = 10)
  > 
  > try(m2a <- Colr(y ~ x, data = d,
  +                   support = c(1L, as.numeric(max(d$y))),
  +                   bounds = c(1L, Inf),
  +                   log_first = TRUE,
  +                   order = 10))
  Error in mlt(model = structure(list(model = structure(list(bresponse = structure(function (data,  : 
    check(responsevar, data) is not TRUE
  > 
  > try(m2b <- Colr(y ~ x, data = d,
  +                   support = c(1L, as.numeric(max(d$y))),
  +                   bounds = c(0L, Inf),
  +                   log_first = TRUE,
  +                   order = 10))
  Error in mlt(model = structure(list(model = structure(list(bresponse = structure(function (data,  : 
    check(responsevar, data) is not TRUE
  > 
  > m3 <- Colr(y.p1 ~ x, data = d,
  +                   support = c(1L, as.numeric(max(d$y.p1))),
  +                   bounds = c(1L, Inf),
  +                   log_first = TRUE,
  +                   order = 10)
  > 
  > l1 <- m1$logliki(coef(as.mlt(m1)), rep(1, nrow(d)))
  > l3 <- m3$logliki(coef(as.mlt(m3)), rep(1, nrow(d)))
  > stopifnot(cor(l1, l3) > .96)
  > 
  > ### normal CDF approximation by Matic et al (2018)
  > x <- -600:600 / 100
  > p1 <- .Call("R_pnormMRS", x)
  > p2 <- pnorm(x)
  > stopifnot(max(abs(p1 - p2)) < 5.79e-6) ### max absolute error
  > 
  > ### 0.3-0
  > data("GBSG2", package = "TH.data")
  > ### this model included an additional intercept
  > m1 <- Coxph(Surv(time, cens) | menostat:tgrade ~ horTh, data = GBSG2)
  > m2 <- Coxph(Surv(time, cens) | 0 + menostat:tgrade ~ horTh, data = GBSG2)
  > stopifnot(max(abs(coef(as.mlt(m1)) - coef(as.mlt(m2)))) < 
  +           sqrt(.Machine$double.eps))
  > ### interaction term with stratum
  > m3 <- Coxph(Surv(time, cens) | horTh ~ menostat + menostat:horTh, 
  +             data = GBSG2)
  > ci <- confint(m3)["menostatPre:horThyes",]
  > stopifnot(all(!is.finite(ci)))
  > 
  > ### problems with responses of class "R", spotted by Balint Tamasi
  > data("wine", package = "ordinal")
  > erating <- wine$rating
  > lrating <- erating
  > rrating <- erating
  > l9 <- lrating[wine$judge == 9]
  > l9[l9 > 1] <- levels(l9)[unclass(l9[l9 > 1]) - 1]
  > r9 <- rrating[wine$judge == 9]
  > r9[r9 < 5] <- levels(r9)[unclass(r9[r9 < 5]) + 1]
  > lrating[wine$judge != 9] <- rrating[wine$judge != 9] <- NA
  > erating[wine$judge == 9] <- NA
  > lrating[wine$judge == 9] <- l9
  > rrating[wine$judge == 9] <- r9
  > which(wine$judge == 9)
  [1] 65 66 67 68 69 70 71 72
  > wine$crating <- R(erating, cleft = lrating, cright = rrating)
  > ### gave an error
  > m <- Polr(crating ~ temp, data = wine, method = "probit")
  > 
  > ### another one
  > data("GBSG2", package = "TH.data")
  > tmp <- GBSG2
  > tmp$y <- R(with(GBSG2, Surv(time, cens)))
  > m1 <- Coxph(Surv(time, cens) ~ horTh, data = tmp)
  > m2 <- Coxph(y ~ horTh, data = tmp)
  > stopifnot(all.equal(coef(as.mlt(m1)), coef(as.mlt(m2)), 
  +                     check.attributes = FALSE))
  > 
  > ### check probabilistic index <-> log-odds ratio conversion
  > stopifnot(max(abs(PI(prob = PI(-15:15)) - (-15:15))) < 1e5)
  > 
  > ### check updating with permutations
  > data("BostonHousing2", package = "mlbench")
  > m <- Colr(Surv(cmedv, cmedv < 50) ~ chas + crim, data = BostonHousing2)
  > mm <- as.mlt(m)
  > p <- sample(1:NROW(BostonHousing2))
  > ## model with crim permuted
  > m1 <- update(mm, perm = "crim", permutation = p)
  > ## permute data directly
  > tmp <- BostonHousing2
  > tmp[, "crim"] <- tmp[p, "crim"]
  > m2 <- Colr(Surv(cmedv, cmedv < 50) ~ chas + crim, data = tmp)
  > stopifnot(all.equal(coef(m1), coef(as.mlt(m2)), tol = 1e-3))
  > stopifnot(all.equal(logLik(m1), logLik(m2), tol = 1e-6))
  > 
  > ### contraints, by Lucas Kook
  > data("GBSG2", package = "TH.data")
  > # gave an error
  > m <- Survreg(Surv(time, cens) ~ horTh + age, data = GBSG2, constraints = c("age >= 0"))
  > 
  > ### mtram with interval censoring, spotted by Sandra Siegfried
  > dir <- system.file("rda", package = "TH.data")
  > load(file.path(dir, "Primary_endpoint_data.rda"))
  > trt <- "randarm5-FU + Oxaliplatin"
  > ### convert "exact" event dates to interval-censoring (+/- one day)
  > tmp <- CAOsurv$iDFS
  > exact <- tmp[,3] == 1 
  > tmp[exact,2] <- tmp[exact,1] + 2
  > tmp[exact,1] <- pmax(tmp[exact,1] - 2, 0)
  > tmp[exact,3] <- 3
  > CAOsurv$iDFS2 <- tmp
  > CAO_SR <- Survreg(iDFS2 ~ randarm, data = CAOsurv, support = c(1, 1700), bounds = c(0, Inf))
  > CAO_SR_mtram <- mtram(CAO_SR, ~ (1 | Block), data = CAOsurv)
  > CAO_Cox <- Coxph(iDFS2 ~ randarm, data = CAOsurv, support = c(1, 1700), bounds = c(0, Inf), log_first = TRUE, order = 1)
  > CAO_Cox_mtram <- mtram(CAO_Cox, ~ (1 | Block), data = CAOsurv)
  > ### wasn't equal
  > stopifnot(isTRUE(all.equal(logLik(CAO_SR_mtram), logLik(CAO_Cox_mtram), tol = 1e-5)))
  > 
  > ## produces negative variances in shift-scale model
  > set.seed(100)
  > N <- 100
  > sc <- -5.5
  > sh <- 3
  >  
  > FZ <- pnorm
  > FZi <- qnorm
  > h2y <- function(y) log(-FZ(y, lower.tail = FALSE, log.p = TRUE))
  >  
  > x <- rep(x0 <- gl(2, 1), each = N)
  > xx <- (0:1)[x]
  >  
  > U <- runif(length(x))
  > scale <- sqrt(exp(sc * xx))
  > shift <- sh * xx
  > d <- data.frame(y = h2y( shift + FZi(U) / scale),
  +                 x = x)
  > m <- BoxCox(y ~ x | x, data = d, scale_shift = TRUE)
  > coef(m)
         x2    scl_x2 
  -4.562880 -5.168216 
  > try(diag(vcov(m)))
  Error in vcov.mlt(as.mlt(object)) : Hessian is not invertible
  > 
  > ### constraints in remove intercept were incorrect
  > ### spotted by Balint
  > stopifnot(isTRUE(all.equal(
  +   coef(as.mlt(Lm(dist ~ 1, data = cars, remove_intercept = FALSE)))[2:1] * c(-1, 1),
  +   coef(as.mlt(Lm(dist ~ 1, data = cars))), check.attributes = FALSE)))
  > stopifnot(isTRUE(all.equal(
  +   coef(as.mlt(Survreg(dist ~ 1, data = cars, remove_intercept = FALSE)))[2:1] * c(-1, 1),
  +   coef(as.mlt(Survreg(dist ~ 1, data = cars))), check.attributes = FALSE)))
  > 
  > ### mmlt with shift-scale had incorrect gradients
  > m1 <- Lm(Sepal.Width ~ Petal.Length | Petal.Width, data = iris)
  > m2 <- Lm(Sepal.Length ~ Petal.Length | Petal.Width, data = iris)
  > m0 <- mmlt(m1, m2, data = iris, formula = ~ 1, domargins = FALSE)
  > cf <- coef(m0)
  > m <- mmlt(m1, m2, data = iris, formula = ~ 1, dofit = FALSE)
  > stopifnot(isTRUE(all.equal(c(m$score(cf)), c(grad(m$ll, cf)), check.attributes = FALSE)))
  > 
  > m1 <- Lm(Sepal.Width ~ Petal.Length | Petal.Width, data = iris, scale_shift = TRUE)
  > m2 <- Lm(Sepal.Length ~ Petal.Length | Petal.Width, data = iris, scale_shift = TRUE)
  > m0 <- mmlt(m1, m2, data = iris, formula = ~ 1, domargins = FALSE)
  > cf <- coef(m0)
  > m <- mmlt(m1, m2, data = iris, formula = ~ 1, dofit = FALSE)
  > stopifnot(isTRUE(all.equal(c(m$score(cf)), c(grad(m$ll, cf)), check.attributes = FALSE)))
  > 
  > 
  > proc.time()
     user  system elapsed 
    34.93    4.54   39.54 
Running the tests in 'tests/mmlt-Ex.R' failed.
Complete output:
  > library("tram")
  Loading required package: mlt
  Loading required package: basefun
  Loading required package: variables
  Loading required package: mvtnorm
  > library("mvtnorm")
  > library("multcomp")
  Loading required package: survival
  Loading required package: TH.data
  Loading required package: MASS
  
  Attaching package: 'TH.data'
  
  The following object is masked from 'package:MASS':
  
      geyser
  
  > library("sandwich")
  > 
  > options(digits = 2)
  > 
  > set.seed(25)
  > chk <- function(...) all.equal(..., tol = 1e-3, check.attributes = FALSE)
  > 
  > OR <- 1
  > 
  > J <- 4
  > N <- 100
  > S <- cov2cor(tcrossprod(matrix(runif(J * J), ncol = J)))
  > y <- rmvnorm(N, sigma = S)
  > u <- as.data.frame(plogis(y))
  > x <- runif(N)
  > d <- cbind(u, x)
  > un <- colnames(d)[1:J]
  > 
  > m <- lapply(un, function(i)
  +     BoxCox(as.formula(paste(i, "~ x")), data = d, bounds = c(0, 1), support = c(0, 1), order = OR))
  > m$data <- d
  > m$formula <- ~ 1
  > mm <- do.call("mmlt", m)
  > 
  > chk(c(logLik(mm)), sum(predict(mm, newdata = d, type = "density", log = TRUE)))
  [1] TRUE
  > L <- as.array(coef(mm, type = "Lambda"))[,,1]
  > chk(as.array(coef(mm, type = "Lambdainv"))[,,1], solve(L))
  [1] TRUE
  > chk(as.array(coef(mm, type = "Sigma"))[,,1], tcrossprod(solve(L)))
  [1] TRUE
  > chk(as.array(coef(mm, type = "Cor"))[,,1], cov2cor(tcrossprod(solve(L))))
  [1] TRUE
  > 
  > chk(colSums(estfun(mm)), mm$score(coef(mm, type = "all")))
  [1] TRUE
  > 
  > ### marginal normal
  > m$conditional <- FALSE
  > mmN <- do.call("mmlt", m)
  > 
  > chk(logLik(mm), logLik(mmN))
  [1] TRUE
  > chk(c(logLik(mmN)), sum(predict(mmN, newdata = d, type = "density", log = TRUE)))
  [1] TRUE
  > 
  > cf1 <- do.call("c", lapply(m[1:J], function(x) coef(as.mlt(x))))
  > cf2 <- coef(mmN)[1:length(cf1)]
  > # cbind(cf1, cf2)
  > 
  > sd1 <- sqrt(do.call("c", lapply(m[1:J], function(x) diag(vcov(as.mlt(x))))))
  > sd2 <- sqrt(diag(vcov(mmN)))[1:length(sd1)]
  > 
  > # cbind(sd1, sd2)
  > vcov(mmN)["V1.x", "V4.x"]
  [1] 0.054
  > 
  > 
  > chk(as.array(coef(mm, type = "Lambda"))[,,1], 
  +     as.array(coef(mmN, type = "Lambda"))[,,1])
  [1] TRUE
  > chk(as.array(coef(mm, type = "Lambdainv"))[,,1], 
  +     as.array(coef(mmN, type = "Lambdainv"))[,,1])
  [1] TRUE
  > chk(as.array(coef(mm, type = "Sigma"))[,,1], 
  +     as.array(coef(mmN, type = "Sigma"))[,,1])
  [1] TRUE
  > chk(as.array(coef(mm, type = "Spearman"))[,,1], 
  +     as.array(coef(mmN, type = "Spearman"))[,,1])
  [1] TRUE
  > 
  > chk(predict(mm, newdata = d, type = "density", log = TRUE), 
  +     predict(mmN, newdata = d, type = "density", log = TRUE))
  [1] TRUE
  > chk(predict(mm, newdata = d, type = "distribution", log = TRUE), 
  +     predict(mmN, newdata = d, type = "distribution", log = TRUE))
  [1] TRUE
  > 
  > chk(predict(mm, newdata = d, margins = 1:2, type = "density", log = TRUE), 
  +     predict(mmN, newdata = d, margins = 1:2, type = "density", log = TRUE))
  [1] TRUE
  > chk(predict(mm, newdata = d, margins = 1:2, type = "distribution", log = TRUE), 
  +     predict(mmN, newdata = d, margins = 1:2, type = "distribution", log = TRUE))
  [1] TRUE
  > chk(predict(mm, newdata = d, margins = 1:3, type = "density", log = TRUE), 
  +     predict(mmN, newdata = d, margins = 1:3, type = "density", log = TRUE))
  [1] TRUE
  > chk(predict(mm, newdata = d, margins = 1:3, type = "distribution", log = TRUE), 
  +     predict(mmN, newdata = d, margins = 1:3, type = "distribution", log = TRUE))
  [1] TRUE
  > 
  > chk(sapply(1:J, function(i) predict(mm, margins = i, newdata = d, type = "density", log = TRUE)),
  +     sapply(1:J, function(i) predict(mmN, margins = i, newdata = d, type = "density", log = TRUE)))
  [1] TRUE
  > 
  > ### check marginal predictions
  > m1 <- m[[1]]
  > m2 <- do.call("mmlt", m[-(3:4)])
  > m3 <- do.call("mmlt", m[-4])
  > 
  > ### we expect differences here
  > if (FALSE) {
  + chk(c(predict(m1, newdata = d, type = "density", log = TRUE)), 
  +     c(predict(mmN, newdata = d, margins = 1, type = "density", log = TRUE)))
  + chk(predict(m2, newdata = d, margins = 1:2, type = "density", log = TRUE), 
  +     predict(mmN, newdata = d, margins = 1:2, type = "density", log = TRUE))
  + chk(predict(m2, newdata = d, margins = 1:2, type = "distribution", log = TRUE), 
  +     predict(mmN, newdata = d, margins = 1:2, type = "distribution", log = TRUE))
  + chk(predict(m3, newdata = d, margins = 1:3, type = "density", log = TRUE), 
  +     predict(mmN, newdata = d, margins = 1:3, type = "density", log = TRUE))
  + chk(predict(m3, newdata = d, margins = 1:3, type = "distribution", log = TRUE), 
  +     predict(mmN, newdata = d, margins = 1:3, type = "distribution", log = TRUE))
  + }
  > 
  > ### marginal normal, implemented differently
  > for (j in 1:J) m[[j]]$todistr$name <- "CarlFriedrich"
  > mmN <- do.call("mmlt", m)
  > 
  > chk(logLik(mm), logLik(mmN))
  [1] TRUE
  > chk(c(logLik(mmN)), sum(predict(mmN, newdata = d, type = "density", log = TRUE)))
  [1] TRUE
  > 
  > cf1 <- do.call("c", lapply(m[1:J], function(x) coef(as.mlt(x))))
  > cf2 <- coef(mmN)[1:length(cf1)]
  > # cbind(cf1, cf2)
  > 
  > sd1 <- sqrt(do.call("c", lapply(m[1:J], function(x) diag(vcov(as.mlt(x))))))
  > sd2 <- sqrt(diag(vcov(mmN)))[1:length(sd1)]
  > 
  > # cbind(sd1, sd2)
  > vcov(mmN)["V1.x", "V4.x"]
  [1] 0.054
  > 
  > chk(as.array(coef(mm, type = "Lambda"))[,,1], 
  +     as.array(coef(mmN, type = "Lambda"))[,,1])
  [1] TRUE
  > chk(as.array(coef(mm, type = "Lambdainv"))[,,1], 
  +     as.array(coef(mmN, type = "Lambdainv"))[,,1])
  [1] TRUE
  > chk(as.array(coef(mm, type = "Sigma"))[,,1], 
  +     as.array(coef(mmN, type = "Sigma"))[,,1])
  [1] TRUE
  > chk(as.array(coef(mm, type = "Spearman"))[,,1], 
  +     as.array(coef(mmN, type = "Spearman"))[,,1])
  [1] TRUE
  > 
  > chk(predict(mm, newdata = d, type = "density", log = TRUE), 
  +     predict(mmN, newdata = d, type = "density", log = TRUE))
  [1] TRUE
  > chk(predict(mm, newdata = d, type = "distribution", log = TRUE), 
  +     predict(mmN, newdata = d, type = "distribution", log = TRUE))
  [1] TRUE
  > 
  > chk(predict(mm, newdata = d, margins = 1:2, type = "density", log = TRUE), 
  +     predict(mmN, newdata = d, margins = 1:2, type = "density", log = TRUE))
  [1] TRUE
  > chk(predict(mm, newdata = d, margins = 1:2, type = "distribution", log = TRUE), 
  +     predict(mmN, newdata = d, margins = 1:2, type = "distribution", log = TRUE))
  [1] TRUE
  > chk(predict(mm, newdata = d, margins = 1:3, type = "density", log = TRUE), 
  +     predict(mmN, newdata = d, margins = 1:3, type = "density", log = TRUE))
  [1] TRUE
  > chk(predict(mm, newdata = d, margins = 1:3, type = "distribution", log = TRUE), 
  +     predict(mmN, newdata = d, margins = 1:3, type = "distribution", log = TRUE))
  [1] TRUE
  > 
  > chk(sapply(1:J, function(i) predict(mm, margins = i, newdata = d, type = "density", log = TRUE)),
  +     sapply(1:J, function(i) predict(mmN, margins = i, newdata = d, type = "density", log = TRUE)))
  [1] TRUE
  > 
  > ### marginal Colr models
  > m <- lapply(un, function(i)
  +     Colr(as.formula(paste(i, "~ x")), data = d, bounds = c(0, 1), support = c(0, 1), order = OR))
  > m$data <- d
  > m$formula <- ~ 1
  > mmC <- do.call("mmlt", m)
  > 
  > chk(c(logLik(mmC)), sum(predict(mmC, newdata = d, type = "density", log = TRUE)))
  [1] TRUE
  > logLik(mmC)
  'log Lik.' 261 (df=18)
  > 
  > ### conditional models
  > m <- lapply(un, function(i)
  +     BoxCox(as.formula(paste(i, "~ x")), data = d, bounds = c(0, 1), support = c(0, 1), order = OR))
  > m$data <- d
  > m$formula <- ~ x
  > mm <- do.call("mmlt", m)
  > 
  > chk(c(logLik(mm)), sum(predict(mm, newdata = d, type = "density", log = TRUE)))
  [1] TRUE
  > L <- as.array(coef(mm, newdata = d, type = "Lambda"))[,,1]
  > chk(as.array(coef(mm, newdata = d, type = "Lambdainv"))[,,1], solve(L))
  [1] TRUE
  > chk(as.array(coef(mm, newdata = d, type = "Sigma"))[,,1], tcrossprod(solve(L)))
  [1] TRUE
  > chk(as.array(coef(mm, newdata = d, type = "Cor"))[,,1], cov2cor(tcrossprod(solve(L))))
  [1] TRUE
  > 
  > ### with marginal parameterisation
  > m$conditional <- FALSE
  > mmN <- do.call("mmlt", m)
  > 
  > chk(logLik(mm), logLik(mmN))
  [1] TRUE
  > chk(c(logLik(mmN)), sum(predict(mmN, newdata = d, type = "density", log = TRUE)))
  [1] TRUE
  > 
  > chk(as.array(coef(mm, newdata = d, type = "Lambda"))[,,1], 
  +     as.array(coef(mmN, newdata = d, type = "Lambda"))[,,1])
  [1] TRUE
  > chk(as.array(coef(mm, newdata = d, type = "Lambdainv"))[,,1], 
  +     as.array(coef(mmN, newdata = d, type = "Lambdainv"))[,,1])
  [1] TRUE
  > chk(as.array(coef(mm, newdata = d, type = "Sigma"))[,,1], 
  +     as.array(coef(mmN, newdata = d, type = "Sigma"))[,,1])
  [1] TRUE
  > chk(as.array(coef(mm, newdata = d, type = "Spearman"))[,,1], 
  +     as.array(coef(mmN, newdata = d, type = "Spearman"))[,,1])
  [1] TRUE
  > 
  > chk(predict(mm, newdata = d, type = "density", log = TRUE), 
  +     predict(mmN, newdata = d, type = "density", log = TRUE))
  [1] TRUE
  > chk(predict(mm, newdata = d, type = "distribution", log = TRUE), 
  +     predict(mmN, newdata = d, type = "distribution", log = TRUE))
  [1] TRUE
  > 
  > chk(predict(mm, newdata = d, margins = 1:2, type = "density", log = TRUE), 
  +     predict(mmN, newdata = d, margins = 1:2, type = "density", log = TRUE))
  [1] TRUE
  > chk(predict(mm, newdata = d, margins = 1:2, type = "distribution", log = TRUE), 
  +     predict(mmN, newdata = d, margins = 1:2, type = "distribution", log = TRUE))
  [1] TRUE
  > chk(predict(mm, newdata = d, margins = 1:3, type = "density", log = TRUE), 
  +     predict(mmN, newdata = d, margins = 1:3, type = "density", log = TRUE))
  [1] TRUE
  > chk(predict(mm, newdata = d, margins = 1:3, type = "distribution", log = TRUE), 
  +     predict(mmN, newdata = d, margins = 1:3, type = "distribution", log = TRUE))
  [1] TRUE
  > 
  > chk(sapply(1:J, function(i) predict(mm, margins = i, newdata = d, type = "density", log = TRUE)),
  +     sapply(1:J, function(i) predict(mmN, margins = i, newdata = d, type = "density", log = TRUE)))
  [1] TRUE
  > 
  > ### implemented differently
  > for (j in 1:J) m[[j]]$todistr$name <- "CarlFriedrich"
  > mmN <- do.call("mmlt", m)
  > 
  > chk(logLik(mm), logLik(mmN))
  [1] TRUE
  > chk(c(logLik(mmN)), sum(predict(mmN, newdata = d, type = "density", log = TRUE)))
  [1] TRUE
  > 
  > chk(as.array(coef(mm, newdata = d, type = "Lambda"))[,,1], 
  +     as.array(coef(mmN, newdata = d, type = "Lambda"))[,,1])
  [1] TRUE
  > chk(as.array(coef(mm, newdata = d, type = "Lambdainv"))[,,1], 
  +     as.array(coef(mmN, newdata = d, type = "Lambdainv"))[,,1])
  [1] TRUE
  > chk(as.array(coef(mm, newdata = d, type = "Sigma"))[,,1], 
  +     as.array(coef(mmN, newdata = d, type = "Sigma"))[,,1])
  [1] TRUE
  > chk(as.array(coef(mm, newdata = d, type = "Spearman"))[,,1], 
  +     as.array(coef(mmN, newdata = d, type = "Spearman"))[,,1])
  [1] TRUE
  > 
  > chk(predict(mm, newdata = d, type = "density", log = TRUE), 
  +     predict(mmN, newdata = d, type = "density", log = TRUE))
  [1] TRUE
  > chk(predict(mm, newdata = d, type = "distribution", log = TRUE), 
  +     predict(mmN, newdata = d, type = "distribution", log = TRUE))
  [1] TRUE
  > 
  > chk(predict(mm, newdata = d, margins = 1:2, type = "density", log = TRUE), 
  +     predict(mmN, newdata = d, margins = 1:2, type = "density", log = TRUE))
  [1] TRUE
  > chk(predict(mm, newdata = d, margins = 1:2, type = "distribution", log = TRUE), 
  +     predict(mmN, newdata = d, margins = 1:2, type = "distribution", log = TRUE))
  [1] TRUE
  > chk(predict(mm, newdata = d, margins = 1:3, type = "density", log = TRUE), 
  +     predict(mmN, newdata = d, margins = 1:3, type = "density", log = TRUE))
  [1] TRUE
  > chk(predict(mm, newdata = d, margins = 1:3, type = "distribution", log = TRUE), 
  +     predict(mmN, newdata = d, margins = 1:3, type = "distribution", log = TRUE))
  [1] TRUE
  > 
  > chk(sapply(1:J, function(i) predict(mm, margins = i, newdata = d, type = "density", log = TRUE)),
  +     sapply(1:J, function(i) predict(mmN, margins = i, newdata = d, type = "density", log = TRUE)))
  [1] TRUE
  > 
  > ### conditional Colr
  > m <- lapply(un, function(i)
  +     Colr(as.formula(paste(i, "~ x")), data = d, bounds = c(0, 1), support = c(0, 1), order = OR))
  > m$data <- d
  > m$formula <- ~ x
  > mmC <- do.call("mmlt", m)
  > 
  > chk(c(logLik(mmC)), sum(predict(mmC, newdata = d, type = "density", log = TRUE)))
  [1] TRUE
  > logLik(mmC)
  'log Lik.' 265 (df=24)
  > 
  > ##### FIRST SCENARIO: CONSTANT LAMBDA #####
  > set.seed(290875)
  > ll <- numeric(50)
  > p <- 3
  > X <- matrix(runif(N * p), ncol = p)
  > m1 <- 1 + X %*% c(2, 1, 1)
  > m2 <- 1 + X %*% c(1, 2, 1)
  > lb <- (off <- 0.5) + X %*% (cf <- c(0, 2, 0))
  > d <- data.frame(X)
  > Y <- matrix(NA, nrow = N, ncol = 2)
  > colnames(Y) <- c("Y1", "Y2")
  > 
  > cr <- numeric(N)
  > for (i in 1:N) {
  +   Si <- diag(2)
  +   Si[1,2] <- Si[2,1] <- .5
  +   cr[i] <- cov2cor(Si)[2,1]
  +   
  +   Y[i,] <- rmvnorm(1, mean = c(m1[i], m2[i]), sigma = Si)
  + }
  > 
  > 
  > ##### only BoxCox margins: ##### 
  > d <- cbind(d, Y)
  > b1 <- as.mlt(Lm(Y1 ~ X1 + X2 + X3, data = d))
  > b2 <- as.mlt(Lm(Y2 ~ X1 + X2 + X3, data = d))
  > 
  > ## constant correlations. expect identical logliks and lambda parameters
  > mm01 <- mmlt(b1, b2, formula = ~ 1, data = d)
  > mm02 <- mmlt(b2, b1, formula = ~ 1, data = d)
  > 
  > chk(logLik(mm01), logLik(mm02))
  [1] TRUE
  > 
  > chk(c(coef(mm01)["Y2.Y1.(Intercept)"]), 
  +     c(coef(mm02)["Y1.Y2.(Intercept)"]))
  [1] TRUE
  > 
  > ## checking gradients
  > chk(c(numDeriv::grad(mm01$ll, mm02$par)),
  +     c(mm01$sc(mm02$par)))
  [1] TRUE
  > 
  > ## predicting marginal distributions and comparing across models with constant lambda
  > chk(predict(mm01, newdata = d[1:5,], q = -2:2, 
  +         margins = 1, type = "distribution"),
  +     predict(mm02, newdata = d[1:5,], q = -2:2, 
  +         margins = 2, type = "distribution"))
  [1] TRUE
  > 
  > ## expect correlations to be the same for the model with constant lambdas
  > chk(c(coef(mm01, newdata = d[1:5,], type = "Cor")), 
  +     c(coef(mm02, newdata = d[1:5,], type = "Cor")))
  [1] TRUE
  > 
  > 
  > ##### mix of BoxCox and Colr margins: ##### 
  > d$Y1 <- (d$Y1 - min(d$Y1))/(max(d$Y1) - min(d$Y1))
  > 
  > b1 <- as.mlt(Colr(Y1 ~ X1 + X2 + X3, data = d, order = OR))
  > b2 <- as.mlt(Lm(Y2 ~ X1 + X2 + X3, data = d))
  > 
  > mm01 <- mmlt(b1, b2, formula = ~ 1, data = d)
  > mm02 <- mmlt(b2, b1, formula = ~ 1, data = d)
  > 
  > chk(logLik(mm01), logLik(mm02))
  [1] TRUE
  > 
  > ### model is marginal, so expect same marginal coeff
  > cf1 <- coef(mm01)
  > cf1 <- cf1[-length(cf1)]
  > cf2 <- coef(mm02)
  > cf2 <- cf2[names(cf1)]
  > chk(cf1, cf2)
  [1] TRUE
  > 
  > ## checking gradient
  > chk(c(numDeriv::grad(mm01$ll, coef(mm01))), c(mm01$sc(coef(mm01))))
  [1] TRUE
  > chk(c(numDeriv::grad(mm02$ll, coef(mm02))), c(mm02$sc(coef(mm02)))) 
  [1] TRUE
  > 
  > ##### SECOND SCENARIO: COVARIATE DEPENDENT LAMBDA #####
  > set.seed(290875)
  > ll <- numeric(50)
  > X <- matrix(runif(N * p), ncol = p)
  > m1 <- 1 + X %*% c(2, 1, 1)
  > m2 <- 1 + X %*% c(1, 2, 1)
  > lb <- (off <- 0.5) + X %*% (cf <- c(0, 2, 0))
  > d <- data.frame(X)
  > Y <- matrix(NA, nrow = N, ncol = 2)
  > colnames(Y) <- c("Y1", "Y2")
  > 
  > cr <- numeric(N)
  > for (i in 1:N) {
  +   L <- diag(2)
  +   L[2,1] <- lb[i]
  +   Si <- solve(L) %*% t(solve(L))
  +   cr[i] <- cov2cor(Si)[2,1]
  +   
  +   Y[i,] <- rmvnorm(1, mean = c(m1[i], m2[i]), sigma = Si)
  + }
  > 
  > 
  > ##### only BoxCox margins: ##### 
  > d <- cbind(d, Y)
  > b1 <- as.mlt(Lm(Y1 ~ X1 + X2 + X3, data = d))
  > b2 <- as.mlt(Lm(Y2 ~ X1 + X2 + X3, data = d))
  > 
  > ## constant correlations. expect identical logliks and lambda parameters
  > mm01 <- mmlt(b1, b2, formula = ~ 1, data = d)
  > mm02 <- mmlt(b2, b1, formula = ~ 1, data = d)
  > 
  > chk(logLik(mm01), logLik(mm02))
  [1] TRUE
  > 
  > ## checking gradients
  > chk(c(numDeriv::grad(mm01$ll, mm02$par)),c(mm01$sc(mm02$par)))
  [1] TRUE
  > 
  > ## x-dependent correlations. expect slightly different logliks when
  > ## conditional = TRUE
  > mm1 <- mmlt(b1, b2, formula = ~ X1 + X2 + X3, data = d, conditional = TRUE)
  Warning message:
  In mmlt(b1, b2, formula = ~X1 + X2 + X3, data = d, conditional = TRUE) :
    Conditional models with covariate-dependent correlations are order-dependent
  > mm2 <- mmlt(b2, b1, formula = ~ X1 + X2 + X3, data = d, conditional = TRUE)
  Warning message:
  In mmlt(b2, b1, formula = ~X1 + X2 + X3, data = d, conditional = TRUE) :
    Conditional models with covariate-dependent correlations are order-dependent
  > 
  > logLik(mm1)
  'log Lik.' -266 (df=14)
  > logLik(mm2)
  'log Lik.' -282 (df=14)
  > 
  > ### BUT: identical models when conditional = FALSE
  > mm1 <- mmlt(b1, b2, formula = ~ X1 + X2 + X3, data = d, conditional = FALSE)
  > mm2 <- mmlt(b2, b1, formula = ~ X1 + X2 + X3, data = d, conditional = FALSE)
  > 
  > chk(logLik(mm1), logLik(mm2))
  [1] TRUE
  > 
  > ## predicting marginal distributions and comparing across models with constant lambda
  > x <- 0:4 / 4
  > nd <- expand.grid(X1 = x, X2 = x, X3 = x)
  > chk(predict(mm01, newdata = nd[1:5,], q = -2:2, 
  +         margins = 1, type = "distribution"),
  +     predict(mm02, newdata = nd[1:5,], q = -2:2, 
  +         margins = 2, type = "distribution"))
  [1] TRUE
  > 
  > ## predicting marginal distributions and comparing across models with
  > ## x-dependent lambda and conditional = FALSE
  > chk(predict(mm1, newdata = nd[1:5,], q = -2:2, 
  +         margins = 1, type = "distribution"),
  +     predict(mm2, newdata = nd[1:5,], q = -2:2, 
  +         margins = 2, type = "distribution"))
  [1] TRUE
  > 
  > ## expect correlations to be the same for the model with constant lambdas
  > chk(c(coef(mm01, newdata = nd[1:5,], type = "Cor")), 
  +     c(coef(mm02, newdata = nd[1:5,], type = "Cor")))
  [1] TRUE
  > 
  > ## correlations for models with x-dependent lambda
  > chk(c(coef(mm1, newdata = nd[1:5,], type = "Cor")),
  +     c(coef(mm2, newdata = nd[1:5,], type = "Cor")))
  [1] TRUE
  > 
  > 
  > ##### mix of BoxCox and Colr margins: ##### 
  > d$Y1 <- (d$Y1 - min(d$Y1))/(max(d$Y1) - min(d$Y1))
  > 
  > b1 <- as.mlt(Colr(Y1 ~ X1 + X2 + X3, data = d, order = OR))
  > b2 <- as.mlt(Lm(Y2 ~ X1 + X2 + X3, data = d))
  > 
  > mm01 <- mmlt(b1, b2, formula = ~ 1, data = d)
  > mm02 <- mmlt(b2, b1, formula = ~ 1, data = d)
  > 
  > chk(logLik(mm01), logLik(mm02))
  [1] TRUE
  > 
  > coef(b1)
  Bs1(Y1) Bs2(Y1)      X1      X2      X3 
     0.52    6.09   -4.35   -1.43   -0.66 
  > coef(b2)
  (Intercept)          Y2          X1          X2          X3 
        -0.18        0.53        0.34        1.06        1.20 
  > coef(mm01)
         Y1.Bs1(Y1)        Y1.Bs2(Y1)             Y1.X1             Y1.X2 
               0.50              6.08             -4.25             -1.44 
              Y1.X3    Y2.(Intercept)             Y2.Y2             Y2.X1 
              -0.66             -0.17              0.52              0.41 
              Y2.X2             Y2.X3 Y2.Y1.(Intercept) 
               1.01              1.17              1.51 
  > coef(mm02)
     Y2.(Intercept)             Y2.Y2             Y2.X1             Y2.X2 
              -0.17              0.52              0.41              1.01 
              Y2.X3        Y1.Bs1(Y1)        Y1.Bs2(Y1)             Y1.X1 
               1.17              0.50              6.08             -4.25 
              Y1.X2             Y1.X3 Y1.Y2.(Intercept) 
              -1.44             -0.66              1.51 
  > # remember that: lb <- (off <- 0.5) + X %*% (cf <- c(0, 2, 0))
  > 
  > 
  > ## checking gradient
  > chk(c(numDeriv::grad(mm01$ll, coef(mm01))),c(mm01$sc(coef(mm01))))
  [1] TRUE
  > chk(c(numDeriv::grad(mm02$ll, coef(mm02))),c(mm02$sc(coef(mm02))))
  [1] TRUE
  > 
  > ## covariate-dependent Lambda
  > mm1 <- mmlt(b1, b2, formula = ~ X1 + X2 + X3, data = d)
  > mm2 <- mmlt(b2, b1, formula = ~ X1 + X2 + X3, data = d)
  > logLik(mm1)
  'log Lik.' -88 (df=14)
  > logLik(mm2)
  'log Lik.' -88 (df=14)
  > 
  > coef(b1)
  Bs1(Y1) Bs2(Y1)      X1      X2      X3 
     0.52    6.09   -4.35   -1.43   -0.66 
  > coef(b2)
  (Intercept)          Y2          X1          X2          X3 
        -0.18        0.53        0.34        1.06        1.20 
  > coef(mm1)
         Y1.Bs1(Y1)        Y1.Bs2(Y1)             Y1.X1             Y1.X2 
               0.43              6.29             -4.43             -1.28 
              Y1.X3    Y2.(Intercept)             Y2.Y2             Y2.X1 
              -0.81             -0.10              0.50              0.32 
              Y2.X2             Y2.X3 Y2.Y1.(Intercept)          Y2.Y1.X1 
               1.13              1.11              1.92             -1.13 
           Y2.Y1.X2          Y2.Y1.X3 
               1.76             -1.09 
  > coef(mm2)
     Y2.(Intercept)             Y2.Y2             Y2.X1             Y2.X2 
              -0.10              0.50              0.32              1.13 
              Y2.X3        Y1.Bs1(Y1)        Y1.Bs2(Y1)             Y1.X1 
               1.11              0.43              6.29             -4.43 
              Y1.X2             Y1.X3 Y1.Y2.(Intercept)          Y1.Y2.X1 
              -1.28             -0.81              1.92             -1.13 
           Y1.Y2.X2          Y1.Y2.X3 
               1.76             -1.09 
  > # remember that: lb <- (off <- 0.5) + X %*% (cf <- c(0, 2, 0))
  > 
  > ## checking gradient for diag = TRUE
  > chk(c(numDeriv::grad(mm1$ll, coef(mm1))),c(mm1$sc(coef(mm1))))
  [1] TRUE
  > chk(c(numDeriv::grad(mm2$ll, coef(mm2))),c(mm2$sc(coef(mm2))))
  [1] TRUE
  > 
  > ### very simple checks with marginal Lm models
  > set.seed(290875)
  > J <- 4
  > S <- cov2cor(tcrossprod(matrix(runif(J * J), ncol = J)))
  > x <- matrix(runif(N*2), ncol = 2)
  > 
  > y <- x %*% matrix(c(1, -1, -.5, .5, -.2, .2, .3, -.3), nrow = 2) + rmvnorm(N, sigma = S)
  > d <- data.frame(y = y, x = x)
  > 
  > m1 <- Lm(y.1 ~ x.1 + x.2, data = d)
  > m2 <- Lm(y.2 ~ x.1 + x.2, data = d)
  > m3 <- Lm(y.3 ~ x.1 + x.2, data = d)
  > m4 <- Lm(y.4 ~ x.1 + x.2, data = d)
  > 
  > ## simple formula
  > mc01 <- mmlt(m1, m2, m3, m4, formula = ~ 1, data = d, conditional = FALSE)
  > 
  > cf <- coef(mc01)
  > vr <- diag(vcov(mc01))
  > i <- grep("x", names(cf))
  > 
  > ### same results
  > ret <- cbind(c(coef(m1), coef(m2), coef(m3), coef(m4)),
  +              cf[i],
  +              c(diag(vcov(m1)), diag(vcov(m2)), diag(vcov(m3)), diag(vcov(m4))),
  +              vr[i])
  > ret
        [,1]   [,2] [,3] [,4]
  x.1  1.036  1.036 0.13 0.13
  x.2 -1.184 -1.184 0.14 0.14
  x.1 -0.702 -0.702 0.12 0.12
  x.2  0.743  0.743 0.13 0.13
  x.1 -0.585 -0.585 0.12 0.12
  x.2  0.156  0.156 0.13 0.13
  x.1  0.038  0.038 0.12 0.12
  x.2 -0.420 -0.420 0.13 0.13
  > 
  > vc <- vcov(mc01)
  > i <- grep("x.1", colnames(vc))
  > vc[i,i]
          y.1.x.1 y.2.x.1 y.3.x.1 y.4.x.1
  y.1.x.1   0.127   0.029   0.085   0.072
  y.2.x.1   0.029   0.124   0.087   0.070
  y.3.x.1   0.085   0.087   0.124   0.111
  y.4.x.1   0.072   0.070   0.111   0.122
  > 
  > summary(g1 <- glht(mmm(m1 = as.mlt(m1), m2 = as.mlt(m2), m3 = as.mlt(m3), m4 = as.mlt(m4)), mlf("x.1 = 0")))
  
  	 Simultaneous Tests for General Linear Hypotheses
  
  Linear Hypotheses:
               Estimate Std. Error z value Pr(>|z|)  
  m1: x.1 == 0   1.0365     0.3569    2.90    0.012 *
  m2: x.1 == 0  -0.7019     0.3528   -1.99    0.126  
  m3: x.1 == 0  -0.5846     0.3517   -1.66    0.241  
  m4: x.1 == 0   0.0383     0.3493    0.11    1.000  
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  (Adjusted p values reported -- single-step method)
  
  > 
  > summary(g2 <- glht(mc01, c("y.1.x.1 = 0", "y.2.x.1 = 0", "y.3.x.1 = 0", "y.4.x.1 = 0")))
  
  	 Simultaneous Tests for General Linear Hypotheses
  
  Fit: mmlt(m1, m2, m3, m4, formula = ~1, data = d, conditional = FALSE)
  
  Linear Hypotheses:
               Estimate Std. Error z value Pr(>|z|)  
  y.1.x.1 == 0   1.0365     0.3569    2.90    0.011 *
  y.2.x.1 == 0  -0.7019     0.3528   -1.99    0.126  
  y.3.x.1 == 0  -0.5846     0.3517   -1.66    0.241  
  y.4.x.1 == 0   0.0383     0.3493    0.11    1.000  
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  (Adjusted p values reported -- single-step method)
  
  > 
  > vcov(g1)
          m1: x.1 m2: x.1 m3: x.1 m4: x.1
  m1: x.1   0.127   0.029   0.086   0.078
  m2: x.1   0.029   0.124   0.086   0.065
  m3: x.1   0.086   0.086   0.124   0.111
  m4: x.1   0.078   0.065   0.111   0.122
  > vcov(g2)
          y.1.x.1 y.2.x.1 y.3.x.1 y.4.x.1
  y.1.x.1   0.127   0.029   0.085   0.072
  y.2.x.1   0.029   0.124   0.087   0.070
  y.3.x.1   0.085   0.087   0.124   0.111
  y.4.x.1   0.072   0.070   0.111   0.122
  > 
  > 
  > #### check density
  > Shat <- as.array(coef(mc01, type = "Cor"))[,,1]
  > 
  > int <- cf[paste("y", 1:J, "(Intercept)", sep = ".")]
  > fct <- cf[paste("y", 1:J, "y", 1:J, sep = ".")]
  > 
  > d1 <- sapply(1:N, function(i) dmvnorm(int + fct * y[i,], mean = x[i,,drop = FALSE] %*% matrix(ret[,1], nrow = 2), sigma = Shat, log = TRUE))
  > d2 <- predict(mc01, newdata = d, type = "density", log = TRUE)
  > 
  > chk(c(d1), c(d2))
  [1] "Mean relative difference: 0.076"
  > 
  > chk(c(logLik(mmlt(m1, m2, m3, m4, formula = ~ 1, data = d))),
  +     c(logLik(mc01)))
  [1] TRUE
  > chk(c(logLik(mc01)),
  +     sum(d2))
  [1] TRUE
  > 
  > ### check if newdata works in logLik
  > chk(logLik(mc01), logLik(mc01, newdata = d))
  [1] TRUE
  > 
  > ### check user interface
  > dgp <- function(N = 100, J = 5, lambda = 0.5773503) {
  + 
  +     Jp <- J * (J - 1) / 2
  +     X <- c(lambda, rep(0, Jp - 1))
  +     L <- ltMatrices(X)
  +     L <- standardize(invchol = L)
  +     Z <- matrix(rnorm(N * J), ncol = N)
  +     ret <- solve(L, Z)
  +     ret <- as.data.frame(t(ret))
  +     colnames(ret) <- paste0("Y", 1:J)
  +     ret
  + }
  > 
  > N <- 100
  > J <- 5
  > 
  > Y <- dgp(N = N, J = J)
  > 
  > m0 <- lapply(colnames(Y)[1:J], function(v) {
  +     fm <- as.formula(paste(v, " ~ 1"))
  +     BoxCox(fm, data = Y, order = OR)
  + })
  > 
  > TF <- c(TRUE, FALSE)
  > args <- expand.grid(scale  = TF, domargins = TF, dofit = TF, theta = TF, fixed = TF, conditional = TF)
  > args <- subset(args, !(conditional & !domargins))
  > 
  > m0$data <- Y
  > m0$conditional <- TRUE
  > m1 <- do.call("mmlt", m0)
  > theta <- coef(m1)
  > CR <- coef(m1, type = "Cor")
  > 
  > fx <- c("Y5.Y3.(Intercept)" = 0, "Y5.Y4.(Intercept)" = 0)
  > 
  > for (i in 1:nrow(args)) {
  +     print(i)
  +     m0$scale <- args$scale[i]
  +     m0$dofit <- args$dofit[i]
  +     m0$domargins <- args$domargins[i]
  +     m0$conditional <- args$conditional[i]
  +     m0$theta <- NULL
  +     if (args$theta[i])
  +         m0$theta <- theta
  +     m0$fixed <- NULL
  +     if (args$fixed[i])
  +         m0$fixed <- fx
  +     m1 <- try(do.call("mmlt", m0))
  +     if (inherits(m1, "mmlt")) {
  +         print(logLik(m1))
  +         print(isTRUE(chk(coef(m1, type = "Cor"), CR)))
  +     } else {
  +         print(m1$ll(theta))
  +     }
  + }
  [1] 1
  'log Lik.' -678 (df=18)
  [1] FALSE
  [1] 2
  'log Lik.' -678 (df=18)
  [1] FALSE
  [1] 3
  'log Lik.' -678 (df=18)
  [1] FALSE
  [1] 4
  'log Lik.' -678 (df=18)
  [1] FALSE
  [1] 5
  'log Lik.' -678 (df=18)
  [1] FALSE
  [1] 6
  'log Lik.' -678 (df=18)
  [1] FALSE
  [1] 7
  [1] 678
  [1] 8
  [1] 678
  [1] 9
  'log Lik.' -678 (df=20)
  [1] TRUE
  [1] 10
  'log Lik.' -678 (df=20)
  [1] TRUE
  [1] 11
  'log Lik.' -678 (df=20)
  [1] TRUE
  [1] 12
  'log Lik.' -678 (df=20)
  [1] TRUE
  [1] 13
  'log Lik.' -678 (df=20)
  [1] TRUE
  [1] 14
  'log Lik.' -678 (df=20)
  [1] TRUE
  [1] 15
  [1] 678
  [1] 16
  [1] 678
  [1] 17
  'log Lik.' -678 (df=18)
  [1] FALSE
  [1] 18
  'log Lik.' -678 (df=18)
  [1] FALSE
  [1] 19
  'log Lik.' -678 (df=18)
  [1] FALSE
  [1] 20
  'log Lik.' -678 (df=18)
  [1] FALSE
  [1] 21
  'log Lik.' -681 (df=18)
  [1] FALSE
  [1] 22
  'log Lik.' -681 (df=18)
  [1] FALSE
  [1] 23
  'log Lik.' -678 (df=18)
  [1] FALSE
  [1] 24
  'log Lik.' -678 (df=18)
  [1] FALSE
  [1] 25
  'log Lik.' -678 (df=18)
  [1] FALSE
  [1] 26
  'log Lik.' -678 (df=18)
  [1] FALSE
  [1] 27
  'log Lik.' -678 (df=18)
  [1] FALSE
  [1] 28
  'log Lik.' -678 (df=18)
  [1] FALSE
  [1] 29
  [1] 681
  [1] 30
  [1] 681
  [1] 31
  'log Lik.' -678 (df=18)
  [1] FALSE
  [1] 32
  'log Lik.' -678 (df=18)
  [1] FALSE
  [1] 33
  'log Lik.' -678 (df=20)
  [1] TRUE
  [1] 34
  'log Lik.' -678 (df=20)
  [1] TRUE
  [1] 35
  'log Lik.' -678 (df=20)
  [1] TRUE
  [1] 36
  'log Lik.' -678 (df=20)
  [1] TRUE
  [1] 37
  'log Lik.' -681 (df=20)
  [1] TRUE
  [1] 38
  'log Lik.' -681 (df=20)
  [1] TRUE
  [1] 39
  'log Lik.' -678 (df=20)
  [1] TRUE
  [1] 40
  'log Lik.' -678 (df=20)
  [1] TRUE
  [1] 41
  'log Lik.' -678 (df=20)
  [1] TRUE
  [1] 42
  'log Lik.' -678 (df=20)
  [1] TRUE
  [1] 43
  'log Lik.' -678 (df=20)
  [1] TRUE
  [1] 44
  'log Lik.' -678 (df=20)
  [1] TRUE
  [1] 45
  [1] 681
  [1] 46
  [1] 681
  [1] 47
  'log Lik.' -678 (df=20)
  [1] TRUE
  [1] 48
  'log Lik.' -678 (df=20)
  [1] TRUE
  > 
  > ### names
  > J <- 5
  > N <- 50
  > df <- as.data.frame(matrix(rnorm(J * N), ncol = J))
  > colnames(df) <- paste0("X", 1:J)
  > 
  > mltargs <- lapply(1:ncol(df), function(j) {
  +   fm <- as.formula(paste0("X", j, "~1"))
  +   BoxCox(fm, data = df, order = OR)
  + })
  > mltargs$data <- df
  > 
  > fx <- c("X3.X1.(Intercept)" = 0, "X4.X1.(Intercept)" = 0, "X5.X1.(Intercept)" = 0,
  +         "X3.X2.(Intercept)" = 0, "X4.X2.(Intercept)" = 0, "X5.X2.(Intercept)" = 0,
  +         "X5.X4.(Intercept)" = 0)
  > tmp <- do.call("mmlt", c(mltargs, list(fixed = fx)))
  > 
  > mltargs$conditional <- TRUE
  > tmp <- do.call("mmlt", c(mltargs, list(fixed = fx)))
  > 
  > cf <- coef(tmp)
  > cf <- cf[grep("Intercept", names(cf))]
  > names(cf) <- substr(names(cf), 1, 5)
  > chk(unclass(coef(tmp, type = "Lambda"))[names(cf),], cf)
  [1] TRUE
  > 
  > ### check discrete models
  > J <- 2
  > N <- 100
  > S <- cov2cor(tcrossprod(matrix(runif(J * J), ncol = J)))
  > y <- rmvnorm(N, sigma = S)
  > u <- as.data.frame(plogis(y))
  > x <- runif(N)
  > d <- cbind(u, x)
  > un <- colnames(d)[1:J]
  > d[1:J] <- lapply(d[1:J], function(x) 
  +     cut(x, breaks = c(-Inf, quantile(x, prob = 1:3 / 4), Inf), ordered_result = TRUE))
  > 
  > m <- lapply(un, function(i)
  +     Polr(as.formula(paste(i, "~ x")), data = d, method = "probit"))
  > m$data <- d
  > m$formula <- ~ 1
  > m$args <- list(seed = 1, M = 100)
  > mm <- do.call("mmlt", m)
  > 
  > L <- as.array(coef(mm, type = "Lambda"))[,,1]
  > chk(as.array(coef(mm, type = "Lambdainv"))[,,1], solve(L))
  [1] TRUE
  > chk(as.array(coef(mm, type = "Sigma"))[,,1], tcrossprod(solve(L)))
  [1] TRUE
  > chk(as.array(coef(mm, type = "Cor"))[,,1], cov2cor(tcrossprod(solve(L))))
  [1] TRUE
  > chk(colSums(estfun(mm)), mm$score(coef(mm, type = "all")))
  [1] TRUE
  > 
  > for (j in 1:J) m[[j]]$todistr$name <- "CarlFriedrich"
  > 
  > mmN <- do.call("mmlt", m)
  > 
  > chk(logLik(mm), logLik(mmN))
  [1] TRUE
  > chk(coef(mm), coef(mmN))
  [1] TRUE
  > chk(diag(vcov(mm)), diag(vcov(mmN)))
  [1] TRUE
  > chk(as.array(coef(mm, type = "Lambda"))[,,1], 
  +     as.array(coef(mmN, type = "Lambda"))[,,1])
  [1] TRUE
  > chk(as.array(coef(mm, type = "Lambdainv"))[,,1], 
  +     as.array(coef(mmN, type = "Lambdainv"))[,,1])
  [1] TRUE
  > chk(as.array(coef(mm, type = "Sigma"))[,,1], 
  +     as.array(coef(mmN, type = "Sigma"))[,,1])
  [1] TRUE
  > chk(as.array(coef(mm, type = "Spearman"))[,,1], 
  +     as.array(coef(mmN, type = "Spearman"))[,,1])
  [1] TRUE
  > 
  > 
  > proc.time()
     user  system elapsed 
    109.8     2.3   112.1 
Running the tests in 'tests/mmlt_ll_sc.R' failed.
Complete output:
  > 
  > library("tram")
  Loading required package: mlt
  Loading required package: basefun
  Loading required package: variables
  Loading required package: mvtnorm
  > set.seed(29)
  > 
  > chk <- function(...) 
  +     stopifnot(isTRUE(all.equal(..., tol = 1e-6, check.attributes = FALSE)))
  > 
  > thischeck <- expression({
  +   ltM <- function(x) ltMatrices(x, diag = FALSE, byrow = TRUE)
  +   ltD <- function(x) ltMatrices(x, diag = TRUE, byrow = TRUE)
  +   prm <- matrix(runif(J * (J - 1) / 2 * N), ncol = N)
  +   L <- ltM(prm)
  + 
  +   obs <- matrix(rnorm(J * N), ncol = N)
  +   lwr <- -2 -abs(obs)
  +   upr <- 2 + abs(obs)
  + 
  +   fun <- tram:::.ll(c(J, 0), standardize = FALSE)
  +   sum(fun$logLik(obs, L))
  +   s <- fun$score(obs, L)
  + 
  +   f <- function(obs = obs, L = L) 
  +     sum(fun$logLik(obs, ltM(L)))
  +   S <- matrix(grad(f, unclass(L), obs = obs), ncol = N)
  +   chk(S, Lower_tri(s$Lambda))
  +   S <- matrix(grad(f, obs, L = L), ncol = N)
  +   chk(S, s$obs)
  + 
  +   fun <- tram:::.ll(c(J, 0), standardize = TRUE)
  +   sum(fun$logLik(obs, L))
  + 
  +   LD <- invcholD(L)
  +   sum(-colSums(Mult(LD, obs)^2 - log(diagonals(LD)^2)))
  + 
  +   C <- solve(L)
  +   CCt <- Tcrossprod(C, diag_only = TRUE)
  +   LD <- invcholD(L, D = sqrt(CCt))
  +   sum(-colSums(Mult(LD, obs)^2 - log(CCt)))
  + 
  +   f <- function(obs, a) {
  +     LD <- ltMatrices(a, diag = TRUE, byrow = TRUE)
  +     sum(-colSums(Mult(LD, obs)^2 - log(diagonals(LD)^2)))
  +   }
  +   f(obs = obs, LD)
  + 
  +   sLD <- ltMatrices(matrix(grad(f, unclass(LD), obs = obs), ncol = N), 
  +                     diag = TRUE, byrow = TRUE)
  +   sobs <- matrix(grad(f, obs, a = unclass(LD)), ncol = N)
  + 
  +   sLDfun <- function(obs, LD) {
  +     cJ <- dim(LD)[2L]
  +     Y <- matrix(obs, byrow = TRUE, nrow = cJ, ncol = N * cJ)
  +     tmp <- -2 * Mult(LD, Mult(LD, obs), transpose = TRUE)
  +     ret <- - 2 * matrix(Mult(LD, obs)[, rep(1:N, each = cJ)] * Y, ncol = N)
  +     M <- matrix(1:(cJ^2), nrow = cJ, byrow = FALSE)
  +     ret <- ltMatrices(ret[M[lower.tri(M, diag = TRUE)],,drop = FALSE], 
  +                           diag = TRUE, byrow = FALSE)
  +     ret <- ltMatrices(ret, 
  +                       diag = TRUE, byrow = TRUE)
  +     diagonals(ret) <- diagonals(ret) + 2 / diagonals(LD)
  +     return(list(Lambda = ret, obs = tmp))
  +   }
  + 
  +   s <- sLDfun(obs, LD)
  +   chk(s$Lambda, sLD)
  +   chk(s$obs, sobs)
  + 
  +   s <- fun$score(obs, L)
  + 
  +   f <- function(obs = obs, L = L) 
  +     sum(fun$logLik(obs, ltM(L)))
  +   S <- matrix(grad(f, unclass(L), obs = obs), ncol = N)
  +   chk(S, Lower_tri(s$Lambda))
  +   S <- matrix(grad(f, obs, L = L), ncol = N)
  +   chk(S, s$obs)
  + 
  +   w <- matrix(runif((J - 1) * M), ncol = M)
  +   fun <- tram:::.ll(c(0, J), standardize = FALSE, list(w = w))
  +   sum(fun$logLik(lower = lwr, upper = upr, Lambda = L))
  +   s <- fun$score(lower = lwr, upper = upr, Lambda = L)
  + 
  +   f <- function(lwr = lwr, upr = upr, L = L) 
  +     sum(fun$logLik(lower = lwr, upper = upr, Lambda = ltM(L)))
  +   S <- grad(f, unclass(L), lwr = lwr, upr = upr)
  +   chk(S, c(Lower_tri(s$Lambda)))
  + 
  +   S <- matrix(grad(f, lwr, upr = upr, L = L), ncol = N)
  +   chk(S, s$lower)
  +   S <- matrix(grad(f, upr, lwr = lwr, L = L), ncol = N)
  +   chk(S, s$upper)
  + 
  +   fun <- tram:::.ll(c(0, J), standardize = TRUE, list(w = w))
  +   sum(fun$logLik(lower = lwr, upper = upr, Lambda = L))
  +   s <- fun$score(lower = lwr, upper = upr, Lambda = L)
  + 
  +   f <- function(lwr = lwr, upr = upr, L = L) 
  +     sum(fun$logLik(lower = lwr, upper = upr, Lambda = ltM(L)))
  +   S <- matrix(grad(f, unclass(L), lwr = lwr, upr = upr), ncol = N)
  +   chk(S, Lower_tri(s$Lambda))
  +   S <- matrix(grad(f, lwr, upr = upr, L = L), ncol = N)
  +   chk(S, s$lower)
  +   S <- matrix(grad(f, upr, lwr = lwr, L = L), ncol = N)
  +   chk(S, s$upper)
  + 
  +   w <- matrix(runif((dJ - 1) * M), ncol = M)
  +   fun <- tram:::.ll(c(cJ, dJ), standardize = FALSE, list(w = w))
  +   sum(fun$logLik(obs[1:cJ,,drop = FALSE], lwr[-(1:cJ),,drop = FALSE], upr[-(1:cJ),,drop = FALSE], L))
  +   s <- fun$score(obs[1:cJ,,drop = FALSE], lwr[-(1:cJ),,drop = FALSE], upr[-(1:cJ),,drop = FALSE], L)
  +   f <- function(obs = obs[1:cJ,,drop = FALSE], lwr = lwr[-(1:cJ),,drop = FALSE], upr = upr[-(1:cJ),,drop = FALSE], L = L) 
  +     sum(fun$logLik(obs, lwr, upr, ltM(L)))
  +   S <- matrix(grad(f, unclass(L), obs = obs[1:cJ,,drop = FALSE], lwr = lwr[-(1:cJ),,drop = FALSE], 
  +                  upr = upr[-(1:cJ),,drop = FALSE]), ncol = N)
  +   chk(S, Lower_tri(s$Lambda))
  +   S <- matrix(grad(f, obs[1:cJ,,drop = FALSE], lwr = lwr[-(1:cJ),,drop = FALSE], 
  +                     upr = upr[-(1:cJ),,drop = FALSE], L = L), ncol = N)
  +   chk(S, s$obs)
  +   S <- matrix(grad(f, lwr[-(1:cJ),,drop = FALSE], obs = obs[1:cJ,,drop = FALSE], 
  +                     upr = upr[-(1:cJ),,drop = FALSE], L = L), ncol = N)
  +   chk(S, s$lower)
  +   S <- matrix(grad(f, upr[-(1:cJ),,drop = FALSE], obs = obs[1:cJ,,drop = FALSE], 
  +                     lwr = lwr[-(1:cJ),,drop = FALSE], L = L), ncol = N)
  +   chk(S, s$upper)
  + 
  +   fun <- tram:::.ll(c(cJ, dJ), standardize = TRUE, list(w = w))
  +   sum(fun$logLik(obs[1:cJ,,drop = FALSE], lwr[-(1:cJ),,drop = FALSE], upr[-(1:cJ),,drop = FALSE], L))
  +   s <- fun$score(obs[1:cJ,,drop = FALSE], lwr[-(1:cJ),,drop = FALSE], 
  +                  upr[-(1:cJ),,drop = FALSE], L)
  + 
  +   f <- function(obs = obs[1:cJ,,drop = FALSE], lwr = lwr[-(1:cJ),,drop = FALSE], upr = upr[-(1:cJ),,drop = FALSE], L = L) 
  +     sum(fun$logLik(obs, lwr, upr, ltM(L)))
  +   S <- matrix(grad(f, unclass(L), obs = obs[1:cJ,,drop = FALSE], lwr = lwr[-(1:cJ),,drop = FALSE], 
  +                  upr = upr[-(1:cJ),,drop = FALSE]), ncol = N)
  +   chk(S, Lower_tri(s$Lambda))
  + 
  +   S <- matrix(grad(f, obs[1:cJ,,drop = FALSE], lwr = lwr[-(1:cJ),,drop = FALSE], 
  +                     upr = upr[-(1:cJ),,drop = FALSE], L = L), ncol = N)
  +   chk(S, s$obs)
  +   S <- matrix(grad(f, lwr[-(1:cJ),,drop = FALSE], obs = obs[1:cJ,,drop = FALSE], 
  +                     upr = upr[-(1:cJ),,drop = FALSE], L = L), ncol = N)
  +   chk(S, s$lower)
  +   S <- matrix(grad(f, upr[-(1:cJ),,drop = FALSE], obs = obs[1:cJ,,drop = FALSE], 
  +                     lwr = lwr[-(1:cJ),,drop = FALSE], L = L), ncol = N)
  +   chk(S, s$upper)
  + 
  + })
  > 
  > 
  > if (require("numDeriv", quietly = TRUE)) {
  + 
  + J <- (cJ <- 5) + (dJ <- 6)
  + N <- 3
  + M <- 10
  + 
  + eval(thischeck)
  + 
  + J <- (cJ <- 1) + (dJ <- 1)
  + 
  + eval(thischeck)
  + 
  + J <- (cJ <- 1) + (dJ <- 4)
  + 
  + eval(thischeck)
  + 
  + J <- (cJ <- 4) + (dJ <- 1)
  + 
  + eval(thischeck)
  + 
  + }
  > 
  > proc.time()
     user  system elapsed 
    27.28    1.20   28.53