glm_joinffx.m

toolbox of BVQX, This is the access between BV and matlab. It will help you to analysis data from BV

function hfile3 = glm_JoinFFX(hfile, hfile2)
% GLM::JoinFFX  - join two fixed effects GLMs
%
% FORMAT:       combined = glm1.JoinFFX(glm2);
%
% Input fields:
%
%       glm2        second FFX GLM to be added to glm1
%
% Joins the Fixed Effects GLM results of glm1 and glm2 to one
% single structure.

% Version:  v0.7b
% Build:    7083022
% Date:     Aug-30 2007, 10:26 PM CEST
% Author:   Jochen Weber, Brain Innovation, B.V., Maastricht, NL
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% argument check
if nargin < 2 || ...
    numel(hfile) ~= 1 || ...
    numel(hfile2) ~= 1 || ...
   ~isBVQXfile(hfile, 'glm') || ...
   ~isBVQXfile(hfile2, 'glm')
    error( ...
        'BVQXfile:BadArgument', ...
        'Invalid call to %s.', ...
        mfilename ...
    );
end
bc1 = bvqxfile_getcont(hfile.L);
bc2 = bvqxfile_getcont(hfile2.L);
if bc1.ProjectTypeRFX ~= 0 || ...
    bc2.ProjectTypeRFX ~= 0 || ...
    bc1 .SeparatePredictors < 1
    error( ...
        'BVQXfile:BadArgument', ...
        'Invalid call to %s.', ...
        mfilename ...
    );
end
% check type and sizes of GLMs
if isfield(bc1.GLMData, 'RFXGlobalMap') || ...
    isfield(bc2.GLMData, 'RFXGlobalMap') || ...
    bc1.ProjectType ~= bc2.ProjectType || ...
    bc1.FileVersion ~= bc2.FileVersion || ...
    bc1.SeparatePredictors ~= bc2.SeparatePredictors || ...
    bc1.TransformationType ~= bc2.TransformationType || ...
    bc1.Resolution ~= bc2.Resolution || ...
    bc1.XStart ~= bc2.XStart || ...
    bc1.XEnd ~= bc2.XEnd || ...
    bc1.YStart ~= bc2.YStart || ...
    bc1.YEnd ~= bc2.YEnd || ...
    bc1.ZStart ~= bc2.ZStart || ...
    bc1.ZEnd ~= bc2.ZEnd
    error( ...
        'BVQXfile:BadObject', ...
        'Invalid object(s) given in call. Crucial fields mismatch.' ...
    );
end

% build new file
DM1 = bc1.DesignMatrix;
DM2 = bc2.DesignMatrix;
NP1 = bc1.NrOfPredictors;
NP2 = bc2.NrOfPredictors;
NPS = NP1 + NP2;
NS1 = bc1.NrOfStudies;
NS2 = bc2.NrOfStudies;
SPS = bc1.SeparatePredictors;
NSS = NS1 + NS2;
NTP = bc1.NrOfTimePoints + bc2.NrOfTimePoints;
VA1 = (bc1.GLMData.MCorrSS > 0);
VA2 = (bc2.GLMData.MCorrSS > 0);
VAC = (VA1 & VA2);
hfile3 = aft_CopyObject(hfile);
sc3 = bvqxfile_getscont(hfile3.L);
sc3.F = '';
bc3 = sc3.C;
bc3.NrOfTimePoints = NTP;
bc3.NrOfPredictors = NPS;
bc3.NrOfStudies = NSS;
bc3.NrOfVoxelsForBonfCorrection = sum(VAC(:));
bc3.Study(1:end) = [];
bc3.Predictor(1:end) = [];
bc3.DesignMatrix = [];
bc3.iXX = [];

% concatenate design matrices
try
    bc3.DesignMatrix = [ ...
        DM1(:, 1:end-NS1), ...
        zeros(size(DM1, 1), size(DM2, 2) - NS2), ...
        DM1(:, end+1-NS1:end), ...
        zeros(size(DM1, 1), NS2); ...
        ...
        zeros(size(DM2, 1), size(DM1, 2) - NS1), ...
        DM2(:, 1:end-NS2), ...
        zeros(size(DM2, 1), NS1), ...
        DM2(:, end+1-NS2:end)];
    bc3.iXX = pinv(bc3.DesignMatrix' * bc3.DesignMatrix);
catch
    error( ...
        'BVQXfile:BadObjects', ...
        'Error pseudo-inverting concatenated design matrix.' ...
    );
end

% copy studys
bc3.Study = [bc1.Study(:)', bc2.Study(:)'];

% copy predictors
tpc = 1;
for pc = 1:(NP1 - NS1)
    bc3.Predictor(tpc) = bc1.Predictor(pc);
    tpc = tpc + 1;
end
for pc = 1:(NP2 - NS2)
    bc3.Predictor(tpc) = bc2.Predictor(pc);
    bc3.Predictor(tpc).Name1 = sprintf('Predictor: %d', tpc);
    if SPS == 1
        prname = bc3.Predictor(tpc).Name2;
        [pm_m{1:3}] = regexpi(prname, '^study\s+(\d+)\:\s+(.*)$');
        if isempty(pm_m{1}) || ...
            isempty(pm_m{3}) || ...
            ~all(size(pm_m{3}{1}) == 2)
            error( ...
                'BVQXfile:InternalError', ...
                'Error extracting study predictor name.' ...
            );
        end
        pm_m = pm_m{3};
        bc3.Predictor(tpc).Name2 = sprintf('Study %d: %s', ...
            str2double(prname(pm_m{1}(1, 1):pm_m{1}(1, 2))) + NS1, ...
            prname(pm_m{1}(2, 1):pm_m{1}(2, 2)));
    end
    tpc = tpc + 1;
end
for pc = (NP1 + 1 - NS1):NP1
    bc3.Predictor(tpc) = bc1.Predictor(pc);
    bc3.Predictor(tpc).Name1 = sprintf('Predictor: %d', tpc);
    tpc = tpc + 1;
end
for pc = (NP2 + 1 - NS2):NP2
    bc3.Predictor(tpc) = bc2.Predictor(pc);
    bc3.Predictor(tpc).Name1 = sprintf('Predictor: %d', tpc);
    tpc = tpc + 1;
end

% create GLMData fields
bc3.GLMData.MultipleRegressionR = sqrt( ...
    (NS1 * bc1.GLMData.MultipleRegressionR .^ 2+ ...
     NS2 * bc1.GLMData.MultipleRegressionR .^ 2) / (NSS));
bc3.GLMData.MultipleRegressionR(~VAC) = 0;
bc3.GLMData.MCorrSS(~VAC) = 0;
bc3.GLMData.MCorrSS(VAC) = NTP;
maporder = [ ...
     1:1:(NP1 - NS1), ...
    -1:-1:-(NP2 - NS2), ...
     (NP1 + 1 - NS1):1:NP1, ...
    -(NP2 + 1 - NS2):-1:-NP2];
maps1 = find(maporder > 0);
maps2 = find(maporder < 0);
if bc1.ProjectType ~= 2
    bc3.GLMData.BetaMaps(:, :, :, maps1) =  ...
        bc1.GLMData.BetaMaps(:, :, :, maporder(maps1));
    bc3.GLMData.BetaMaps(:, :, :, maps2) = ...
        bc2.GLMData.BetaMaps(:, :, :, -maporder(maps2));
    bc3.GLMData.XY(:, :, :, maps1) = ...
        bc1.GLMData.XY(:, :, :, maporder(maps1));
    bc3.GLMData.XY(:, :, :, maps2) = ...
        bc2.GLMData.XY(:, :, :, -maporder(maps2));
else
    bc3.GLMData.BetaMaps(:, maps1) = ...
        bc1.GLMData.BetaMaps(:, maporder(maps1));
    bc3.GLMData.BetaMaps(:, maps2) = ...
        bc2.GLMData.BetaMaps(:, -maporder(maps2));
    bc3.GLMData.XY(:, maps1) = ...
        bc1.GLMData.XY(:, maporder(maps1));
    bc3.GLMData.XY(:, maps2) = ...
        bc2.GLMData.XY(:, -maporder(maps2));
end
betasz = size(bc3.GLMData.BetaMaps);
numvox = prod(betasz) / NPS;
bc3.GLMData.BetaMaps = reshape(bc3.GLMData.BetaMaps, [numvox, NPS]);
bc3.GLMData.BetaMaps(~VAC(:), :) = 0;
bc3.GLMData.BetaMaps = reshape(bc3.GLMData.BetaMaps, betasz);
bc3.GLMData.XY = reshape(bc3.GLMData.XY, [numvox, NPS]);
bc3.GLMData.XY(~VAC(:), :) = 0;
bc3.GLMData.XY = reshape(bc3.GLMData.XY, betasz);
bc3.GLMData.TimeCourseMean = ...
    (NS1 * bc1.GLMData.TimeCourseMean + NS2 * bc2.GLMData.TimeCourseMean) / NSS;
bc3.GLMData.TimeCourseMean(~VAC) = 0;

% clean up fields
ct = eps;
bc3.DesignMatrix(abs(bc3.DesignMatrix) < ct) = 0;
bc3.iXX(abs(bc3.iXX) < ct) = 0;
bc3.GLMData.MultipleRegressionR(abs(bc3.GLMData.MultipleRegressionR) < ct) = 0;
bc3.GLMData.MCorrSS(abs(bc3.GLMData.MCorrSS) < ct) = 0;
bc3.GLMData.BetaMaps(abs(bc3.GLMData.BetaMaps) < ct) = 0;
bc3.GLMData.XY(abs(bc3.GLMData.XY) < ct) = 0;
bc3.GLMData.TimeCourseMean(abs(bc3.GLMData.TimeCourseMean) < ct) = 0;

% put back
sc3.C = bc3;
bvqxfile_setscont(hfile3.L, sc3);