I was sent a copy of the data, and this is what I get on a different machine:
fit <- clogit(cc ~ addContr(A) + addContr(C) + addContr(A.C) + strata(set),
data=pscc)
Warning messages:
1: In fitter(X, Y, strats, offset, init, control, weights = weights, :
Loglik converged before variable 1,2,3,4 ; beta may be infinite.
2: In coxph(formula = Surv(rep(1, 13110L), cc) ~ addContr(A) + addContr(C) + :
X matrix deemed to be singular; variable 5
fit
coef exp(coef) se(coef) z p
addContr(A)2 -0.14250 0.867 110812 -1.29e-06 1
addContr(C)2 0.00525 1.005 110812 4.74e-08 1
addContr(A.C)1-2 -0.30097 0.740 110812 -2.72e-06 1
addContr(A.C)2-1 -0.47712 0.621 110812 -4.31e-06 1
addContr(A.C)2-2 NA NA 0 NA NA
xmat <- model.matrix(fit)
svd(xmat)$d
[1] 1.932097e+02 2.700101e+01 1.624731e+01 6.049630e-15 2.031334e-15
The primary issue is that the covariates matrix is singular, having rank 3
instead of rank 5.
The coxph routine prints two warning messages that things are not good about the matrix.
Warning messages should not be ignored! The insane se(coef) values in the printed result
are an even bigger clue that the model fit is suspect. Unfortunately, some small change in
the iteration path or numerics has put this data set over the edge from being seen as rank
3 (old run) to rank 4 (new run). Moral: coxph does pretty well at detecting redundat
variables, but if you know of some it never hurts to help the routine out by removing them
before the fit.
Singularity of the X matrix in a Cox model is very difficult to detect reliably; the
current threshold is the result of long experience and experiment to give as few false
messages as possible. (The RMS package in particular used truncated power basis functions
for the splines, which lead to X matrices that look almost singular numerically, but are
not.) Setting a little less stringent threshold for declaring singularity in the cholesky
decompostion sufficies for this data set.
fit2 <- clogit(cc ~ addContr(A) + addContr(C) + addContr(A.C) + strata(set),
data=pscc, toler.chol=1e-10)
I'll certainly add this to my list of test problems that I use to tune those
constants.
Terry Therneau
On 12/11/2013 09:30 PM, Hin-Tak Leung wrote:
Here is a rather long discussion etc about freetype 2.5.2, problem with the
survival
package, and build R 2.15.x with gcc 4.8.x. Please feel free to skip forward.
- freetype 2.5.2:
the fix to cope with one of the Mac OS X's system fonts just before the release
of
freetype 2.5.1 caused a regression, crashing over one of Microsoft windows'
system fonts.
So there is a 2.5.2. There are new 2.5.2 bundles for windows & Mac OS X. The
official
win/mac binaries of R were built statically with 2+-years-old freetype with a
few known
problems. Most should upgrade/rebuild.
http://sourceforge.net/projects/outmodedbonsai/files/R/
- problem with the survival package:
Trying to re-run a vignette to get the same result as two years ago
reveal a strange change. I went and bisected it down to
r11513 and r11516 of the survival package.
-------------- r11513 --------------------
clogit(cc ~ addContr(A) + addContr(C) + addContr(A.C) + strata(set))
coef exp(coef) se(coef) z p
addContr(A)2 -0.620 0.538 0.217 -2.86 0.0043
addContr(C)2 0.482 1.620 0.217 2.22 0.0270
addContr(A.C)1-2 -0.778 0.459 0.275 -2.83 0.0047
addContr(A.C)2-1 NA NA 0.000 NA NA
addContr(A.C)2-2 NA NA 0.000 NA NA
Likelihood ratio test=26 on 3 df, p=9.49e-06 n= 13110, number of events= 3524
------------------------------------------
------------- r11516 ---------------------
clogit(cc ~ addContr(A) + addContr(C) + addContr(A.C) + strata(set))
coef exp(coef) se(coef) z p
addContr(A)2 -0.14250 0.867 110812 -1.29e-06 1
addContr(C)2 0.00525 1.005 110812 4.74e-08 1
addContr(A.C)1-2 -0.30097 0.740 110812 -2.72e-06 1
addContr(A.C)2-1 -0.47712 0.621 110812 -4.31e-06 1
addContr(A.C)2-2 NA NA 0 NA NA
Likelihood ratio test=26 on 4 df, p=3.15e-05 n= 13110, number of events= 3524
------------------------------------------
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
