Search code examples
rsurvival-analysisrate

Generate a crude incidence rate table (stratified by a factor variable) from a Lexis Model

I am using the 'Epi' package in R to model follow-up data from a study. I am having no issues with declaring the Lexis model or running Poisson and (combined with the survival package) Cox regressions.

As part of the initial data review I want to find a simple way to make a table of crude unadjusted incidence/event rates from data in a lexis model in R (pre-fitting any poisson/cox models).

I have found a coded approach which allows me to do this and to stratify by a variable as part of exploratory data analysis:

    #Generic Syntax Example
total <-cbind(tapply(lexis_model$lex.Xst,lexis_model$stratifying_var,sum),tapply(lexis_model$lex.dur,lexis_model$stratifying_var,sum))
    #Add up the number of events within the stratifying variable
    #Add up the amount of follow-up time within the stratifying the variable
    rates <- tapply(lexis_model$lex.Xst,lexis_model$stratifying_var,sum)/tapply(lexis_model$lex.dur,lexis_model$stratifying_var,sum)*10^3
    #Given rates per 1,000 person years
    ratetable <- (cbind(totals,rates))

#Specific Example based on the dataset
totals <-cbind(tapply(lexis_model$lex.Xst,lexis_model$grade,sum),tapply(lexis_model$lex.dur,lexis_model$grade,sum))
rates <- tapply(lexis_model$lex.Xst,lexis_model$grade,sum)/tapply(lexis_model$lex.dur,lexis_model$grade,sum)*10^3
ratetable <- (cbind(totals,rates))
ratetable

                 rates
1 90 20338.234 4.4251630
2 64  7265.065 8.8092811
#Shows number of events, years follow-up, number of events per 1000 years follow-up, stratified by the stratifying variable

Note this is crude unadjusted/absolute rates - not the output of a Poisson model. Whilst I appreciate that the code above does indeed produce the desired output (and is pretty straightforward) I wanted to see if people were aware of a command which can take a lexis dataset and output this. I've had a look at the available commands in the Epi and epitools package - may have missed somthing but could not see an obvious way to do this.

As this is a quite common thing to want to do I wondered if anyone was aware of a package/function that could do this by specifying the simply the lexis dataset and the stratification variable (or indeed a single function to the steps above in a single go).

Ideally the output would look something like the below (which is taken from STATA which I am trying to move away from in favour of R!):enter image description here

A copy of the first twenty rows or so of the actual data is here (the data has already been put in to a lexis model using Epi package so all relevant lexis variables are there): https://www.dropbox.com/s/yjyz1kzusysz941/rate_table_data.xlsx?dl=0

Solution

  • I would do this simply using the tidyverse R package as such:

    library(tidyverse)
lexis_model %>%
  group_by(grade) %>%
  summarise(sum_Xst = sum(lex.Xst), sum_dur = sum(lex.dur)) %>%
  mutate(rate = sum_Xst/sum_dur*10^3) -> rateable
rateable

#  A tibble: 2 x 4
#   grade sum_Xst   sum_dur     rate
#   <dbl>   <int>     <dbl>    <dbl>
# 1     1       2 375.24709 5.329821
# 2     2       0  92.44079 0.000000

    And you could wrap this into a function yourself:

    rateFunc <- function(data, strat_var)
{
  lexis_model %>%
    group_by_(strat_var) %>%
    summarise(sum_Xst = sum(lex.Xst), sum_dur = sum(lex.dur)) %>%
    mutate(rate = sum_Xst/sum_dur*10^3)
}

    which you would then call:

    rateFunc(lexis_model, "grade")

    This is useful because, using the combination of tidyverse summarise and mutate it is very easy to add more summary statistics to the table.

    EDIT: After clarification on the question, this can be done using the popEpi package using the rate command:

    popEpi::rate(lexis_model, obs = lex.Xst, pyrs = lex.dur, print = grade)

# Crude rates and 95% confidence intervals:
#    grade lex.Xst  lex.dur       rate     SE.rate     rate.lo   rate.hi
# 1:     1       2 375.2472 0.00532982 0.003768752 0.001332942 0.0213115
# 2:     2       0  92.4408 0.00000000 0.000000000 0.000000000       NaN