% Laboratory in Oceanography - Data and Methods
% SMAST, UMass Dartmouth
%
% Examples from statistics toolbox:
%   ANOVA analysis of light vs. DO 
%
% NOTE: light and DO data not included online due to proprietary data.
%
% Original code written by Miles A. Sundermeyer, 4/10/08
% Last modifid for Lab in Oceanography 3/9/09

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for n=1:8
  figure(n)
  clf;
end

combineflag = 1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[NUMERIC,TXT,RAW] = xlsread('Light_DO_data_for_Brian.xls','raw data');

% NUMERIC has columns:
% ------------------
% 1 = nan
%
% 2 = Jul 06 light
% 3 = Jul 06 % high DO
%
% 4 = nan
%
% 5 = Aug 06 light
% 6 = Aug 06 % high DO
%
% 7 = nan
%
% 8 = Both 06 light
% 9 = Both 06 % high DO
%
% 10 = nan
%
% 11 = Jul 06 light
% 12 = Jul 06 % low DO
%
% 13 = nan
%
% 14 = Aug 06 light
% 15 = Aug 06 % low DO
%
% 16 = nan
%
% 17 = Both 06 light
% 18 = Both 06 % low DO

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% clear some variables

light = nan*ones(size(NUMERIC,1),6);
pctHigh = nan*ones(size(NUMERIC,1),6);
pctLow = nan*ones(size(NUMERIC,1),6);

lightind = [2 5 8];
highind = [3 6 9];
lowind = [12 15 18];

for n=[1:3]
  light(:,n) = NUMERIC(:,lightind(n));
  pcthigh(:,n) = NUMERIC(:,highind(n));
  pctlow(:,n) = NUMERIC(:,lowind(n));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1)
clf

plot(light(:,1),pcthigh(:,1),'ro')
hold on
plot(light(:,2),pcthigh(:,2),'rs')

xlabel('Daily Avg. Light (day * uM/m^2/s)')
ylabel('% Avail Locations w/ events')
title('High DO events (Jul=o, Aug=[])')

set(gca,'xlim',[50 750]);
set(gca,'ylim',[0 100]);

grid on

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(2)
clf

plot(light(:,1),pctlow(:,1),'bo')
hold on
plot(light(:,2),pctlow(:,2),'bs')

xlabel('Daily Avg. Light (day * uM/m^2/s)')
ylabel('% Avail Locations w/ events')
title('Low DO events (Jul=o, Aug=[])')

set(gca,'xlim',[50 750]);
set(gca,'ylim',[0 100]);

grid on

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Run a one-way ANOVA - see: Matlab web help
% First put the data in a right format, with columns representing bins of data on the x-axis

nbins = 6;
binsize = 100;
thisdataJH = nan*ones(20,nbins);
thisdataAH = nan*ones(20,nbins);
thisdataJL = nan*ones(20,nbins);
thisdataAL = nan*ones(20,nbins);
thisdataBH = nan*ones(30,nbins);
thisdataBL = nan*ones(30,nbins);

for n=1:nbins		% fit bins from 100 - 700 light level
  lightbins(n) = n*binsize;

  % Jul
  thisind = find((light(:,1) >= n*binsize) & (light(:,1) < (n+1)*binsize));
  if ~isempty(thisind)
    thisdataJH(1:length(thisind),n) = pcthigh(thisind,1);
    thisdataJL(1:length(thisind),n) = pctlow(thisind,1);
  end

  % Aug
  thisind = find((light(:,2) >= n*binsize) & (light(:,2) < (n+1)*binsize));
  if ~isempty(thisind)
    thisdataAH(1:length(thisind),n) = pcthigh(thisind,2);
    thisdataAL(1:length(thisind),n) = pctlow(thisind,2);
  end

  % both
  thisind = find((light(:,3) >= n*binsize) & (light(:,3) < (n+1)*binsize));
  if ~isempty(thisind)
    thisdataBH(1:length(thisind),n) = pcthigh(thisind,3);
    thisdataBL(1:length(thisind),n) = pctlow(thisind,3);
  end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
meanJH = nanmean(thisdataJH);
meanJL = nanmean(thisdataJL);
meanAH = nanmean(thisdataAH);
meanAL = nanmean(thisdataAL);
meanBH = nanmean(thisdataBH);
meanBL = nanmean(thisdataBL);

stddevJH = nanstd(thisdataJH);
stddevJL = nanstd(thisdataJL);
stddevAH = nanstd(thisdataAH);
stddevAL = nanstd(thisdataAL);
stddevBH = nanstd(thisdataBH);
stddevBL = nanstd(thisdataBL);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
hold on
if combineflag
  figure(1)
  h = fill([lightbins fliplr(lightbins)],[meanBH+stddevBH fliplr(meanBH-stddevBH)],[1 0.75 0.75]);
  set(h,'edgecolor','none');
  plot(lightbins,meanBH,'r.-')

  figure(2)
  h = fill([lightbins fliplr(lightbins)],[meanBL+stddevBL fliplr(meanBL-stddevBL)],[0.75 0.75 1]);
  set(h,'edgecolor','none');
  plot(lightbins,meanBL,'b.-')
else
  figure(1)
  h = fill([lightbins fliplr(lightbins)],[meanJH+stddevJH fliplr(meanJH-stddevJH)],[1 0.75 0.75]);
  set(h,'edgecolor','none');
  h = fill([lightbins fliplr(lightbins)],[meanAH+stddevAH fliplr(meanAH-stddevAH)],[1 0.75 0.75]);
  set(h,'edgecolor','none');
  plot(lightbins,meanJH,'r.-')
  plot(lightbins,meanAH,'b.-')

  figure(2)
  h = fill([lightbins fliplr(lightbins)],[meanJL+stddevJL fliplr(meanJL-stddevJL)],[0.75 0.75 1]);
  set(h,'edgecolor','none');
  h = fill([lightbins fliplr(lightbins)],[meanAL+stddevAL fliplr(meanAL-stddevAL)],[0.75 0.75 1]);
  set(h,'edgecolor','none');
  plot(lightbins,meanJL,'r.-')
  plot(lightbins,meanAL,'b.-')
end

% replot the original data
figure(1)
plot(light(:,1),pcthigh(:,1),'ro')
plot(light(:,2),pcthigh(:,2),'rs')

figure(2)
plot(light(:,1),pctlow(:,1),'ro')
plot(light(:,2),pctlow(:,2),'rs')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% now run Matlab's ANOVA on the full data set
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% first the highs
[p,tlb,stats] = anova1(thisdataBH);

% and do a multiple comparison test
figure(3)
[c,m] = multcompare(stats)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% do the same for the lows
[p,tlb,stats] = anova1(thisdataBL);

% and do a multiple comparison test
figure(4)
[c,m] = multcompare(stats)