Using the code below, I have managed to create a univariate survival analysis from ggsurvplot using the survminer package.
I have been able to generate a univariate analysis with cox regression and plot using ggsurvplot, but when I run my multivariate analysis (that includes the covariables: age + sex + smoking_status + match_RF_HF + match_RF_AMI + match_RF_T2D + match_RF_HTA + match_RF_MI, I get the same table risk table with the same values, whilst the risk table should be different between the univariate and the multivariate analyses.
Here is the code to get the univariate analysis:
# your custom theme; unchanged
custom_theme <- function(){
theme_survminer() %+replace%
theme(
legend.background = element_rect(fill = "white", color = "black"),
plot.title=element_text(hjust=0.4)
)
}
# Apply filters to create the filtered dataset
#filtered_data <- df_ukb %>%
# filter(is.na(time_diff_MACE_months) | time_diff_MACE_months >= 0) %>% # remove time diff negative
#filter(time_diff_MACE_months <= N * 12) %>% ## remove time diff more than N follow up
#filter(time_diff_MACE_months >= 1) %>% ## remove AF proc within 1 month after AF
#filter(dateAF_plusNy <= Sys.Date()) ##remove N more than today
filtered_data <- df_ukb %>%
filter(is.na(time_diff_MACE_months) | time_diff_MACE_months >= 0) %>% # remove time diff negative
filter(time_diff_MACE_months <= N * 12) %>% ## remove time diff more than N follow up
filter(time_diff_MACE_months >= 1) %>% ## remove AF proc within 1 month after AF
filter(dateAF_plusNy <= Sys.Date())
## specific order
scoreMVPA_order <- c("no AF proc + inactive", "no AF proc + active", "AF proc + inactive", "AF proc + active")
filtered_data <-
filtered_data %>%
mutate(scoreMVPA = fct_relevel(scoreMVPA,scoreMVPA_order))
# Create the formula
filtered_data$event <- filtered_data$match_MACE
filtered_data$event_followup <- filtered_data$time_diff_MACE_months
formula1 <- as.formula("Surv(event_followup, event) ~ scoreMVPA")
# Fit Cox proportional hazards model with reference group
cox1 <- coxph(Surv(event_followup, event) ~ scoreMVPA, data = filtered_data)
cox1 <- coxph(formula1, data = filtered_data)
cox1
coxP1 <- data.frame(summary(cox1)$coefficients)[,5]
coxP1
# Generate the confidence intervals and p-values with renamed scores
coxConf1 <- data.frame(summary(cox1)$conf.int) %>%
rownames_to_column() %>%
mutate(p = coxP1,
p2 = case_when(
round(p, 3) > p ~ '=',
round(p, 3) < p ~ '=',
round(p, 3) == p ~ '='
),
p3 = ifelse(round(p, 2) == 1, T, F),
tag = paste0(rowname %>% gsub("(\\D)(\\d)", "\\1 \\2", .),
". HR ", exp.coef. %>% sprintf(fmt = "%.2f", .),
" (95% CI ",
lower..95 %>% sprintf(fmt = "%.2f", .),
"-", upper..95 %>% sprintf(fmt = "%.2f", .),
"), ",
ifelse(p < 0.001, "P<0.001",
ifelse(p >= 0.01 & p < 0.05, "P<0.05",
ifelse(p >= 0.05, paste0("P", p2, sprintf(fmt = "%.2f", round(p, 2))),
paste0("P", p2, sprintf(fmt = "%.3f", p))))))) %>%
select(tag) %>%
rename(new_variable_name = tag)
coxConf1$new_variable_name <- gsub("scoreMVPA", "", coxConf1$new_variable_name)
# Remove leading and trailing spaces
coxConf1$new_variable_name <- trimws(coxConf1$new_variable_name)
# Validate the result
coxConf1
# model to plot-
###----
fit1 <- survfit(Surv(event_followup, event) ~ scoreMVPA, data = filtered_data)
# Generate the plot with the reordered legends
uni <- ggsurvplot(fit1,
data = uni_data,
test.for.trend = TRUE,
conf.int = FALSE,
censor = FALSE,
ggtheme = custom_theme(),
legend = c(0.4, 0.25),
font.legend = c(size = 12),
legend.title = element_blank(),
legend.labs = c("no AF proc + inactive. ", unlist(coxConf1)),
xlab = "Follow-up (months)",
font.x = c(size = 12),
risk.table = T,
risk.table.side = "left",
risk.table.y.text = F,
risk.table.title = "Risk Table",
risk.table.col = 1,
ylab = "Free from MACE",
font.y = c(size = 12),
ylim = c(0, 1.05),
xlim = c(0, N * 12),
break.y.by = 0.2,
break.x.by = 6,
axes.offset = TRUE)
And the code for the multivariate analysis:
##MULTIVARIATE ANALYSIS
##covariables
# Define covariates: "age", "sex", "smoking_status", "match_RF_HF", "match_RF_AMI", "match_RF_T2D", "match_RF_HTA", "match_RF_MI"
covariates <- c("age", "sex", "smoking_status", "match_RF_HF", "match_RF_AMI", "match_RF_T2D", "match_RF_HTA", "match_RF_MI")
# Covariate formula
covariate_formula <- as.formula(paste("Surv(event_followup, event) ~ scoreMVPA +", paste(covariates, collapse = "+"))) ###I have tried to combine my covariates but I am not able to
# Apply filters in Cox proportional hazards model
cox2 <- coxph(Surv(event_followup, event) ~ age + sex + smoking_status + match_RF_HF + match_RF_AMI + match_RF_T2D + match_RF_HTA + match_RF_MI + strata(scoreMVPA),
data=filtered_data)
coxP2 <- data.frame(summary(cox2)$coefficients)[,5]
#create theme
coxConf2 <- data.frame(summary(cox2)$conf.int) %>%
rownames_to_column() %>%
mutate(p = coxP2,
p2 = case_when(
round(p, 3) > p ~ '=',
round(p, 3) < p ~ '=',
round(p, 3) == p ~ '='
),
p3 = ifelse(round(p, 2) == 1, T, F),
tag = paste0(rowname %>% gsub("(\\D)(\\d)", "\\1 \\2", .),
". HR ", exp.coef. %>% sprintf(fmt = "%.2f", .),
" (95% CI ",
lower..95 %>% sprintf(fmt = "%.2f", .),
"-", upper..95 %>% sprintf(fmt = "%.2f", .),
"), ",
ifelse(p < 0.001, "P<0.001",
ifelse(p >= 0.01 & p < 0.05, "P<0.05",
ifelse(p >= 0.05, paste0("P", p2, sprintf(fmt = "%.2f", round(p, 2))),
paste0("P", p2, sprintf(fmt = "%.3f", p))))))) %>%
select(tag) %>%
rename(new_variable_name = tag)
coxConf2$new_variable_name <- gsub("scoreMVPA", "", coxConf2$new_variable_name)
# Remove leading and trailing spaces
coxConf2$new_variable_name <- trimws(coxConf2$new_variable_name)
coxConf2 <- coxConf2 %>% slice(1:3)
##
survfit2 <- survfit(cox2)
# Generate the plot with the reordered legends
multi <- ggsurvplot(survfit2,
data = multi_data,
test.for.trend = TRUE,
conf.int = FALSE,
censor = FALSE,
ggtheme = custom_theme(),
legend = c(0.4, 0.25),
font.legend = c(size = 12),
legend.title = element_blank(),
legend.labs = c("no AF proc + inactive. ", unlist(coxConf2)),
xlab = "Follow-up (months)",
font.x = c(size = 12),
risk.table = TRUE,
risk.table.side = "left",
risk.table.y.text = FALSE,
risk.table.title = "Risk Table",
risk.table.col = 1,
ylab = "Free from MACE",
font.y = c(size = 12),
ylim = c(0, 1.05),
xlim = c(0, N * 12),
break.y.by = 0.2,
break.x.by = 6,
axes.offset = TRUE)
Here is a sample of my dataset:
structure(list(sex = c("Male", "Female", "Male", "Male", "Male",
"Male", "Male", "Male", "Male", "Male", "Female", "Male", "Male",
"Male", "Male", "Male", "Male", "Male", "Male", "Male", "Male",
"Male", "Female", "Male", "Male", "Male", "Male", "Male", "Male",
"Female", "Male", "Male", "Male", "Male", "Female", "Female",
"Male", "Female", "Female", "Male", "Male", "Female", "Female",
"Female", "Male", "Male", "Male", "Male", "Male", "Male", "Female",
"Male", "Male", "Female", "Female", "Male", "Male", "Female",
"Female", "Female", "Male", "Female", "Female", "Male", "Male",
"Female", "Male", "Female", "Male", "Male", "Male", "Male", "Male",
"Male", "Female", "Male", "Male", "Male", "Male", "Male", "Male",
"Female", "Male", "Male", "Male", "Male", "Male", "Male", "Male",
"Female", "Female", "Female", "Male", "Male", "Male", "Female",
"Male", "Male", "Female", "Male", "Male", "Female", "Male", "Male",
"Male", "Male", "Male", "Female", "Female", "Male", "Male", "Female",
"Female", "Male", "Male", "Female", "Male", "Male", "Female",
"Female", "Female", "Male", "Male", "Male", "Male", "Male", "Male",
"Male", "Female", "Male", "Male", "Female", "Male", "Female",
"Male", "Male", "Male", "Female", "Male", "Male", "Male", "Male",
"Female", "Female", "Male", "Female", "Male", "Male", "Male",
"Male", "Male", "Male", "Female", "Female", "Male", "Male", "Male",
"Male", "Female", "Male", "Female", "Male", "Male", "Female",
"Female", "Male", "Male", "Female", "Male", "Male", "Female",
"Female", "Male", "Male", "Male", "Male", "Male", "Male", "Female",
"Female", "Male", "Male", "Male", "Male", "Male", "Male", "Female",
"Male", "Male", "Female", "Male", "Male", "Female", "Female",
"Male", "Male", "Female", "Female", "Female", "Female"), age = c(41,
58, 63, 65, 62, 66, 68, 48, 68, 59, 63, 67, 64, 66, 52, 52, 69,
64, 61, 67, 65, 67, 61, 62, 58, 54, 53, 66, 67, 67, 64, 65, 64,
69, 63, 51, 67, 65, 68, 64, 65, 58, 62, 65, 67, 63, 60, 52, 69,
69, 54, 53, 55, 68, 61, 46, 67, 67, 57, 55, 61, 65, 61, 66, 48,
47, 69, 55, 58, 69, 64, 53, 65, 64, 62, 68, 64, 65, 60, 58, 54,
43, 51, 69, 60, 69, 70, 64, 59, 65, 61, 69, 67, 63, 69, 58, 55,
61, 63, 64, 65, 69, 62, 69, 61, 65, 60, 53, 41, 56, 61, 65, 65,
68, 65, 51, 69, 52, 61, 54, 44, 63, 55, 62, 62, 68, 54, 63, 41,
69, 67, 50, 65, 60, 59, 67, 66, 61, 67, 47, 63, 66, 65, 64, 57,
52, 65, 65, 60, 44, 48, 67, 43, 63, 46, 59, 63, 66, 67, 62, 62,
59, 70, 64, 65, 63, 62, 62, 42, 60, 62, 61, 64, 61, 60, 68, 64,
56, 66, 61, 69, 62, 58, 59, 60, 57, 69, 47, 69, 51, 43, 63, 63,
64, 68, 64, 60, 70, 66, 62), smoking_status = c(0, 0, 0, 0, 0,
1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0,
0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0,
0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 1, 0, 0, 0, 0), match_MACE = c(0, 0, 1, 1, 0, 0, 0, 0, 1,
1, 0, 0, 1, 1, 1, 0, 1, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0,
0, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 1, 0, 1, 0, 0,
0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 1, 0, 1, 0, 0,
1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 1, 1, 0, 1, 0,
0, 1, 0, 0, 1, 1, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 1, 1, 0, 0, 0,
0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 1, 0, 0, 0,
1, 1, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 1, 0, 0, 0, 0, 1, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1,
0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 1, 1, 0, 1, 1, 0, 0,
0, 0), time_diff_MACE_months = c(72, 72, 28.4, 59.4, 72, 72,
72, 72, 68.9, 11.4, 72, 72, 11.7, 7.8, 1.5, 72, 1.7, 16.8, 72,
72, 72, 26.1, 72, 72, 62.3, 72, 72, 72, 38.9, 72, 72, 49.8, 72,
72, 13, 26.4, 72, 72, 72, 72, 72, 18.8, 24, 49.9, 25.1, 72, 16.4,
72, 27, 72, 72, 72, 72, 72, 72, 72, 3.6, 72, 72, 72, 72, 40.8,
72, 1.6, 72, 72, 69, 3.4, 72, 28.1, 72, 72, 52.8, 72, 39.6, 72,
59.9, 72, 72, 72, 72, 72, 72, 6.8, 72, 37.9, 72, 9.3, 18.7, 5.9,
72, 37.1, 72, 72, 5.9, 72, 72, 20.4, 45.6, 68.5, 72, 72, 25.4,
72, 72, 18.5, 72, 72, 72, 10.1, 24, 72, 72, 72, 72, 72, 11.3,
72, 72, 72, 72, 41.6, 72, 72, 72, 66.6, 24.1, 72, 72, 8, 72,
1.2, 72, 72, 72, 2.1, 30.6, 72, 10, 72, 55.5, 20.9, 67, 11.7,
4, 72, 18, 14, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 7.7, 72,
72, 72, 72, 64.3, 56.8, 2.6, 72, 3.9, 72, 72, 72, 72, 72, 72,
72, 37.9, 22.6, 72, 72, 30.4, 72, 72, 72, 72, 34.5, 72, 8, 72,
24.3, 72, 13.5, 9.3, 12, 72, 3, 1.3, 72, 72, 72, 72), match_RF_stroke = c(0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0), match_RF_HF = c(0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0), match_RF_AMI = c(0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0), match_RF_T2D = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), match_RF_HTA = c(0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0), match_RF_MI = c(0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0), scoreMVPA = structure(c(1L, 1L, 2L, 1L, 3L,
2L, 2L, 2L, 2L, 2L, 1L, 2L, 2L, 1L, 1L, 2L, 1L, 2L, 1L, 1L, 2L,
1L, 2L, 2L, 1L, 3L, 2L, 1L, 3L, 2L, 1L, 2L, 2L, 2L, 1L, 1L, 2L,
2L, 2L, 1L, 2L, 1L, 1L, 2L, 1L, 2L, 2L, 4L, 2L, 1L, 1L, 1L, 2L,
2L, 1L, 1L, 2L, 2L, 1L, 2L, 1L, 2L, 3L, 2L, 2L, 1L, 1L, 1L, 2L,
1L, 2L, 2L, 4L, 1L, 1L, 1L, 1L, 2L, 1L, 2L, 4L, 2L, 2L, 1L, 1L,
1L, 1L, 2L, 2L, 2L, 2L, 3L, 2L, 1L, 2L, 2L, 2L, 2L, 1L, 1L, 2L,
1L, 1L, 2L, 2L, 1L, 1L, 2L, 1L, 1L, 2L, 2L, 1L, 2L, 2L, 1L, 1L,
1L, 2L, 2L, 2L, 3L, 4L, 2L, 1L, 2L, 3L, 1L, 2L, 2L, 2L, 2L, 1L,
2L, 2L, 2L, 2L, 1L, 3L, 1L, 2L, 2L, 1L, 2L, 1L, 2L, 2L, 1L, 1L,
2L, 4L, 2L, 1L, 1L, 2L, 1L, 1L, 2L, 1L, 2L, 1L, 1L, 3L, 1L, 1L,
2L, 2L, 4L, 2L, 1L, 1L, 2L, 1L, 2L, 1L, 1L, 1L, 2L, 1L, 1L, 2L,
2L, 2L, 2L, 2L, 1L, 1L, 2L, 1L, 1L, 2L, 2L, 1L, 2L, 2L, 1L, 1L,
1L, 1L, 3L), levels = c("no AF proc + inactive", "no AF proc + active",
"AF proc + inactive", "AF proc + active"), class = "factor")), row.names = c(NA,
-200L), class = c("tbl_df", "tbl", "data.frame"))
I hope it makes sense.
It is possible to obtain a marginal adjusted survival curve using a cox model. You could either use the ggadjustedcurves function from the survminer package, or you could use the adjustedCurves package if you also need confidence intervals to go with it.
Creating a survival curve using this method is absolutely valid, but it is impossible to get an "adjusted risk table" then. There is no such thing. It is not a Kaplan-Meier estimate anymore and as such does not use numbers at risk etc.
The only form of adjusted risk table that does exist is if you use inverse probability of treatment weighting to obtain the adjusted survival curves (as described in https://robindenz1.github.io/adjustedCurves/reference/surv_iptw_km.html). That method actually weights the number of events and number of risks, making it possible to obtain the "adjusted risk table". If you use the adjustedCurves package, that table will be included in the output object but it will still not be possible to directly plot it. To me it does not really seem all that important though.