% script m-file StudentDataAnalyses.m. Analyses of data provided by
% students in the Advanced Statistics course

%% Natalie Obrect's data

NatalieData=1; % Creating a variable into which we then paste the data
% from the Excel spreadsheet Natalie provided. There are 11 columns of
% data. Natalie provided the following guide to those data:

% DEFINING DATA FILE COLUMNS:
% 
% 1 = condition (1 = control, 2 = conexp, 3 = exp)
% 2 = order (1 or 2)
% 3 = subject number (1-337)
% 
% /*
% The condition numbers are: 
% 1)+++, 2)++-, 3)+-+, 4)-++, 5)--+, 6)-+-, 7)+--, 8)---

% Wxperiment: Each subject made a confidence judgment about 8 cases,
% that is, eight different comparisons of consumer ratings of a pair of
% products. What varied from comparison to comparison was the power level
% of the sample size (small sample vs large sample), the power level of the
% mean difference (large difference between the means, small difference)
% and the power level of the variability (low variability w/i each sample
% vs high variability).
% 
% where position 1 = sample size, position 2 = mean difference, and postion 3 = std dev.
% and + = high power and - = low power
% 
% order 1: 7)+--, 4)-++, 2)++-, 8)---, 3)+-+, 1)+++, 6)-+-, 5)--+ 
% order 2: 2)++-, 3)+-+, 8)---, 1)+++, 7)+--, 5)--+, 4)-++, 6)-+-

% +-- to means that there was a large sample size, but a
% small difference between the means and large variability.  -++
% means that the sample sizes were small, but there was a large
% difference between the means and low within-sample variabilty. And so on

% The numbers in Columns 4-11 of the data array are the confidence ratings
% that the subject gave for each of the 8 comparisons she was asked to make
% */
% 
% 4 = cond1 (+++)
% 5 = cond2 (++-)
% 6 = cond3 (+-+)
% 7 = cond4 (-++)
% 8 = cond5 (--+)
% 9 = cond6 (-+-)
% 10= cond7 (+--)
% 11= cond8 (---)
% 
% 
% WHAT WE WANT TO KNOW:
% 
% Ss are asked to compare pairs of products that were rated by consumers and decide
% how confident they are that the product with the higher mean rating is really
% better than the comparison product.  The product data are set up so that if Ss
% judgments are in line with statistical power/t-statistics/LRs they should treat 
% N ,mean diff, and sd equally.
% 
% Some Ss were: 
% -only given statistical summaries of the data (N, mean, sd for each product);
% this is called the control or stat-only group
% -others were given this statistical information plus the corresponding raw data; 
% there are the conexp or stat+data group and finally
% -the last group were given only the raw data; called exp or data-only group
% 
% We want to know 1) how subjects weight these 3 factors and 2) whether the way that
% the data are provided affect how these factors are used.
%               *****************************************

NatalieData(1,:)=[]; % First row were column numbers; this deletes that row
% from the array, leaving only the data

% To do an analysis of variance (and for other purposes as well), the most
% convenient way to arrange the data is to have it all in a single column,
% with other columns of the same array specifying the conditions under
% which the datum in that row was obtained. Each of these
% condition-specifying columns corresponds to what is called a factor. The
% column contains different identifying numbers (indices) for each
% different level of that factor. In Natalie's case, we have the following
% six factors (with the number of levels of each factor indicated in
% parentheses):
%   Condition (3); 'Stat Only', 'Stat+Data', 'Data Only'

%   Order (2); this is a "nuisance variable" meaning that we are not
%   interested in its effects & have counterbalanced order so as to
%   distribute whatever effects this variable may have across the cells we
%   are interested in
%
%   Subject (337); indexed by number

%   N (number in samples) (2); 'L', 'S'

%   M (the size of the difference btwn the means) (2); 'L', 'S'

%   S (for sigma, the within-sample variability) (2); 'L','S'



% However, the data are organized with the levels of the first three
% factors varying from row to row, while the levels of the last three
% factors vary from column to column (Columns 4-11). So we begin by
% rearranging the data array to get all the data in a single column:
%%
NatDat=reshape(NatalieData(:,4:11),8*337,1); % Creates the single data
% column 2696 rows long. Now, we need to make appropriately indexed
% columns for each of the factors

Icols=repmat(NatalieData(:,1:3),8,1); % The indices in NatalieData form
% a 337x3 matrix. This replicates that matrix 8 times and places the 8
% replications one below the other to create a matrix 2696x3 matrix of
% indices for the Condition, Order and Subject factors. Now we need to
% create index columns for the three other factors (Mean, Variability & N).
% Columns 4-6 of the original matrix contained data from the
% large-mean-difference condition; all other columns were from the
% low-mean-difference level
Icols([1:1011 2023:2359],4)=1; % This creates a 4th column in Icols and puts 1's in the
% the rows where there was a large N . Matlab fills in 0's in the remaining
% rows.
Icols([1:674 1012:1348 1686:2022],5)=1; % Creates a
% 5th column and fills in 1s in the rows where there was a large mean diff.
% Matlab fills in 0s in all the other rows

Icols([1:337 675:1685],6)=1; % Creates a 6th column in Icols and puts 1s in the
% rows where there was a small std dev

Icols(:,7)=sort(repmat([1;2;3;4;5;6;7;8],337,1)); % This adds a 7th index column
% corresponding to the 8 data columns in the original array. These indices
% are redundant in that they could be inferred from the pattern in the
% three index columns to the left. They will, however, prove useful later
% on when we want to do operations subtracting one of those original
% columns from another

% Now we create the cell array that identifies our factors and their
% levels. The first row of the cell matrix gives the name of the factor.
% This is for our benefit only. It will not be used in the computations.
% The second row of the cell array gives in each cell either a numerical or
% a character vector 2696 rows long, giving an index (level) for each datum
% for that factor. In other words, identifying the level of the factor at
% which the datum was obtained
%%
Factors{1,1}='Form of Data';
Factors{1,2}='Order';
Factors{1,3}='Subject';
Factors{1,4}='N';
Factors{1,5}='MeanDiff';
Factors{1,6}='Sigma';

Factors{2,1}(find(Icols(:,1)==1),1)='S';
Factors{2,1}(find(Icols(:,1)==2),1)='B';
Factors{2,1}(find(Icols(:,1)==3),1)='D';
% These commands create a cell with a character vector 2696 rows long, with
% each row containing an 'S' if it was a Stats-Only level of the Form
% factor, a 'B' if it was both stats & data, and a 'D' if it was only the
% data
Factors{2,2}=Icols(:,2); % Cell containing the order indices (1 or 2)

Factors{2,3}=Icols(:,3); % Cell containing subject-number vector

Factors{2,4}=Icols(:,4); % Cell with 1's whereever power of N was high

Factors{2,5}=Icols(:,5); % Cell with 1's whereever power of mean diff high

Factors{2,6}=Icols(:,6); % Cell with 1's wherever power of sigma is high
% (i.e., variability is low)
%%
% Now we are ready to start looking at the data graphically. This is the
% most important step in the analysis. You should thoroughly understand
% your data by means of graphic visualization BEFORE attempting a
% statistical analysis.

% First, let's look at whether the power has any effect. There are three
% factors contributing to the power: N (size of samples), M (size of mean
% differences) and Sigma (variability). The indices of these columns are in
% columns 4-6 of Icols. There are 1's wherever the level of one of those
% factors was such as to increase the power. Therefore, the sum of those
% three columns gives three levels of power.
%%
[Means,SEM,Counts,Gname] = grpstats(NatDat,sum(Icols(:,4:6),2),.05);

% The first argument of grpstats is the data vector, the second argument is
% the vector that determines the breakdown of the data into groups. In this
% case, we are grouping the data on the basis of the sum of the power
% factors. This sum will have four different values (0, 1, 2, 3), so there
% will be four groups

xlabel('Total Power')
ylabel('Mean Confidence Ratings (w. 95% conf intervals)')
ylim([2.2 3.6])

% This looks reasonable; as we increase the total power (put more and more
% relevant factors at their high-power levels, subjects' confidence grows

% Now, let's break this down by the Form factor, that is, look at the
% relation at each of the three different levels of the Form factor

%%
figure % new figure window-- figure(2)

Rows1=find(Icols(:,1)==1); % The rows where subjects saw only the stats

Rows2=find(Icols(:,1)==2); % The rows where subjects saw stats & data
Rows3=find(Icols(:,1)==3); % The rows where subjects saw only data

[Means1,SEM1,Counts1,Gname] = grpstats(NatDat(Rows1),sum(Icols(Rows1,4:6),2),.05);
hold on % & used editor to color black
[Means2,SEM2,Counts2,Gname] = grpstats(NatDat(Rows2),sum(Icols(Rows2,4:6),2),.05);
% use editor to color this curve red

[Means3,SEM3,Counts3,Gname] = grpstats(NatDat(Rows3),sum(Icols(Rows3,4:6),2),.05);
% used editor to color this one green
xlabel('Total Power')

ylabel('Confidence')

legend('Stats Only','Stats&Data','Data Only')
% The u-shaped dependence of confidence on total power appears with all
% three forms of presentation ('Stats Only', 'Stats & Data' and 'Data
% Only'). but it is most pronounced in the Data Only case, and the bottom
% of the U shifts to Power Level 2 in that case
%%

%%
% Now, let's look at the effects of the different stats (N, M and S)
figure  % third figure

[MeansN,SEMN,CountsN,Gname] = grpstats(NatDat,Icols(:,4),.05);
% Col 4 of Icols indexes the two levels of N. The high levels come from
% Columns 4, 5, 6 and 10 of the original data array
hold on  % leave this colored blue
%%
[MeansM,SEMM,CountsM,Gname] = grpstats(NatDat,Icols(:,5),.05);
% Col 5 of Icols indexes the two levels of M, the mean difference. The high
% levels come from Cols 4, 5, 6 & 9 of the original data array (use editor
% to color this curve brown)
%%
[MeansS,SEMS,CountsS,Gname] = grpstats(NatDat,Icols(:,6),.05);
% Col 6 of Icols indexes the two levels of S, the variability. The
% high-power (low variability levels) come from Cols 4, 6, 7 and 8 of the
% original data array

legend('N','M','S')
xlabel('Power')
ylabel('Confidence')
title('Effects of the Power of N, M and S')

% We have a classic crossing interaction: The effect of power depends on
% which parameter one manipulates: increasing power by increasing N has no
% effect, whereas increasing power by increasing M (mean diff) has a large
% effect and increasing it by decreasing variability(S) has a more modest
% effect. Crossing interactions are prized because they cannot be caused by
% non-linear scaling of the measurement axis (to be explained in class)


%%
% Now we further break these down by Form (looking for 3rd order
% interactions)

figure % 4th figure

subplot(3,1,1) % a separate plot for each of the three stats

[MeansN1,SEMN1,CountsN1,Gname] = grpstats(NatDat(Rows1),Icols(Rows1,4),.05);
% the "Rows1 in "NatDat(Rows1) limits the data being considered to those in
% which the presentation form was Stats Only. The Icols(Rows1,4) specifies
% the indices in the same rows of the column containing the indices that
% distinguish the levels of N.

% leave this curve colored blue
hold on
%%
[MeansN2,SEMN2,CountsN2,Gname] = grpstats(NatDat(Rows2),Icols(Rows2,4),.05);
%%
[MeansN3,SEMN3,CountsN3,Gname] = grpstats(NatDat(Rows3),Icols(Rows3,4),.05);
title('Effect of N, under different forms of presentation')
ylabel('Confidence')
%%
legend('Stats Only','Stats&Data','Data Only')
xlabel('Power')

% We have a 2d-order interaction: the effect of increasing N's power,
% depends on the Form of presentation: increasing N increases confidence
% (moderately) when Form is Stats Only, but (oddly) decreases it when Form
% is Data-Only. Adding Data to Stats eliminates positive effect of
% increasing power by increasing N
%%
subplot(3,1,2) % a separate plot for each of the three stats

[MeansM1,SEMM1,CountsM1,Gname] = grpstats(NatDat(Rows1),Icols(Rows1,5),.05);
% NatDat(Rows1) is the data from the stats-only condition. Icols(Rows1,5)
% specifies the indices in the same rows of Col 5 of Icols, which is the
% column that has the indices for the M factor. See farther below for a
% check on the values returned into MeansM1. Something is wrong
% leave blue
hold on
%%
[MeansM2,SEMM2,CountsM2,Gname] = grpstats(NatDat(Rows2),Icols(Rows2,5),.05);
%%
[MeansM3,SEMM3,CountsM3,Gname] = grpstats(NatDat(Rows3),Icols(Rows3,5),.05);

title('Effect of Mean Difference, under different forms of presentation')

ylabel('Confidence')

% Now we see a 3rd order interaction, because the 2nd-order interaction in
% the top panel goes away when we increase power by increasing the
% difference between the means

%%
subplot(3,1,3) % a separate plot for each of the three stats

[MeansM1,SEMM1,CountsM1,Gname] = grpstats(NatDat(Rows1),Icols(Rows1,6),.05);

hold on
%%
[MeansM2,SEMM2,CountsM2,Gname] = grpstats(NatDat(Rows2),Icols(Rows2,6),.05);
%%
[MeansM3,SEMM3,CountsM3,Gname] = grpstats(NatDat(Rows3),Icols(Rows3,6),.05);

title('Effect of Variability, under different forms of presentation')

ylabel('Confidence')

xlabel('Low-High Power')
%%
subplot(3,1,1);ylim([2 3.5]) % Making vert scales all the same

subplot(3,1,2);ylim([2 3.5])

subplot(3,1,3);ylim([2 3.5])

% More 3rd-order intervaction: the 2nd-order interaction seen in the top
% panel reverses in the bottom panel: the Data-Only Form of presentation
% endows the manipulation of power of S with a positive effect; whereas it
% had a negative effect when the parameter varied was N
%%
RowsStsOnly=find(NatalieData(:,1)==1);

RowsStsData=find(NatalieData(:,1)==2);

RowsDatOnly=find(NatalieData(:,1)==3);
%%


% Summary of what has been revealed by graphic examination of the data
% There are clear interactions, so main effects must be viewed with
% caution. The effects of increasing the power of one of the parameters
% depend on the Form in which the subject observes these parameters and on
% which parameter is being varied. Increasing the power of a parameter
% (e.g., N) may increase confidence under one Form of presentation and
% decrease it under another (top panel, fig 4) but with another parameter,
% this interaction may be absent or even reversed (bottom panel of fig 4).

% We might hazard an ANOVA at this point, looking only at interactions of
% 3rd order and less (the higher the order of the interactions, the more of
% them there are, the more computationally burdensome the ANOVA, and the
% harder it is to make sense of the output.) This is a complex ANOVA. There
% are 6 factors:
% From(3 levels)xSubject(337)xOrder(2)xN(2)xM(2)xS(2). Subject is a random
% factor, that is, the subjects were chosen more or less randomly from a
% much larger pool of potential subjects, and we are very unlikely to rerun
% those particular subjects in any subsequent replication. The other levels
% are fixed. Subjects are nested within levels of Form, that is, each level
% of Form has its own subjects. The observations at 2 different levels of
% the N, M and S factors are, on the other hand, repeated within each
% subject. It may have an effect but we do not want to think about what
% that might mean. We have counterbalanced order across subjects within
% each level of Form so as to spread its effects across the levels of the
% factors we are interested in. Thus, we can drop the Order factor, so now
% we have only a 5-factor ANOVA. Still intimidating.

% To use anovan, we have to have our indices in cell arrays: the cell array
% must have one cell for each factor and that cell must contain an index vector
% equal in length to the data vector. Each row of an index vector specifies
% the level of that factor for the datum in the corresponding row of the
% data vector. Our index vectors are the columns of Icols. We created the
% requisite cell array way back at the beginning, but now we drop the cell
% for the Order factor
Factors(:,2)=[]

[P,T,Stats,Terms]=anovan(NatDat,Factors(2,:),'random',[2],'model',3,...
    'varnames',{['Form'] ['Subjects'] ['N'] ['M'] ['S']})
% The first argument (NatDat) is the data; the second (Factors(2,:)) is the cell array
% with the index variables.  The
% next two arguments--'random' & [2]--indicate that Subjects is a random factor.
% The next two
% arguments--'model' & 3--specify that our model includes 3rd-order
% interactions. The final pair of
% arguments--'varnames' & {['Form'] ['Subjects'] ['N'] ['M'] ['S']}--causes
% the output to use our names in referring to the factors.

% This command generates an out-of-memory error after about a minute. I
% told you this was a complex ANOVA! My laptop has 1 gigabyte of RAM! This
% is the only time I have encountered an out-of-memory error.

%%
% Since we cannot run the omnibus F at the moment, let's look for
% interactions between the three power factors, N, M and S. Does the effect
% of N depend on whether M and/or S is set at a high or low level of power?

% This question is most easily approached using the original data format in
% that, it involves simply subtracting columns. For example, if we subtract
% the 4th data column (-++) from the 1st (+++), we get the effect of the
% change in the power of N when the level of the other two power factors is
% high. If we now subtract the 8th column (---) from the 7th (+--), we get
% the effect when the other two are low. Of course, we are seeing this
% effect under three different levels of the Form factor, so we will need
% to see wether that  matters. (One begins to see that there are a very
% large number of higher order interactions, which is why the full ANOVA is
% so computationally burdensome.)
%%
% Now let's look at whether the levels of the other parameter factors
% affect the effect of manipulating the power of one parameter.
% Conceptually, we want to compare different columns in the original data
% format, because in the last 8 columns of that original array, the
% combinations of the different parameter levels were varied
% systematically. It was to make this easy that we added a 7th column to
% Icols, the array of indices. That column will allow us to pluck out the
% rows that were in a given column. The first of the 8 columns was +++ (all
% parameters high). The 4th of the 8 columns was -++; it differed from the
% first, only in the level of N; the other two factors were both high in
% both columns. In creating meaningful variable names for what follows, I
% use 'p' for '+' and 'n' for '-' and R for rows.
%%
pppR=find(Icols(:,7)==1); % This variable now contains
% the numbers of the rows in the data vector that contain the data from
% Natalie's +++ column

ppnR=find(Icols(:,7)==2); % ++- column

ppnR=find(Icols(:,7)==2); % etc

pnpR=find(Icols(:,7)==3);

nppR=find(Icols(:,7)==4);

nnpR=find(Icols(:,7)==5);

npnR=find(Icols(:,7)==6);

pnnR=find(Icols(:,7)==7);

nnnR=find(Icols(:,7)==8);
%%
R1=find(Icols(:,1)==1); % row numbers of StatsOnly dat
R2=find(Icols(:,1)==2);
R3=find(Icols(:,1)==3);

% The intersection of, for example, pppR and R1, is a vector containing the
% numbers of the data rows that were in the +++ column and at the Stats-Only
% level of the Form factor. Thus, these intersections give the greatest
% possible breakdown of the data for the factors of interest (Form and the
% 3 power parameters. They break the data down to the cell level. In
% drawing contrasts between such cells, e.g, R1ppp vs R1npp, (varying the
% power of N when both other power factors are high and subjects saw stats
% only, one wants to look at the differences between data that came from
% the same subject, thereby taking out the between-subject variance
%%

figure % Fifth figure
subplot(4,3,1) % this figure will have 12 plots arranged in 3 columns of
% 4 plots each

plot([1 2],[mean(NatDat(Icols(find(Icols(nppR,4)==0)))) ...
    mean(NatDat(Icols(find(Icols(pppR,4)==1))))],'ko-')

xlim([.5 2.5])

hold on

plot([1 2],[mean(NatDat(Icols(find(Icols(nnnR,4)==0)))) ...
    mean(NatDat(Icols(find(Icols(pnnR,4)==1))))],'rx-')