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.

> #Check fit4
> #
> # The first model has random effects of 1 + trt | inst.
> # The second form is 1 | inst/trt.  The test is to make them
> #   the same using variance matrices for the second.
> #
> # The first has 18 dummy variables i1-i9 and t1-t9, the second has
> #   g1-g18, one per group.
> # Arithmetic on the dummys shows that g1= i1-t1, g2 = t1, g3=i2-t2, g4=t2,
> #   etc.  So 
> #    var(g1) =     var(i1) + var(t1) - 2*cov(i1,t1)
> #    var(g2) =     var(t1)
> #    cov(g1,g2) =  cov(i1,t1) - var(t1)
> #
> # So if we model g, we get a covariance that is a set of 2x2 blocks
> #        a+b-c  c-b
> #        c-b    b
> # which can be written as a*mat1 + b*mat2 + c*mat3 for 3 constructed matrices
> #
> 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
> 
> contr.none <- function(n,contrasts=T) {
+         if(is.numeric(n) && length(n) == 1.)
+                 levs <- 1.:n
+         else {
+                 levs <- n
+                 n <- length(n)
+         }
+         contr <- array(0., c(n, n), list(levs, levs))
+         contr[seq(1., n^2., n + 1.)] <- 1.
+         contr
+         }
> options(contrasts=c('contr.none', 'contr.poly'))
> igchol <- function(x) {
+     dd <- diag(x)
+     ll <- as.matrix(x)
+     ll %*% diag(dd) %*% t(ll)
+     }
> 
> # The basic building blocks
> tempx <- matrix(0., nrow(simdata), ncol=18)
> for (i in 1:9) {
+     tempx[,i] <- 1*(simdata$inst==i)
+     tempx[,i+9] <- tempx[,i] * simdata$trt
+     }
> cox1<- coxph(Surv(time, status) ~ tempx + age + trt, simdata,
+              iter=0, x=T)
> dt1 <- coxph.detail(cox1)
> u1 <- apply(dt1$score, 2, sum)
> imat1 <- apply(dt1$imat, 1:2, sum)
> 
> group <- strata(simdata$inst, simdata$trt, shortlabel=T, sep='/')
> cox2 <- coxph(Surv(time, status) ~ group + age + trt, simdata,
+               iter=0, x=T)
> dt2 <- coxph.detail(cox2)
> u2 <- apply(dt2$score, 2, sum)
> imat2 <- apply(dt2$imat, 1:2, sum)
> 
> map <- matrix(0, 20,20)
> for (i in 1:9) {
+     map[i, 2*i -1] <- 1
+     map[i+9, 2*i -1] <- -1
+     map[i+9, 2*i] <- 1
+     }
> map[19,19] <- 1
> map[20,20] <- 1
> imap <- round(solve(map))  #inverse map is integers as well
> 
> aeq(u1 %*% map, u2)
[1] TRUE
> aeq(u1, u2 %*% imap)
[1] TRUE
> aeq(t(map) %*% imat1 %*% map, imat2)
[1] TRUE
> aeq(imat1, t(imap) %*% imat2 %*% imap)
[1] TRUE
> 
> 
> # Now for a fit of the first form of the model
> fit0a <- coxme(Surv(time, status) ~ age + trt + (1 +trt |inst), simdata,
+                iter=0, vfixed=c(.2, .1, .3))
> 
> aeq(fit0a$u, u1)
[1] TRUE
> 
> pen1 <- matrix(0., 20, 20)
> pen1[cbind(1:9, 1:9)] <- .2
> pen1[cbind(10:18, 10:18)] <- .3
> pen1[cbind(1:9, 10:18)] <- .1*sqrt(.2* .3)
> pen1[cbind(10:18, 1:9)] <- .1*sqrt(.2* .3)
> ipen1 <- solve(gchol(pen1)) #generalized inverse
> 
> aeq(imat1 + ipen1, igchol(fit0a$hmat))
[1] TRUE
> step1 <- solve(fit0a$hmat, fit0a$u)
> 
> # iteration 1
> fit1a <- coxme(Surv(time, status) ~ age + trt + (1 +trt |inst), simdata,
+                iter=1, vfixed=c(.2, .1, .3))
> aeq(step1, c(unlist(ranef(fit1a)), fixef(fit1a)))
[1] TRUE
> 
> cox1.1<- coxph(Surv(time, status) ~ tempx + age + trt, simdata,
+              iter=0, x=T, init=step1)
> 
> dt1.1 <- coxph.detail(cox1.1)
> aeq(apply(dt1.1$score, 2, sum) - step1 %*% ipen1, fit1a$u)
[1] TRUE
> aeq(apply(dt1.1$imat, 1:2,sum) + ipen1, igchol(fit1a$hmat))
[1] TRUE
> step2 <- solve(fit1a$hmat, fit1a$u)
> 
> #iteration 2
> fit2a <- coxme(Surv(time, status) ~ age + trt + (1 +trt |inst), simdata,
+                iter=2, vfixed=c(.2, .1, .3))
> aeq(step1+step2, c(unlist(ranef(fit2a)), fixef(fit2a)))
[1] TRUE
> 
> 
> 
> # Now try the solution using method #2
> vfix <- c(.2, .3, .1*sqrt(.2*.3))  #var, var, covar
> gname <- levels(group)
> 
> mat1 <- bdsmatrix(rep(c(1,0,0,0), 9), blocksize=rep(2,9),
+                   dimnames=list(gname, gname))
> mat2 <- bdsmatrix(rep(c(1,-1,-1,1), 9), blocksize=rep(2,9),
+                   dimnames=list(gname, gname))
> mat3 <- bdsmatrix(rep(c(-2,1,1,0), 9), blocksize=rep(2,9),
+                   dimnames=list(gname, gname))
> mat1 <- as.matrix(mat1)  #force non-sparse
> 
> pen2 <- matrix(0, 20, 20)
> pen2[1:18, 1:18]  <- vfix[1]*mat1 + vfix[2]*as.matrix(mat2) + 
+                     vfix[3]*as.matrix(mat3)
> ipen2 <- solve(gchol(pen2))
> aeq( t(map) %*% pen1 %*% map, pen2)
[1] TRUE
> aeq( ipen1, map %*% ipen2 %*% t(map))
[1] TRUE
> 
> fit0b <- coxme(Surv(time, status) ~ age + trt + (1|inst/trt), simdata,
+                iter=0, vfixed=vfix,
+               varlist=coxmeMlist(list(mat1,mat2,mat3), pdcheck=F, rescale=F))
>                
> aeq(u2, fit0b$u)
[1] TRUE
> aeq(imat2 + ipen2, igchol(fit0b$hmat))
[1] TRUE
> step1b <- solve(fit0b$hmat, fit0b$u)
> 
> # Note that step1b != step1 %*% map, or t(map) or imap or t(imap),
> #  since U and imat transform by map, but ipen by inverse-map = imap
> 
> # Iteration 1
> fit1b <- coxme(Surv(time, status) ~ age + trt + (1|inst/trt), simdata,
+                iter=1, vfixed=vfix,
+               varlist=coxmeMlist(list(mat1,mat2,mat3), pdcheck=F, rescale=F))
> aeq(step1b, c(unlist(ranef(fit1b)), fixef(fit1b)))
[1] TRUE
> 
> 
> # And now the full fit
> #  After more looks at this data, I realized that the likelihood
> # wrt the correlation term is very very flat. Tiny changes in
> # starting point or solution path move it a lot, but not the final
> # loglik.
> fita <- coxme(Surv(time, status) ~ age + trt + (1+trt | inst), simdata)
>  
> # fitb has a hard time, by way of wandering into bad solutions
> # vtemp <- VarCorr(fita)[[1]][c(1,4,2)]
> # vtemp[3] <- vtemp[3] * sqrt(vtemp[1] * vtemp[2])
> #fitb <- coxme(Surv(time, status) ~ age + trt + (1|inst/trt), simdata,
> #              vinit=vtemp,
> #              varlist=coxmeMlist(list(mat1,mat2,mat3), pdcheck=F, rescale=F,
> #                                  positive=F))
> 
> # So create our own variance function, which expects v1, v2, cor
> #  and matrices that are already in the right form
> myvar <- function(varlist) {
+     varlist <- varlist
+     init <- function(vinit, fixed, intercept, G, X, sparse, ...) {
+         ngroup <- length(G)
+         n <- nrow(G)
+     
+         if (length(vinit) ==3) theta <- vinit 
+         else theta <- c(.1, .2, 0)
+         theta[3] <- (1+theta[3])/(1-theta[3])
+         
+     
+         G <- rev(expand.nested(G))
+         imap <- as.matrix(as.numeric(G[,1]))
+         rname <- names(G)[1]
+         bname <- levels(G[[1]])
+ 
+         list(theta=log(theta), imap=imap, X=NULL, xmap=NULL, 
+              parms=list(varlist=varlist, bname=bname, rname=rname,
+                         vname=names(varlist)))
+         }
+     generate <- function(newtheta, parms) {
+         theta <- exp(newtheta)
+         theta[3] <- (theta[3]-1)/(theta[3] +1)  #correlation
+         theta[3] <- theta[3] * sqrt(theta[1]*theta[2]) #covar
+         
+         theta[1]*parms$varlist[[1]] + theta[2]*parms$varlist[[2]] +
+             theta[3] * parms$varlist[[3]]
+         }
+     wrapup <- function(newtheta, b, parms) {
+         theta <- exp(newtheta)
+         theta[3] <- (theta[3]-1)/(theta[3] +1)  #correlation
+         
+         mname <- c("Intercept","Slope")
+         tmat <- matrix(theta[c(1,3,3,2)], 2,
+                        dimnames=list(mname, mname))
+         tmat <- list(tmat)
+         names(tmat) <- "Inst"
+         list(theta=tmat, b=b)
+         }
+     out <- list(initialize=init, generate=generate, wrapup=wrapup)
+     class(out) <- 'coxmevar'
+     out
+     }
> 
> fitc <- coxme(Surv(time, status) ~ age + trt + (1|inst/trt), simdata,
+               varlist=myvar(list(mat1, mat2, mat3)))
> 
> aeq(fitc$log, fita$log, tol=1e-5)
[1] TRUE
> aeq(fixef(fita), fixef(fitc), tol=1e-4)
[1] TRUE
> 
> vtemp <- VarCorr(fita)[[1]][c(1,4,2)]
> fitc2 <- coxme(Surv(time, status) ~ age + trt + (1|inst/trt), simdata,
+               vfixed=vtemp, varlist=myvar(list(mat1, mat2, mat3)))
> aeq(unlist(ranef(fita)),  map[1:18, 1:18] %*% unlist(ranef(fitc2)), tol=1e5)
[1] TRUE
> aeq(fitc2$log, fita$log, tol=1e-5)
[1] TRUE
> 
> proc.time()
   user  system elapsed 
   6.81    0.31    7.12