function [pattern,model,graph] = AIClim(block, time, treat, response, factors);
% AIClim - fit a set of candidate models to a bioassay experiment with 
% 2 limiting nutrients added as replicated factorial design and sampled 
% repeatedly over time, and select the model with the lowest AIC
%
% Inputs -----------------------------------------------------------
% block    : Column vector of experiment block indicators as integers 1, 2, ...
% time     : Column vector of sampling times as successive integers 0, 1, 2, ...
% treat    : 2-column binary (0, 1) matrix indicating nutrient additions
% response : Column vector of measured response in each experimental unit
% factors  : 2-element string of liming factor names (default 'PN')
%
% Outputs ----------------------------------------------------------
% pattern  : 5-character string representing the limitation pattern with lowest AIC
%          : e.g., 'XP(2)' = Exclusive P limitation with quadratic time effect
% model    : Structure with the following elements:
%            AICmin   : AIC for the selected model
%            lim_m    : Limitation pattern effect matrix
%            time_m   : Orthogonal time effect matrix
%            reorder  : Standard reordering of elements of the response
%            design_m : General linear model design matrix for fitting
%                     : (b = design_m \ response(reorder))
% graph    : Handle to a figure showing observations and predictions
%
%   Tom Andersen 2005

if nargin < 4,
    error('Too few input arguments');
end;

if nargin < 5,
    factors = 'PN'; % Default factor names
end;

[Nsamp,Nblock] = checkinput(block,time,treat,response);

% Reorder rows of input data to standard ordering scheme
%   - treatments are sorted within sampling times
%   - sampling times sorted within replication blocks
code = treat * [1 2]'; % Collape binary matrix to treatment code vector
[dummy, stdorder] = sortrows([block,time,code]);
y = response(stdorder);

% Generate candidate contrast matrices
[Contrast,cname] = limcontrast(factors);

% Time effects as orthogonal polynomials matrix
Time = cell(Nsamp-1,1);
for i = 1:(Nsamp-1),
    Time(i) = {orthpoly(Nsamp,i+1)};
end;

Z = cell(length(Time),length(Contrast));
for i = 1:length(Time),
    for j = 1:length(Contrast),
        Ti = cell2mat(Time(i));
        Cj = cell2mat(Contrast(j));
        Zij = kron(Ti,Cj);
        X(i,j) = {repmat(Zij,Nblock,1)};
    end;
end;
[order, cnum] = ind2sub(size(X),[1:prod(size(X))]');
X = X(:);

a = zeros(length(X),1);
for j = 1:length(X),
    Xj = cell2mat(X(j));
    a(j) = AIClm(Xj, y);
end;
[amin, k] = min(a); % Select model with min. AIC

pattern = [cname(cnum(k),:) '(' num2str(order(k)) ')'];

if nargout > 1,
    model = struct('AICmin',amin,...
        'lim_m',cell2mat(Contrast(cnum(k))),...
        'time_m',cell2mat(Time(order(k))),...
        'reorder',stdorder,'design_m',cell2mat(X(k)));
end;

if nargout > 2,
    Xk = cell2mat(X(k));
    Zk = Xk(1:(4*Nsamp),:); % 1st block of X for predicting
    zk = Zk * (Xk \ y);     % Predicted y
    
    tfit = sort(unique(time));
    tobs = repmat(tfit,Nblock,1);
    yfit = exp(reshape(zk,4,Nsamp)');      % Backtransform -> linear
    yobs = exp(reshape(y,4,Nsamp*Nblock)');% Backtransform -> linear
    
    graph = plot(tobs,yobs(:,2),'og',tobs,yobs(:,3),'oy',...
        tobs,yobs(:,1),'ob',tobs,yobs(:,4),'or',...
        tfit,yfit(:,2),'-g',tfit,yfit(:,3),'-y',...
        tfit,yfit(:,1),':b',tfit,yfit(:,4),':r',...
        'markersize',5,'linewidth',4);
    set(gca,'Xlim',[-0.2 (Nsamp-0.8)],'XTick',[0:(Nsamp-1)]);
    xlabel('Time'); ylabel('Response');
    [leg, obj] = legend(['00';[factors(1) '0'];['0' factors(2)]; factors],0);
end;

%% End - AIClim %%


function [Nsamp,Nblock] = checkinput(block,time,treat,response);
% Checkinput - Check consistency of input and give error messages if necessary
%
% Nsamp  : Number of repeated samples from each unit
% Nblock : Number of replicated treatment blocks
%
%   Tom Andersen 2005

N = length(time);
n = [length(time) max(size(treat)) length(block) length(response)];
if any(n  ~= N),
    error('Inputs time, treat, block, and response must have the same number of rows');
end;

b = hist(block, unique(block));
if any(b  ~= b(1)),
    error('Unequal number of replication block samples');
end;

Nblock = length(unique(block));
if (min(block) > 1) | (mod(N,Nblock) ~= 0) | (max(block) > Nblock),
    error('Block must be successive integers from 1');
end;

b = hist(time, unique(time));
if any(b  ~= b(1)),
    error('Unequal number of time samples');
end;

Nsamp = length(unique(time));
if (min(time) > 0) |  (mod(N,Nsamp) ~= 0) | (max(time) > (Nsamp - 1)),
    error('Time must be successive integers from 0');
end;

if (min(size(treat)) ~= 2),
    error('Treatment matrix must have 2 columns');
end;

if any(treat ~= (treat > 0)),
    error('Treatments must be coded as binary');
end;

code = treat * [1 2]'; % Collape binary matrix to treatment code vector
b = hist(code, unique(code));
if any(b  ~= b(1)),
    error('Unequal number of treatments');
end;

if ~ all(isfinite(response)),
    error('Missing response valuse not allowed');
end;

%% End - checkinput %%


function [contrast,cname] = limcontrast(factors);
%
% limcontrast - Generate the 7 candidate contrast matrices 
% corresponding to interpretable limitation patterns;
%
% contrast  : 7-element cell array of matrices
% cname     : string matrix of limitation pattern codes
%
%   Tom Andersen 2005

M1 = [1 -1; 1 1];   % Elementary orthogonal contrast matrix
M2 = kron(M1,M1);   % Full factorial design contrast matrix

contrast = cell(7,1);
cname = char(double(' ') * ones(7,2));

% Pattern 1: No treatment effect (00 == 10 == 01 == 11)
contrast(1) = {[M2(:,1)]};
cname(1,:) = '00';

% Pattern 2: Exclusive factor 1 effect ((00 == 01) <> (10 == 11))
contrast(2) = {[M2(:,1) M2(:,2)]};
cname(2,:) = ['X' factors(1)];

% Pattern 3: Primary factor 1 effect ((00 == 01) <> 10 <> 11)
contrast(3) = {[M2(:,1) M2(:,2) (M2(:,3) + M2(:,4))]};
cname(3,:) = [factors(1) '1'];

% Pattern 4: Exclusive factor 2 effect ((00 == 10) <> (01 == 11))
contrast(4) = {[M2(:,1) M2(:,3)]};
cname(4,:) = ['X' factors(2)];

% Pattern 5: Primary factor 2 effect ((00 == 10) <> 01 <> 11)
contrast(5) = {[M2(:,1) M2(:,3) (M2(:,2) + M2(:,4))]};
cname(5,:) = [factors(2) '1'];

% Pattern 6: Exclusive factor 1+2 effect ((00 == 10 == 01) <> 11)
contrast(6) = {[M2(:,1)  (M2(:,2) + M2(:,3) + M2(:,4))]};
cname(6,:) = 'XC';

% Pattern 7: Combined factor 1+2 effect (00 <> 10 <> 01 <> 11)
contrast(7) = {M2};
cname(7,:) = 'C1';

%% End - limcontrast %%


function T = orthpoly(N,M);
% orthpoly - Orthogonal polynomials up to order M-1
% on a regular grid with N points (x = 0:N-1)
% (Chebyshev polynomials of a discrete variable)
%
% Computed by the recurrence relation on p. 791 in
% Abramowitz & Stegun "Handbook of Mathematical functions"
%
%   Tom Andersen 1995

T = zeros(N,M);
T(:,1) = ones(N,1);
T(:,2) = T(:,1) - 2 * [0:(N-1)]' / (N-1);

for i = 3:M,
	a1 = ((2*i - 3)  * (N - 1))  / ((i - 1) * (N - i + 1));
	a2 = ((i - 2) * (N + i - 2)) / ((i - 1) * (N - i + 1));
	T(:,i) = a1 * (T(:,2) .* T(:,i-1)) - a2 * T(:,i-2);
end;

%% End - orthpoly %%


function aic = AIClm(X,y);
% AIClm - Akaike's (1973) information criterion for a linear model y = X * b
% with correction for overfitting on small sample data sets, as suggested by
% Bedrick & Tsai (1994)
%
% Akaike, H. (1973) Information theory and an extension of the maximum
% likelihood principle. pp 267-281 in B. N. Petrov and F. Csaki (eds.)
% "2nd international symposium on Information theory", Akademia Kiado, Budapest
%
% Bedrick, E. J. & Tsai, C. L. (1994) Model selection for multivariate
% regression in small samples. Biometrics 50: 226-231
%
%   Tom Andersen 2003

n = length(y);  % Number of observations
p = size(X,2);  % Number of parameters

sigma = y' * (eye(n,n) - X * inv(X' * X) * X') * y;
aic   = log(det(sigma)) + ((n + p) / (n - p - 2));    % Small sample correction

%% End - AIClm %%