tests/test0.R In coxme: Mixed Effects Cox Models

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

# Really simple dataset -- covariate x1 is our old friend from the
#  validation section at the back of the book.  (Well, as yet only
#  my online copy has the examples with weights).
#
tdata0 <- data.frame(time  =c(5,4,1,1,2,2,2,2,3),
status=c(0,1,1,0,1,1,1,0,0),
x1    =c(0,1,2,0,1,1,0,1,0),
wt    =c(1,2,1,2,3,4,3,2,1),
x2    =c(1,3,5,2,3,6,4,3,1),
grp   =c(1,1,2,2,1,1,2,2,1))

# these 3 functions give the loglik, and the u/imat results for variable
#   x1
lfun <- function(beta, efron=T) {
r <- exp(beta)
a <- 7*r +3
b <- 4*r +2
temp1 <- 11*beta - (log(r^2 + 11*r +7) + 2*log(2*r +1))
if (efron) temp2 <- (10/3)*(log(a+b) + log(a*2/3 +b) + log(a/3 +b))
else       temp2 <- 10 * log(a+b)
temp1 - temp2
}

ufun <- function(beta, efron=T) {
r <- exp(beta)
a <- 7*r +3
b <- 4*r +2
xbar1 <- (2*r^2 + 11*r)/(r^2 + 11*r + 7)
xbar2 <- 11*r/(11*r +5)
xbar3 <- 2*r/(2*r +1)
xbar2b <- (7*r*2/3 + 4*r)/(a*2/3 +b)
xbar2c <- (7*r/3 + 4*r)/(a/3 + b)

temp1 <- 11 - (xbar1 + 2*xbar3)
if (efron) temp2 <- (10/3)*(xbar2 + xbar2b + xbar2c)
else       temp2 <- 10*xbar2
temp1 - temp2
}

ifun <- function(beta, efron=T) {
r <- exp(beta)
a <- 7*r +3
b <- 4*r +2
xbar1 <- (2*r^2 + 11*r)/(r^2 + 11*r + 7)
xbar2 <- 11*r/(11*r +5)
xbar3 <- 2*r/(2*r +1)
xbar2b <- (7*r*2/3 + 4*r)/(a*2/3 +b)
xbar2c <- (7*r/3 + 4*r)/(a/3 + b)

temp1 <- (4*r^2 + 11*r)/(r^2 + 11*r +7) - xbar1^2
if (efron) temp2 <- (10/3)*((xbar2 - xbar2^2) + (xbar2b - xbar2b^2) +
(xbar2c- xbar2c^2))
else       temp2 <- 10*(xbar2 - xbar2^2)
temp1 + temp2 + 2*(xbar3 - xbar3^2)
}

tfit <- coxme(Surv(time, status) ~ x1 + x2 + (1|grp), data=tdata0,
vfixed=.5, weight=wt, iter=0,
ties='breslow', sparse.calc=0)

aeq(tfit\$loglik[1], lfun(0,F))
aeq(tfit\$u[3], ufun(0,F))
aeq((solve(tfit\$var))[3,3], ifun(0,F))

tfit1 <- coxme(Surv(time, status) ~ x1 + x2 + (1|grp), data=tdata0,
vfixed=.5, weight=wt, iter=0,
ties='breslow', sparse.calc=1)
aeq(tfit\$u, tfit1\$u)
all.equal(tfit\$var, tfit1\$var)
aeq(tfit\$loglik, tfit1\$loglik)

# Do the matrix form, using coxmeMlist
dmat <- diag(2)
dimnames(dmat) <- list(1:2, 1:2)
tfit2 <- coxme(Surv(time, status) ~ x1 + x2 + (1|grp), data=tdata0,
vfixed=.5, weight=wt, iter=0,
ties='breslow', varlist=dmat)
aeq(tfit\$u, tfit2\$u)
aeq(as.matrix(tfit\$var), as.matrix(tfit2\$var))
aeq(tfit\$loglik, tfit2\$loglik)

#Now the Efron approx
tfit <- coxme(Surv(time, status) ~ x1 + x2 + (1|grp), data=tdata0,
vfixed=.5, weight=wt, iter=0,
ties='efron')

aeq(tfit\$loglik[3], lfun(0,T))
aeq(tfit\$u[3], ufun(0,T))
aeq((solve(tfit\$var))[3,3], ifun(0,T))

#An initial value other than 0--
tfit <- coxme(Surv(time, status) ~ x1 + x2 + (1|grp), data=tdata0,
vfixed=.5, weight=wt, iter=0,
ties='breslow', init=c(pi,0), sparse.calc=0)

aeq(tfit\$loglik[3], lfun(pi,F))
aeq(tfit\$u[3], ufun(pi,F))
aeq((solve(tfit\$var))[3,3], ifun(pi,F))

# Check out that the old style code -- use of a separate random statement --
#  still works.  One day we will drop this test and the backwards
#  compatablility.
#
tfit1 <- coxme(Surv(time, status) ~ x1 + x2, data=tdata0,
random= ~1|grp, vfixed=.5, weight=wt, iter=0,
ties='breslow', init=c(pi,0), sparse.calc=1)
aeq(tfit\$u, tfit1\$u)
all.equal(tfit\$var, tfit1\$var)
aeq(tfit\$loglik, tfit1\$loglik)

# Efron approximation
tfit <- coxme(Surv(time, status) ~ x1 + x2 + (1|grp), data=tdata0,
vfixed=.5, weight=wt, iter=0,
ties='efron', init=c(pi,0), sparse.calc=0)
aeq(tfit\$loglik[3], lfun(pi,T))
aeq(tfit\$u[3], ufun(pi,T))
aeq((solve(tfit\$var))[3,3], ifun(pi,T))

tfit1 <- coxme(Surv(time, status) ~ x1 + x2 + (1|grp), data=tdata0,
vfixed=.5, weight=wt, iter=0,
ties='efron', init=c(pi,0), sparse.calc=1)
aeq(tfit\$u, tfit1\$u)
all.equal(tfit\$var, tfit1\$var)
aeq(tfit\$loglik, tfit1\$loglik)

# Use (start, stop] style input
dummy <- rep(0, nrow(tdata0))
tfit <- coxme(Surv(dummy, time, status) ~ x1 + x2 + (1|grp), data=tdata0,
vfixed=.5, weight=wt, iter=0,
ties='efron', init=c(pi,0), sparse.calc=0)
aeq(tfit\$loglik[3], lfun(pi,T))
aeq(tfit\$u[3], ufun(pi,T))
aeq((solve(tfit\$var))[3,3], ifun(pi,T))

tfit1 <- coxme(Surv(dummy, time, status) ~ x1 + x2 + (1|grp), data=tdata0,
vfixed=.5, weight=wt, iter=0,
ties='efron', init=c(pi,0), sparse.calc=1)
aeq(tfit\$u, tfit1\$u)
all.equal(tfit\$var, tfit1\$var)
aeq(tfit\$loglik, tfit1\$loglik)

#
# a copy of tdata0, but with several subjects broken into multiple pieces,
#  exercises the "add in & take out" portions of the code
#
rcnt <- c(3,2,1,1,1,2,2,1,3)  # rep count
tdata0b <- data.frame(time2  = c(3,4,5, 2,4, 1,1,2, 1,2, .5,2, 2, 1,2,3),
time1  = c(0,3,4, 0,2, 0,0,0, 0,1, 0,.5, 0, 0,1,2),
status = c(0,0,0, 0,1, 1,0,1, 0,1, 0, 1, 0, 0,0,0),
x1     = rep(tdata0\$x1, rcnt),
wt     = rep(tdata0\$wt, rcnt),
x2     = rep(tdata0\$x2, rcnt),
grp    = rep(tdata0\$grp, rcnt))
tfit <- coxme(Surv(time1, time2, status) ~ x1 + x2 + (1|grp), data=tdata0b,
vfixed=.5, weight=wt, iter=0,
ties='efron', init=c(pi,0), sparse.calc=0)
aeq(tfit\$loglik[3], lfun(pi,T))
aeq(tfit\$u[3], ufun(pi,T))
aeq((solve(tfit\$var))[3,3], ifun(pi,T))

tfit1 <- coxme(Surv(time1, time2, status) ~ x1 + x2 + (1|grp), data=tdata0b,
vfixed=.5, weight=wt, iter=0,
ties='efron', init=c(pi,0), sparse.calc=1)
aeq(tfit\$u, tfit1\$u)
all.equal(tfit\$var, tfit1\$var)
aeq(tfit\$loglik, tfit1\$loglik)

tfit <- coxme(Surv(time1, time2, status) ~ x1 + x2 + (1|grp), data=tdata0b,
vfixed=.5, weight=wt, iter=0,
ties='breslow', init=c(pi,0), sparse.calc=0)
aeq(tfit\$loglik[3], lfun(pi,F))
aeq(tfit\$u[3], ufun(pi,F))
aeq((solve(tfit\$var))[3,3], ifun(pi,F))

tfit1 <- coxme(Surv(time1, time2, status) ~ x1 + x2 + (1|grp), data=tdata0b,
vfixed=.5, weight=wt, iter=0,
ties='breslow', init=c(pi,0), sparse.calc=1)
aeq(tfit\$u, tfit1\$u)
all.equal(tfit\$var, tfit1\$var)
aeq(tfit\$loglik, tfit1\$loglik)

rm (rcnt, tfit, tfit1, dummy)
```

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.