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

> 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 
   2.32    0.21    2.48