%
% Preparation
%
%%% 1) Create a structure array specifying details of data
ACEfit_Par.P_nm = 'HCP_imgs.txt'; % A list of images or other data
% specification; see FileFormats.txt
ACEfit_Par.InfMx = 'KinInf.csv'; % Kinship information matrix of 4
% columns with headers
% ('SubjectID','MotherID',
% 'FatherID','Zygosity'). MotherID
% and FatherID must be numeric,
% and Zygosity is one of 'MZ',
% 'NotMZ' and 'NotTwin'. The order
% of subjects must match the order
% of image paths or subjects in
% ACEfit_Par.P_nm.
ACEfit_Par.ResDir = '/my/path/ResDir';
%%% The rest are optional; omit to use default values
ACEfit_Par.Pmask = 'HCP_mask.nii'; % Brain mask image (default: whole volume)
ACEfit_Par.Dsnmtx = ''; % Design matrix (default: all-ones vector)
ACEfit_Par.Nlz = 1; % Inverse Gaussian normalisation options
% (applied to each phenotype voxel/element)
% 0 - None
% 1 - Gaussian normalisation *before*
% forming residuals
% 2 - Gaussian normalisation *after*
% forming residuals; note, though,
% that after this Gaussian
% normalisation the residuals
% will be formed again to ensure
% orthognality of residuals to
% the design matrix.
ACEfit_Par.AggNlz = 0; % Aggregate heritability normalisation options
% (applied to each phenotype voxel/element)
% 0 - Default, meaning that only
% determined by Nlz option.
% This always entails
% de-meaning and removal of any effects
% in Dsnmtx, but variance at each
% voxel/element unchanged.
% 1 - Same as 0, but with variance
% normalisation; this is full
% 'studentization'.
% 2 - *Undo* mean centering of
% residual formation. (If Dsnmtx
% is empty or default value of
% an all-ones vector, this is
% equivalent to using the raw
% input data. If Dsnmtx
% contains nuisance regressors
% the residuals will be formed
% and the mean added back in.
% 3 - Same as 2, but with
% variance normalisation;
% note that for each
% voxel/element, mean is
% first added back and then
% data are divided by stdev.
ACEfit_Par.ContSel = []; % Select a single contrast (a volume
% in a 4D Nifti file supplied for
% each subject, or at the last
% dimension in a cifti image; NOT
% compatibile with a single file
% containing all subjects' data.)
ACEfit_Par.NoImg = 0; % If 1, suppress image-wise inference,
% and only compute summaries.
%%% 2) Updata 'ACEfit_Par' with the input data information
ACEfit_Par = PrepData(ACEfit_Par);
%%% 3) Run the original data once
ACEfit_Par.alpha_CFT = []; % Cluster-forming threshold (default: 0.05)
ACEfit_Par = ACEfit(ACEfit_Par);
%%% 4) Add permutation and bootstrapping information, and save "ACEfit_Par.mat"
ACEfit_Par.nPerm = 1000; % Number of permutations
ACEfit_Par.nBoot = 1000; % Number of bootstrap replicates
nParallel = []; % Number of parallel runs (default: 1, without parallelization)
PrepParallel(ACEfit_Par,nParallel);
%
% Permutations
%
if ACEfit_Par.nPerm>0
%%% 1)
%%%%% The following code can be scripted or parallelized as you wish.
%%%%% Please refer to "README_APACE_intro.pdf" for the example snippets
%%%%% for parallelization.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load(fullfile(ACEfit_Par.ResDir,'ACEfit_Par.mat'));
RunID = 1;
ACEfit_Perm_Parallel(ACEfit_Par,RunID);
% Inside ResDir, it will create result sets, ACEfit_Parallel_XXXX.mat,
% one for each RunID.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% 2) This will merge together all the results from the specified RunID's
load(fullfile(ACEfit_Par.ResDir,'ACEfit_Par.mat'));
ACEfit_Perm_Parallel_Results(ACEfit_Par);
%%% 3)
FWEalpha = 0.05;
FDRalpha = 0.05;
load(fullfile(ACEfit_Par.ResDir,'ACEfit_Par.mat'));
ACEfit_Results(ACEfit_Par,FWEalpha,FDRalpha);
end
%
% Bootstrapping
%
if ACEfit_Par.nBoot>0
%%% 1)
%%%%% The following code can be scripted or parallelized as you wish.
%%%%% Please refer to "README_APACE_intro.pdf" for the example snippets
%%%%% for parallelization.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load(fullfile(ACEfit_Par.ResDir,'ACEfit_Par.mat'));
RunID = 1;
ACEfit_Boot_Parallel(ACEfit_Par,RunID);
% Inside ResDir, it will create result sets, BootCI_Parallel_XXXX.mat,
% one for each RunID.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% 2) This will merge together all the results from the specified RunID's
load(fullfile(ACEfit_Par.ResDir,'ACEfit_Par.mat'));
ACEfit_Boot_Parallel_Results(ACEfit_Par);
%%% 3)
Balpha = 0.05;
load(fullfile(ACEfit_Par.ResDir,'ACEfit_Par.mat'));
Boot_CIs(ACEfit_Par,Balpha)
end
%
% Aggregate heritability (aka "Steve's method"), with P-values via
% permutation and CI's via boostrapping.
%
% Note that permuation (steps 1-3) and bootstrapping (steps 1-3) can be
% skipped by setting ACEfit_Par.nPerm=0 and ACEfit_Par.nBoot=0 separately;
% the following code can be run immediately after "PrepParallel".
%
load(fullfile(ACEfit_Par.ResDir,'ACEfit_Par.mat'));
Palpha = 0.05; % Significance threhold for permutations (plotting only)
Balpha = 0.05; % Confidence level, where CI's have level 100*(1-alpha)
AgHe_Method(ACEfit_Par,Palpha,Balpha)
% % Once with no variance normalisation
% ACEfit_Par.AggNlz = 0; % de-meaning only
% AgHe_Method(ACEfit_Par,Palpha,Balpha,'_NoNorm')
%
% % Now, again, but with Variance normalisation
% ACEfit_Par.AggNlz = 1; % de-meaning and scaling to have stdev of 1.0
% AgHe_Method(ACEfit_Par,Palpha,Balpha,'_Norm')
%
% Generate summary file
%
APACEsummary(ACEfit_Par,'ResultSummary')