Using categories as new variables?

Question

I am teaching myself R since some days and I am stuck with a cox regression analysis.

I managed to divide each of the the continuous variables I have into 2 categorical groups using the cut() function. Now I wonder if I can combine these categories seperately in a cox regression analysis.

For example: I have 3 variables (A,B,C), each consisting of 2 categories: A-,A+ & B-,B+ & C-,C+. Now if I run the coxph() with naming the variables I only get the results for the "+" categories (which as I understand is because the "-" categories are used as a reference group). However, since C- seems to have a negative effect on survival, I would be more interested in recieving the result for that category.

Also I wonder if there is a way I can define each category as a new group/variable and combine them individually to see the impact of each on the survival? Or is this unnecessary?

Edit: I'll try giving a more specific example, I hope it works :)

#example data:
test<-structure(list(A = c(8, 6, 42, 97, 55, 1, 5, 7, 55, 4), B = c(93, 9, 65, 2, 51, 89, 1, 1, 5, 62), C = c(58, 99, 100, 98, 92, 100, 99, 95, 81, 67), time = c(1.6, 34.6, 1.5, 35.8, 7.7, 38.6, 40.2, 4.7, 37.6, 8.6), event= c(1, 0, 0, 0, 1, 0, 0, 1, 0, 1)))
mydata<-as.data.frame(test)

It should look like so:

 mydata
    A  B   C time status
1   8 93  58  1.6      1
2   6  9  99 34.6      0
3  42 65 100  1.5      0
4  97  2  98 35.8      0
5  55 51  92  7.7      1
6   1 89 100 38.6      0
7   5  1  99 40.2      0
8   7  1  95  4.7      1
9  55  5  81 37.6      0
10  4 62  67  8.6      1

And concerning the questions above, this is what I have done so far:

#load survival package
library("survival")

#define variables
A <- c(mydata[1:10,1])
B <- c(mydata[1:10,2])
C <- c(mydata[1:10,3])
OS <- c(mydata[1:10,4])
Event <- c(mydata[1:10,5])

# dependent and independent variables
y <- Surv(OS, Event)
x <- cbind(A, B, C)
mydata <- data.frame(cbind(x,y))

#Cox proportional hazard model, with the "raw data"
coxph1 <- coxph(y~x,data=mydata, method="breslow")
summary(coxph1)

#categorising the variables

CA=cut(mydata$A, br=c(-1,20,101), labels = c("[A-]", "[A+]"))
CB=cut(mydata$B, br=c(-1,20,101), labels = c("[B-]", "[B+]"))
CC=cut(mydata$C, br=c(-1,96,101), labels = c("[C-]", "[C+]"))

#Cox regression, combined with cut intervals
coxph2=coxph(y~CA+CB+CC, data=mydata, method="breslow")
summary(coxph2)

And the expected output is:

coxph(formula = y ~ x, data = mydata, method = "breslow")

  n= 10, number of events= 4 

         coef  exp(coef)   se(coef)      z Pr(>|z|)
xA  0.0001443  1.0001443  0.0238329  0.006    0.995
xB  0.0104826  1.0105378  0.0211830  0.495    0.621
xC -0.0497490  0.9514682  0.0383305 -1.298    0.194

   exp(coef) exp(-coef) lower .95 upper .95
xA    1.0001     0.9999    0.9545     1.048
xB    1.0105     0.9896    0.9694     1.053
xC    0.9515     1.0510    0.8826     1.026

Concordance= 0.769  (se = 0.167 )
Rsquare= 0.29   (max possible= 0.799 )
Likelihood ratio test= 3.43  on 3 df,   p=0.33
Wald test            = 3.3  on 3 df,   p=0.3476
Score (logrank) test = 4.24  on 3 df,   p=0.2364


coxph(formula = y ~ CA + CB + CC, data = mydata, method = "breslow")

  n= 10, number of events= 4 

             coef  exp(coef)   se(coef)      z Pr(>|z|)
CA[A+] -1.036e+00  3.549e-01  1.262e+00 -0.821    0.412
CB[B+]  4.294e-01  1.536e+00  1.274e+00  0.337    0.736
CC[C+] -2.162e+01  4.095e-10  2.094e+04 -0.001    0.999

       exp(coef) exp(-coef) lower .95 upper .95
CA[A+] 3.549e-01  2.818e+00    0.0299     4.213
CB[B+] 1.536e+00  6.509e-01    0.1266    18.653
CC[C+] 4.095e-10  2.442e+09    0.0000       Inf

Concordance= 0.904  (se = 0.165 )
Rsquare= 0.542   (max possible= 0.799 )
Likelihood ratio test= 7.8  on 3 df,   p=0.05031
Wald test            = 1.15  on 3 df,   p=0.7653
Score (logrank) test = 6.42  on 3 df,   p=0.09288

Answer 1

First off, on several occasions there is unnecessarily convoluted code. The test data can be written:

mydata <- data.frame(A = c(8, 6, 42, 97, 55, 1, 5, 7, 55, 4), 
                     B = c(93, 9, 65, 2, 51, 89, 1, 1, 5, 62), 
                     C = c(58, 99, 100, 98, 92, 100, 99, 95, 81, 67), 
                     time = c(1.6, 34.6, 1.5, 35.8, 7.7, 38.6, 40.2, 4.7, 37.6, 8.6), 
                     event= c(1, 0, 0, 0, 1, 0, 0, 1, 0, 1))

Instead of OS <- c(mydata[1:10,4]) you can write OS <- mydata$time but there is no need to take them our of the mydata. As you can see, the model is the same as the first one in your question:

> (coxph1 <- coxph(Surv(time, event) ~ ., data=mydata, method="breslow"))
Call:
coxph(formula = Surv(time, event) ~ ., data = mydata, method = "breslow")

       coef exp(coef)  se(coef)     z    p
A  0.000144  1.000144  0.023833  0.01 1.00
B  0.010483  1.010538  0.021183  0.49 0.62
C -0.049749  0.951468  0.038331 -1.30 0.19

Likelihood ratio test=3.43  on 3 df, p=0.33
n= 10, number of events= 4 

Regarding your question of combining the covariates individually - the ~. means using all other covariates. you could specify them specifically with ~ A + B + C or any other combination.

Regarding changing the reference category - this is only needed with >2 categories. The meaning of the coefficient is the difference between the category and the reference category. If only 2 categories are present, than changing the reference will give the same coefficient, bit with a '-' sign. To change a reference category in a Factor use the relevel function:

 mydata$CA <- cut(mydata$A, br=c(-1,20,101), labels = c("[A-]", "[A+]"))
 mydata$CB <- cut(mydata$B, br=c(-1,20,101), labels = c("[B-]", "[B+]"))
 mydata$CC <- cut(mydata$C, br=c(-1,96,101), labels = c("[C-]", "[C+]"))

mydata$CA <- relevel(mydata$CA, 2)
> (coxph1 <- coxph(Surv(time, event) ~ CA, data=mydata, method="breslow"))
Call:
coxph(formula = Surv(time, event) ~ CA, data = mydata, method = "breslow")

        coef exp(coef) se(coef)    z    p
CA[A-] 0.559     1.749    1.158 0.48 0.63

Hope this helps :)