Use of Covariance Matrices from R

[Philipp Riedel]

Dear Christopher,

Thanks for putting your brainGraph library on CRAN. I currently try to work with the library, but encountered a problem when importing the correlation matrices with create_mats():

myCM <- list.files(pattern = "*RestingState_CM*")

braingraph_input <- create_mats(myCM, modality = "fmri", divisor = "none", div.files = NULL, threshold.by = "consensus", sub.thresh = 0.3)

Error in scan(x, what = numeric(0), n = Nv * ncols, quiet = TRUE) : 

scan() expected 'a real', got '""'

I think it has to do with the format of the connection matrices. Instead of using a "ROICorrelation.txt" or a similar file, I actually used cor() in R to get a correlation matrix from ROI time courses. What are the requirements for the format of the connection matrices?

I'd appreciate any help, but I also know you're probably very busy.

Thanks,

Philipp

################################

Hi Philipp,

First I would like to suggest that you re-submit this question either to the GitHub page (if you have an account), or to the Google Group mailing list, if you don't mind. That way, other users will be able to see it and there will be a history of it along with other issues. If you prefer not to, that is OK.

How did you save the text files of the "cor" output? I suspect there may be some character entries in the files. I think saving with:

write.table(my_result, file="RestingState_CMX', quote=FALSE)

will work. It *might* be necessary to also include "row.names=FALSE" and "col.names=FALSE".

If you provide more details regarding your workflow, I can add that information to the User Guide if others would like to use that method.

Chris

################################

Thank you for your response! I just thought that the row and column names might be the problem and fixed it. It worked fine afterwards. I'll try to post the conversation on one of the two platforms.

[Chris Watson]

Hi Philipp, thanks for forwarding your question and my response to the list.

It looks like the issue was, indeed, because there were row and/or
column names, or there was a character value somewhere. My quoted
response explains this.

> --
> You received this message because you are subscribed to the Google Groups
> "brainGraph-help" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to brainGraph-he...@googlegroups.com.
> To post to this group, send email to brainGr...@googlegroups.com.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/brainGraph-help/541ee805-bb6f-4b0e-beab-d90f865dde32%40googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.
>

[Philipp Riedel]

Hi Chris,

I actually have some follow up questions regarding the covariance matrices (CM):

1) Do I need to replace the correlation of a ROI with itself (diagonal 1's)? If so, to 0 (e.g., CM[CM == 1] = 0) or to NA?

2) As I understand it from a different post in this group ("Adjacency matrix - does it need to be all positive values"), I need to somehow handle the negative correlation coefficients in my covariance matrices. I decided to just set them to 0 for now (e.g., CM[CM < 0] <- 0). Does that seem ok? Or should I use something like CM <- abs(CM) instead?

3) In further steps, there seem to be some, somehow different errors when I use very high (> 0.7) or very low (<0.3) thresholds (e.g., rev(seq(0.1, 0.8, 0.05))). Is that a common problem? I don't know if it matters, but I am using a customised Destrieux atlas in those further steps, that is:

destrieux.edit <- destrieux

destrieux.edit$index <- c(1:148)

destrieux.edit <- destrieux.edit[-c(24, 31, 32, 64, 98, 105, 106, 138, 8, 82, 62, 136, 4, 78, 21, 23, 33, 43, 50, 95, 97, 107, 117, 124), ]

vtp.destrieux <- setDT(destrieux.edit)

atlas <- "custom" <- 'vtp.destrieux'

Thanks for your help,

Philipp

[Philipp Riedel]

To add some information with regard to questions 1 and 2:

g.tmp <- graph_from_adjacency_matrix(opt.norm.sub[[j]][, , inds[[i]][k]], mode = 'undirected', diag = F, weighted = T)

g.tmp <- set_brainGraph_attr(g.tmp, atlas = atlas, modality = modality, weighting = 'pearson', threshold = thresholds[j], subject = covars[groups[i], ID[k]], group = groups[i], use.parallel=FALSE, A = opt.norm.sub[[j]][, , inds[[i]][k]])

On Tuesday, July 9, 2019 at 8:39:16 AM UTC-7, Philipp Riedel wrote:

[Philipp Riedel]

To add some information with regard to question 3:

thresholds <- rev(seq(0.2, 0.8, 0.05))

sub.thresh <- 0

…

> ## STEP 1 of 7.3 Graph Creation ## OPT

> 

> setwd(savedirA)

> 

> g.group <- g <- fnames <- vector('list', length=length(groups))

> 

> for (i in seq_along(groups)) {

+ for (j in seq_along(thresholds)) {

+     print(paste0('Threshold ', j, '/', length(thresholds), '; group ', i, '; ', format(Sys.time(), '%H:%M:%S')))

+ for (k in seq_along(inds[[i]])) {

+ g.tmp <- graph_from_adjacency_matrix(opt.norm.sub[[j]][, , inds[[i]][k]], mode = 'undirected', diag = F, weighted = T)

+       g.tmp <- set_brainGraph_attr(g.tmp, atlas = atlas, modality = modality, weighting = 'pearson', threshold = thresholds[j], subject = covars[groups[i], ID[k]], group = groups[i], use.parallel=FALSE, A = opt.norm.sub[[j]][, , inds[[i]][k]])

+       saveRDS(g.tmp, file = paste0(savedirA, sprintf('g%i_thr%02i_subj%03i%s', i, j, k, '.rds')))

+ }

+ }

+ 

+ # Group mean weighted graphs

+ print(paste0('Group', i, '; ', format(Sys.time(), '%H:%M:%S')))

+ g.group[[i]] <- lapply(seq_along(thresholds), function(x) graph_from_adjacency_matrix(opt.norm.mean[[x]][[i]], mode = 'undirected', diag = F, weighted = T))

+ g.group[[i]] <- llply(seq_along(thresholds), function(x) set_brainGraph_attr(g.group[[i]][[x]], atlas = atlas, modality = modality, weighting = 'pearson', threshold = thresholds[x], group = groups[i], A = opt.norm.mean[[x]][[i]], use.parallel = FALSE), .parallel = TRUE)

+ }

[1] "Threshold 1/13; group 1; 16:28:25"

[1] "Threshold 2/13; group 1; 16:31:35"

[1] "Threshold 3/13; group 1; 16:35:17"

[1] "Threshold 4/13; group 1; 16:39:08"

[1] "Threshold 5/13; group 1; 16:42:58"

[1] "Threshold 6/13; group 1; 16:46:49"

Error in if (nrow(adjmatrix) != ncol(adjmatrix)) { : 

  argument is of length zero

The code only works when I set sub.thresh to at least 0.5; that seems unusually high for fMRI data though.

My guess is that the ROIs are too highly correlated and will be completely connected above a certain (lower) threshold when I use a sub.thresh < 0.5

Thanks for your help

On Tuesday, July 9, 2019 at 8:39:16 AM UTC-7, Philipp Riedel wrote:

[Chris Watson]

Hi Philipp,

1. No, not necessary. This is because "diag=FALSE" is used throughout the package.

2. You only "need" to if you use "set_brainGraph_attrs"; in a future update, I
will handle such cases without exiting with an error but for the time being you
will have to remove negative values. This is more of a theoretical issue, and
(as far as I am aware) is an open question in the literature. I think setting
them equal to 0 is most common, but I am sure an argument could be made for
using the absolute value. It depends on whether you think that a strong negative
correlation between 2 regions implies there is a strong functional connection
between them. One issue w/ the absolute value is that you would no longer be
able to distinguish which edges had been negative and which had been positive.

3. I don't know what the issue there is. A few suggestions/points:
* Why do you change the "index" variable in the data.table? Is that how your
data (connectivity matrix) are ordered?
* The "setDT" line is not necessary; also I am not sure if you can assign to the
output of that function. You don't need a new variable/name, either; you can
just use "destrieux.edit"
* The final line should only be "atlas <- 'destrieux.edit'"

* Make sure "nrow(destrieux.edit) == nrow(A)", where "A" is one of the adjacency matrices
* After it throws the error, can you type "traceback()" and tell me what it
says? Can you also tell me the values of "i" and "j" are? It may be necessary to
look at the actual data, so if you want, you can send me a ".rda" file
containing all variables needed to reproduce the issue.

Chris

On Wed, Jul 10, 2019 at 12:18 PM, Philipp Riedel <philipp.rie...@gmail.com> wrote:

> from: Philipp Riedel <philipp.rie...@gmail.com>
> date: Wed, Jul 10 10:18 AM -07:00 2019
> to: brainGraph-help <brainGr...@googlegroups.com>
> reply-to: brainGr...@googlegroups.com
> subject: [brainGraph-help] Re: Use of Covariance Matrices from R

> --
> You received this message because you are subscribed to the Google Groups "brainGraph-help" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to brainGraph-he...@googlegroups.com.
> To post to this group, send email to brainGr...@googlegroups.com.

> To view this discussion on the web visit https://groups.google.com/d/msgid/brainGraph-help/5ae675a9-328c-4caf-b008-6eb279ef203b%40googlegroups.com.

[Philipp Riedel]

Hi Chris,

Thanks four your reply.

1) I thought so. Good to be sure.

2) Okay.

3A) Yes, I re-ordered the index variable to match the ordering in my correlation matrix. That shouldn't be a problem, should it? My edited Destrieux atlas fits the matrices, i.e. nrow(destrieux.edit) == nrow(A).

3B) This is what it says:

[1] "Threshold 1/11; group 1; 16:45:38"

[1] "Threshold 2/11; group 1; 16:49:22"

[1] "Threshold 3/11; group 1; 16:53:11"

[1] "Threshold 4/11; group 1; 16:57:07"

[1] "Threshold 5/11; group 1; 17:01:03"

Error in if (nrow(adjmatrix) != ncol(adjmatrix)) { : 

  argument is of length zero

> 

> traceback()

4: graph.adjacency.dense(adjmatrix, mode = mode, weighted = weighted, 

       diag = diag)

3: graph_from_adjacency_matrix(A[X[[i]], X[[i]]], mode = "undirected", 

       weighted = weighted)

2: efficiency(g, type = "local", use.parallel = use.parallel, A = A)

1: set_brainGraph_attr(g.tmp, atlas = atlas, modality = modality, 

       weighting = "pearson", threshold = thresholds[j], subject = covars[groups[i], 

           ID[k]], group = groups[i], use.parallel = FALSE, A = rest.norm.sub[[j]][, 

           , inds[[i]][k]])

> 

Could I maybe send you a workspace datafile with a subset of the data of list (e.g., via Email)?

Thanks,

Philipp 

On Tuesday, July 9, 2019 at 8:39:16 AM UTC-7, Philipp Riedel wrote:

[Philipp Riedel]

I forgot...

3B)

> i

[1] 1

> j

[1] 5

> thresholds[5]

[1] 0.55

On Tuesday, July 9, 2019 at 8:39:16 AM UTC-7, Philipp Riedel wrote:

[Chris Watson]

So it looks like an issue when calculating local efficiency. I have seen a
similar error from other users recently. If you send me the data I will look
into it. I will try to push a fix to GitHub this weekend if it is
straightforward enough.

On Thu, Jul 11, 2019 at 07:11 PM, Philipp Riedel <philipp.rie...@gmail.com> wrote:

> from: Philipp Riedel <philipp.rie...@gmail.com>

> date: Thu, Jul 11 05:11 PM -07:00 2019

> to: brainGraph-help <brainGr...@googlegroups.com>
> reply-to: brainGr...@googlegroups.com
> subject: [brainGraph-help] Re: Use of Covariance Matrices from R
>

> --
> You received this message because you are subscribed to the Google Groups "brainGraph-help" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to brainGraph-he...@googlegroups.com.
> To post to this group, send email to brainGr...@googlegroups.com.

> To view this discussion on the web visit https://groups.google.com/d/msgid/brainGraph-help/a8dc571f-cb5b-44fe-83ed-994012f8da33%40googlegroups.com.

[Chris Watson]

Hi Philipp, it looks like the graphs are in fact "too dense" so that an error is
thrown when calculating local efficiency. This occurs when calculating the
neighborhood of each vertex; in a fully connected graph, the neighborhood of
each vertex is the whole graph, and the local efficiency is identical for all
vertices. In this case, the value for local efficiency is equal to the global
efficiency.

It first happens at threshold 9 because 1 of the subject graphs is fully
connected, although several others are very high. I think you will have to use
higher thresholds. Even for the first threshold, the average density is ~0.74.
Another option is to increase "sub.thresh".

Chris

On Thu, Jul 11, 2019 at 11:17 PM, Chris Watson <chris.w...@gmail.com> wrote:

> from: Chris Watson <chris.w...@gmail.com>
> date: Thu, Jul 11 11:17 PM -05:00 2019

[Philipp Riedel]

Hi Chris, Thank you for your help and reply.  That's what I was afraid of.  I'm still wondering why all of the 124 ROIs are so highly correlated.  This can somehow not be explained by smoothing (5mm) alone.  What do you think about adding global signal regression to preprocessing in this case?  Otherwise, I'll go with your suggestion to use higher thresholds.  Is the 'consensus' sub.thresh option actually as common for fMRI as it is for DTI (de Reus & van den Heuvel, 2013)?  Thanks, Philipp

On Friday, July 12, 2019 at 9:33:41 PM UTC-7, Chris Watson wrote:

Hi Philipp, it looks like the graphs are in fact "too dense" so that an error is
thrown when calculating local efficiency. This occurs when calculating the
neighborhood of each vertex; in a fully connected graph, the neighborhood of
each vertex is the whole graph, and the local efficiency is identical for all
vertices. In this case, the value for local efficiency is equal to the global
efficiency.

It first happens at threshold 9 because 1 of the subject graphs is fully
connected, although several others are very high. I think you will have to use
higher thresholds. Even for the first threshold, the average density is ~0.74.
Another option is to increase "sub.thresh".

Chris

On Thu, Jul 11, 2019 at 11:17 PM, Chris Watson <chris....@gmail.com> wrote:

> from: Chris Watson <chris....@gmail.com>

> date: Thu, Jul 11 11:17 PM -05:00 2019
> to: brainGr...@googlegroups.com
> subject: [brainGraph-help] Re: Use of Covariance Matrices from R
>
> So it looks like an issue when calculating local efficiency. I have seen a
> similar error from other users recently. If you send me the data I will look
> into it. I will try to push a fix to GitHub this weekend if it is
> straightforward enough.
>

> On Thu, Jul 11, 2019 at 07:11 PM, Philipp Riedel <philipp.ri...@gmail.com> wrote:
>
>> from: Philipp Riedel <philipp.ri...@gmail.com>

>> To unsubscribe from this group and stop receiving emails from it, send an email to brainGr...@googlegroups.com.

[Chris Watson]

I am not sure about doing global signal regression. I don't read too closely the
literature for rs-fMRI preprocessing, but I have seen quite a few papers in the
past few years that look at these kinds of issues specifically. See for example:
* Chen et al., 2018, Human Brain Mapping 39(11)
* Andellini et al., 2015, J Neurosci Methods
* Bright et al., 2017, NeuroImage 154
* Murphy et al., 2017, NeuroImage 154

and references therein. I am not sure if consensus thresholding is common for
fMRI but I think it is a good idea to use one; maybe not 0.5, but something
greater than 0.

Chris

On Sat, Jul 13, 2019 at 03:49 PM, Philipp Riedel <philipp.rie...@gmail.com> wrote:

> from: Philipp Riedel <philipp.rie...@gmail.com>

> date: Sat, Jul 13 01:49 PM -07:00 2019

> to: brainGraph-help <brainGr...@googlegroups.com>
> reply-to: brainGr...@googlegroups.com

> subject: Re: [brainGraph-help] Re: Use of Covariance Matrices from R

>> <javascript:>> wrote:
>>
>> > from: Chris Watson <chris....@gmail.com <javascript:>>

>> > date: Thu, Jul 11 11:17 PM -05:00 2019

>> > to: brainGr...@googlegroups.com <javascript:>

>> > subject: [brainGraph-help] Re: Use of Covariance Matrices from R
>> >
>> > So it looks like an issue when calculating local efficiency. I have seen
>> a
>> > similar error from other users recently. If you send me the data I will
>> look
>> > into it. I will try to push a fix to GitHub this weekend if it is
>> > straightforward enough.
>> >
>> > On Thu, Jul 11, 2019 at 07:11 PM, Philipp Riedel <

>> philipp.ri...@gmail.com <javascript:>> wrote:
>> >
>> >> from: Philipp Riedel <philipp.ri...@gmail.com <javascript:>>

>> >> date: Thu, Jul 11 05:11 PM -07:00 2019

>> >> to: brainGraph-help <brainGr...@googlegroups.com <javascript:>>
>> >> reply-to: brainGr...@googlegroups.com <javascript:>

>> an email to brainGr...@googlegroups.com <javascript:>.

>> >> To post to this group, send email to brainGr...@googlegroups.com

>> <javascript:>.

>> >> To view this discussion on the web visit
>> https://groups.google.com/d/msgid/brainGraph-help/a8dc571f-cb5b-44fe-83ed-994012f8da33%40googlegroups.com.
>>
>> >> For more options, visit https://groups.google.com/d/optout.
>>
>>
>

> --
> You received this message because you are subscribed to the Google Groups "brainGraph-help" group.

> To unsubscribe from this group and stop receiving emails from it, send an email to brainGraph-he...@googlegroups.com.

> To post to this group, send email to brainGr...@googlegroups.com.

> To view this discussion on the web visit https://groups.google.com/d/msgid/brainGraph-help/59703af5-2de0-4156-a14c-8dff1ec12adb%40googlegroups.com.

[Philipp Riedel]

Thanks for your suggestions and comments!

On Saturday, July 13, 2019 at 9:50:27 PM UTC-7, Chris Watson wrote:

I am not sure about doing global signal regression. I don't read too closely the
literature for rs-fMRI preprocessing, but I have seen quite a few papers in the
past few years that look at these kinds of issues specifically. See for example:
* Chen et al., 2018, Human Brain Mapping 39(11)
* Andellini et al., 2015, J Neurosci Methods
* Bright et al., 2017, NeuroImage 154
* Murphy et al., 2017, NeuroImage 154

and references therein. I am not sure if consensus thresholding is common for
fMRI but I think it is a good idea to use one; maybe not 0.5, but something
greater than 0.

Chris

> To unsubscribe from this group and stop receiving emails from it, send an email to brainGr...@googlegroups.com.

[Philipp Riedel]

Hi Chris,

This is just an add-on for anyone who is interested in the topic.  I didn't think the consensus thresholding was the best solution for my analyses, so I read the very informative discussion you had with Sean in this group ("Threshold selection and group graph creation"; https://groups.google.com/d/msg/braingraph-help/aaxNHf7_BJg/QWTDrDTVCAAJ).

Now, instead of using the create_mats() function, I am using the below code for individual thresholding by a set of rho values.  I did not have to change the rest of your example code.  In my case, I rearranged my input files so that they are in separate directories for each group.  It worked out pretty well for me.

A1 <- brainGraph:::read.array(CM.list)

Nv <- nrow(A1)

A1[is.nan(A1)] <- 0

A.norm.sub <- lapply(seq_along(thresholds), function(z) {

  #cat("Current threshold: ", thresholds[z], "\n")

  lapply(seq_along(inds), function(x) {

    array(sapply(inds[[x]], function(y) {

      #cat(" --- current inds: ", inds[[x]][y], "\n")

      ifelse(A1[, , y] > thresholds[z], A1[, , y], 0) 

    })

    ,dim = dim(A1[, , inds[[x]]]))

  })

})

A.norm.sub <- lapply(A.norm.sub, function(x) do.call(abind::abind, x))

Thanks to you and Sean!

Best,

Philipp

On Saturday, July 13, 2019 at 9:50:27 PM UTC-7, Chris Watson wrote:

I am not sure about doing global signal regression. I don't read too closely the
literature for rs-fMRI preprocessing, but I have seen quite a few papers in the
past few years that look at these kinds of issues specifically. See for example:
* Chen et al., 2018, Human Brain Mapping 39(11)
* Andellini et al., 2015, J Neurosci Methods
* Bright et al., 2017, NeuroImage 154
* Murphy et al., 2017, NeuroImage 154

and references therein. I am not sure if consensus thresholding is common for
fMRI but I think it is a good idea to use one; maybe not 0.5, but something
greater than 0.

Chris

> To unsubscribe from this group and stop receiving emails from it, send an email to brainGr...@googlegroups.com.

[Chris Watson]

Hi Philipp,
Thanks for reporting back, and for bringing that thread up. I meant to add that
feature to the function; I will do so for the next release. It will likely be
called via "threshold.by='raw'", and the "sub.thresh" value will be ignored.
This would then have the behavior you expect; i.e., keep all edges for a given
subject if they exceed the threshold, irrespective of other subjects.

If you have any other problems in your analyses, please let me know.
Chris

On Thu, Jul 18, 2019 at 02:17 PM, Philipp Riedel <philipp.rie...@gmail.com> wrote:

> from: Philipp Riedel <philipp.rie...@gmail.com>

> date: Thu, Jul 18 12:17 PM -07:00 2019

> to: brainGraph-help <brainGr...@googlegroups.com>
> reply-to: brainGr...@googlegroups.com
> subject: Re: [brainGraph-help] Re: Use of Covariance Matrices from R
>
> Hi Chris,
>

>> <javascript:>> wrote:
>>
>> > from: Philipp Riedel <philipp.ri...@gmail.com <javascript:>>

>> > date: Sat, Jul 13 01:49 PM -07:00 2019

>> https://groups.google.com/d/msgid/brainGraph-help/59703af5-2de0-4156-a14c-8dff1ec12adb%40googlegroups.com.
>>
>> > For more options, visit https://groups.google.com/d/optout.
>>
>>
>

> --
> You received this message because you are subscribed to the Google Groups "brainGraph-help" group.

> To unsubscribe from this group and stop receiving emails from it, send an email to brainGraph-he...@googlegroups.com.

> To post to this group, send email to brainGr...@googlegroups.com.

> To view this discussion on the web visit https://groups.google.com/d/msgid/brainGraph-help/56cd6163-844b-4ec1-ad3b-7abbcc886e1d%40googlegroups.com.