function [m] = model(d)

% gm = grandmean
gm = matmean(d);

% main effects
e = effects(d);

% first-order interactions

[ad, ap, am, dp, dm, pm] = interact1(d);

% second-order interactions

[adp, adm, apm, dpm] = interact2(d);

% third-order interactions
[adpm] = interact3(d);

for i = 1:3
    for j = 1:3
        for k = 1:3
            for l = 1:3
             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);
             third = adpm(i,j,k,l);
       % full model
       %      m(i,j,k,l) = gm + main + first + second + third;
       
       % simplified model
              m(i,j,k,l) = gm + main + dm(j,l);
            end
        end
    end
end