# tests/twoterm.R In coxme: Mixed Effects Cox Models

```library(coxme)
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)
aeq(unlist(ranef(fit1)), unlist(ranef(fit2)))
aeq(fixef(fit1), fixef(fit2))

# 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))
aeq(fit1\$loglik, fit1b\$loglik)

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)
all.equal(coef(fit2b), coef(fit2))
```

## Try the coxme package in your browser

Any scripts or data that you put into this service are public.

coxme documentation built on July 4, 2019, 5:05 p.m.