R Under development (unstable) (2023-08-12 r84939 ucrt) -- "Unsuffered Consequences" Copyright (C) 2023 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. > options(na.action=na.exclude) # preserve missings > options(contrasts=c('contr.treatment', 'contr.poly')) #ensure constrast type > library(survival) > > # > # A test to exercise the "infinity" check on 2 variables > # > test3 <- data.frame(futime=1:12, fustat=c(1,0,1,0,1,0,0,0,0,0,0,0), + x1=rep(0:1,6), x2=c(rep(0,6), rep(1,6))) > > # This will produce a warning message, which is the point of the test. > # The variance is close to singular and gives different answers > # on different machines > fit3 <- coxph(Surv(futime, fustat) ~ x1 + x2, test3, iter=25) Warning message: In coxph.fit(X, Y, istrat, offset, init, control, weights = weights, : Loglik converged before variable 1,2 ; coefficient may be infinite. > > all(fit3$coef < -22) [1] TRUE > all.equal(round(fit3$log, 4),c(-6.8669, -1.7918)) [1] TRUE > > # > # Actual solution > # time 1, 12 at risk, 3 each of x1/x2 = 00, 01, 10, 11 > # time 2, 10 at risk, 2, 3, 2 , 3 > # time 5, 8 at risk, 1, 3, 1, 3 > # Let r1 = exp(beta1), r2= exp(beta2) > # loglik = -log(3 + 3r1 + 3r2 + 3 r1*r2) - log(2 + 2r1 + 3r2 + 3 r1*r2) - > # log(1 + r1 + 3r2 + 3 r1*r2) > true <- function(beta) { + r1 <- exp(beta[1]) + r2 <- exp(beta[2]) + loglik <- -log(3*(1+ r1+ r2+ r1*r2)) - log(2+ 2*r1 + 3*r2 + 3*r1*r2) - + log(1 + r1 + 3*r2 + 3*r1*r2) + loglik + } > > all.equal(fit3$loglik[2], true(fit3$coef), check.attributes=FALSE) [1] TRUE > > proc.time() user system elapsed 1.18 0.07 1.26