Running 'brfit.R' [4s]
  Running 'lmekin1.R' [3s]
  Comparing 'lmekin1.Rout' to 'lmekin1.Rout.save' ... OK
  Running 'lmekin2.R' [166s]
  Running 'nest1.R' [2s]
  Comparing 'nest1.Rout' to 'nest1.Rout.save' ... OK
  Running 'nest2.R' [2s]
  Comparing 'nest2.Rout' to 'nest2.Rout.save' ... OK
  Running 'refine.R' [2s]
  Comparing 'refine.Rout' to 'refine.Rout.save' ... OK
  Running 'refine1.R' [145s]
  Running 'refine2.R' [155s]
  Running 'refine3.R' [130s]
  Running 'shrink.R' [1s]
  Comparing 'shrink.Rout' to 'shrink.Rout.save' ... OK
  Running 'slope1.R' [116s]
  Running 'slope2.R' [100s]
  Running 'slope3.R' [2s]
  Comparing 'slope3.Rout' to 'slope3.Rout.save' ... OK
  Running 'slope4.R' [5s]
  Comparing 'slope4.Rout' to 'slope4.Rout.save' ... OK
  Running 'strata1.R' [83s]
  Running 'test0.R' [1s]
  Comparing 'test0.Rout' to 'test0.Rout.save' ... OK
  Running 'test1.R' [1s]
  Comparing 'test1.Rout' to 'test1.Rout.save' ... OK
  Running 'test2.R' [1s]
  Comparing 'test2.Rout' to 'test2.Rout.save' ... OK
  Running 'test3.R' [92s]
  Running 'test4.R' [91s]
  Running 'twoterm.R' [2s]
  Comparing 'twoterm.Rout' to 'twoterm.Rout.save' ... OK
Running the tests in 'tests/lmekin2.R' failed.
Complete output:
  > library(coxme)
  Loading required package: survival
  Loading required package: bdsmatrix
  
  Attaching package: 'bdsmatrix'
  
  The following object is masked from 'package:base':
  
      backsolve
  
  > aeq <- function(x, y, ...) all.equal(as.vector(x), as.vector(y), ...)
  > #
  > # A set of tests using user-defined matrices
  > #  The lung data is handy, though models of time that ignore status aren't
  > #  very sensible.
  > #
  > fit1 <- lmekin(time ~ age + (1|ph.ecog), lung, vfixed=.5)
  > 
  > # Solve using the extended normal equations
  > #  The loglik is sum(y - hat y)^2 + b'P b/2 where the penalty matrix
  > #  is all zeros for the fixed effects, and A-inverse for the random
  > #  where b is the coefficient vector and A is the variance of the random effects
  > #  In our case A-inverse/2 = identity matrix.  
  > # If the variance is known this equation is easy to solve.
  > efit <- function(p, pmat) {
  +     tlung <- na.omit(lung[,c("time", "age", "ph.ecog")])
  +     cmat <- cbind(0,0, rbind(0,0, pmat))
  +     x <- with(tlung, cbind(1, age, ph.ecog==0, ph.ecog==1,
  +                           ph.ecog==2, ph.ecog==3))
  +     as.vector(solve(t(x) %*% x + p*cmat, t(x) %*% tlung$time))
  + }
  > efit1 <- efit(2, diag(4))
  > 
  > aeq(fixef(fit1), efit1[1:2])
  [1] TRUE
  > aeq(unlist(ranef(fit1)), efit1[3:6])
  [1] TRUE
  > 
  > # Now put in a more complex pmat (linear trend constraint)
  > # Since coxme is using the inverse variance as the penalty, we
  > #  need to give the inverse penalty as the "variance".
  > # (This is an lmekin test with non-zero off diagonal elements 
  > #  in the variance for the random effect, which is key to working
  > #  with kinship matrices.)
  > pmat <- matrix(c(2,-1,0,0, -1,2,-1,0, 0, -1, 2, -1, 0,0,-1,2),4)
  > dimnames(pmat) <- list(0:3, 0:3)
  > 
  > fit2 <- lmekin(time ~ age + (1|ph.ecog), lung, vfixed=.5,
  +                varlist= coxmeMlist(solve(pmat), rescale=F))
  > efit2 <- efit(2, pmat)
  > aeq(fixef(fit2), efit2[1:2])
  [1] TRUE
  > aeq(unlist(ranef(fit2)), efit2[3:6])
  [1] TRUE
  > 
  > 
  > proc.time()
     user  system elapsed 
     1.87    0.25    2.12 
Running the tests in 'tests/refine1.R' failed.
Complete output:
  > library(coxme)
  Loading required package: survival
  Loading required package: bdsmatrix
  
  Attaching package: 'bdsmatrix'
  
  The following object is masked from 'package:base':
  
      backsolve
  
  > # Test of refine.n option
  > #  Matches the first example in the laplace vignette
  > #
  > aeq <- function(x,y, ...) all.equal(as.vector(x), as.vector(y), ...)
  > nsim <- 100
  > set.seed(1953)  #an auspicious birth year :-)
  > mkdata <- function(n, beta=c(.4, .1), sitehaz=c(.5,1.5, 2,1)) {
  +     nsite <- length(sitehaz)
  +     site <- rep(1:nsite, each=n)
  +     trt1 <- rep(0:1, length=n*nsite)
  +     hazard <- sitehaz[site] + beta[1]*trt1 + beta[2]*trt1 * (site-mean(site))
  +     stime <- rexp(n*nsite, exp(hazard))
  +     q80 <- quantile(stime, .8)
  +     data.frame(site=site,
  +                trt = trt1,
  +                trt2 = 1-trt1,
  +                futime= pmin(stime, q80),
  +                status= ifelse(stime>q80, 0, 1),
  +                hazard=hazard
  +                )
  + }
  > trdata <- mkdata(150)  #150 enrolled per site
  > 
  > set.seed(12345)
  > 
  > fit1 <- coxme(Surv(futime, status) ~  trt + (1| site/trt), trdata,
  +               refine.n=nsim, refine.detail=TRUE)
  > 
  > hmatb <- fit1$hmat[1:12, 1:12]
  > hmat.inv <- as.matrix(solve(hmatb))
  > 
  > chidf <- coxme.control()$refine.df
  > set.seed(12345)
  > bmat <- matrix(rnorm(12*nsim), ncol=nsim)
  > 
  > b2 <- backsolve(hmatb, bmat, upper=TRUE)
  > htran <- as(hmatb, "dtCMatrix")  #verify that backsolve works correctly
  > all.equal(as.matrix(htran %*% b2), bmat, check.attributes=FALSE)
  [1] TRUE
  > 
  > b2 <- scale(b2, center=F, scale=sqrt(rchisq(nsim, df=chidf)/chidf))
  > b3 <- b2 + unlist(ranef(fit1)) 
  > aeq(b3, fit1$refine.detail$bmat)
  [1] TRUE
  > 
  > # Check the fitted Cox models
  > logvec <- double(nsim)
  > indx <- -1 + trdata$trt + 2*trdata$site  #random treatment effect index
  > for (i in 1:nsim) {
  +     eta <- fixef(fit1) * trdata$trt + b3[indx, i] +
  +             b3[trdata$site+8,i] 
  +     tfit <- coxph(Surv(futime, status) ~ offset(eta), data= trdata)
  +     logvec[i] <- tfit$loglik[1] 
  + }
  > aeq(logvec, fit1$refine.detail$loglik)
  [1] TRUE
  > 
  > # The penalties
  > temp <- unlist(VarCorr(fit1))
  > Ainv <- diag(rep(1/temp, c(8,4)))
  > p1 <- diag(t(b3) %*% Ainv %*% b3)/2
  > aeq(p1, fit1$refine.detail$penalty1)
  [1] TRUE
  > 
  > p2 <- mahalanobis(t(b3), center=unlist(ranef(fit1)), cov=hmat.inv)
  > p3 <- mahalanobis(t(b2), center=rep(0,12), cov=hmat.inv)
  > aeq(p2, p3)
  [1] TRUE
  > aeq(p2/2, fit1$refine.detail$penalty2)
  [1] TRUE
  > 
  > # The log-density of the t
  > require(mvtnorm)
  Loading required package: mvtnorm
  > tdens <- dmvt(t(b3), delta=unlist(ranef(fit1)), sigma=hmat.inv, df=chidf)
  > aeq(tdens, fit1$refine.detail$tdens)
  [1] TRUE
  >  
  > # The normalization constant for the Gaussian penalty
  > #  Determinants are easy for a diagonal matrix
  > gnorm <- -(6*log(2*pi) - .5*sum(log(diag(Ainv))))
  > 
  > # My errhat vector
  > t1 <- logvec + gnorm - (tdens + p1 + fit1$loglik[2])
  > t2 <- fit1$loglik[3] + gnorm - (tdens + p2/2 + fit1$loglik[2])
  > errhat <- exp(t1) - exp(t2)
  > aeq(errhat, fit1$refine.detail$errhat)
  [1] TRUE
  > 
  > proc.time()
     user  system elapsed 
     2.01    0.31    2.31 
Running the tests in 'tests/refine2.R' failed.
Complete output:
  > #
  > # Check out the code on a simulated data set.  On first cut it appeared to
  > #  do very badly, this turned out to be a bug in coxfit6d.c, when there
  > #  were random slopes it miscalculated the linear predictor.
  > #
  > library(coxme)
  Loading required package: survival
  Loading required package: bdsmatrix
  
  Attaching package: 'bdsmatrix'
  
  The following object is masked from 'package:base':
  
      backsolve
  
  > aeq <- function(x,y) all.equal(as.vector(x), as.vector(y))
  > 
  > set.seed(1979)
  > mkdata <- function(n, beta=c(.4, .1), sitehaz=c(.5,1.5, 2,1)) {
  +     nsite <- length(sitehaz)
  +     site <- rep(1:nsite, each=n)
  +     trt1 <- rep(0:1, length=n*nsite)
  +     hazard <- sitehaz[site] + beta[1]*trt1 + beta[2]*trt1 * (site-mean(site))
  +     stime <- rexp(n*nsite, exp(hazard))
  +     q80 <- quantile(stime, .8)
  +     data.frame(site=site,
  +                trt1 = trt1,
  +                trt2 = 1-trt1,
  +                futime= pmin(stime, q80),
  +                status= ifelse(stime>q80, 0, 1),
  +                hazard=hazard
  +                )
  +     }
  > 
  > trdata <- mkdata(150)  #150 enrolled per site
  > #fixf <- coxph(Surv(futime, status) ~ factor(site)*trt1, trdata)
  > #pfit <- pyears(Surv(futime, status) ~ site + trt1, trdata)
  > #(pfit$event/sum(pfit$event))/ (pfit$pyears/sum(pfit$pyears))
  > 
  > set.seed(50)
  > #nsim <- 500
  > nsim <- 20  # speedup for CRAN
  > 
  > fit <- coxme(Surv(futime, status) ~ trt2 + (1 + trt2 | site), trdata,
  +              refine.n=nsim, refine.detail=TRUE)
  > debug <- fit$refine.detail
  > 
  > # Recreate the variance-covariance and sim data that was used
  > set.seed(50)
  > bmat <- matrix(rnorm(8*nsim), nrow=8)
  > hmatb <- fit$hmat[1:8, 1:8]
  > hmat.inv <- as.matrix(solve(hmatb))
  > 
  > chidf <- coxme.control()$refine.df
  > 
  > b2 <- backsolve(hmatb, bmat, upper=TRUE)
  > htran <- as(hmatb, "dtCMatrix")  #verify that backsolve works correctly
  > all.equal(as.matrix(htran %*% b2), bmat, check.attributes=FALSE)
  [1] TRUE
  > 
  > b2 <- scale(b2, center=F, scale=sqrt(rchisq(nsim, df=chidf)/chidf))
  > b3 <- b2 + unlist(ranef(fit)) 
  > aeq(b3, debug$bmat)
  [1] TRUE
  > 
  > 
  > temp <- VarCorr(fit)[[1]]
  > sigma <- diag(c(rep(temp[1],4), rep(temp[4],4)))
  > sigma[cbind(1:4,5:8)] <- temp[3]* sqrt(temp[1] * temp[4])
  > sigma[cbind(5:8,1:4)] <- temp[3]* sqrt(temp[1] * temp[4])
  > 
  > if (!is.null(debug))
  +     all.equal(as.matrix(gchol(sigma), ones=F), as.matrix(debug$gkmat, ones=F))
  [1] TRUE
  > 
  > coxll <- double(nsim)
  > for (i in 1:nsim) {
  +     lp <- trdata$trt2 * fixef(fit) + b3[trdata$site,i] +
  +            b3[trdata$site +4,i] * trdata$trt2
  +     tfit <- coxph(Surv(futime, status) ~ offset(lp), trdata)
  +     coxll[i] <- tfit$loglik
  +     }
  > if (!is.null(debug)) all.equal(coxll, debug$log)
  [1] TRUE
  > 
  > # How does the direct average compare to the IPL?
  > # (This is not guarranteed to be accurate, just curious)
  > constant <- .5*(log(2*pi) + sum(log(diag(gchol(sigma)))))
  > fit$log[2] + c(IPL=0, sim=log(mean(exp(coxll-fit$log[2]))) - constant)
        IPL       sim 
  -2715.550 -2689.821 
  > 
  > # Compute the Taylor series for the IPL
  > bhat <- unlist(ranef(fit))
  > b.sig <- t(b2) %*% fit$hmat[1:8, 1:8]  # b times sqrt(H)
  > taylor <- rowSums(b.sig^2)/2   # vector of b'H b/2
  > aeq(taylor, debug$penalty2)
  [1] TRUE
  > 
  > # Look at how well the Taylor series does
  > pen <- solve(sigma)
  > simpen <- colSums(b3 *(pen%*%b3))/2 #penalty for each simulated b
  > aeq(simpen, debug$penalty1)
  [1] TRUE
  > 
  > #plot(coxll- simpen, fit$log[3] - taylor, 
  > #     xlab="Actual PPL for simulated points", 
  > #     ylab="Taylor series approximation")
  > #abline(0,1)
  > #
  > 
  > #And now I need the Gaussian integration contant and the t-density
  > require(mvtnorm)
  Loading required package: mvtnorm
  > tdens <- dmvt(t(b3), delta=unlist(ranef(fit)), sigma=hmat.inv, df=chidf)
  > aeq(tdens, debug$tdens)
  [1] TRUE
  >  
  > # The normalization constant for the Gaussian penalty
  > #  
  > temp <- determinant(sigma, log=TRUE)
  > gnorm <- -(4*log(2*pi) + .5*temp$sign *temp$modulus)
  > aeq(gnorm, debug$gdens)
  [1] TRUE
  > 
  > m2 <- fit$loglik[2]
  > errhat <- exp(coxll+gnorm -(simpen + tdens +m2)) - 
  +          exp(fit$log[3]+ gnorm -(taylor + tdens + m2))
  > if (!is.null(debug)) {
  +     aeq(errhat, debug$errhat)
  +     }
  [1] TRUE
  > 
  > proc.time()
     user  system elapsed 
     1.56    0.18    1.71 
Running the tests in 'tests/refine3.R' failed.
Complete output:
  > library(coxme)
  Loading required package: survival
  Loading required package: bdsmatrix
  
  Attaching package: 'bdsmatrix'
  
  The following object is masked from 'package:base':
  
      backsolve
  
  > aeq <- function(x,y) all.equal(as.vector(x), as.vector(y))
  > 
  > #
  > # Check out the code on a simulated data set, 
  > #  arose out a concern that the first derivatives were wrong
  > #  when the model had random slopes and intercepts
  > # (The code was ok, but another test case is always good.)
  > # The heart of the test is to re-compute  quantities using
  > #   a combination of coxph or matrix algebra, and check these with 
  > #   the coxme answers.
  > #
  > approx <- "efron"
  > set.seed(9291978)
  > mkdata <- function(n, beta=c(.4, .1), sitehaz=c(.5,1.5, 2,1)) {
  +     nsite <- length(sitehaz)
  +     site <- rep(1:nsite, each=n)
  +     trt1 <- rep(0:1, length=n*nsite)
  +     hazard <- sitehaz[site] + beta[1]*trt1 + beta[2]*trt1 * (site-mean(site))
  +     stime <- rexp(n*nsite, exp(hazard))
  +     q80 <- quantile(stime, .8)
  +     data.frame(site=site,
  +                trt1 = trt1,
  +                trt2 = 1-trt1,
  +                futime= pmin(stime, q80),
  +                status= ifelse(stime>q80, 0, 1),
  +                hazard=hazard
  +                )
  +     }
  > 
  > trdata <- mkdata(150)  #150 enrolled per site
  > fixf <- coxph(Surv(futime, status) ~ factor(site)*trt1, trdata,
  +               method=approx)
  > 
  > fit1 <- coxme(Surv(futime, status) ~ trt1 + (1 + trt1 | site), trdata,
  +               ties=approx)
  > fit2 <- coxme(Surv(futime, status) ~ trt2 + (1 + trt2 | site), trdata,
  +               ties=approx)
  > # Fit1 and fit2 are actually different models, with different 
  > #  iteration paths and maxima.
  > 
  > # Test out loglik
  > eta <- fit1$linear
  > fit1b <- coxph(Surv(futime, status) ~ eta, trdata, init=1, iter=0,
  +                method=approx)
  > aeq(fit1b$loglik[2], fit1$loglik[3]+fit1$penal)
  [1] TRUE
  > 
  > # First derivatives
  > temp <- VarCorr(fit1)[[1]]
  > sigma <- diag(c(rep(temp[1],4), rep(temp[4],4)))
  > sigma[cbind(1:4,5:8)] <- temp[3]* sqrt(temp[1] * temp[4])
  > sigma[cbind(5:8,1:4)] <- temp[3]* sqrt(temp[1] * temp[4])
  > pen <- matrix(0., 9,9)
  > pen[1:8, 1:8] <- solve(sigma)
  > pcoef <- c(unlist(ranef(fit1)), fixef(fit1))
  > aeq(pcoef %*% pen %*% pcoef/2, fit1$penal)
  [1] TRUE
  > 
  > xx <- 1* with(trdata, cbind(site==1, site==2, site==3, site==4))
  > xx <- cbind(xx, xx*trdata$trt1)
  > 
  > cfit1 <- coxph(Surv(futime, status) ~ xx + trt1, iter=0,
  +                method=approx, init=pcoef, trdata)
  > aeq(cfit1$log, fit1b$log) #check my typing
  [1] TRUE
  > 
  > # Using coxph.detail, I find the first and second derivatives
  > #  via a completely different code path than coxme
  > dt1 <- coxph.detail(cfit1)
  > first <- colSums(dt1$score) - pcoef%*%pen
  > aeq(first, fit1$u)  #do they agree?
  [1] TRUE
  > 
  > i2 <- apply(dt1$imat, 1:2, sum) + pen  #information is -1*second deriv
  > aeq(as.matrix(fit1$hmat, one=F), as.matrix(gchol(i2), ones=F))
  [1] TRUE
  > 
  > proc.time()
     user  system elapsed 
     1.96    0.10    2.06 
Running the tests in 'tests/slope1.R' failed.
Complete output:
  > library(coxme)
  Loading required package: survival
  Loading required package: bdsmatrix
  
  Attaching package: 'bdsmatrix'
  
  The following object is masked from 'package:base':
  
      backsolve
  
  > options(na.action='na.exclude', contrasts=c('contr.treatment', 'contr.poly'))
  > aeq <- function(x,y) all.equal(as.vector(x), as.vector(y))
  > 
  > #
  > # Test of fitting random slopes
  > #
  > # Simulation data with 9 institutions, strong age effects
  > #  and a random treatment effect
  > #      
  > set.seed(56)
  > n.subject <- seq(180, by=21, length=9) # number of subjects
  > slope <- sort(-.5 + rnorm(9, sd=.5))         # true treament effects
  > 
  > inst <- rep(1:9, n.subject)
  > n <- length(inst)
  > simdata <- data.frame(id=1:n, inst=inst,
  +                       trt= rep(0:1, length=n),
  +                       age= runif(n, 40, 70))
  > #risk goes up 30%/decade of age
  > simdata$hazard <- .8* exp(simdata$trt * rep(slope, n.subject) +
  +                           (simdata$age-55) * .03)
  > 
  > rtime <- function(hazard, censor=c(1,2)) {
  +     stime <- rexp(length(hazard), rate=hazard)
  +     ctime <- runif(length(hazard), censor[1], censor[2])
  +     list(time= pmin(stime, ctime), status=1*(stime <=ctime))
  +     }
  > temp <- rtime(simdata$hazard)
  > simdata$time <- temp$time
  > simdata$status <- temp$status
  > 
  > coef0 <- matrix(0., 2,9)
  > for (i in 1:9) {
  +     fit0 <- coxph(Surv(time, status) ~ age + trt, simdata,
  +                   subset=(inst==i))
  +     coef0[,i] <- fit0$coef
  +     }
  > 
  > # Several of these fits will differ in the last few digits on a 
  > #   64 bit vs 32 bit Intel processer.  The loglike is very
  > #   flat on top so tiny changes in the compute path lead to a small
  > #   change in the final solution.  Hence the "digits" argument.
  > fit0 <- coxph(Surv(time, status) ~ age + trt, simdata)
  > fit1 <- coxme(Surv(time, status) ~ age + trt + (1|inst), simdata)
  > print(fit1, rcoef=TRUE, digits=4)
  Cox mixed-effects model fit by maximum likelihood
    Data: simdata
    events, n = 1504, 2376
    Iterations= 8 35 
                      NULL Integrated    Fitted
  Log-likelihood -10892.32  -10826.57 -10819.37
  
                    Chisq   df p   AIC    BIC
  Integrated loglik 131.5 3.00 0 125.5 109.55
   Penalized loglik 145.9 7.38 0 131.1  91.88
  
  Model:  Surv(time, status) ~ age + trt + (1 | inst) 
  Fixed and penalized coefficients
              coef exp(coef) se(coef)     z       p
  age     0.029948    1.0304 0.003023  9.91 0.0e+00
  trt    -0.249769    0.7790 0.051899 -4.81 1.5e-06
  inst.1 -0.122696    0.8845 0.078439              
  inst.2 -0.007691    0.9923 0.075521              
  inst.3 -0.043795    0.9572 0.074255              
  inst.4 -0.120717    0.8863 0.073124              
  inst.5  0.032154    1.0327 0.070791              
  inst.6  0.091249    1.0955 0.069919              
  inst.7 -0.017063    0.9831 0.068940              
  inst.8  0.024342    1.0246 0.067913              
  inst.9  0.164217    1.1785 0.066641              
  
  Random effects
   Group Variable  Std Dev Variance
   inst  Intercept 0.11329 0.01283 
  > 
  > fit2 <- coxme(Surv(time, status) ~ age + trt + (1|inst/trt), simdata)
  > print(fit2, rcoef=TRUE, digits=3)
  Cox mixed-effects model fit by maximum likelihood
    Data: simdata
    events, n = 1504, 2376
    Iterations= 14 60 
                      NULL Integrated Fitted
  Log-likelihood -10892.32  -10814.83 -10796
  
                    Chisq   df p AIC   BIC
  Integrated loglik   155  4.0 0 147 125.7
   Penalized loglik   193 14.2 0 164  88.8
  
  Model:  Surv(time, status) ~ age + trt + (1 | inst/trt) 
  Fixed and penalized coefficients
                    coef exp(coef) se(coef)     z      p
  age           0.030902     1.031  0.00304 10.18 0.0000
  trt          -0.303468     0.738  0.10867 -2.79 0.0052
  inst/trt.1/0 -0.037941     0.963  0.12209             
  inst/trt.1/1 -0.262439     0.769  0.13045             
  inst/trt.2/0  0.093971     1.099  0.11654             
  inst/trt.2/1 -0.124835     0.883  0.12498             
  inst/trt.3/0  0.111253     1.118  0.11498             
  inst/trt.3/1 -0.220075     0.802  0.12246             
  inst/trt.4/0 -0.123843     0.884  0.11411             
  inst/trt.4/1 -0.155038     0.856  0.12008             
  inst/trt.5/0  0.043624     1.045  0.11024             
  inst/trt.5/1  0.027408     1.028  0.11527             
  inst/trt.6/0  0.154576     1.167  0.10811             
  inst/trt.6/1  0.042069     1.043  0.11446             
  inst/trt.7/0 -0.035562     0.965  0.10859             
  inst/trt.7/1 -0.000864     0.999  0.11119             
  inst/trt.8/0 -0.149839     0.861  0.10780             
  inst/trt.8/1  0.227553     1.256  0.10864             
  inst/trt.9/0 -0.056239     0.945  0.10533             
  inst/trt.9/1  0.466220     1.594  0.10656             
  inst.1       -0.001390     0.999  0.01364             
  inst.2       -0.000143     1.000  0.01364             
  inst.3       -0.000504     0.999  0.01364             
  inst.4       -0.001290     0.999  0.01364             
  inst.5        0.000329     1.000  0.01364             
  inst.6        0.000910     1.001  0.01364             
  inst.7       -0.000169     1.000  0.01364             
  inst.8        0.000360     1.000  0.01364             
  inst.9        0.001897     1.002  0.01364             
  
  Random effects
   Group    Variable    Std Dev  Variance
   inst/trt (Intercept) 0.201142 0.040458
   inst     (Intercept) 0.013682 0.000187
  > 
  > # And so will this one
  > fit3 <- coxme(Surv(time, status) ~ age + trt + (1|inst) + (trt|inst),simdata)
  > print(fit3, rcoef=TRUE, digits=3)
  Cox mixed-effects model fit by maximum likelihood
    Data: simdata
    events, n = 1504, 2376
    Iterations= 25 104 
                      NULL Integrated    Fitted
  Log-likelihood -10892.32  -10812.96 -10798.58
  
                    Chisq   df p AIC BIC
  Integrated loglik   159  4.0 0 151 129
   Penalized loglik   188 10.8 0 166 109
  
  Model:  Surv(time, status) ~ age + trt + (1 | inst) + (trt | inst) 
  Fixed and penalized coefficients
                 coef exp(coef) se(coef)     z      p
  age         0.03073     1.031  0.00303 10.14 0.0000
  trt        -0.30570     0.737  0.10695 -2.86 0.0043
  inst.1     -0.02143     0.979  0.05618             
  inst.2      0.02234     1.023  0.05563             
  inst.3      0.02398     1.024  0.05541             
  inst.4     -0.04352     0.957  0.05443             
  inst.5      0.01729     1.017  0.05414             
  inst.6      0.05591     1.058  0.05403             
  inst.7     -0.01037     0.990  0.05348             
  inst.8     -0.04290     0.958  0.05284             
  inst.9     -0.00131     0.999  0.05249             
  inst.1:trt -0.29599     0.744  0.16082             
  inst.2:trt -0.15830     0.854  0.15333             
  inst.3:trt -0.26950     0.764  0.15102             
  inst.4:trt -0.13695     0.872  0.14825             
  inst.5:trt  0.02349     1.024  0.14239             
  inst.6:trt  0.00704     1.007  0.14141             
  inst.7:trt  0.01545     1.016  0.13845             
  inst.8:trt  0.29617     1.345  0.13549             
  inst.9:trt  0.51860     1.680  0.13306             
  
  Random effects
   Group Variable  Std Dev Variance
   inst  Intercept 0.06266 0.00393 
   inst  trt       0.27847 0.07754 
  > 
  > fit4 <- coxme(Surv(time, status) ~ age + trt + (1 +trt |inst), simdata)
  > 
  > sfit0 <- coxph(Surv(time, status) ~ age + trt + strata(inst), simdata)
  > sfit1 <- coxme(Surv(time, status) ~ age + trt + (trt|inst) + strata(inst),
  +                simdata)
  > print(sfit1, rcoef=TRUE, digits=4)
  Cox mixed-effects model fit by maximum likelihood
    Data: simdata
    events, n = 1504, 2376
    Iterations= 7 31 
                      NULL Integrated    Fitted
  Log-likelihood -7636.963  -7564.201 -7553.573
  
                    Chisq  df p   AIC   BIC
  Integrated loglik 145.5 3.0 0 139.5 123.6
   Penalized loglik 166.8 8.3 0 150.2 106.1
  
  Model:  Surv(time, status) ~ age + trt + (trt | inst) + strata(inst) 
  Fixed and penalized coefficients
                 coef exp(coef) se(coef)     z      p
  age         0.03097    1.0315 0.003047 10.16 0.0000
  trt        -0.31038    0.7332 0.116564 -2.66 0.0077
  inst.1:trt -0.23574    0.7900 0.190686             
  inst.2:trt -0.23583    0.7899 0.179152             
  inst.3:trt -0.34250    0.7100 0.176387             
  inst.4:trt -0.03709    0.9636 0.174728             
  inst.5:trt -0.01858    0.9816 0.165983             
  inst.6:trt -0.09182    0.9123 0.163471             
  inst.7:trt  0.03376    1.0343 0.162333             
  inst.8:trt  0.39190    1.4798 0.159318             
  inst.9:trt  0.53589    1.7090 0.155067             
  
  Random effects
   Group Variable Std Dev Variance
   inst  trt      0.31042 0.09636 
  > 
  > # Check that the start,stop code does the same
  > dummy <- runif(nrow(simdata), -4, -1)  #all start times before first event
  > fit4b <- coxme(Surv(dummy, time, status) ~ age + trt + (1 +trt |inst), simdata)
  > all.equal(fit4b$loglik, fit4$loglik)
  [1] TRUE
  > all.equal(fit4b$coef, fit4$coef, tolerance=1e-7) # different order of internal
  [1] TRUE
  >                                                # sums => tiny difference
  > 
  > #Comparison plot
  > y <- cbind(slope, NA, coef0[2,], fixef(sfit1)[2] + unlist(ranef(sfit1)),
  +            fixef(fit3)[2] + ranef(fit3)[[2]],
  +            fixef(fit4)[2] + ranef(fit4)[[1]][,2])
  > matplot(c(1, 1.5, 2:5), t(y), type='b', xaxt='n', xlab="Simulation", 
  +         ylab="Treatment coefficient", lty=1)
  > axis(1, 1:5, c("Sim", "Separate", "Strata", "Uncor", "Corr"))
  > 
  > 
  > #
  > # Now compute some things exactly
  > #
  > contr.none <- function(n,contrasts=T) {
  +         if(is.numeric(n) && length(n) == 1.)
  +                 levs <- 1.:n
  +         else {
  +                 levs <- n
  +                 n <- length(n)
  +         }
  +         contr <- array(0., c(n, n), list(levs, levs))
  +         contr[seq(1., n^2., n + 1.)] <- 1.
  +         contr
  +         }
  > options(contrasts=c('contr.none', 'contr.poly'))
  > igchol <- function(x) {
  +     dd <- diag(x)
  +     ll <- as.matrix(x)
  +     ll %*% diag(dd) %*% t(ll)
  +     }
  > 
  > # For fit2
  > vtemp <- unlist(VarCorr(fit2))
  > names(vtemp) <- names(VarCorr(fit2))
  > fit2a <- coxme(Surv(time, status) ~ age + trt + (1|inst/trt), simdata,
  +                iter=0, vfixed=vtemp)
  > temp <- strata(simdata$inst, simdata$trt, sep='/', shortlabel=TRUE)
  > cfit <- coxph(Surv(time, status) ~ factor(temp) +factor(inst) +age + trt,
  +               simdata, iter=0, x=T)
  > dt2 <- coxph.detail(cfit)
  > u2 <- apply(dt2$score, 2, sum)
  > aeq(u2, fit2a$u)
  [1] TRUE
  > imat2 <- apply(dt2$imat, 1:2, sum) + diag(c(rep(1/vtemp, c(18,9)),0,0))
  > aeq(imat2, as.matrix(igchol(fit2a$hmat)))
  [1] TRUE
  > 
  > # For fit3
  > vtemp <- as.vector(unlist(VarCorr(fit3)))  #name not needed
  > fit3a <- coxme(Surv(time, status) ~ age + trt + (1|inst) + (trt|inst),
  +                simdata, iter=0, vfixed=as.list(vtemp))
  > cfit <- coxph(Surv(time, status) ~ factor(inst) * trt + age, simdata,
  +               iter=0, x=T)
  > dt3 <- coxph.detail(cfit)
  > u3 <- apply(dt3$score, 2, sum)
  > indx <- c(1:9, 12:20, 11, 10)
  > aeq(u3[indx], fit3a$u)
  [1] TRUE
  > imat2 <- apply(dt3$imat, 1:2, sum)[indx,indx] + 
  +     diag(c(rep(1/vtemp, c(9,9)),0,0))
  > aeq(imat2, as.matrix(igchol(fit3a$hmat)))
  [1] TRUE
  > 
  > fit3b <- coxme(Surv(time, status) ~ age + trt + (trt|inst) +(1|inst),
  +                simdata, iter=0, vfixed=as.list(rev(vtemp)))
  > aeq(fit3a$u, fit3b$u)
  [1] TRUE
  > aeq(fit3b$imat, fit3b$imat)
  [1] TRUE
  > 
  > #For sfit1
  > vtemp <- .0966
  > fit <- coxme(Surv(time, status) ~ age + trt + strata(inst) + (trt|inst),
  +                simdata, iter=0, vfixed=vtemp)
  > cfit <- coxph(Surv(time, status) ~ factor(inst):trt + trt+ age+ strata(inst),
  +               simdata, iter=0, x=T)
  > dt3 <- coxph.detail(cfit)
  > u3 <- apply(dt3$score, 2, sum)
  > indx <- c(3:11,2,1)
  > aeq(u3[indx], fit$u)
  [1] TRUE
  > 
  > imat3 <- apply(dt3$imat,1:2, sum)[indx,indx] + diag(c(rep(1/vtemp,9),0,0))
  > aeq(imat3, as.matrix(igchol(fit$hmat)))
  [1] TRUE
  > 
  > 
  > proc.time()
     user  system elapsed 
     7.25    0.12    7.35 
Running the tests in 'tests/slope2.R' failed.
Complete output:
  > library(coxme)
  Loading required package: survival
  Loading required package: bdsmatrix
  
  Attaching package: 'bdsmatrix'
  
  The following object is masked from 'package:base':
  
      backsolve
  
  > options(na.action='na.exclude', contrasts=c('contr.treatment', 'contr.poly'))
  > aeq <- function(x,y, ...) all.equal(as.vector(x), as.vector(y), ...)
  > 
  > #
  > # Test of fitting random slopes
  > #
  > # Same data set as slope1, refit using variance matrices
  > #
  > set.seed(56)
  > n.subject <- seq(180, by=21, length=9) # number of subjects
  > slope <- sort(-.5 + rnorm(9, sd=.5))         # true treament effects
  > 
  > inst <- rep(1:9, n.subject)
  > n <- length(inst)
  > simdata <- data.frame(id=1:n, inst=inst,
  +                       trt= rep(0:1, length=n),
  +                       age= runif(n, 40, 70))
  > #risk goes up 30%/decade of age
  > simdata$hazard <- .8* exp(simdata$trt * rep(slope, n.subject) +
  +                           (simdata$age-55) * .03)
  > 
  > rtime <- function(hazard, censor=c(1,2)) {
  +     stime <- rexp(length(hazard), rate=hazard)
  +     ctime <- runif(length(hazard), censor[1], censor[2])
  +     list(time= pmin(stime, ctime), status=1*(stime <=ctime))
  +     }
  > temp <- rtime(simdata$hazard)
  > simdata$time <- temp$time
  > simdata$status <- temp$status
  > 
  > #Check fit1
  > fit1 <- coxme(Surv(time, status) ~ age + trt + (1|inst), simdata)
  > fit1b <- coxme(Surv(time, status) ~ age + trt + (1|inst), simdata,
  +                varlist=diag(9))
  > aeq(fit1$log, fit1b$log)
  [1] TRUE
  > aeq(as.matrix(fit1$var), as.matrix(fit1b$var))
  [1] TRUE
  > aeq(fixef(fit1), fixef(fit1b))
  [1] TRUE
  > 
  > # Check fit2
  > #  To get the same iteration path you need to force the same
  > #  starting estimates
  > idlist <- sort(outer(1:9, 0:1, paste, sep='/'))
  > fit2 <- coxme(Surv(time, status) ~ age + trt + (1|inst/trt), simdata,
  +                varlist=coxmeFull(collapse=TRUE), vinit=c(.1, .1))
  > 
  > mat1 <- matrix(diag(18), 18, dimnames=list(idlist, idlist))
  > mat2 <-  bdsBlock(idlist, rep(1:9, each=2))
  > fit2b <- coxme(Surv(time, status) ~ age + trt + (1|inst/trt), simdata,
  +                varlist=list(mat1, mat2), vinit=c(.1, .1))
  > aeq(fit2$log, fit2b$log)
  [1] TRUE
  > aeq(as.matrix(fit2$var), as.matrix(fit2b$var), tol=1e-7)
  [1] TRUE
  > aeq(fixef(fit2), fixef(fit2b))
  [1] TRUE
  > 
  > # Check fit3
  > # This takes two different paths to the solution, due to use of a
  > #  slope coefficient rather than nested.  So results are a bit
  > fit3 <- coxme(Surv(time, status) ~ age + trt + (1|inst) + (trt|inst),simdata,
  +               vinit=list(.1, .2))
  > mat3 <- diag(rep(0:1, 9))
  > dimnames(mat3) <- list(idlist, idlist)
  > fit3b <-  coxme(Surv(time, status) ~ age + trt + (1|inst/trt), simdata,
  +                varlist=coxmeMlist(list(mat2, mat3), rescale=F, pdcheck=F),
  +                vinit=c(.1, .2))
  > 
  > aeq(fit3$log, fit3b$log)
  [1] TRUE
  > aeq(fixef(fit3), fixef(fit3b))
  [1] TRUE
  > all.equal(unlist(VarCorr(fit3)), unlist(VarCorr(fit3b)), tolerance=1e-7, 
  +           check.attributes=FALSE)
  [1] TRUE
  > # This function should map between coefficients
  > map <- matrix(0, 20,20)
  > for (i in 1:9) {
  +     map[i, 2*i -1] <- 1
  +     map[i+9, 2*i -1] <- -1
  +     map[i+9, 2*i] <- 1
  +     }
  > map[19,19] <- 1
  > map[20,20] <- 1
  > 
  > # Some of the random effects are very close to zero
  > aeq(unlist(ranef(fit3)), c(map[1:18,1:18] %*% unlist(ranef(fit3b))),
  +     tol=1e-7)
  [1] TRUE
  > aeq(as.matrix(fit3$var), map %*% as.matrix(fit3b$var) %*% t(map),
  +     tol=1e-7)
  [1] TRUE
  > 
  > proc.time()
     user  system elapsed 
     2.73    0.09    2.81 
Running the tests in 'tests/strata1.R' failed.
Complete output:
  > library(coxme)
  Loading required package: survival
  Loading required package: bdsmatrix
  
  Attaching package: 'bdsmatrix'
  
  The following object is masked from 'package:base':
  
      backsolve
  
  > aeq <- function(x, y) all.equal(as.vector(x), as.vector(y))
  > options(na.action=na.exclude)
  > 
  > #
  > # trivial test of strata
  > #
  > tdata0 <- data.frame(time  =c(5,4,1,1,2,2,2,2,3), 
  +                      status=c(0,1,1,0,1,1,1,0,0),
  +                      x1    =c(0,1,2,0,1,1,0,1,0),
  +                      wt    =c(1,2,1,2,3,4,3,2,1),
  +                      x2    =c(1,3,5,2,3,6,4,3,1),
  +                      grp   =c(1,1,2,2,1,1,2,2,1))
  > 
  > rcnt <- c(3,2,1,1,1,2,2,1,3)  # rep count
  > tdata0b <- data.frame(time2  = c(3,4,5, 2,4, 1,1,2, 1,2, .5,2, 2, 1,2,3), 
  + 		      time1  = c(0,3,4, 0,2, 0,0,0, 0,1, 0,.5, 0, 0,1,2),
  + 		      status = c(0,0,0, 0,1, 1,0,1, 0,1, 0, 1, 0, 0,0,0),
  +                       x1     = rep(tdata0$x1, rcnt),
  +                       wt     = rep(tdata0$wt, rcnt),
  +                       x2     = rep(tdata0$x2, rcnt),
  +                       grp    = rep(tdata0$grp, rcnt))
  > 
  > # Undo a gchol: given gchol(x), returns x
  > igchol <- function(x) {
  +     dd <- diag(x)
  +     ll <- as.matrix(x)
  +     ll %*% diag(dd) %*% t(ll)
  +     }
  > 
  > tfit <- coxme(Surv(time, status) ~ x1 + x2 + (1|grp), data=tdata0,
  +               vfixed=.5, weight=wt, iter=0,
  +               ties='breslow', init=c(pi,0), sparse.calc=0)
  > 
  > sdata0 <- rbind(tdata0, tdata0)
  > sdata0$strat <- rep(1:2, each=nrow(tdata0))
  > 
  > sfit0 <- coxme(Surv(time,status) ~ x1 + x2 + strata(strat) + (1|grp), 
  +                data=sdata0, ties='breslow', iter=0, sparse.calc=0,
  + 	       vfixed=0.5, weight=wt, init=c(pi,0))
  > 
  > sfit1 <- coxme(Surv(time,status) ~ x1 + x2 + strata(strat) + (1|grp), 
  +                data=sdata0, ties='breslow', iter=0, sparse.calc=1,
  + 	       vfixed=0.5, weight=wt, init=c(pi,0))
  > 
  > aeq(sfit0$u, sfit1$u)
  [1] TRUE
  > all.equal(sfit0$hmat, sfit1$hmat)
  [1] TRUE
  > aeq(sfit0$u/2, tfit$u)
  [1] TRUE
  > tpen <- diag(c(2,2,0,0))
  > aeq(igchol(sfit0$hmat)-tpen, 2*(igchol(tfit$hmat)-tpen))
  [1] TRUE
  > 
  > #
  > # Repeat for start/stop data
  > #   The subset makes the test more rigorous (someone exits a risk
  > #   set without entering it immediatly again).
  > tfit <- coxme(Surv(time1, time2, status) ~ x1 + x2 + (1|grp), data=tdata0b,
  +               vfixed=0.6, weight=wt, iter=0,
  +               ties='efron', init=c(pi,0), subset=(-1))
  > 
  > sdata0b <- rbind(tdata0b[-1,], tdata0b[-1,])
  > sdata0b$strat <- rep(1:2, each=nrow(tdata0b)-1)
  > 
  > sfit0 <- coxme(Surv(time1,time2,status) ~ x1 + x2 + strata(strat) +(1|grp),
  + 	       data=sdata0b, iter=0, sparse.calc=0,
  + 	       vfixed=0.6, weight=wt, init=c(pi,0))
  > sfit1 <- coxme(Surv(time1,time2,status) ~ x1 + x2 + strata(strat) + (1|grp),
  + 	       data=sdata0b, iter=0, sparse.calc=1,
  + 	       vfixed=0.6, weight=wt, init=c(pi,0))
  > aeq(tfit$u, sfit0$u/2)
  [1] TRUE
  > aeq(sfit0$u, sfit1$u)
  [1] TRUE
  > tpen <- diag(c(1/.6, 1/.6, 0,0))
  > aeq(sfit0$imat, sfit1$imat)
  [1] TRUE
  > aeq(igchol(sfit0$hmat)-tpen, 2*(igchol(tfit$hmat)-tpen))
  [1] TRUE
  > 
  > proc.time()
     user  system elapsed 
     1.09    0.12    1.20 
Running the tests in 'tests/test3.R' failed.
Complete output:
  > library(coxme)
  Loading required package: survival
  Loading required package: bdsmatrix
  
  Attaching package: 'bdsmatrix'
  
  The following object is masked from 'package:base':
  
      backsolve
  
  > options(na.action='na.exclude', contrasts=c('contr.treatment', 'contr.poly'))
  > aeq <- function(x,y) all.equal(as.vector(x), as.vector(y))
  > 
  > # $Id: test3.s,v 1.4 2003/08/21 21:24:40 therneau Exp $
  > #
  > # Check out a mixed sparse/non-sparse cacluation
  > #
  > # The choice of "sparse" below forces group 1 to be non-sparse and
  > #   group 2-4 to be sparse.  It will also force the coefficients to be
  > #   in 2,3,4,1 order
  > # The test is much more complete if the grouping variable has 4 levels instead
  > #   of just 2.
  > #
  > tdata1 <- data.frame(time  =c(5,4,1,1,2,2,2,2,3, 1:8), 
  +                      status=c(0,1,1,0,1,1,1,0,0, rep(1:0,4)),
  +                      x1    =c(0,1,2,0,1,1,0,1,0, rep(4:1,2)),
  +                      wt    =c(1,2,1,2,3,4,3,2,1, rep(1:2,4)),
  +                      x2    =c(1,3,5,2,3,6,4,3,1, rep(1:2,4)),
  +                      grp   =c(1,1,2,2,1,1,2,2,1, rep(3,4), rep(4,4)))
  > 
  > theta=.77
  > zero <- rep(0, nrow(tdata1))
  > fit0 <- coxme(Surv(zero, time, status) ~ x1 + x2 + (1|grp), data=tdata1,
  +               vfixed=theta, weight=wt, iter=0,
  +               sparse=c(2, .25), ties='breslow')
  > 
  > tfit <- coxph(Surv(time, status) ~ I(grp==2) + I(grp==3) +I(grp==4) + 
  +               I(grp==1) + x1 + x2,
  +               data=tdata1, x=T, weight=wt, iter=0, method='breslow')
  > dt0 <- coxph.detail(tfit)
  > aeq(apply(dt0$score,2,sum), fit0$u)
  [1] TRUE
  > 
  > h0 <- apply(dt0$imat,1:2,sum) + diag(c(1/theta, 1/theta,1/theta, 1/theta,0,0))
  > h0[1,2:3] <- h0[2:3,1] <- 0  #sparse part
  > h0[2,3] <- h0[3,2] <- 0
  > 
  > aeq(as.matrix(gchol(h0)), as.matrix(fit0$hmat))
  [1] TRUE
  > aeq(diag(gchol(h0)), diag(fit0$hmat))
  [1] TRUE
  > hinv <- solve(h0)
  > hinv[1,2:3] <-hinv[2:3,1] <- 0
  > hinv[2,3] <- hinv[3,2] <- 0
  > aeq(hinv, as.matrix(fit0$var))  
  [1] TRUE
  > 
  > fit1 <- coxme(Surv(zero, time, status) ~ x1 + x2 + (1|grp), data=tdata1,
  +               vfixed=theta, weight=wt, iter=1,
  +               sparse=c(2, .25), ties='breslow')
  > 
  > temp <- solve(fit0$hmat, fit0$u)
  > temp[1:4] <- temp[1:4] - mean(temp[1:4])
  > aeq(temp, c(unlist(ranef(fit1)), fixef(fit1)))
  [1] TRUE
  > 
  > # Now test out Breslow/cox
  > fit0 <- coxme(Surv(time, status) ~ x1 + x2 + (1|grp), data=tdata1,
  +               vfixed=theta, weight=wt, iter=0,
  +               sparse=c(2, .25), ties='breslow')
  > aeq(apply(dt0$score,2,sum), fit0$u)
  [1] TRUE
  > aeq(as.matrix(gchol(h0)), as.matrix(fit0$hmat))
  [1] TRUE
  > aeq(diag(gchol(h0)), diag(fit0$hmat))
  [1] TRUE
  > aeq(hinv, as.matrix(fit0$var))  
  [1] TRUE
  > 
  > #Efron/Cox
  > fit0 <- coxme(Surv(time, status) ~ x1 + x2 + (1|grp), data=tdata1,
  +               vfixed=theta, weight=wt, iter=0,
  +               sparse=c(2, .25))
  > 
  > tfit <- coxph(Surv(time, status) ~ I(grp==2) + I(grp==3) +I(grp==4) + 
  +               I(grp==1) + x1 + x2,
  +               data=tdata1, x=T, weight=wt, iter=0)
  > dt0 <- coxph.detail(tfit)
  > aeq(apply(dt0$score,2,sum), fit0$u)
  [1] TRUE
  > 
  > h0 <- apply(dt0$imat,1:2,sum) + diag(c(1/theta, 1/theta,1/theta, 1/theta,0,0))
  > h0[1,2:3] <- h0[2:3,1] <- 0  #sparse part
  > h0[2,3] <- h0[3,2] <- 0
  > 
  > aeq(as.matrix(gchol(h0)), as.matrix(fit0$hmat))
  [1] TRUE
  > aeq(diag(gchol(h0)), diag(fit0$hmat))
  [1] TRUE
  > hinv <- solve(h0)
  > hinv[1,2:3] <-hinv[2:3,1] <- 0
  > hinv[2,3] <- hinv[3,2] <- 0
  > aeq(hinv, as.matrix(fit0$var))  
  [1] TRUE
  > 
  > # Efron/ag
  > fit0 <- coxme(Surv(zero,time, status) ~ x1 + x2 + (1|grp), data=tdata1,
  +               vfixed=theta, weight=wt, iter=0,
  +               sparse=c(2, .25))
  > aeq(apply(dt0$score,2,sum), fit0$u)
  [1] TRUE
  > aeq(as.matrix(gchol(h0)), as.matrix(fit0$hmat))
  [1] TRUE
  > aeq(diag(gchol(h0)), diag(fit0$hmat))
  [1] TRUE
  > aeq(hinv, as.matrix(fit0$var))  
  [1] TRUE
  > 
  > rm(tfit, dt0, fit0, h0, hinv, fit1, theta, zero)
  > 
  > proc.time()
     user  system elapsed 
     1.29    0.14    1.42 
Running the tests in 'tests/test4.R' failed.
Complete output:
  > library(coxme)
  Loading required package: survival
  Loading required package: bdsmatrix
  
  Attaching package: 'bdsmatrix'
  
  The following object is masked from 'package:base':
  
      backsolve
  
  > options(na.action='na.exclude', contrasts=c('contr.treatment', 'contr.poly'))
  > aeq <- function(x,y) all.equal(as.vector(x), as.vector(y))
  > 
  > #
  > # A further test of mixed sparse/ non-sparse computation
  > #  Make sure all the terms are added in and taken out properly,
  > #  when people enter and leave the risk set. 
  > #
  > # The choice of "sparse" below forces group 1 to be non-sparse and
  > #   group 2-4 to be sparse.  It will also force the coefficients to be
  > #   in 2,3,4,1 order
  > 
  > tdata1 <- data.frame(time  =c(5,4,1,1,2,2,2,2,3, 1:8), 
  +                      status=c(0,1,1,0,1,1,1,0,0, rep(1:0,4)),
  +                      x1    =c(0,1,2,0,1,1,0,1,0, rep(4:1,2)),
  +                      wt    =c(1,2,1,2,3,4,3,2,1, rep(1:2,4)),
  +                      x2    =c(1,3,5,2,3,6,4,3,1, rep(1:2,4)),
  +                      grp   =c(1,1,2,2,1,1,2,2,1, rep(3,4), rep(4,4)))
  > temp1 <-  5*tdata1$time -3  # each subject spends 4 days not at risk
  > tdata2 <- data.frame(t1 = c(rep(0, length(temp1)), temp1+4),
  + 		     t2 = c(temp1, 10*tdata1$time),
  +                      status=c(tdata1$status, rev(tdata1$status)),
  + 		     x1 = rep(tdata1$x1, 2),
  + 		     wt = rep(tdata1$wt, 2),
  + 		     x2 = rep(tdata1$x2, 2),
  + 		     grp= rep(tdata1$grp, 2),
  + 		     id = rep(1:17, 2))
  > 
  > # What's it look like?
  > #plot(c(0,80), c(0,17), type='n')
  > #segments(tdata2$t1, tdata2$id, tdata2$t2, tdata2$id)
  > #points(tdata2$t2[tdata2$status==1], tdata2$id[tdata2$status==1], pch=1)
  > 
  > theta <- .83
  > fit1 <- coxme(Surv(t1, t2, status) ~ x1 + x2 +(1|grp), data=tdata2,
  + 	      vfixed=theta, ties='breslow', 
  +               sparse=c(2, .25))
  > 
  > # This fit will complain when it tries to invert the information matrix--
  > #  ignore the complaint
  > fit2 <- coxph(Surv(t1, t2, status) ~ I(grp==2) + I(grp==3) + I(grp==4) +
  + 	      I(grp==1) + x1 + x2, data=tdata2, method='breslow',
  + 	      init=c(unlist(ranef(fit1)), fixef(fit1)),
  + 	      iter=0)
  > 
  > lp <- c(cbind(tdata2$x1, tdata2$x2) %*% fixef(fit1) + 
  + 	unlist(ranef(fit1))[c(4,1,2,3)[tdata2$grp]])
  > aeq(lp, fit1$linear)
  [1] TRUE
  > fit3 <- coxph(Surv(t1, t2, status) ~ offset(lp), data=tdata2, method='breslow')
  > 
  > aeq(fit1$loglik[3] + fit1$penalty, fit3$loglik)
  [1] TRUE
  > aeq(fit3$loglik, fit2$loglik[2])
  [1] TRUE
  > 
  > 
  > #
  > # And now, do it mostly by hand as the definitive check
  > #
  > lp2 <- fit1$linear - sum(fit1$mean * fixef(fit1))
  > temp.risk <- exp(lp2)
  > dtimes <- sort(unique(tdata2$t2[tdata2$status==1]))
  > nd <- length(dtimes)
  > denom <- double(nd)
  > xbar  <- matrix(0., nrow=nd, ncol=6)
  > xtemp <- cbind(1*(tdata2$grp==2), 1*(tdata2$grp==3), 1*(tdata2$grp==4),
  + 	       1*(tdata2$grp==1), tdata2$x1, tdata2$x2)
  > 
  > for (i in 1:nd) {
  +     who <- (tdata2$t1 <dtimes[i] & tdata2$t2 >= dtimes[i])
  +     denom[i] <- sum(temp.risk[who])
  +     xbar[i,] <- colSums(temp.risk[who] * xtemp[who,])/denom[i]
  +     }
  > temp.log <- sum(lp2[tdata2$status==1] -
  + 		log(denom[match(tdata2$t2[tdata2$status==1], dtimes)]))
  > aeq(fit1$loglik[3] + fit1$pen, temp.log)
  [1] TRUE
  > aeq(fit1$penalty, sum(unlist(ranef(fit1))^2)/(2*theta))
  [1] TRUE
  > # The linear predictors will differ by a constant, since fit2 will have
  > #   subtracted means from the factor terms.
  > aeq(diff(range(fit1$linear-fit2$linear)),0)
  [1] TRUE
  > 
  > # Score statistic
  > temp.score <- xtemp[tdata2$status==1,] - 
  + 	      xbar[match(tdata2$t2[tdata2$status==1], dtimes),]
  > aeq(colSums(temp.score) - c(unlist(ranef(fit1)),0,0)/theta, fit1$u)
  [1] TRUE
  > 
  > 
  > dt <- coxph.detail(fit2, riskmat=T)
  > aeq(dt$x, xtemp[as.numeric(dimnames(dt$x)[[1]]),])
  [1] TRUE
  > nevent <- table(tdata2$t2, tdata2$status)[,2]
  > nevent <- nevent[nevent>0]
  > aeq(dt$nevent, nevent)
  [1] TRUE
  > aeq(dt$means, xbar)
  [1] TRUE
  > aeq(colSums(temp.score), colSums(dt$score))
  [1] TRUE
  > 
  > # Check out the information matrix.
  > #  The coxme model has extra added to the diagonal, and due to
  > #  the sparseness the off-diagonal elements for vars 1-3 are zero.
  > temp1 <- as.matrix(fit1$hmat) %*% diag(sqrt(diag(fit1$hmat)))
  > temp.imat <- temp1 %*% t(temp1)
  > dt.imat <- apply(dt$imat,1:2, sum)
  > aeq(diag(temp.imat), diag(dt.imat) + c(rep(1/theta,4),0,0))
  [1] TRUE
  > aeq(temp.imat[4:6, 3], dt.imat[4:6,3])
  [1] TRUE
  > aeq(temp.imat[5:6, 4:6], dt.imat[5:6,4:6])
  [1] TRUE
  > 
  > 
  > proc.time()
     user  system elapsed 
     1.09    0.14    1.21