% script m-file SungHoDataAnalysis. Analysis of SungHo's A' data. For
% description of experiment, see SungHoDataDescription.pdf on Sakai website
% in Resource folder

% SungHoData=1; % Creating variable into which to paste the data array copied
% from the Excel spreadsheet
%%
format short

DataColHeaders={['FHole2'] ['FHole4'] ['FHole6'] ['EHole2'] ['EHole4'] ['EHole6'] ...
    ['FBay2'] ['FBay4'] ['FBay6'] ['EBay2'] ['EBay4'] ['EBay6']}

% These measurements are all within subject, that is, there is an A' value
% for each subject for each condition (each column of SungHoData
%%
x=[2 4 6];

plot(x,[mean(SungHoData(:,1)) mean(SungHoData(:,2)) mean(SungHoData(:,3))],'ro-')
% Figure 1. Plot of the filled-hole data

hold on

plot(x,[mean(SungHoData(:,4)) mean(SungHoData(:,5)) mean(SungHoData(:,6))],'gs--')
% Plot of the empty-hole data. According one model, these A's 
% should be substantially lower than for the full-hole stimuli. According
% to another model, these A's come from the same distribution as the full
% ones

plot(x,[mean(SungHoData(:,7)) mean(SungHoData(:,8)) mean(SungHoData(:,9))],'bo-')
% Plot of the full-bay data

plot(x,[mean(SungHoData(:,10)) mean(SungHoData(:,11)) mean(SungHoData(:,12))],'ms--')
% Plot of the empty-bay data. Ditto above comment re models

xlim([1 7])

xlabel('Number of Stimuli in Displays')

ylabel('Mean A')

title('Sung-Ho Data')

legend('Full Hole','Empty Holes','Full Bays','Empty Bays')

% The A' decreases linearly with the number of stimuli in the displays. The
% hole stimuli are more difficult than the bay stimuli (their A's are
% lower), but it does not appear that the filled-vs empty manipulation has
% any effect
%%
% These data are within-subject. Therefore, they should be analyzed by
% looking at the differences within subjects. Let's look at the
% distributions of those differences. We want to look first at the key
% theoretical differences, which are those between stimuli that differ only
% with regard to whether the holes or bays are filled or not. There are six
% columns of these differences (FHoles2-EHoles2, FHoles4-EHoles4, etc

Diffs=[SungHoData(:,1)-SungHoData(:,4) SungHoData(:,2)-SungHoData(:,5) ...
    SungHoData(:,3)-SungHoData(:,6) SungHoData(:,7)-SungHoData(:,10) ...
    SungHoData(:,8)-SungHoData(:,11) SungHoData(:,9)-SungHoData(:,12)]

DiffColHdrs={['FH2-EH2'] ['FH4-EH4'] ['FH6-EH6'] ['FB2-EB2'] ['FB4-EB4'] ['FB6-EB6']}


%%
figure % Second figure. A plot of the histograms of the key differences
% for different numbers of stimuli
subplot(3,2,1)
hist(Diffs(:,1))
subplot(3,2,2);hist(Diffs(:,4))
subplot(3,2,3);hist(Diffs(:,2))
subplot(3,2,4);hist(Diffs(:,5))
subplot(3,2,5);hist(Diffs(:,3))
subplot(3,2,6);hist(Diffs(:,6))
subplot(3,2,1);xlim([-.6 .6])
for i=1:6
    subplot(3,2,i)
    xlim([-.6 .6])
end
% Creating common x-axis scale. Notice, again the importance of making the
% scales the same whenever one wants people (including oneself!!) to
% compare the plots in different panels on the same sheet.

subplot(3,2,5);xlabel('Filled-A - Empty-A')
subplot(3,2,6);xlabel('Filled-A - Empty-A')
subplot(3,2,1);title('Histograms of Filled A - Empty A for Holes')
subplot(3,2,2);title('Histograms of Filled A - Empty A for Bays')
subplot(3,2,1); ylabel('Instances (# stimuli=2')
subplot(3,2,3); ylabel('Instances (# stimuli=4')
subplot(3,2,5); ylabel('Instances (# stimuli=6')
% Using axis labels and titles to make it clear what each histogram shows.
% Titles above the top two panels label the columns. Ylabels at left of
% each row label the rows

% It is apparent that there is not homogeneity of variance within each cell
% of the prospective ANOVA table, because there is clearly more variability
% in the differences in the high-number-of-stimuli conditions than in the
% low-number-of-stimuli conditions. I BELIEVE (I'm not 100% sure of this)
% that this is a serious problem for the standard ANOVA because it averages
% the within cell variability across levels of the factors. That clearly is
% unwarranted in this case. Another major problem with the traditional
% approach is that it cannot give us a measure of the extent to which the
% data actually favor one of the two hypotheses we are entertaining,
% because one of them predicts a null

% We can do a preliminary check that the seemingly small differences we see
% in the graph are in fact not significant

[Hfe,Pfe]=ttest(Diffs)
% Result is:
% Hfe =
% 
%      0     0     0     0     0     0
% 
% Pfe =
% 
%     0.7510    0.7120    0.6855    0.0615    0.1205    0.1097

% So the small differences we see do not approach signficance except in one
% case (.0615). However, these numbers do not tell us that the data
% actually favor the null hypothesis (no difference). For that, we need
% Bayes Factors

%% Computing Bayes Factors
% The easy way is to use the Schwartz criterion, which avoids integrating
% over the priors. When you use this approach, you give each competing
% model its best shot, that is, you find how likely that model makes the
% data for the maximum-likelihood values of its parameters, that is, the
% values that make the data maximmally likely. Then, you penalize models
% for the number of free(!!) parameters that they employ, that is, for their
% complexity, because the more parameters in a model that are estimated
% from the data (that's what we mean by "free" parameters), the more likely
% they will make the data. The question is whether the increase in the
% likelihood of the data is commensurate with the increase in the
% complexity of the model.

% We can see from looking at the histograms of the differences that these
% differences are approximately normally distributed, and that is what we
% expect from the central limit theorem, so our model assumes normal
% distributions. The no-difference model says that the A' values for the
% filled stimuli are drawn from the same distributions as the A' values for
% the empty stimuli. Thus, the expected value of the A' difference between
% the filled and empty cases is 0. (This is a key step in the argument; be
% sure you understand it.) In this model, the only parameter to be
% estimated from the data is the variability (sigma), because the model
% says that mu_diff (the mean of the differences) is 0. The other model
% says that the mean of the difference is some non-zero value that is to be
% estimated from the data. Thus, in this model, there is a 2nd free
% parameter, namely the mean. This specification of the second model is not
% really accurate, because, in fact, the second model specifies the
% direction (sign) of the mean difference: it says that empty is harder
% than filled, so the average A' in the empty case should be lower than the
% average A' in the filled case, so if empty values are subtracted from
% corresponding(!!) filled values, the mean of the differences should be
% positive. (Corresponding means from the same subject and for displays with
% the same number of stimuli. However, to test this specification, we cannot
% use the baby Bayesian (shortcut) approach we will now describe. We have to
% integrate over the priors. That is computationally much more difficult, and,
% also, conceptually, somewhat of a can of worms, because you have to
% specify a probability distribution for possible values of the parameters
% ("specify the priors").

% We can see from the histograms in the second figure that we must
% separately analyze the cases for different numbers of stimuli in the
% displays, because the variability of the distributions of the differences
% clearly increases as the A's get lower.

for i=1:6 % stepping through the 6 columns of the Diffs array
    
    SigmaM1(i)=sqrt(sum(Diffs(:,i).^2)/length(Diffs(:,1)));
    % For a normally distributed data of specified mean, the value of sigma
    % that maximizes the likelihood of those data is the mean square
    % deviation of those data from the specified mean. In this case, the
    % specified mean is 0, so the deviations simply are the differences.
    % Thus, the estimate of sigma that maximizes the likelihood of the data
    % in a column is simply the mean of the squares of the data in that
    % column
    
    Params(i,:)=mle(Diffs(:,i)); % Maximum likelihood values for the mean
    % and sigma of each column under Model2, which assumes that the means
    % have some value that is to be estimated from the data. Notice that
    % Model 1 is nested within Model 2, because it stipulates that the mean
    % of the differences is 0, which is one among the infinity of possible
    % values for this mean in Model 2. Thus, Model 1 is a special case of
    % Model 2. Params has two columns: the first contains the maximum
    % likelihood value for the mean; the second, the maximum likelihood
    % value for sigma
    
    MLLM1(i)=sum(log(normpdf(Diffs(:,i),0,SigmaM1(i)))); % The maximum log
    % likelihood for Model 1
    
    MLLM2(i)=sum(log(normpdf(Diffs(:,i),Params(i,1),Params(i,2)))); % The maximum
    % log likelihood for Model 2
    
    BF(i)=MLLM1(i)-MLLM2(i)-.5*(1-2)*log(length(Diffs(:,i))); % The Schwarz criterion
    % says that the log of the Bayes Factor is the difference in the
    % maximum log likelihoods of the two models [MLLM1(i)-MLLM2(i)]
    % minus a penalty for complexity.
    % The penalty for complexity (extra free parameters) is:
    % -.5*(1-2)*log(length(Diffs(:,i))). Notice that it depends on both the
    % difference in the number of free parameters between the two models
    % [(1-2)] AND the log of the number of data being "explained" The
    % Bayes Factor is the factor by which the data alter the prior odds.
    % That is, the posterior odds (the odds in favor of one model or the
    % other AFTER the data have been considered) are the prior odds (the
    % odds BEFORE the data are considered) times the Bayes Factor. The
    % Bayes Factor is the antilog of the result returned by the calculation
    % specified in this command line.
    
end
    
BF % displaying the logs of the Bayes factors
% Result is:
% BF =
% 
%     1.5353    1.5164    1.5016    1.5741    1.5375    1.5465

% In every comparison, the Bayes Factor favors Model 1. The overall
% likelihood of Model 1 is obtained by summing these (natural) logs and
% taking the antilog of the sum:

format bank

exp(sum(BF))

% Result is:

% ans =
% 
%       10010.08

% The Bayes Factor favors Model 1 by 10,000:1, and this calculation
% seriously underestimates the posterior odds in favor of Model 1, because
% we specified Model 2 in such a way that negative mean differences counted
% in its favor, when, in fact, only positive mean differences would count
% in its favor. However, adding up the log likelihoods from the different
% comparisons (= multiplying the Bayes Factors) is valid only if subjects'
% performance under one condition is independent of their performance under
% the other conditions, which is unlikely to be true. We'll check on that
% later

% It is instructive to look at the simple differences in the maximum log
% likelihoods:

format short

MLLM1-MLLM2

% ans =
% 
%  -0.0537   -0.0727   -0.0874   -0.0149   -0.0515   -0.0425

% which translate into odds of:
format bank

1/exp(sum(MLLM1-MLLM2))

% ans =
% 
%           1.38    (again, assuming independence, which is doubtful)

% The differences in maximum log likelihood are all negative, meaning that,
% before correction for complexity, the odds favor the more complex model.
% They almost always do, as we learned when we did our first relative
% likelihood analysis. In this case, they only favor it very slightly,
% because the difference in complexity is modest (only a factor of 2 in the
% number of free parameters); whereas, in the previous case, the difference
% in complexity was huge

%% The difference between holes and bays
% From Figure 1, it appears that the displays with holes may be more difficult
% to discriminate (have lower A' scores) than the displays with bays. We
% apply the same methods to compute the relative likelihood of this model
% (the mean hole-bay difference is non-zero) and the model that assumes
% that there is no difference between hole stimuli and bay stimuli

HBDiffs=[SungHoData(:,1)-SungHoData(:,7) SungHoData(:,2)-SungHoData(:,8) ...
    SungHoData(:,3)-SungHoData(:,9) SungHoData(:,4)-SungHoData(:,10) ...
    SungHoData(:,5)-SungHoData(:,11) SungHoData(:,6)-SungHoData(:,12)]

% As a first check on whether we should believe the differences between
% holes and bays that we see in the graph is to use ttest to test whether
% the means of the differences differ significantly from 0
format short

[H,P]=ttest(HBDiffs) % This t-test commands tests whether the mean differs
% from 0. 

% Result is:
% H =
% 
%      0     0     0     0     0     0
% 
% P =
% 
%     0.0639    0.2475    0.3168    0.0694    0.0571    0.1075

% so none of the differences is significant


%%
figure % Third figure. A plot of the histograms of the hole-bay differences
% for different numbers of stimulus elements
subplot(3,2,1)
hist(HBDiffs(:,1)) % this time we normalize histograms
subplot(3,2,2);hist(HBDiffs(:,4))
subplot(3,2,3);hist(HBDiffs(:,2))
subplot(3,2,4);hist(HBDiffs(:,5))
subplot(3,2,5);hist(HBDiffs(:,3))
subplot(3,2,6);hist(HBDiffs(:,6))
subplot(3,2,1);xlim([-.6 .6])
for i=1:6
    subplot(3,2,i)
    xlim([-.6 .6])
end
% Creating common x-axis scale. Notice, again the importance of making the
% scales the same whenever one wants people (including oneself!!) to
% compare the plots in different panels on the same sheet.

subplot(3,2,5);xlabel('FHA'' - FBA''')
subplot(3,2,6);xlabel('EHA'' - EBA''')
subplot(3,2,1);title('A''for filled holes - A''for filled bays')
subplot(3,2,2);title('A''for empty holes - A''for empty bays')
subplot(3,2,1); ylabel('Rel Freq (# stimuli=2')
subplot(3,2,3); ylabel('Rel Freq (# stimuli=4')
subplot(3,2,5); ylabel('Rel Freq (# stimuli=6')

% As we expect from our first graph, these histograms of the within-subject
% hole-bay differences appear to have their means somewhat on the negative
% side of 0 (we subtracted the bay A's (the higher ones) from the hole A's)
%%
for i=1:6 % stepping through the 6 columns of the HBDiffs array
    
    SigHBM1(i)=sqrt(sum(HBDiffs(:,i).^2)/length(HBDiffs(:,1)));
    % For a normally distributed data of specified mean, the value of sigma
    % that maximizes the likelihood of those data is the mean square
    % deviation of those data from the specified mean. In this case, the
    % specified mean is 0, so the deviations simply are the differences.
    % Thus, the estimate of sigma that maximizes the likelihood of the data
    % in a column is simply the mean of the squares of the data in that
    % column
    
    ParamsHB(i,:)=mle(HBDiffs(:,i)); % Maximum likelihood values for the mean
    % and sigma of each column under Model2, which assumes that the means
    % have some value that is to be estimated from the data. Notice that
    % Model 1 is nested within Model 2, because it stipulates that the mean
    % of the differences is 0, which is one among the infinity of possible
    % values for this mean in Model 2. Thus, Model 1 is a special case of
    % Model 2. Params has two columns: the first contains the maximum
    % likelihood value for the mean; the second, the maximum likelihood
    % value for sigma
    
    MLLM1HB(i)=sum(log(normpdf(HBDiffs(:,1),0,SigHBM1(i)))); % The maximum log
    % likelihood for Model 1
    
    MLLM2HB(i)=sum(log(normpdf(HBDiffs(:,1),ParamsHB(i,1),ParamsHB(i,2)))); % The maximum
    % log likelihood for Model 2
    
    BFHB(i)=MLLM1HB(i)-MLLM2HB(i)-.5*(1-2)*log(length(HBDiffs(:,i))); % The Schwarz criterion
    % says that the log of the Bayes Factor is the difference in the
    % maximum log likelihoods of the two models [MLLM1(i)-MLLM2(i)]
    % minus a penalty for complexity.
    % The penalty for complexity (extra free parameters) is:
    % -.5*(1-2)*log(length(Diffs(:,i))). Notice that it depends on both the
    % difference in the number of free parameters between the two models
    % [(1-2)] AND the log of the number of data being "explained" The
    % Bayes Factor is the factor by which the data alter the prior odds.
    % That is, the posterior odds (the odds in favor of one model or the
    % other AFTER the data have been considered) are the prior odds (the
    % odds BEFORE the data are considered) times the Bayes Factor. The
    % Bayes Factor is the antilog of the result returned by the calculation
    % specified in this command line.
    
end
    
BFHB % displaying the logs of the Bayes factors
% Result is:

% BFHB =
% 
%    -0.2417    0.8777    1.1934   -0.2196    0.4959    1.0589

% The Bayes Factor favors the no-difference model in some cases and the
% difference model in others, but the odds either way are low:

exp(BFHB) % This gives us the odds for each comparison (odds less than one
% are odds in favor of Model 2 (that there is a difference)

% ans =
% 
%     0.7853    2.4053    3.2984    0.8028    1.6419    2.8831

% Combining the odds across comparisons (which assumes independence):

exp(sum(BFHB))

% ans =
% 
%    23.6782

% gives strongish odds in favor of the no-difference model, but whether
% combining evidence across conditions in this way is legit or not depends
% on whether subjects' behavior in one condition is independent of their
% behavior in the other

%% Correlations

figure % 4th figure
subplot(2,1,1) % first of two plots. This one will look at the extent to
% which a subject's filled-hole-empty hole difference predicts their
% filled-bay-empty-bay differrence. It superposes scatter plots from the
% three levels, using color and symbol shape to distinguish them.

plot(Diffs(:,1),Diffs(:,4),'ko') % scatter plot of the FB-EB diff versus the
% FH-EH diff when there are only 2 elements in the displays
hold on
plot(Diffs(:,2),Diffs(:,5),'rs') % scatter plot of the FB-EB diff versus the
% FH-EH diff when there are only 2 elements in the displays

plot(Diffs(:,3),Diffs(:,6),'gv') % scatter plot of the FB-EB diff versus the
% FH-EH diff when there are only 2 elements in the displays

xlabel('Filled vs Empty Holes (FH-EH)')
ylabel('Filled vs Empty Bays (FB-EB)')

title('Correlation btwn Filled-Empty Differences for Holes and Bays')

legend('2-element displays','4-element displays','6-element displays')

% These correlations are not significant:

[R14,P14]=corrcoef(Diffs(:,[1 4]))
[R25,P25]=corrcoef(Diffs(:,[2 5]))
[R36,P36]=corrcoef(Diffs(:,[3 6]))
% which gives
% R14 =
% 
%     1.0000    0.2103
%     0.2103    1.0000
% 
% P14 =
% 
%     1.0000    0.3239
%     0.3239    1.0000
% 
% R25 =
% 
%     1.0000   -0.0352
%    -0.0352    1.0000
% 
% P25 =
% 
%     1.0000    0.8704
%     0.8704    1.0000
% 
% R36 =
% 
%     1.0000    0.0832
%     0.0832    1.0000
% 
% P36 =
% 
%     1.0000    0.6990
%     0.6990    1.0000

% Just look at the values in the upper right corner of each 2x2 array. The
% arrays beginning w 'R' give the correlation coeeficients; those beginning
% w 'P' give the significance of the correlation. The correlations are
% small and insignificant, which is good news.

subplot(2,1,2) % Now we look at the extent to which a subject's hole-bay
% difference in the filled condition predicts their hole-bay difference in
% the empty condition

plot(HBDiffs(:,1),HBDiffs(:,4),'ko') % scatter plot of the FB-EB diff versus the
% FH-EH diff when there are only 2 elements in the displays
hold on
plot(HBDiffs(:,2),HBDiffs(:,5),'rs') % scatter plot of the FB-EB diff versus the
% FH-EH diff when there are only 2 elements in the displays

plot(HBDiffs(:,3),HBDiffs(:,6),'gv') % scatter plot of the FB-EB diff versus the
% FH-EH diff when there are only 2 elements in the displays

xlabel('Filled Holes vs Filled Bays (FH-FB)')
ylabel('Empty Holes vs Empty Bays (EH-EB)')

title('Correlation btwn Hole-Bay Differences in Filled & Empty Conditions')

% 
[RHB14,PHB14]=corrcoef(HBDiffs(:,[1 4]))
[RHB25,PHB25]=corrcoef(HBDiffs(:,[2 5]))
[RHB36,PHB36]=corrcoef(HBDiffs(:,[3 6]))
% which gives
% RHB14 =
% 
%     1.0000    0.1300
%     0.1300    1.0000
% 
% PHB14 =
% 
%     1.0000    0.5449
%     0.5449    1.0000
% 
% RHB25 =
% 
%     1.0000    0.1123
%     0.1123    1.0000
% 
% PHB25 =
% 
%     1.0000    0.6013
%     0.6013    1.0000
% 
% RHB36 =
% 
%     1.0000    0.0092
%     0.0092    1.0000
% 
% PHB36 =
% 
%     1.0000    0.9659
%     0.9659    1.0000

% Again, the correlations are small and insignificant. This tends to
% legitimize our having combined the Bayes Factors from the different
% comparisons into one summary Bayes Factor

%% Running the ANOVA
RspdData=reshape(SungHoData,12*24,1); % Putting data in a single col vector

Factors(1,:)={['Subject'] ['F-E'] ['H-B'] ['Elements']};
% Cell array w factor indices
Factors{2,1}=repmat((1:24)',12,1); % Subject indices
Factors{2,2}=repmat([repmat('f',3*24,1);repmat('e',3*24,1)],2,1) % f-e
Factors{2,3}=[repmat('h',6*24,1);repmat('b',6*24,1)] % hole-bay
Factors{2,4}=repmat([2*ones(24,1);4*ones(24,1);6*ones(24,1)],4,1) % # elements

[P,T,Stats,Terms]=anovan(RspdData,Factors(2,:),'random',[1],'model',2,...
    'varnames',Factors(1,:)) % The n-way anova command. The first argument
% is the data vector; the 2nd is the cell array with the index vectors for
% the factors. Each cell contains the indices for one of the 4 factors
% (Subjects, F-E, H-B, #Elements). The third & fourth arguments tell Matlab
% that the Subject factor is random. The 4th & 5th arguments tell it to
% compute interactions only up to the 2nd order. The 5th & 6th arguments
% give it the names of the factors, which names it will use in the ANOVA
% tables it returns
% 
% % Result (from the ANOVA table, which appears in its own window):
% 
% Analysis of Variance
%   Source	 Sum Sq.	d.f.	Mean Sq.	  F	        Prob>F
%   Subject	0.21229	    23	  0.0092299	   1.9138	    0.22827
%   F-E	  0.00012535	1	  0.00012535	0.022969	0.88086
%   H-B	  0.058084	    1	  0.058084	   10.4268	    0.0037116
%   Elements	1.3524	    2	  0.67621	  118.0414	    0
%   Subject*F-E	0.12552	23	  0.0054572	    0.9146	    0.57917
%   Subject*H-B	0.12812	23	  0.0055706	   0.93361	    0.55411
%   Subject*Elements	0.26351	46	0.0057286	0.96008	    0.55063
%   F-E*H-B	0.0022781	1	  0.0022781	    0.3818	    0.5375
%   F-E*Elements	0.00059653	2	0.00029826	0.049988    0.95126
%   H-B*Elements	0.012517	2	0.0062587	1.0489	    0.35267
%   Error	  0.97258	   163	  0.0059668		
%   Total	  3.128	       287			
% Constrained (Type III) sums of squares.

% There is a hugely significant effect of the number of
% elements in the display (F=118) and a significant effect of holes vs bays
% but the effect of filled vs empty does not approach significance and none
% of the interactions is significant. The significant effect of holes vs
% bays is not congruent with what we learned from the relative likelihood
% analyses of the differences. It would be good to get to the bottom of
% this. One question concerns the fact that the within-cell variance is not
% constant across levels of the Number-of-Elements factor; as the mean A's
% get lower, their variability increases, and so, naturally, does the
% variability in their differences. By averaging the variability from the
% low-#-elements conditions together with the variability from the
% high-#-elements conditions, the ANOVA may come up with an a too-low
% estimate of the error variance against which the F for the largish
% difference in mean A's at the high-#-elements end is pitted. I don't THINK
% the ANOVA takes this into account, but I am not sure. Recall that we also
% did t-tests for whether for any given number of stimulus elements and a
% given value for filled vs empty, the hole-bay differences were
% significant. In no case were they, although some of these comparisons
% were nearly significant. All in all, I would not be inclined to make much
% of the hole-bay difference apparent in Figure 1.
