Hi Christian,<div><br></div><div>Thank you - you have been very helpful. I have a few more questions regarding the PCA based variance normalization especially in regards to multi subject analysis.</div><div><br></div><div>
(1) When you are doing concat ICA, does the PCA based variance normalization operate separately for each subject? (my guess is yes). If yes, then this will require latent dimensionality estimation for each subject. If so, is the average latent dimensionality over single subjects used for concatenated data latent dimension?</div>
<div><br></div><div>(2) What is the sequence of steps for concat ICA? My guess is as follows:</div><div><br></div><div>(a) For each subject, estimate latent dimension and hat_sigma_i for variance normalization.</div><div>
<br></div><div>(b) Variance normalize data for each subject.</div><div><br></div><div>(c) Concatenate variance normalized data for all subjects and choose latent dim = mean latent dim over single subjects.</div><div><br></div>
<div>(d) Proceed as if analyzing single subject data.</div><div><br></div><div>Is this accurate? If not, I am very curious to know the true sequence of steps.</div><div><br></div><div>Thanks,</div><div>Hans.</div><div><br>
</div><div><br><div class="gmail_quote">On Thu, Aug 20, 2009 at 4:44 PM, Christian F. Beckmann <span dir="ltr">&lt;<a href="mailto:[log in to unmask]"><a href="/cgi-bin/wa-jisc.exe?LOGON=A3%3Dind0908%26L%3DFSL%26E%3Dquoted-printable%26P%3D1630085%26B%3D--0016e644b9e0a47c4f0471aa060a%26T%3Dtext%252Fhtml%3B%2520charset%3Dwindows-1252%26pending%3D" target="_parent" >[log in to unmask]</a></a>&gt;</span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex;">
<div style="word-wrap:break-word">Hi Hans,<div><br></div><div>Sorry to not have been more clear in the section of the paper you quote. In section IV A we go into more details and illustrate the entire methodlogy (including the re-estimation of \Sigma_i from the noise estimate) in fig 4. Specifically (page 146) says</div>
<div><i><br></i></div><div><div><i>We can re-estimate \Sigma_I</i></div><div><i>from the residuals and iterate the entire cycle. In practice, </i></div><div><i>the output results do not suggest a strong dependency on the </i></div>
<div><i>form of \Sigma_i and preliminary results suggest that it is sufficient </i></div><div><i>to iterate these steps only once.</i></div><div><br></div><div>That is, we do not re-estimate again (which one could do in theory) as the estimate does not seem to improve after the main problem (different tissue types) has been addressed. </div>
<div><br></div><div>The section you refer to relates to the very first (initial) variance normalisation _before_ calculating the first PCA basis (that&#39;s why at that stage it is described as pre-processing). </div><div>
This is similar to e.g. one deals with modeling temporal autocorrelation: ignore it first, then get residuals and estimate AR parameters from these....</div><div>In practice, the initial variance normalisation (based on x_i only)  becomes irrelevant once if you re-estimate the noise variance, so melodic does not do this initial step and simply does a PCA on the de-meaned data in order to get the first set of residuals.</div>
<div><br></div><div>Note, however, that if one wasn&#39;t to iterate at all (i.e. only does the simple vn) then I still wouldn&#39;t call this &#39;wrong&#39; - the estimator is not optimal and not very robust but correct under the null....tissue-type induced differences are by far the largest in almost all cases (unless your data is real rubbish), so in practice even the simple variance normalisation (deviding by the standard deviation of x_i only) does a pretty good job - in fact, the improvements I see in simulations when doing the more complicated variance normalisation are small compared to the improvements we already see with the simple variance normalisation...  it&#39;s only &#39;wrong&#39; in the sense that such an approach has reduced sensitivity...</div>
<div><br></div><div>hth</div><div>Christian</div><div><div></div><div class="h5"><div><br></div><div><br></div><div><br></div><div><div>On 20 Aug 2009, at 19:20, Hans Tissot wrote:</div><br><blockquote type="cite">Hi Christian,<div>
<br></div><div>Thanks for your explanation and for sending that figure. I fully agree with you that noise standard deviation is in fact not isotropic and so we need to process the data so that our model assumptions are more closely satisfied. It is also clear that there are two approaches to variance normalization:</div>
 <div><br></div><div>(A) Dividing each timecourse by its standard deviation </div><div>(perhaps used in previous MELODIC versions and also suggested in your IEEE paper)</div><div><br></div><div>(B) Division of each timecourse by hat_sigma_i (estimated after an initial PCA). </div>
 <div>(perhaps the newer MELODIC uses this approach)</div><div><br></div><div>Let x_i = demeaned time courses.</div><div><br></div><div>(1) </div><div>Yes, under the null the two approaches are identical since if the time course is pure Gaussian noise then the non-Gaussian components s_i will be near 0 and [ x_i - mu - A * s_i will be approximately equal to x_i - mu ] . Hence std(x_i - mu) will give the correct standard deviation hat_sigma_i. </div>
 <div><br></div><div>However, if the signal contains *any* non-Gaussian component (s_i is non-zero) then std(x_i - mu) will NOT give hat_sigma_i. In this case we have to resort to std( x_i - mu - A * s_i ) to get hat_sigma_i where the product A * s_i can be computed by an initial PCA step only.<br>
</div><div><br></div><div><br></div><div>(2) </div><div>You call the newer approach to get hat_sigma_i (based on x_i - mu - A * s_i) to be a &quot;robust estimate&quot;. I have to disagree with you here. I feel that this newer approach is the &quot;correct&quot; estimate and the former based on std(x_i - mu) and ignoring A*s_i is &quot;wrong&quot;. </div>
 <div><br></div><div>The older approach will not achieve equalization of voxelwise standard deviation, in fact it might introduce unwanted differences across voxels based on s_i!! Since your IEEE paper talks about using std(x_i) approach for variance normalization, my guess is that previous versions of MELODIC did in fact use this approach. </div>
 <div><br></div><div>This deeply concerns me - since this would mean that previous MELODIC versions produced invalid results? Please let me know your thoughts.</div><div><br></div><div>Thanks,</div><div>Hans.</div><div><br>
<div>On Thu, Aug 20, 2009 at 1:13 PM, Christian F. Beckmann &lt;<a href="mailto:[log in to unmask]" target="_blank"><a href="/cgi-bin/wa-jisc.exe?LOGON=A3%3Dind0908%26L%3DFSL%26E%3Dquoted-printable%26P%3D1630085%26B%3D--0016e644b9e0a47c4f0471aa060a%26T%3Dtext%252Fhtml%3B%2520charset%3Dwindows-1252%26pending%3D" target="_parent" >[log in to unmask]</a></a>&gt; wrote:<br><div>Hi Hans,<div><br></div><div>Under the null assumption these two approaches are identical, the only difference being that in the case of the more complicated variance normalisation the variance is based on an estimate from a truncated Gaussian (i.e. you normalise to e.g. 0.95 globally rather than 1 globally but this is irrelavant given the later Z-transform). Under the alternative hypothesis of there being signal somewhere, the PCA based variance normalisation is basically a robust estimate for the the noise standard deviation. </div>
 </div><div> <br></div><div><div>Given the parameterisation we can flexibly move between both extremes (no variance normalisation to full variance normalisation). The details are not in the paper, the overall concept of having to do such a normalising step is.</div>
 <div><br></div><div>Given typical sparseness and low SNR of BOLD signals of interest the difference between the two cases is largely negligible - where it really matters is wrt non-BOLD signal variations, e.g. slice dropouts or spiking, which can induce changes in variance _much_ above the level you get from signals due to activation.</div>
 <div>The reason why it matters is that in such extreme cases the noise standard deviation does not correspond to an isotropic case which then messes with the mixture model fit used for inference in that voxels which do not contain any signal (the background noise class) cannot simply be modelled using a single Gaussian (as the corresponding PDF for the background noise will be a mixture of Gaussians just as in the case where no variance normalisation is carried out)  - it is for this reason that we now use a robust estimator of the voxel-wise variances in the software.</div>
 <div>The thresholding simply is used temporarily to get a robust estimate of the noise which is then applied to the non-thresholded data. </div><div>Practically, apart from the extreme cases of scanner-indiced variation, this normalisation deals with the very simple problem of tissue-dependent noise variance. Below is a plot of this (on log-scale!) highlighting that if no variance normalisation is carried out then the a single CSF voxel can easily weigh as much as 10 GM voxels in the calculation of the data covariance - one reason for not doing what you suggest (i.e. only use this normalisation in the dim-estimation step)</div>
 <div><br></div><div>hth</div><div>Christian</div><div><br></div><div><br></div><div></div></div><br><div><div></div><div><br></div><div><br></div><div><br><div><div>On 20 Aug 2009, at 17:02, Hans Tissot wrote:</div><br><blockquote type="cite">
Hi Christian,<div><br></div><div>Sorry for jumping into this thread but I am very confused. I have 2 questions:</div><div><br></div> <div>(1) I found your description in the previous e-mail different from that described in your paper. Here&#39;s the text from your IEEE paper:</div>
 <div><div><br></div><div>&quot;An immediate consequence of the fact that we are operating</div> <div>under an isotropic noise model is that as an initial preprocessing</div> <div>step we will modify the original data time courses to be normalized</div>
 <div>to zero mean and unit variance. This appears to be a sensible</div> <div>step in that on the one hand we know that the voxel-wise standard</div> <div>deviation of resting state data varies significantly over the</div>
 <div>brain but on the other hand, all voxels’ time courses are assumed</div> <div>to be generated from the same noise process. This variance-normalization</div> <div>preconditions the data under the “null hypotheses” of</div>
 <div>purely Gaussian noise, i.e., in the absence of any signal: the data</div> <div>matrix X is identical up to second-order statistics to a simple set</div> <div>of realizations from a N(0,1) noise process. Any signal component</div>
 <div>contained in X will have to reveal itself via its deviation</div> <div>from Gaussianity. This will turn out to be of prime importance</div> <div>both for the estimation of the number of sources and the final</div> <div>
inferential steps.</div> <div><br></div><div>After a voxel-wise normalization of variance, two voxels with</div><div>comparable noise level that are modulated by the same signal</div><div>time course, aj, say, but by different amounts will have the same</div>
<div>regression coefficient upon regression against aj. The difference</div><div>in the original amount of modulation is, therefore, contained</div><div>in the standard deviation of the residual noise. Forming</div><div>voxel-wise statistics, i.e., dividing the PICA maps by the estimated</div>
<div>standard deviation of eta, thus, is invariant under the initial</div><div>variance-normalization.&quot;</div> </div><div><br></div><div>Based on this description it seems that each time course is indeed normalized to unit variance (and this is a critical step as per the description above). Could you please clarify.</div>
 <div><br></div><div>(2) </div><div>A different approach that you described in the previous e-mail seems to divide each time course x_i by hat_sigma_i to equalize noise variance across voxels as assumed in the model. I could not find anything in the paper about thresholding the loadings at p .05 and doing the division by hat_sigma_i only on those voxels which pass the threshold. If I missed it, could you please direct me to the correct section? Why is this thresholding required anyways, why not just divide all voxels by their sigma_hat_i, once the latent dimensionality is estimated?</div>
 <div><br></div><div><br></div><div>Thanks,</div><div>Hans Tissot.</div><div><br><div>On Wed, Aug 19, 2009 at 4:37 PM, Christian F. Beckmann &lt;<a href="mailto:[log in to unmask]" target="_blank"><a href="/cgi-bin/wa-jisc.exe?LOGON=A3%3Dind0908%26L%3DFSL%26E%3Dquoted-printable%26P%3D1630085%26B%3D--0016e644b9e0a47c4f0471aa060a%26T%3Dtext%252Fhtml%3B%2520charset%3Dwindows-1252%26pending%3D" target="_parent" >[log in to unmask]</a></a>&gt; wrote:<br>
<div>Hi<div><br></div><div>Neither - Melodic variance normalisation attempts to normalise the residual _noise_ and therefore scales by the inverse of a robust estimate of the noise, rather than normalise by signal + noise (which, as you say, would remove important amplitude information). More specifically, melodic runs a full PCA on the data, thresholds loadings at a p .05 level (uncorrected) in order to indentify signals which potentially will later on pass thresholding. Melodic then regresses this out of the data and uses the residuals from this regression to identify voxel-wise standard deviations. Each voxel then gets normalised by the inverse of this estimate. This then renders the data closer to the model assumptions of the probabilistic ICA model (which assumes isotropic noise) and improves signal detection and gives a much improved mixture model fit for inference - see Beckmann and Smith (2004) IEEE TMI for details.</div>
 <div>hth</div><div>Christian</div><div><br></div><div><br></div><div><div><div></div><div><br><div><div>On 19 Aug 2009, at 21:17, MCLAREN, Donald wrote:</div><br><blockquote type="cite"><div>I&#39;ve been trying to figure out how the variance normalization works in FSL. </div>
 <div> </div> <div>Is the variance of each voxel set to 1 or is the overall variance of all voxels set to 1?</div> <div> </div> <div>In the context of resting fMRI, if the former is true, isn&#39;t information about the amplitude of the flucuations being removed? This would seem to be important if your comparing data from different scans were done in two different states (e.g. eyes open versus eyes closed) since the amplitude of the flucuations how been shown to be different in these states.</div>
 <div> </div> <div>Any thoughts would be appreciated.<br clear="all"></div> <div></div><br>-- <br>Best Regards, Donald McLaren<br>=====================<br>Support the Alzheimer&#39;s Association Memory Walk 2009 ~~~<br>Join the Wisconsin Alzheimer&#39;s Disease Research Center team or make a donation.<br>
<a href="http://madison.kintera.org/2009/donaldmclaren" target="_blank">http://madison.kintera.org/2009/donaldmclaren</a><br><br><br>=================<br>D.G. McLaren<br>University of Wisconsin - Madison<br>Neuroscience Training Program<br>
 Office: (608) 265-9672<br>Lab: (608) 256-1901 ext 12914<br>=====================<br>This e-mail contains CONFIDENTIAL INFORMATION which may contain PROTECTED HEALTHCARE INFORMATION and may also be LEGALLY PRIVILEGED and which is intended only for the use of the individual or entity named above. If the reader of the e-mail is not the intended recipient or the employee or agent responsible for delivering it to the intended recipient, you are hereby notified that you are in possession of confidential and privileged information. Any unauthorized use, disclosure, copying or the taking of any action in reliance on the contents of this information is strictly prohibited and may be unlawful. If you have received this e-mail unintentionally, please immediately notify the sender via telephone at (608) 265-9672 or email.<br>
</blockquote></div><br></div></div><div><div><div>_______________________________________________<br>Christian F. Beckmann, DPhil<br>Senior Lecturer, Clinical Neuroscience Department<br>Division of Neuroscience and Mental Health<br>
Imperial College London, Hammersmith Campus<br> Rm 419, Burlington Danes Bldg, Du Cane Road, London W12 0NN, UK<br>Tel.: +44 (0)20 7594 6685   ---   Fax: +44 (0)20 7594 6548<br>Email: <a href="mailto:[log in to unmask]" target="_blank"><a href="/cgi-bin/wa-jisc.exe?LOGON=A3%3Dind0908%26L%3DFSL%26E%3Dquoted-printable%26P%3D1630085%26B%3D--0016e644b9e0a47c4f0471aa060a%26T%3Dtext%252Fhtml%3B%2520charset%3Dwindows-1252%26pending%3D" target="_parent" >[log in to unmask]</a></a><br>
<a href="http://www.imperial.ac.uk/medicine/people/c.beckmann/" target="_blank">http://www.imperial.ac.uk/medicine/people/c.beckmann/</a><br><br>Senior Research Fellow, FMRIB Centre<br>University of Oxford<br>JR Hospital - Oxford OX3 9DU<br>
 Tel.: +44 (0)1865 222551 --- Fax: +44 (0)1865 222717<br>Email: <a href="mailto:[log in to unmask]" target="_blank"><a href="/cgi-bin/wa-jisc.exe?LOGON=A3%3Dind0908%26L%3DFSL%26E%3Dquoted-printable%26P%3D1630085%26B%3D--0016e644b9e0a47c4f0471aa060a%26T%3Dtext%252Fhtml%3B%2520charset%3Dwindows-1252%26pending%3D" target="_parent" >[log in to unmask]</a></a><br><a href="http://www.fmrib.ox.ac.uk/~beckmann" target="_blank">http://www.fmrib.ox.ac.uk/~beckmann</a><br>
<br><br></div></div> </div><br></div></div></div><br></div></blockquote></div><br><div><div><div>_______________________________________________<br>Christian F. Beckmann, DPhil<br>Senior Lecturer, Clinical Neuroscience Department<br>
Division of Neuroscience and Mental Health<br>Imperial College London, Hammersmith Campus<br> Rm 419, Burlington Danes Bldg, Du Cane Road, London W12 0NN, UK<br>Tel.: +44 (0)20 7594 6685   ---   Fax: +44 (0)20 7594 6548<br>
Email: <a href="mailto:[log in to unmask]" target="_blank"><a href="/cgi-bin/wa-jisc.exe?LOGON=A3%3Dind0908%26L%3DFSL%26E%3Dquoted-printable%26P%3D1630085%26B%3D--0016e644b9e0a47c4f0471aa060a%26T%3Dtext%252Fhtml%3B%2520charset%3Dwindows-1252%26pending%3D" target="_parent" >[log in to unmask]</a></a><br><a href="http://www.imperial.ac.uk/medicine/people/c.beckmann/" target="_blank">http://www.imperial.ac.uk/medicine/people/c.beckmann/</a><br>
<br>Senior Research Fellow, FMRIB Centre<br>University of Oxford<br>JR Hospital - Oxford OX3 9DU<br> Tel.: +44 (0)1865 222551 --- Fax: +44 (0)1865 222717<br>Email: <a href="mailto:[log in to unmask]" target="_blank"><a href="/cgi-bin/wa-jisc.exe?LOGON=A3%3Dind0908%26L%3DFSL%26E%3Dquoted-printable%26P%3D1630085%26B%3D--0016e644b9e0a47c4f0471aa060a%26T%3Dtext%252Fhtml%3B%2520charset%3Dwindows-1252%26pending%3D" target="_parent" >[log in to unmask]</a></a><br>
<a href="http://www.fmrib.ox.ac.uk/~beckmann" target="_blank">http://www.fmrib.ox.ac.uk/~beckmann</a><br><br><br></div></div> </div><br></div></div><br></div><br></div></blockquote></div><br><div><div><div>_______________________________________________<br>
Christian F. Beckmann, DPhil<br>Senior Lecturer, Clinical Neuroscience Department<br>Division of Neuroscience and Mental Health<br>Imperial College London, Hammersmith Campus<br>Rm 419, Burlington Danes Bldg, Du Cane Road, London W12 0NN, UK<br>
Tel.: +44 (0)20 7594 6685   ---   Fax: +44 (0)20 7594 6548<br>Email: <a href="mailto:[log in to unmask]" target="_blank"><a href="/cgi-bin/wa-jisc.exe?LOGON=A3%3Dind0908%26L%3DFSL%26E%3Dquoted-printable%26P%3D1630085%26B%3D--0016e644b9e0a47c4f0471aa060a%26T%3Dtext%252Fhtml%3B%2520charset%3Dwindows-1252%26pending%3D" target="_parent" >[log in to unmask]</a></a><br><a href="http://www.imperial.ac.uk/medicine/people/c.beckmann/" target="_blank">http://www.imperial.ac.uk/medicine/people/c.beckmann/</a><br>
<br>Senior Research Fellow, FMRIB Centre<br>University of Oxford<br>JR Hospital - Oxford OX3 9DU<br>Tel.: +44 (0)1865 222551 --- Fax: +44 (0)1865 222717<br>Email: <a href="mailto:[log in to unmask]" target="_blank"><a href="/cgi-bin/wa-jisc.exe?LOGON=A3%3Dind0908%26L%3DFSL%26E%3Dquoted-printable%26P%3D1630085%26B%3D--0016e644b9e0a47c4f0471aa060a%26T%3Dtext%252Fhtml%3B%2520charset%3Dwindows-1252%26pending%3D" target="_parent" >[log in to unmask]</a></a><br>
<a href="http://www.fmrib.ox.ac.uk/~beckmann" target="_blank">http://www.fmrib.ox.ac.uk/~beckmann</a><br><br><br></div></div> </div><br></div></div></div><div><div></div><div class="h5"><br><br><div> <span style="border-collapse:separate;color:rgb(0, 0, 0);font-family:Helvetica;font-size:12px;font-style:normal;font-variant:normal;font-weight:normal;letter-spacing:normal;line-height:normal;text-align:auto;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px"><div style="word-wrap:break-word">
<div>_______________________________________________<br>Christian F. Beckmann, DPhil<br>Senior Lecturer, Clinical Neuroscience Department<br>Division of Neuroscience and Mental Health<br>Imperial College London, Hammersmith Campus<br>
Rm 419, Burlington Danes Bldg, Du Cane Road, London W12 0NN, UK<br>Tel.: +44 (0)20 7594 6685   ---   Fax: +44 (0)20 7594 6548<br>Email: <a href="mailto:[log in to unmask]" target="_blank"><a href="/cgi-bin/wa-jisc.exe?LOGON=A3%3Dind0908%26L%3DFSL%26E%3Dquoted-printable%26P%3D1630085%26B%3D--0016e644b9e0a47c4f0471aa060a%26T%3Dtext%252Fhtml%3B%2520charset%3Dwindows-1252%26pending%3D" target="_parent" >[log in to unmask]</a></a><br>
<a href="http://www.imperial.ac.uk/medicine/people/c.beckmann/" target="_blank">http://www.imperial.ac.uk/medicine/people/c.beckmann/</a><br><br>Senior Research Fellow, FMRIB Centre<br>University of Oxford<br>JR Hospital - Oxford OX3 9DU<br>
Tel.: +44 (0)1865 222551 --- Fax: +44 (0)1865 222717<br>Email: <a href="mailto:[log in to unmask]" target="_blank"><a href="/cgi-bin/wa-jisc.exe?LOGON=A3%3Dind0908%26L%3DFSL%26E%3Dquoted-printable%26P%3D1630085%26B%3D--0016e644b9e0a47c4f0471aa060a%26T%3Dtext%252Fhtml%3B%2520charset%3Dwindows-1252%26pending%3D" target="_parent" >[log in to unmask]</a></a><br><a href="http://www.fmrib.ox.ac.uk/~beckmann" target="_blank">http://www.fmrib.ox.ac.uk/~beckmann</a><br>
<br><br></div></div></span> </div><br></div></div></div></blockquote></div><br></div>