function [ad, ap, am, dp, dm, pm] = interact1(d)
% gm = grandmean
gm = matmean(d);

% compute main effects
e = effects(d);

% compute first-order interactions

% interactions for AD
for i =1:3
    for j = 1:3
    m=matmean(d(i,j,:,:));
    ad(i,j) = m - gm - e(1,i) - e(2,j);
    end
end
% interactions for AP
for i =1:3
    for k = 1:3
    m=matmean(d(i,:,k,:));
    ap(i,k) = m - gm - e(1,i) - e(3,k);
    end
end
% interactions for AM
for i =1:3
    for l = 1:3
    m=matmean(d(i,:,:,l));
    am(i,l) = m - gm - e(1,i) - e(4,l);
    end
end
% interactions for DP
for j =1:3
    for k = 1:3
    m=matmean(d(:,j,k,:));
    dp(j,k) = m - gm - e(2,j) - e(3,k);
    end
end
% interactions for DM
for j =1:3
    for l = 1:3
    m=matmean(d(:,j,:,l));
    dm(j,l) = m - gm - e(2,j) - e(4,l);
    end
end
% interactions for PM
for k =1:3
    for l = 1:3
    m=matmean(d(:,:,k,l));
    pm(k,l) = m - gm - e(3,k) - e(4,l);
    end
end