Event History Analysis Ass2

load library

gc()

##           used  (Mb) gc trigger  (Mb) limit (Mb) max used  (Mb)
## Ncells 2904526 155.2    5601877 299.2         NA  3999272 213.6
## Vcells 4873666  37.2   10146329  77.5      16384  6981735  53.3

library(ipumsr)

## 
## Attaching package: 'ipumsr'

## The following objects are masked from 'package:sjlabelled':
## 
##     as_factor, zap_labels

gc()

##           used  (Mb) gc trigger  (Mb) limit (Mb) max used  (Mb)
## Ncells 2908488 155.4    5601877 299.2         NA  3999272 213.6
## Vcells 4877932  37.3   10146329  77.5      16384  6981735  53.3

setwd("~/Downloads")

Read the IPUMS NHIS data

gc()

##           used  (Mb) gc trigger  (Mb) limit (Mb) max used  (Mb)
## Ncells 2908484 155.4    5601877 299.2         NA  3999272 213.6
## Vcells 4878089  37.3   10146329  77.5      16384  6981735  53.3

ddi <- read_ipums_ddi("nhis_00008.xml")
data <- read_ipums_micro(ddi)

## Use of data from IPUMS NHIS is subject to conditions including that users
## should cite the data appropriately. Use command `ipums_conditions()` for more
## details.

data<- haven::zap_labels(data)

Filter data to only mortality eligible and young adult age 18-35

data <- data %>%
  filter(MORTELIG == 1)

data <- data %>%
filter(AGE ==18)

data <- data%>%
filter(complete.cases(.))

1. Define your event variable

Answer: Young Adult mortality in the United State(Age 18-35),

This analysis uses NHIS data to see the mortality event among young adult age 18-35 in the United state. Due to data limitation the data used 2001-2018 variable as those who are eligible in 2001 will be 18years and their maximum age in 2018 will be 35years 2001-2018

2. Define a duration or time variable.

I would be considering time of young adult death

3. Define a censoring indicator

1 means dead, 0 means alive

4. Estimate the survival function for your outcome and plot it : Age at Death

Cross tabulation of mortality status and year

table(data$MORTSTAT, data$YEAR)

##    
##     2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015
##   1   27   26   23   21   23   16    6    8   22   15    9    7   12    3    1
##   2 1227 1087 1111 1166 1132  960 1019 1048 1259 1234 1285 1386 1260 1411  321
##    
##     2016 2017 2018
##   1    2    0    0
##   2  396  306  204

data <- data %>%
  mutate(death_age = ifelse( MORTSTAT ==1, 
                             MORTDODY - (YEAR - AGE), 
                             2018 - (YEAR - AGE)), 
         d.event = ifelse(MORTSTAT == 1, 1, 0))


library(survival)

age_fit <- survfit(Surv(death_age, d.event) ~ 1, 
                   data = data)

library(ggsurvfit)

age_fit %>%
  ggsurvfit() +
  add_confidence_interval(type = "ribbon") +
  add_quantile() 

5a. Kaplan-Meier survival analysis of the outcome( time to death)

Time to death

data <- data  %>%
  mutate(death_time = ifelse( MORTSTAT ==1, 
                             MORTDODY - YEAR , 
                             2018 - YEAR ), 
         d5.event = ifelse(MORTSTAT == 1 & death_time <= 5, 1, 0))

library(survival)

time_fit <- survfit(Surv(death_time, d5.event) ~ 1, 
                  conf.type = "log",
                   data = data )


time_fit %>%
  ggsurvfit()+
  xlim(-1, 6) +
  add_censor_mark() +
  add_confidence_interval()+
  add_quantile(y_value = .975) 

## Warning: Removed 12 row(s) containing missing values (geom_path).

5 b. Define a grouping variable, this can be dichotomous or categorical

Gender is my grouping variable

5c Do you have a research hypothesis about the survival patterns for the levels of the categorical variable? State it.

My Assumption is that males will have higher mortality compared to male young adult(Age18-35)

#d.5d Comparison of Kaplan-Meier survival across grouping variables in your data. Interpret your results.

# Testing the hypothesis to see if there is a significant difference in the Pharmacy survival chances between the two groups 

survdiff(Surv(death_time, d.event) ~ SEX, 
                   data = data)

## Call:
## survdiff(formula = Surv(death_time, d.event) ~ SEX, data = data)
## 
##          N Observed Expected (O-E)^2/E (O-E)^2/V
## SEX=1 9332      175      114      32.1        67
## SEX=2 8701       46      107      34.5        67
## 
##  Chisq= 67  on 1 degrees of freedom, p= 3e-16

The result shows that there is a significant difference between males and females, as male was higher as hypothesized.

5e. Plot the survival function for the analysis for each level of the group variable

# Estimating the survival function for the groups
group_fit <- survfit(Surv(death_time, d.event) ~ SEX, 
                   data = data)

#summary(group_fit)
# Plotting the Survival functions for the two groups
group_fit%>%
  ggsurvfit() +
   add_confidence_interval() +
   add_risktable() +
  labs(title = "Young adult mortality Survival Plot by Gender",
       caption = "Source: NHIS 2001-2018 \n Calculations by Joseph Jaiyeola", x= "time",
      y = "Survival Probability")+  theme_minimal() +
  theme(legend.position = "bottom",axis.line = element_line(linetype = "solid"),
    axis.ticks = element_line(linetype = "blank"),
    panel.grid.major = element_line(colour = "gray80"),
    panel.grid.minor = element_line(colour = "gray94",
        linetype = "blank"), axis.title = element_text(family = "serif",
        face = "bold"), axis.text = element_text(face = "bold"),
    plot.title = element_text(family = "serif",
        size = 14, face = "bold") , legend.text = element_text(face = "bold",
        family = "serif"), legend.title = element_text(face = "bold",
        family = "serif"), panel.background = element_rect(fill = "gray97"),
    legend.key = element_rect(fill = NA,
        linetype = "solid"), legend.background = element_rect(linetype = "solid")) +
  guides(color = guide_legend(ncol = 1)) +
    coord_cartesian(xlim = c(0, 30)) +
  scale_y_continuous(
    limits = c(.95, 1),
    labels = scales::percent, 
    expand = c(0.01, 0))+
  scale_color_manual(values = c('#54738E', '#82AC7C')) +
  scale_fill_manual(values = c('#54738E', '#82AC7C')) 

## Warning: Removed 1 row(s) containing missing values (geom_path).