MATLAB Examples

Contents

function regress_plot_twocategories_allslope(input_x1,input_y1,input_x2,input_y2,desc,x_label,y_label)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Name - regress_plot_twocategories_allslope
% Creation Date - 18th Oct 2011
% Author: Soumya Banerjee
% Website: https://sites.google.com/site/neelsoumya/
%
% Acknowledgements: Dr. Melanie Moses
%
% Description: Takes four column vectors, a description, x label and y label and
% plots the data and outputs all the statistics (r2, OLS slope, RMA slope, 95% CI
% intervals)
% Same syntax as plot(x,y) : x-axis data 1st followed by y-axis data
% Plots the slopes of both the categories as well as combined
%
% input_x1: category 1 column vector x
% input_y1: category 1 column vector y
% input_x2: category 2 column vector x
% input_y2: category 2 column vector y
% desc:     description of plot (goes in the title of the plot)
% x_label:  label for x-axis of plot
% y_label:  label for y-axis of plot
%
% Example usage:
% regress_plot_twocategories_allslope([1 1 1 1]', [1 1.1 1.2 1.3]', [2 2 2 2 2]',[2.2 2.1 2.3 2 2.4]',...
% 'analysis of all individuals mcmc burn-in thinning','log_1_0(M)','log_1_0(delta)')
%
% License - BSD
%
% Change History -
%                   18th Oct 2011 - Creation by Soumya Banerjee
%                   5th  Jan 2013 - Modified by Soumya Banerjee
%                                       Plots slopes for both categories
%                   2nd June 2014 - Modified by Soumya Banerjee
%                                       Automaticaly combines two
%                                       categories of data (fewer input
%                                       parameters needed)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

combined

input_x = [input_x1; input_x2];
input_y = [input_y1; input_y2];

msapop_ones = [ones(size(input_x)), input_x];
[b1 bint r rint stats] = regress(input_y,msapop_ones);
display('RMA Slope')
rmaslope = b1(2,1) / stats(1,1) ^ 0.5
display('bint')
bint
display('95% CI on intercept')
bint(1,:)
display('95% CI on slope')
bint(2,:)
display('r2')
stats(1,1)
display('OLS Slope')
b1(2,1)
olsslope = b1(2,1);
display('p value')
stats(1,3)
display('b1')
b1
RMA Slope

rmaslope =

    1.0817

bint

bint =

   -0.2792    0.4792
    0.8178    1.2822

95% CI on intercept

ans =

   -0.2792    0.4792

95% CI on slope

ans =

    0.8178    1.2822

r2

ans =

    0.9423

OLS Slope

ans =

    1.0500

p value

ans =

   1.3735e-05

b1

b1 =

    0.1000
    1.0500

first category

msapop_ones_1 = [ones(size(input_x1)), input_x1];
[b1_1 bint_1 r_1 rint_1 stats_1] = regress(input_y1,msapop_ones_1);
display('RMA Slope category 1')
rmaslope_1 = b1_1(2,1) / stats_1(1,1) ^ 0.5
display('bint category 1')
bint_1
display('95% CI on intercept category 1')
bint_1(1,:)
display('95% CI on slope category 1')
bint_1(2,:)
display('r2 category 1')
stats_1(1,1)
display('OLS Slope category 1')
b1_1(2,1)
olsslope_1 = b1_1(2,1);
display('p value category 1')
stats_1(1,3)
display('b1 category 1')
b1_1
Warning: X is rank deficient to within machine precision. 
RMA Slope category 1

rmaslope_1 =

   NaN

bint category 1

bint_1 =

    0.9446    1.3554
         0         0

95% CI on intercept category 1

ans =

    0.9446    1.3554

95% CI on slope category 1

ans =

     0     0

r2 category 1

ans =

     0

OLS Slope category 1

ans =

     0

p value category 1

ans =

   NaN

b1 category 1

b1_1 =

    1.1500
         0

second category

msapop_ones_2 = [ones(size(input_x2)), input_x2];
[b1_2 bint_2 r_2 rint_2 stats_2] = regress(input_y2,msapop_ones_2);
display('RMA Slope category 2')
rmaslope_2 = b1_2(2,1) / stats_2(1,1) ^ 0.5
display('bint category 2')
bint_2
display('95% CI on intercept category 2')
bint_2(1,:)
display('95% CI on slope category 2')
bint_2(2,:)
display('r2 category 2')
stats_2(1,1)
display('OLS Slope category 2')
b1_2(2,1)
olsslope_2 = b1_2(2,1);
display('p value category 2')
stats_2(1,3)
display('b1 category 2')
b1_2
Warning: X is rank deficient to within machine precision. 
RMA Slope category 2

rmaslope_2 =

   Inf

bint category 2

bint_2 =

         0         0
    1.0018    1.1982

95% CI on intercept category 2

ans =

     0     0

95% CI on slope category 2

ans =

    1.0018    1.1982

r2 category 2

ans =

     0

OLS Slope category 2

ans =

    1.1000

p value category 2

ans =

   NaN

b1 category 2

b1_2 =

         0
    1.1000

plot

figure
plot(input_x1,input_y1,'ks','MarkerEdgeColor','k','MarkerFaceColor','k','MarkerSize',10)
xlabel(x_label,'fontsize',23);
ylabel(y_label,'fontsize',23);
title(strcat(strcat(strcat(strcat(desc,strcat(strcat(strcat(strcat('r^2 value = ',num2str(stats(1,1))), strcat('OLSslope = ',num2str(olsslope))), 'RMAslope = '),num2str(rmaslope)))),'p value = '),num2str(stats(1,3))),'fontsize',18);
hold on
plot(input_x2,input_y2,'ro','MarkerEdgeColor','k','MarkerFaceColor','r','MarkerSize',10)
%  axis([-2 0.55 0.25 0.75])
% axis([1 11 1 17])
plot(input_x,olsslope .* input_x + b1(1,1),'-b');
hold on
% plot([-2 -1 0 1 2 3],[5 5 5 5 5 5],'-r')
% plot(plot_x,rmaslope * plot_x + 1,'-g');
plot(input_x1,olsslope_1 .* input_x1 + b1_1(1,1),'-k');
hold on
plot(input_x2,olsslope_2 .* input_x2 + b1_2(1,1),'-r');
hold on
legend('Category 1','Category 2','Combined','Location','NorthEast')

hold off

display('b1')
b1
b1

b1 =

    0.1000
    1.0500

end