Issue using FengConv

[pwighton]

I'm getting an error in aif-fengconv. Perhaps it's an environment issue? I'm using the container from PR #9.

Running bloodstream on ds004230 using the following config file:

{
  "Subsets": {
    "sub": [""],
    "ses": [""],
    "rec": [""],
    "task": [""],
    "run": [""],
    "TracerName": [""],
    "ModeOfAdministration": [""],
    "InstitutionName": [""],
    "PharmaceuticalName": [""]
  },
  "Model": {
    "ParentFraction": {
      "Method": ["Fit Hierarchically: HGAM"],
      "set_ppf0": [true],
      "starttime": [0],
      "endtime": ["Inf"],
      "gam_k": ["6"],
      "hgam_formula": ["s(time, k=8) + s(time, sub, bs='fs', k=5) + s(time, pet, bs='fs', k=5)"]
    },
    "BPR": {
      "Method": ["Fit Hierarchically: HGAM"],
      "starttime": [0],
      "endtime": ["Inf"],
      "gam_k": [6],
      "hgam_formula": ["s(time, k=8) + s(time, sub, bs='fs', k=5) + s(time, pet, bs='fs', k=5)"]
    },
    "AIF": {
      "Method": ["Fit Individually: FengConv"],
      "starttime": [0],
      "endtime": ["Inf"],
      "expdecay_props": ["NA", "NA"],
      "inftime": ["NA"],
      "spline_kb": [""],
      "spline_ka_m": [""],
      "spline_ka_a": [""]
    },
    "WholeBlood": {
      "Method": ["Interpolation"],
      "dispcor": [false],
      "starttime": [0],
      "endtime": ["Inf"],
      "spline_kb": [""],
      "spline_ka_m": [""],
      "spline_ka_a": [""]
    }
  }
}

Generates the following error:

Quitting from lines 1054-1110 [aif-fengconv] (template.rmd)

Error in `mutate()`:
ℹ In argument: `aif_fit = map(aif, ~blmod_fengconv(time = .x$time, activity = .x$aif, Method = .x$Method))`.
Caused by error in `map()`:
ℹ In index: 1.
Caused by error in `if (startpars$t0 < 0) ...`:
! missing value where TRUE/FALSE needed
Traceback:

1. bloodstream(studypath = studypath, configpath = configfile)
2. rmarkdown::render(input = paste0(system.file(package = "bloodstream"), 
 .     "/rmd/template.rmd"), output_file = paste0(studypath, "/derivatives/bloodstream", 
 .     config_suffix, "/", "bloodstream_report_config", ifelse(stringr::str_length(config_suffix) > 
 .         0, yes = "-", no = ""), config_suffix, ".html"), params = list(configpath = configpath, 
 .     studypath = studypath), knit_root_dir = studypath)
3. knitr::knit(knit_input, knit_output, envir = envir, quiet = quiet)
4. process_file(text, output)
5. withCallingHandlers(withCallingHandlers(if (tangle) process_tangle(group) else process_group(group), 
 .     error = function(e) if (xfun::pkg_available("rlang", "1.0.0")) rlang::entrace(e)), 
 .     error = function(e) {
 .         setwd(wd)
 .         write_utf8(res, output %n% stdout())
 .         message("\nQuitting from lines ", paste(current_lines(i), 
 .             collapse = "-"), if (labels[i] != "") 
 .             sprintf(" [%s]", labels[i]), sprintf(" (%s)", knit_concord$get("infile")))
 .     })
6. withCallingHandlers(if (tangle) process_tangle(group) else process_group(group), 
 .     error = function(e) if (xfun::pkg_available("rlang", "1.0.0")) rlang::entrace(e))
7. process_group(group)
8. process_group.block(group)
9. call_block(x)
10. block_exec(params)
11. eng_r(options)
12. in_input_dir(evaluate(code, envir = env, new_device = FALSE, 
  .     keep_warning = if (is.numeric(options$warning)) TRUE else options$warning, 
  .     keep_message = if (is.numeric(options$message)) TRUE else options$message, 
  .     stop_on_error = if (is.numeric(options$error)) options$error else {
  .         if (options$error && options$include) 
  .             0L
  .         else 2L
  .     }, output_handler = knit_handlers(options$render, options)))
13. in_dir(input_dir(), expr)
14. evaluate(code, envir = env, new_device = FALSE, keep_warning = if (is.numeric(options$warning)) TRUE else options$warning, 
  .     keep_message = if (is.numeric(options$message)) TRUE else options$message, 
  .     stop_on_error = if (is.numeric(options$error)) options$error else {
  .         if (options$error && options$include) 
  .             0L
  .         else 2L
  .     }, output_handler = knit_handlers(options$render, options))
15. evaluate::evaluate(...)
16. evaluate_call(expr, parsed$src[[i]], envir = envir, enclos = enclos, 
  .     debug = debug, last = i == length(out), use_try = stop_on_error != 
  .         2L, keep_warning = keep_warning, keep_message = keep_message, 
  .     log_echo = log_echo, log_warning = log_warning, output_handler = output_handler, 
  .     include_timing = include_timing)
17. timing_fn(handle(ev <- withCallingHandlers(withVisible(eval_with_user_handlers(expr, 
  .     envir, enclos, user_handlers)), warning = wHandler, error = eHandler, 
  .     message = mHandler)))
18. handle(ev <- withCallingHandlers(withVisible(eval_with_user_handlers(expr, 
  .     envir, enclos, user_handlers)), warning = wHandler, error = eHandler, 
  .     message = mHandler))
19. withCallingHandlers(withVisible(eval_with_user_handlers(expr, 
  .     envir, enclos, user_handlers)), warning = wHandler, error = eHandler, 
  .     message = mHandler)
20. withVisible(eval_with_user_handlers(expr, envir, enclos, user_handlers))
21. eval_with_user_handlers(expr, envir, enclos, user_handlers)
22. eval(expr, envir, enclos)
23. eval(expr, envir, enclos)
24. aif_data %>% mutate(aif_fit = map(aif, ~blmod_fengconv(time = .x$time, 
  .     activity = .x$aif, Method = .x$Method)))   # at line 20-22 of file <text>
25. mutate(., aif_fit = map(aif, ~blmod_fengconv(time = .x$time, 
  .     activity = .x$aif, Method = .x$Method)))
26. mutate.data.frame(., aif_fit = map(aif, ~blmod_fengconv(time = .x$time, 
  .     activity = .x$aif, Method = .x$Method)))
27. mutate_cols(.data, dplyr_quosures(...), by)
28. withCallingHandlers(for (i in seq_along(dots)) {
  .     poke_error_context(dots, i, mask = mask)
  .     context_poke("column", old_current_column)
  .     new_columns <- mutate_col(dots[[i]], data, mask, new_columns)
  . }, error = dplyr_error_handler(dots = dots, mask = mask, bullets = mutate_bullets, 
  .     error_call = error_call, error_class = "dplyr:::mutate_error"), 
  .     warning = dplyr_warning_handler(state = warnings_state, mask = mask, 
  .         error_call = error_call))
29. mutate_col(dots[[i]], data, mask, new_columns)
30. mask$eval_all_mutate(quo)
31. eval()
32. map(aif, ~blmod_fengconv(time = .x$time, activity = .x$aif, Method = .x$Method))
33. map_("list", .x, .f, ..., .progress = .progress)
34. with_indexed_errors(i = i, names = names, error_call = .purrr_error_call, 
  .     call_with_cleanup(map_impl, environment(), .type, .progress, 
  .         n, names, i))
35. withCallingHandlers(expr, error = function(cnd) {
  .     if (i == 0L) {
  .     }
  .     else {
  .         message <- c(i = "In index: {i}.")
  .         if (!is.null(names) && !is.na(names[[i]]) && names[[i]] != 
  .             "") {
  .             name <- names[[i]]
  .             message <- c(message, i = "With name: {name}.")
  .         }
  .         else {
  .             name <- NULL
  .         }
  .         cli::cli_abort(message, location = i, name = name, parent = cnd, 
  .             call = error_call, class = "purrr_error_indexed")
  .     }
  . })
36. call_with_cleanup(map_impl, environment(), .type, .progress, 
  .     n, names, i)
37. .f(.x[[i]], ...)
38. blmod_fengconv(time = .x$time, activity = .x$aif, Method = .x$Method)
39. blmod_feng_startpars(time, activity, expdecay_props)
40. .handleSimpleError(function (cnd) 
  . {
  .     if (i == 0L) {
  .     }
  .     else {
  .         message <- c(i = "In index: {i}.")
  .         if (!is.null(names) && !is.na(names[[i]]) && names[[i]] != 
  .             "") {
  .             name <- names[[i]]
  .             message <- c(message, i = "With name: {name}.")
  .         }
  .         else {
  .             name <- NULL
  .         }
  .         cli::cli_abort(message, location = i, name = name, parent = cnd, 
  .             call = error_call, class = "purrr_error_indexed")
  .     }
  . }, "missing value where TRUE/FALSE needed", base::quote(if (startpars$t0 < 
  .     0) {
  .     bloodrise <- bloodrise[blood$activity >= 0.1 * startpars$peakval, 
  .         ]
  .     rise_lm <- lm(activity ~ time, data = bloodrise)
  .     rise_coef <- coef(rise_lm)
  .     startpars$t0 <- as.numeric(-1 * (rise_coef[1]/rise_coef[2]))
  . }))
41. h(simpleError(msg, call))
42. cli::cli_abort(message, location = i, name = name, parent = cnd, 
  .     call = error_call, class = "purrr_error_indexed")
43. rlang::abort(message, ..., call = call, use_cli_format = TRUE, 
  .     .frame = .frame)
44. signal_abort(cnd, .file)
45. signalCondition(cnd)
46. (function (cnd) 
  . {
  .     local_error_context(dots, i = frame[[i_sym]], mask = mask)
  .     if (inherits(cnd, "dplyr:::internal_error")) {
  .         parent <- error_cnd(message = bullets(cnd))
  .     }
  .     else {
  .         parent <- cnd
  .     }
  .     message <- c(cnd_bullet_header(action), i = if (has_active_group_context(mask)) cnd_bullet_cur_group_label())
  .     abort(message, class = error_class, parent = parent, call = error_call)
  . })(structure(list(message = c(i = "In index: 1."), trace = structure(list(
  .     call = list(aif_data %>% mutate(aif_fit = map(aif, ~blmod_fengconv(time = .x$time, 
  .         activity = .x$aif, Method = .x$Method))), mutate(., aif_fit = map(aif, 
  .         ~blmod_fengconv(time = .x$time, activity = .x$aif, Method = .x$Method))), 
  .         mutate.data.frame(., aif_fit = map(aif, ~blmod_fengconv(time = .x$time, 
  .             activity = .x$aif, Method = .x$Method))), mutate_cols(.data, 
  .             dplyr_quosures(...), by), withCallingHandlers(for (i in seq_along(dots)) {
  .             poke_error_context(dots, i, mask = mask)
  .             context_poke("column", old_current_column)
  .             new_columns <- mutate_col(dots[[i]], data, mask, 
  .                 new_columns)
  .         }, error = dplyr_error_handler(dots = dots, mask = mask, 
  .             bullets = mutate_bullets, error_call = error_call, 
  .             error_class = "dplyr:::mutate_error"), warning = dplyr_warning_handler(state = warnings_state, 
  .             mask = mask, error_call = error_call)), mutate_col(dots[[i]], 
  .             data, mask, new_columns), mask$eval_all_mutate(quo), 
  .         eval(), map(aif, ~blmod_fengconv(time = .x$time, activity = .x$aif, 
  .             Method = .x$Method)), map_("list", .x, .f, ..., .progress = .progress), 
  .         with_indexed_errors(i = i, names = names, error_call = .purrr_error_call, 
  .             call_with_cleanup(map_impl, environment(), .type, 
  .                 .progress, n, names, i)), withCallingHandlers(expr, 
  .             error = function(cnd) {
  .                 if (i == 0L) {
  .                 }
  .                 else {
  .                   message <- c(i = "In index: {i}.")
  .                   if (!is.null(names) && !is.na(names[[i]]) && 
  .                     names[[i]] != "") {
  .                     name <- names[[i]]
  .                     message <- c(message, i = "With name: {name}.")
  .                   }
  .                   else {
  .                     name <- NULL
  .                   }
  .                   cli::cli_abort(message, location = i, name = name, 
  .                     parent = cnd, call = error_call, class = "purrr_error_indexed")
  .                 }
  .             }), call_with_cleanup(map_impl, environment(), .type, 
  .             .progress, n, names, i), .f(.x[[i]], ...), blmod_fengconv(time = .x$time, 
  .             activity = .x$aif, Method = .x$Method), blmod_feng_startpars(time, 
  .             activity, expdecay_props), .handleSimpleError(`<fn>`, 
  .             "missing value where TRUE/FALSE needed", base::quote(if (startpars$t0 < 
  .                 0) {
  .                 bloodrise <- bloodrise[blood$activity >= 0.1 * 
  .                   startpars$peakval, ]
  .                 rise_lm <- lm(activity ~ time, data = bloodrise)
  .                 rise_coef <- coef(rise_lm)
  .                 startpars$t0 <- as.numeric(-1 * (rise_coef[1]/rise_coef[2]))
  .             })), h(simpleError(msg, call)), cli::cli_abort(message, 
  .             location = i, name = name, parent = cnd, call = error_call, 
  .             class = "purrr_error_indexed"), rlang::abort(message, 
  .             ..., call = call, use_cli_format = TRUE, .frame = .frame)), 
  .     parent = c(0L, 0L, 0L, 3L, 4L, 4L, 6L, 7L, 0L, 9L, 10L, 11L, 
  .     10L, 10L, 14L, 15L, 0L, 17L, 18L, 19L), visible = c(TRUE, 
  .     TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, 
  .     TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE), 
  .     namespace = c(NA, "dplyr", "dplyr", "dplyr", "base", "dplyr", 
  .     NA, "dplyr", "purrr", "purrr", "purrr", "base", "purrr", 
  .     "bloodstream", "kinfitr", "kinfitr", "base", "purrr", "cli", 
  .     "rlang"), scope = c(NA, "::", ":::", ":::", "::", ":::", 
  .     NA, "local", "::", ":::", ":::", "::", ":::", "local", "::", 
  .     "::", "::", "local", "::", "::"), error_frame = c(FALSE, 
  .     FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, TRUE, FALSE, 
  .     FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, 
  .     FALSE)), row.names = c(NA, -20L), version = 2L, class = c("rlang_trace", 
  . "rlib_trace", "tbl", "data.frame")), parent = structure(list(
  .     message = "missing value where TRUE/FALSE needed", call = if (startpars$t0 < 
  .         0) {
  .         bloodrise <- bloodrise[blood$activity >= 0.1 * startpars$peakval, 
  .             ]
  .         rise_lm <- lm(activity ~ time, data = bloodrise)
  .         rise_coef <- coef(rise_lm)
  .         startpars$t0 <- as.numeric(-1 * (rise_coef[1]/rise_coef[2]))
  .     }), class = c("simpleError", "error", "condition")), location = 1L, 
  .     name = NULL, rlang = list(inherit = TRUE), call = map(aif, 
  .         ~blmod_fengconv(time = .x$time, activity = .x$aif, Method = .x$Method)), 
  .     use_cli_format = TRUE), class = c("purrr_error_indexed", 
  . "rlang_error", "error", "condition")))
47. abort(message, class = error_class, parent = parent, call = error_call)
48. signal_abort(cnd, .file)

[mathesong]

I've seen this error message before, but didn't previously bump into it when processing this particular dataset - or at least an earlier version. I think it's running into an issue with creating approximate starting parameters while trying to figure out where the AIF rises. I don't think this is likely to be caused by the dockerisation...

When I encountered this before, it was, if I recall correctly, when I'd forgotten to include the automatic blood samples, and the model was confused by the sparseness of the measurements at the start. Are you sure that the data is complete for all individuals? Are there any individuals missing their automatic samples for instance?

If not, could you possible remove the FengConv modelling of the AIF and replace it with either a spline (preferably) or linear interpolation, and send me the html output? That way I could take a closer look at the points themselves and see if I can figure out what's going wrong.

Regardless, I should probably implement a few more checks to give a more informative error message, or else at least come up with some kind of edge-case approximation.

[pwighton]

Thanks for the help, @mathesong!

There are indeed some subjects that are missing *autosampler_blood.tsv. However even after excluding those subjects, I'm getting the same error.

Bloodstream reports are included in the attached zip file (I couldn't attach html files directly). bloodstream_report_config--aif-interp--subject-subset.html was generated using config--aif-interp--subject-subset.json. For reference bloodstream_report_config-2022-09-27_id-T3gL.html is also included, which is what I'm trying to reproduce.

bloodstream-issue.zip

[mathesong]

Looking at the html file, you can see that there are no AIF (red) datapoints before the peak. That presents a problem for all the AIF models, because they're trying to model both the ascent and descent of the AIF curve. So when they're trying to prepare for fitting the model and orient around when everything happens, they run into issues (and I should really go ahead and improve that error message).

That said, you mentioned that there were autosampler samples for these measurements? They don't appear to be being parsed. I can see several possibilities here:

• You might be subsetting your data the wrong way around (i.e. including only those without automatic samples, rather than only those with the automatic samples)
  • Could you confirm that PS19 and PS20 do have automatic samples?
• (More likely) There is an issue with parsing the data, and the tool is somehow not picking up the automatic samples. I actually recently updated this to handle some funny edge cases, but I can't think of how it could be having issues here.
  • Could you make sure that the BIDS file naming is correct, and that they completely match the manual sample filenames except for the recording-autosampler key-value pair?
  • If that's all fine, then if you boot up your docker, could you run the following and send me the outputs? That way, I can see how the parser could be failing.

library(tidyverse)
library(kinfitr)

studyfolder <- *insert study folder*

parsed_files <- kinfitr::bids_parse_files(studyfolder)

parsed_files

parsed_files %>% 
  arrange(sub, ses) %>% 
  slice(1:4) %>% 
  unnest(filedata) %>% 
  print(n=20)

parsed_files %>% 
  filter(sub %in% c("PS19", "PS20")) %>% 
  unnest(filedata) %>% 
  print(n=20)

list.files(studyfolder, pattern="_recording-autosampler_blood.tsv", recursive = T)
list.files(studyfolder, pattern="_recording-manual_blood.tsv", recursive = T)

• Lastly, if it's still unclear, I can send you a zipped up version of a minituarised version of this dataset that I made for testing bloodstream, and you can see if you get the same odd result with that. That would suggest that the dockerisation is really messing with the code function somehow - but I think that seems like an extremely remote possibility.

[mathesong]

Update:

I wanted to try out your docker container, and have been busy re-purposing an old laptop with a fresh Linux install anyhow, so I decided to give it a go. I could build the docker container without a hitch, and pulled the data using datalad (which is much more convenient because it can download the miniaturized dataset, and that's all we need for the blood). I fired up the docker image, loaded in the data, and it appears to be picking up all the automatic data just fine...

So that implies that it's probably not the docker image, nor the dataset online (through datalad at least), nor does it appear to be a code issue. The only remaining thing that I can think of is that your local dataset is compromised in some way, either through downloading through a different method and getting different data (are you also using datalad?), or else could you have an old version of the dataset?

Maybe you could try downloading the (mini) data through datalad (it's <1MB), and running bloodstream on that to see what happens?

[pwighton]

Well, that was embarrassing.. I was indeed working with an old version of the data. I downloaded the most recent version and everything is working well!

Thanks for the support @mathesong, I'll go ahead and close this issue.