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))
> 
> # Really simple dataset -- covariate x1 is our old friend from the 
> #  validation section at the back of the book.  (Well, as yet only
> #  my online copy has the examples with weights).
> #
> 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))
> 
> # these 3 functions give the loglik, and the u/imat results for variable
> #   x1
> lfun <- function(beta, efron=T) {
+     r <- exp(beta)
+     a <- 7*r +3
+     b <- 4*r +2
+     temp1 <- 11*beta - (log(r^2 + 11*r +7) + 2*log(2*r +1))
+     if (efron) temp2 <- (10/3)*(log(a+b) + log(a*2/3 +b) + log(a/3 +b))
+     else       temp2 <- 10 * log(a+b)
+     temp1 - temp2
+     }
> 
> ufun <- function(beta, efron=T) {
+     r <- exp(beta)
+     a <- 7*r +3
+     b <- 4*r +2
+     xbar1 <- (2*r^2 + 11*r)/(r^2 + 11*r + 7)
+     xbar2 <- 11*r/(11*r +5)
+     xbar3 <- 2*r/(2*r +1)
+     xbar2b <- (7*r*2/3 + 4*r)/(a*2/3 +b)
+     xbar2c <- (7*r/3 + 4*r)/(a/3 + b)
+ 
+     temp1 <- 11 - (xbar1 + 2*xbar3)
+     if (efron) temp2 <- (10/3)*(xbar2 + xbar2b + xbar2c)
+     else       temp2 <- 10*xbar2
+     temp1 - temp2
+     }
> 
> ifun <- function(beta, efron=T) {
+     r <- exp(beta)
+     a <- 7*r +3
+     b <- 4*r +2
+     xbar1 <- (2*r^2 + 11*r)/(r^2 + 11*r + 7)
+     xbar2 <- 11*r/(11*r +5)
+     xbar3 <- 2*r/(2*r +1)
+     xbar2b <- (7*r*2/3 + 4*r)/(a*2/3 +b)
+     xbar2c <- (7*r/3 + 4*r)/(a/3 + b)
+ 
+     temp1 <- (4*r^2 + 11*r)/(r^2 + 11*r +7) - xbar1^2
+     if (efron) temp2 <- (10/3)*((xbar2 - xbar2^2) + (xbar2b - xbar2b^2) +
+                                 (xbar2c- xbar2c^2))
+     else       temp2 <- 10*(xbar2 - xbar2^2)
+     temp1 + temp2 + 2*(xbar3 - xbar3^2)
+     }
> 
> tfit <- coxme(Surv(time, status) ~ x1 + x2 + (1|grp), data=tdata0,
+               vfixed=.5, weight=wt, iter=0,
+               ties='breslow', sparse.calc=0)
> 
> aeq(tfit$loglik[1], lfun(0,F))
[1] TRUE
> aeq(tfit$u[3], ufun(0,F))
[1] TRUE
> aeq((solve(tfit$var))[3,3], ifun(0,F))
[1] TRUE
> 
> tfit1 <- coxme(Surv(time, status) ~ x1 + x2 + (1|grp), data=tdata0,
+               vfixed=.5, weight=wt, iter=0,
+               ties='breslow', sparse.calc=1)
> aeq(tfit$u, tfit1$u)
[1] TRUE
> all.equal(tfit$var, tfit1$var)
[1] TRUE
> aeq(tfit$loglik, tfit1$loglik)
[1] TRUE
> 
> # Do the matrix form, using coxmeMlist
> dmat <- diag(2)
> dimnames(dmat) <- list(1:2, 1:2)
> tfit2 <- coxme(Surv(time, status) ~ x1 + x2 + (1|grp), data=tdata0,
+               vfixed=.5, weight=wt, iter=0,
+               ties='breslow', varlist=dmat)
> aeq(tfit$u, tfit2$u)
[1] TRUE
> aeq(as.matrix(tfit$var), as.matrix(tfit2$var))
[1] TRUE
> aeq(tfit$loglik, tfit2$loglik)
[1] TRUE
> 
> #Now the Efron approx
> tfit <- coxme(Surv(time, status) ~ x1 + x2 + (1|grp), data=tdata0,
+               vfixed=.5, weight=wt, iter=0,
+               ties='efron')
> 
> aeq(tfit$loglik[3], lfun(0,T))
[1] TRUE
> aeq(tfit$u[3], ufun(0,T))
[1] TRUE
> aeq((solve(tfit$var))[3,3], ifun(0,T))
[1] TRUE
> 
> 
> #An initial value other than 0--
> 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)
> 
> aeq(tfit$loglik[3], lfun(pi,F))
[1] TRUE
> aeq(tfit$u[3], ufun(pi,F))
[1] TRUE
> aeq((solve(tfit$var))[3,3], ifun(pi,F))
[1] TRUE
> 
> # Check out that the old style code -- use of a separate random statement --
> #  still works.  One day we will drop this test and the backwards
> #  compatablility.
> #
> tfit1 <- coxme(Surv(time, status) ~ x1 + x2, data=tdata0,
+               random= ~1|grp, vfixed=.5, weight=wt, iter=0,
+               ties='breslow', init=c(pi,0), sparse.calc=1)
Warning message:
In coxme(Surv(time, status) ~ x1 + x2, data = tdata0, random = ~1 |  :
  The random argument of coxme is depreciated
> aeq(tfit$u, tfit1$u)
[1] TRUE
> all.equal(tfit$var, tfit1$var)
[1] TRUE
> aeq(tfit$loglik, tfit1$loglik)
[1] TRUE
> 
> # Efron approximation
> tfit <- coxme(Surv(time, status) ~ x1 + x2 + (1|grp), data=tdata0,
+               vfixed=.5, weight=wt, iter=0,
+               ties='efron', init=c(pi,0), sparse.calc=0)
> aeq(tfit$loglik[3], lfun(pi,T))
[1] TRUE
> aeq(tfit$u[3], ufun(pi,T))
[1] TRUE
> aeq((solve(tfit$var))[3,3], ifun(pi,T))
[1] TRUE
> 
> tfit1 <- coxme(Surv(time, status) ~ x1 + x2 + (1|grp), data=tdata0,
+               vfixed=.5, weight=wt, iter=0,
+               ties='efron', init=c(pi,0), sparse.calc=1)
> aeq(tfit$u, tfit1$u)
[1] TRUE
> all.equal(tfit$var, tfit1$var)
[1] TRUE
> aeq(tfit$loglik, tfit1$loglik)
[1] TRUE
> 
> # Use (start, stop] style input
> dummy <- rep(0, nrow(tdata0))
> tfit <- coxme(Surv(dummy, time, status) ~ x1 + x2 + (1|grp), data=tdata0,
+               vfixed=.5, weight=wt, iter=0,
+               ties='efron', init=c(pi,0), sparse.calc=0)
> aeq(tfit$loglik[3], lfun(pi,T))
[1] TRUE
> aeq(tfit$u[3], ufun(pi,T))
[1] TRUE
> aeq((solve(tfit$var))[3,3], ifun(pi,T))
[1] TRUE
> 
> tfit1 <- coxme(Surv(dummy, time, status) ~ x1 + x2 + (1|grp), data=tdata0,
+               vfixed=.5, weight=wt, iter=0,
+               ties='efron', init=c(pi,0), sparse.calc=1)
> aeq(tfit$u, tfit1$u)
[1] TRUE
> all.equal(tfit$var, tfit1$var)
[1] TRUE
> aeq(tfit$loglik, tfit1$loglik)
[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))
> tfit <- coxme(Surv(time1, time2, status) ~ x1 + x2 + (1|grp), data=tdata0b,
+               vfixed=.5, weight=wt, iter=0,
+               ties='efron', init=c(pi,0), sparse.calc=0)
> aeq(tfit$loglik[3], lfun(pi,T))
[1] TRUE
> aeq(tfit$u[3], ufun(pi,T))
[1] TRUE
> aeq((solve(tfit$var))[3,3], ifun(pi,T))
[1] TRUE
> 
> tfit1 <- coxme(Surv(time1, time2, status) ~ x1 + x2 + (1|grp), data=tdata0b,
+               vfixed=.5, weight=wt, iter=0,
+               ties='efron', init=c(pi,0), sparse.calc=1)
> aeq(tfit$u, tfit1$u)
[1] TRUE
> all.equal(tfit$var, tfit1$var)
[1] TRUE
> aeq(tfit$loglik, tfit1$loglik)
[1] TRUE
> 
> tfit <- coxme(Surv(time1, time2, status) ~ x1 + x2 + (1|grp), data=tdata0b,
+               vfixed=.5, weight=wt, iter=0,
+               ties='breslow', init=c(pi,0), sparse.calc=0)
> aeq(tfit$loglik[3], lfun(pi,F))
[1] TRUE
> aeq(tfit$u[3], ufun(pi,F))
[1] TRUE
> aeq((solve(tfit$var))[3,3], ifun(pi,F))
[1] TRUE
> 
> tfit1 <- coxme(Surv(time1, time2, status) ~ x1 + x2 + (1|grp), data=tdata0b,
+               vfixed=.5, weight=wt, iter=0,
+               ties='breslow', init=c(pi,0), sparse.calc=1)
> aeq(tfit$u, tfit1$u)
[1] TRUE
> all.equal(tfit$var, tfit1$var)
[1] TRUE
> aeq(tfit$loglik, tfit1$loglik)
[1] TRUE
> 
> rm (rcnt, tfit, tfit1, dummy)
> 
> proc.time()
   user  system elapsed 
   1.73    0.29    2.03