help > Clarification on contrasts in CONN 2nd-level multivariate analysis
Showing 1-17 of 17 posts
Oct 4, 2018 09:10 AM | Martyn McFarquhar
Clarification on contrasts in CONN 2nd-level multivariate analysis
Dear Alfonso,
I am currently working with some colleagues on a data analysis using CONN. They have sent me the individual subject component maps from an MVPA analysis, and have reported 2nd-level results using these components in CONN, but I have been unable to replicate their results using my own MANOVA implementation. I was just wondering whether you could clarify exactly what CONN is doing so that I can hopefully replicate this.
For context, the analysis is based on the 2 x 2 x 2 mixed ANOVA example on page 26 of the CONN manual. There are pre- and post-treatment scans and so each subject has a set of components for the pre and a set for the post. Based on this I was assuming that CONN was running some sort of omnibus test across components, so the contrasts would be something like (using 3 components as an example):
C = [1 -1 -1 1]
M = [1 -1 0 0 0 0
0 0 1 -1 0 0
0 0 0 0 1 -1]
such that you would get a multivariate test across components on the difference between pre and post, however, this does not seem to be producing the same results in my own implementation (using Wilks Lambda as the test statistic).
Could you please clarify how the omnibus test from across components is formulated in CONN and whether I have misunderstood anything about the CONN implementation.
Best wishes,
Martyn McFarquhar
I am currently working with some colleagues on a data analysis using CONN. They have sent me the individual subject component maps from an MVPA analysis, and have reported 2nd-level results using these components in CONN, but I have been unable to replicate their results using my own MANOVA implementation. I was just wondering whether you could clarify exactly what CONN is doing so that I can hopefully replicate this.
For context, the analysis is based on the 2 x 2 x 2 mixed ANOVA example on page 26 of the CONN manual. There are pre- and post-treatment scans and so each subject has a set of components for the pre and a set for the post. Based on this I was assuming that CONN was running some sort of omnibus test across components, so the contrasts would be something like (using 3 components as an example):
C = [1 -1 -1 1]
M = [1 -1 0 0 0 0
0 0 1 -1 0 0
0 0 0 0 1 -1]
such that you would get a multivariate test across components on the difference between pre and post, however, this does not seem to be producing the same results in my own implementation (using Wilks Lambda as the test statistic).
Could you please clarify how the omnibus test from across components is formulated in CONN and whether I have misunderstood anything about the CONN implementation.
Best wishes,
Martyn McFarquhar
Oct 5, 2018 10:10 AM | Alfonso Nieto-Castanon - Boston University
RE: Clarification on contrasts in CONN 2nd-level multivariate analysis
Dear Martyn,
That looks exactly correct. If you have one within-subject factor (pre- vs. post- scans; coded as two conditions in CONN), three components from a MVPA analyses, and two additional between-subject factors (coded as four covariates, each one encoding each single-cell in a 2x2 design), then if you look at the three factor interaction the second-level analysis will be exactly as you describe (in CONN you would select the four covariates in the 'subject effects' list and enter a [1 -1 -1 1] contrast; select the two conditions in the 'conditions' list and enter a [1 -1] contrast, and select the three MVPA components in the 'sources' list and leave the default [1 0 0;0 1 0;0 0 1] contrast there). In order to have CONN produce the exact same results as your own implementation please make sure to select the 'non-parametric' statistics options in the results explorer window (so that it uses the Wilks lambda statistics there instead of SPM's reml procedure). If that is still producing different results perhaps click on the "design" button in CONN's gui to look at your second-level analysis design and C/M contrasts to make sure everything looks identical to your own analyses. Note that CONN will have reduced the design so that it enters into the second-level model 3 maps directly containing the differences between pre- and post- in each of the three MVPA components and then uses a M=[1 0 0;0 1 0; 0 0 1] contrast matrix (instead of entering the 6 original maps and then entering a M=[1 -1 0 0 0 0;0 0 1 -1 0 0;0 0 0 0 1 -1] contrast matrix), but those two approaches will lead to exactly the same Wilks lambda statistics either way, so there must be some other difference between your two versions of the analyses.
Hope this helps
Alfonso
Originally posted by Martyn McFarquhar:
That looks exactly correct. If you have one within-subject factor (pre- vs. post- scans; coded as two conditions in CONN), three components from a MVPA analyses, and two additional between-subject factors (coded as four covariates, each one encoding each single-cell in a 2x2 design), then if you look at the three factor interaction the second-level analysis will be exactly as you describe (in CONN you would select the four covariates in the 'subject effects' list and enter a [1 -1 -1 1] contrast; select the two conditions in the 'conditions' list and enter a [1 -1] contrast, and select the three MVPA components in the 'sources' list and leave the default [1 0 0;0 1 0;0 0 1] contrast there). In order to have CONN produce the exact same results as your own implementation please make sure to select the 'non-parametric' statistics options in the results explorer window (so that it uses the Wilks lambda statistics there instead of SPM's reml procedure). If that is still producing different results perhaps click on the "design" button in CONN's gui to look at your second-level analysis design and C/M contrasts to make sure everything looks identical to your own analyses. Note that CONN will have reduced the design so that it enters into the second-level model 3 maps directly containing the differences between pre- and post- in each of the three MVPA components and then uses a M=[1 0 0;0 1 0; 0 0 1] contrast matrix (instead of entering the 6 original maps and then entering a M=[1 -1 0 0 0 0;0 0 1 -1 0 0;0 0 0 0 1 -1] contrast matrix), but those two approaches will lead to exactly the same Wilks lambda statistics either way, so there must be some other difference between your two versions of the analyses.
Hope this helps
Alfonso
Originally posted by Martyn McFarquhar:
Dear Alfonso,
I am currently working with some colleagues on a data analysis using CONN. They have sent me the individual subject component maps from an MVPA analysis, and have reported 2nd-level results using these components in CONN, but I have been unable to replicate their results using my own MANOVA implementation. I was just wondering whether you could clarify exactly what CONN is doing so that I can hopefully replicate this.
For context, the analysis is based on the 2 x 2 x 2 mixed ANOVA example on page 26 of the CONN manual. There are pre- and post-treatment scans and so each subject has a set of components for the pre and a set for the post. Based on this I was assuming that CONN was running some sort of omnibus test across components, so the contrasts would be something like (using 3 components as an example):
C = [1 -1 -1 1]
M = [1 -1 0 0 0 0
0 0 1 -1 0 0
0 0 0 0 1 -1]
such that you would get a multivariate test across components on the difference between pre and post, however, this does not seem to be producing the same results in my own implementation (using Wilks Lambda as the test statistic).
Could you please clarify how the omnibus test from across components is formulated in CONN and whether I have misunderstood anything about the CONN implementation.
Best wishes,
Martyn McFarquhar
I am currently working with some colleagues on a data analysis using CONN. They have sent me the individual subject component maps from an MVPA analysis, and have reported 2nd-level results using these components in CONN, but I have been unable to replicate their results using my own MANOVA implementation. I was just wondering whether you could clarify exactly what CONN is doing so that I can hopefully replicate this.
For context, the analysis is based on the 2 x 2 x 2 mixed ANOVA example on page 26 of the CONN manual. There are pre- and post-treatment scans and so each subject has a set of components for the pre and a set for the post. Based on this I was assuming that CONN was running some sort of omnibus test across components, so the contrasts would be something like (using 3 components as an example):
C = [1 -1 -1 1]
M = [1 -1 0 0 0 0
0 0 1 -1 0 0
0 0 0 0 1 -1]
such that you would get a multivariate test across components on the difference between pre and post, however, this does not seem to be producing the same results in my own implementation (using Wilks Lambda as the test statistic).
Could you please clarify how the omnibus test from across components is formulated in CONN and whether I have misunderstood anything about the CONN implementation.
Best wishes,
Martyn McFarquhar
Oct 5, 2018 11:10 AM | Martyn McFarquhar
RE: Clarification on contrasts in CONN 2nd-level multivariate analysis
Hi Alfonso,
That's excellent, thank you very much for clarification and advice. I will forward this on to my collaborators and hopefully we can get to the bottom of the discrepancy soon. Thanks again.
Best wishes,
- Martyn
That's excellent, thank you very much for clarification and advice. I will forward this on to my collaborators and hopefully we can get to the bottom of the discrepancy soon. Thanks again.
Best wishes,
- Martyn
Oct 5, 2018 04:10 PM | Martyn McFarquhar
RE: Clarification on contrasts in CONN 2nd-level multivariate analysis
Hi Alfonso,
Just on the back of this, could you clarify what CONN would be using instead of Wilks' Lambda and how the SPM ReML approach (presumably for dealing with heterogeneity of variance) is implemented for multiple dependant variables? I can't think how this is possible given that you would get different heterogeneity estimates per-DV and so couldn't filter the design matrix in the same fashion (unless you used multiple design matrices)?
Best wishes,
- Martyn
Just on the back of this, could you clarify what CONN would be using instead of Wilks' Lambda and how the SPM ReML approach (presumably for dealing with heterogeneity of variance) is implemented for multiple dependant variables? I can't think how this is possible given that you would get different heterogeneity estimates per-DV and so couldn't filter the design matrix in the same fashion (unless you used multiple design matrices)?
Best wishes,
- Martyn
Oct 5, 2018 10:10 PM | Alfonso Nieto-Castanon - Boston University
RE: Clarification on contrasts in CONN 2nd-level multivariate analysis
Hi Martyn,
You are actually on the right track, it is very much like having multiple design matrices. To be precise, if a standard multivariate test with three dependent variables uses a model of the form:
[Y1 Y2 Y3] = X * [B1 B2 B3] + E
with Y# your Nx1 dependent variables, X your NxM design matrix and E a Nx3 noise term with N independent samples, zero mean, and arbitrary 3x3 noise covariance matrix S, SPM's ReML approach would use instead a model of the form:
[Y1 = [X 0 0 * [B1 + E
Y2 0 X 0 B2
Y3] 0 0 X] B3]
with the same NxM design matrix X (now repeated in a block-diagonal manner three times into a new 3Nx3M design matrix) and a different 3Nx1 noise term E which is now univariate with zero mean, arbitrary scalar covariance s, and non-independent samples modeled by a 3Nx3N sample-to-sample covariance of the form kron(S,I) (to be precise in SPM the scalar covariance s is allowed to be different for each voxel, while the sample-to-sample covariance kron(S,I) is assumed constant -up to the previous scaling factor- across the entire brain)
In CONN both approaches are implemented (for seed-to-voxel or voxel-to-voxel analyses). You can see the differences in the corresponding design matrices and contrasts by clicking the "design" button in the GUI and then switching in the new window in the bottom dropdown menu between "univariate model (SPM)" and "multivariate model". In the "results explorer" window, CONN will use the "multivariate model" approach for non-parametric analyses, and the "univariate model (SPM)" approach for parametric analyses. Typically, if you have a single dependent variable (or even with multiple dependent variable if you are using a vector between-conditions and between-sources contrasts) then both approaches are actually identical and produce exactly the same statistics, but when you have multiple dependents (e.g. a between-conditions contrast matrix instead of a vector) then the two models will produce (slightly) different results (mostly due to the difference in the assumption regarding spatial homogeneity of the residual covariance matrix).
Hope this helps
Alfonso
Originally posted by Martyn McFarquhar:
You are actually on the right track, it is very much like having multiple design matrices. To be precise, if a standard multivariate test with three dependent variables uses a model of the form:
[Y1 Y2 Y3] = X * [B1 B2 B3] + E
with Y# your Nx1 dependent variables, X your NxM design matrix and E a Nx3 noise term with N independent samples, zero mean, and arbitrary 3x3 noise covariance matrix S, SPM's ReML approach would use instead a model of the form:
[Y1 = [X 0 0 * [B1 + E
Y2 0 X 0 B2
Y3] 0 0 X] B3]
with the same NxM design matrix X (now repeated in a block-diagonal manner three times into a new 3Nx3M design matrix) and a different 3Nx1 noise term E which is now univariate with zero mean, arbitrary scalar covariance s, and non-independent samples modeled by a 3Nx3N sample-to-sample covariance of the form kron(S,I) (to be precise in SPM the scalar covariance s is allowed to be different for each voxel, while the sample-to-sample covariance kron(S,I) is assumed constant -up to the previous scaling factor- across the entire brain)
In CONN both approaches are implemented (for seed-to-voxel or voxel-to-voxel analyses). You can see the differences in the corresponding design matrices and contrasts by clicking the "design" button in the GUI and then switching in the new window in the bottom dropdown menu between "univariate model (SPM)" and "multivariate model". In the "results explorer" window, CONN will use the "multivariate model" approach for non-parametric analyses, and the "univariate model (SPM)" approach for parametric analyses. Typically, if you have a single dependent variable (or even with multiple dependent variable if you are using a vector between-conditions and between-sources contrasts) then both approaches are actually identical and produce exactly the same statistics, but when you have multiple dependents (e.g. a between-conditions contrast matrix instead of a vector) then the two models will produce (slightly) different results (mostly due to the difference in the assumption regarding spatial homogeneity of the residual covariance matrix).
Hope this helps
Alfonso
Originally posted by Martyn McFarquhar:
Hi Alfonso,
Just on the back of this, could you clarify what CONN would be using instead of Wilks' Lambda and how the SPM ReML approach (presumably for dealing with heterogeneity of variance) is implemented for multiple dependant variables? I can't think how this is possible given that you would get different heterogeneity estimates per-DV and so couldn't filter the design matrix in the same fashion (unless you used multiple design matrices)?
Best wishes,
- Martyn
Just on the back of this, could you clarify what CONN would be using instead of Wilks' Lambda and how the SPM ReML approach (presumably for dealing with heterogeneity of variance) is implemented for multiple dependant variables? I can't think how this is possible given that you would get different heterogeneity estimates per-DV and so couldn't filter the design matrix in the same fashion (unless you used multiple design matrices)?
Best wishes,
- Martyn
Oct 8, 2018 11:10 AM | Martyn McFarquhar
RE: Clarification on contrasts in CONN 2nd-level multivariate analysis
Hi Alfonso,
Thank you for the detailed response. So the approach is basically to form the univariate expression of the multivariate model? So for the example of 3 components with a pre a post you can have
y = reshape(Y', [n 1]);
G = kron(X, eye(6));
such that
y ~ N(G*beta, V)
and
V = blkdiag(Sigma(1), Sigma(2), Sigma(3), Sigma(4))
where each Sigma(i) is 6 x 6 covariance matrix of the dependancy structure between the original dependent variables (so the 3 components x 2 time-points)?
The SPM ReML procedure would then estimate a single covariance matrix from a pool of voxels from an initial model fit, which is then scaled by a voxel-specific noise term. So identical to a regular 2nd-level SPM analysis?
Presumably this is not suitable for DVs on different scales, but works fine for the case of the component?
The contrasts would then be something like?
C = kron(L,M)
Putting aside issues of having a global covariance estimate (which is par for the course in SPM), my question is then about the correct estimation of the variance-covariance matrix of the parameters. If you are using essentially a block-diagonal version of the original cell means design matrix, how are you able to derive the correct error terms for each contrast without inclusion of subject as a block in the design matrix?
For the regular multivariate statistics, this is implicit in the derivation of the SSCPh and SSCPe matrices, and for marginal models (such as the sandwich-estimator) this is formed from the design matrix and the estimate of V. For the regular SPM F-statistic, you have to force the correct error term due to the implicit use of the model residuals to form the singular variance estimate. This involves multiple models and the inclusion of subject (and all interactions) in the model. As you are essentially forming a univariate repeated-measures model in this approach, this is surely the only way to make sure the test statistic is correct. Is this what CONN is doing, or is there some other statistic being used that is able to accommodate the correlation structure in the derivation of the test statistic?
Best wishes,
- Martyn
Thank you for the detailed response. So the approach is basically to form the univariate expression of the multivariate model? So for the example of 3 components with a pre a post you can have
y = reshape(Y', [n 1]);
G = kron(X, eye(6));
such that
y ~ N(G*beta, V)
and
V = blkdiag(Sigma(1), Sigma(2), Sigma(3), Sigma(4))
where each Sigma(i) is 6 x 6 covariance matrix of the dependancy structure between the original dependent variables (so the 3 components x 2 time-points)?
The SPM ReML procedure would then estimate a single covariance matrix from a pool of voxels from an initial model fit, which is then scaled by a voxel-specific noise term. So identical to a regular 2nd-level SPM analysis?
Presumably this is not suitable for DVs on different scales, but works fine for the case of the component?
The contrasts would then be something like?
C = kron(L,M)
Putting aside issues of having a global covariance estimate (which is par for the course in SPM), my question is then about the correct estimation of the variance-covariance matrix of the parameters. If you are using essentially a block-diagonal version of the original cell means design matrix, how are you able to derive the correct error terms for each contrast without inclusion of subject as a block in the design matrix?
For the regular multivariate statistics, this is implicit in the derivation of the SSCPh and SSCPe matrices, and for marginal models (such as the sandwich-estimator) this is formed from the design matrix and the estimate of V. For the regular SPM F-statistic, you have to force the correct error term due to the implicit use of the model residuals to form the singular variance estimate. This involves multiple models and the inclusion of subject (and all interactions) in the model. As you are essentially forming a univariate repeated-measures model in this approach, this is surely the only way to make sure the test statistic is correct. Is this what CONN is doing, or is there some other statistic being used that is able to accommodate the correlation structure in the derivation of the test statistic?
Best wishes,
- Martyn
Oct 8, 2018 06:10 PM | Alfonso Nieto-Castanon - Boston University
RE: Clarification on contrasts in CONN 2nd-level multivariate analysis
Hi Martyn,
Yes, you are right the "univariate (SPM)" approach for multiple DVs is simply a univariate repeated measures model in SPM where both the conditions (pre- post-) and the components (MVPA) are treated as within-subject factors (with levels assumed non-independent with unequal variance). You are also right that, due to the approach chosen in CONN, a different model will be used/implemented for every different contrast that you may want to estimate.
Just for reference / additional info, in general CONN uses SPM's "two-stage" approach (described in more detail in ref attached), both for multivariate and repeated-measures univariate statistics. This means that, if you specify a contrast matrix M with K rows, then a separate image is computed for each subject and for each row of M, and the resulting NxK images are then entered into a new second-level model now with M=eye(K). In your example with 3-components (from MVPA) and 2-conditions (pre- and post-), and M=[1 -1 0 0 0 0;0 0 1 -1 0 0;0 0 0 0 1 -1] a new image will be computed separately for each subject and component estimating the difference pre- and post- in the corresponding MVPA component (so subject effects are already taken into account by this first-level contrast estimation). So, for the univariate repeated-measures approach, you would have something like:
y = reshape(Y*M', 3*N,1);
G = kron(eye(3),X);
y = G*beta + e;
with:
V = E{e*e'} = kron(Sigma,eye(N))
where Y is N-by-6 data matrix, V is a 3*N-by-3*N matrix modelling the repeated-measures structure of our data (non-independence between pairs of DVs within the same subject), the 3-by-3 matrix Sigma is estimated jointly across all voxels (note that this also allows DVs to have different scales, as long as the ratio between those scales is relatively constant across voxels), X in our example was a N-by-4 design matrix, the vector beta (in our example a 12 element vector) is estimated separately for each voxel, and the hypothesis tested is kron(eye(3),C)*beta = [0 0 0]' (M does not appear here as it has already been subsumed in the data y, and C in our original example was [1 -1 -1 1] modelling the interaction with the two between-subject effects)
In contrast, for the multivariate analyses on the same data, you would have something like:
y = Y*M';
y = X*beta + e;
with:
Sigma = E{e'*e}
where the 3-by-3 matrix Sigma, as well as beta (in our example now a 4-by-3 matrix), are estimated separately for each voxel, and the hypothesis tested now is C*beta = [0 0 0] (again with X, Y, M, and C just the same as before)
Hope this helps
Alfonso
Originally posted by Martyn McFarquhar:
Yes, you are right the "univariate (SPM)" approach for multiple DVs is simply a univariate repeated measures model in SPM where both the conditions (pre- post-) and the components (MVPA) are treated as within-subject factors (with levels assumed non-independent with unequal variance). You are also right that, due to the approach chosen in CONN, a different model will be used/implemented for every different contrast that you may want to estimate.
Just for reference / additional info, in general CONN uses SPM's "two-stage" approach (described in more detail in ref attached), both for multivariate and repeated-measures univariate statistics. This means that, if you specify a contrast matrix M with K rows, then a separate image is computed for each subject and for each row of M, and the resulting NxK images are then entered into a new second-level model now with M=eye(K). In your example with 3-components (from MVPA) and 2-conditions (pre- and post-), and M=[1 -1 0 0 0 0;0 0 1 -1 0 0;0 0 0 0 1 -1] a new image will be computed separately for each subject and component estimating the difference pre- and post- in the corresponding MVPA component (so subject effects are already taken into account by this first-level contrast estimation). So, for the univariate repeated-measures approach, you would have something like:
y = reshape(Y*M', 3*N,1);
G = kron(eye(3),X);
y = G*beta + e;
with:
V = E{e*e'} = kron(Sigma,eye(N))
where Y is N-by-6 data matrix, V is a 3*N-by-3*N matrix modelling the repeated-measures structure of our data (non-independence between pairs of DVs within the same subject), the 3-by-3 matrix Sigma is estimated jointly across all voxels (note that this also allows DVs to have different scales, as long as the ratio between those scales is relatively constant across voxels), X in our example was a N-by-4 design matrix, the vector beta (in our example a 12 element vector) is estimated separately for each voxel, and the hypothesis tested is kron(eye(3),C)*beta = [0 0 0]' (M does not appear here as it has already been subsumed in the data y, and C in our original example was [1 -1 -1 1] modelling the interaction with the two between-subject effects)
In contrast, for the multivariate analyses on the same data, you would have something like:
y = Y*M';
y = X*beta + e;
with:
Sigma = E{e'*e}
where the 3-by-3 matrix Sigma, as well as beta (in our example now a 4-by-3 matrix), are estimated separately for each voxel, and the hypothesis tested now is C*beta = [0 0 0] (again with X, Y, M, and C just the same as before)
Hope this helps
Alfonso
Originally posted by Martyn McFarquhar:
Hi Alfonso,
Thank you for the detailed response. So the approach is basically to form the univariate expression of the multivariate model? So for the example of 3 components with a pre a post you can have
y = reshape(Y', [n 1]);
G = kron(X, eye(6));
such that
y ~ N(G*beta, V)
and
V = blkdiag(Sigma(1), Sigma(2), Sigma(3), Sigma(4))
where each Sigma(i) is 6 x 6 covariance matrix of the dependancy structure between the original dependent variables (so the 3 components x 2 time-points)?
The SPM ReML procedure would then estimate a single covariance matrix from a pool of voxels from an initial model fit, which is then scaled by a voxel-specific noise term. So identical to a regular 2nd-level SPM analysis?
Presumably this is not suitable for DVs on different scales, but works fine for the case of the component?
The contrasts would then be something like?
C = kron(L,M)
Putting aside issues of having a global covariance estimate (which is par for the course in SPM), my question is then about the correct estimation of the variance-covariance matrix of the parameters. If you are using essentially a block-diagonal version of the original cell means design matrix, how are you able to derive the correct error terms for each contrast without inclusion of subject as a block in the design matrix?
For the regular multivariate statistics, this is implicit in the derivation of the SSCPh and SSCPe matrices, and for marginal models (such as the sandwich-estimator) this is formed from the design matrix and the estimate of V. For the regular SPM F-statistic, you have to force the correct error term due to the implicit use of the model residuals to form the singular variance estimate. This involves multiple models and the inclusion of subject (and all interactions) in the model. As you are essentially forming a univariate repeated-measures model in this approach, this is surely the only way to make sure the test statistic is correct. Is this what CONN is doing, or is there some other statistic being used that is able to accommodate the correlation structure in the derivation of the test statistic?
Best wishes,
- Martyn
Thank you for the detailed response. So the approach is basically to form the univariate expression of the multivariate model? So for the example of 3 components with a pre a post you can have
y = reshape(Y', [n 1]);
G = kron(X, eye(6));
such that
y ~ N(G*beta, V)
and
V = blkdiag(Sigma(1), Sigma(2), Sigma(3), Sigma(4))
where each Sigma(i) is 6 x 6 covariance matrix of the dependancy structure between the original dependent variables (so the 3 components x 2 time-points)?
The SPM ReML procedure would then estimate a single covariance matrix from a pool of voxels from an initial model fit, which is then scaled by a voxel-specific noise term. So identical to a regular 2nd-level SPM analysis?
Presumably this is not suitable for DVs on different scales, but works fine for the case of the component?
The contrasts would then be something like?
C = kron(L,M)
Putting aside issues of having a global covariance estimate (which is par for the course in SPM), my question is then about the correct estimation of the variance-covariance matrix of the parameters. If you are using essentially a block-diagonal version of the original cell means design matrix, how are you able to derive the correct error terms for each contrast without inclusion of subject as a block in the design matrix?
For the regular multivariate statistics, this is implicit in the derivation of the SSCPh and SSCPe matrices, and for marginal models (such as the sandwich-estimator) this is formed from the design matrix and the estimate of V. For the regular SPM F-statistic, you have to force the correct error term due to the implicit use of the model residuals to form the singular variance estimate. This involves multiple models and the inclusion of subject (and all interactions) in the model. As you are essentially forming a univariate repeated-measures model in this approach, this is surely the only way to make sure the test statistic is correct. Is this what CONN is doing, or is there some other statistic being used that is able to accommodate the correlation structure in the derivation of the test statistic?
Best wishes,
- Martyn
Oct 9, 2018 07:10 AM | Martyn McFarquhar
RE: Clarification on contrasts in CONN 2nd-level multivariate analysis
Hi Alfonso,
Thanks again for the detailed response, I promise the questions will stop soon!
Just to contextualise, we are still trying to resolve the discrepancy in results and my collaborator informs me that they didn't select the "non-parametric" option in CONN and as such I am assuming that the results they are reporting come from the univariate SPM approach. As such, I am just trying to understand this method fully so we can decide what to do.
I just want to check about the implementation of the "two-stage" approach as the guidance given in the attached paper (and by Will Penny on the SPM Wiki) do not actually give the correct F-statistics or degrees of freedom when you have a within-subject factor with > 2 levels. I have been in discussions with Guillaume Flandin at the FIL about this and I can send you scripts and examples to show this if you want.
For now I just want to clarify how the model works in CONN because in order to get an equivalent to the multivariate test across components (i.e. when the M contrast matrix is the identity) you would still need multiple components in the same model and thus would still require the inclusion of the subject factor in the design matrix in order to get the correct degrees of freedom and compute the correct error term for the F-tests. Is this being done in CONN?
Best wishes,
Martyn
Thanks again for the detailed response, I promise the questions will stop soon!
Just to contextualise, we are still trying to resolve the discrepancy in results and my collaborator informs me that they didn't select the "non-parametric" option in CONN and as such I am assuming that the results they are reporting come from the univariate SPM approach. As such, I am just trying to understand this method fully so we can decide what to do.
I just want to check about the implementation of the "two-stage" approach as the guidance given in the attached paper (and by Will Penny on the SPM Wiki) do not actually give the correct F-statistics or degrees of freedom when you have a within-subject factor with > 2 levels. I have been in discussions with Guillaume Flandin at the FIL about this and I can send you scripts and examples to show this if you want.
For now I just want to clarify how the model works in CONN because in order to get an equivalent to the multivariate test across components (i.e. when the M contrast matrix is the identity) you would still need multiple components in the same model and thus would still require the inclusion of the subject factor in the design matrix in order to get the correct degrees of freedom and compute the correct error term for the F-tests. Is this being done in CONN?
Best wishes,
Martyn
Oct 9, 2018 11:10 PM | Alfonso Nieto-Castanon - Boston University
RE: Clarification on contrasts in CONN 2nd-level multivariate analysis
Hi Martyn,
Thank you for the clarification (and no need to stop the questions, I appreciate the interesting points your are bringing up)
Regarding the easy portion of your question, in this particular case (MVPA analyses) the best way to proceed, in my opinion, would be to use the multivariate (non-parametric) analysis version in CONN. Beyond the point regarding the correct degrees of freedom (more on this below) the main problem for univariate SPM analyses in the context of MVPA is that the components entered into your second-level analysis are being computed using an entirely different Principal Component decomposition for each voxel, which means that, even if for nearby voxels the subspace spanned by those components is expected to be similar, the identity of the components is not (i.e. the same component across two adjacent voxels cannot be guaranteed to represent the same or similar aspect of the pattern of connectivity with the rest of the brain), so the standard SPM assumptions in second-level analyses (including homogeneity of covariance across voxels as well as random field smoothness assumptions) are not all that reasonable in this case. Hence the recommendation here is to use multivariate analyses together with non-parametric cluster statistics (i.e simply selecting "non-parametric stats" in the results explorer window).
Now, coming back to the more difficult but perhaps also more interesting portion of your question, CONN does NOT include a subject factor explicitly in the second-level design matrix, even in the context of these repeated-measure analyses. I do not believe it would be appropriate to do so, for the following reasons (but please feel free to correct me and/or let me know your thoughts): 1) in the example analysis that you describe, subject effects (average MVPA component values for each component pre- and post- intervention) are already explicitly excluded by computing the "post - pre" differences in component scores, separately for each subject, and entering those differences directly (instead of entering separately both the pre- and post- scores) into the second-level analyses; 2) the same is not done across MVPA components (i.e. entering the difference in scores between components) because we are not interested in testing whether there are any differences between the 3 MVPA component scores, but rather whether there are any effects when pooling across the 3 components; and 3) as far as I understand the ReML pooled estimation of the covariance between the DVs, when the within-subject factor is modeled as having non-independent levels, is followed by a whitening step, which (I guess assuming that the covariance is not rank-defficient) should effectively remove the dependencies between the DVs (implicitly making the degrees of freedom of the F statistic used in these analyses correct).
To illustrate the latter point (and for me to double-check what I am saying here, since I could be mistaken) I am including below some test code that simulates random data for 10 subjects and 3 non-independent DVs (e.g. simulating the three "differences in MVPA scores" variables in our example analyses), and then it runs the "univariate (SPM)"-style second-level repeated-measures analysis without the subjects factor included in the design matrix (just like CONN does when using the "univariate (SPM)" option). The point here is to illustrate that, without the whitening operation (i.e. if you were to change "W=..." in the code below to "W=eye(L*N)", in order to skip whitening) you would be exactly right that the F- statistics computed there would be incorrect and have the incorrect degrees of freedom (this is shown by the code producing non-uniform distribution of p-values in this case, and it is, as far as I can tell, simply the result of the F-test incorrect assumption about the 30 data samples in there being independent), but with the whitening operation the distribution of p-values is now correct (i.e. uniform distribution under the null hypothesis) indicating that the degrees of freedom being used in these analyses (F(3,27), rather than the F(3,9) degrees of freedom that you would get for a multivariate version of this same analysis or with the inclusion of subject factors) are correct.
In any way, please let me know your thoughts, and/or if I am misinterpreting your question (and of course I would be more than happy to hear more / learn about why the two-stage approach in SPM might not provide the correct F-statistics or degrees of freedom in the context of a within-subject factor with more than two levels), and hope this helps
Alfonso
-----------------------------------
N = 10; % number of subjects
L = 3; % number of DVs/components
% prepares design and data
X = kron(eye(L),ones(N,1)); % design matrix
C = eye(L); % contrast of interest
A = randn(L); A=chol(eye(L)+A'*A); % cholesky factor (dependence between DVs)
W = kron(sqrtm(inv(A'*A)),eye(N)); % whitening matrix
X = W*X; % whitened design matrix
% runs simulations (second-level analyses with null data)
dofe = size(X,1)-rank(X); % denominator degrees of freedom
Nc0 = rank(X*C'); % numerator degrees of freedom
iX = pinv(X'*X);
ir = pinv(C*iX*C');
iXX = iX*X';
FF = [];
for niter=1:1e6, % number of simulations
Y = randn(N,L)*A; % simulated data
Y = W*reshape(Y,[],1); % data concatenated across DVs and whitened
B = iXX*Y; % model parameters
E = Y-X*B; % model error
EE = sum(abs(E).^2,1);
h = C*B;
BB = sum(h.*(ir*h),1);
F = real(BB./max(eps,EE))*dofe/Nc0; % F-statistic
FF(niter) = F;
end
% display null distribution of p-values (expected uniform distribution if valid analysis)
p=1-spm_Fcdf(FF,Nc0,dofe);
hist(p,100)
-------------------------------------------------
Originally posted by Martyn McFarquhar:
Thank you for the clarification (and no need to stop the questions, I appreciate the interesting points your are bringing up)
Regarding the easy portion of your question, in this particular case (MVPA analyses) the best way to proceed, in my opinion, would be to use the multivariate (non-parametric) analysis version in CONN. Beyond the point regarding the correct degrees of freedom (more on this below) the main problem for univariate SPM analyses in the context of MVPA is that the components entered into your second-level analysis are being computed using an entirely different Principal Component decomposition for each voxel, which means that, even if for nearby voxels the subspace spanned by those components is expected to be similar, the identity of the components is not (i.e. the same component across two adjacent voxels cannot be guaranteed to represent the same or similar aspect of the pattern of connectivity with the rest of the brain), so the standard SPM assumptions in second-level analyses (including homogeneity of covariance across voxels as well as random field smoothness assumptions) are not all that reasonable in this case. Hence the recommendation here is to use multivariate analyses together with non-parametric cluster statistics (i.e simply selecting "non-parametric stats" in the results explorer window).
Now, coming back to the more difficult but perhaps also more interesting portion of your question, CONN does NOT include a subject factor explicitly in the second-level design matrix, even in the context of these repeated-measure analyses. I do not believe it would be appropriate to do so, for the following reasons (but please feel free to correct me and/or let me know your thoughts): 1) in the example analysis that you describe, subject effects (average MVPA component values for each component pre- and post- intervention) are already explicitly excluded by computing the "post - pre" differences in component scores, separately for each subject, and entering those differences directly (instead of entering separately both the pre- and post- scores) into the second-level analyses; 2) the same is not done across MVPA components (i.e. entering the difference in scores between components) because we are not interested in testing whether there are any differences between the 3 MVPA component scores, but rather whether there are any effects when pooling across the 3 components; and 3) as far as I understand the ReML pooled estimation of the covariance between the DVs, when the within-subject factor is modeled as having non-independent levels, is followed by a whitening step, which (I guess assuming that the covariance is not rank-defficient) should effectively remove the dependencies between the DVs (implicitly making the degrees of freedom of the F statistic used in these analyses correct).
To illustrate the latter point (and for me to double-check what I am saying here, since I could be mistaken) I am including below some test code that simulates random data for 10 subjects and 3 non-independent DVs (e.g. simulating the three "differences in MVPA scores" variables in our example analyses), and then it runs the "univariate (SPM)"-style second-level repeated-measures analysis without the subjects factor included in the design matrix (just like CONN does when using the "univariate (SPM)" option). The point here is to illustrate that, without the whitening operation (i.e. if you were to change "W=..." in the code below to "W=eye(L*N)", in order to skip whitening) you would be exactly right that the F- statistics computed there would be incorrect and have the incorrect degrees of freedom (this is shown by the code producing non-uniform distribution of p-values in this case, and it is, as far as I can tell, simply the result of the F-test incorrect assumption about the 30 data samples in there being independent), but with the whitening operation the distribution of p-values is now correct (i.e. uniform distribution under the null hypothesis) indicating that the degrees of freedom being used in these analyses (F(3,27), rather than the F(3,9) degrees of freedom that you would get for a multivariate version of this same analysis or with the inclusion of subject factors) are correct.
In any way, please let me know your thoughts, and/or if I am misinterpreting your question (and of course I would be more than happy to hear more / learn about why the two-stage approach in SPM might not provide the correct F-statistics or degrees of freedom in the context of a within-subject factor with more than two levels), and hope this helps
Alfonso
-----------------------------------
N = 10; % number of subjects
L = 3; % number of DVs/components
% prepares design and data
X = kron(eye(L),ones(N,1)); % design matrix
C = eye(L); % contrast of interest
A = randn(L); A=chol(eye(L)+A'*A); % cholesky factor (dependence between DVs)
W = kron(sqrtm(inv(A'*A)),eye(N)); % whitening matrix
X = W*X; % whitened design matrix
% runs simulations (second-level analyses with null data)
dofe = size(X,1)-rank(X); % denominator degrees of freedom
Nc0 = rank(X*C'); % numerator degrees of freedom
iX = pinv(X'*X);
ir = pinv(C*iX*C');
iXX = iX*X';
FF = [];
for niter=1:1e6, % number of simulations
Y = randn(N,L)*A; % simulated data
Y = W*reshape(Y,[],1); % data concatenated across DVs and whitened
B = iXX*Y; % model parameters
E = Y-X*B; % model error
EE = sum(abs(E).^2,1);
h = C*B;
BB = sum(h.*(ir*h),1);
F = real(BB./max(eps,EE))*dofe/Nc0; % F-statistic
FF(niter) = F;
end
% display null distribution of p-values (expected uniform distribution if valid analysis)
p=1-spm_Fcdf(FF,Nc0,dofe);
hist(p,100)
-------------------------------------------------
Originally posted by Martyn McFarquhar:
Hi Alfonso,
Thanks again for the detailed response, I promise the questions will stop soon!
Just to contextualise, we are still trying to resolve the discrepancy in results and my collaborator informs me that they didn't select the "non-parametric" option in CONN and as such I am assuming that the results they are reporting come from the univariate SPM approach. As such, I am just trying to understand this method fully so we can decide what to do.
I just want to check about the implementation of the "two-stage" approach as the guidance given in the attached paper (and by Will Penny on the SPM Wiki) do not actually give the correct F-statistics or degrees of freedom when you have a within-subject factor with > 2 levels. I have been in discussions with Guillaume Flandin at the FIL about this and I can send you scripts and examples to show this if you want.
For now I just want to clarify how the model works in CONN because in order to get an equivalent to the multivariate test across components (i.e. when the M contrast matrix is the identity) you would still need multiple components in the same model and thus would still require the inclusion of the subject factor in the design matrix in order to get the correct degrees of freedom and compute the correct error term for the F-tests. Is this being done in CONN?
Best wishes,
Martyn
Thanks again for the detailed response, I promise the questions will stop soon!
Just to contextualise, we are still trying to resolve the discrepancy in results and my collaborator informs me that they didn't select the "non-parametric" option in CONN and as such I am assuming that the results they are reporting come from the univariate SPM approach. As such, I am just trying to understand this method fully so we can decide what to do.
I just want to check about the implementation of the "two-stage" approach as the guidance given in the attached paper (and by Will Penny on the SPM Wiki) do not actually give the correct F-statistics or degrees of freedom when you have a within-subject factor with > 2 levels. I have been in discussions with Guillaume Flandin at the FIL about this and I can send you scripts and examples to show this if you want.
For now I just want to clarify how the model works in CONN because in order to get an equivalent to the multivariate test across components (i.e. when the M contrast matrix is the identity) you would still need multiple components in the same model and thus would still require the inclusion of the subject factor in the design matrix in order to get the correct degrees of freedom and compute the correct error term for the F-tests. Is this being done in CONN?
Best wishes,
Martyn
Oct 11, 2018 10:10 AM | Martyn McFarquhar
RE: Clarification on contrasts in CONN 2nd-level multivariate analysis
Hi Alfonso,
Thank you for the clarification and recommendations on the 2nd-level MVPA approach. I will pass this information on to my collaborators.
I should apologise first of all as my use of the word "correct" in terms of the F-statistic is not really accurate as you are indeed right, the whitening will guarantee (when all other assumptions are valid) that the test statistic will follow an F-distribution under the null. This is confirmed by your simulations.
The question is really one about the most appropriate F-statistic to use to test specific hypotheses, which comes down to the issue of partitioned vs pooled errors. The majority of statistical packages use the traditional "split-plot" approach for testing hypotheses in repeated-measures models because the tests are more sensitive. As such, we would do well in neuroimaging to take the same approach. I would hazard that this is the expectation of most researchers, namely that their models can match the models implemented in other packages.
Unfortunately, as far as I can tell, the two-stage approach advocated in the SPM chapter you indicated, and one the SPM Wiki (https://en.wikibooks.org/wiki/SPM/Group_...) doesn't agree with other stats packages when the within-subject factor has > 2 levels. Correct partitioned errors seem to be only achievable through the use of the flexible factorial approach (using a full over-parameterised design) or through careful combination of first-level contrasts and first-level averages.
I have attached a script and some data demonstrating this. The models looking at within-subject effects with 2-levels can actually be corrected by using two-sample t-test models for all comparisons (rather than the one-sample models advocated on the Wiki), however there is no such simple fix for those with > 3 levels, as far as I can tell. The models given for the portioned error examples are actually simplifications of the more complex over-parameterised models that I discuss in a recent preprint (https://psyarxiv.com/a5469/) but the principles are still the same.
In terms of the approach for CONN, my feeling would be that partitioned-errors are still necessary when looking across components at the second-level. Taking through differential effects of the second within-subject factor is akin to the approach advocated on the wiki, which does not result in the same statistics as produced by standard packages. The fact that you are looking across components rather than between components I don't believe should make a difference to how the data are modelled. My feeling is that the same error term that you would use for looking between components should be used when looking across components, as you are still looking at a within-subject effect (although I'd be interested to hear your take on this).
This is still a thorny issue in the community as there are questions about it almost weekly on the SPM list. I'd be interested to hear your take on this, although I realise we've gone beyond the original remit of my questions!
Best wishes,
- Martyn
Thank you for the clarification and recommendations on the 2nd-level MVPA approach. I will pass this information on to my collaborators.
I should apologise first of all as my use of the word "correct" in terms of the F-statistic is not really accurate as you are indeed right, the whitening will guarantee (when all other assumptions are valid) that the test statistic will follow an F-distribution under the null. This is confirmed by your simulations.
The question is really one about the most appropriate F-statistic to use to test specific hypotheses, which comes down to the issue of partitioned vs pooled errors. The majority of statistical packages use the traditional "split-plot" approach for testing hypotheses in repeated-measures models because the tests are more sensitive. As such, we would do well in neuroimaging to take the same approach. I would hazard that this is the expectation of most researchers, namely that their models can match the models implemented in other packages.
Unfortunately, as far as I can tell, the two-stage approach advocated in the SPM chapter you indicated, and one the SPM Wiki (https://en.wikibooks.org/wiki/SPM/Group_...) doesn't agree with other stats packages when the within-subject factor has > 2 levels. Correct partitioned errors seem to be only achievable through the use of the flexible factorial approach (using a full over-parameterised design) or through careful combination of first-level contrasts and first-level averages.
I have attached a script and some data demonstrating this. The models looking at within-subject effects with 2-levels can actually be corrected by using two-sample t-test models for all comparisons (rather than the one-sample models advocated on the Wiki), however there is no such simple fix for those with > 3 levels, as far as I can tell. The models given for the portioned error examples are actually simplifications of the more complex over-parameterised models that I discuss in a recent preprint (https://psyarxiv.com/a5469/) but the principles are still the same.
In terms of the approach for CONN, my feeling would be that partitioned-errors are still necessary when looking across components at the second-level. Taking through differential effects of the second within-subject factor is akin to the approach advocated on the wiki, which does not result in the same statistics as produced by standard packages. The fact that you are looking across components rather than between components I don't believe should make a difference to how the data are modelled. My feeling is that the same error term that you would use for looking between components should be used when looking across components, as you are still looking at a within-subject effect (although I'd be interested to hear your take on this).
This is still a thorny issue in the community as there are questions about it almost weekly on the SPM list. I'd be interested to hear your take on this, although I realise we've gone beyond the original remit of my questions!
Best wishes,
- Martyn
Oct 11, 2018 03:10 PM | Ali Amad
RE: Clarification on contrasts in CONN 2nd-level multivariate analysis
Dear Alfonso,
first of all, thank you for help on this topic and for all the clarifications.
I have worked with Martyn's team on these analysis. Following MVPA analysis, we have indeed some significant results by using parametric statistics and nothing with non-parametric stats. However, by using the results from parametric stats as seeds to run seed to voxel analysis we found several interesting results.
Do you think that all these results are worthless and cannot be used at all ? Or can these results can be discussed with caution ?
Many thanks for all you help.
Best wishes,
Ali
Originally posted by Alfonso Nieto-Castanon:
first of all, thank you for help on this topic and for all the clarifications.
I have worked with Martyn's team on these analysis. Following MVPA analysis, we have indeed some significant results by using parametric statistics and nothing with non-parametric stats. However, by using the results from parametric stats as seeds to run seed to voxel analysis we found several interesting results.
Do you think that all these results are worthless and cannot be used at all ? Or can these results can be discussed with caution ?
Many thanks for all you help.
Best wishes,
Ali
Originally posted by Alfonso Nieto-Castanon:
Hi
Martyn,
Thank you for the clarification (and no need to stop the questions, I appreciate the interesting points your are bringing up)
Regarding the easy portion of your question, in this particular case (MVPA analyses) the best way to proceed, in my opinion, would be to use the multivariate (non-parametric) analysis version in CONN. Beyond the point regarding the correct degrees of freedom (more on this below) the main problem for univariate SPM analyses in the context of MVPA is that the components entered into your second-level analysis are being computed using an entirely different Principal Component decomposition for each voxel, which means that, even if for nearby voxels the subspace spanned by those components is expected to be similar, the identity of the components is not (i.e. the same component across two adjacent voxels cannot be guaranteed to represent the same or similar aspect of the pattern of connectivity with the rest of the brain), so the standard SPM assumptions in second-level analyses (including homogeneity of covariance across voxels as well as random field smoothness assumptions) are not all that reasonable in this case. Hence the recommendation here is to use multivariate analyses together with non-parametric cluster statistics (i.e simply selecting "non-parametric stats" in the results explorer window).
Now, coming back to the more difficult but perhaps also more interesting portion of your question, CONN does NOT include a subject factor explicitly in the second-level design matrix, even in the context of these repeated-measure analyses. I do not believe it would be appropriate to do so, for the following reasons (but please feel free to correct me and/or let me know your thoughts): 1) in the example analysis that you describe, subject effects (average MVPA component values for each component pre- and post- intervention) are already explicitly excluded by computing the "post - pre" differences in component scores, separately for each subject, and entering those differences directly (instead of entering separately both the pre- and post- scores) into the second-level analyses; 2) the same is not done across MVPA components (i.e. entering the difference in scores between components) because we are not interested in testing whether there are any differences between the 3 MVPA component scores, but rather whether there are any effects when pooling across the 3 components; and 3) as far as I understand the ReML pooled estimation of the covariance between the DVs, when the within-subject factor is modeled as having non-independent levels, is followed by a whitening step, which (I guess assuming that the covariance is not rank-defficient) should effectively remove the dependencies between the DVs (implicitly making the degrees of freedom of the F statistic used in these analyses correct).
To illustrate the latter point (and for me to double-check what I am saying here, since I could be mistaken) I am including below some test code that simulates random data for 10 subjects and 3 non-independent DVs (e.g. simulating the three "differences in MVPA scores" variables in our example analyses), and then it runs the "univariate (SPM)"-style second-level repeated-measures analysis without the subjects factor included in the design matrix (just like CONN does when using the "univariate (SPM)" option). The point here is to illustrate that, without the whitening operation (i.e. if you were to change "W=..." in the code below to "W=eye(L*N)", in order to skip whitening) you would be exactly right that the F- statistics computed there would be incorrect and have the incorrect degrees of freedom (this is shown by the code producing non-uniform distribution of p-values in this case, and it is, as far as I can tell, simply the result of the F-test incorrect assumption about the 30 data samples in there being independent), but with the whitening operation the distribution of p-values is now correct (i.e. uniform distribution under the null hypothesis) indicating that the degrees of freedom being used in these analyses (F(3,27), rather than the F(3,9) degrees of freedom that you would get for a multivariate version of this same analysis or with the inclusion of subject factors) are correct.
In any way, please let me know your thoughts, and/or if I am misinterpreting your question (and of course I would be more than happy to hear more / learn about why the two-stage approach in SPM might not provide the correct F-statistics or degrees of freedom in the context of a within-subject factor with more than two levels), and hope this helps
Alfonso
-----------------------------------
N = 10; % number of subjects
L = 3; % number of DVs/components
% prepares design and data
X = kron(eye(L),ones(N,1)); % design matrix
C = eye(L); % contrast of interest
A = randn(L); A=chol(eye(L)+A'*A); % cholesky factor (dependence between DVs)
W = kron(sqrtm(inv(A'*A)),eye(N)); % whitening matrix
X = W*X; % whitened design matrix
% runs simulations (second-level analyses with null data)
dofe = size(X,1)-rank(X); % denominator degrees of freedom
Nc0 = rank(X*C'); % numerator degrees of freedom
iX = pinv(X'*X);
ir = pinv(C*iX*C');
iXX = iX*X';
FF = [];
for niter=1:1e6, % number of simulations
Y = randn(N,L)*A; % simulated data
Y = W*reshape(Y,[],1); % data concatenated across DVs and whitened
B = iXX*Y; % model parameters
E = Y-X*B; % model error
EE = sum(abs(E).^2,1);
h = C*B;
BB = sum(h.*(ir*h),1);
F = real(BB./max(eps,EE))*dofe/Nc0; % F-statistic
FF(niter) = F;
end
% display null distribution of p-values (expected uniform distribution if valid analysis)
p=1-spm_Fcdf(FF,Nc0,dofe);
hist(p,100)
-------------------------------------------------
Originally posted by Martyn McFarquhar:
Thank you for the clarification (and no need to stop the questions, I appreciate the interesting points your are bringing up)
Regarding the easy portion of your question, in this particular case (MVPA analyses) the best way to proceed, in my opinion, would be to use the multivariate (non-parametric) analysis version in CONN. Beyond the point regarding the correct degrees of freedom (more on this below) the main problem for univariate SPM analyses in the context of MVPA is that the components entered into your second-level analysis are being computed using an entirely different Principal Component decomposition for each voxel, which means that, even if for nearby voxels the subspace spanned by those components is expected to be similar, the identity of the components is not (i.e. the same component across two adjacent voxels cannot be guaranteed to represent the same or similar aspect of the pattern of connectivity with the rest of the brain), so the standard SPM assumptions in second-level analyses (including homogeneity of covariance across voxels as well as random field smoothness assumptions) are not all that reasonable in this case. Hence the recommendation here is to use multivariate analyses together with non-parametric cluster statistics (i.e simply selecting "non-parametric stats" in the results explorer window).
Now, coming back to the more difficult but perhaps also more interesting portion of your question, CONN does NOT include a subject factor explicitly in the second-level design matrix, even in the context of these repeated-measure analyses. I do not believe it would be appropriate to do so, for the following reasons (but please feel free to correct me and/or let me know your thoughts): 1) in the example analysis that you describe, subject effects (average MVPA component values for each component pre- and post- intervention) are already explicitly excluded by computing the "post - pre" differences in component scores, separately for each subject, and entering those differences directly (instead of entering separately both the pre- and post- scores) into the second-level analyses; 2) the same is not done across MVPA components (i.e. entering the difference in scores between components) because we are not interested in testing whether there are any differences between the 3 MVPA component scores, but rather whether there are any effects when pooling across the 3 components; and 3) as far as I understand the ReML pooled estimation of the covariance between the DVs, when the within-subject factor is modeled as having non-independent levels, is followed by a whitening step, which (I guess assuming that the covariance is not rank-defficient) should effectively remove the dependencies between the DVs (implicitly making the degrees of freedom of the F statistic used in these analyses correct).
To illustrate the latter point (and for me to double-check what I am saying here, since I could be mistaken) I am including below some test code that simulates random data for 10 subjects and 3 non-independent DVs (e.g. simulating the three "differences in MVPA scores" variables in our example analyses), and then it runs the "univariate (SPM)"-style second-level repeated-measures analysis without the subjects factor included in the design matrix (just like CONN does when using the "univariate (SPM)" option). The point here is to illustrate that, without the whitening operation (i.e. if you were to change "W=..." in the code below to "W=eye(L*N)", in order to skip whitening) you would be exactly right that the F- statistics computed there would be incorrect and have the incorrect degrees of freedom (this is shown by the code producing non-uniform distribution of p-values in this case, and it is, as far as I can tell, simply the result of the F-test incorrect assumption about the 30 data samples in there being independent), but with the whitening operation the distribution of p-values is now correct (i.e. uniform distribution under the null hypothesis) indicating that the degrees of freedom being used in these analyses (F(3,27), rather than the F(3,9) degrees of freedom that you would get for a multivariate version of this same analysis or with the inclusion of subject factors) are correct.
In any way, please let me know your thoughts, and/or if I am misinterpreting your question (and of course I would be more than happy to hear more / learn about why the two-stage approach in SPM might not provide the correct F-statistics or degrees of freedom in the context of a within-subject factor with more than two levels), and hope this helps
Alfonso
-----------------------------------
N = 10; % number of subjects
L = 3; % number of DVs/components
% prepares design and data
X = kron(eye(L),ones(N,1)); % design matrix
C = eye(L); % contrast of interest
A = randn(L); A=chol(eye(L)+A'*A); % cholesky factor (dependence between DVs)
W = kron(sqrtm(inv(A'*A)),eye(N)); % whitening matrix
X = W*X; % whitened design matrix
% runs simulations (second-level analyses with null data)
dofe = size(X,1)-rank(X); % denominator degrees of freedom
Nc0 = rank(X*C'); % numerator degrees of freedom
iX = pinv(X'*X);
ir = pinv(C*iX*C');
iXX = iX*X';
FF = [];
for niter=1:1e6, % number of simulations
Y = randn(N,L)*A; % simulated data
Y = W*reshape(Y,[],1); % data concatenated across DVs and whitened
B = iXX*Y; % model parameters
E = Y-X*B; % model error
EE = sum(abs(E).^2,1);
h = C*B;
BB = sum(h.*(ir*h),1);
F = real(BB./max(eps,EE))*dofe/Nc0; % F-statistic
FF(niter) = F;
end
% display null distribution of p-values (expected uniform distribution if valid analysis)
p=1-spm_Fcdf(FF,Nc0,dofe);
hist(p,100)
-------------------------------------------------
Originally posted by Martyn McFarquhar:
Hi Alfonso,
Thanks again for the detailed response, I promise the questions will stop soon!
Just to contextualise, we are still trying to resolve the discrepancy in results and my collaborator informs me that they didn't select the "non-parametric" option in CONN and as such I am assuming that the results they are reporting come from the univariate SPM approach. As such, I am just trying to understand this method fully so we can decide what to do.
I just want to check about the implementation of the "two-stage" approach as the guidance given in the attached paper (and by Will Penny on the SPM Wiki) do not actually give the correct F-statistics or degrees of freedom when you have a within-subject factor with > 2 levels. I have been in discussions with Guillaume Flandin at the FIL about this and I can send you scripts and examples to show this if you want.
For now I just want to clarify how the model works in CONN because in order to get an equivalent to the multivariate test across components (i.e. when the M contrast matrix is the identity) you would still need multiple components in the same model and thus would still require the inclusion of the subject factor in the design matrix in order to get the correct degrees of freedom and compute the correct error term for the F-tests. Is this being done in CONN?
Best wishes,
Martyn
Thanks again for the detailed response, I promise the questions will stop soon!
Just to contextualise, we are still trying to resolve the discrepancy in results and my collaborator informs me that they didn't select the "non-parametric" option in CONN and as such I am assuming that the results they are reporting come from the univariate SPM approach. As such, I am just trying to understand this method fully so we can decide what to do.
I just want to check about the implementation of the "two-stage" approach as the guidance given in the attached paper (and by Will Penny on the SPM Wiki) do not actually give the correct F-statistics or degrees of freedom when you have a within-subject factor with > 2 levels. I have been in discussions with Guillaume Flandin at the FIL about this and I can send you scripts and examples to show this if you want.
For now I just want to clarify how the model works in CONN because in order to get an equivalent to the multivariate test across components (i.e. when the M contrast matrix is the identity) you would still need multiple components in the same model and thus would still require the inclusion of the subject factor in the design matrix in order to get the correct degrees of freedom and compute the correct error term for the F-tests. Is this being done in CONN?
Best wishes,
Martyn
Oct 12, 2018 01:10 AM | Alfonso Nieto-Castanon - Boston University
RE: Clarification on contrasts in CONN 2nd-level multivariate analysis
Hi Martyn,
Thank you for the additional clarifications and scripts, you continue to raise rather interesting points. I would like to follow-up on your assertion that the two-stage approach is not optimal (it has lower sensitivity than an differently-framed partitioned errors approach) when a within-subjects factor has more than two levels. I say here "differently-framed" because, as far as I can tell, the two-stage approach in the Henson&Penny paper is simply an implementation of the same partitioned-errors approach as the one SPSS and many other software packages use (and I will try to convince you of this below).
First, regarding the scripts that you sent me, there were a couple of minor glitches there, if I am not mistaken, so just to get those out of the way (and please feel free to correct me if I am wrong):
1) in the computation of the two-level "Texture x Drink" statistics, I believe you have incorrectly defined the design matrix there. In particular, instead of:
X6 = blkdiag(ones(14,1), ones(14,1), ones(16,1), ones(16,1));
I believe it should read:
X6 = kron(blkdiag(ones(14,1),ones(16,1)),eye(2));
2) similarly, in the computation of the two-level "Location x Texture" statistics, instead of:
X7 = blkdiag(ones(28,1),ones(32,1));
I believe it should read:
X7 = kron(ones(30,1),eye(2));
3) and last, in the computation of the two-level "Location x Texture x Drink" statistics, instead of:
X8 = blkdiag(ones(14,1), ones(14,1), ones(16,1), ones(16,1));
I believe it should read:
X8 = kron(blkdiag(ones(14,1),ones(16,1)),eye(2));
With those changes, the differences between the statistics reported by the two methods (SPSS-like implementation and CONN-like implementation) are not really as dramatic as they appeared originally. Out of the 7 effects tested, the two methods report as strongly significant only the two-way within-subjects "Location x Texture" interaction, and the reported F- values across these 7 tests are also rather similar (two are identical, three are higher with the two-level approach are two are higher with SPSS's approach). Because of this I remain unconvinced that the SPSS-like approach would actually be more sensitive than the two-stage approach that CONN uses (at least given the current evidence). I also believe that the differences that we are seeing in this example are actually only due to the data not being really white/independent, so they really reflect differences in how the two approaches handle violations in their sphericity assumptions rather than intrinsic differences in the sensitivity/power of the two approaches. To illustrate that point, I will use as example the last three-way interaction tested (Location x Texture x Drink). Following your script notation, the SPSS-like approach would use something like:
% SPSS-like
Y8 = kron(eye(30),kron(eye(3),[1 -1])) * Y;
X8 = [kron(blkdiag(ones(14,1),ones(16,1)),eye(3)) kron(eye(30),ones(3,1))];
[t8, df8] = spm_ancova(X8,eye(90),Y8,[1 -1 0 -1 1 0 zeros(1,30);0 1 -1 0 -1 1 zeros(1,30)]')
which reports F=1.1671, while a CONN-like approach would use something like:
% CONN-like
Y8 = kron(eye(30),kron([-1 1 0;0 -1 1],[1 -1])) * Y;
X8 = kron(blkdiag(ones(14,1),ones(16,1)),eye(2));
[t8,df8] = spm_ancova(X8,eye(60),Y8,kron([1 -1],eye(2))')
which reports F=0.7480. The difference between these two approaches stems solely, in my opinion, from CONN's having reduced explicitly the measures across the three-level to two between-level differences before entering those to spm_ancova. To be precise, a similar two-stage approach that does not reduce the dimensionality of the 3-levels "Texture" effect would be:
% CONN-like (but with 3 dimensions)
Y8 = kron(eye(30),kron([2 -1 -1;-1 2 -1;-1 -1 2],[1 -1])) * Y;
X8 = kron(blkdiag(ones(14,1),ones(16,1)),eye(3));
[t8,df8] = spm_ancova(X8,eye(90),Y8,kron([1 -1],eye(3))')
which, unsurprisingly, reports F=1.1671 (just like the SPSS-like approach). A multivariate analysis of this 3-way interaction test would be exactly invariant to precisely this sort of differences between the two CONN-like approaches (the reduction from three linearly-dependent dimensions to two dimensions) so the reason why we are seeing those differences here must be related to the remaining dependencies between DVs. But clearly, one limitation/problem with all of these simulations here is that we have not considered that, in a real fMRI analysis scenario, the data would have been actually whitened before entering into the "spm_ancova" computation. If we go ahead and do that (of course in a real-world scenario this whitening comes from the estimation of the residual covariance across multiple voxels rather than a single-voxel, as we are considering in this example, but anyway, just for illustration purposes), we would have something like:
% SPSS-like (with actually-white residuals)
Y8 = kron(eye(30),kron(eye(3),[1 -1])) * Y;
X8 = [kron(blkdiag(ones(14,1),ones(16,1)),eye(3)) kron(eye(30),ones(3,1))];
y = Y8-X8*(X8\Y8); y=reshape(y,[],30)'; C=y'*y; W=kron(eye(30),real(sqrtm(pinv(C)))); X8=W*X8;Y8=W*Y8;
[t8, df8] = spm_ancova(X8,kron(eye(30),C),Y8,[1 -1 0 -1 1 0 zeros(1,30);0 1 -1 0 -1 1 zeros(1,30)]')
% CONN-like (with actually-white residuals)
Y8 = kron(eye(30),kron([-1 1 0;0 -1 1],[1 -1])) * Y;
X8 = kron(blkdiag(ones(14,1),ones(16,1)),eye(2));
y = Y8-X8*(X8\Y8); y=reshape(y,[],30)'; C=y'*y; W=kron(eye(30),real(sqrtm(pinv(C)))); X8=W*X8;Y8=W*Y8;
[t8,df8] = spm_ancova(X8,kron(eye(30),C),Y8,kron([1 -1],eye(2))')
% CONN-like (but with 3 dimensions with actually-white residuals)
Y8 = kron(eye(30),kron([2 -1 -1;-1 2 -1;-1 -1 2],[1 -1])) * Y;
X8 = kron(blkdiag(ones(14,1),ones(16,1)),eye(3));
y = Y8-X8*(X8\Y8); y=reshape(y,[],30)'; C=y'*y; W=kron(eye(30),real(sqrtm(pinv(C)))); X8=W*X8;Y8=W*Y8;
[t8,df8] = spm_ancova(X8,kron(eye(30),C),Y8,kron([1 -1],eye(3))')
and now all three approaches report exactly the same F=1.1006 value.
So, summarizing, SPM's "two-stage" approach used by CONN is identical, as far as I can tell, as the standard partitioned-errors approach once the data has been whitened in order to account for the repeated-measures dependency structure (as SPM&CONN do). In addition, residual differences between the approaches likely reflect remaining departures from truly-white residuals, and as such are, in my opinion, unlikely to be systematic (i.e. sensitivity should not be consistently higher for one approach vs. the other). Please let me know your thoughts/comments/suggestions
I should also probably mention that I am of course aware that this is a somewhat "thorny"/complicated issue in the SPM list. SPM gives users the flexibility to implement an incredibly wide variety of analyses. This flexibility also means that users have to make a lot of perhaps-not-obvious choices when defining their analyses using SPM (particularly in the context of full- and flexible factorial models). Unavoidably, some of those choices lead to incorrectly-defined analyses. In contrast CONN, by design, offers a much more limited range of possible analyses that users can run. The advantage is that, for those limited analyses only, we have worked out what (hopefully) is a reasonable/general/appropriate way to implement them, so the chance of potential issues or mis-specifications is also considerably smaller. That is also why I am particularly interested in making sure that the multivariate and the two-stage univariate approaches that CONN offers for mixed-design models remain the best choices available, so I am always happy to contemplate alternatives as well as to follow-up on interesting thoughts such as the ones you have been raising in this thread.
Best
Alfonso
Originally posted by Martyn McFarquhar:
Thank you for the additional clarifications and scripts, you continue to raise rather interesting points. I would like to follow-up on your assertion that the two-stage approach is not optimal (it has lower sensitivity than an differently-framed partitioned errors approach) when a within-subjects factor has more than two levels. I say here "differently-framed" because, as far as I can tell, the two-stage approach in the Henson&Penny paper is simply an implementation of the same partitioned-errors approach as the one SPSS and many other software packages use (and I will try to convince you of this below).
First, regarding the scripts that you sent me, there were a couple of minor glitches there, if I am not mistaken, so just to get those out of the way (and please feel free to correct me if I am wrong):
1) in the computation of the two-level "Texture x Drink" statistics, I believe you have incorrectly defined the design matrix there. In particular, instead of:
X6 = blkdiag(ones(14,1), ones(14,1), ones(16,1), ones(16,1));
I believe it should read:
X6 = kron(blkdiag(ones(14,1),ones(16,1)),eye(2));
2) similarly, in the computation of the two-level "Location x Texture" statistics, instead of:
X7 = blkdiag(ones(28,1),ones(32,1));
I believe it should read:
X7 = kron(ones(30,1),eye(2));
3) and last, in the computation of the two-level "Location x Texture x Drink" statistics, instead of:
X8 = blkdiag(ones(14,1), ones(14,1), ones(16,1), ones(16,1));
I believe it should read:
X8 = kron(blkdiag(ones(14,1),ones(16,1)),eye(2));
With those changes, the differences between the statistics reported by the two methods (SPSS-like implementation and CONN-like implementation) are not really as dramatic as they appeared originally. Out of the 7 effects tested, the two methods report as strongly significant only the two-way within-subjects "Location x Texture" interaction, and the reported F- values across these 7 tests are also rather similar (two are identical, three are higher with the two-level approach are two are higher with SPSS's approach). Because of this I remain unconvinced that the SPSS-like approach would actually be more sensitive than the two-stage approach that CONN uses (at least given the current evidence). I also believe that the differences that we are seeing in this example are actually only due to the data not being really white/independent, so they really reflect differences in how the two approaches handle violations in their sphericity assumptions rather than intrinsic differences in the sensitivity/power of the two approaches. To illustrate that point, I will use as example the last three-way interaction tested (Location x Texture x Drink). Following your script notation, the SPSS-like approach would use something like:
% SPSS-like
Y8 = kron(eye(30),kron(eye(3),[1 -1])) * Y;
X8 = [kron(blkdiag(ones(14,1),ones(16,1)),eye(3)) kron(eye(30),ones(3,1))];
[t8, df8] = spm_ancova(X8,eye(90),Y8,[1 -1 0 -1 1 0 zeros(1,30);0 1 -1 0 -1 1 zeros(1,30)]')
which reports F=1.1671, while a CONN-like approach would use something like:
% CONN-like
Y8 = kron(eye(30),kron([-1 1 0;0 -1 1],[1 -1])) * Y;
X8 = kron(blkdiag(ones(14,1),ones(16,1)),eye(2));
[t8,df8] = spm_ancova(X8,eye(60),Y8,kron([1 -1],eye(2))')
which reports F=0.7480. The difference between these two approaches stems solely, in my opinion, from CONN's having reduced explicitly the measures across the three-level to two between-level differences before entering those to spm_ancova. To be precise, a similar two-stage approach that does not reduce the dimensionality of the 3-levels "Texture" effect would be:
% CONN-like (but with 3 dimensions)
Y8 = kron(eye(30),kron([2 -1 -1;-1 2 -1;-1 -1 2],[1 -1])) * Y;
X8 = kron(blkdiag(ones(14,1),ones(16,1)),eye(3));
[t8,df8] = spm_ancova(X8,eye(90),Y8,kron([1 -1],eye(3))')
which, unsurprisingly, reports F=1.1671 (just like the SPSS-like approach). A multivariate analysis of this 3-way interaction test would be exactly invariant to precisely this sort of differences between the two CONN-like approaches (the reduction from three linearly-dependent dimensions to two dimensions) so the reason why we are seeing those differences here must be related to the remaining dependencies between DVs. But clearly, one limitation/problem with all of these simulations here is that we have not considered that, in a real fMRI analysis scenario, the data would have been actually whitened before entering into the "spm_ancova" computation. If we go ahead and do that (of course in a real-world scenario this whitening comes from the estimation of the residual covariance across multiple voxels rather than a single-voxel, as we are considering in this example, but anyway, just for illustration purposes), we would have something like:
% SPSS-like (with actually-white residuals)
Y8 = kron(eye(30),kron(eye(3),[1 -1])) * Y;
X8 = [kron(blkdiag(ones(14,1),ones(16,1)),eye(3)) kron(eye(30),ones(3,1))];
y = Y8-X8*(X8\Y8); y=reshape(y,[],30)'; C=y'*y; W=kron(eye(30),real(sqrtm(pinv(C)))); X8=W*X8;Y8=W*Y8;
[t8, df8] = spm_ancova(X8,kron(eye(30),C),Y8,[1 -1 0 -1 1 0 zeros(1,30);0 1 -1 0 -1 1 zeros(1,30)]')
% CONN-like (with actually-white residuals)
Y8 = kron(eye(30),kron([-1 1 0;0 -1 1],[1 -1])) * Y;
X8 = kron(blkdiag(ones(14,1),ones(16,1)),eye(2));
y = Y8-X8*(X8\Y8); y=reshape(y,[],30)'; C=y'*y; W=kron(eye(30),real(sqrtm(pinv(C)))); X8=W*X8;Y8=W*Y8;
[t8,df8] = spm_ancova(X8,kron(eye(30),C),Y8,kron([1 -1],eye(2))')
% CONN-like (but with 3 dimensions with actually-white residuals)
Y8 = kron(eye(30),kron([2 -1 -1;-1 2 -1;-1 -1 2],[1 -1])) * Y;
X8 = kron(blkdiag(ones(14,1),ones(16,1)),eye(3));
y = Y8-X8*(X8\Y8); y=reshape(y,[],30)'; C=y'*y; W=kron(eye(30),real(sqrtm(pinv(C)))); X8=W*X8;Y8=W*Y8;
[t8,df8] = spm_ancova(X8,kron(eye(30),C),Y8,kron([1 -1],eye(3))')
and now all three approaches report exactly the same F=1.1006 value.
So, summarizing, SPM's "two-stage" approach used by CONN is identical, as far as I can tell, as the standard partitioned-errors approach once the data has been whitened in order to account for the repeated-measures dependency structure (as SPM&CONN do). In addition, residual differences between the approaches likely reflect remaining departures from truly-white residuals, and as such are, in my opinion, unlikely to be systematic (i.e. sensitivity should not be consistently higher for one approach vs. the other). Please let me know your thoughts/comments/suggestions
I should also probably mention that I am of course aware that this is a somewhat "thorny"/complicated issue in the SPM list. SPM gives users the flexibility to implement an incredibly wide variety of analyses. This flexibility also means that users have to make a lot of perhaps-not-obvious choices when defining their analyses using SPM (particularly in the context of full- and flexible factorial models). Unavoidably, some of those choices lead to incorrectly-defined analyses. In contrast CONN, by design, offers a much more limited range of possible analyses that users can run. The advantage is that, for those limited analyses only, we have worked out what (hopefully) is a reasonable/general/appropriate way to implement them, so the chance of potential issues or mis-specifications is also considerably smaller. That is also why I am particularly interested in making sure that the multivariate and the two-stage univariate approaches that CONN offers for mixed-design models remain the best choices available, so I am always happy to contemplate alternatives as well as to follow-up on interesting thoughts such as the ones you have been raising in this thread.
Best
Alfonso
Originally posted by Martyn McFarquhar:
Hi Alfonso,
Thank you for the clarification and recommendations on the 2nd-level MVPA approach. I will pass this information on to my collaborators.
I should apologise first of all as my use of the word "correct" in terms of the F-statistic is not really accurate as you are indeed right, the whitening will guarantee (when all other assumptions are valid) that the test statistic will follow an F-distribution under the null. This is confirmed by your simulations.
The question is really one about the most appropriate F-statistic to use to test specific hypotheses, which comes down to the issue of partitioned vs pooled errors. The majority of statistical packages use the traditional "split-plot" approach for testing hypotheses in repeated-measures models because the tests are more sensitive. As such, we would do well in neuroimaging to take the same approach. I would hazard that this is the expectation of most researchers, namely that their models can match the models implemented in other packages.
Unfortunately, as far as I can tell, the two-stage approach advocated in the SPM chapter you indicated, and one the SPM Wiki (https://en.wikibooks.org/wiki/SPM/Group_...) doesn't agree with other stats packages when the within-subject factor has > 2 levels. Correct partitioned errors seem to be only achievable through the use of the flexible factorial approach (using a full over-parameterised design) or through careful combination of first-level contrasts and first-level averages.
I have attached a script and some data demonstrating this. The models looking at within-subject effects with 2-levels can actually be corrected by using two-sample t-test models for all comparisons (rather than the one-sample models advocated on the Wiki), however there is no such simple fix for those with > 3 levels, as far as I can tell. The models given for the portioned error examples are actually simplifications of the more complex over-parameterised models that I discuss in a recent preprint (https://psyarxiv.com/a5469/) but the principles are still the same.
In terms of the approach for CONN, my feeling would be that partitioned-errors are still necessary when looking across components at the second-level. Taking through differential effects of the second within-subject factor is akin to the approach advocated on the wiki, which does not result in the same statistics as produced by standard packages. The fact that you are looking across components rather than between components I don't believe should make a difference to how the data are modelled. My feeling is that the same error term that you would use for looking between components should be used when looking across components, as you are still looking at a within-subject effect (although I'd be interested to hear your take on this).
This is still a thorny issue in the community as there are questions about it almost weekly on the SPM list. I'd be interested to hear your take on this, although I realise we've gone beyond the original remit of my questions!
Best wishes,
- Martyn
Thank you for the clarification and recommendations on the 2nd-level MVPA approach. I will pass this information on to my collaborators.
I should apologise first of all as my use of the word "correct" in terms of the F-statistic is not really accurate as you are indeed right, the whitening will guarantee (when all other assumptions are valid) that the test statistic will follow an F-distribution under the null. This is confirmed by your simulations.
The question is really one about the most appropriate F-statistic to use to test specific hypotheses, which comes down to the issue of partitioned vs pooled errors. The majority of statistical packages use the traditional "split-plot" approach for testing hypotheses in repeated-measures models because the tests are more sensitive. As such, we would do well in neuroimaging to take the same approach. I would hazard that this is the expectation of most researchers, namely that their models can match the models implemented in other packages.
Unfortunately, as far as I can tell, the two-stage approach advocated in the SPM chapter you indicated, and one the SPM Wiki (https://en.wikibooks.org/wiki/SPM/Group_...) doesn't agree with other stats packages when the within-subject factor has > 2 levels. Correct partitioned errors seem to be only achievable through the use of the flexible factorial approach (using a full over-parameterised design) or through careful combination of first-level contrasts and first-level averages.
I have attached a script and some data demonstrating this. The models looking at within-subject effects with 2-levels can actually be corrected by using two-sample t-test models for all comparisons (rather than the one-sample models advocated on the Wiki), however there is no such simple fix for those with > 3 levels, as far as I can tell. The models given for the portioned error examples are actually simplifications of the more complex over-parameterised models that I discuss in a recent preprint (https://psyarxiv.com/a5469/) but the principles are still the same.
In terms of the approach for CONN, my feeling would be that partitioned-errors are still necessary when looking across components at the second-level. Taking through differential effects of the second within-subject factor is akin to the approach advocated on the wiki, which does not result in the same statistics as produced by standard packages. The fact that you are looking across components rather than between components I don't believe should make a difference to how the data are modelled. My feeling is that the same error term that you would use for looking between components should be used when looking across components, as you are still looking at a within-subject effect (although I'd be interested to hear your take on this).
This is still a thorny issue in the community as there are questions about it almost weekly on the SPM list. I'd be interested to hear your take on this, although I realise we've gone beyond the original remit of my questions!
Best wishes,
- Martyn
Oct 12, 2018 10:10 AM | Martyn McFarquhar
RE: Clarification on contrasts in CONN 2nd-level multivariate analysis
Hi Alfonso,
This is turning into a very interesting conversation!
Apologies for the error in the code, you are indeed correct on all counts. However, I should clarify what I meant as I was simply saying that partitioned-errors are more sensitive than pooled errors and that the two-stage approach, although claiming to implement partitioned errors, does not agree with the partitioned error approaches in other software. My point was just that we need to get the partitioned-errors right.
Even with the changes to the code, the approaches do not agree. I think we should be aiming for better than "close enough" with these methods, and the two-stage approach (based on Will Penny's recommendations) is still not estimating the F-statistic correctly nor the degrees of freedom for the within-subject models with > 2 levels. This can't be due to differences in the treatment of sphericity because the SPSS results I've reported have not adjusted for non-sphericity and nor have the SPM results. As such, the results should be identical.
My thinking is that the non-sphericity correction should not be necessary in order to achieve the "correct" F-statistic when sphericity is assumed. This speaks to a larger issue in terms of the GLM implementation in other packages that do assume sphericity (such as FSL) where the model we use should allow us to achieve the classical results, under assumptions of sphericity. If one wished to assume sphericity then the two-level models (given in the script I sent, not the ones you have indicated) do not actually achieve the classical results under assumptions of sphericity. My aim would be to get the models to agree under assumptions of sphericity and then apply the non-sphericity correction, as I don't think the correction should be seen as an integral part of the model, but rather an adjunct control for liberal/conservative conclusions when the assumptions are not adequately met.
Additionally, although the code you provided does give the same F-statistic it does not appear to give the same degrees of freedom (SPSS-like = 1.9781,55.3860 [which agrees with the Greenhouse-Geisser correction], CONN-like [2-dimensions] = 1.5814, 44.2780, CONN-like [3-dimensions] = 1.9781,55.3860 [again, Greenhouse-Geisser]). And so the 2-dimension model will give a different p-value to the SPSS and 3-dimension model. I suppose the question in then whether this is to be expected or not as to my mind all the models are testing the same hypothesis with the same data and so the conclusions should be invariant to this sort of choice? This is the claim made by the two-stage approach, but as far as I can tell it still does not agree.
Furthermore, even if the degrees of freedom agree with the Greenhouse-Geisser correction, this correction should be applied to the original F-statistic of 1.167, rather than the F-statistic of 1.1006 produced by the code. As such, the SPSS Greenhouse-Geisser corrected results of F(1.98,55.39) = 1.167, p = 0.319 is still not matched by the non-sphericity-corrected results of F(1.98,55.39) = 1.1009, p = 0.340.
I know this seems like splitting hairs, but we need to be convinced that the models we are using are appropriate and I still think the best way to do so is make sure we can achieve equivalency with the classical results in the textbooks (and implemented in other software) before adding on additional corrections (particularly when these corrections are not available in all software). This is easy enough when dealing with 2-levels, but > 2-levels is still a sticking point for me.
Looking forward to hearing your thoughts!
Best wishes,
- Martyn
This is turning into a very interesting conversation!
Apologies for the error in the code, you are indeed correct on all counts. However, I should clarify what I meant as I was simply saying that partitioned-errors are more sensitive than pooled errors and that the two-stage approach, although claiming to implement partitioned errors, does not agree with the partitioned error approaches in other software. My point was just that we need to get the partitioned-errors right.
Even with the changes to the code, the approaches do not agree. I think we should be aiming for better than "close enough" with these methods, and the two-stage approach (based on Will Penny's recommendations) is still not estimating the F-statistic correctly nor the degrees of freedom for the within-subject models with > 2 levels. This can't be due to differences in the treatment of sphericity because the SPSS results I've reported have not adjusted for non-sphericity and nor have the SPM results. As such, the results should be identical.
My thinking is that the non-sphericity correction should not be necessary in order to achieve the "correct" F-statistic when sphericity is assumed. This speaks to a larger issue in terms of the GLM implementation in other packages that do assume sphericity (such as FSL) where the model we use should allow us to achieve the classical results, under assumptions of sphericity. If one wished to assume sphericity then the two-level models (given in the script I sent, not the ones you have indicated) do not actually achieve the classical results under assumptions of sphericity. My aim would be to get the models to agree under assumptions of sphericity and then apply the non-sphericity correction, as I don't think the correction should be seen as an integral part of the model, but rather an adjunct control for liberal/conservative conclusions when the assumptions are not adequately met.
Additionally, although the code you provided does give the same F-statistic it does not appear to give the same degrees of freedom (SPSS-like = 1.9781,55.3860 [which agrees with the Greenhouse-Geisser correction], CONN-like [2-dimensions] = 1.5814, 44.2780, CONN-like [3-dimensions] = 1.9781,55.3860 [again, Greenhouse-Geisser]). And so the 2-dimension model will give a different p-value to the SPSS and 3-dimension model. I suppose the question in then whether this is to be expected or not as to my mind all the models are testing the same hypothesis with the same data and so the conclusions should be invariant to this sort of choice? This is the claim made by the two-stage approach, but as far as I can tell it still does not agree.
Furthermore, even if the degrees of freedom agree with the Greenhouse-Geisser correction, this correction should be applied to the original F-statistic of 1.167, rather than the F-statistic of 1.1006 produced by the code. As such, the SPSS Greenhouse-Geisser corrected results of F(1.98,55.39) = 1.167, p = 0.319 is still not matched by the non-sphericity-corrected results of F(1.98,55.39) = 1.1009, p = 0.340.
I know this seems like splitting hairs, but we need to be convinced that the models we are using are appropriate and I still think the best way to do so is make sure we can achieve equivalency with the classical results in the textbooks (and implemented in other software) before adding on additional corrections (particularly when these corrections are not available in all software). This is easy enough when dealing with 2-levels, but > 2-levels is still a sticking point for me.
Looking forward to hearing your thoughts!
Best wishes,
- Martyn
Oct 12, 2018 07:10 PM | Alfonso Nieto-Castanon - Boston University
RE: Clarification on contrasts in CONN 2nd-level multivariate analysis
Hi Martyn,
Fair enough, you are not convinced until seeing exactly the same results using both approaches, which sounds like healthy (and useful) scepticism so I am going to give it a try. I have taken a closer look at your scripts and these are, I believe, the additional changes there that should make all 7 statistics exactly equal across the three methods (without the need to whiten the data):
1) in the "Main effect of Location" and "Main effect of Texture" sections change the design matrix to include the between-subject factor. Basically, from:
X3 = ones(30,1);
[t3, df3] = spm_ancova(X3,[],Y3,1);
X5 = blkdiag(ones(28,1),ones(32,1));
[t5, df5] = spm_ancova(X5,[],Y5,eye(2));
to:
X3 = blkdiag(ones(14,1),ones(16,1));
[t3, df3] = spm_ancova(X3,[],Y3,[1 1]');
X5 = kron(blkdiag(ones(14,1),ones(16,1)),eye(2));
[t5, df5] = spm_ancova(X5,[]Y5,[1 0 1 0;0 1 0 1]');
2) in the "Main effect of Texture", "Texture x Drink", "Location x Texture", and "Location x Texture x Drink" sections, orthogonalize your contrast matrices (within-subject contrasts) before computing your reduced data Y, basically changing from:
Y5 = kron(eye(30),[1 1 -1 -1 0 0; 0 0 1 1 -1 -1]) * Y;
Y6 = kron(eye(30),[1 1 -1 -1 0 0; 0 0 1 1 -1 -1]) * Y;
Y7 = kron(eye(30),[1 -1 -1 1 0 0; 0 0 1 -1 -1 1]) * Y;
Y8 = kron(eye(30),[1 -1 -1 1 0 0; 0 0 1 -1 -1 1]) * Y;
to:
Y5 = kron(eye(30),orth([1 1 -1 -1 0 0; 0 0 1 1 -1 -1]')') * Y;
Y6 = kron(eye(30),orth([1 1 -1 -1 0 0; 0 0 1 1 -1 -1]')') * Y;
Y7 = kron(eye(30),orth([1 -1 -1 1 0 0; 0 0 1 -1 -1 1]')') * Y;
Y8 = kron(eye(30),orth([1 -1 -1 1 0 0; 0 0 1 -1 -1 1]')') * Y;
That should produce exactly the same F- values (as well as exactly the same degrees of freedom) across the three methods reported in your scripts (partitioned-errrors, two-stage, and SPSS). Please let me know otherwise. The reason why (1) is needed is because your alternative statistics there use a full-factorial model, so even when you are testing the presence of within-subject effects you should still include the between-subject factor if you want to get exactly the same results across the three methods (and then in your contrast specification simply define a contrast that will aggregate across all levels of that factor). The reason why (2) is needed is because none of the statistics here include any form of whitening so they are not invariant to non-rank-reducing linear transformations of your DVs (and orthogonalizing the contrast just makes sure that we are not modifying the scaling of the data that is being entered into spm_ancova, which, again, we need if we want to get exactly the same results across the three methods). To be precise CONN does not use this form of contrast orthogonalization when computing the intermediate volumes as part of CONN's two-stage implementation, which brings the good point that perhaps I should actually have CONN do that (when multiple within-subject contrasts are entered, orthogonalize the within-subject contrast matrix before computing subject-level intermediate volumes). As far as I can tell this should not make a difference in any of the two approaches used by CONN, namely multivariate analyses (since those explicitly compute the data covariance, so they are exactly invariant to non-rank-reducing linear transformations of the data) nor in univariate SPM analyses using ReML covariance modeling (since again the whitening step there will guarantee invariance to these sort of linear transformations), but there might be some subtleties involved, so I will have to think about this a bit more. Any thoughts/comments/suggestions are welcome.
So, at the very least do you think this should be enough to convince you that the two-stage procedure advocated in the Henson&Penny paper is in fact a valid implementation of M-way ANOVAs with partitioned-errors (irrespective of number of levels of a factor)?
Hope this helps, and let me know your thoughts
Alfonso
Originally posted by Martyn McFarquhar:
Fair enough, you are not convinced until seeing exactly the same results using both approaches, which sounds like healthy (and useful) scepticism so I am going to give it a try. I have taken a closer look at your scripts and these are, I believe, the additional changes there that should make all 7 statistics exactly equal across the three methods (without the need to whiten the data):
1) in the "Main effect of Location" and "Main effect of Texture" sections change the design matrix to include the between-subject factor. Basically, from:
X3 = ones(30,1);
[t3, df3] = spm_ancova(X3,[],Y3,1);
X5 = blkdiag(ones(28,1),ones(32,1));
[t5, df5] = spm_ancova(X5,[],Y5,eye(2));
to:
X3 = blkdiag(ones(14,1),ones(16,1));
[t3, df3] = spm_ancova(X3,[],Y3,[1 1]');
X5 = kron(blkdiag(ones(14,1),ones(16,1)),eye(2));
[t5, df5] = spm_ancova(X5,[]Y5,[1 0 1 0;0 1 0 1]');
2) in the "Main effect of Texture", "Texture x Drink", "Location x Texture", and "Location x Texture x Drink" sections, orthogonalize your contrast matrices (within-subject contrasts) before computing your reduced data Y, basically changing from:
Y5 = kron(eye(30),[1 1 -1 -1 0 0; 0 0 1 1 -1 -1]) * Y;
Y6 = kron(eye(30),[1 1 -1 -1 0 0; 0 0 1 1 -1 -1]) * Y;
Y7 = kron(eye(30),[1 -1 -1 1 0 0; 0 0 1 -1 -1 1]) * Y;
Y8 = kron(eye(30),[1 -1 -1 1 0 0; 0 0 1 -1 -1 1]) * Y;
to:
Y5 = kron(eye(30),orth([1 1 -1 -1 0 0; 0 0 1 1 -1 -1]')') * Y;
Y6 = kron(eye(30),orth([1 1 -1 -1 0 0; 0 0 1 1 -1 -1]')') * Y;
Y7 = kron(eye(30),orth([1 -1 -1 1 0 0; 0 0 1 -1 -1 1]')') * Y;
Y8 = kron(eye(30),orth([1 -1 -1 1 0 0; 0 0 1 -1 -1 1]')') * Y;
That should produce exactly the same F- values (as well as exactly the same degrees of freedom) across the three methods reported in your scripts (partitioned-errrors, two-stage, and SPSS). Please let me know otherwise. The reason why (1) is needed is because your alternative statistics there use a full-factorial model, so even when you are testing the presence of within-subject effects you should still include the between-subject factor if you want to get exactly the same results across the three methods (and then in your contrast specification simply define a contrast that will aggregate across all levels of that factor). The reason why (2) is needed is because none of the statistics here include any form of whitening so they are not invariant to non-rank-reducing linear transformations of your DVs (and orthogonalizing the contrast just makes sure that we are not modifying the scaling of the data that is being entered into spm_ancova, which, again, we need if we want to get exactly the same results across the three methods). To be precise CONN does not use this form of contrast orthogonalization when computing the intermediate volumes as part of CONN's two-stage implementation, which brings the good point that perhaps I should actually have CONN do that (when multiple within-subject contrasts are entered, orthogonalize the within-subject contrast matrix before computing subject-level intermediate volumes). As far as I can tell this should not make a difference in any of the two approaches used by CONN, namely multivariate analyses (since those explicitly compute the data covariance, so they are exactly invariant to non-rank-reducing linear transformations of the data) nor in univariate SPM analyses using ReML covariance modeling (since again the whitening step there will guarantee invariance to these sort of linear transformations), but there might be some subtleties involved, so I will have to think about this a bit more. Any thoughts/comments/suggestions are welcome.
So, at the very least do you think this should be enough to convince you that the two-stage procedure advocated in the Henson&Penny paper is in fact a valid implementation of M-way ANOVAs with partitioned-errors (irrespective of number of levels of a factor)?
Hope this helps, and let me know your thoughts
Alfonso
Originally posted by Martyn McFarquhar:
Hi Alfonso,
This is turning into a very interesting conversation!
Apologies for the error in the code, you are indeed correct on all counts. However, I should clarify what I meant as I was simply saying that partitioned-errors are more sensitive than pooled errors and that the two-stage approach, although claiming to implement partitioned errors, does not agree with the partitioned error approaches in other software. My point was just that we need to get the partitioned-errors right.
Even with the changes to the code, the approaches do not agree. I think we should be aiming for better than "close enough" with these methods, and the two-stage approach (based on Will Penny's recommendations) is still not estimating the F-statistic correctly nor the degrees of freedom for the within-subject models with > 2 levels. This can't be due to differences in the treatment of sphericity because the SPSS results I've reported have not adjusted for non-sphericity and nor have the SPM results. As such, the results should be identical.
My thinking is that the non-sphericity correction should not be necessary in order to achieve the "correct" F-statistic when sphericity is assumed. This speaks to a larger issue in terms of the GLM implementation in other packages that do assume sphericity (such as FSL) where the model we use should allow us to achieve the classical results, under assumptions of sphericity. If one wished to assume sphericity then the two-level models (given in the script I sent, not the ones you have indicated) do not actually achieve the classical results under assumptions of sphericity. My aim would be to get the models to agree under assumptions of sphericity and then apply the non-sphericity correction, as I don't think the correction should be seen as an integral part of the model, but rather an adjunct control for liberal/conservative conclusions when the assumptions are not adequately met.
Additionally, although the code you provided does give the same F-statistic it does not appear to give the same degrees of freedom (SPSS-like = 1.9781,55.3860 [which agrees with the Greenhouse-Geisser correction], CONN-like [2-dimensions] = 1.5814, 44.2780, CONN-like [3-dimensions] = 1.9781,55.3860 [again, Greenhouse-Geisser]). And so the 2-dimension model will give a different p-value to the SPSS and 3-dimension model. I suppose the question in then whether this is to be expected or not as to my mind all the models are testing the same hypothesis with the same data and so the conclusions should be invariant to this sort of choice? This is the claim made by the two-stage approach, but as far as I can tell it still does not agree.
Furthermore, even if the degrees of freedom agree with the Greenhouse-Geisser correction, this correction should be applied to the original F-statistic of 1.167, rather than the F-statistic of 1.1006 produced by the code. As such, the SPSS Greenhouse-Geisser corrected results of F(1.98,55.39) = 1.167, p = 0.319 is still not matched by the non-sphericity-corrected results of F(1.98,55.39) = 1.1009, p = 0.340.
I know this seems like splitting hairs, but we need to be convinced that the models we are using are appropriate and I still think the best way to do so is make sure we can achieve equivalency with the classical results in the textbooks (and implemented in other software) before adding on additional corrections (particularly when these corrections are not available in all software). This is easy enough when dealing with 2-levels, but > 2-levels is still a sticking point for me.
Looking forward to hearing your thoughts!
Best wishes,
- Martyn
This is turning into a very interesting conversation!
Apologies for the error in the code, you are indeed correct on all counts. However, I should clarify what I meant as I was simply saying that partitioned-errors are more sensitive than pooled errors and that the two-stage approach, although claiming to implement partitioned errors, does not agree with the partitioned error approaches in other software. My point was just that we need to get the partitioned-errors right.
Even with the changes to the code, the approaches do not agree. I think we should be aiming for better than "close enough" with these methods, and the two-stage approach (based on Will Penny's recommendations) is still not estimating the F-statistic correctly nor the degrees of freedom for the within-subject models with > 2 levels. This can't be due to differences in the treatment of sphericity because the SPSS results I've reported have not adjusted for non-sphericity and nor have the SPM results. As such, the results should be identical.
My thinking is that the non-sphericity correction should not be necessary in order to achieve the "correct" F-statistic when sphericity is assumed. This speaks to a larger issue in terms of the GLM implementation in other packages that do assume sphericity (such as FSL) where the model we use should allow us to achieve the classical results, under assumptions of sphericity. If one wished to assume sphericity then the two-level models (given in the script I sent, not the ones you have indicated) do not actually achieve the classical results under assumptions of sphericity. My aim would be to get the models to agree under assumptions of sphericity and then apply the non-sphericity correction, as I don't think the correction should be seen as an integral part of the model, but rather an adjunct control for liberal/conservative conclusions when the assumptions are not adequately met.
Additionally, although the code you provided does give the same F-statistic it does not appear to give the same degrees of freedom (SPSS-like = 1.9781,55.3860 [which agrees with the Greenhouse-Geisser correction], CONN-like [2-dimensions] = 1.5814, 44.2780, CONN-like [3-dimensions] = 1.9781,55.3860 [again, Greenhouse-Geisser]). And so the 2-dimension model will give a different p-value to the SPSS and 3-dimension model. I suppose the question in then whether this is to be expected or not as to my mind all the models are testing the same hypothesis with the same data and so the conclusions should be invariant to this sort of choice? This is the claim made by the two-stage approach, but as far as I can tell it still does not agree.
Furthermore, even if the degrees of freedom agree with the Greenhouse-Geisser correction, this correction should be applied to the original F-statistic of 1.167, rather than the F-statistic of 1.1006 produced by the code. As such, the SPSS Greenhouse-Geisser corrected results of F(1.98,55.39) = 1.167, p = 0.319 is still not matched by the non-sphericity-corrected results of F(1.98,55.39) = 1.1009, p = 0.340.
I know this seems like splitting hairs, but we need to be convinced that the models we are using are appropriate and I still think the best way to do so is make sure we can achieve equivalency with the classical results in the textbooks (and implemented in other software) before adding on additional corrections (particularly when these corrections are not available in all software). This is easy enough when dealing with 2-levels, but > 2-levels is still a sticking point for me.
Looking forward to hearing your thoughts!
Best wishes,
- Martyn
Oct 12, 2018 11:10 PM | Alfonso Nieto-Castanon - Boston University
RE: Clarification on contrasts in CONN 2nd-level multivariate analysis
Dear Ali,
Sorry I had missed this message. Yes, I would definitely report those results and simply mention that the conclusions should be treated with caution until replicated given that the results do not survive a more conservative non-parametric randomization test.
Hope this helps
Alfonso
Originally posted by Ali Amad:
Sorry I had missed this message. Yes, I would definitely report those results and simply mention that the conclusions should be treated with caution until replicated given that the results do not survive a more conservative non-parametric randomization test.
Hope this helps
Alfonso
Originally posted by Ali Amad:
Dear
Alfonso,
first of all, thank you for help on this topic and for all the clarifications.
I have worked with Martyn's team on these analysis. Following MVPA analysis, we have indeed some significant results by using parametric statistics and nothing with non-parametric stats. However, by using the results from parametric stats as seeds to run seed to voxel analysis we found several interesting results.
Do you think that all these results are worthless and cannot be used at all ? Or can these results can be discussed with caution ?
Many thanks for all you help.
Best wishes,
Ali
Originally posted by Alfonso Nieto-Castanon:
first of all, thank you for help on this topic and for all the clarifications.
I have worked with Martyn's team on these analysis. Following MVPA analysis, we have indeed some significant results by using parametric statistics and nothing with non-parametric stats. However, by using the results from parametric stats as seeds to run seed to voxel analysis we found several interesting results.
Do you think that all these results are worthless and cannot be used at all ? Or can these results can be discussed with caution ?
Many thanks for all you help.
Best wishes,
Ali
Originally posted by Alfonso Nieto-Castanon:
Hi
Martyn,
Thank you for the clarification (and no need to stop the questions, I appreciate the interesting points your are bringing up)
Regarding the easy portion of your question, in this particular case (MVPA analyses) the best way to proceed, in my opinion, would be to use the multivariate (non-parametric) analysis version in CONN. Beyond the point regarding the correct degrees of freedom (more on this below) the main problem for univariate SPM analyses in the context of MVPA is that the components entered into your second-level analysis are being computed using an entirely different Principal Component decomposition for each voxel, which means that, even if for nearby voxels the subspace spanned by those components is expected to be similar, the identity of the components is not (i.e. the same component across two adjacent voxels cannot be guaranteed to represent the same or similar aspect of the pattern of connectivity with the rest of the brain), so the standard SPM assumptions in second-level analyses (including homogeneity of covariance across voxels as well as random field smoothness assumptions) are not all that reasonable in this case. Hence the recommendation here is to use multivariate analyses together with non-parametric cluster statistics (i.e simply selecting "non-parametric stats" in the results explorer window).
Now, coming back to the more difficult but perhaps also more interesting portion of your question, CONN does NOT include a subject factor explicitly in the second-level design matrix, even in the context of these repeated-measure analyses. I do not believe it would be appropriate to do so, for the following reasons (but please feel free to correct me and/or let me know your thoughts): 1) in the example analysis that you describe, subject effects (average MVPA component values for each component pre- and post- intervention) are already explicitly excluded by computing the "post - pre" differences in component scores, separately for each subject, and entering those differences directly (instead of entering separately both the pre- and post- scores) into the second-level analyses; 2) the same is not done across MVPA components (i.e. entering the difference in scores between components) because we are not interested in testing whether there are any differences between the 3 MVPA component scores, but rather whether there are any effects when pooling across the 3 components; and 3) as far as I understand the ReML pooled estimation of the covariance between the DVs, when the within-subject factor is modeled as having non-independent levels, is followed by a whitening step, which (I guess assuming that the covariance is not rank-defficient) should effectively remove the dependencies between the DVs (implicitly making the degrees of freedom of the F statistic used in these analyses correct).
To illustrate the latter point (and for me to double-check what I am saying here, since I could be mistaken) I am including below some test code that simulates random data for 10 subjects and 3 non-independent DVs (e.g. simulating the three "differences in MVPA scores" variables in our example analyses), and then it runs the "univariate (SPM)"-style second-level repeated-measures analysis without the subjects factor included in the design matrix (just like CONN does when using the "univariate (SPM)" option). The point here is to illustrate that, without the whitening operation (i.e. if you were to change "W=..." in the code below to "W=eye(L*N)", in order to skip whitening) you would be exactly right that the F- statistics computed there would be incorrect and have the incorrect degrees of freedom (this is shown by the code producing non-uniform distribution of p-values in this case, and it is, as far as I can tell, simply the result of the F-test incorrect assumption about the 30 data samples in there being independent), but with the whitening operation the distribution of p-values is now correct (i.e. uniform distribution under the null hypothesis) indicating that the degrees of freedom being used in these analyses (F(3,27), rather than the F(3,9) degrees of freedom that you would get for a multivariate version of this same analysis or with the inclusion of subject factors) are correct.
In any way, please let me know your thoughts, and/or if I am misinterpreting your question (and of course I would be more than happy to hear more / learn about why the two-stage approach in SPM might not provide the correct F-statistics or degrees of freedom in the context of a within-subject factor with more than two levels), and hope this helps
Alfonso
-----------------------------------
N = 10; % number of subjects
L = 3; % number of DVs/components
% prepares design and data
X = kron(eye(L),ones(N,1)); % design matrix
C = eye(L); % contrast of interest
A = randn(L); A=chol(eye(L)+A'*A); % cholesky factor (dependence between DVs)
W = kron(sqrtm(inv(A'*A)),eye(N)); % whitening matrix
X = W*X; % whitened design matrix
% runs simulations (second-level analyses with null data)
dofe = size(X,1)-rank(X); % denominator degrees of freedom
Nc0 = rank(X*C'); % numerator degrees of freedom
iX = pinv(X'*X);
ir = pinv(C*iX*C');
iXX = iX*X';
FF = [];
for niter=1:1e6, % number of simulations
Y = randn(N,L)*A; % simulated data
Y = W*reshape(Y,[],1); % data concatenated across DVs and whitened
B = iXX*Y; % model parameters
E = Y-X*B; % model error
EE = sum(abs(E).^2,1);
h = C*B;
BB = sum(h.*(ir*h),1);
F = real(BB./max(eps,EE))*dofe/Nc0; % F-statistic
FF(niter) = F;
end
% display null distribution of p-values (expected uniform distribution if valid analysis)
p=1-spm_Fcdf(FF,Nc0,dofe);
hist(p,100)
-------------------------------------------------
Originally posted by Martyn McFarquhar:
Thank you for the clarification (and no need to stop the questions, I appreciate the interesting points your are bringing up)
Regarding the easy portion of your question, in this particular case (MVPA analyses) the best way to proceed, in my opinion, would be to use the multivariate (non-parametric) analysis version in CONN. Beyond the point regarding the correct degrees of freedom (more on this below) the main problem for univariate SPM analyses in the context of MVPA is that the components entered into your second-level analysis are being computed using an entirely different Principal Component decomposition for each voxel, which means that, even if for nearby voxels the subspace spanned by those components is expected to be similar, the identity of the components is not (i.e. the same component across two adjacent voxels cannot be guaranteed to represent the same or similar aspect of the pattern of connectivity with the rest of the brain), so the standard SPM assumptions in second-level analyses (including homogeneity of covariance across voxels as well as random field smoothness assumptions) are not all that reasonable in this case. Hence the recommendation here is to use multivariate analyses together with non-parametric cluster statistics (i.e simply selecting "non-parametric stats" in the results explorer window).
Now, coming back to the more difficult but perhaps also more interesting portion of your question, CONN does NOT include a subject factor explicitly in the second-level design matrix, even in the context of these repeated-measure analyses. I do not believe it would be appropriate to do so, for the following reasons (but please feel free to correct me and/or let me know your thoughts): 1) in the example analysis that you describe, subject effects (average MVPA component values for each component pre- and post- intervention) are already explicitly excluded by computing the "post - pre" differences in component scores, separately for each subject, and entering those differences directly (instead of entering separately both the pre- and post- scores) into the second-level analyses; 2) the same is not done across MVPA components (i.e. entering the difference in scores between components) because we are not interested in testing whether there are any differences between the 3 MVPA component scores, but rather whether there are any effects when pooling across the 3 components; and 3) as far as I understand the ReML pooled estimation of the covariance between the DVs, when the within-subject factor is modeled as having non-independent levels, is followed by a whitening step, which (I guess assuming that the covariance is not rank-defficient) should effectively remove the dependencies between the DVs (implicitly making the degrees of freedom of the F statistic used in these analyses correct).
To illustrate the latter point (and for me to double-check what I am saying here, since I could be mistaken) I am including below some test code that simulates random data for 10 subjects and 3 non-independent DVs (e.g. simulating the three "differences in MVPA scores" variables in our example analyses), and then it runs the "univariate (SPM)"-style second-level repeated-measures analysis without the subjects factor included in the design matrix (just like CONN does when using the "univariate (SPM)" option). The point here is to illustrate that, without the whitening operation (i.e. if you were to change "W=..." in the code below to "W=eye(L*N)", in order to skip whitening) you would be exactly right that the F- statistics computed there would be incorrect and have the incorrect degrees of freedom (this is shown by the code producing non-uniform distribution of p-values in this case, and it is, as far as I can tell, simply the result of the F-test incorrect assumption about the 30 data samples in there being independent), but with the whitening operation the distribution of p-values is now correct (i.e. uniform distribution under the null hypothesis) indicating that the degrees of freedom being used in these analyses (F(3,27), rather than the F(3,9) degrees of freedom that you would get for a multivariate version of this same analysis or with the inclusion of subject factors) are correct.
In any way, please let me know your thoughts, and/or if I am misinterpreting your question (and of course I would be more than happy to hear more / learn about why the two-stage approach in SPM might not provide the correct F-statistics or degrees of freedom in the context of a within-subject factor with more than two levels), and hope this helps
Alfonso
-----------------------------------
N = 10; % number of subjects
L = 3; % number of DVs/components
% prepares design and data
X = kron(eye(L),ones(N,1)); % design matrix
C = eye(L); % contrast of interest
A = randn(L); A=chol(eye(L)+A'*A); % cholesky factor (dependence between DVs)
W = kron(sqrtm(inv(A'*A)),eye(N)); % whitening matrix
X = W*X; % whitened design matrix
% runs simulations (second-level analyses with null data)
dofe = size(X,1)-rank(X); % denominator degrees of freedom
Nc0 = rank(X*C'); % numerator degrees of freedom
iX = pinv(X'*X);
ir = pinv(C*iX*C');
iXX = iX*X';
FF = [];
for niter=1:1e6, % number of simulations
Y = randn(N,L)*A; % simulated data
Y = W*reshape(Y,[],1); % data concatenated across DVs and whitened
B = iXX*Y; % model parameters
E = Y-X*B; % model error
EE = sum(abs(E).^2,1);
h = C*B;
BB = sum(h.*(ir*h),1);
F = real(BB./max(eps,EE))*dofe/Nc0; % F-statistic
FF(niter) = F;
end
% display null distribution of p-values (expected uniform distribution if valid analysis)
p=1-spm_Fcdf(FF,Nc0,dofe);
hist(p,100)
-------------------------------------------------
Originally posted by Martyn McFarquhar:
Hi Alfonso,
Thanks again for the detailed response, I promise the questions will stop soon!
Just to contextualise, we are still trying to resolve the discrepancy in results and my collaborator informs me that they didn't select the "non-parametric" option in CONN and as such I am assuming that the results they are reporting come from the univariate SPM approach. As such, I am just trying to understand this method fully so we can decide what to do.
I just want to check about the implementation of the "two-stage" approach as the guidance given in the attached paper (and by Will Penny on the SPM Wiki) do not actually give the correct F-statistics or degrees of freedom when you have a within-subject factor with > 2 levels. I have been in discussions with Guillaume Flandin at the FIL about this and I can send you scripts and examples to show this if you want.
For now I just want to clarify how the model works in CONN because in order to get an equivalent to the multivariate test across components (i.e. when the M contrast matrix is the identity) you would still need multiple components in the same model and thus would still require the inclusion of the subject factor in the design matrix in order to get the correct degrees of freedom and compute the correct error term for the F-tests. Is this being done in CONN?
Best wishes,
Martyn
Thanks again for the detailed response, I promise the questions will stop soon!
Just to contextualise, we are still trying to resolve the discrepancy in results and my collaborator informs me that they didn't select the "non-parametric" option in CONN and as such I am assuming that the results they are reporting come from the univariate SPM approach. As such, I am just trying to understand this method fully so we can decide what to do.
I just want to check about the implementation of the "two-stage" approach as the guidance given in the attached paper (and by Will Penny on the SPM Wiki) do not actually give the correct F-statistics or degrees of freedom when you have a within-subject factor with > 2 levels. I have been in discussions with Guillaume Flandin at the FIL about this and I can send you scripts and examples to show this if you want.
For now I just want to clarify how the model works in CONN because in order to get an equivalent to the multivariate test across components (i.e. when the M contrast matrix is the identity) you would still need multiple components in the same model and thus would still require the inclusion of the subject factor in the design matrix in order to get the correct degrees of freedom and compute the correct error term for the F-tests. Is this being done in CONN?
Best wishes,
Martyn
Oct 15, 2018 08:10 AM | Martyn McFarquhar
RE: Clarification on contrasts in CONN 2nd-level multivariate analysis
Hi Alfonso,
This is excellent work, thank you!
Indeed, the necessity of orthogonalising the contrast weights when the number of levels > 2 is not something I had previously realised. I'm intrigued by this and certainly has given me some food for thought. I do not believe that SPM orthogonalises contrast weights by default and so this is something I'm going to dig into more to try and gain some intuition.
Rather frustratingly, I notice that Penny and Henson do actually mention this in their SPM chapter, however it is relegated to a footnote and a short explanation in the appendices. I had missed this previously, hence the issue with implementing their approach. Perhaps more frustrating is the fact that Will Penny does not mention or implement this in the example on the SPM Wiki. Perhaps the assumption is that the non-sphericity correction will always be used? Nevertheless, this seems like such an important point that I've not seen made anywhere else previously.
Anyway, thank you for all the help and useful discussion. I hope you got something out of this as well and hasn't been seen as a waste of your time. Hopefully other people will find the discussion instructive as well.
Best wishes,
- Martyn
This is excellent work, thank you!
Indeed, the necessity of orthogonalising the contrast weights when the number of levels > 2 is not something I had previously realised. I'm intrigued by this and certainly has given me some food for thought. I do not believe that SPM orthogonalises contrast weights by default and so this is something I'm going to dig into more to try and gain some intuition.
Rather frustratingly, I notice that Penny and Henson do actually mention this in their SPM chapter, however it is relegated to a footnote and a short explanation in the appendices. I had missed this previously, hence the issue with implementing their approach. Perhaps more frustrating is the fact that Will Penny does not mention or implement this in the example on the SPM Wiki. Perhaps the assumption is that the non-sphericity correction will always be used? Nevertheless, this seems like such an important point that I've not seen made anywhere else previously.
Anyway, thank you for all the help and useful discussion. I hope you got something out of this as well and hasn't been seen as a waste of your time. Hopefully other people will find the discussion instructive as well.
Best wishes,
- Martyn
Oct 19, 2018 08:10 AM | Ali Amad
RE: Clarification on contrasts in CONN 2nd-level multivariate analysis
Many thanks for your answer and all
your fantastic work Alfonso.
Best wishes
Ali
Originally posted by Alfonso Nieto-Castanon:
Best wishes
Ali
Originally posted by Alfonso Nieto-Castanon:
Dear
Ali,
Sorry I had missed this message. Yes, I would definitely report those results and simply mention that the conclusions should be treated with caution until replicated given that the results do not survive a more conservative non-parametric randomization test.
Hope this helps
Alfonso
Originally posted by Ali Amad:
Sorry I had missed this message. Yes, I would definitely report those results and simply mention that the conclusions should be treated with caution until replicated given that the results do not survive a more conservative non-parametric randomization test.
Hope this helps
Alfonso
Originally posted by Ali Amad:
Dear
Alfonso,
first of all, thank you for help on this topic and for all the clarifications.
I have worked with Martyn's team on these analysis. Following MVPA analysis, we have indeed some significant results by using parametric statistics and nothing with non-parametric stats. However, by using the results from parametric stats as seeds to run seed to voxel analysis we found several interesting results.
Do you think that all these results are worthless and cannot be used at all ? Or can these results can be discussed with caution ?
Many thanks for all you help.
Best wishes,
Ali
Originally posted by Alfonso Nieto-Castanon:
first of all, thank you for help on this topic and for all the clarifications.
I have worked with Martyn's team on these analysis. Following MVPA analysis, we have indeed some significant results by using parametric statistics and nothing with non-parametric stats. However, by using the results from parametric stats as seeds to run seed to voxel analysis we found several interesting results.
Do you think that all these results are worthless and cannot be used at all ? Or can these results can be discussed with caution ?
Many thanks for all you help.
Best wishes,
Ali
Originally posted by Alfonso Nieto-Castanon:
Hi
Martyn,
Thank you for the clarification (and no need to stop the questions, I appreciate the interesting points your are bringing up)
Regarding the easy portion of your question, in this particular case (MVPA analyses) the best way to proceed, in my opinion, would be to use the multivariate (non-parametric) analysis version in CONN. Beyond the point regarding the correct degrees of freedom (more on this below) the main problem for univariate SPM analyses in the context of MVPA is that the components entered into your second-level analysis are being computed using an entirely different Principal Component decomposition for each voxel, which means that, even if for nearby voxels the subspace spanned by those components is expected to be similar, the identity of the components is not (i.e. the same component across two adjacent voxels cannot be guaranteed to represent the same or similar aspect of the pattern of connectivity with the rest of the brain), so the standard SPM assumptions in second-level analyses (including homogeneity of covariance across voxels as well as random field smoothness assumptions) are not all that reasonable in this case. Hence the recommendation here is to use multivariate analyses together with non-parametric cluster statistics (i.e simply selecting "non-parametric stats" in the results explorer window).
Now, coming back to the more difficult but perhaps also more interesting portion of your question, CONN does NOT include a subject factor explicitly in the second-level design matrix, even in the context of these repeated-measure analyses. I do not believe it would be appropriate to do so, for the following reasons (but please feel free to correct me and/or let me know your thoughts): 1) in the example analysis that you describe, subject effects (average MVPA component values for each component pre- and post- intervention) are already explicitly excluded by computing the "post - pre" differences in component scores, separately for each subject, and entering those differences directly (instead of entering separately both the pre- and post- scores) into the second-level analyses; 2) the same is not done across MVPA components (i.e. entering the difference in scores between components) because we are not interested in testing whether there are any differences between the 3 MVPA component scores, but rather whether there are any effects when pooling across the 3 components; and 3) as far as I understand the ReML pooled estimation of the covariance between the DVs, when the within-subject factor is modeled as having non-independent levels, is followed by a whitening step, which (I guess assuming that the covariance is not rank-defficient) should effectively remove the dependencies between the DVs (implicitly making the degrees of freedom of the F statistic used in these analyses correct).
To illustrate the latter point (and for me to double-check what I am saying here, since I could be mistaken) I am including below some test code that simulates random data for 10 subjects and 3 non-independent DVs (e.g. simulating the three "differences in MVPA scores" variables in our example analyses), and then it runs the "univariate (SPM)"-style second-level repeated-measures analysis without the subjects factor included in the design matrix (just like CONN does when using the "univariate (SPM)" option). The point here is to illustrate that, without the whitening operation (i.e. if you were to change "W=..." in the code below to "W=eye(L*N)", in order to skip whitening) you would be exactly right that the F- statistics computed there would be incorrect and have the incorrect degrees of freedom (this is shown by the code producing non-uniform distribution of p-values in this case, and it is, as far as I can tell, simply the result of the F-test incorrect assumption about the 30 data samples in there being independent), but with the whitening operation the distribution of p-values is now correct (i.e. uniform distribution under the null hypothesis) indicating that the degrees of freedom being used in these analyses (F(3,27), rather than the F(3,9) degrees of freedom that you would get for a multivariate version of this same analysis or with the inclusion of subject factors) are correct.
In any way, please let me know your thoughts, and/or if I am misinterpreting your question (and of course I would be more than happy to hear more / learn about why the two-stage approach in SPM might not provide the correct F-statistics or degrees of freedom in the context of a within-subject factor with more than two levels), and hope this helps
Alfonso
-----------------------------------
N = 10; % number of subjects
L = 3; % number of DVs/components
% prepares design and data
X = kron(eye(L),ones(N,1)); % design matrix
C = eye(L); % contrast of interest
A = randn(L); A=chol(eye(L)+A'*A); % cholesky factor (dependence between DVs)
W = kron(sqrtm(inv(A'*A)),eye(N)); % whitening matrix
X = W*X; % whitened design matrix
% runs simulations (second-level analyses with null data)
dofe = size(X,1)-rank(X); % denominator degrees of freedom
Nc0 = rank(X*C'); % numerator degrees of freedom
iX = pinv(X'*X);
ir = pinv(C*iX*C');
iXX = iX*X';
FF = [];
for niter=1:1e6, % number of simulations
Y = randn(N,L)*A; % simulated data
Y = W*reshape(Y,[],1); % data concatenated across DVs and whitened
B = iXX*Y; % model parameters
E = Y-X*B; % model error
EE = sum(abs(E).^2,1);
h = C*B;
BB = sum(h.*(ir*h),1);
F = real(BB./max(eps,EE))*dofe/Nc0; % F-statistic
FF(niter) = F;
end
% display null distribution of p-values (expected uniform distribution if valid analysis)
p=1-spm_Fcdf(FF,Nc0,dofe);
hist(p,100)
-------------------------------------------------
Originally posted by Martyn McFarquhar:
Thank you for the clarification (and no need to stop the questions, I appreciate the interesting points your are bringing up)
Regarding the easy portion of your question, in this particular case (MVPA analyses) the best way to proceed, in my opinion, would be to use the multivariate (non-parametric) analysis version in CONN. Beyond the point regarding the correct degrees of freedom (more on this below) the main problem for univariate SPM analyses in the context of MVPA is that the components entered into your second-level analysis are being computed using an entirely different Principal Component decomposition for each voxel, which means that, even if for nearby voxels the subspace spanned by those components is expected to be similar, the identity of the components is not (i.e. the same component across two adjacent voxels cannot be guaranteed to represent the same or similar aspect of the pattern of connectivity with the rest of the brain), so the standard SPM assumptions in second-level analyses (including homogeneity of covariance across voxels as well as random field smoothness assumptions) are not all that reasonable in this case. Hence the recommendation here is to use multivariate analyses together with non-parametric cluster statistics (i.e simply selecting "non-parametric stats" in the results explorer window).
Now, coming back to the more difficult but perhaps also more interesting portion of your question, CONN does NOT include a subject factor explicitly in the second-level design matrix, even in the context of these repeated-measure analyses. I do not believe it would be appropriate to do so, for the following reasons (but please feel free to correct me and/or let me know your thoughts): 1) in the example analysis that you describe, subject effects (average MVPA component values for each component pre- and post- intervention) are already explicitly excluded by computing the "post - pre" differences in component scores, separately for each subject, and entering those differences directly (instead of entering separately both the pre- and post- scores) into the second-level analyses; 2) the same is not done across MVPA components (i.e. entering the difference in scores between components) because we are not interested in testing whether there are any differences between the 3 MVPA component scores, but rather whether there are any effects when pooling across the 3 components; and 3) as far as I understand the ReML pooled estimation of the covariance between the DVs, when the within-subject factor is modeled as having non-independent levels, is followed by a whitening step, which (I guess assuming that the covariance is not rank-defficient) should effectively remove the dependencies between the DVs (implicitly making the degrees of freedom of the F statistic used in these analyses correct).
To illustrate the latter point (and for me to double-check what I am saying here, since I could be mistaken) I am including below some test code that simulates random data for 10 subjects and 3 non-independent DVs (e.g. simulating the three "differences in MVPA scores" variables in our example analyses), and then it runs the "univariate (SPM)"-style second-level repeated-measures analysis without the subjects factor included in the design matrix (just like CONN does when using the "univariate (SPM)" option). The point here is to illustrate that, without the whitening operation (i.e. if you were to change "W=..." in the code below to "W=eye(L*N)", in order to skip whitening) you would be exactly right that the F- statistics computed there would be incorrect and have the incorrect degrees of freedom (this is shown by the code producing non-uniform distribution of p-values in this case, and it is, as far as I can tell, simply the result of the F-test incorrect assumption about the 30 data samples in there being independent), but with the whitening operation the distribution of p-values is now correct (i.e. uniform distribution under the null hypothesis) indicating that the degrees of freedom being used in these analyses (F(3,27), rather than the F(3,9) degrees of freedom that you would get for a multivariate version of this same analysis or with the inclusion of subject factors) are correct.
In any way, please let me know your thoughts, and/or if I am misinterpreting your question (and of course I would be more than happy to hear more / learn about why the two-stage approach in SPM might not provide the correct F-statistics or degrees of freedom in the context of a within-subject factor with more than two levels), and hope this helps
Alfonso
-----------------------------------
N = 10; % number of subjects
L = 3; % number of DVs/components
% prepares design and data
X = kron(eye(L),ones(N,1)); % design matrix
C = eye(L); % contrast of interest
A = randn(L); A=chol(eye(L)+A'*A); % cholesky factor (dependence between DVs)
W = kron(sqrtm(inv(A'*A)),eye(N)); % whitening matrix
X = W*X; % whitened design matrix
% runs simulations (second-level analyses with null data)
dofe = size(X,1)-rank(X); % denominator degrees of freedom
Nc0 = rank(X*C'); % numerator degrees of freedom
iX = pinv(X'*X);
ir = pinv(C*iX*C');
iXX = iX*X';
FF = [];
for niter=1:1e6, % number of simulations
Y = randn(N,L)*A; % simulated data
Y = W*reshape(Y,[],1); % data concatenated across DVs and whitened
B = iXX*Y; % model parameters
E = Y-X*B; % model error
EE = sum(abs(E).^2,1);
h = C*B;
BB = sum(h.*(ir*h),1);
F = real(BB./max(eps,EE))*dofe/Nc0; % F-statistic
FF(niter) = F;
end
% display null distribution of p-values (expected uniform distribution if valid analysis)
p=1-spm_Fcdf(FF,Nc0,dofe);
hist(p,100)
-------------------------------------------------
Originally posted by Martyn McFarquhar:
Hi Alfonso,
Thanks again for the detailed response, I promise the questions will stop soon!
Just to contextualise, we are still trying to resolve the discrepancy in results and my collaborator informs me that they didn't select the "non-parametric" option in CONN and as such I am assuming that the results they are reporting come from the univariate SPM approach. As such, I am just trying to understand this method fully so we can decide what to do.
I just want to check about the implementation of the "two-stage" approach as the guidance given in the attached paper (and by Will Penny on the SPM Wiki) do not actually give the correct F-statistics or degrees of freedom when you have a within-subject factor with > 2 levels. I have been in discussions with Guillaume Flandin at the FIL about this and I can send you scripts and examples to show this if you want.
For now I just want to clarify how the model works in CONN because in order to get an equivalent to the multivariate test across components (i.e. when the M contrast matrix is the identity) you would still need multiple components in the same model and thus would still require the inclusion of the subject factor in the design matrix in order to get the correct degrees of freedom and compute the correct error term for the F-tests. Is this being done in CONN?
Best wishes,
Martyn
Thanks again for the detailed response, I promise the questions will stop soon!
Just to contextualise, we are still trying to resolve the discrepancy in results and my collaborator informs me that they didn't select the "non-parametric" option in CONN and as such I am assuming that the results they are reporting come from the univariate SPM approach. As such, I am just trying to understand this method fully so we can decide what to do.
I just want to check about the implementation of the "two-stage" approach as the guidance given in the attached paper (and by Will Penny on the SPM Wiki) do not actually give the correct F-statistics or degrees of freedom when you have a within-subject factor with > 2 levels. I have been in discussions with Guillaume Flandin at the FIL about this and I can send you scripts and examples to show this if you want.
For now I just want to clarify how the model works in CONN because in order to get an equivalent to the multivariate test across components (i.e. when the M contrast matrix is the identity) you would still need multiple components in the same model and thus would still require the inclusion of the subject factor in the design matrix in order to get the correct degrees of freedom and compute the correct error term for the F-tests. Is this being done in CONN?
Best wishes,
Martyn