I have run a mass-univariate spatiotemporal model using 'lme_mass_fit_EMinit' and then applied an FDR correction across both hemispheres (syntax provided below). On completion I realised that most of the corrected p-values for my main and interaction effects were well above the alpha of 0.05. The corrected p-value thresholds in log base 10 format are included in the attached figure and summarised below.

Sample

Our sample includes 112 controls subjects, and 106 at-risk subjects (individuals with a first degree relative  of BD but not personal history of BD). All subjects are aged between 12 and 30 years. 153 or these subjects have 2 time-points (77 control, 76 at-risk), while the remaining 65 individuals have only 1 MRI time-point. We also have mixed ethnicites - 165 Caucasians, 23 Asians, and 30 mixed (Asians-Caucasians).

Analysis Design

The QDEC file contains the following 6 variables:

1. Y (years between scans)
2. A (baseline age)
3. G (group, 1=case, 0=control)
4. S (sex, 1=female, 0=male)
5. E1 (Ethnicity 1, 1=asian, 0=other)
6. E2 (Ethnicity 1, 1=mixed-asian-caucasian, 0=other)

Design matrix: [ones(length(M),1) M M(:,1).*M(:,3)]
                       i.e. main effects for each of the qdec variables + an interaction term for years X group

DVs:     Cortical thickness, area and volume

Model: Mass-univariate spatiotemporal model using 'lme_mass_fit_EMinit'        

                CODE:

                % Read in surface files

                [Y,mri] = fs_read_Y(mgh);

                % Read in qdec file

                Qdec = fReadQdec(qdec);
                % Remove fsid from qdec
                Qdec = rmQdecCol(Qdec,1);
                % Store col 1 (fsid-base) in sID variable
                sID = Qdec(2:end,1);
                % Remove col 1 (fsid-base) from Qdec array
                Qdec = rmQdecCol(Qdec,1);
                % Convert Qdec to numeric matrix M
                M = Qdec2num(Qdec);

                % Sort data and evaluate design matrix

                [M,Y,ni] = sortData(M,1,Y,sID);
                X = eval([ones(length(M),1) M M(:,1).*M(:,3)])

                % Compute vertex-wise temporal covariance estimates.

                [Th0, Re] = lme_mass_fit_EMinit(X,[1],Y,ni,cortex,3);

                %Segmentation and model fitting.

                [Rgs, RgMeans, stats] = fit(Th0, Re, [1], sphere, cortex, X, Y, ni);  

                %Check surfaces.

                surfcomp(Th0, RgMeans, sphere, fig1, fig2)         

Correction for multiple comparisons: FDR across both hemispheres

                CODE:

                 P = [ F_lhstats.pval(lhcortex) F_rhstats.pval(rhcortex) ];
                 G = [ F_lhstats.sgn(lhcortex) F_rhstats.sgn(rhcortex) ];
                 [detvtx, sided_pval, pth] = lme_mass_FDR2(P,G,[],0.05,0);
                 pcor = -log10(pth);
                 thrlh(r,2) = {pcor};
                 thrrh(r,2) = {pcor};
                 [~,dc] = size(detvtx);
                 dvtx(r,2) = {dc};