R Under development (unstable) (2024-03-03 r86036 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. > 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)) > > # > # Test of fitting random slopes > # > # Simulation data with 9 institutions, strong age effects > # and a random treatment effect > # > 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 > > coef0 <- matrix(0., 2,9) > for (i in 1:9) { + fit0 <- coxph(Surv(time, status) ~ age + trt, simdata, + subset=(inst==i)) + coef0[,i] <- fit0$coef + } > > # Several of these fits will differ in the last few digits on a > # 64 bit vs 32 bit Intel processer. The loglike is very > # flat on top so tiny changes in the compute path lead to a small > # change in the final solution. Hence the "digits" argument. > fit0 <- coxph(Surv(time, status) ~ age + trt, simdata) > fit1 <- coxme(Surv(time, status) ~ age + trt + (1|inst), simdata) > print(fit1, rcoef=TRUE, digits=4) Cox mixed-effects model fit by maximum likelihood Data: simdata events, n = 1504, 2376 Iterations= 8 35 NULL Integrated Fitted Log-likelihood -10892.32 -10826.57 -10819.37 Chisq df p AIC BIC Integrated loglik 131.5 3.00 0 125.5 109.55 Penalized loglik 145.9 7.38 0 131.1 91.88 Model: Surv(time, status) ~ age + trt + (1 | inst) Fixed and penalized coefficients coef exp(coef) se(coef) z p age 0.029948 1.0304 0.003023 9.91 0.0e+00 trt -0.249769 0.7790 0.051899 -4.81 1.5e-06 inst.1 -0.122696 0.8845 0.078439 inst.2 -0.007691 0.9923 0.075521 inst.3 -0.043795 0.9572 0.074255 inst.4 -0.120717 0.8863 0.073124 inst.5 0.032154 1.0327 0.070791 inst.6 0.091249 1.0955 0.069919 inst.7 -0.017063 0.9831 0.068940 inst.8 0.024342 1.0246 0.067913 inst.9 0.164217 1.1785 0.066641 Random effects Group Variable Std Dev Variance inst Intercept 0.11329 0.01283 > > fit2 <- coxme(Surv(time, status) ~ age + trt + (1|inst/trt), simdata) > print(fit2, rcoef=TRUE, digits=3) Cox mixed-effects model fit by maximum likelihood Data: simdata events, n = 1504, 2376 Iterations= 14 60 NULL Integrated Fitted Log-likelihood -10892.32 -10814.83 -10796 Chisq df p AIC BIC Integrated loglik 155 4.0 0 147 125.7 Penalized loglik 193 14.2 0 164 88.8 Model: Surv(time, status) ~ age + trt + (1 | inst/trt) Fixed and penalized coefficients coef exp(coef) se(coef) z p age 0.030902 1.031 0.00304 10.18 0.0000 trt -0.303468 0.738 0.10867 -2.79 0.0052 inst/trt.1/0 -0.037941 0.963 0.12209 inst/trt.1/1 -0.262439 0.769 0.13045 inst/trt.2/0 0.093971 1.099 0.11654 inst/trt.2/1 -0.124835 0.883 0.12498 inst/trt.3/0 0.111253 1.118 0.11498 inst/trt.3/1 -0.220075 0.802 0.12246 inst/trt.4/0 -0.123843 0.884 0.11411 inst/trt.4/1 -0.155038 0.856 0.12008 inst/trt.5/0 0.043624 1.045 0.11024 inst/trt.5/1 0.027408 1.028 0.11527 inst/trt.6/0 0.154576 1.167 0.10811 inst/trt.6/1 0.042069 1.043 0.11446 inst/trt.7/0 -0.035562 0.965 0.10859 inst/trt.7/1 -0.000864 0.999 0.11119 inst/trt.8/0 -0.149839 0.861 0.10780 inst/trt.8/1 0.227553 1.256 0.10864 inst/trt.9/0 -0.056239 0.945 0.10533 inst/trt.9/1 0.466220 1.594 0.10656 inst.1 -0.001390 0.999 0.01364 inst.2 -0.000143 1.000 0.01364 inst.3 -0.000504 0.999 0.01364 inst.4 -0.001290 0.999 0.01364 inst.5 0.000329 1.000 0.01364 inst.6 0.000910 1.001 0.01364 inst.7 -0.000169 1.000 0.01364 inst.8 0.000360 1.000 0.01364 inst.9 0.001897 1.002 0.01364 Random effects Group Variable Std Dev Variance inst/trt (Intercept) 0.201142 0.040458 inst (Intercept) 0.013682 0.000187 > > # And so will this one > fit3 <- coxme(Surv(time, status) ~ age + trt + (1|inst) + (trt|inst),simdata) > print(fit3, rcoef=TRUE, digits=3) Cox mixed-effects model fit by maximum likelihood Data: simdata events, n = 1504, 2376 Iterations= 25 104 NULL Integrated Fitted Log-likelihood -10892.32 -10812.96 -10798.58 Chisq df p AIC BIC Integrated loglik 159 4.0 0 151 129 Penalized loglik 188 10.8 0 166 109 Model: Surv(time, status) ~ age + trt + (1 | inst) + (trt | inst) Fixed and penalized coefficients coef exp(coef) se(coef) z p age 0.03073 1.031 0.00303 10.14 0.0000 trt -0.30570 0.737 0.10695 -2.86 0.0043 inst.1 -0.02143 0.979 0.05618 inst.2 0.02234 1.023 0.05563 inst.3 0.02398 1.024 0.05541 inst.4 -0.04352 0.957 0.05443 inst.5 0.01729 1.017 0.05414 inst.6 0.05591 1.058 0.05403 inst.7 -0.01037 0.990 0.05348 inst.8 -0.04290 0.958 0.05284 inst.9 -0.00131 0.999 0.05249 inst.1:trt -0.29599 0.744 0.16082 inst.2:trt -0.15830 0.854 0.15333 inst.3:trt -0.26950 0.764 0.15102 inst.4:trt -0.13695 0.872 0.14825 inst.5:trt 0.02349 1.024 0.14239 inst.6:trt 0.00704 1.007 0.14141 inst.7:trt 0.01545 1.016 0.13845 inst.8:trt 0.29617 1.345 0.13549 inst.9:trt 0.51860 1.680 0.13306 Random effects Group Variable Std Dev Variance inst Intercept 0.06266 0.00393 inst trt 0.27847 0.07754 > > fit4 <- coxme(Surv(time, status) ~ age + trt + (1 +trt |inst), simdata) > > sfit0 <- coxph(Surv(time, status) ~ age + trt + strata(inst), simdata) > sfit1 <- coxme(Surv(time, status) ~ age + trt + (trt|inst) + strata(inst), + simdata) > print(sfit1, rcoef=TRUE, digits=4) Cox mixed-effects model fit by maximum likelihood Data: simdata events, n = 1504, 2376 Iterations= 7 31 NULL Integrated Fitted Log-likelihood -7636.963 -7564.201 -7553.573 Chisq df p AIC BIC Integrated loglik 145.5 3.0 0 139.5 123.6 Penalized loglik 166.8 8.3 0 150.2 106.1 Model: Surv(time, status) ~ age + trt + (trt | inst) + strata(inst) Fixed and penalized coefficients coef exp(coef) se(coef) z p age 0.03097 1.0315 0.003047 10.16 0.0000 trt -0.31038 0.7332 0.116564 -2.66 0.0077 inst.1:trt -0.23574 0.7900 0.190686 inst.2:trt -0.23583 0.7899 0.179152 inst.3:trt -0.34250 0.7100 0.176387 inst.4:trt -0.03709 0.9636 0.174728 inst.5:trt -0.01858 0.9816 0.165983 inst.6:trt -0.09182 0.9123 0.163471 inst.7:trt 0.03376 1.0343 0.162333 inst.8:trt 0.39190 1.4798 0.159318 inst.9:trt 0.53589 1.7090 0.155067 Random effects Group Variable Std Dev Variance inst trt 0.31042 0.09636 > > # Check that the start,stop code does the same > dummy <- runif(nrow(simdata), -4, -1) #all start times before first event > fit4b <- coxme(Surv(dummy, time, status) ~ age + trt + (1 +trt |inst), simdata) > all.equal(fit4b$loglik, fit4$loglik) [1] TRUE > all.equal(fit4b$coef, fit4$coef, tolerance=1e-7) # different order of internal [1] TRUE > # sums => tiny difference > > #Comparison plot > y <- cbind(slope, NA, coef0[2,], fixef(sfit1)[2] + unlist(ranef(sfit1)), + fixef(fit3)[2] + ranef(fit3)[[2]], + fixef(fit4)[2] + ranef(fit4)[[1]][,2]) > matplot(c(1, 1.5, 2:5), t(y), type='b', xaxt='n', xlab="Simulation", + ylab="Treatment coefficient", lty=1) > axis(1, 1:5, c("Sim", "Separate", "Strata", "Uncor", "Corr")) > > > # > # Now compute some things exactly > # > 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) + } > > # For fit2 > vtemp <- unlist(VarCorr(fit2)) > names(vtemp) <- names(VarCorr(fit2)) > fit2a <- coxme(Surv(time, status) ~ age + trt + (1|inst/trt), simdata, + iter=0, vfixed=vtemp) > temp <- strata(simdata$inst, simdata$trt, sep='/', shortlabel=TRUE) > cfit <- coxph(Surv(time, status) ~ factor(temp) +factor(inst) +age + trt, + simdata, iter=0, x=T) > dt2 <- coxph.detail(cfit) > u2 <- apply(dt2$score, 2, sum) > aeq(u2, fit2a$u) [1] TRUE > imat2 <- apply(dt2$imat, 1:2, sum) + diag(c(rep(1/vtemp, c(18,9)),0,0)) > aeq(imat2, as.matrix(igchol(fit2a$hmat))) [1] TRUE > > # For fit3 > vtemp <- as.vector(unlist(VarCorr(fit3))) #name not needed > fit3a <- coxme(Surv(time, status) ~ age + trt + (1|inst) + (trt|inst), + simdata, iter=0, vfixed=as.list(vtemp)) > cfit <- coxph(Surv(time, status) ~ factor(inst) * trt + age, simdata, + iter=0, x=T) > dt3 <- coxph.detail(cfit) > u3 <- apply(dt3$score, 2, sum) > indx <- c(1:9, 12:20, 11, 10) > aeq(u3[indx], fit3a$u) [1] TRUE > imat2 <- apply(dt3$imat, 1:2, sum)[indx,indx] + + diag(c(rep(1/vtemp, c(9,9)),0,0)) > aeq(imat2, as.matrix(igchol(fit3a$hmat))) [1] TRUE > > fit3b <- coxme(Surv(time, status) ~ age + trt + (trt|inst) +(1|inst), + simdata, iter=0, vfixed=as.list(rev(vtemp))) > aeq(fit3a$u, fit3b$u) [1] TRUE > aeq(fit3b$imat, fit3b$imat) [1] TRUE > > #For sfit1 > vtemp <- .0966 > fit <- coxme(Surv(time, status) ~ age + trt + strata(inst) + (trt|inst), + simdata, iter=0, vfixed=vtemp) > cfit <- coxph(Surv(time, status) ~ factor(inst):trt + trt+ age+ strata(inst), + simdata, iter=0, x=T) > dt3 <- coxph.detail(cfit) > u3 <- apply(dt3$score, 2, sum) > indx <- c(3:11,2,1) > aeq(u3[indx], fit$u) [1] TRUE > > imat3 <- apply(dt3$imat,1:2, sum)[indx,indx] + diag(c(rep(1/vtemp,9),0,0)) > aeq(imat3, as.matrix(igchol(fit$hmat))) [1] TRUE > > > proc.time() user system elapsed 7.06 0.23 7.28