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

%%%%%%%%%%%%%%%%%%%%%%%%%  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);
mod1 = randn(200,1);
mod2 = randn(200,1);
e    = rand(200,1);
dv   = 8.7*idv +  -1.49*mod1 +  2.2*mod2 + 1.23*idv.*mod1 + 3.1*idv.*mod2 ...
       + -1.3*idv.*mod1.*mod2  + e*30;
%%%%%%%%%%%%%%%%%%%%%%%%%  END OF TRIAL-RUN DATA  %%%%%%%%%%%%%%%%%%%%%%%%%


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

% Specify the number of chunks for MOD1
chunkMOD1 = 5 ;

% Specify the number of chunks for MOD1
chunkMOD2 = 5 ;

% 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 MOD1 in an optimal design
lowMOD1  = 1;

% Enter the highest value for MOD1 in an optimal design
highMOD1 = 5;

% Enter the lowest value for MOD1 in an optimal design
lowMOD2  = 1;

% Enter the highest value for MOD1 in an optimal design
highMOD2 = 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 3-way Interaction Between Continuous Variables:']);

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

% centering the IDV & MOD1 & MOD2
idv  = idv  - mean(idv);
mod1 = mod1 - mean(mod1);
mod2 = mod2 - mean(mod2);

% 3-D min, max, increm
IDVmax  = max(idv)  + .00000001;
IDVmin  = min(idv)  - .00000001;
MOD1max = max(mod1) + .00000001;
MOD1min = min(mod1) - .00000001;
MOD2max = max(mod2) + .00000001;
MOD2min = min(mod2) - .00000001;
increm1 = (IDVmax - IDVmin) / chunkIDV;
increm2 = (MOD1max - MOD1min) / chunkMOD1;
increm3 = (MOD2max - MOD2min) / chunkMOD2;

% Moderated Regression
x1 = idv  .* mod1;
x2 = idv  .* mod2;
x3 = mod1 .* mod2;
x4 = idv .* mod1 .* mod2;
cr = corrcoef([ idv mod1 mod2 x1 x2 x3 x4 dv ]);
beta = inv(cr(1:7,1:7)) * cr(1:7,8);
sd = std( [ idv mod1 mod2 x1 x2 x3 x4 dv ] );
mn = mean( [ idv mod1 mod2 x1 x2 x3 x4 dv ] );
b = (sd(1,8) ./ sd(1,1:7))   .*  beta';
a = mn(1,8) - ( sum( (mn(1,1:7) .* b ),2) );
r2all  = beta' * cr(1:7,8);
r2main =  ( inv(cr(1:6,1:6))*cr(1:6,8))'  *  cr(1:6,8);
r2chXn = r2all - r2main;
F = (r2all-r2main) / ((1-r2all)/(bigN-7-1));
dferror = bigN - 7 - 1;
pF = 1 - fcdf(F,1,dferror);

% Moderated Regression including quadratic terms
idvq = idv .* idv;
mod1q = mod1 .* mod1;
mod2q = mod2 .* mod2;
crq = corrcoef([ idv mod1 mod2 idvq mod1q mod2q x1 x2 x3 x4 dv ]);
betaq = inv(crq(1:10,1:10)) * crq(1:10,11);
sdq = std( [ idv mod1 mod2 idvq mod1q mod2q x1 x2 x3 x4 dv ] );
mnq = mean([ idv mod1 mod2 idvq mod1q mod2q x1 x2 x3 x4 dv ] );
bq = (sdq(1,11) ./ sdq(1,1:10))   .*  betaq';
aq = mnq(1,11) - ( sum( (mnq(1,1:10) .* bq ),2) );
r2allq  = betaq' * crq(1:10,11);
r2mainq =  ( inv(crq(1:9,1:9))*crq(1:9,11))'  *  crq(1:9,11);
r2chXnq = r2allq - r2mainq;
Fq = (r2allq-r2mainq) / ((1-r2allq)/(bigN-10-1));
dferrorq = bigN - 10 - 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 mod1 mod2 x1 x2 x3 x4  ];
sse = dv' * dv  -  [ a; b']' * xx' * dv;
mse = sse / ( bigN - 7 - 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:7,1:7)); 
SExn = ( sd(1,8) ./ sd(1,7) ) * sqrt( (1 - r2all) / ( bigN - 7 - 1) ) * sqrt( Rijm1(7,7) ); 
Vxzxz2 = mse / ( ( bigN - 7 - 1) * SExn.^2 );

% method 1: using the MSE that results from regressing XZ on X & Z approach
cr2 = corrcoef([ idv mod1 mod2 x1 x2 x3 x4 ]);
beta2 = inv(cr2(1:6,1:6)) * cr2(1:6,7);
sd2 = std(  [ idv mod1 mod2 x1 x2 x3 x4 ] );
mn2 = mean( [ idv mod1 mod2 x1 x2 x3 x4 ] );
b2 = (sd2(1,7) ./ sd2(1,1:6))   .*  beta2';
a2 = mn2(1,7) - ( sum( (mn2(1,1:6) .* b2 ),2) );
xx2 = [ ones(bigN,1) idv mod1 mod2 x1 x2 x3 ];
sse2 = x4' * x4  -  [ a2; b2']' * xx2' * x4;
% mse  Darlington p 121
mse2 = sse2 / ( bigN - 6 -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;
mod1top = -9999;
mod1bot = -9999;
mod2top = -9999;
mod2bot = -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  mod1(luper,1) > (MOD1max - increm2); mod1top = [ mod1top ; mod1(luper,1) ];end
if  mod1(luper,1) < (MOD1min + increm2); mod1bot = [ mod1bot ; mod1(luper,1) ];end
if  mod2(luper,1) > (MOD2max - increm3); mod2top = [ mod2top ; mod2(luper,1) ];end
if  mod2(luper,1) < (MOD2min + increm3); mod2bot = [ mod2bot ; mod2(luper,1) ];end
end
highIDV   = mean(idvtop(2:size(idvtop,1),1));
lowIDV    = mean(idvbot(2:size(idvbot,1),1));
highMOD1  = mean(mod1top(2:size(mod1top,1),1));
lowMOD1   = mean(mod1bot(2:size(mod1bot,1),1));
highMOD2  = mean(mod2top(2:size(mod2top,1),1));
lowMOD2   = mean(mod2bot(2:size(mod2bot,1),1));
end

% maximim value of Vxzxz (p. 381, 383)
maxVxzxz = (( (highIDV - lowIDV) / 2 ) .^2 )  *  (( (highMOD1 - lowMOD1) / 2 ) .^2 )  *  (( (highMOD2 - lowMOD2) / 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,7).^2 * Vxzxz1) ) );

% PRE for an optimal design p 384
PRE2 = 1 / ( 1+ ( mse / (b(1,7).^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(['MOD1 Mean:                                      ' num2str(mean(mod1)) ]);
disp(['MOD1 Standard Deviation:                        ' num2str(std(mod1)) ]);
disp(['MOD1 Lowest Score:                              ' num2str(MOD1min) ]);
disp(['MOD1 Highest Score:                             ' num2str(MOD1max) ]);disp([' ']);
disp(['MOD2 Mean:                                      ' num2str(mean(mod2)) ]);
disp(['MOD2 Standard Deviation:                        ' num2str(std(mod2)) ]);
disp(['MOD2 Lowest Score:                              ' num2str(MOD2min) ]);
disp(['MOD2 Highest Score:                             ' num2str(MOD2max) ]);disp([' ']);
disp(['Specified # of chunks for IDV:                  ' num2str(chunkIDV) ]);
disp(['Specified # of chunks for MOD1:                 ' num2str(chunkMOD1) ]);
disp(['Specified # of chunks for MOD2:                 ' num2str(chunkMOD2) ]);disp([' ']);
disp(['IDV Low  Value for an Optimal Design:           ' num2str(lowIDV) ]);
disp(['IDV High Value for an Optimal Design:           ' num2str(highIDV) ]);disp([' ']);
disp(['MOD1 Low  Value for an Optimal Design:          ' num2str(lowMOD1) ]);
disp(['MOD1 High Value for an Optimal Design:          ' num2str(highMOD1) ]);disp([' ']);
disp(['MOD2 Low  Value for an Optimal Design:          ' num2str(lowMOD2) ]);
disp(['MOD2 High Value for an Optimal Design:          ' num2str(highMOD2) ]);disp([' ']);disp([' ']);
disp(['Vxzxz1 (residual variance of the product):      ' num2str(Vxzxz1) ]);disp([' ']);
disp(['Vxzxz2 (residual variance of the product):      ' 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 IDV:                      ' num2str(b(1,1)) ]);disp([' ']);
disp(['Slope coefficient for MOD1:                     ' num2str(b(1,2)) ]);disp([' ']);
disp(['Slope coefficient for MOD2:                     ' num2str(b(1,3)) ]);disp([' ']);
disp(['Slope coefficient for IDV * MOD1:               ' num2str(b(1,4)) ]);disp([' ']);
disp(['Slope coefficient for IDV * MOD2:               ' num2str(b(1,5)) ]);disp([' ']);
disp(['Slope coefficient for MOD1 * MOD2:              ' num2str(b(1,6)) ]);disp([' ']);
disp(['Slope coefficient for IDV * MOD1 * MOD2:        ' num2str(b(1,7)) ]);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 mod1 mod2 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;
k4 = fix( (bigN - luper + 1) * rand(1) + 1 )  + luper - 1;
d1 = permdata(luper,1);
d2 = permdata(luper,2);
d3 = permdata(luper,3);
d4 = permdata(luper,4);
permdata(luper,1) = permdata(k1,1);
permdata(luper,2) = permdata(k2,2);
permdata(luper,3) = permdata(k3,3);
permdata(luper,4) = permdata(k4,4);
permdata(k1,1) = d1;
permdata(k2,2) = d2;
permdata(k3,3) = d3;
permdata(k3,4) = d4;
end

idv12 = permdata(:,1);
mod12 = permdata(:,2);
mod22 = permdata(:,3);
dv12  = permdata(:,4);

% Moderated Regression
x12 = idv12  .* mod12;
x22 = idv12  .* mod22;
x32 = mod12  .* mod22;
x42 = idv12  .* mod12 .* mod22;
cr2 = corrcoef([ idv12 mod12 mod22 x12 x22 x32 x42 dv12 ]);
beta2 = inv(cr2(1:7,1:7)) * cr2(1:7,8);
r2all2  = beta2' * cr2(1:7,8);
r2main2 =  ( inv(cr2(1:6,1:6))*cr2(1:6,8))'  *  cr2(1:6,8);
F2 = (r2all2-r2main2) / ((1-r2all2)/(bigN-7-1));
if F2 >= F;  counter = counter + 1;end

% Moderated Regression including quadratic terms
idvq2 = idv12 .* idv12;
mod1q2 = mod12 .* mod12;
mod2q2 = mod22 .* mod22;
crq2 = corrcoef([ idv12 mod12 mod22 idvq2 mod1q2 mod2q2 x12 x22 x32 x42 dv12 ]);
betaq2 = inv(crq2(1:10,1:10)) * crq2(1:10,11);
r2allq2  = betaq2' * crq2(1:10,11);
r2mainq2 =  ( inv(crq2(1:9,1:9))*crq2(1:9,11))'  *  crq2(1:9,11);
Fq2 = (r2allq2-r2mainq2) / ((1-r2allq2)/(bigN-10-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 IDV:                      ' num2str(bq(1,1)) ]);disp([' ']);
disp(['Slope coefficient for MOD1:                     ' num2str(bq(1,2)) ]);disp([' ']);
disp(['Slope coefficient for MOD2:                     ' num2str(bq(1,3)) ]);disp([' ']);
disp(['Slope coefficient for the IDV quadratic term:   ' num2str(bq(1,4)) ]);disp([' ']);
disp(['Slope coefficient for the MOD1 quadratic term:  ' num2str(bq(1,5)) ]);disp([' ']);
disp(['Slope coefficient for the MOD2 quadratic term:  ' num2str(bq(1,6)) ]);disp([' ']);
disp(['Slope coefficient for IDV * MOD1:               ' num2str(bq(1,7)) ]);disp([' ']);
disp(['Slope coefficient for IDV * MOD2:               ' num2str(bq(1,8)) ]);disp([' ']);
disp(['Slope coefficient for MOD1 * MOD2:              ' num2str(bq(1,9)) ]);disp([' ']);
disp(['Slope coefficient for IDV * MOD1 * MOD2:        ' num2str(bq(1,10)) ]);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([' '])
