search

appliedepi / epiRhandbook_eng

The repository for the English version of the Epidemiologist R Handbook
Other
97 stars 56 forks source link

3_External_review #282

Closed jarvisc1 closed 3 weeks ago

arranhamlet commented 1 month ago

@nsbatra @aspina7 

dsr <- function(data, event, fu, subgroup, ..., refdata, mp, method = "normal", sig = 0.95, decimals ) {    

     subgroup <- enquo(subgroup)
     event <- enquo(event)
     fu <- enquo(fu)
     stdvars <- quos(...)

     #Compute crude and standardized rates and variances
     all_data_st = data %>%
          left_join(refdata) %>%
          group_by(!!subgroup) %>%
          mutate(n = sum(!!event, na.rm = T),
                 d = sum(!!fu),
                 cr_rate = n/d,
                 cr_var = n/d^2,
                 wts = pop/sum(pop, na.rm = T),
                 st_rate = sum(wts*(!!event/!!fu), na.rm = T),
                 st_var = sum(as.numeric((wts^2)*(!!event/(!!fu )^2)), na.rm = T)) %>%
          distinct(!!subgroup, .keep_all=TRUE) %>%
          select(!!subgroup, n, d, cr_rate, cr_var, st_rate, st_var)

     #Compute Confidence Intervals (CI) according to method. The default is 'normal'
     if(method=="gamma") {

          tmp1 =   all_data_st %>%
               mutate(
                    c_rate = mp*cr_rate,
                    c_lower = mp*qgamma((1-sig)/2, shape = cr_rate^2/(cr_var))/(cr_rate/cr_var),
                    c_upper = mp*qgamma(1-((1-sig)/2), shape = 1+cr_rate^2/(cr_var))/(cr_rate/cr_var),
                    s_rate = mp*st_rate,
                    s_lower = mp*qgamma((1-sig)/2, shape = st_rate^2/st_var)/(st_rate/st_var),
                    s_upper = mp*qgamma(1-((1-sig)/2), shape = 1+(st_rate^2/st_var))/(st_rate/st_var)) %>%
               select(!!subgroup, n, d, c_rate, c_lower, c_upper, s_rate, s_lower, s_upper)

     } else if (method == "normal") {

          tmp1 =   all_data_st %>%
               mutate(
                    c_rate = mp*cr_rate,
                    c_lower = mp*(cr_rate+qnorm((1-sig)/2)*sqrt(cr_var)),
                    c_upper = mp*(cr_rate-qnorm((1-sig)/2)*sqrt(cr_var)),
                    s_rate = mp*st_rate,
                    s_lower = mp*(st_rate+qnorm((1-sig)/2)*sqrt(st_var)),
                    s_upper = mp*(st_rate-qnorm((1-sig)/2)*sqrt(st_var))) %>%
               select(!!subgroup, n, d, c_rate, c_lower, c_upper, s_rate, s_lower, s_upper)

     } else if (method=="lognormal") {

          tmp1 =   all_data_st %>%
               mutate(
                    c_rate = mp*cr_rate,
                    c_lower = mp*exp((log(cr_rate)+qnorm((1-sig)/2)*sqrt(cr_var)/(cr_rate))),
                    c_upper = mp*exp((log(cr_rate)-qnorm((1-sig)/2)*sqrt(cr_var)/(cr_rate))),
                    s_rate = mp*st_rate,
                    s_lower = mp*exp((log(st_rate)+qnorm((1-sig)/2)*sqrt(st_var)/(st_rate))),
                    s_upper = mp*exp((log(st_rate)-qnorm((1-sig)/2)*sqrt(st_var)/(st_rate)))) %>%
               select(!!subgroup, n, d, c_rate, c_lower, c_upper, s_rate, s_lower, s_upper)

     }

     #Clean up and output

     tmp1$c_rate  <- round(tmp1$c_rate,  digits=decimals)
     tmp1$c_lower <- round(tmp1$c_lower, digits=decimals)
     tmp1$c_upper <- round(tmp1$c_upper, digits=decimals)
     tmp1$s_rate  <- round(tmp1$s_rate, digits=decimals)
     tmp1$s_lower <- round(tmp1$s_lower, digits=decimals)
     tmp1$s_upper <- round(tmp1$s_upper, digits=decimals)

     colnames(tmp1) <-  c('Subgroup', 'Numerator','Denominator',
                          paste('Crude Rate (per ', mp,')', sep = ''),
                          paste(sig*100,'% LCL (Crude)', sep = ''),
                          paste(sig*100,'% UCL (Crude)', sep = ''),
                          paste('Std Rate (per ', mp,')', sep = ''),
                          paste(sig*100,'% LCL (Std)', sep = ''),
                          paste(sig*100,'% UCL (Std)', sep = ''))

     tmp1 <- as.data.frame(tmp1)

}
aspina7 commented 1 month ago

There's a few parts to this. I think it was originally written by someone at UKHSA (back then PHE). Someone (maybe Liza Coyer?) at some point was going to reach out and see if they were using an alternative.

We could just recreate it and put it in {epikit} - but leads to a few questions. 1) probably for @jarvisc1 : are there more modern ways of doing this? 2) is this really done often enough in practice that we feel it merits us managing a package (personally I think probably not). And almost certainly not before the end of the year. 3) the existing code is.... clean enough.... But are there any tests at all? Assuming the package setup is probably very old school. If we were going to recreate it would probably need quite a bit of tearing apart and restarting.

aspina7 commented 1 month ago

Is this the new version?

https://cran.r-project.org/web/packages/PHEindicatormethods/vignettes/DSR-vignette.html

https://github.com/ukhsa-collaboration/PHEindicatormethods

nsbatra commented 1 month ago

As I recall there are two packages highlighted in the chapter, {dsr} and {pheindicatormethods}. There was some difference in functionality.

As Arran found, the first is deprecated.

If the PHE one becomes the only one in the Handbook, are there functionality gaps?

If gaps, we can try to revise {dsr} code in a new function (agree with Alex, not a priority) and maybe just post that code in Forum and link to it, before actually putting in Handbook?

We've gotten about 5 emails over the 3 years asking for help because {dsr} is gone. Not much, but not nothing.

My 2 cents

On Tue, Sep 10, 2024, 13:39 Alex Spina @.***> wrote:

Is this the new version?

https://cran.r-project.org/web/packages/PHEindicatormethods/vignettes/DSR-vignette.html

— Reply to this email directly, view it on GitHub https://github.com/appliedepi/epiRhandbook_eng/issues/282#issuecomment-2340425495, or unsubscribe https://github.com/notifications/unsubscribe-auth/AMN4O77QW4NA66U7IA72GITZV3LATAVCNFSM6AAAAABNSFRCH2VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDGNBQGQZDKNBZGU . You are receiving this because you were mentioned.Message ID: @.***>

nsbatra commented 1 month ago

Ah I see. Alex is saying PHE should cover us

On Tue, Sep 10, 2024, 14:12 Neale Batra @.***> wrote:

As I recall there are two packages highlighted in the chapter, {dsr} and {pheindicatormethods}. There was some difference in functionality.

As Arran found, the first is deprecated.

If the PHE one becomes the only one in the Handbook, are there functionality gaps?

If gaps, we can try to revise {dsr} code in a new function (agree with Alex, not a priority) and maybe just post that code in Forum and link to it, before actually putting in Handbook?

We've gotten about 5 emails over the 3 years asking for help because {dsr} is gone. Not much, but not nothing.

My 2 cents

On Tue, Sep 10, 2024, 13:39 Alex Spina @.***> wrote:

Is this the new version?

https://cran.r-project.org/web/packages/PHEindicatormethods/vignettes/DSR-vignette.html

— Reply to this email directly, view it on GitHub https://github.com/appliedepi/epiRhandbook_eng/issues/282#issuecomment-2340425495, or unsubscribe https://github.com/notifications/unsubscribe-auth/AMN4O77QW4NA66U7IA72GITZV3LATAVCNFSM6AAAAABNSFRCH2VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDGNBQGQZDKNBZGU . You are receiving this because you were mentioned.Message ID: @.***>

arranhamlet commented 1 month ago

Okay awesome, let's just go with the PHEindicatormethods DSR approach