## Copyright (C) 2004  Rafael Laboissiere
##
## This program is free software. You can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 2, or (at your option)
## any later version.
##
## This program is distributed in the hope that it will be useful, but
## WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
## General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this program; see the file COPYING.  If not, write to the
## Free Software Foundation, 59 Temple Place - Suite 330, Boston, MA
## 02111-1307, USA.

## -*- texinfo -*-
## @deftypefn {Function File} {[@var{threshold}, @var{slope}, @var{chisquare_fitting}, @var{dof_fitting}, @var{pvalue_fitting}, @var{threshold_indep}, @var{slope_indep}, @var{chisquare_parallel}, @var{dof_parallel}, @var{pvalue_parallel}, @var{t_threshold}, @var{pvalue_threshold}] =} probit_analysis (@var{control_var}, @var{batch_size}, @var{success_count})
## @end deftypefn
##
## @paragraphindent 0
## Do a probit analysis of quantal response data as presented by
## Finney[1].  The data is grouped according to different conditions and
## cumulative normal curves are fitted for each conditions.  The slopes
## of the curves at P = 0.5 are forced to be the same accross all the
## conditions.  The following statistical tests are executed:
##
## @enumerate
## @item
## Goodness of fit (chi square test)
## @item
## Parallelism of the fitted curves (chi square test)
## @item
## Separation of parallel curves at P = 0.5 (Student t test)
## @end enumerate
##
## The data is provided in the three input arguments @var{control_var},
## @var{batch_size}, and @var{success_count}.  These arguments must be
## cell variables with the same length, with the elements
## @code{@var{control_var}@{i@}}, @code{@var{batch_size}@{i@}}, and
## @code{@var{success_count}@{i@}} representing the data for condition
## @code{i}.  These elements for condition @code{i} must be vectors of
## the same length whose meanings are:
##
## @table @code
## @item @var{control_var}@{i@}
## The controlled variable of the experiment.  This is usually shown as
## the horizontal axis in plots.  Notice that in some cases it may be
## better to work with the logarithm of the controlled variables, in
## particular for case where the underlying sensitivity function is
## asymetric (see [1]).  The needed transformations on @var{control_var}
## are let to the user.
## @item @var{batch_size}@{i@}
## The size of the batch in number of individuals or trials.
## @item @var{success_count}@{i@}
## The number of individuals or trials that responded positively in the
## experiment under condition @code{i}.
## @end table
##
## The meaning of the output arguments are as follows:
##
## @table @code
## @item @var{threshold}
## Vector with the same length as @var{control_var} containing the fitted
## thresholds for each condition, assuming the parallelism hypothesis
## (null hypothesis).  A threshold is the value of the controlled variable
## for which the fitted curve gives P = 0.5.
## @item @var{slope}
## Scalar with the fitted common slope for all conditions.
## @item @var{chisquare_fitting}
## Chi square value for use in the goodness of fit test.  Large values of
## @var{chisquare_fitting} suggest rejection of the probit regression fit.
## @item @var{dof_fitting}
## Number of degrees of freedom for the goodness of fit test.
## @item @var{pvalue_fitting}
## p-value for the goodness of fit test.
## @item @var{threshold_indep}
## Same as @var{threshold}, but for each condition being fitted separately.
## @item @var{slope_indep}
## Vector containing the slopes for the curves fitted individually for each
## condition.
## @item @var{chisquare_parallel}
## Chi square value for the parallel component of the model.  It is
## computed as the difference between the chi square value for the parallel
## fitting minus the total chi square for the individual fittings.  This is
## a valid procedure because the two models considered are nested. A large
## value of @var{chisquare_parallel} indicates that the null hypothesis
## (i,e. the parallelism of the curves) can be rejected.
## @item @var{dof_parallel}
## Number of degrees of freedom for the parallelism test.
## @item @var{pvalue_parallel}
## p-value of the parallelism component of the model.
## @item @var{t_threshold}
## A matrix containing the Student-t values for the fiducial intervals.
## These values are computed for each two-by-two combinations of the
## experimental conditions.  @code{@var{t_threshold}(i,j)} is the
## t value for the fiducial interval of the difference in the thresholds
## in conditions @code{i} and @code{j} being equal to zero.
## @item @var{pvalue_threshold}
## p-value for the Student-t values in @var{t_threshold}.
## @end table
##
## CAVEAT: When heterogeneity is present in the data, which arises from small
## batch sizes or from correlations among individual responses (see [1],
## section 4.6, page 71), the chi square statistics are inflated and must be
## corrected.  This is not done in the current implementation of this
## function.  Also, when either the goodness of fit or the parallelism
## null hypothesis can be rejected, the t test of difference in thresholds
## is meaningless.
##
## [1] Finney, D.J. (1971) "Probit Analysis", Cambridge University Press,
## 3rd ed.
##
## @seealso{probit_likelihood_function, probit_finney_examples.m}

## Author: Rafael Laboissiere <rafael@laboissiere.net>
## Description: Finney's probit analysis
## $Id$

function [threshold, slope, chisquare_fitting, dof_fitting, \
          pvalue_fitting, threshold_indep, slope_indep, \
          chisquare_parallel, dof_parallel, pvalue_parallel, \
          t_threshold, pvalue_threshold] \
      = probit_analysis (control_var, batch_size, success_count)

  global probit_analysis_data

  ## Input argument sanity checks

  if nargin != 3
    usage ("probit_analysis (control_var, batch_size, success_count)");
  endif

  if ! iscell (control_var) || ! iscell (batch_size) \
        || ! iscell (success_count)
    if ! is_vector (control_var) || ! is_vector (batch_size) \
          || ! is_vector (success_count)

      error ("input arguments must be all cells or all vectors");

    else
      control_var = {control_var};
      batch_size = {batch_size};
      success_count = {success_count};
    endif

  endif

  if all (size (control_var) != 1) || all (size (batch_size) != 1) \
        || all (size (success_count) != 1)
        error ("input arguments must be vector cells");
  endif

  nc = length (control_var);
  nb = length (batch_size);
  ns = length (success_count);

  if nc != nb || nc != ns || nb != ns
    error ("input arguments must be cells of same length");
  endif

  ## Fit individual conditions

  chisquare_total = 0;
  dof_total = 0;
  slope_indep = zeros (1, nc);
  threshold_indep = zeros (1, nc);
  mcv = zeros (1, nc);

  for i = 1 : nc

    ci = control_var {i};
    bi = batch_size {i};
    si = success_count {i};

    if ! is_vector (ci) || ! is_vector (bi) || ! is_vector (si)
      error (sprintf ("vectors not found for condition %d", i));
    endif

    mc = length (ci);
    mb = length (bi);
    ms = length (si);

    if mc != mb || mc != ms || mb != ms
      error (sprintf ("in condition %d, vectors have different lengths", i));
    endif

    dof_total = dof_total + mc - 2;

    mcv (i) = mean (control_var {i});
    x = fmins ("probit_likelihood_function", [(mcv (i)), 0],
               [], [], {control_var {i}}, {batch_size {i}}, \
               {success_count {i}});
    [L, P, Q, y, t, s, chisq] \
        = probit_likelihood_function (x, {control_var {i}}, \
               {batch_size {i}}, {success_count {i}});
    threshold_indep (i) = t;
    slope_indep (i) = s;
    chisquare_total = chisquare_total + chisq;

  endfor

  ## Global fit with parallelism assumption

  x = fmins ("probit_likelihood_function", \
             [threshold_indep, (mean (slope_indep))], [], [], \
             control_var, batch_size, success_count);
  [L, P, Q, y, threshold, slope, chisquare_fitting] \
      = probit_likelihood_function (x, control_var, batch_size, \
                                    success_count);

  dof_fitting = dof_total + nc - 1;
  pvalue_fitting = 1 - chisquare_cdf (chisquare_fitting, dof_fitting);

  dof_parallel = dof_fitting - dof_total;
  chisquare_parallel = chisquare_fitting - chisquare_total;
  pvalue_parallel = 1 - chisquare_cdf (chisquare_parallel, dof_parallel);

  ## Fiducial limits

  t_threshold = zeros (nc, nc);

  for i = 1 : nc

    z = control_var {i};
    n = batch_size {i};

    w = (1 / (2 * pi)) * exp (- y {i} .^ 2) ./ (P {i} .* Q {i});
    Sw = sum (w);

    m1 = threshold (i);
    zm1 = sum (w .* z) / Sw;

    Vb1 = sum (n .* w .* (z - zm1) .^ 2);
    Vy1 = 1 / sum (n .* w);

    for j = i + 1 : nc

      z = control_var {j};
      n = batch_size {j};

      w = (1 / (2 * pi)) * exp (- y {j} .^ 2) ./ (P {j} .* Q {j});
      Sw = sum (w);

      m2 = threshold (j);
      zm2 = sum (w .* z) / Sw;

      Vb2 = sum (n .* w .* (z - zm2) .^ 2);
      Vy2 = 1 / sum (n .* w);

      Vb = 1 / (Vb1 + Vb2);
      Vy = Vy1 + Vy2;

      if m1 < m2
        [m1, m2] = swap (m1, m2);
        [zm1, zm2] = swap (zm1, zm2);
      endif

      probit_analysis_data.m = m1 - m2;
      probit_analysis_data.zm = zm1 - zm2;
      probit_analysis_data.b = slope;
      probit_analysis_data.Vb = Vb;
      probit_analysis_data.Vy = Vy;

      t_threshold (i, j) = fsolve ("probit_fiducial_interval", 0);
      t_threshold (j, i) = t_threshold (i, j);

    endfor
  endfor

  ## The t variable under the null hypothesis is a normal deviate with
  ## one degree of freedom.

  pvalue_threshold = 1 - erf (t_threshold / sqrt(2));

endfunction

%!demo
%! type probit_finney_example
%! probit_finney_example
%! fprintf ("Results:\n");
%! [threshold, slope, chisquare_fitting, dof_fitting, \
%!          pvalue_fitting, threshold_indep, slope_indep, \
%!          chisquare_parallel, dof_parallel, pvalue_parallel, \
%!          t_threshold, pvalue_threshold] \
%!      = probit_analysis (control_var, batch_size, success_count)
%! n = length (control_var);
%! minx = Inf;
%! maxx = -Inf;
%! for i = 1 : n
%!   minx = min ([minx, (control_var {i})]);
%!   maxx = max ([maxx, (control_var {i})]);
%! endfor
%! dx = (maxx - minx) / 10;
%! x = linspace (minx - dx, maxx + dx, 100);
%! y = linspace (0, 1, 50);
%! figure;
%! gset key right bottom
%! eval (sprintf ("gset xrange [%f:%f]", minx - dx, maxx + dx))
%! color = {"r", "b", "g", "k"};
%! subplot (2, 1, 1);
%! hold off
%! for i = 1 : n
%!   plot (control_var {i}, success_count {i} ./ batch_size {i}, \
%!         sprintf ("%s*;;", color {i}));
%!   hold on
%!   [L, P] = probit_likelihood_function ([(threshold_indep (i)), \
%!                                         (slope_indep (i))], {x});
%!   plot (x, P {1}, sprintf ("%s;;", color {i}));
%!   plot (threshold_indep (i) * ones (size (y)), y, \
%!         sprintf ("%s.;;", color {i}));
%! endfor
%! title ("Individual fit");
%! xlabel ("log dose");
%! ylabel ("Prob");
%! X = cell (1, n);
%! for i = 1 : n
%!   X {i} = x;
%! endfor
%! [L, P] = probit_likelihood_function ([threshold , slope], X);
%! subplot (2, 1, 2);
%! hold off
%! for i = 1 : n
%!   plot (control_var {i}, success_count {i} ./ batch_size {i}, \
%!         sprintf ("%s*;%s;", color {i}, condition_label {i}));
%!   hold on
%!   plot (x, P {i}, sprintf ("%s;;", color {i}));
%!   plot (threshold (i) * ones (size (y)), y, sprintf ("%s.;;", color {i}));
%! endfor
%! title ("Parallel fit");
%! xlabel ("log dose");
%! ylabel ("Prob");
