function [adp, adm, apm, dpm] = interact2(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

% interactions for ADP
for i =1:3
    for j = 1:3
         for k = 1:3
             m=matmean(d(i,j,k,:));
             adp(i,j,k)=m-gm-e(1,i)-e(2,j)-e(3,k)-ad(i,j)-ap(i,k)-dp(j,k);
         end
    end
end
% interactions for ADM
for i =1:3
    for j = 1:3
         for l = 1:3
             m=matmean(d(i,j,:,l));
             adm(i,j,l)=m-gm-e(1,i)-e(2,j)-e(4,l)-ad(i,j)-am(i,l)-dm(j,l);
         end
    end
end
% interactions for APM
for i =1:3
    for k = 1:3
         for l = 1:3
             m=matmean(d(i,:,k,l));
             apm(i,k,l)=m-gm-e(1,i)-e(3,k)-e(4,l)-ap(i,k)-am(i,l)-pm(k,l);
         end
    end
end
% interactions for DPM
for j =1:3
    for k = 1:3
         for l = 1:3
             m=matmean(d(:,j,k,l));
             dpm(j,k,l)=m-gm-e(2,j)-e(3,k)-e(4,l)-dp(j,k)-dm(j,l)-pm(k,l);
         end
    end
end
