I have a time to event dataset were I’m looking at the time until individuals perform a specific action. I’m monitoring these individuals for a certain time and if that action is not performed during the monitoring time the individual gets censored in the sense that the monitored trait is neither “yes or no”.
I was thinking that a cox proportional hazards model is perfect to analyse my data and I’ve been able to fit the model. However, my problem is that the output of the model skips some of the independent variables and there might be something that I don’t grasp when it comes to the output. Let me explain this a bit further and I’m sure someone out there have an answer to my problem.
I’m working with R. The first line of code looks the following:
coxph_interaction <-coxph(Surv(Mating, No_mating) ~ A * B, data = dat)
Where “Mating” is the dependent variable, “No_mating” is the censored data and A and B are the independent variables where I’m also interested in their interaction. Both A and B consists of high, low and metabolite which gives us the following combinations:
A B
High High
Low Low
Metabolite Metabolite
With 9 different combinations of the two (3x3=9) When I’m running the model and then check:
Summary(coxph_interaction)
I get the following output:
> summary(coxph_interaction)
Call:
coxph(formula = Surv(Mating, No_mating) ~
A * B, data = dat)
n= 187, number of events= 187
(92 observations deleted due to missingness)
coef exp(coef) se(coef) z Pr(>|z|)
A high 0.1400 1.1502 0.2784 0.503 0.61518
A low -0.4099 0.6637 0.2813 -1.457 0.14509
B high -0.2511 0.7780 0.3169 -0.792 0.42824
B low -0.9005 0.4063 0.3001 -3.001 0.00269 **
A high : B high -0.1475 0.8628 0.4509 -0.327 0.74349
A low : B high 0.5164 1.6760 0.4315 1.197 0.23134
A high : B low 0.4075 1.5030 0.4269 0.954 0.33988
A low : B low 0.3094 1.3625 0.4413 0.701 0.48327
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
exp(coef) exp(-coef) lower .95 upper .95
A high 1.1502 0.8694 0.6665 1.9852
A low 0.6637 1.5066 0.3824 1.1519
B high 0.7780 1.2854 0.4180 1.4479
B low 0.4063 2.4610 0.2257 0.7317
A high : B high 0.8628 1.1590 0.3566 2.0878
A low : B high 1.6760 0.5967 0.7195 3.9043
A high : B low 1.5030 0.6653 0.6510 3.4704
A low : B low 1.3625 0.7339 0.5738 3.2357
Concordance= 0.619 (se = 0.027 )
Likelihood ratio test= 22.58 on 8 df, p=0.004
Wald test = 20.9 on 8 df, p=0.007
Score (logrank) test = 21.88 on 8 df, p=0.005
What I can’t wrap my head around is why “metabolite” is missing. Are the other ones compared to that one? I would expect there to be:
A high
A low
A metabolite
B high
B low
B metabolite
A high : B high
A low : B high
A metabolite : B high
A high : B low
A low : B low
A metabolite : B low
A high : b metabolite
A low : B metabolite
A Metabolite : B metabolite
Could someone give me an explanation to this. I would love to have output for the metabolite A, B and interactions (9 of them) as well.
library(survival)
library(lubridate)
library(ggsurvfit)
library(gtsummary)
library(tidycmprsk)
library(condsurv)
library(survminer)
These are the libraries loaded in R.
Here's the head of my dataset head(dat)
A B ` Mating No_mating
<chr> <chr> <chr> <dbl>
1 High High NA 0
2 High High 2 1
3 High High 4 1
4 High High 2 1
5 High High 6 1
6 High High 2 1
If mating did not occur within the timeframe there's a NA, I then put "No_mating" as 0 if mating did not occur to censor the data and 1 if it did occur.
EDIT:
Made the changes in Mating from NA to last event time (9) as suggested and now I get the following output:
> summary(coxph_interaction)
Call:
coxph(formula = Surv(Mating, No_mating) ~
A * B, data = dat)
n= 279, number of events= 187
coef exp(coef) se(coef) z Pr(>|z|)
A High 0.11011 1.11640 0.27724 0.397 0.69125
A Low -0.01703 0.98312 0.27718 -0.061 0.95102
B High -0.60628 0.54537 0.31657 -1.915 0.05547 .
B Low -0.63799 0.52835 0.29368 -2.172 0.02983 *
A High : B High 0.12763 1.13613 0.44984 0.284 0.77663
A Low : B High 1.12667 3.08536 0.43139 2.612 0.00901 **
A High : B Low 0.24185 1.27361 0.42380 0.571 0.56822
A Low : B Low 0.11071 1.11707 0.43813 0.253 0.80051
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
exp(coef) exp(-coef) lower .95 upper .95
A High 1.1164 0.8957 0.6484 1.9222
A Low 0.9831 1.0172 0.5710 1.6926
B High 0.5454 1.8336 0.2932 1.0143
B Low 0.5284 1.8927 0.2971 0.9395
A High: B High 1.1361 0.8802 0.4705 2.7437
A Low : B High 3.0854 0.3241 1.3247 7.1863
A High: B Low 1.2736 0.7852 0.5550 2.9226
A Low : B Low 1.1171 0.8952 0.4733 2.6365
Concordance= 0.621 (se = 0.022 )
Likelihood ratio test= 22.92 on 8 df, p=0.003
Wald test = 23.83 on 8 df, p=0.002
Score (logrank) test = 25.14 on 8 df, p=0.001
Mating
column. With the NAs the software is just ignoring those cases. You then should have a complete data set. $\endgroup$time = 9
. Note that you could also include individuals who dropped out before the end of the study this way, if you simply use each individual's own last observation time in that column. That might not be an issue for your study, but it's critical in many other applications. $\endgroup$