function [adpm] = interact3(d)
% gm = grandmean
gm = matmean(d);

% compute main effects
e = effects(d);

% compute first order interactions
[ad, ap, am, dp, dm, pm] = interact1(d);

% compute second-order interactions
[adp, adm, apm, dpm] = interact2(d);

adpm=zeros(3,3,3,3);

% compute third-order interaction

for i =1:3
    for j = 1:3
        for k = 1:3
            for l = 1:3
             m=d(i,j,k,l);
             main = e(1,i)+e(2,j)+e(3,k)+e(4,l);
             first = ad(i,j)+ap(i,k)+am(i,l)+dp(j,k)+dm(j,l)+pm(k,l);
             second = adp(i,j,k)+adm(i,j,l)+apm(i,k,l)+dpm(j,k,l);
             adpm(i,j,k,l)=m-gm-main-first-second;
            end
        end
    end
end