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

> # 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.29    0.25    2.53