clear;tic;
% Moderated Regression For A 2-way Interaction Between Continuous Variables



% Load or create column vectors of data with the following names: idv, mod, dv


%%%%%%%%%%%%%%%%%%%%%%%%%  START OF TRIAL-RUN DATA  %%%%%%%%%%%%%%%%%%%%%%%%%
%  The following commands generate artificial data that can be used for
%  a trial-run of the program.  Just run this whole file. 
idv = randn(200,1);
mod = randn(200,1);
e   = rand(200,1);
dv  = 8.7*idv +  -1.49*mod +  1.23*(idv .* mod)  + e*30;
%%%%%%%%%%%%%%%%%%%%%%%%%  END OF TRIAL-RUN DATA  %%%%%%%%%%%%%%%%%%%%%%%%%


% Specify the number of chunks for the IDV
chunkIDV = 7 ;

% Specify the number of chunks for the MOD
chunkMOD = 7 ;

% Specify the method of determining highest & lowest values for the IDV & MOD
% for an optimal design
% Enter 1 to let the program do it automatically
% Enter 2 to use your own preferred values
optdes = 1;

% If you entered 2 on the above "optdes =" statement, then enter your preferred
% values on the next four statements:

% Enter the lowest value for the IDV in an optimal design
lowIDV = 1;

% Enter the highest value for the IDV in an optimal design
highIDV = 5;

% Enter the lowest value for the MOD in an optimal design
lowMOD = 1;

% Enter the highest value for the MOD in an optimal design
highMOD = 5;

% Specify the # of randomized data sets for the randomization
% test of statistical significance.  Recommendations: 
% use 100 data sets for a trial run, but use 1000 or more for final results
permutes = 100;

% End of required user specifications


disp([' ']);disp([' ']);disp([' ']);disp([' ']);
disp(['Moderated Regression For A 2-way Interaction Between Continuous Variables:']);

% sample size
bigN = size(idv,1);

% centering the IDV & MOD
idv = idv - mean(idv);
mod = mod - mean(mod);

% 2-D frequency distribution.
IDVmax = max(idv) + .00000001;
IDVmin = min(idv) - .00000001;
MODmax = max(mod) + .00000001;
MODmin = min(mod) - .00000001;
increm1 = (IDVmax - IDVmin) / chunkIDV;
increm2 = (MODmax - MODmin) / chunkMOD;
freqs = repmat(0, chunkIDV, chunkMOD);
for luper = 1 : bigN;
IDVchunk =  ceil( ( idv(luper,1) - IDVmin ) / increm1 );
MODchunk =  ceil( ( mod(luper,1) - MODmin ) / increm2 );
freqs( IDVchunk, MODchunk ) = freqs( IDVchunk, MODchunk ) + 1;
end

% Moderated Regression
x = idv .* mod;
cr = corrcoef([ idv mod x dv ]);
beta = inv(cr(1:3,1:3)) * cr(1:3,4);
sd = std( [ idv mod x dv ] );
mn = mean( [ idv mod x dv ] );
b = (sd(1,4) ./ sd(1,1:3))   .*  beta';
a = mn(1,4) - ( sum( (mn(1,1:3) .* b ),2) );
r2all  = beta' * cr(1:3,4);
r2main =  ( inv(cr(1:2,1:2))*cr(1:2,4))'  *  cr(1:2,4);
r2chXn = r2all - r2main;
F = (r2all-r2main) / ((1-r2all)/(bigN-3-1));
dferror = bigN - 3 - 1;
pF = 1 - fcdf(F,1,dferror);

% Moderated Regression including quadratic terms
idvq = idv .* idv;
modq = mod .* mod;
crq = corrcoef([ idv mod idvq modq x dv ]);
betaq = inv(crq(1:5,1:5)) * crq(1:5,6);
sdq = std( [ idv mod idvq modq x dv ] );
mnq = mean([ idv mod idvq modq x dv ] );
bq = (sdq(1,6) ./ sdq(1,1:5))   .*  betaq';
aq = mnq(1,6) - ( sum( (mnq(1,1:5) .* bq ),2) );
r2allq  = betaq' * crq(1:5,6);
r2mainq =  ( inv(crq(1:4,1:4))*crq(1:4,6))'  *  crq(1:4,6);
r2chXnq = r2allq - r2mainq;
Fq = (r2allq-r2mainq) / ((1-r2allq)/(bigN-5-1));
dferrorq = bigN - 5 - 1;
pFq = 1 - fcdf(Fq,1,dferror);

% f-squared = the proportion of systematic variance accounted for by the effect 
%             relative to unexplained variance in the criterion (A & W, 1991, p. 157)
fsquare  = (r2all - r2main) / (1 - r2all);
fsquareq = (r2allq - r2mainq) / (1 - r2allq);

% mse  Darlington p 121
xx = [ ones(bigN,1) idv mod x  ];
sse = dv' * dv  -  [ a; b']' * xx' * dv;
mse = sse / ( bigN - 3 - 1);

% Vxzxz = the residual variance of the product XZ after controlling for X & Z (p. 378);
%         2 ways of computing Vxzxz
%         method 1 = the MSE that results from regressing XZ on X & Z (p. 381)
%         method 2 = using the full mod reg equation, = 
%         MSE / (n * stand. error of the estimate for the product.^2) (p. 381)

% method 2: the full mod reg equation approach (yields slightly diff results in unusual sits)
% Cohen & Cohen, 2003, p. 632
Rijm1 = inv(cr(1:3,1:3)); 
SExn = ( sd(1,4) ./ sd(1,3) ) * sqrt( (1 - r2all) / ( bigN - 3 - 1) ) * sqrt( Rijm1(3,3) ); 
Vxzxz2 = mse / ( ( bigN - 3 - 1) * SExn.^2 );

% method 1: using the MSE that results from regressing XZ on X & Z approach
cr2 = corrcoef([ idv mod x ]);
beta2 = inv(cr2(1:2,1:2)) * cr2(1:2,3);
sd2 = std( [ idv mod x ] );
mn2 = mean( [ idv mod x ] );
b2 = (sd2(1,3) ./ sd2(1,1:2))   .*  beta2';
a2 = mn2(1,3) - ( sum( (mn2(1,1:2) .* b2 ),2) );
xx2 = [ ones(bigN,1) idv mod ];
sse2 = x' * x  -  [ a2; b2']' * xx2' * x;
% mse  Darlington p 121
mse2 = sse2 / ( bigN - 2 -1);
Vxzxz1 = mse2;


% PRE, & Relative efficiency / optimal design values

if optdes == 1;
% Using the mean scores in the top and bottoms chunks, based on the specified chunk sizes
idvtop = -9999;
idvbot = -9999;
modtop = -9999;
modbot = -9999;
for luper = 1 : size(idv,1)
if  idv(luper,1) > (IDVmax - increm1); idvtop = [ idvtop ; idv(luper,1) ];end
if  idv(luper,1) < (IDVmin + increm1); idvbot = [ idvbot ; idv(luper,1) ];end
if  mod(luper,1) > (MODmax - increm2); modtop = [ modtop ; mod(luper,1) ];end
if  mod(luper,1) < (MODmin + increm2); modbot = [ modbot ; mod(luper,1) ];end
end
highIDV  = mean(idvtop(2:size(idvtop,1),1));
lowIDV   = mean(idvbot(2:size(idvbot,1),1));
highMOD  = mean(modtop(2:size(modtop,1),1));
lowMOD   = mean(modbot(2:size(modbot,1),1));
end

% maximim value of Vxzxz (p. 381, 383)
maxVxzxz = (( (highIDV - lowIDV) / 2 ) .^2 )  *  (( (highMOD - lowMOD) / 2 ) .^2 );

% Relative Efficiency 
releffic = Vxzxz1 / maxVxzxz;

% PRE squared partial correlation = proportional reduction in error (p 377, 384)
%   = the model improvement due to adding the product term
PRE = 1 / ( 1+ ( mse / (b(1,3).^2 * Vxzxz1) ) );

% PRE for an optimal design p 384
PRE2 = 1 / ( 1+ ( mse / (b(1,3).^2 * maxVxzxz) ) );


disp([' ']);
disp(['Sample Size:                                      ' num2str(bigN) ]);disp([' ']);
disp(['IDV Mean:                                         ' num2str(mean(idv)) ]);
disp(['IDV Standard Deviation:                           ' num2str(std(idv)) ]);
disp(['IDV Lowest Score:                                 ' num2str(IDVmin) ]);
disp(['IDV Highest Score:                                ' num2str(IDVmax) ]);disp([' ']);
disp(['MOD Mean:                                         ' num2str(mean(mod)) ]);
disp(['MOD Standard Deviation:                           ' num2str(std(mod)) ]);
disp(['MOD Lowest Score:                                 ' num2str(MODmin) ]);
disp(['MOD Highest Score:                                ' num2str(MODmax) ]);disp([' ']);
disp(['Specified # of chunks for the IDV:                ' num2str(chunkIDV) ]);
disp(['Specified # of chunks for the MOD:                ' num2str(chunkMOD) ]);disp([' ']);
disp(['IDV Low  Value for an Optimal Design:             ' num2str(lowIDV) ]);
disp(['IDV High Value for an Optimal Design:             ' num2str(highIDV) ]);disp([' ']);
disp(['MOD Low  Value for an Optimal Design:             ' num2str(lowMOD) ]);
disp(['MOD High Value for an Optimal Design:             ' num2str(highMOD) ]);disp([' ']);
disp([' ']);disp(['Joint Frequencies of the IDV & MOD:   ']);disp([' ']);disp(freqs);
bar3(freqs)
set(gca,'FontName','Helvetica', 'FontSize',14, 'fontweight', 'normal' );
xlabel('Predictor');ylabel('Moderator');zlabel('Frequency');
title ('Joint Distribution of the Predictor & Moderator Variables:');disp([' ']);
disp(['Vxzxz1 (residual variance of the product term):    ' num2str(Vxzxz1) ]);disp([' ']);
disp(['Vxzxz2 (residual variance of the product term):    ' num2str(Vxzxz2) ]);disp([' ']);
disp(['maxVxzxz (maximum possible value of Vxzxz):        ' num2str(maxVxzxz) ]);disp([' ']);
disp(['Relative Efficiency (Vxzxz / maxVxzxz):            ' num2str(releffic) ]);disp([' ']);
disp(['PRE (Proportional Reduction in Error):             ' num2str(PRE) ]);disp([' ']);
disp(['PRE for the optimal design:                        ' num2str(PRE2) ]);disp([' ']);disp([' ']);
disp(['Regression results for the equation that DOES NOT include quadratic terms:']);disp([' ']);
disp(['Slope coefficient for the IDV:                     ' num2str(b(1,1)) ]);disp([' ']);
disp(['Slope coefficient for the MOD:                     ' num2str(b(1,2)) ]);disp([' ']);
disp(['Slope coefficient for the product term:            ' num2str(b(1,3)) ]);disp([' ']);
disp(['Intercept:                                         ' num2str(a) ]);disp([' ']);
disp(['R-squared change for the interaction term:          ' num2str(r2chXn) ]);disp([' ']);
disp(['f-squared effect size for the interaction term:     ' num2str(fsquare) ]);disp([' ']);
disp(['F value for the interaction term:                  ' num2str(F) ]);disp([' ']);
disp(['Significance level of F (tabled value):            ' num2str(pF) ]);disp([' ']);
disp(['Specified number of randomized data sets:          ' num2str(permutes) ]);disp([' ']);


% Analyses of the randomized data
counter  = 0;
counterq = 0;
for nperms = 1:permutes

%  The raw data permutations are based on column-wise random shufflings
%  of the values in the raw data matrix using Castellan's (1992, 
%  BRMIC, 24, 72-77) algorithm; The distributions of the original 
%  raw variables are exactly preserved in the shuffled versions used
%  in the parallel analyses.

% data matrix to permute
permdata = [ idv mod dv ];
for luper = 1:(bigN -1)
k1 = fix( (bigN - luper + 1) * rand(1) + 1 )  + luper - 1;
k2 = fix( (bigN - luper + 1) * rand(1) + 1 )  + luper - 1;
k3 = fix( (bigN - luper + 1) * rand(1) + 1 )  + luper - 1;
d1 = permdata(luper,1);
d2 = permdata(luper,2);
d3 = permdata(luper,3);
permdata(luper,1) = permdata(k1,1);
permdata(luper,2) = permdata(k2,2);
permdata(luper,3) = permdata(k3,3);
permdata(k1,1) = d1;
permdata(k2,2) = d2;
permdata(k3,3) = d3;
end

idv2 = permdata(:,1);
mod2 = permdata(:,2);
dv2  = permdata(:,3);

% Moderated Regression
x2 = idv2 .* mod2;
cr2 = corrcoef([ idv2 mod2 x2 dv2 ]);
beta2 = inv(cr2(1:3,1:3)) * cr2(1:3,4);
r2all2  = beta2' * cr2(1:3,4);
r2main2 =  ( inv(cr2(1:2,1:2))*cr2(1:2,4))' * cr2(1:2,4);
F2 = (r2all2-r2main2) / ((1-r2all2)/(bigN-3-1));
if F2 >= F;  counter = counter + 1;end

% Moderated Regression including quadratic terms
idvq2 = idv2 .* idv2;
modq2 = mod2 .* mod2;
crq2 = corrcoef([ idv2 mod2 idvq2 modq2 x2 dv2 ]);
betaq2 = inv(crq2(1:5,1:5)) * crq2(1:5,6);
r2allq2  = betaq2' * crq2(1:5,6);
r2mainq2 =  ( inv(crq2(1:4,1:4))*crq2(1:4,6))'  *  crq2(1:4,6);
Fq2 = (r2allq2-r2mainq2) / ((1-r2allq2)/(bigN-5-1));
if Fq2 >= Fq;  counterq = counterq + 1;end

end

% significance level computation from Noreen (1989, p. 56)
siglevel  = (counter  + 1) / (permutes + 1);
siglevelq = (counterq + 1) / (permutes + 1);

disp(['Significance level of F (randomization test):      ' num2str(siglevel) ]);disp([' ']);disp([' ']);

disp(['Regression results for the equation that INCLUDES quadratic terms:']);disp([' ']);
disp(['Slope coefficient for the IDV:                     ' num2str(bq(1,1)) ]);disp([' ']);
disp(['Slope coefficient for the MOD:                     ' num2str(bq(1,2)) ]);disp([' ']);
disp(['Slope coefficient for the IDV quadratic term:      ' num2str(bq(1,3)) ]);disp([' ']);
disp(['Slope coefficient for the MOD quadratic term:      ' num2str(bq(1,4)) ]);disp([' ']);
disp(['Slope coefficient for the product term:            ' num2str(bq(1,5)) ]);disp([' ']);
disp(['Intercept:                                         ' num2str(aq) ]);disp([' ']);
disp(['R-squared change for the interaction term:         ' num2str(r2chXnq) ]);disp([' ']);
disp(['f-squared effect size for the interaction term:    ' num2str(fsquareq) ]);disp([' ']);
disp(['F value for the interaction term:                  ' num2str(Fq) ]);disp([' ']);
disp(['Significance level of F (tabled value):            ' num2str(pFq) ]);disp([' ']);
disp(['Specified number of randomized data sets:          ' num2str(permutes) ]);disp([' ']);
disp(['Significance level of F (randomization test):      ' num2str(siglevelq) ]);disp([' ']);disp([' ']);

disp(['Time for this problem (in seconds):                ', num2str(toc) ]); disp([' '])

