# empir_ARL2_profile.m

Final plot of the empirical Avg Response Levels w/ 56-trial blocks. See Figures 4 and 5 in ProtoPaper/expmt.tex. Compare with model_ARL2_predictions.m

    File: work/CMUExper/fdbk1/data/empir_ARL2_profile.m
    Date: 2007-08-12
    Alexander Petrov, http://alexpetrov.com

## Contents

- Load Fdbk1 human data
- Group descriptors
- Collapse across groups
- Group the points in pairs
- Save the ARL2 profiles for subsequent fits
- Plot the ARL profile. Same LineStyles as model_ARL2_predictions.
- Export ARL2 to SPSS in ASCII format
- Export same ARL2 to SSPS, rearranging for within-sbj ANOVA
- Define Assimilation := ARL(N)-ARL(P) and plot Assim profile
- Multiple Regression (ANCOVA) on the Assim group data
- Full model
- Model without interaction terms
- Significance of the Time covariate
- Significance of the Fdbk factor
- Significance of the Order factor
- Clean up

## Load Fdbk1 human data

    ARL is calculated as in the CatRat1 experiment
    See work/CMUExper/CatRat1/data/trend/

    clear all ;
    cd(fullfile(work_pathstr,'CMUExper','fdbk1','data')) ;
    load('S.mat') ;
    idx = find([S(:).sbj_number]>=99) ;   % eliminate pilot runs
    S(idx) = [] ;
    ARL = [S(:).ARL] ;    % 17x55

## Group descriptors

    gr=[S(:).group]' ;
    xtab1(gr)
    fdbk_first_p = [1 1 1 0 0 0] ;    % exp schedule by group
    %gr_descr = {'Ufb','Pfb','Nfb','xxx','Pno','Nno' } ;
    gr_descr = {'U1','L1','H1','xxx','L2','H2' } ;

       Value   Count  Percent  Cum_cnt  Cum_pct
    -------------------------------------------
           1      11    20.00       11    20.00
           2      11    20.00       22    40.00
           3      11    20.00       33    60.00
           5      11    20.00       44    80.00
           6      11    20.00       55   100.00
    -------------------------------------------

## Collapse across groups

    grARL = zeros(size(ARL,1),max(gr)) ;   % 17x6
    for k = 1:max(gr)
        idx = find(gr==k) ;
        if (~isempty(idx))
            grARL(:,k) = mean(ARL(:,idx),2) ;
        end
    end
    grARL

    grARL =
        4.3358    4.1354    4.2075         0    4.1789    4.2197
        4.0411    4.0838    4.1610         0    4.4570    4.3693
        4.0639    4.0350    4.1723         0    4.8891    4.0287
        4.1521    4.0537    4.2286         0    4.8043    4.1585
        4.2236    3.8352    4.3046         0    4.7330    4.0194
        4.3239    4.3174    4.4587         0    4.0544    4.2597
        4.5233    4.5207    4.4153         0    3.9155    4.4074
        4.3777    4.2732    4.4184         0    4.0007    4.3081
        4.3067    4.2592    4.4738         0    3.8678    4.3894
        4.3178    3.8628    4.4319         0    3.9486    4.2385
        4.1629    3.7666    4.3952         0    4.2545    4.3144
        4.0674    3.7753    4.3084         0    4.1976    4.4453
        4.1249    3.7122    4.4734         0    4.2139    4.4834
        4.2917    3.9530    4.2607         0    3.9207    4.3365
        4.2568    3.9123    4.4973         0    3.8637    4.4364
        4.3501    3.6305    4.5398         0    3.8915    4.3348
        4.3486    3.7241    4.5097         0    3.6673    4.4448

## Group the points in pairs

The difference between this and re-fitting the ARLs over 56-trial windows is less than 1%. I checked it on synthetic data. The problem is that data_by_sbj.m does it with 28-trial blocks and takes care of missing values. There is no need to redo this, given that the outcome will be virtually the same.

    idx = [2:2:16]' ;
    grARL2 = (grARL(idx,:)+grARL(idx+1,:))./2 ;
    grARL2 = [grARL(1,:) ; grARL2]               % 9x6

    ARL2 = (ARL(idx,:)+ARL(idx+1,:))./2 ;
    ARL2 = [ARL(1,:) ; ARL2] ;                   % 9x55

    grARL2 =
        4.3358    4.1354    4.2075         0    4.1789    4.2197
        4.0525    4.0594    4.1666         0    4.6730    4.1990
        4.1878    3.9444    4.2666         0    4.7687    4.0890
        4.4236    4.4190    4.4370         0    3.9850    4.3336
        4.3422    4.2662    4.4461         0    3.9343    4.3488
        4.2404    3.8147    4.4136         0    4.1016    4.2765
        4.0962    3.7438    4.3909         0    4.2058    4.4644
        4.2742    3.9327    4.3790         0    3.8922    4.3864
        4.3494    3.6773    4.5247         0    3.7794    4.3898

## Save the ARL2 profiles for subsequent fits

    S2 = S ;
    for k = 1:length(S)
        S2(k).ARL2 = ARL2(:,k) ;
    end
    save S2 S2

## Plot the ARL profile. Same LineStyles as model_ARL2_predictions.

    ls = {'r-' 'b-' 'k<-' 'r--' 'b--' 'k<--'} ;   % group line styles
    lw = [1 2 2  1 2 2] ;                         % group line widths
    %lw = [.5 1.5 1.5  .5 1.5 1.5] ;               % group line widths

    x2 = [14 , 56:56:448]' ;
    xtick = [28 , 140:112:476] ;
    ytick = [3.6:.2:4.8] ;
    ax = [0 476 3.6 4.8] ;

    clf ; hold on ;
    for k = unique(gr)'
        hh = plot(x2,grARL2(:,k),ls{k}) ;
        set(hh,'LineWidth',lw(k),'MarkerSize',4) ;
    end
    hold off ; box on ;
    axis(ax) ; set(gca,'xtick',xtick,'ytick',ytick,'xgrid','on') ;
    hh = xlabel('Trial') ; set(hh,'FontSize',14) ;
    hh = ylabel('Average response level') ; set(hh,'FontSize',14) ;
    hh = legend('U1','L1','H1','L2','H2',3) ; set(hh,'FontSize',14) ;
    % Saved as work/papers/ProtoPaper/fig/ARL-empir.pdf

## Export ARL2 to SPSS in ASCII format

Write a text file with one row per data point and the following columns:

- 1. Sbj -- 1-55
- 2. Group -- 1=U1, 2=L1, 3=H1, 5=L2, 6=H2
- 3. Context -- 1=U, 2=L, 3=H
- 4. Order -- 1=Feedback_first, 2=No_feedback_first
- 5. Block -- 1-9. Each is 56 trials (except the first, which is 28).
- 6. Feedback -- 0=None, 1=Fdbk
- 7. ARL2

    % Design matrix for each group
    Context = [1 2 3 1 2 3] ;
    Order = [1 1 1 2 2 2] ;
    Feedback = [1 1 1 0 0 1 1 0 0 ;  ...    % when order=1
                1 0 0 1 1 0 0 1 1 ]' ;      % when order=2
    Block = [1:9]' ;
    Design = zeros([9 5 6]) ;     % 9 blocks x 5 variables x 6 groups
    D = zeros(9,5) ;
    for k = 1:6
        D(:,1) = k ;
        D(:,2) = Context(k) ;
        D(:,3) = Order(k) ;
        D(:,4) = Block ;
        D(:,5) = Feedback(:,Order(k)) ;
        Design(:,:,k) = D ;
    end

    % Data matrix with one row per point
    data1 = zeros(0,7) ;
    D = zeros(9,7) ;
    for k = 1:length(S) ;
        D(:,1) = S(k).sbj_number ;
        D(:,2:6) = Design(:,:,S(k).group) ;
        D(:,7) = ARL2(:,k) ;
        data1 = [data1 ; D] ;
    end

    % Write to an ASCII file
    fid = fopen('ARL2-bn-Ss.dat','w') ;
    fprintf(fid,'Sbj Group Context Order Block Feedback ARL2\n') ;
    for k = 1:size(data1,1)
        fprintf(fid,'%02d %1d %1d %1d %02d %1d %5.3f\n',data1(k,:)) ;
    end
    fclose(fid) ;

## Export same ARL2 to SSPS, rearranging for within-sbj ANOVA

Write a text file with one row per subject and the following columns:

- 1. Sbj -- 1-55
- 2. Group -- 1=U1, 2=L1, 3=H1, 5=L2, 6=H2
- 3. Context -- 1=U, 2=L, 3=H
- 4. Order -- 1=Feedback_first, 2=No_feedback_first
- 5. W -- ARL2 value for the "warm-up" initial block
- 6-9. F1-F4 -- ARL2 values during the four feedback blocks
- 10-13. N1-N4 -- ARL2 values during the four no-feedback blocks

Note that this data format throws away the temporal order of observations. They are rearranged depending on Order so that in the resulting ASCII file all feedback ARLs are listed in columns F1-F4 and all no-feedback ARLs are listed in columns N1-N4.

Results of SPSS analyses available in work/CMUExper/fdbk1/data/SPSS/

    rearrange = [1, 2 3 6 7, 4 5 8 9 ; ...     % when order=1
                 1, 4 5 8 9, 2 3 6 7 ]' ;      % when order=2

    data2 = zeros(length(S),13) ;
    for k = 1:length(S)
        data2(k,1) = S(k).sbj_number ;
        hh = S(k).group ;
        data2(k,2:4) = Design(1,1:3,hh) ;
        data2(k,5:13) = ARL2(rearrange(:,Order(hh)),k)' ;
    end

    % Write to an ASCII file
    fid = fopen('ARL2-wn-Ss.dat','w') ;
    fprintf(fid,'Sbj Group Context Order W F1 F2 F3 F4 N1 N2 N3 N4\n') ;
    for k = 1:size(data2,1)
        fprintf(fid,...
          '%02d %1d %1d %1d %5.3f %5.3f %5.3f %5.3f %5.3f %5.3f %5.3f %5.3f %5.3f\n',...
          data2(k,:)) ;
    end
    fclose(fid) ;

## Define Assimilation := ARL(N)-ARL(P) and plot Assim profile

See CMUExper/fdbk1/data/assimil_profiles.m

- 1 = Feedback first = Group 3 - Group 2
- 2 = No-feedback first = Group 6 - Group 5

    Assim = grARL2(:,[3 6]) - grARL2(:,[2 5])   % <-- Sic!

    ax = [0 476 -0.85 1.15] ;
    ytick1 = [-.8:.2:1] ;
    idx1 = [1 2 3 6 7] ;    % feedback blocks for the feedback-first groups
    idx2 = [1 4 5 8 9] ;    % feedback blocks for the no-feedback-first grps

    clf ;
    hh = plot(x2,Assim(:,1),'r-',x2,Assim(:,2),'b--',...
              x2(idx1),Assim(idx1,1),'ro',x2(idx2),Assim(idx2,2),'bs') ;
        set(hh,'LineWidth',2) ;
        axis(ax) ;
        hh=refline(0,0) ;
        set(hh,'Color','k','LineStyle','-','LineWidth',1) ;
    box on ;
    axis(ax) ; set(gca,'xtick',xtick,'ytick',ytick1,'xgrid','on') ;
    hh = xlabel('Trial') ; set(hh,'FontSize',14) ;
    hh = ylabel('ARL(H) - ARL(L)') ; set(hh,'FontSize',14) ;
    hh = legend('Feedback first','No-feedback first',4) ; set(hh,'FontSize',14) ;
    % Saved as work/papers/ProtoPaper/fig/Assim-empir.pdf

    Assim =
        0.0720    0.0408
        0.1073   -0.4740
        0.3222   -0.6797
        0.0180    0.3486
        0.1799    0.4145
        0.5989    0.1749
        0.6471    0.2586
        0.4463    0.4943
        0.8474    0.6104

## Multiple Regression (ANCOVA) on the Assim group data

Ignore the "warm-up" block 1. Regress the remaining 16 Assim points w.r.t the following predictors:

- 1. Time: 2-9
- 2. Feedback: -1=None, +1=Fdbk
- 3. Order: -1=No-fdbk-first, +1=Feedback-first
- 4. Fdbk X Order

Cf. David Howell (1992). "Statistical Methods for Psychology" (3rd Ed.). Chap. 16: ANOVA and ANCOVA as general linear models. Design matrix, p.553

    y = reshape(Assim(2:9,:),16,1) ;    % dependent variable, N=16

    X_Time = [2:9 , 2:9]' ;
    X_Fdbk = [+1 +1 -1 -1 +1 +1 -1 -1 , -1 -1 +1 +1 -1 -1 +1 +1]' ;
    X_Order = [ones(8,1) ; -ones(8,1)] ;
    X_FxO = X_Fdbk .* X_Order ;

    % Design matrix X
    X = ones(16,1) ;                                   % Regression constant
    X = [X (X_Time - mean(X_Time))./std(X_Time,1)] ;   % Standardize X_Time
    X = [X X_Fdbk X_Order X_FxO]                       % Already standardized

    X =
        1.0000   -1.5275    1.0000    1.0000    1.0000
        1.0000   -1.0911    1.0000    1.0000    1.0000
        1.0000   -0.6547   -1.0000    1.0000   -1.0000
        1.0000   -0.2182   -1.0000    1.0000   -1.0000
        1.0000    0.2182    1.0000    1.0000    1.0000
        1.0000    0.6547    1.0000    1.0000    1.0000
        1.0000    1.0911   -1.0000    1.0000   -1.0000
        1.0000    1.5275   -1.0000    1.0000   -1.0000
        1.0000   -1.5275   -1.0000   -1.0000    1.0000
        1.0000   -1.0911   -1.0000   -1.0000    1.0000
        1.0000   -0.6547    1.0000   -1.0000   -1.0000
        1.0000   -0.2182    1.0000   -1.0000   -1.0000
        1.0000    0.2182   -1.0000   -1.0000    1.0000
        1.0000    0.6547   -1.0000   -1.0000    1.0000
        1.0000    1.0911    1.0000   -1.0000   -1.0000
        1.0000    1.5275    1.0000   -1.0000   -1.0000

## Full model

Use REGRESS in the STATS toolbox to do multiple regression. The regression coeff of Fx0 is not stat. diff. from 0:

    [b_full,bint_full,r,rint,stats_full] = regress(y,X) ;
    [b_full bint_full]'    % [Const Time Fdbk Order FxO]
    stats_full             % [R2  F  p  MSE]

    ans =
        0.2697    0.2738    0.1732    0.1262   -0.0308
        0.1778    0.1717    0.0814    0.0344   -0.1328
        0.3615    0.3759    0.2651    0.2180    0.0713
    stats_full =
        0.8710   18.5597    0.0001    0.0278

## Model without interaction terms

    [b_main,bint_main,r,rint,stats_main] = regress(y,X(:,1:4)) ;
    [b_main bint_main]'    % [Const Time Fdbk Order]
    stats_main             % [R2  F  p  MSE]
    rmse = sqrt(stats_main(4))

    ans =
        0.2697    0.2872    0.1732    0.1262
        0.1809    0.1985    0.0845    0.0375
        0.3584    0.3760    0.2620    0.2150
    stats_main =
        0.8658   25.8036    0.0000    0.0265
    rmse =
        0.1629

## Significance of the Time covariate

    [b,bint,r,rint,stats] = regress(y,X(:,[1 3 4])) ;
    R2_noTime = stats(1)
    [F_Time,p] = R2signif(16,stats_main(1),3,R2_noTime,2)

    R2_noTime =
        0.3097
    Warning: Could not find an exact (case-sensitive) match for 'Fcdf'.
    /Applications/MATLAB74/toolbox/stats/fcdf.m is a case-insensitive match and will be used instead.
    You can improve the performance of your code by using exact
    name matches and we therefore recommend that you update your
    usage accordingly. Alternatively, you can disable this warning using
    warning('off','MATLAB:dispatcher:InexactMatch').
    F_Time =
       49.7215
    p =
       1.3350e-05

## Significance of the Fdbk factor

    [b,bint,r,rint,stats] = regress(y,X(:,[1 2 4])) ;
    R2_noFdbk = stats(1)
    [F_Fdbk,p] = R2signif(16,stats_main(1),3,R2_noFdbk,2)

    R2_noFdbk =
        0.6635
    F_Fdbk =
       18.0880
    p =
        0.0011

## Significance of the Order factor

    [b,bint,r,rint,stats] = regress(y,X(:,[1 2 3])) ;
    R2_noOrder = stats(1)
    [F_Order,p] = R2signif(16,stats_main(1),3,R2_noOrder,2)

    R2_noOrder =
        0.7584
    F_Order =
        9.6011
    p =
        0.0092

## Clean up

    clear hh k idx D b bint r rint stats ;