Plotting BayesGLM Results

[smeisler]

Hi,

Continuing the conversation from #17 since this is now a different issue.

I have tried a few ways to plot results and get different errors. I just downloaded updated 1.9 branch.

1

plot_BayesGLM_slice(results) -->

Error in in_binary_mask[, 2:1]: incorrect number of dimensions
Traceback:

1. plot_BayesGLM_slice(results)

2

plot_slice(results$beta_estimates$'run-1') -->

ERROR while rich displaying an object: Error: geom_raster requires the following missing aesthetics: x and y

Traceback:
1. FUN(X[[i]], ...)
2. tryCatch(withCallingHandlers({
 .     if (!mime %in% names(repr::mime2repr)) 
 .         stop("No repr_* for mimetype ", mime, " in repr::mime2repr")
 .     rpr <- repr::mime2repr[[mime]](obj)
 .     if (is.null(rpr)) 
 .         return(NULL)
 .     prepare_content(is.raw(rpr), rpr)
 . }, error = error_handler), error = outer_handler)
3. tryCatchList(expr, classes, parentenv, handlers)
4. tryCatchOne(expr, names, parentenv, handlers[[1L]])
5. doTryCatch(return(expr), name, parentenv, handler)
6. withCallingHandlers({
 .     if (!mime %in% names(repr::mime2repr)) 
 .         stop("No repr_* for mimetype ", mime, " in repr::mime2repr")
 .     rpr <- repr::mime2repr[[mime]](obj)
 .     if (is.null(rpr)) 
 .         return(NULL)
 .     prepare_content(is.raw(rpr), rpr)
 . }, error = error_handler)
7. repr::mime2repr[[mime]](obj)
8. repr_text.default(obj)
9. paste(capture.output(print(obj)), collapse = "\n")
10. capture.output(print(obj))
11. withVisible(...elt(i))
12. print(obj)
13. print.ggplot(obj)
14. ggplot_build(x)
15. ggplot_build.ggplot(x)
16. by_layer(function(l, d) l$compute_geom_1(d))
17. f(l = layers[[i]], d = data[[i]])
18. l$compute_geom_1(d)
19. f(..., self = self)
20. check_required_aesthetics(self$geom$required_aes, c(names(data), 
  .     names(self$aes_params)), snake_class(self$geom))
21. abort(glue("{name} requires the following missing aesthetics: ", 
  .     glue_collapse(lapply(missing_aes, glue_collapse, sep = ", ", 
  .         last = " and "), sep = " or ")))
22. signal_abort(cnd)

3

plot(results) -->

Error in in_binary_mask[, 2:1]: incorrect number of dimensions
Traceback:

1. plot(results)
2. plot(results)
3. plot.BayesGLM(results)
4. plot_BayesGLM_slice(x, ...)

Thank you, Steven

[damondpham]

Thanks Steven!

For now, I would plot "xifti" objects directly. They are nested in the BayesGLM output. For this code:

q <- BayesGLM_cifti(
    cifti_fname = fnames_ts[ss],
    ...
)

q is a "BayesGLM_cifti" object and its entries q$betas_Bayesian and q$betas_classical are "xifti" objects. So you can visualize those coefficients with e.g. plot(q$betas_Bayesian). (plot is an S3 method for view_xifti_surface in this case, so you can look at more options for the plot with ?view_xifti_surface.)

In the future, I will be helping to make it possible for users to do plot(q) directly.

Let me know if you have further questions or suggestions!

---

EDIT: I said "xifti" objects are nested in the BayesGLM output. Actually, they are nested in the BayesGLM_cifti output, and are not returned by BayesGLM. Sorry for the confusion this caused!

[smeisler]

Thank you for the response. However, my BayesGLM_cifti object does not seem to have the xifti objects. The output of names(results) returns: 'INLA_result''beta_estimates''result_classical''mesh''mesh_orig''mask''design''session_names''beta_names''theta_posteriors''mu.theta''Q.theta''y''X''prewhiten_info''call'

Would one of these suffice?

I'm also happy to wait until the plotting functionality and output documentation have been more developed.

Steven

[damondpham]

Could you write your code which was used to get 'results', and/or attach this object?

[smeisler]

The object may be downloaded from https://drive.google.com/file/d/19oB8BGmE7XwzVkB7N6gMHMIOftIJ_ilz/view?usp=sharing

Code:

results <- BayesGLM_cifti(
    fnames_ts[1],
    surfL_fname = fname_gifti_left,
    surfR_fname = fname_gifti_right,
    brainstructures = c("left", "right"),
    onsets = events_1,
    TR = TR,
    nuisance = confounds_1,
    nuisance_include = c("drift", "dHRF"),
    scale_BOLD = TRUE,
    scale_design = TRUE,
    Bayes = TRUE,
    ar_order = 6,
    ar_smooth = 5,
    resamp_res = 10000,
    num.threads = parallel::detectCores() - 2,
    verbose = TRUE,
    outfile = outfile,
    return_INLA_result = TRUE,
    avg_sessions = TRUE,
    session_names = c('run-1'),
    meanTol = 0.1,
    varTol = 1e-06,
    trim_INLA = TRUE
)

[damondpham]

Hello Steven,

Would you have happened to do something like results <- results$GLMs_Bayesian$cortexL In between getting results and trying to plot it? The direct result of the call to BayesGLM_cifti should be a "BayesGLM_cifti" object with information for both hemispheres, but the attached file is a "BayesGLM" object with only the left hemisphere result.

I'm able to plot your results with the following code:

z <- results$beta_estimates$`run-1`
z <- ciftiTools::as.xifti(cortexL=z, cortexL_mwall=!is.na(z[,1]))
plot(z)

This code extracts the beta estimates as a (vertices x coefs) numeric matrix. It converts it to a "xifti" object with as.xifti; note how the medial wall information is given by the rows with NA values in the data matrix. Th constructed "xifti" object is then plotted.

But if you didn't do something like results <- results$GLMs_Bayesian$cortexL, we still have a mystery at hand. So please let me know if that's the case!

[mandymejia]

To built on Damon's message:

The BayesGLM_cifti function basically calls the BayesGLM function twice (once per hemisphere) then puts the resulting estimates together in xifti format. Each call to BayesGLM returns a "BayesGLM" object class. The call to BayesGLM_cifti returns a "BayesGLM_cifti" object, which is a list that includes the BayesGLM objects for each hemisphere, the xifti objects of the beta estimates, and more.

So if your object is of class BayesGLM, it must be the result of a call to the BayesGLM function, either directly or via the BayesGLM_cifti function.

On Sat, Jan 29, 2022, 1:05 AM Damon Pham @.***> wrote:

Hello Steven,

Would you have happened to do something like results <- results$GLMs_Bayesian$cortexL In between getting results and trying to plot it? The direct result of the call to BayesGLM_cifti should be a "BayesGLM_cifti" object with information for both hemispheres, but the attached file is a "BayesGLM" object with only the left hemisphere result.

I'm able to plot your results with the following code:

z <- results$beta_estimates$run-1 z <- ciftiTools::as.xifti(cortexL=z, cortexL_mwall=!is.na(z[,1])) plot(z)

This code extracts the beta estimates as a (vertices x coefs) numeric matrix. It converts it to a "xifti" object with as.xifti; note how the medial wall information is given by the rows with NA values in the data matrix. Th constructed "xifti" object is then plotted.

But if you didn't do something like results <- results$GLMs_Bayesian$cortexL, we still have a mystery at hand. So please let me know if that's the case!

— Reply to this email directly, view it on GitHub https://github.com/mandymejia/BayesfMRI/issues/20#issuecomment-1024844473, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABHVKUXQQZMYIQZH7HXZRK3UYN7SVANCNFSM5NCLY65Q . Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub.

You are receiving this because you are subscribed to this thread.Message ID: @.***>

[smeisler]

Ah, I loaded the results file with results <—RDSload(outfile_lefthemi), so I guess that is similar to indexing the left cortex results when defining results. Without rerunning the GLM again, is there a way I can use load previously computed outputs in a way that will be useful for plotting?

[smeisler]

Using the code you provided:

z <- results$beta_estimates$`run-1`
z <- ciftiTools::as.xifti(cortexL=z, cortexL_mwall=!is.na(z[,1]))
plot(z)

I am seeing results

Given that I had two columns of tasks, how should I interpret this single image?

[damondpham]

view_xifti_surface automatically plots the first column, so you can see the second column (the second set of coefficients) with plot(z, idx=2).

More information can be found at ?view_xifti_surface.

Also, you can print a quick overview with summary(z) or equivalently just z.

[smeisler]

Understood thank you. I am having trouble getting results to display when I move to RStudio (I was previously working in Jupyter). The plotting code plot(z)proceeds with no errors (just a warning about data limits since zlim is undefined), but nothing appears in the plot window or console. Simple ggplots tend to work fine. Is this something you have experienced before?

[damondpham]

If you are working interactively in RStudio, by default, view_xifti_surface should open a separate window using OpenGL (Windows and Linux) or XQuartz (Mac; you have to download XQuartz). Are you on your Linux machine still? If not, what OS are you using? There may be some idiosyncrasies related to your OS and software versions that would require some digging.

You can try circumventing this issue with the argument widget=TRUE. In RStudio, this will create an interactive htmlwidget in the Plots window.

You can also try the argument fname to save your plots to PNG files and view those.

What happens when you try to knit an RMD document which generates a plot? Try knitting something like this (replace the straight quotes with backticks):

'''{r}
library(rgl)
rgl::setupKnitr()

# Sometimes the first OpenGL window doesn't render properly.
rgl::rgl.open(); rgl::rgl.close()

my_xifti <- [...] # create the xifti object
'''

'''{r, fig.cap="My Plot", rgl=TRUE, format="png", fig.height=4.2, fig.width=5}
plot(my_xifti, idx=1, zlim=c(-5, 5), title="My Title")
'''

[smeisler]

plot(z,widget=TRUE) did the trick (popped up in the "viewer" tab, not "plot")!

To make matters confusing OS-wise, I am running Rstudio-server within a singularity container on my linux machine, and tunneling that via SSH to my local machine (Mac) 🙃. Sometimes working on university HPCs are nice. This is not one of those times.

[damondpham]

Ahh right, it should be in that tab.

Haha, that sounds confusing! Well, feel free to reach out when you encounter any more roadblocks.

[mandymejia]

Steven, do these patterns of activation make sense for your task? It might help to set zlim to (-1,1) if it's not already.

On Sat, Jan 29, 2022, 1:21 PM Damon Pham @.***> wrote:

Ahh right, it should be in that tab.

Haha, that sounds confusing! Well, feel free to reach out when you encounter any more roadblocks.

— Reply to this email directly, view it on GitHub https://github.com/mandymejia/BayesfMRI/issues/20#issuecomment-1024971363, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABHVKUUJL4WRSOOU3H7LYY3UYQ4Z5ANCNFSM5NCLY65Q . Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub.

You are receiving this because you commented.Message ID: @.***>

[smeisler]

Hello,

Preliminarily, yes, this looks like it could be right, but I would want to see a group context before making more definitive judgements. Is there an API for doing this within BayesfMRI or ciftiTools?

[mandymejia]

Hi Steven,

We do have some options for group GLM analysis in BayesfMRI. I'm not sure that's what you mean though?

Mandy

On Sun, Jan 30, 2022, 12:36 PM Steven Meisler @.***> wrote:

Hello,

Preliminarily, yes, this looks like it could be right, but I would want to see a group context before making more definitive judgements. Is there an API for doing this within BayesfMRI or ciftiTools?

— Reply to this email directly, view it on GitHub https://github.com/mandymejia/BayesfMRI/issues/20#issuecomment-1025191886, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABHVKUSF7EZBOG4FNOS6CG3UYVZITANCNFSM5NCLY65Q . Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub.

You are receiving this because you commented.Message ID: @.***>

[smeisler]

Hello,

Yes, a group GLM analysis would be great!

Steven

[smeisler]

Hi,

Just following up on this. It is unclear to me which function would perform the group analyses. Any guidance would be appreciated.

Thanks, Steven

[danieladamspencer]

Hi Steven,

The BayesGLM2 function can perform the group-level analysis, given that you have performed your subject-level analysis. You can feed in a list of BayesGLM or BayesGLM_cifti objects into the function (the documentation needs to be updated to reflect that BayesGLM_cifti objects can also be used). Take a look and let us know if you run into trouble!

All the best, Dan

On Wed, Feb 9, 2022 at 7:04 PM Steven Meisler @.***> wrote:

Hi,

Just following up on this. It is unclear to me which function would perform the group analyses. Any guidance would be appreciated.

Thanks, Steven

— Reply to this email directly, view it on GitHub https://github.com/mandymejia/BayesfMRI/issues/20#issuecomment-1034333585, or unsubscribe https://github.com/notifications/unsubscribe-auth/ADFBPK6H4VG5Y445ASIX5KTU2L6H5ANCNFSM5NCLY65Q . Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub.

You are receiving this because you are subscribed to this thread.Message ID: @.***>

[smeisler]

Got it, I will give it a shot thanks! I see there has been activity on the master branch, should I use that or keep on using 1.9?

[danieladamspencer]

Hey Steven,

I think 1.9 should be your go-to for now. We're doing a lot of testing on different branches, but 1.9 should be mostly stable.

Good luck! Dan

On Thu, Feb 10, 2022 at 12:01 PM Steven Meisler @.***> wrote:

Got it, I will give it a shot thanks! I see there has been activity on the master branch, should I use that or keep on using 1.9?

— Reply to this email directly, view it on GitHub https://github.com/mandymejia/BayesfMRI/issues/20#issuecomment-1035168833, or unsubscribe https://github.com/notifications/unsubscribe-auth/ADFBPKYPVGDCJPWTNB5I323U2PVN7ANCNFSM5NCLY65Q . Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub.

You are receiving this because you commented.Message ID: @.***>

[smeisler]

Sounds good, closing this for now so it doesn't get too off-topic.