tests/slope4.R

In coxme: Mixed Effects Cox Models

#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)
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)
aeq(u1, u2 %*% imap)
aeq(t(map) %*% imat1 %*% map, imat2)
aeq(imat1, t(imap) %*% imat2 %*% imap)


# 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)

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))
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)))

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)
aeq(apply(dt1.1$imat, 1:2,sum) + ipen1, igchol(fit1a$hmat))
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)))



# 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)
aeq( ipen1, map %*% ipen2 %*% t(map))

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)
aeq(imat2 + ipen2, igchol(fit0b$hmat))
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)))


# 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)
aeq(fixef(fita), fixef(fitc), tol=1e-4)

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)
aeq(fitc2$log, fita$log, tol=1e-5)