Matlab simulation demonstrating how squared distance matrices can be linearly combined - probably requires my pilab and matlab-plotting repos at least. Generates a nice figure at the end.

test_multirsa_again.m

% demonstration of how squared distances can be combined linearly in
% multiple regression RSA, while ordinary distances cannot. In this case
% the distortions from doing multiple regression are relatively subtle but
% perhaps the difference could be made more dramatic by other means.
%
% Actually, the less subtle effect is on parameter estimates - the fits
% from model_multi2 below are [.5 .5], which is absolutely correct since
% the predictors were already scaled to match the data. By contrast the
% non-squared fit comes up with [.41 .65] - regardless of how minor the
% distortions are in the overall fitted RDM, the parameter estimates are
% wrong!
%
function test_multirsa_again

% let's start with a true grid representation. We will attempt to use
% multiple regression to reconstruct the full XY matrix from one regressor
% representing X and one Y.
p = [1:5]';
[x,y] = meshgrid(p,p);
xy = [x(:) y(:)];
% just to complicate things, scale and stretch the distances a bit.
xy(:,1) = xy(:,1) .^2;
xy(:,2) = xy(:,2) * 3;

% finally rotate the matrix (this shouldn't actually matter)
% xy_rot = rotatepoints(xy',1/3*pi)';

% now our predictors in RDM space are
% full distances
rdm_euc = pdist(xy)';
% distances in X
rdm_x = pdist(xy(:,1))';
% distances in Y
rdm_y = pdist(xy(:,2))';

% do multidimensional scaling on the predicted (fitted) RDM from each
mdwrap = @(thismodel)mdscale(asrdmmat(predictY(thismodel)),2,...
    'criterion','metricstress');
mdwraproot = @(thismodel)mdscale(asrdmmat(sqrt(predictY(thismodel))),2,...
    'criterion','metricstress');

plotfields.original = xy;

% fit to the original - just to confirm that this does recover the RDM
model_org = RSA(rdm_euc,rdm_euc);
bfields.original = [1,1];

% ordinary multiple regression - tends to produce weird bendy lines in the
% MDscale representation
model_multi = RSA([rdm_x rdm_y],rdm_euc);
plotfields.fit_multi = mdwrap(model_multi);
rfields.fit_multi = corr(predictY(model_multi),rdm_euc);
bfields.fit_multi = fit(model_multi);

% squared distance multiple regression - after taking the square root of
% the fitted RDM we recover the original RDM perfectly.
model_multi2 = RSA([rdm_x rdm_y].^2,rdm_euc.^2);
plotfields.fit_squared = mdwraproot(model_multi2);
rfields.fit_squared = corr(sqrt(predictY(model_multi2)),rdm_euc);
bfields.fit_squared = fit(model_multi2);

% plot all these
plotgroups = fieldnames(plotfields);
nc = 2;
nr = numel(plotgroups);
f = figure(111);
clf(f);
ax = [];
for p = 1:numel(plotgroups)
    ax(end+1) = subplot(nr,nc,numel(ax)+1);
    bar(bfields.(plotgroups{p}));
    set(gca,'xtick',[1:2],'xticklabel',{'x','y'});
    xlabel('regressor');
    if any(strfind(plotgroups{p},'original'))
        ylabel('ground truth');
    else
        ylabel('parameter estimate');
    end
    box off;
    ax(end+1) = subplot(nr,nc,numel(ax)+1);
    plot(plotfields.(plotgroups{p})(:,1),...
        plotfields.(plotgroups{p})(:,2),'o');
    set(ax(end),'dataaspectratio',[1 1 1]);
    tstr = stripbadcharacters(plotgroups{p},' ');
    if p>1
        tstr = [tstr sprintf(' r=%.2f',rfields.(plotgroups{p}))];
    end
    axis off;
    titlebetter(tstr,ax(end-1:end),false);
end
[fundir,funname,ext] = fileparts(mfilename('fullpath'));
printstandard(fullfile(fundir,funname));
keyboard;