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

> aeq <- function(x,y) all.equal(as.vector(x), as.vector(y))
> #
> # Ensure that a call with two terms and a call with one term but zero
> #  correlation give the same answer.  These are identical mathmatically,
> #  but follow two different indexing paths in the code.
> #
> approx <- "efron"
> set.seed(3101953)
> 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
+                )
+     }
> 
> 
> # The true model has per site coefficients of
> #    site 1:  .5 + .25*trt1
> #    site 2: 1.5 + .35*trt1
> #    site 3: 2.0 + .45*trt1
> #    site 4: 1.0 + .55*trt1
> smdata <- mkdata(150)  # 150 per site
> 
> # Variance function that sets the parameters to a fixed value
> myvar <- function(var=c(.25, .025, .078)) {
+     initialize <- function(vinit, fixed, intercept, G, X, sparse, ...) {
+         imap <- as.matrix(as.numeric(G[[1]]))
+         list(theta=NULL, imap=imap, X=X, xmap=imap+4,
+              parms=list(theta=var, fixed=rep(TRUE,8),
+                         xname=dimnames(X)[[2]], nvar=ncol(X)))
+     }
+     generate <- function(newtheta, parms) {
+         theta <- parms$theta
+         varmat <- diag(rep(theta[1:2], each=4))
+         for (i in 1:4) varmat[i,i+4] <- varmat[i+4,i] <- theta[3]
+         varmat
+     }
+     
+     wrapup <- function(newtheta, b, parms) {
+         theta <- parms$theta
+         rtheta <- list(site=c(var1=theta[1], var2=theta[2], covar=theta[3]))
+         b <- matrix(b, ncol=2)
+         dimnames(b) <- list(1:4, c("Intercept", "Slope"))
+         list(theta=rtheta, b=b)
+     }
+     out <- list(initialize=initialize, generate=generate, wrapup=wrapup)
+     class(out) <- 'coxmevar'
+     out
+     }
> 
> fit1<-  coxme(Surv(futime, status) ~ trt1 + (1 | site) + (trt1|site), smdata,
+               ties=approx)
> fit2 <- coxme(Surv(futime, status) ~ trt1 + (1+trt1 | site), smdata,
+               ties=approx, varlist=myvar(c(unlist(VarCorr(fit1)), 0)))
> 
> aeq(fit1$log, fit2$log)
[1] TRUE
> aeq(unlist(ranef(fit1)), unlist(ranef(fit2)))
[1] TRUE
> aeq(fixef(fit1), fixef(fit2))
[1] TRUE
> 
> # Same models, but fit with the start,stop part of the code
> dummy <- runif(nrow(smdata), -4, -1)  #all start times before any stop times
> fit1b <- coxme(Surv(dummy, futime, status) ~ trt1 + (1 | site) + (trt1|site), 
+                smdata, ties=approx)
> all.equal(ranef(fit1b), ranef(fit1))
[1] TRUE
> aeq(fit1$loglik, fit1b$loglik)
[1] TRUE
> 
> fit2b <- coxme(Surv(dummy, futime, status) ~ trt1 + (1+trt1 | site), smdata,
+               ties=approx, varlist=myvar(c(unlist(VarCorr(fit1)), 0)))
> all.equal(fit2b$loglik, fit2$loglik)
[1] TRUE
> all.equal(coef(fit2b), coef(fit2))
[1] TRUE
> 
> proc.time()
   user  system elapsed 
   1.35    0.06    1.42