I am trying to fit a Cox regression (package needed:survival), and have run into a problem. When I try to fit a "regular" Cox regression model with my data, I receive the error message "X matrix deemed to be singular; variable 9" (and if I remove variable 9, the problem becomes variable 8). As far as I understand the problem, this happens because too many patients with these variables have the same event (I believe in another question this was called "perfect classification")
That's why I tried to fit a Cox model with coxphf function (package of the same name), as this should take care of the problem by using "Firth’s penalized maximum likelihood bias reduction method" for the Cox regression. But this also doesn't seem to work, until I increase the "maxit" from the default 50 to 1000 and remove the "undefined" variable from the equation. But if I remove the undefined variable from my dataset (it is only 1 person), the model doesn't seem to work again.
So my question is, how can I solve this? Is it even appropriate/necessary to remove the whole variable (a therefore that 1 person) from the dataset? Probably made some very obvious mistakes, but please bear with me, since I have absolutely no background in statistics. Thank you very much already in advance.
I included the following sample data, as well as my attempts to fit the Cox model. So this is how I managed to get the model to work, by leaving out the "Undefined" variable from the model:
# load packages
library("survival")
library ("coxphf")
example<-structure(list(Pat.nr. = c(1L, 2L, 5L, 7L, 8L, 10L, 13L, 14L,
15L, 17L, 18L, 19L, 20L, 21L, 22L, 23L, 26L, 28L, 29L, 30L, 31L,
32L, 33L, 34L, 35L, 36L, 37L, 39L, 41L, 42L, 44L, 45L, 46L, 47L,
48L, 49L, 50L, 52L, 53L, 56L, 57L, 58L, 59L, 60L, 61L, 62L, 63L,
67L, 68L, 69L), OS = c(1.6, 34.6, 1.5, 35.8, 7.7, 38.6, 37.6,
8.6, 0.6, 5.7, 0.6, 43.9, 25.8, 7.3, 28.1, 43.8, 12.8, 18.5,
36.1, 43.1, 15.4, 37.6, 8.6, 2.7, 10.2, 8.1, 37.3, 25.3, 3.7,
26.1, 41.2, 5.9, 15.5, 56.8, 29.5, 52.1, 5.4, 54.8, 53.5, 16.6,
49.2, 53.8, 8.5, 56, 7.4, 28, 3.3, 38, 55.7, 0.4), Event = c(1L,
0L, 0L, 0L, 1L, 0L, 0L, 1L, 1L, 1L, 1L, 0L, 0L, 1L, 1L, 0L, 1L,
0L, 0L, 0L, 1L, 1L, 1L, 1L, 1L, 0L, 0L, 0L, 1L, 0L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 0L, 0L, 1L, 0L, 0L, 1L, 0L, 1L, 0L, 1L, 0L, 1L,
1L), Age = c(68.41, 54.9, 44.44, 64.14, 68.86, 62.93, 40.76,
31.06, 42.97, 69.16, 47.39, 60.14, 27.9, 56.57, 19.63, 47.75,
45.58, 66.22, 43.73, 45.34, 38.83, 54.46, 48.91, 70.3, 60.51,
68.55, 63.18, 55.89, 68.27, 57.25, 56.17, 60.83, 74.42, 71.3,
40.36, 50.85, 59.61, 50.14, 45.77, 19.34, 56.32, 53.38, 70.7,
55.25, 56.05, 44.06, 51.36, 69.37, 69.71, 75.44), Favorable = c(0L,
0L, 1L, 0L, 0L, 1L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 1L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 1L,
0L, 0L, 1L, 1L, 0L, 1L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L,
0L), Intermediate = c(0L, 2L, 0L, 0L, 2L, 0L, 0L, 0L, 0L, 2L,
0L, 2L, 2L, 0L, 2L, 0L, 2L, 2L, 0L, 0L, 0L, 0L, 2L, 0L, 0L, 0L,
2L, 0L, 0L, 2L, 2L, 2L, 0L, 2L, 2L, 0L, 0L, 2L, 0L, 2L, 0L, 0L,
0L, 2L, 2L, 0L, 2L, 0L, 2L, 2L), Adverse = c(1L, 0L, 0L, 1L,
0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 0L, 1L, 0L,
1L, 0L, 0L, 1L, 1L, 1L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 0L, 1L, 0L, 0L, 0L, 0L), Undefined = c(0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L)), .Names = c("Pat.nr.", "OS", "Event", "Age", "Favorable",
"Intermediate", "Adverse", "Undefined"), row.names = c(NA, 50L
), class = "data.frame")
#row&columns
n_row <- dim(example)[1]
n_col <- dim(example)[2]
#variables:
OS <- c(example[1:n_row,2])
Event <- c(example[1:n_row,3])
age <- c(example[1:n_row,4])
Favorable <- c(example[1:n_row,5])
Intermediate <- c(example[1:n_row,6])
Adverse <- c(example[1:n_row,7])
Undefined <- c(example[1:n_row,8])
# dependent and independent variables
y <- Surv(OS, Event)
x <- cbind(age, Favorable, Intermediate, Adverse, Undefined)
example <- data.frame(cbind(x,y))
#coxphf with Firth's Penalized Likelihood --> which doesn't seem to work
cox2<-coxphf(data=example,y~x, firth=TRUE, pl=TRUE, maxit=1000)
summary(cox2)
#coxphf with Firth's Penalized Likelihood (without Variable "Undefined") --> this works
cox2<-coxphf(data=example,y~age+Favorable+Intermediate+Adverse, firth=TRUE, pl=TRUE, maxit=1000)
summary(cox2)
And here, I have modified the dataset to not include the undefined variable (and the model doesn't work anymore):
example1<-structure(list(Pat.nr. = c(1L, 2L, 5L, 7L, 8L, 10L, 13L, 14L,
15L, 17L, 19L, 20L, 21L, 22L, 23L, 26L, 28L, 29L, 30L, 31L, 32L,
33L, 34L, 35L, 36L, 37L, 39L, 41L, 42L, 44L, 45L, 46L, 47L, 48L,
49L, 50L, 52L, 53L, 56L, 57L, 58L, 59L, 60L, 61L, 62L, 63L, 67L,
68L, 69L, 72L), OS = c(1.6, 34.6, 1.5, 35.8, 7.7, 38.6, 37.6,
8.6, 0.6, 5.7, 43.9, 25.8, 7.3, 28.1, 43.8, 12.8, 18.5, 36.1,
43.1, 15.4, 37.6, 8.6, 2.7, 10.2, 8.1, 37.3, 25.3, 3.7, 26.1,
41.2, 5.9, 15.5, 56.8, 29.5, 52.1, 5.4, 54.8, 53.5, 16.6, 49.2,
53.8, 8.5, 56, 7.4, 28, 3.3, 38, 55.7, 0.4, 2.8), Event = c(1L,
0L, 0L, 0L, 1L, 0L, 0L, 1L, 1L, 1L, 0L, 0L, 1L, 1L, 0L, 1L, 0L,
0L, 0L, 1L, 1L, 1L, 1L, 1L, 0L, 0L, 0L, 1L, 0L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 0L, 0L, 1L, 0L, 0L, 1L, 0L, 1L, 0L, 1L, 0L, 1L, 1L,
1L), Age = c(68.41, 54.9, 44.44, 64.14, 68.86, 62.93, 40.76,
31.06, 42.97, 69.16, 60.14, 27.9, 56.57, 19.63, 47.75, 45.58,
66.22, 43.73, 45.34, 38.83, 54.46, 48.91, 70.3, 60.51, 68.55,
63.18, 55.89, 68.27, 57.25, 56.17, 60.83, 74.42, 71.3, 40.36,
50.85, 59.61, 50.14, 45.77, 19.34, 56.32, 53.38, 70.7, 55.25,
56.05, 44.06, 51.36, 69.37, 69.71, 75.44, 71.05), Favorable = c(0L,
0L, 1L, 0L, 0L, 1L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 1L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 1L, 0L,
0L, 1L, 1L, 0L, 1L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L,
1L), Intermediate = c(0L, 2L, 0L, 0L, 2L, 0L, 0L, 0L, 0L, 2L,
2L, 2L, 0L, 2L, 0L, 2L, 2L, 0L, 0L, 0L, 0L, 2L, 0L, 0L, 0L, 2L,
0L, 0L, 2L, 2L, 2L, 0L, 2L, 2L, 0L, 0L, 2L, 0L, 2L, 0L, 0L, 0L,
2L, 2L, 0L, 2L, 0L, 2L, 2L, 0L), Adverse = c(1L, 0L, 0L, 1L,
0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 0L, 1L, 0L, 1L,
0L, 0L, 1L, 1L, 1L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
0L, 0L, 0L, 1L, 0L, 1L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L)), .Names = c("Pat.nr.",
"OS", "Event", "Age", "Favorable", "Intermediate", "Adverse"), row.names = c(NA,
50L), class = "data.frame")
#row&columns
n_row <- dim(example1)[1]
n_col <- dim(example1)[2]
#variables:
OS <- c(example1[1:n_row,2])
Event <- c(example1[1:n_row,3])
age <- c(example1[1:n_row,4])
Favorable <- c(example1[1:n_row,5])
Intermediate <- c(example1[1:n_row,6])
Adverse <- c(example1[1:n_row,7])
# dependent and independent variables
y <- Surv(OS, Event)
x <- cbind(age, Favorable, Intermediate, Adverse)
example <- data.frame(cbind(x,y))
# dependent and independent variables
y <- Surv(OS, Event)
x <- cbind(age, Favorable, Intermediate, Adverse)
example1 <- data.frame(cbind(x,y))
#coxphf with Firth's Penalized Likelihood
cox2<-coxphf(data=example,y~x, firth=TRUE, pl=TRUE, maxit=1000)
summary(cox2)