3
$\begingroup$

I need to calculate the linear predictor of a Cox PH model by hand.

I can get continuous and binary variables to match the output of predict.coxph (specifying 'lp') but I can't seem to figure out how to calculate it for categorical variables with more than 2 levels.

My aim is to assess calibration of a published model in my own data-I only have coefficients so need to be able to do this by hand.

This post on stackoverflow describes how to calculate for continuous variables...

(https://stackoverflow.com/questions/42843239/coxph-predictions-dont-match-the-coefficients)

Any advice would be appreciated! Thanks

R example...

URL   <- "http://socserv.mcmaster.ca/jfox/Books/Companion/data/Rossi.txt"
Rossi <- read.table(URL, header=TRUE) 


summary(Rossi[,c("week", "arrest", "fin")])
 #      week           arrest        fin     
 # Min.   : 1.00   Min.   :0.0000   no :216  
 # 1st Qu.:50.00   1st Qu.:0.0000   yes:216  
 # Median :52.00   Median :0.0000            
 # Mean   :45.85   Mean   :0.2639            
 # 3rd Qu.:52.00   3rd Qu.:1.0000            
 # Max.   :52.00   Max.   :1.0000 

library(survival)     

#for binary variable
fitCPH <- coxph(Surv(week, arrest) ~ fin, data=Rossi)    #Cox-PH model
(coefCPH <- coef(fitCPH))                               # estimated coefficients
#   finyes 
#-0.3690692 

head(predict(fitCPH,type="lp"))
#[1]  0.1845346  0.1845346  0.1845346 -0.1845346  0.1845346  0.1845346

head(((as.numeric(Rossi$fin) - 1) - mean(as.numeric(Rossi$fin) - 1)) * coef(fitCPH))
#[1]  0.1845346  0.1845346  0.1845346 -0.1845346  0.1845346  0.1845346

#for categorical variable
set.seed(170981)
Rossi$categorical.example <-    as.factor(sample(1:3,nrow(Rossi),replace = TRUE))
summary(Rossi[,c("week", "arrest", "categorical.example")])
 #      week           arrest       categorical.example
 # Min.   : 1.00   Min.   :0.0000   1:156              
 # 1st Qu.:50.00   1st Qu.:0.0000   2:138              
 # Median :52.00   Median :0.0000   3:138              
 # Mean   :45.85   Mean   :0.2639                      
 # 3rd Qu.:52.00   3rd Qu.:1.0000                      
 # Max.   :52.00   Max.   :1.0000               

fitCPH2 <- coxph(Surv(week, arrest) ~ categorical.example, data=Rossi)    #Cox-PH model
(coefCPH2 <- coef(fitCPH2))  
#categorical.example2 categorical.example3 
#           -0.181790            -0.103019 

head(predict(fitCPH2,type="lp"))
#[1]  0.09098066 -0.01203832 -0.01203832  0.09098066 -0.09080938 -0.09080938

#How to calculate by hand??
$\endgroup$
1
  • $\begingroup$ It seems that neither riskRegression:::coxLP.coxph nor predict(fitCPH,type="lp") center the categorical and binary variables anymore. Any confirmation? $\endgroup$ Commented Jan 12, 2023 at 23:51

2 Answers 2

3
$\begingroup$

You need to know the how the dummy variables are generated when you fit the model. Otherwise, it is impossible to get linear predictor.

For example, suppose categorical variable X has 3 levels. Generally, 2 dummy variables will be generated. One coding scheme is as following:

                                dummy variables
       -------------------------------------------
          X level             X1          x2      
       -------------------------------------------
            1                  0           0      
            2                  1           0      
            3                  0           1      
       --------------------------------------------

Then you can get the linear predictor for subjects in different X level by apply different regression coefficients.

For your example, calculate a shift first. Total number of objects = 432. Among them 138 with categorical.example = 2, and 138 with categorical.example =3. So averages are 138/432 and 138/432. The shift is $s = -(-0.181790 \times 138/432 + -0.103019 \times 138/432) = 0.09098066$.

Then for each individual, the linear predictor is

$X_1 \beta_1 + X_2\beta_2 +s$

For categorical.example = 1, we have $X_1 = 0$ and $X_2=0$, so the linear predictor = $s = 0.09098066$

For categorical.example = 2, we have $X_1 = 1$ and $X_2=0$, so the linear predictor = $\beta_1 + s = -0.181790 + 0.09098066 = -0.9080938$

For categorical.example = 2, we have $X_1 = 1$ and $X_2=0$, so the linear predictor = $\beta_2 + s = -0.103019 + 0.09098066 = -0.01203832$

$\endgroup$
1
  • $\begingroup$ Thank you-could you explain this using the R example added to original question? Thanks in advance for your help! $\endgroup$
    – Michelle
    Commented Dec 5, 2018 at 10:34
1
$\begingroup$

Thanks @user158565 for your help, I figured it out from the dummy variable clue you posted!

Code below if anyone interested!

dv1 <- as.numeric(Rossi$categorical.example)-1    #make 0,1,2 rather than 1,2,3
dv1[dv1==2] <- 0

dv2 <- as.numeric(Rossi$categorical.example)-2    #make -1,0,1 rather than 1,2,3
dv2[dv2==-1] <- 0

meandv1  <- mean(dv1)   
meandv2  <- mean(dv2) 

head(((dv1-meandv1)*coefCPH2 [1])+((dv2-meandv2)*coefCPH2 [2]))
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.