R Under development (unstable) (2024-03-03 r86036 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), ...)
> 
> #
> # 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.32    0.09    1.42