plotting a mean response curve

Apologies - I am a Matlab newbie and have spent a few hours looking for online help with this. All I want to be able to do is to pull in (or generate) numerous "dose-response" datasets (each one a pair of vectors, one "x", the other "R" (the response). The x vector is common to all. I want to: 1. Plot all curves 2. Plot the mean response curve 3. Show standard deviations at each x point

I have tried:

x=[0, 50, 200, 500]  
R1=[10, 30, 90, 200];
R2=[10, 15, 60, 90];
R3=[1, 5, 40, 75];
Rave = mean([R1 R2 R3],2); % assuming Ys are column vectors and defines mean of the response vectors 
plot(x,R1,'b')
hold on
plot(x,R2,'b')
hold on
plot(x,R3,'b')
plot(x, Rave,'-s','MarkerSize',5,'MarkerEdgeColor','red','MarkerFaceColor',[1 .6 .6]) %// mean of all blue curves) %plots the 3 vectors in blue

This is not working (in the sense that the means are not right). I am sure this is simple so apologies in advance but any help hugely welcome

Answer 1

[Star Strider]

‘assuming Ys are column vectors and defines mean of the response vectors’

They are not, at least in the code that you posted. If you want them as column vectors, use semicolons ‘;’ that will concatenate the elements of the vectors vertically instead of commas ‘,’ that will concatenate them horizontally.

This is probably what you want:

x=[0, 50, 200, 500]  
R1=[10; 30; 90; 200];
R2=[10; 15; 60; 90];
R3=[1; 5; 40; 75];
Rave = mean([R1 R2 R3],2);

figure(1)
plot(x,R1,'b')
hold on
plot(x,R2,'b')
hold on
plot(x,R3,'b')
plot(x, Rave,'-sr','MarkerSize',5,'MarkerEdgeColor','red','MarkerFaceColor',[1 .6 .6]) %// mean of all blue curves) %plots the 3 vectors in blue

The other change I made was to plot the line for the mean entirely in red to make it easier to see.

[Jason]

Many thanks for this. - appreciated.

[Jason]

Does anyone know how to then generate the associated standard deviation around the mean?

[Star Strider]

My pleasure.

Use the std function to calculate the standard deviation:

Rstd = std([R1 R2 R3],[],2);

although the most appropriate statistic is the ‘standard error of the mean’, from which you can calculate the confidence intervals if you want:

Rmtx = [R1 R2 R3];
Rsem = std(Rmtx,[],2)/sqrt(size(Rmtx,2));

Use the tinv function with degrees-of-freedom given by:

v = size(Rmtx,2)-1;

since you are calculating only one parameter.

The full code is then:

x=[0, 50, 200, 500]  
R1=[10; 30; 90; 200];
R2=[10; 15; 60; 90];
R3=[1; 5; 40; 75];
Rmtx = [R1 R2 R3];
Rave = mean(Rmtx,2);
Rsem = std(Rmtx,[],2)/sqrt(size(Rmtx,2));
v = size(Rmtx,2)-1;
tstat = tinv(0.95,v);                                       % For 95% Confidence Intervals
Rci = Rave + bsxfun(@times, [-1 1], Rsem*tstat);

figure(1)
plot(x,R1,'b')
hold on
plot(x,R2,'b')
plot(x,R3,'b')
plot(x, Rave,'-sr','MarkerSize',5,'MarkerEdgeColor','red','MarkerFaceColor',[1 .6 .6]) %// mean of all blue curves) %plots the 3 vectors in blue
plot(x, Rci, '--sr','MarkerSize',5,'MarkerEdgeColor','red','MarkerFaceColor',[1 .6 .6])
hold off

I did this from memory, so it would be worthwhile to check it to be sure the calculations are correct.

If my Answer helped solve your problem, please Accept it!

[Jason]

Really appreciate that many thanks. Would you mind explain what the Rmtx command does and also what the line: Rci = Rave + bsxfun(@times, [-1 1], Rsem*tstat); is doing? I don't know what bsxfun, @times mean.

[Jason]

Hi - sorry for the delay - I had to install a toolbox for this to run. It actually doesn't do quite what it wanted. The Rci line seems to plot two dashed red lines spanning the raw data plots and the mean plot. What I actually wanted was a vertical "error bar" shown at each of the mean data points. Maybe I was completely unclear to start off with

[Star Strider]

My pleasure.

The ‘Rmtx’ assignment simply concatenates your ‘R’ vectors to form a matrix. This makes it easier to plot them and to calculate statistics from them.

This assignment:

Rci = Rave + bsxfun(@times, [-1 1], Rsem*tstat);

multiplies the standard error vector ‘Rsem’ by the scalar critical value ‘tstat’, and then expands the ‘[-1 1]’ vector to the correct dimensions to multiply it element-wise. This is what bsxfun does. (The bsxfun function expands the arguments to the appropriate dimensions to do whatever operation is defined in the first argument, here being element-wise multiplication using the ‘times’ function. It is necessary to refer to it by its function handle designated by the ‘@’ sign, so @times.) This creates the 95% confidence intervals, that are then added to the mean vector ‘Rave’ to create the ‘Rci’ matrix for the plot.

See the documentation on bsxfun (link) for a more detailed description of what it can do.

[Jason]

Hi, thanks again. Can I use the product of Rsem and stat to generate a plotted error bar on each of the Rmean values?

[Star Strider]

My pleasure.

To plot errorbars using the errorbar function, the ‘Rci’ calculation changes slightly, and the entire code becomes:

x=[0, 50, 200, 500];  
R1=[10; 30; 90; 200];
R2=[10; 15; 60; 90];
R3=[1; 5; 40; 75];
Rmtx = [R1 R2 R3];
Rave = mean(Rmtx,2);
Rsem = std(Rmtx,[],2)/sqrt(size(Rmtx,2));
v = size(Rmtx,2)-1;
tstat = tinv(0.95,v);                                       % For 95% Confidence Intervals
Rci = bsxfun(@times, [-1 1], Rsem*tstat);

figure(1)
plot(x,R1,'b')
hold on
plot(x,R2,'b')
plot(x,R3,'b')
plot(x, Rave,'-sr','MarkerSize',5,'MarkerEdgeColor','red','MarkerFaceColor',[1 .6 .6])
errorbar(x, Rave, Rci(:,1), Rci(:,2), '-r')             % Plot Error Bars
hold off
set(gca, 'XLim', [-10 510])                             % Expand X-Axis Slightly To Show Errorbars

To include a legend:

x=[0, 50, 200, 500];  
R1=[10; 30; 90; 200];
R2=[10; 15; 60; 90];
R3=[1; 5; 40; 75];
Rmtx = [R1 R2 R3];
Rave = mean(Rmtx,2);
Rsem = std(Rmtx,[],2)/sqrt(size(Rmtx,2));
v = size(Rmtx,2)-1;
tstat = tinv(0.95,v);                                       % For 95% Confidence Intervals
Rci = bsxfun(@times, [-1 1], Rsem*tstat);

figure(1)
h1 = plot(x,R1,'b');
hold on
plot(x,R2,'b')
plot(x,R3,'b')
plot(x, Rave,'-sr','MarkerSize',5,'MarkerEdgeColor','red','MarkerFaceColor',[1 .6 .6])
he = errorbar(x, Rave, Rci(:,1), Rci(:,2), '-r');           % Plot Error Bars
hold off
legend([h1 he], 'Data', 'Mean ± 95% CI', 'Location','NW')
set(gca, 'XLim', [-10 510])                                 % Expand X-Axis Slightly To Show Errorbars

Experiment to get it to look the way you want. See the documentation on errorbar (link) and legend (link) for details.

[Jason]

You are an absolute star - thank you very much.

[Star Strider]

Thank you! As always, my pleasure.