Running 'brfit.R' [4s]
  Running 'lmekin1.R' [3s]
  Comparing 'lmekin1.Rout' to 'lmekin1.Rout.save' ... OK
  Running 'lmekin2.R' [145s]
  Running 'nest1.R' [136s]
  Running 'nest2.R' [2s]
  Comparing 'nest2.Rout' to 'nest2.Rout.save' ... OK
  Running 'refine.R' [122s]
  Running 'refine1.R' [92s]
  Running 'refine2.R' [2s]
  Comparing 'refine2.Rout' to 'refine2.Rout.save' ... OK
  Running 'refine3.R' [2s]
  Comparing 'refine3.Rout' to 'refine3.Rout.save' ... OK
  Running 'shrink.R' [2s]
  Comparing 'shrink.Rout' to 'shrink.Rout.save' ... OK
  Running 'slope1.R' [7s]
  Comparing 'slope1.Rout' to 'slope1.Rout.save' ... OK
  Running 'slope2.R' [3s]
  Comparing 'slope2.Rout' to 'slope2.Rout.save' ... OK
  Running 'slope3.R' [84s]
  Running 'slope4.R' [5s]
  Comparing 'slope4.Rout' to 'slope4.Rout.save' ... OK
  Running 'strata1.R' [1s]
  Comparing 'strata1.Rout' to 'strata1.Rout.save' ... OK
  Running 'test0.R' [1s]
  Comparing 'test0.Rout' to 'test0.Rout.save' ... OK
  Running 'test1.R' [70s]
  Running 'test2.R' [88s]
  Running 'test3.R' [1s]
  Comparing 'test3.Rout' to 'test3.Rout.save' ... OK
  Running 'test4.R' [81s]
  Running 'twoterm.R' [1s]
  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.89    0.15    2.03 
Running the tests in 'tests/nest1.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))
  > #
  > # Very simple data set, to test out nested factors
  > #
  > # The first version of coxme put f1 before f1/f2, the current does
  > #  the opposite
  > # old: indx <- 1:10
  > indx <- c(4:9, 1:3, 10)
  > 
  > # Undo a gchol: given gchol(x), returns x
  > #
  > igchol <- function(x) {
  +     dd <- diag(x)
  +     ll <- as.matrix(x)
  +     ll %*% diag(dd) %*% t(ll)
  +     }
  > 
  > simple <- data.frame(time=c(9:3,1,1), status=c(rep(0,7),1,1),
  +                      f1=rep(1:3,3), f2=c(rep(1,6), rep(2,3)), x=1:9)
  > sfit <- coxme(Surv(time, status) ~ x + (1| f1/f2), data=simple, ties='breslow',
  +                vfixed=c(1,2), varlist=list(coxmeFull(collapse=F)), iter=0)
  > 
  > ta <- 2/9
  > tb <- -1/9
  > tc <- 4/27
  > td <- 14/81 
  > te <- 8/81 
  > tf <- -2/27
  > tg <- 2/27
  > th <- -1/27
  > ti <- -2/81
  > tj <- -4/81
  > tk <- -1/81
  > tx <- c(-3,0,3,-5,2,-3,3,-1,4,60)/9
  > itrue <- bdsmatrix(c(ta, tb, tb, tc, tg, tf, th, tf, th,
  +                          ta, tb, tf, th, tc, tg, tf, th,
  +                              ta, tf, th, tf, th, tc, tg,
  +                                  td, ti, tj, ti, tj, ti,
  +                                      te, ti, tk, ti, tk,
  +                                          td, ti, tj, ti,
  +                                              te, ti, tk,
  +                                                  td, ti, te),
  +                    blocksize=9, rmat=as.matrix(tx))
  > ibreslow <- 2*itrue[indx, indx]
  > diag(ibreslow) <- diag(ibreslow) + rep(c(1,1/2,0), c(6,3,1))
  > 
  > aeq(as.matrix(igchol(sfit$hmat)), as.matrix(ibreslow))
  [1] TRUE
  > 
  > 
  > sfit2 <- coxme(Surv(time, status) ~ x + (1| f1/f2), 
  +                data=simple, ties='breslow',  varlist=coxmeFull(collapse=F),
  +               vfixed=c(1,2), iter=0, sparse.calc=1)
  > aeq(as.matrix(igchol(sfit2$hmat)), as.matrix(ibreslow))
  [1] TRUE
  > 
  > # Check the collapse option
  > sfit2 <- coxme(Surv(time, status) ~ x + (1| f1/f2), 
  +                data=simple, ties='breslow',  varlist=coxmeFull(collapse=F),
  +               vfixed=c(1,2))
  > sfit3 <- coxme(Surv(time, status) ~ x + (1| f1/f2), 
  +                data=simple, ties='breslow',  varlist=coxmeFull(collapse=T),
  +               vfixed=c(1,2))
  > aeq(sfit2$log, sfit3$log)
  [1] TRUE
  > aeq(fixef(sfit2), fixef(sfit3))
  [1] TRUE
  > aeq(unlist(ranef(sfit3)), ranef(sfit2)[[1]] +ranef(sfit2)[[2]][c(1,1,2,2,3,3)])
  [1] TRUE
  > 
  > # Now for the Efron approx
  > sfit <- coxme(Surv(time, status) ~ x + (1| f1/f2), data=simple, 
  +               ties='efron', varlist=coxmeFull(collapse=FALSE),
  +               vfixed=c(1,2), iter=0)
  > 
  > # the matrix for the first death, where each of the last 2 obs has weight
  > #  1/2
  > tx <- c(-54, -13, 67, -132, 78, -68, 55, -4, 71, 1471)
  > i2 <- bdsmatrix(c( 60, -30, -30,  40,  20, -24,  -6, -24,  -6, 
  +                         55, -25, -20, -10,  44,  11, -20,  -5, 
  +                              55, -20, -10, -20,  -5,  44,  11, 
  +                                   48,  -8, -16,  -4, -16,  -4, 
  +                                        28,  -8,  -2,  -8,  -2, 
  +                                             48,  -4, -16,  -4, 
  +                                                  15,  -4,  -1, 
  +                                                       48,  -4,     
  + 		                                           15),
  +                    blocksize=9, rmat=as.matrix(tx))
  > iefron <- (itrue + i2/256)[indx, indx]
  > diag(iefron) <- diag(iefron) +  rep(c(1,1/2,0), c(6,3,1)) #add penalty
  > aeq(as.matrix(igchol(sfit$hmat)), as.matrix(iefron))
  [1] TRUE
  > 
  > sfit2 <- coxme(Surv(time, status) ~ x + (1| f1/f2), data=simple, ties='efron', 
  +               varlist=coxmeFull(collapse=FALSE), vfixed=c(1,2), 
  +                iter=0, sparse.calc=1)
  > aeq(as.matrix(igchol(sfit2$hmat)), as.matrix(iefron))
  [1] TRUE
  > 
  > proc.time()
     user  system elapsed 
     1.42    0.15    1.57 
Running the tests in 'tests/refine.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), ...)
  > 
  > #
  > # Same data set as slope1
  > #
  > 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
  > 
  > #
  > # Test out the refine.n code, using the simdata
  > #  A simple diagonal variance
  > #
  > # For original testing we had nsim=100, changed to 10 for a CRAN speedup
  > nsim <- 10
  > var  <- .3   #sizeable
  > 
  > set.seed(20)
  > fit1 <- coxme(Surv(time, status) ~ age + trt + (trt|inst) + strata(inst),
  +               vfixed=.3, simdata, refine.n=nsim, refine.detail=TRUE)
  > 
  > debug <- fit1$refine.detail 
  > 
  > nfrail <- length(unlist(ranef(fit1)))
  > hmatbb <- fit1$hmat[1:nfrail, 1:nfrail]
  > bhat <- unlist(ranef(fit1))  #random coefs
  > set.seed(20)
  > rdf <- coxme.control()$refine.df
  > bmat <- matrix(rnorm(nfrail*nsim), nfrail)  # replicate the simulations
  > bmat <- backsolve(hmatbb, bmat) /
  +     rep(sqrt(rchisq(nsim, rdf)/rdf), each=nfrail)
  > bmat <- bmat + bhat
  > if (!is.null(debug)) all.equal(bmat, debug$bmat)
  [1] TRUE
  > 
  > clog <- double(nsim)
  > Xmat <- scale(as.matrix(simdata[,c('age', 'trt')]), fit1$means, FALSE)
  > 
  > # Part 1, loglik for a set of nearby Cox models
  > fix.lp <- Xmat %*% fixef(fit1)
  > for (i in 1:nsim) {
  +     lp <- fix.lp + bmat[simdata$inst,i]*simdata$trt
  +     tfit <- coxph(Surv(time, status) ~ offset(lp) + strata(inst), simdata)
  +     clog[i] <- tfit$loglik
  +     }
  > if (!is.null(debug)) aeq(clog, debug$loglik)
  [1] TRUE
  > 
  > # Part 2: Taylor series for the PPL
  > b.sig <- t(bmat-bhat) %*% hmatbb  #b time sqrt(H)
  > taylor <- rowSums(b.sig^2)/2
  > 
  > temp2 <- cbind(clog-colSums(bmat^2)/.6 , fit1$log[3] - taylor)
  > m2 <- fit1$log[2]
  > errhat <- exp(temp2[,1]-m2) - exp(temp2[,2]-m2)
  > 
  > require(mvtnorm)
  Loading required package: mvtnorm
  > tdens <- dmvt(t(bmat), delta=bhat, sigma=as.matrix(solve(hmatbb)), df=rdf)
  > gdens <- -(log(2*pi) + log(.3))*nfrail/2
  > 
  > errhat <- errhat * exp(gdens-tdens)  # account for importance sampling dis
  > if (!is.null(debug)) {
  +     aeq(errhat, debug$errhat)
  +     }
  [1] TRUE
  > mtemp <- mean(errhat)
  > aeq(c(log(1+ mean(errhat)), sqrt(var(errhat)/nsim)/(1+mtemp)), fit1$refine)
  [1] TRUE
  > 
  > 
  > proc.time()
     user  system elapsed 
     1.54    0.20    1.73 
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 
     1.71    0.17    1.87 
Running the tests in 'tests/slope3.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), ...)
  > 
  > #
  > # Same data set as slope1
  > #
  > 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
  > 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)
  +     }
  > 
  > 
  > #Hammer through fit3 an iteration at a time
  > # Iteration 0, first form
  > cox1 <- coxph(Surv(time, status) ~ factor(inst)*trt + age, simdata, 
  +               iter=0, x=T)
  > indx <- c(1:9, 12:20, 11,10) #order of variables in coxph
  > cx1 <- cox1$x[,indx]
  > cox1 <- coxph(Surv(time, status) ~ cx1, simdata, iter=0, x=T)
  > dt1 <- coxph.detail(cox1)
  > u1 <- apply(dt1$score, 2, sum)
  > imat1 <- apply(dt1$imat,1:2, sum)
  > 
  > fit3a <- coxme(Surv(time, status) ~ age + trt + (1|inst) + (trt|inst),
  +                simdata, vfixed=list(.1,.3), iter=0)
  > aeq(u1, fit3a$u)
  [1] TRUE
  > pen1 <- diag(c(rep(1/.1,9), rep(1/.3,9), 0,0))
  > aeq(imat1+pen1, as.matrix(igchol(fit3a$hmat)))
  [1] TRUE
  > 
  > # Iteration 0, second form
  > idlist <- sort(outer(1:9, 0:1, paste, sep='/'))
  > mat1 <- matrix(diag(18), 18, dimnames=list(idlist, idlist))
  > mat2 <-  bdsBlock(idlist, rep(1:9, each=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),
  +                vfixed=c(.1,.3), iter=0)
  > 
  > group <- strata(simdata$inst, simdata$trt, shortlabel=T, sep='/')
  > cox2 <- coxph(Surv(time, status) ~ group + age + trt, simdata,
  +               iter=0, x=T)
  > dt2 <- coxph.detail(cox2)
  > u2 <- apply(dt2$score, 2, sum)
  > aeq(fit3b$u, u2)
  [1] TRUE
  > 
  > imat2 <- apply(dt2$imat, 1:2, sum)
  > pen2 <- matrix(0., 20,20)
  > pen2[1:18, 1:18] <- solve(.1*mat2 + .3*mat3)
  > aeq(imat2 + pen2, as.matrix(igchol(fit3b$hmat)))
  [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
  > aeq(u1 %*%map, u2)     #scores map
  [1] TRUE
  > aeq(pen2, t(map) %*% pen1 %*% map)  #penalty matrices map
  [1] TRUE
  > aeq(imat2, t(map) %*% imat1 %*% map) # Cox variances map
  [1] TRUE
  > 
  > # Now iteration 1
  > step1 <- solve(imat1 + pen1, u1)
  > step2 <- solve(imat2 + pen2, u2)
  > aeq(solve(fit3a$hmat, fit3a$u), step1)
  [1] TRUE
  > aeq(solve(fit3b$hmat, fit3b$u), step2)
  [1] TRUE
  > 
  > 
  > fit3a.1 <- coxme(Surv(time, status) ~ age + trt + (1|inst) + (trt|inst),
  +                simdata, vfixed=list(.1,.3), iter=1)
  > aeq(c(unlist(ranef(fit3a.1)), fixef(fit3a.1)), step1)
  [1] TRUE
  > 
  > fit3b.1 <-  coxme(Surv(time, status) ~ age + trt + (1|inst/trt), simdata,
  +                varlist=coxmeMlist(list(mat2, mat3), rescale=F, pdcheck=F),
  +                vfixed=c(.1,.3), iter=1)
  > aeq(c(unlist(ranef(fit3b.1)), fixef(fit3b.1)), step2)
  [1] TRUE
  > 
  > proc.time()
     user  system elapsed 
     2.04    0.10    2.15 
Running the tests in 'tests/test1.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))
  > 
  > #
  > # Now do iterations on tdata0, checking it out with coxph.detail
  > #   Note that the coxph calls give some "singular matrix" warnings
  > #   which can be ignored
  > #
  > 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))
  > theta <- .53
  > fit0 <- coxme(Surv(time, status) ~ x1 + x2 + (1|grp), data=tdata0,
  +               vfixed=theta, weight=wt, iter=0)
  > 
  > tfit <- coxph(Surv(time, status) ~ I(grp==1) + I(grp==2) + x1 + x2,
  + 	      data=tdata0, 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,0,0))
  > aeq(as.matrix(gchol(h0)), as.matrix(fit0$hmat))
  [1] TRUE
  > aeq(as.matrix(solve(fit0$var, full=F)), h0)
  [1] TRUE
  > 
  > # Now iteration 1
  > fit1 <- coxme(Surv(time, status) ~ x1 + x2 +(1|grp), data=tdata0,
  +               vfixed=theta, weight=wt, iter=1)
  > aeq(fit0$u %*% fit0$var, c(unlist(ranef(fit1)), fixef(fit1)))
  [1] TRUE
  > tfit <- coxph(Surv(time, status) ~ I(grp==1) + I(grp==2) + x1 + x2,
  + 	      data=tdata0, x=T, weight=wt, iter=0,
  + 	      init=c(unlist(ranef(fit1)), fixef(fit1)))
  > dt1 <- coxph.detail(tfit)
  > 
  > aeq(apply(dt1$score,2,sum)- c(unlist(ranef(fit1)), 0,0)/theta, fit1$u)
  [1] TRUE
  > h1 <- apply(dt1$imat,1:2,sum) + diag(c(1/theta, 1/theta,0,0))
  > aeq(as.matrix(gchol(h1)), as.matrix(fit1$hmat))
  [1] TRUE
  > aeq(as.matrix(solve(fit1$var, full=F)), h1)
  [1] TRUE
  > 
  > # And iteration 2
  > fit2 <- coxme(Surv(time, status) ~ x1 + x2 + (1|grp), data=tdata0,
  +               vfixed=theta, weight=wt, iter=2)
  > aeq(solve(fit1$hmat, fit1$u), 
  +     c(unlist(ranef(fit2)), fixef(fit2)) - 
  +     c(unlist(ranef(fit1)), fixef(fit1)))
  [1] TRUE
  > aeq(solve(fit1$hmat, fit1$u), fit1$variance %*% fit1$u)
  [1] TRUE
  > 
  > #
  > # a copy of tdata0, but with several subjects broken into multiple pieces,
  > #  exercises the "add in & take out" portions of the code
  > #
  > 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))
  > fit0 <- coxme(Surv(time1, time2, status) ~ x1 + x2 + (1|grp), data=tdata0b,
  +               vfixed=theta, weight=wt, iter=0)
  > 
  > tfit <- coxph(Surv(time1, time2, status) ~ I(grp==1) + I(grp==2) + x1 + x2,
  + 	      data=tdata0b, 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,0,0))
  > aeq(as.matrix(gchol(h0)), as.matrix(fit0$hmat))
  [1] TRUE
  > aeq(as.matrix(solve(fit0$var, full=F)), h0)
  [1] TRUE
  > 
  > # Now iteration 1
  > fit1 <- coxme(Surv(time1, time2, status) ~ x1 + x2 + (1|grp), data=tdata0b,
  +               vfixed=theta, weight=wt, iter=1)
  > aeq(fit0$u %*% fit0$var, c(unlist(ranef(fit1)), fixef(fit1)))
  [1] TRUE
  > tfit <- coxph(Surv(time1, time2, status) ~ I(grp==1) + I(grp==2) + x1 + x2,
  + 	      data=tdata0b, x=T, weight=wt, iter=0,
  + 	      init=c(unlist(ranef(fit1)), fixef(fit1)))
  > dt1 <- coxph.detail(tfit)
  > 
  > aeq(apply(dt1$score,2,sum)- c(unlist(ranef(fit1)), 0,0)/theta, fit1$u)
  [1] TRUE
  > h1 <- apply(dt1$imat,1:2,sum) + diag(c(1/theta, 1/theta,0,0))
  > aeq(as.matrix(gchol(h1)), as.matrix(fit1$hmat))
  [1] TRUE
  > aeq(as.matrix(solve(fit1$var, full=F)), h1)
  [1] TRUE
  > 
  > # And iteration 2
  > fit2 <- coxme(Surv(time1, time2, status) ~ x1 + x2 + (1|grp), data=tdata0b,
  +               vfixed=theta, weight=wt, iter=2)
  > aeq(fit1$u %*% fit1$var, c(unlist(ranef(fit2)), fixef(fit2)) - 
  +                          c(unlist(ranef(fit1)), fixef(fit1)))
  [1] TRUE
  > 
  > rm(fit0, fit1, dt0, dt1, tfit, h0, h1, theta)
  > rm(fit2)
  > 
  > proc.time()
     user  system elapsed 
     1.15    0.12    1.25 
Running the tests in 'tests/test2.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')
  > aeq <- function(x,y) all.equal(as.vector(x), as.vector(y))
  > 
  > #
  > # Similar to test1.s, but with the rats data.  This has enough groups to
  > #  force sparse matrix computations.
  > #
  > femrat <- rats[rats$sex=='f',] 
  > 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'))
  > 
  > theta <- pi/2
  > 
  > # First with no sparse
  > fit0 <- coxme(Surv(time, status) ~ rx + (1|litter), data=femrat,
  + 	      vfixed=theta, iter=0, sparse=c(100, .001))
  > tfit <- coxph(Surv(time, status) ~ factor(litter) + rx,
  +               data=femrat, x=T, 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(rep(1/theta, 50),0))
  > aeq(as.matrix(gchol(h0)), as.matrix(fit0$hmat))
  [1] TRUE
  > aeq(diag(gchol(h0)), diag(fit0$hmat))
  [1] TRUE
  > aeq(diag(fit0$var), diag(solve(h0)))
  [1] TRUE
  > 
  > 
  > # then sparse
  > fit0 <- coxme(Surv(time, status) ~ rx + (1|litter), data=femrat,
  + 	      vfixed=theta, iter=0, sparse=c(20, .1))
  > 
  > h0[1:50,1:50] <- diag(diag(h0)[1:50])
  > aeq(as.matrix(gchol(h0)), as.matrix(fit0$hmat))
  [1] TRUE
  > aeq(diag(gchol(h0)), diag(fit0$hmat))
  [1] TRUE
  > aeq(diag(fit0$var), diag(solve(h0)))
  [1] TRUE
  > 
  > # Now iteration 1
  > fit1 <- coxme(Surv(time, status) ~ rx + (1|litter), data=femrat,
  +               vfixed=theta, iter=1, sparse=c(10, .1))
  > update0 <- solve(fit0$hmat, fit0$u)
  > update0[1:50] <- update0[1:50] - mean(update0[1:50])
  > aeq(update0, c(unlist(ranef(fit1)), fixef(fit1)))
  [1] TRUE
  > tfit <- coxph(Surv(time, status) ~ factor(litter) + rx,
  +               data=femrat, x=T, iter=0,
  +               init=c(unlist(ranef(fit1)), fixef(fit1)))
  > dt1 <- coxph.detail(tfit)
  > 
  > aeq(apply(dt1$score,2,sum)- c(unlist(ranef(fit1)), 0)/theta, fit1$u)
  [1] TRUE
  > h1 <- apply(dt1$imat,1:2,sum) + diag(c(rep(1/theta, 50),0))
  > h1[1:50,1:50] <- diag(diag(h1)[1:50])
  > aeq(as.matrix(gchol(h1)), as.matrix(fit1$hmat))
  [1] TRUE
  > aeq(diag(gchol(h1)), diag(fit1$hmat))
  [1] TRUE
  > aeq(diag(fit1$var), diag(solve(h1)))
  [1] TRUE
  > 
  > 
  > # And iteration 2
  > fit2 <- coxme(Surv(time, status) ~ rx + (1|litter), data=femrat,
  +               vfixed=theta, iter=2)
  > 
  > update1 <- solve(fit1$hmat, fit1$u)
  > update1[1:50] <- update1[1:50] - mean(update1[1:50])
  > aeq(update1, c(unlist(ranef(fit2)), fixef(fit2)) -
  +     c(unlist(ranef(fit1)), fixef(fit1)))
  [1] TRUE
  > 
  > #
  > # Same computation, using a specified matrix
  > #
  > fit2b <- coxme(Surv(time, status) ~ rx + (1|litter), data=femrat,
  +               vfixed=theta, iter=2, varlist=bdsI(seq(1, 99, 2)))
  > all.equal(fit2b$u, fit2$u)
  [1] TRUE
  > all.equal(fit2b$variance, fit2$variance)
  [1] TRUE
  > all.equal(fit2b$loglik, fit2$loglik)
  [1] TRUE
  > 
  > 
  > 
  > proc.time()
     user  system elapsed 
     1.18    0.12    1.29 
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.40    0.14    1.51