rbind() warning in coxph() with pspline() terms

[ThomasSoeiro]

Hi again,

Sorry to flood you with issues. I hope at least some of them are useful...

I get this warning with a model that have about 30 covariates, including psplines() terms and strata() terms:

coxph(Surv(time, status) ~ <...>, data = d, subset = group == "A")
# <...>
# Iterations: 9 outer, 26 Newton-Raphson
#      Theta= 0.985 
#      Theta= 0.000000000000000000000000000837 
# Degrees of freedom for terms=  1 13  6  4  4  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  4  1  1  1  1  1  1  1  1  1 
# Likelihood ratio test=3561  on 59 df, p=<0.0000000000000002
# n= 236595, number of events= 4566 
#    (13127 observations deleted due to missingness)
# Warning message:
# In rbind(c(linear, sqrt(lvar1), sqrt(lvar2), test2, 1, pchisq(test2,  :
#   number of columns of result is not a multiple of vector length (arg 2)

Also, I get different results if I use the subset argument (above) or if I supply data that is already subsetted (below):

coxph(Surv(time, status) ~ <...>, data = d[d$group == "A", ])
# <...>
# Iterations: 8 outer, 26 Newton-Raphson
#      Theta= 0.986 
#      Theta= 0.951 
# Degrees of freedom for terms=  1.0 13.0  6.0  4.0  4.1  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  4.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0 
# Likelihood ratio test=3564  on 59.1 df, p=<0.0000000000000002
# n= 236595, number of events= 4566 
#    (13127 observations deleted due to missingness)

Could you please confirm whether the issue is in the data/model or in the software? Thanks!

---

Here is a reproducible example for the warning. I didn't manage to create a reproducible example returning different results with the subset argument or with data that is already subsetted. I can send you some data from my study if needed.

library(survival)

time <- c(1, 33, 136, 150, 28, 16, 28, 68, 88, 31, 16, 6, 68, 46, 98, 12, 1, 96, 48, 99)
status <- c(0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0)
x <- c(1, 2, 2, 1, 1, 2, 2, 1, 1, 2, 2, 4, 3, 2, 3, 1, 2, 2, 3, 4)

coxph(Surv(time, status) ~ pspline(x, 3))
# Call:
# coxph(formula = Surv(time, status) ~ pspline(x, 3))
# 
#                       coef se(coef) se2 Chisq DF p
# pspline(x, 3), linear             0   0        1  
# pspline(x, 3), nonlin                       1     
# 
# Iterations: 10 outer, 40 Newton-Raphson
#      Theta= NaN 
# Degrees of freedom for terms= 2 
# Likelihood ratio test=1.85  on 2 df, p=0.4
# n= 20, number of events= 3 
# Warning message:
# In rbind(c(linear, sqrt(lvar1), sqrt(lvar2), test2, 1, pchisq(test2,  :
#   number of columns of result is not a multiple of vector length (arg 2)