% This file refers to section 5 of:
% Technical Report RP-SOR-0601, Robust Experimental Design for Multivariate Generalized Linear Models, Hovav A. Dror and David M. Steinberg, January 2006.
% Available at http://www.math.tau.ac.il/~dms/GLM_Design

% It compares cluster design to
% Dr. Dave Woods' Example from section 6 (Example 2) of his paper: Woods, D.C.,
% Lewis, S.M., Eccleston, J.A. and Russell, K.G. (2005) Designs for generalized linear models with several variables and model uncertainty, Technometrics, to appear.

% WFC is Woods' Compromise Design
% WF1 through WF4 are the 4 local designs (by woods)
% WF5 is all the 4 local designs together

% Model 1	probit	3.0 + 1.6x1 + 4.1x2
% Model 2	probit	1.2 + 1.7x1 + 5.4x2 -1.7x1x2
% Model 3	comp	3.0 + 1.6x1 + 4.1x2
% Model 4	comp 	1.2 + 1.7x1 + 5.4x2 -1.7x1x2


family='binomial';
beta1=[3 1.6 4.1]; beta3=beta1; 
beta2=[1.2 1.7 5.4 -1.7]; beta4=beta2;
model1='linear'; model2='i'; model3='linear'; model4='i';
link1='probit'; link2='probit'; link3='cll'; link4='cll';
EFF=[];

% Producing Local Designs:
n=12; % n is the number of points in each local D-optimal design
            % Although we need a compromise design of 6 points
            % it is possible that we better cluster bigger Local designs
            % the number is chosen such that each local design reaches
            % its maximum efficiency
WF1=DGLM(beta1,2,n,model1,family,link1,1e-4);
WF2=DGLM(beta2,2,n,model2,family,link2,1e-4);
WF3=DGLM(beta3,2,n,model3,family,link3,1e-4);
WF4=DGLM(beta4,2,n,model4,family,link4,1e-4);

WFC=[-1	-0.102544 % This is Woods' Compromise Design
0.517194	-0.292961
-1	-0.593289
1	-0.999
-1	0.188718
1	-1
];

% Show that there was a reason for re-calculating  the local Designs
% Woods Local Designs:
W1=[-1	-0.619047 ;-1	-0.619048 ;1	-1 ;-1	-0.063877 ;1	-1 ;-1	-0.0638778 ];
W2=[0.671292	-0.297445 ; 1	-1 ; -0.0946724	-0.352366 ; 1	-1 ; -1	-0.0526539 ; -1	0.218024 ];
W3=[-1	-0.66774 ; -1	-0.667742 ; -1	-0.102529 ; 1	-1 ; -1	-0.102528 ; 1	-1 ];
W4=[-0.542084	-0.235569 ; 1	-0.547392 ; -1	0.194956 ; -1	0.19496 ; 0.352757	-0.202024 ; 1	-1 ];

% This compares Woods' 6-point Local D-optimal designs, to the DGLM created
%               n-point local designs
EFF(1)=( det(InfoMtrxGLM(x2fx(W1,model1),beta1,family,link1)) / det(InfoMtrxGLM(x2fx(WF1,model1),beta1,family,link1)) )^(1/3);
EFF(2)=( det(InfoMtrxGLM(x2fx(W2,model2),beta2,family,link2)) / det(InfoMtrxGLM(x2fx(WF2,model2),beta2,family,link2)) )^(1/4);
EFF(3)=( det(InfoMtrxGLM(x2fx(W3,model3),beta3,family,link3)) / det(InfoMtrxGLM(x2fx(WF3,model3),beta3,family,link3)) )^(1/3);
EFF(4)=( det(InfoMtrxGLM(x2fx(W4,model4),beta4,family,link4)) / det(InfoMtrxGLM(x2fx(WF4,model4),beta4,family,link4)) )^(1/4);
EFF=EFF*n/6;
fprintf('Efficiency of Woods Local D-optimal Designs vs. DGLM Local D-optimal designs: \n');
disp(EFF);



% Calculate Local Efficiencies for Compromise design:
    %EFF(1)=( det(InfoMtrxGLM(x2fx(WFC,model1),beta1,family,link1)) / det(InfoMtrxGLM(x2fx(WF1,model1),beta1,family,link1)) )^(1/3);
    %EFF(2)=( det(InfoMtrxGLM(x2fx(WFC,model2),beta2,family,link2)) / det(InfoMtrxGLM(x2fx(WF2,model2),beta2,family,link2)) )^(1/4);
    %EFF(3)=( det(InfoMtrxGLM(x2fx(WFC,model3),beta3,family,link3)) / det(InfoMtrxGLM(x2fx(WF3,model3),beta3,family,link3)) )^(1/3);
    %EFF(4)=( det(InfoMtrxGLM(x2fx(WFC,model4),beta4,family,link4)) / det(InfoMtrxGLM(x2fx(WF4,model4),beta4,family,link4)) )^(1/4);
% efficiency calculated vs. woods local designs
EFF(1)=( det(InfoMtrxGLM(x2fx(WFC,model1),beta1,family,link1)) / det(InfoMtrxGLM(x2fx(W1,model1),beta1,family,link1)) )^(1/3);
EFF(2)=( det(InfoMtrxGLM(x2fx(WFC,model2),beta2,family,link2)) / det(InfoMtrxGLM(x2fx(W2,model2),beta2,family,link2)) )^(1/4);
EFF(3)=( det(InfoMtrxGLM(x2fx(WFC,model3),beta3,family,link3)) / det(InfoMtrxGLM(x2fx(W3,model3),beta3,family,link3)) )^(1/3);
EFF(4)=( det(InfoMtrxGLM(x2fx(WFC,model4),beta4,family,link4)) / det(InfoMtrxGLM(x2fx(W4,model4),beta4,family,link4)) )^(1/4);


WF5= [WF1 ; WF2 ; WF3 ; WF4 ];

WF5=WF5-sign(WF5).*rand(size(WF5)).*0.0001; % Jittering

bestD=-inf;
N=100; % Number of cluster repetitions
for i=1:N
    [IDX,KF,s] =kmeans(WF5,6,'replicates',1,'maxiter',500,'emptyaction','singleton','distance','cityblock','display','off');
        
    thisR=0;
    
    F=x2fx(KF,model1);
    thisR=thisR+log(det(InfoMtrxGLM(F,beta1,family,link1)));
    F=x2fx(KF,model2);
    thisR=thisR+log(det(InfoMtrxGLM(F,beta2,family,link2)));
    F=x2fx(KF,model3);
    thisR=thisR+log(det(InfoMtrxGLM(F,beta3,family,link3)));
    F=x2fx(KF,model4);
    thisR=thisR+log(det(InfoMtrxGLM(F,beta4,family,link4)));
    
    if thisR>bestD % if this cluster improves
        bestD=thisR;
        KFB=KF; % This is the best cluster design
    end;
    
end; % for i

% Now calculate the efficiency of the cluster design   
       % EFF(2,1)=( det(InfoMtrxGLM(x2fx(KFB,model1),beta1,family,link1)) / det(InfoMtrxGLM(x2fx(WF1,model1),beta1,family,link1)) )^(1/3);
        %EFF(2,2)=( det(InfoMtrxGLM(x2fx(KFB,model2),beta2,family,link2)) / det(InfoMtrxGLM(x2fx(WF2,model2),beta2,family,link2)) )^(1/4);
        %EFF(2,3)=( det(InfoMtrxGLM(x2fx(KFB,model3),beta3,family,link3)) / det(InfoMtrxGLM(x2fx(WF3,model3),beta3,family,link3)) )^(1/3);
        %EFF(2,4)=( det(InfoMtrxGLM(x2fx(KFB,model4),beta4,family,link4)) / det(InfoMtrxGLM(x2fx(WF4,model4),beta4,family,link4)) )^(1/4);
%  efficiency vs. Woods' local designs (or - best 6 point local design)
EFF(2,1)=( det(InfoMtrxGLM(x2fx(KFB,model1),beta1,family,link1)) / det(InfoMtrxGLM(x2fx(W1,model1),beta1,family,link1)) )^(1/3);
EFF(2,2)=( det(InfoMtrxGLM(x2fx(KFB,model2),beta2,family,link2)) / det(InfoMtrxGLM(x2fx(W2,model2),beta2,family,link2)) )^(1/4);
EFF(2,3)=( det(InfoMtrxGLM(x2fx(KFB,model3),beta3,family,link3)) / det(InfoMtrxGLM(x2fx(W3,model3),beta3,family,link3)) )^(1/3);
EFF(2,4)=( det(InfoMtrxGLM(x2fx(KFB,model4),beta4,family,link4)) / det(InfoMtrxGLM(x2fx(W4,model4),beta4,family,link4)) )^(1/4);

fprintf('\n\nTop row: Compromise Efficiency vs. the 4 Local D-optimal designs \n');
fprintf('Second row: Cluster Design Efficiency \n');
fprintf('                    Efficiencies are vs. DGLM computed local designs \n');
disp(EFF)