R Under development (unstable) (2024-08-21 r87038 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2024 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

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

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

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

> 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.78    0.15    1.87