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

```library(coxme)
require(nlme)
aeq <- function(x, y, ...) all.equal(as.vector(x), as.vector(y), ...)
tdata <- eortc
tdata\$center2 <- factor(tdata\$center)
tdata\$trt2 <- factor(tdata\$trt)

fit1 <- lme(y ~ trt, random= ~ 1|center2, data=tdata,
method="ML")
fit2 <- lmekin(y ~ trt + (1|center), tdata)

#I get the same loglik but slightly different coefficients.
# (The loglik is flat on top, and we are using different starting
#  estimates.)  So we need to use a toloerance for the test.
vcomp <- function(m1, m2, ...) {
temp1 <- as.numeric(VarCorr(m1)[,2])  #lme gives a matrix of text
#    temp2 <- m2\$sigma^2 * c(unlist(VarCorr(m2)), 1) #lmekin the variance struct
temp2 <- c(unlist(VarCorr(m2)), m2\$sigma^2) #lmekin the variance struct
aeq(temp1, sqrt(temp2), ...)
}

aeq(fit1\$logLik, fit2\$loglik)
aeq(fit1\$sigma, fit2\$sigma, tol=1e-4)
aeq(fixef(fit1), fixef(fit2), tol=1e-4)
vcomp(fit1, fit2, tol=1e-3)

# A second fit, using the matrix form
ucen <- sort(unique(tdata\$center))
kmat <- diag(length(ucen))
dimnames(kmat) <- list(ucen, ucen)
fit3 <- lmekin(y ~ trt + (1|center), tdata, varlist=kmat)
aeq(fit3\$coefficients, fit2\$coefficients)
aeq(unlist(VarCorr(fit3)), unlist(VarCorr(fit2)))

# Force the same coefs.
temp <- as.numeric(unclass(fit1\$modelStruct\$reStruct)[[1]])
fit3 <- lmekin(y~trt + (1|center), tdata,
vfixed= exp(2*temp))
aeq(fit1\$logLik, fit3\$loglik)
aeq(fit1\$sigma, fit3\$sigma)
aeq(fixef(fit1), fixef(fit3))
aeq(ranef(fit1), ranef(fit3), check.attributes=FALSE)
vcomp(fit1, fit3, tol=1e-7)
aeq(residuals(fit1), residuals(fit3))

# Now a model with random slopes and intercepts
fit5 <- lme(y ~ trt, random= ~ trt|center, data=tdata,
method="REML")
fit6 <- lmekin(y~ trt + (1+trt | center), data=tdata,
method="REML")
aeq(fit5\$logLik, fit6\$loglik)
aeq(fit5\$sigma, fit6\$sigma, tol=1e-4)
aeq(fixef(fit5), fixef(fit6), tol=1e-4)

#
# Use an lme data set
#
mystool <- as.data.frame(ergoStool) #get rid of contrast attributes

efit1 <-  lme(effort ~ Type, data=mystool, random= ~1|Subject,
method="ML")
efit2 <- lmekin(effort ~ Type + (1|Subject), mystool)
aeq(efit1\$logLik, efit2\$loglik)
aeq(fixef(efit1), fixef(efit2))
vcomp(efit1, efit2, tol=1e-5)
aeq(efit1\$sigma, efit2\$sigma, tol=1e-5)

efit3 <-lme(effort ~ Type, data=mystool, random= ~1|Subject,
method="REML")
efit4 <- lmekin(effort ~ Type + (1|Subject), mystool, method="REML")

aeq(efit3\$logLik, efit4\$loglik)
aeq(fixef(efit3), fixef(efit4))
vcomp(efit3, efit4, tol=1e-3)
aeq(efit3\$sigma, efit4\$sigma, tol=1e-4)
```

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