diff --git a/R/mc5.R b/R/mc5.R index 2d3b627..1897bf4 100644 --- a/R/mc5.R +++ b/R/mc5.R @@ -122,245 +122,263 @@ mc5 <- function(ae, wr = FALSE) { invisible(rapply(exprs, eval, envir = fenv)) } - # if we're using v3 schema we want to tcplfit2 - dat <- tcplHit2(dat, coff = cutoff) - } else { - # Legacy fitting - - - - ## Apply the model type - dat[ , model_type := model_type] - - ## Determine winning model - dat[ , maic := pmin(cnst_aic, hill_aic, gnls_aic, na.rm = TRUE)] - # Order matters here, because in the case of a tie the simpler model will - # overwrite the more complex model as the winner. - dat[gnls_aic == maic, modl := "gnls"] - dat[hill_aic == maic, modl := "hill"] - dat[cnst_aic == maic, modl := "cnst"] - - ## Make the hitcall - dat[ , hitc := FALSE] - dat[modl == "hill" & hill_tp >= coff & max_med >= coff, hitc := TRUE] - dat[modl == "gnls" & gnls_tp >= coff & max_med >= coff, hitc := TRUE] - - ###--------------------- Bin the Dose-Response Sets -----------------------### - - ## Calculate AC05 and AC95 - dat[ , hill_95 := tcplHillACXX(95, hill_tp, hill_ga, hill_gw)] - dat[ , gnls_95 := tcplHillACXX(95, gnls_tp, gnls_ga, gnls_gw)] - - ## Add a few helper columns - dat[ , coff_upper := 1.2 * coff] - dat[ , coff_lower := 0.8 * coff] - dat[ , rgbl := !is.na(hill)] - dat[ , hill_gu := hill_tp > coff_upper] - dat[ , hill_lu := hill_tp <= coff_upper] - dat[ , hill_gl := hill_tp >= coff_lower] - dat[ , hill_ll := hill_tp < coff_lower] - dat[ , hill_gc := hill_tp >= coff] - dat[ , hill_lc := hill_tp < coff] - dat[ , hill_al := hill_ga <= logc_min] - dat[ , hill_ai := hill_ga > logc_min & hill_95 < logc_max] - dat[ , hill_au := hill_ga > logc_min & hill_95 >= logc_max] - dat[ , gnls_gu := gnls_tp > coff_upper] - dat[ , gnls_lu := gnls_tp <= coff_upper] - dat[ , gnls_gl := gnls_tp >= coff_lower] - dat[ , gnls_ll := gnls_tp < coff_lower] - dat[ , gnls_gc := gnls_tp >= coff] - dat[ , gnls_lc := gnls_tp < coff] - dat[ , gnls_al := gnls_ga <= logc_min] - dat[ , gnls_ai := gnls_ga > logc_min & gnls_95 < logc_max] - dat[ , gnls_au := gnls_ga > logc_min & gnls_95 >= logc_max] - dat[ , cnst_win := modl == "cnst"] - dat[ , hill_win := modl == "hill"] - dat[ , gnls_win := modl == "gnls"] - - ## Initiate fitc column - dat[ , fitc := NA_integer_] - - ## 02: CANNOT DETERMINE - dat[is.na(cnst), fitc := 2L] - ## 04: RESP < BLINE - dat[!rgbl & !is.na(cnst), fitc := 4L] - ## 07: NO TP >= 0.8(COFF) - dat[rgbl & cnst_win & (!hill_gl|is.na(hill_tp)) & (!gnls_gl|is.na(gnls_tp)), - fitc := 7L] - ## 09: NO TP >= COFF - dat[is.na(fitc) & rgbl & cnst_win & - (!hill_gc|is.na(hill_tp)) & (!gnls_gc|is.na(gnls_tp)), - fitc := 9L] - ## 10: ANY TP >= COFF - dat[is.na(fitc) & rgbl & cnst_win & (hill_gc | gnls_gc), fitc := 10L] - ## 54: GNLS DNC - dat[!hitc & hill_win & hill_ll & !gnls, fitc := 54] - ## 13: GNLS < 0.8(COFF) - dat[is.na(fitc) & !hitc & hill_win & hill_ll & gnls_ll, fitc := 13L] - ## 15: GNLS < COFF - dat[is.na(fitc) & !hitc & hill_win & hill_ll & gnls_lc, fitc := 15L] - ## 16: GNLS >= COFF - dat[is.na(fitc) & !hitc & hill_win & hill_ll & gnls_gc, fitc := 16L] - ## 55: GNLS DNC - dat[!hitc & hill_win & hill_gl & !gnls, fitc := 55] - ## 18: GNLS < 0.8(COFF) - dat[is.na(fitc) & !hitc & hill_win & hill_gl & gnls_ll, fitc := 18L] - ## 20: GNLS < COFF - dat[is.na(fitc) & !hitc & hill_win & hill_gl & gnls_lc, fitc := 20L] - ## 21: GNLS >= COFF - dat[is.na(fitc) & !hitc & hill_win & hill_gl & gnls_gc, fitc := 21L] - ## 52: HILL DNC - dat[!hitc & gnls_win & gnls_ll & !hill, fitc := 52] - ## 24: HILL < 0.8(COFF) - dat[is.na(fitc) & !hitc & gnls_win & gnls_ll & hill_ll, fitc := 24L] - ## 26: HILL < COFF - dat[is.na(fitc) & !hitc & gnls_win & gnls_ll & hill_lc, fitc := 26L] - ## 27: HILL >= COFF - dat[is.na(fitc) & !hitc & gnls_win & gnls_ll & hill_gc, fitc := 27L] - ## 53: HILL DNC - dat[!hitc & gnls_win & gnls_gl & !hill, fitc := 53] - ## 29: HILL < 0.8(COFF) - dat[is.na(fitc) & !hitc & gnls_win & gnls_gl & hill_ll, fitc := 29L] - ## 31: HILL < COFF - dat[is.na(fitc) & !hitc & gnls_win & gnls_gl & hill_lc, fitc := 31L] - ## 32: HILL >= COFF - dat[is.na(fitc) & !hitc & gnls_win & gnls_gl & hill_gc, fitc := 32L] - ## 36: GA <= LCONC - dat[hitc & hill_win & hill_lu & hill_al, fitc := 36L] - ## 37: LCONC < GA < HCONC - dat[hitc & hill_win & hill_lu & hill_ai, fitc := 37L] - ## 38: GA >= HCONC - dat[hitc & hill_win & hill_lu & hill_au, fitc := 38L] - ## 40: GA <= LCONC - dat[hitc & hill_win & hill_gu & hill_al, fitc := 40L] - ## 41: LCONC < GA < HCONC - dat[hitc & hill_win & hill_gu & hill_ai, fitc := 41L] - ## 42: GA >= HCONC - dat[hitc & hill_win & hill_gu & hill_au, fitc := 42L] - ## 45: GA <= LCONC - dat[hitc & gnls_win * gnls_lu & gnls_al, fitc := 45L] - ## 46: LCONC < GA < HCONC - dat[hitc & gnls_win * gnls_lu & gnls_ai, fitc := 46L] - ## 47: GA >= HCONC - dat[hitc & gnls_win * gnls_lu & gnls_au, fitc := 47L] - ## 49: GA <= LCONC - dat[hitc & gnls_win * gnls_gu & gnls_al, fitc := 49L] - ## 50: LCONC < GA < HCONC - dat[hitc & gnls_win * gnls_gu & gnls_ai, fitc := 50L] - ## 51: GA >= HCONC - dat[hitc & gnls_win * gnls_gu & gnls_au, fitc := 51L] - - ## Update hit-calls with the cannot determine call - dat[ , hitc := as.integer(hitc)] - dat[is.na(cnst), hitc := -1L] - - ## Add model fields - modl_pars <- c("modl_acb", - "modl_acc", - "modl_ac10", - "modl_er", - "modl_tp", - "modl_ga", - "modl_gw", - "modl_la", - "modl_lw", - "modl_rmse", - "modl_prob") - dat[ , modl_pars := NA_real_] - dat[modl == "cnst", modl_er := cnst_er] - dat[modl == "cnst", modl_rmse := cnst_rmse] - dat[modl == "cnst", modl_prob := cnst_prob] - dat[modl == "hill", modl_er := hill_er] - dat[modl == "hill", modl_tp := hill_tp] - dat[modl == "hill", modl_ga := hill_ga] - dat[modl == "hill", modl_gw := hill_gw] - dat[modl == "hill", modl_rmse := hill_rmse] - dat[modl == "hill", modl_prob := hill_prob] - dat[modl == "hill" & hill_tp >= 3 * bmad, - modl_acb := tcplHillConc(3 * bmad, hill_tp, hill_ga, hill_gw)] - dat[modl == "hill" & hill_tp >= coff, - modl_acc := tcplHillConc(coff, hill_tp, hill_ga, hill_gw)] - dat[modl == "hill", modl_ac10 := tcplHillACXX(10, hill_tp, hill_ga, hill_gw)] - dat[modl == "gnls", modl_er := gnls_er] - dat[modl == "gnls", modl_tp := gnls_tp] - dat[modl == "gnls", modl_ga := gnls_ga] - dat[modl == "gnls", modl_gw := gnls_gw] - dat[modl == "gnls", modl_la := gnls_la] - dat[modl == "gnls", modl_lw := gnls_lw] - dat[modl == "gnls", modl_rmse := gnls_rmse] - dat[modl == "gnls", modl_prob := gnls_prob] - dat[modl == "gnls" & gnls_tp >= 3 * bmad, - modl_acb := tcplHillConc(3 * bmad, gnls_tp, gnls_ga, gnls_gw)] - dat[modl == "gnls" & gnls_tp >= coff, - modl_acc := tcplHillConc(coff, gnls_tp, gnls_ga, gnls_gw)] - dat[modl == "gnls", modl_ac10 := tcplHillACXX(10, gnls_tp, gnls_ga, gnls_gw)] - - ## Add activity probability - dat[ , actp := 1 - cnst_prob] - - ## Complete the loec calculations - if (loec.mthd) { - # Complete the todo list to adjust for the loec method by calling loec.coff in mc5_mthds - - coff <- unique(dat$coff) # coff for aeid - calc_z <-function(resp) { + ## Complete the loec calculations + if (loec.mthd) { - # Original Z-score methodology - #if (length(resp) <= 1) {sdev=0.1} - #else {sdev=sd(resp)} - #mu = mean(resp) - #Z = (mu - coff)/sdev - - # New methodology where all resp > coff - above.coff <- all(resp>coff) # All responses must be greater than coff - if (above.coff == T){ - Z = 1 - } else { - Z = 0 + all_resp_gt_conc <-function(resp) { + # all resp > coff + return(as.integer(all(abs(resp) > cutoff))) # All responses must be greater than coff } + mc3 <- tcplLoadData(3L, fld='aeid', val=ae, type='mc') + + mc3[, loec_coff := lapply(.SD, all_resp_gt_conc), by=.(spid, conc), .SDcols = c("resp")] + suppressWarnings(mc3[, loec := min(conc[loec_coff == 1]), by = spid]) # Define the loec for each SPID + mc3[is.infinite(loec), loec := NA] #convert Inf to NA + mc3[, loec_hitc := max(loec_coff), by = spid] # is there a loec? used for hitc + mc3 <- mc3[dat, mult='first', on='spid', nomatch=0L] + + dat <- dat[mc3[,c("spid","loec", "loec_hitc")],on = "spid"] - return(Z) + } else { + # if we're using v3 schema and not loec method we want to tcplfit2 + dat <- tcplHit2(dat, coff = cutoff) } + } else { + # Legacy fitting + + ## Apply the model type + dat[ , model_type := model_type] + ## Determine winning model + dat[ , maic := pmin(cnst_aic, hill_aic, gnls_aic, na.rm = TRUE)] + # Order matters here, because in the case of a tie the simpler model will + # overwrite the more complex model as the winner. + dat[gnls_aic == maic, modl := "gnls"] + dat[hill_aic == maic, modl := "hill"] + dat[cnst_aic == maic, modl := "cnst"] - tmp.mc3 <- tcplLoadData(3L, fld='aeid', val=ae, type='mc') + ## Make the hitcall + dat[ , hitc := FALSE] + dat[modl == "hill" & hill_tp >= coff & max_med >= coff, hitc := TRUE] + dat[modl == "gnls" & gnls_tp >= coff & max_med >= coff, hitc := TRUE] + ###--------------------- Bin the Dose-Response Sets -----------------------### - tmp.mc3[, Z:=lapply(.SD, calc_z), by=.(spid, logc), .SDcols = c("resp")] - tmp.mc3[Z >= 1, loec_coff :=1] - tmp.mc3[Z < 1, loec_coff :=0] - suppressWarnings(tmp.mc3[, loec := min(logc[loec_coff == 1]), by = spid]) # Define the loec for each SPID - tmp.mc3 <- tmp.mc3[dat, mult='first', on='spid', nomatch=0L] - tmp.mc3[is.infinite(loec), loec_coff :=0] - tmp.mc3[is.finite(loec), loec_coff :=1] - is.na(tmp.mc3$loec) <- !is.finite(tmp.mc3$loec) # change + ## Calculate AC05 and AC95 + dat[ , hill_95 := tcplHillACXX(95, hill_tp, hill_ga, hill_gw)] + dat[ , gnls_95 := tcplHillACXX(95, gnls_tp, gnls_ga, gnls_gw)] + ## Add a few helper columns + dat[ , coff_upper := 1.2 * coff] + dat[ , coff_lower := 0.8 * coff] + dat[ , rgbl := !is.na(hill)] + dat[ , hill_gu := hill_tp > coff_upper] + dat[ , hill_lu := hill_tp <= coff_upper] + dat[ , hill_gl := hill_tp >= coff_lower] + dat[ , hill_ll := hill_tp < coff_lower] + dat[ , hill_gc := hill_tp >= coff] + dat[ , hill_lc := hill_tp < coff] + dat[ , hill_al := hill_ga <= logc_min] + dat[ , hill_ai := hill_ga > logc_min & hill_95 < logc_max] + dat[ , hill_au := hill_ga > logc_min & hill_95 >= logc_max] + dat[ , gnls_gu := gnls_tp > coff_upper] + dat[ , gnls_lu := gnls_tp <= coff_upper] + dat[ , gnls_gl := gnls_tp >= coff_lower] + dat[ , gnls_ll := gnls_tp < coff_lower] + dat[ , gnls_gc := gnls_tp >= coff] + dat[ , gnls_lc := gnls_tp < coff] + dat[ , gnls_al := gnls_ga <= logc_min] + dat[ , gnls_ai := gnls_ga > logc_min & gnls_95 < logc_max] + dat[ , gnls_au := gnls_ga > logc_min & gnls_95 >= logc_max] + dat[ , cnst_win := modl == "cnst"] + dat[ , hill_win := modl == "hill"] + dat[ , gnls_win := modl == "gnls"] - dat <- dat[tmp.mc3[,c("spid","loec","loec_coff")],on = "spid"] - dat[(!cnst_win), modl_acc := loec] - dat[(!cnst_win), modl_acb := loec] - dat[(!cnst_win), modl_ga := loec] - dat[(!cnst_win), fitc := 100L] - dat[(!cnst_win), model_type := 1] - dat[(!cnst_win), hitc := loec_coff] - dat <- dat[,-c("loec","loec_coff")] + ## Initiate fitc column + dat[ , fitc := NA_integer_] + ## 02: CANNOT DETERMINE + dat[is.na(cnst), fitc := 2L] + ## 04: RESP < BLINE + dat[!rgbl & !is.na(cnst), fitc := 4L] + ## 07: NO TP >= 0.8(COFF) + dat[rgbl & cnst_win & (!hill_gl|is.na(hill_tp)) & (!gnls_gl|is.na(gnls_tp)), + fitc := 7L] + ## 09: NO TP >= COFF + dat[is.na(fitc) & rgbl & cnst_win & + (!hill_gc|is.na(hill_tp)) & (!gnls_gc|is.na(gnls_tp)), + fitc := 9L] + ## 10: ANY TP >= COFF + dat[is.na(fitc) & rgbl & cnst_win & (hill_gc | gnls_gc), fitc := 10L] + ## 54: GNLS DNC + dat[!hitc & hill_win & hill_ll & !gnls, fitc := 54] + ## 13: GNLS < 0.8(COFF) + dat[is.na(fitc) & !hitc & hill_win & hill_ll & gnls_ll, fitc := 13L] + ## 15: GNLS < COFF + dat[is.na(fitc) & !hitc & hill_win & hill_ll & gnls_lc, fitc := 15L] + ## 16: GNLS >= COFF + dat[is.na(fitc) & !hitc & hill_win & hill_ll & gnls_gc, fitc := 16L] + ## 55: GNLS DNC + dat[!hitc & hill_win & hill_gl & !gnls, fitc := 55] + ## 18: GNLS < 0.8(COFF) + dat[is.na(fitc) & !hitc & hill_win & hill_gl & gnls_ll, fitc := 18L] + ## 20: GNLS < COFF + dat[is.na(fitc) & !hitc & hill_win & hill_gl & gnls_lc, fitc := 20L] + ## 21: GNLS >= COFF + dat[is.na(fitc) & !hitc & hill_win & hill_gl & gnls_gc, fitc := 21L] + ## 52: HILL DNC + dat[!hitc & gnls_win & gnls_ll & !hill, fitc := 52] + ## 24: HILL < 0.8(COFF) + dat[is.na(fitc) & !hitc & gnls_win & gnls_ll & hill_ll, fitc := 24L] + ## 26: HILL < COFF + dat[is.na(fitc) & !hitc & gnls_win & gnls_ll & hill_lc, fitc := 26L] + ## 27: HILL >= COFF + dat[is.na(fitc) & !hitc & gnls_win & gnls_ll & hill_gc, fitc := 27L] + ## 53: HILL DNC + dat[!hitc & gnls_win & gnls_gl & !hill, fitc := 53] + ## 29: HILL < 0.8(COFF) + dat[is.na(fitc) & !hitc & gnls_win & gnls_gl & hill_ll, fitc := 29L] + ## 31: HILL < COFF + dat[is.na(fitc) & !hitc & gnls_win & gnls_gl & hill_lc, fitc := 31L] + ## 32: HILL >= COFF + dat[is.na(fitc) & !hitc & gnls_win & gnls_gl & hill_gc, fitc := 32L] + ## 36: GA <= LCONC + dat[hitc & hill_win & hill_lu & hill_al, fitc := 36L] + ## 37: LCONC < GA < HCONC + dat[hitc & hill_win & hill_lu & hill_ai, fitc := 37L] + ## 38: GA >= HCONC + dat[hitc & hill_win & hill_lu & hill_au, fitc := 38L] + ## 40: GA <= LCONC + dat[hitc & hill_win & hill_gu & hill_al, fitc := 40L] + ## 41: LCONC < GA < HCONC + dat[hitc & hill_win & hill_gu & hill_ai, fitc := 41L] + ## 42: GA >= HCONC + dat[hitc & hill_win & hill_gu & hill_au, fitc := 42L] + ## 45: GA <= LCONC + dat[hitc & gnls_win * gnls_lu & gnls_al, fitc := 45L] + ## 46: LCONC < GA < HCONC + dat[hitc & gnls_win * gnls_lu & gnls_ai, fitc := 46L] + ## 47: GA >= HCONC + dat[hitc & gnls_win * gnls_lu & gnls_au, fitc := 47L] + ## 49: GA <= LCONC + dat[hitc & gnls_win * gnls_gu & gnls_al, fitc := 49L] + ## 50: LCONC < GA < HCONC + dat[hitc & gnls_win * gnls_gu & gnls_ai, fitc := 50L] + ## 51: GA >= HCONC + dat[hitc & gnls_win * gnls_gu & gnls_au, fitc := 51L] + ## Update hit-calls with the cannot determine call + dat[ , hitc := as.integer(hitc)] + dat[is.na(cnst), hitc := -1L] + ## Add model fields + modl_pars <- c("modl_acb", + "modl_acc", + "modl_ac10", + "modl_er", + "modl_tp", + "modl_ga", + "modl_gw", + "modl_la", + "modl_lw", + "modl_rmse", + "modl_prob") + dat[ , modl_pars := NA_real_] + dat[modl == "cnst", modl_er := cnst_er] + dat[modl == "cnst", modl_rmse := cnst_rmse] + dat[modl == "cnst", modl_prob := cnst_prob] + dat[modl == "hill", modl_er := hill_er] + dat[modl == "hill", modl_tp := hill_tp] + dat[modl == "hill", modl_ga := hill_ga] + dat[modl == "hill", modl_gw := hill_gw] + dat[modl == "hill", modl_rmse := hill_rmse] + dat[modl == "hill", modl_prob := hill_prob] + dat[modl == "hill" & hill_tp >= 3 * bmad, + modl_acb := tcplHillConc(3 * bmad, hill_tp, hill_ga, hill_gw)] + dat[modl == "hill" & hill_tp >= coff, + modl_acc := tcplHillConc(coff, hill_tp, hill_ga, hill_gw)] + dat[modl == "hill", modl_ac10 := tcplHillACXX(10, hill_tp, hill_ga, hill_gw)] + dat[modl == "gnls", modl_er := gnls_er] + dat[modl == "gnls", modl_tp := gnls_tp] + dat[modl == "gnls", modl_ga := gnls_ga] + dat[modl == "gnls", modl_gw := gnls_gw] + dat[modl == "gnls", modl_la := gnls_la] + dat[modl == "gnls", modl_lw := gnls_lw] + dat[modl == "gnls", modl_rmse := gnls_rmse] + dat[modl == "gnls", modl_prob := gnls_prob] + dat[modl == "gnls" & gnls_tp >= 3 * bmad, + modl_acb := tcplHillConc(3 * bmad, gnls_tp, gnls_ga, gnls_gw)] + dat[modl == "gnls" & gnls_tp >= coff, + modl_acc := tcplHillConc(coff, gnls_tp, gnls_ga, gnls_gw)] + dat[modl == "gnls", modl_ac10 := tcplHillACXX(10, gnls_tp, gnls_ga, gnls_gw)] - # ms <- tcplMthdLoad(lvl = 5L, id = ae, type = "mc") - # ms <- ms[mthd_id == 13] - # - # exprs <- lapply(mthd_funcs[ms$mthd], do.call, args = list(dat)) - # fenv <- environment() - # invisible(rapply(exprs, eval, envir = fenv)) + ## Add activity probability + dat[ , actp := 1 - cnst_prob] - } - - outcols <- c("m4id", "aeid", "modl", "hitc", "fitc", - "coff", "actp", "model_type", modl_pars) # Added model_type here - dat <- dat[ , .SD, .SDcols = outcols] + ## Complete the loec calculations + if (loec.mthd) { + # Complete the todo list to adjust for the loec method by calling loec.coff in mc5_mthds + + coff <- unique(dat$coff) # coff for aeid + calc_z <-function(resp) { + + # Original Z-score methodology + #if (length(resp) <= 1) {sdev=0.1} + #else {sdev=sd(resp)} + #mu = mean(resp) + #Z = (mu - coff)/sdev + + # New methodology where all resp > coff + above.coff <- all(resp>coff) # All responses must be greater than coff + if (above.coff == T){ + Z = 1 + } else { + Z = 0 + } + + + return(Z) + } + + + tmp.mc3 <- tcplLoadData(3L, fld='aeid', val=ae, type='mc') + + + tmp.mc3[, Z:=lapply(.SD, calc_z), by=.(spid, logc), .SDcols = c("resp")] + tmp.mc3[Z >= 1, loec_coff :=1] + tmp.mc3[Z < 1, loec_coff :=0] + suppressWarnings(tmp.mc3[, loec := min(logc[loec_coff == 1]), by = spid]) # Define the loec for each SPID + tmp.mc3 <- tmp.mc3[dat, mult='first', on='spid', nomatch=0L] + tmp.mc3[is.infinite(loec), loec_coff :=0] + tmp.mc3[is.finite(loec), loec_coff :=1] + is.na(tmp.mc3$loec) <- !is.finite(tmp.mc3$loec) # change + + + dat <- dat[tmp.mc3[,c("spid","loec","loec_coff")],on = "spid"] + dat[(!cnst_win), modl_acc := loec] + dat[(!cnst_win), modl_acb := loec] + dat[(!cnst_win), modl_ga := loec] + dat[(!cnst_win), fitc := 100L] + dat[(!cnst_win), model_type := 1] + dat[(!cnst_win), hitc := loec_coff] + dat <- dat[,-c("loec","loec_coff")] + + + + + # ms <- tcplMthdLoad(lvl = 5L, id = ae, type = "mc") + # ms <- ms[mthd_id == 13] + # + # exprs <- lapply(mthd_funcs[ms$mthd], do.call, args = list(dat)) + # fenv <- environment() + # invisible(rapply(exprs, eval, envir = fenv)) + + } + + outcols <- c("m4id", "aeid", "modl", "hitc", "fitc", + "coff", "actp", "model_type", modl_pars) # Added model_type here + dat <- dat[ , .SD, .SDcols = outcols] } @@ -371,6 +389,13 @@ mc5 <- function(ae, wr = FALSE) { invisible(rapply(exprs, eval, envir = fenv)) } + # apply loec.coff + if (loec.mthd) { + exprs <- lapply(mthd_funcs[c("loec.coff")], do.call, args = list()) + fenv <- environment() + invisible(rapply(exprs, eval, envir = fenv)) + } + ttime <- round(difftime(Sys.time(), stime, units = "sec"), 2) ttime <- paste(unclass(ttime), units(ttime)) diff --git a/R/mc5_mthds.R b/R/mc5_mthds.R index ff27271..4ab15c6 100644 --- a/R/mc5_mthds.R +++ b/R/mc5_mthds.R @@ -106,8 +106,10 @@ #' \item{maxmed20pct}{Add a cutoff value of 20 percent of the maximum of all endpoint maximal #' average response values (max_med).} #' \item{coff_2.32}{Add a cutoff value of 2.32.} -#' \item{loec.coff}{Method not yet updated for tcpl implementation. Identify the lowest observed -#' effective concentration (loec) compared to baseline.} +#' \item{loec.coff}{Identify the lowest observed effective concentration (loec) where the values +#' of all responses are outside the cutoff band (i.e. abs(resp) > cutoff). If loec exists, assume +#' hit call = 1, fitc = 100, model_type = 1. Winning model is not selected based on curve fits +#' and therefore additional potency estimates are not derived.} #' \item{ow_bidirectional_loss}{Multiply winning model hitcall (hitc) by -1 for models fit in the #' positive analysis direction. Typically used for endpoints where only negative responses are #' biologically relevant.} @@ -350,6 +352,15 @@ mc5_mthds <- function(ae) { e1 <- bquote(coff <- c(coff, 40)) list(e1) + }, + + loec.coff = function() { + + e1 <- bquote(dat[, c("modl", "fitc", "model_type", "hitc") := list("loec", 100L, 1, loec_hitc)]) + e2 <- bquote(dat <- dat |> melt(measure.vars = c("loec"), variable.name = "hit_param", value.name = "hit_val")) + e3 <- bquote(dat <- dat[,c("m4id", "aeid", "modl", "hitc", "fitc", "coff", "model_type", "hit_param", "hit_val")]) + list(e1, e2, e3) + } ) diff --git a/man/MC5_Methods.Rd b/man/MC5_Methods.Rd index bd6dffb..9c47a0b 100644 --- a/man/MC5_Methods.Rd +++ b/man/MC5_Methods.Rd @@ -111,8 +111,10 @@ Log Base 2 \item{maxmed20pct}{Add a cutoff value of 20 percent of the maximum of all endpoint maximal average response values (max_med).} \item{coff_2.32}{Add a cutoff value of 2.32.} - \item{loec.coff}{Method not yet updated for tcpl implementation. Identify the lowest observed - effective concentration (loec) compared to baseline.} + \item{loec.coff}{Identify the lowest observed effective concentration (loec) where the values + of all responses are outside the cutoff band (i.e. abs(resp) > cutoff). If loec exists, assume + hit call = 1, fitc = 100, model_type = 1. Winning model is not selected based on curve fits + and therefore additional potency estimates are not derived.} \item{ow_bidirectional_loss}{Multiply winning model hitcall (hitc) by -1 for models fit in the positive analysis direction. Typically used for endpoints where only negative responses are biologically relevant.}