CODEWITHSUNDEEP
rstatisticsdata-visualization
Noob3000
Noob3000

Reputation: 53

Data Interpretation and Visualisation

I am trying to fit a Cox Proportional Hazard model using mother's age and infant's gender as covariates.

I am having problems with visually representing the data in R.

I'll post my code so far:

These are my packages:

#Getting started:

# load libraries
pkgTest <- function(pkg){
  new.pkg <- pkg[!(pkg %in% installed.packages()[,  "Package"])]
  if (length(new.pkg)) 
    install.packages(new.pkg,  dependencies = TRUE)
  sapply(pkg,  require,  character.only = TRUE)
}

lapply(c("survival", "eha", "tidyverse", "ggfortify", "stargazer"),  pkgTest)

Here is the data I am using and trying to run a test on:


data(infants)

imr <- with(infants, Surv(enter, exit, event))

cox <- coxph(imr ~ sex + age, data = infants)
summary(cox)
drop1(cox, test = "Chisq")
stargazer(cox, type = "text")


cox_fit <- survfit(cox)
autoplot(cox_fit)

newdat <- with(infants, 
               data.frame(
                 sex = c("male", "female"), age="Age"
               )
)

plot(survfit(cox, newdata = newdat), xscale = 12,
     conf.int = T,
     ylim = c(0.6, 1),
     col = c("red", "blue"),
     xlab = "Time",
     ylab = "Survival proportion",
     main = "")
legend("bottomleft",
       legend=c("Male", "Female"),
       lty = 1, 
       col = c("red", "blue"),
       text.col = c("red", "blue"))

# Adding an interaction
cox.int <- coxph(imr ~ sex * age, data = infants)
summary(cox.int)
drop1(cox.int, test = "Chisq")
stargazer(cox.int, type = "text")

plot(survfit(cox.int, newdata = newdat), xscale = 12,
      conf.int = T,
      ylim = c(0.6, 1),
      col = c("male", "female"),
      xlab = "Age",
      ylab = "Survival proportion",
      main = "")

any advice on this would be most helpful! I can't figure out what I am doing wrong in generating a graph from the final data.

Upvotes: 0

Views: 44

Answers (1)

Gnueghoidune
Gnueghoidune

Reputation: 1357

To plot the impact of sex on survival probability, rename the sex levels in your newdat according to the sex levels in infants. Age is fixed to their average value, repeated for the number of different sex levels.

Code

data(infants)

imr <- with(infants, Surv(enter, exit, event))

cox <- coxph(imr ~ sex + age, data = infants)

cox_fit <- survfit(cox)

newdat <- with(infants, 
               data.frame(
                 sex = c("girl", "boy"), age=rep(mean(age, na.rm = TRUE), 2)))

plot(survfit(cox, newdata = newdat), xscale = 12,
     conf.int = T,
     ylim = c(0.6, 1),
     col = c("red", "blue"),
     xlab = "Time",
     ylab = "Survival proportion",
     main = "")
legend("bottomleft",
       legend=c("Male", "Female"),
       lty = 1, 
       col = c("red", "blue"),
       text.col = c("red", "blue"))

Output

enter image description here

I recommend using ggsurvplot instead, which allows for more customization:

ggsurvplot(survfit(cox, newdata = newdat), data = infants,
       legend.labs = c("Male", "Female"),
       legend.title = "Sex")

enter image description here

Upvotes: 1

Related Questions